一類新冠肺炎的流行病學(xué)模型穩(wěn)定性研究

金德泉,陳芃合,王凱明

(1.廣西大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,廣西 南寧 530004;2.廣西大學(xué)廣西應(yīng)用數(shù)學(xué)中心,廣西 南寧 530004;3.長安大學(xué)理學(xué)院數(shù)學(xué)與信息科學(xué)系,陜西 西安 710061)

1 引言

目前,新冠肺炎疫情在全球范圍內(nèi)大規(guī)模流行,對各國的社會經(jīng)濟(jì)生活造成了嚴(yán)重的影響[1-2].2020年3月,世界衛(wèi)生組織(WHO)正式宣布新冠肺炎為大流行病.因此,根據(jù)新冠肺炎傳播的流行病學(xué)特點(diǎn),建立符合新冠肺炎傳播規(guī)律的數(shù)學(xué)模型,刻畫新冠肺炎的內(nèi)在傳播機(jī)制,并進(jìn)一步定性定量的研究其傳播規(guī)律,將對制定和評估防控措施方面具有重要的指導(dǎo)作用.

研究傳染病傳播模型的主要方法包括SIS模型、SIR模型、SEIR模型與SEIS模型等.文獻(xiàn)[3]利用離散時(shí)間的SIS模型,模擬了新冠肺炎變異的傳播過程.文獻(xiàn)[4]利用經(jīng)典的SIR模型與支持向量機(jī)方法,對印度尼西亞新冠肺炎的傳播進(jìn)行了預(yù)測.文獻(xiàn)[5]通過對SEIR模型進(jìn)行改進(jìn),即使用SEIQRDP模型對南非新冠肺炎的發(fā)展趨勢進(jìn)行了擬合.文獻(xiàn)[6]基于廣義SEIR模型,預(yù)測了新冠肺炎在中國的傳播趨勢.文獻(xiàn)[7]研究了具有飽和接觸率的SEIS模型的動力學(xué)性質(zhì),得到了決定疾病絕滅或持續(xù)生存的基本再生數(shù).文獻(xiàn)[8]討論了一類具有常數(shù)輸入且傳染率為非線性的SEIS流行病傳播數(shù)學(xué)模型,給出了疾病消失和成為地方病的基本再生數(shù).

以上研究成果在很多方面成功刻畫了新冠肺炎的一些傳播特點(diǎn),為新冠肺炎的防控給出了理論依據(jù).但是SIS模型與SIR模型沒有充分考慮潛伏期對于新冠肺炎傳播的影響,難以很好地刻畫此次新冠病毒存在潛伏期長的特性.同時(shí),與傳統(tǒng)的天花、流感、麻疹這一類流行病不同,目前已發(fā)現(xiàn)的新冠肺炎患者被治愈后依舊存在再次感染的可能性,而接種疫苗也只能降低新冠肺炎的因病死亡率和重癥率,獲得抗體后依舊可能感染新冠肺炎.另外,新冠肺炎具有很強(qiáng)的變異性,這也顯著削弱了抗體在預(yù)防新冠病毒感染中的作用.但是一般意義下的SEIS模型與SEIR模型難以準(zhǔn)確地描述新冠肺炎的這些特點(diǎn).綜上所述,針對新冠病毒因病死亡率率高、傳染力強(qiáng)、潛伏期長、抗體在預(yù)防感染中作用不顯著等特點(diǎn)[9],需要建立一類能夠反映這些特點(diǎn)的流行病學(xué)模型來恰當(dāng)?shù)孛枋龃舜涡鹿诜窝椎膭恿W(xué)特征.

通過建立一類描述新冠肺炎傳播動力學(xué)機(jī)制的改進(jìn)SEIS模型,證明了新冠肺炎傳播的無病平衡點(diǎn)和地方病平衡點(diǎn)全局漸近穩(wěn)定性的條件,并通過數(shù)值仿真的方法進(jìn)行了分析驗(yàn)證.

2 預(yù)備知識

2.1 常微分方程穩(wěn)定性理論

定義2.1[10]設(shè)平面自治系統(tǒng)

且O(0,0)是系統(tǒng)(1)的平衡點(diǎn),則線性系統(tǒng)

引理2.1[11]設(shè)O(0,0)是幾乎線性系統(tǒng)(1)的初等奇點(diǎn),則當(dāng)O(0,0)是其線性近似系統(tǒng)(2)的鞍點(diǎn)、結(jié)點(diǎn)或焦點(diǎn)時(shí),它也必是系統(tǒng)(1)的鞍點(diǎn)、結(jié)點(diǎn)或焦點(diǎn),且具有相同的穩(wěn)定性.

引理2.2[12]設(shè)高維線性系統(tǒng)

的系數(shù)矩陣 A 的特征根為 λ1,λ2,···,λr,則有

(i)若 λ1,λ2,···,λr均具有負(fù)實(shí)部,則系統(tǒng) (3)的零解是漸近穩(wěn)定的;

(ii)若λ1,λ2,···,λr中至少有一個(gè)具有正實(shí)部,則系統(tǒng) (3)的零解是不穩(wěn)定的.

若 λ1,λ2,···,λr中沒有正實(shí)部的根但是有零根或零實(shí)部的純虛根,則當(dāng)零根或零實(shí)部的純虛根的初級因子都是一次的,系統(tǒng)(3)的零解是穩(wěn)定的.當(dāng)零根或零實(shí)部的純虛根中至少有一個(gè)的初級因子次數(shù)大于一次時(shí),系統(tǒng)(3)的零解是不穩(wěn)定的.

2.2 奇異攝動法

引理2.3[13]對攝動系統(tǒng)

若系統(tǒng)(4)滿足

(i)它的退化系統(tǒng)

其中g(shù)(x,y,ε)=0是一條曲線,即存在函數(shù)y=φ(x),使得g(x,φ(x),0)≡0;

(ii)f(x,y,ε)和 g(x,y,ε)關(guān)于 x,y,ε都有連續(xù)光滑導(dǎo)數(shù);

(v)矩陣gy(x*,y*,0)的特征值的實(shí)部都是負(fù)的,則如果完全系統(tǒng)(4)具有局部穩(wěn)定的平衡點(diǎn),那么當(dāng)ε充分小的時(shí)候該平衡點(diǎn)是全局穩(wěn)定的.

2.3 傳染病模型

SIS傳染病模型將人口分為兩類:易感者類,記為類S,該人群目前無病,且無免疫能力,可以被傳染而患病;患者類,記為類I,該人群已經(jīng)患病,且可以把疾病傳染給類S人群.用S(t)和I(t)分別表示t時(shí)刻上述兩類人群占總?cè)丝诘谋壤?SIS模型的假設(shè)條件為:

(b)人群中兩類成員均勻分布,傳播方式為接觸傳播.單位時(shí)間內(nèi)一個(gè)患者類I的成員和其他類S成員接觸率為常數(shù)λ,與類S成員的接觸總數(shù)為λS,則單位時(shí)間內(nèi)與類S成員接觸總數(shù)為λNSI,即單位時(shí)間內(nèi)類I成員增加的數(shù)量,稱為發(fā)病率.

SIS傳染病模型假設(shè)患者病愈后并未產(chǎn)生免疫能力,仍為S類成員.因此,SIS模型[15]為

由S(t)和I(t)的定義及人口數(shù)為常數(shù)可知S+I=1,則模型(6)可化簡為

將方程(7)改寫為

當(dāng)傳染病具有潛伏期時(shí),可以將處于潛伏期人群比例設(shè)為E;同時(shí),如果治愈或接種疫苗后的人群具有長期免疫力,可將其比例設(shè)為R.由此可以得到SIS模型的改進(jìn)形式,如SEIS和SEIR等流行病模型.

3 SEIS模型及其平衡點(diǎn)的全局穩(wěn)定性

新冠肺炎具有不同于傳統(tǒng)傳染病的一些特點(diǎn).首先,在接種新冠疫苗之前,新冠肺炎發(fā)病導(dǎo)致的死亡率很高;其次,新冠肺炎有一個(gè)很長的潛伏期,且在潛伏期中具有傳染性;第三,個(gè)體在感染新冠和接種了新冠疫苗仍然有可能感染新冠肺炎.對于具有潛伏期,且不具備長期免疫性特點(diǎn)的傳染病,常用的模型是SEIS模型.在傳統(tǒng)的SEIS模型中,E代表處于潛伏期的人群比例.但是,新冠肺炎具有與傳統(tǒng)傳染病不同的流行病學(xué)特點(diǎn).首先,由于大流行背景下檢測條件的限制,潛伏期感染者自愈及轉(zhuǎn)化為患者的比例難以精確統(tǒng)計(jì).其次,新冠肺炎的潛伏期感染者具有一定的傳染性,但難以確定其接觸感染率.最后,新冠肺炎的潛伏期和發(fā)病期的界定仍存在爭議,界限不夠清晰,這些都對建立較為準(zhǔn)確的新冠肺炎傳播的動力學(xué)模型帶來了較大的困難.

為了解決這些問題,該模型對處于新冠肺炎潛伏期的感染者人群進(jìn)行更為細(xì)致的劃分.首先,由于新冠肺炎病毒變異性強(qiáng),同時(shí)潛伏期感染者自愈后所帶抗體作用時(shí)間短,將感染后未發(fā)病、未檢出,不具傳染性的自愈人群歸入類S.其次,將感染了新冠肺炎,通過咽拭子檢測未檢出,但最終發(fā)展為新冠肺炎患者的人群歸入類E.最后,將處于潛伏期,但可以通過咽拭子檢出的感染者人群歸入類I.由此得到如下關(guān)于新冠肺炎傳播的改進(jìn)SEIS模型:

在該模型中,λ的定義修改為單位時(shí)間內(nèi)一個(gè)類I的成員和類S成員接觸感染率,即接觸并造成感染的比例.k是單位時(shí)間內(nèi)E類人群中的個(gè)體由潛伏期的感染者發(fā)展為新冠肺炎患者的比例,即發(fā)病率.d是單位時(shí)間內(nèi)類I的人群因?yàn)樾鹿诜窝滓鹚劳龅谋壤?即因病死亡率.其他都如2.3節(jié)中定義.則由計(jì)算可得,基本再生數(shù)為

Rc0的流行病學(xué)意義是,在病期內(nèi)每個(gè)新冠肺炎患者平均感染其他健康者的人數(shù).為了方便分析模型的動力學(xué)性質(zhì),將模型(9)化簡為

其中a=k+μ,b=γ+μ+d.

在以下的討論中,將結(jié)合奇異攝動理論[14],研究系統(tǒng)(10)平衡點(diǎn)的全局穩(wěn)定性.奇異攝動理論允許一些非攝動系統(tǒng)通過尺度變化為攝動系統(tǒng).下面利用Lyapunov函數(shù)和Dulac函數(shù),以及奇異攝動方法和理論證明系統(tǒng)平衡點(diǎn)的全局穩(wěn)定性.

設(shè)D={(S,E,I)|S＞0,E＞0,I＞0,S+E+I≤ 1}.

命題 3.1D為系統(tǒng)(10)的正向不變集,且全局吸引.

3.1 系統(tǒng)變量的尺度變換

首先對時(shí)間進(jìn)行尺度變換,令τ=kt.再分別對S,E和I尺度變換,令

則模型(10)變?yōu)?/p>

令ε=0,則得到退化系統(tǒng)

下面將對退化系統(tǒng)(14)-(16)進(jìn)行討論.

3.2 退化系統(tǒng)的平衡點(diǎn)性質(zhì)

通過計(jì)算可以得到系統(tǒng)(17)具有兩個(gè)平衡點(diǎn)(μ,0),

根據(jù)流行病學(xué)性質(zhì),前者稱為無病平衡點(diǎn),后者稱為地方病平衡點(diǎn).為了證明系統(tǒng)(17)平衡點(diǎn)的全局性質(zhì),本文先證明系統(tǒng)(17)平衡點(diǎn)的局部性質(zhì).

證明(1)令 x=X+μ,y=Y,則系統(tǒng)(17)變?yōu)?/p>

它的線性近似系統(tǒng)為

系統(tǒng)(19)的系數(shù)矩陣為

因此,它的線性近似系統(tǒng)為

它的系數(shù)矩陣為

特征方程為m2+pm+q=0,其中

當(dāng)Rc0＞1,p和q都是正的.因此系統(tǒng)(17)的無病平衡點(diǎn)是局部漸近穩(wěn)定的.

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

命題3.4當(dāng)Rc0＜1時(shí),系統(tǒng)(10)的無病平衡點(diǎn)是全局漸近穩(wěn)定的.

3.4 系統(tǒng) (10)地方病平衡點(diǎn)的全局穩(wěn)定性

下面證明系統(tǒng)(13)滿足定理2.3.

命題3.5如果Rc0＞1且ε充分小,系統(tǒng)(13)的地方病平衡點(diǎn)是全局漸近穩(wěn)定的.

證明對該模型,有

可得

它的特征方程為m2+pm+q=0.其中當(dāng)Rc0-1＞0時(shí),

所以J當(dāng)Rc0-1只有負(fù)的特征值時(shí),滿足條件(iii).

由命題3.1,可知只要t足夠大,S+E+I≤1,即命題3.3成立,完全系統(tǒng)的地方病平衡點(diǎn)是局部漸近穩(wěn)定的.所以由定理2.3可知,如果Rc0＞1且ε充分小,系統(tǒng)(13)的地方病平衡點(diǎn)是全局漸近穩(wěn)定的.

4 參數(shù)估計(jì)與數(shù)值模擬

4.1 參數(shù)估計(jì)

4.1 數(shù)值模擬

圖1 f=0.01時(shí)x和y的軌線,此時(shí)ε=0.0007,Rc0=0.05

圖2 f=0.05時(shí)x和y的軌線,此時(shí)ε=0.00037,Rc0=0.24

圖3 f=0.2時(shí)x和y的軌線,此時(shí)ε=0.00175,Rc0=0.95

圖4 f=0.3時(shí)x和y的軌線,此時(shí)ε=0.003,Rc0=1.4

由上面的數(shù)值實(shí)驗(yàn)可以看到,當(dāng)Rc0＜1時(shí),系統(tǒng)(13)和系統(tǒng)(17)的軌線都收斂于無病平衡點(diǎn)(0.014,0),當(dāng)ε很小時(shí)符合得相當(dāng)?shù)暮?但隨著ε的增大,兩個(gè)系統(tǒng)的軌線的差距開始增大.當(dāng)取到f=0.3時(shí),軌線不再收斂于無病平衡點(diǎn),而是收斂于地方病平衡點(diǎn)(0.01,0.0039).

為了觀察地方病平衡點(diǎn)的性質(zhì),取 λ=30,這在實(shí)際情況中已經(jīng)是很大的值.當(dāng) f=0.1時(shí),則

此時(shí)地方病平衡點(diǎn)為(0.0099,0.0041).

圖5 f=0.1時(shí)x和y的軌線,此時(shí)ε=0.000778,Rc0=1.42

圖6 f=0.3時(shí)x和y的軌線,此時(shí)ε=0.003,Rc0=4.26

通過上面的實(shí)驗(yàn)可以看到,在上述幾種情況下,Rc0＜1時(shí)模型的無病平衡點(diǎn)和Rc0＞1時(shí)模型的地方病平衡點(diǎn)的穩(wěn)定性是明顯的,退化系統(tǒng)和完全系統(tǒng)在平衡點(diǎn)附近的性質(zhì)是一致的,在一定條件下可以運(yùn)用退化系統(tǒng)來模擬完全系統(tǒng)的性質(zhì).這與該SEIS模型的理論證明是一致的.

5 結(jié)論

疾病的傳染率、潛伏期、因病死亡率是傳染病流行的重要因素,影響著流行病的動力學(xué)性質(zhì),因此對具有類似新冠肺炎傳播特點(diǎn)的動力學(xué)傳染病模型進(jìn)行定性分析具有現(xiàn)實(shí)意義.由以上的證明結(jié)果和仿真實(shí)驗(yàn)可知,要防止新冠肺炎的流行,必須減少基本再生數(shù)Rc0.若Rc0＜1,即一個(gè)患者在平均患病期能傳染的最大人數(shù)小于1時(shí),患者和潛伏期內(nèi)感染者的人數(shù)會逐漸減少,疾病自然逐漸消亡;反之,若Rc0＞1,即一個(gè)患者在平均患病期能傳染的最大人數(shù)大于1時(shí),患者人數(shù)會逐漸增加,但不會無限增加,在一定時(shí)間后,易感染人群比例、潛伏期感染者比例、患者比例會趨于一個(gè)平衡狀態(tài),疾病將始終存在而形成地方病.

為了減小基本再生數(shù),可以采取的方法包括:

(1)減少λ,即患者接觸其他健康者的數(shù)量.對患者進(jìn)行隔離可以達(dá)到此效果;

(2)提高治愈率γ,即縮短治療時(shí)間.這可以通過提高醫(yī)療條件達(dá)到此效果;

(3)減少k,即減少感染者的發(fā)病率;這可以通過研究相應(yīng)的疫苗達(dá)到此效果.

(4)當(dāng)因病死亡率d增大時(shí),基本再生數(shù)Rc0也會減少.這表明,新冠病毒雖然存在出現(xiàn)毒性加強(qiáng)的變異而導(dǎo)致因病死亡率過高的情況,但由于死亡率過高降低了Rc0,該變異類型很可能只能在局部地區(qū)造成疫情加重的情況,難以在大范圍地區(qū)形成大規(guī)模流行.

由前面的討論可以看到,在沒有特效藥的情況下,即使人群普遍接種了疫苗,如果放任人群在無防護(hù)情況下發(fā)生大規(guī)模聚集和流動,依然存在基本再生數(shù)過高,無法根除新冠肺炎流行,從而使其發(fā)展為在一定范圍內(nèi)長期流行的傳染病的風(fēng)險(xiǎn).因此,在國外新冠肺炎仍處于大流行的背景下,在研制出防護(hù)能力較強(qiáng)的高效疫苗和特效藥之前,仍然需要采用包括佩戴口罩、控制跨國人口流動等措施來降低接觸感染率,維持國內(nèi)來之不易的較好的防疫形勢,保證國內(nèi)社會生產(chǎn)和人民生活的正常運(yùn)行.

值得注意的是,該SEIS模型是在簡化了一些情況和條件下得到的,這使得模型的預(yù)測能力受到了一些影響.為此,在進(jìn)一步的工作中,可以對人群進(jìn)行更細(xì)致的分類,進(jìn)一步完善模型,以此提高模型的預(yù)測能力和預(yù)測精確性.