基于SIJRD -系統(tǒng)模型的COVID -19流行研究

李丹陽, 姚斌, 杜雷鳴

(云南師范大學(xué) 物理與電子信息學(xué)院,云南 昆明 650500)

1 背景簡介

2020年1月起，COVID -19逐步在全球范圍內(nèi)流行，已累計(jì)確診40 675 337例，累計(jì)死亡1 119 412人(截至2020年10月20日除中國以外數(shù)據(jù)).宋倩倩等人[1]推測，COVID -19的傳播路徑不是單一的，而且人自身?xiàng)l件的差異，可能導(dǎo)致感染情況不同，所以公開數(shù)據(jù)可能只是確診病例的數(shù)據(jù)，還有一些例如無癥狀感染者自愈及有癥狀感染者非確診死亡等篩查困難的數(shù)據(jù)，稱為隱性數(shù)據(jù)，對(duì)應(yīng)人群稱為隱性患者，這在Yang 等人[2]和Ivorra等人[3]的研究中得到證實(shí)，但他們的研究僅得出無癥狀感染者在病毒傳播中的作用不能忽視，所以本研究利用公開數(shù)據(jù)來估計(jì)隱性數(shù)據(jù)，以探究隱性患者對(duì)整個(gè)疫情的影響，對(duì)未來可能出現(xiàn)的類似疫情以及其他信息傳播的研究提供參考.

2 模型建立

2.1 SIR模型

SIR模型[4-6]是基于微分動(dòng)力學(xué)的最基本的傳染病模型，它將所研究人群分為3類，易感者(S)被感染成為具有傳染性的感染者(I)，感染者(I)死亡或者治愈后移除傳播系統(tǒng)成為移除者(R).設(shè)當(dāng)前環(huán)境總?cè)丝跀?shù)為N(忽略自然出生和死亡、遷入和遷出帶來的影響)，則N=S+I+R.

SIR模型只能描述疫情的大體傳播流行過程，無法描述隱性患者在疫情過程中起到的作用，因此在SIR模型基礎(chǔ)上，基于二叉樹原理，建立SIJRD -系統(tǒng)模型，將人群從4類逐步細(xì)化為8類，增加了4個(gè)求解函數(shù)，提高數(shù)據(jù)擬合度.

2.2 SIJRD -系統(tǒng)模型

SIJRD -系統(tǒng)模型由3個(gè)子模型構(gòu)成，分別是SIRR、SIIR和SIIJRD模型，它們均是以SIR模型作為基礎(chǔ)進(jìn)行延伸的模型.

2.2.1 SIJRD -系統(tǒng)模型架構(gòu)與假設(shè)

將所研究地區(qū)的人群分類(如表1所示).

假設(shè):

(1)當(dāng)前環(huán)境總?cè)丝跀?shù)為N=S+I+R.即整個(gè)研究地區(qū)看作一個(gè)封閉系統(tǒng),忽略自然出生和非COVID -19死亡以及人口遷入、遷出帶來的影響.

(2)移除者對(duì)該病具有終生免疫能力，即可排除在傳染系統(tǒng)外.

(3)感染率為a，即單位時(shí)間內(nèi)一個(gè)感染者傳染給一個(gè)易感者的概率，t時(shí)刻傳染者占總?cè)丝诘腎(t)/N，則t時(shí)刻單位時(shí)間內(nèi)被感染者傳染的易感者的人數(shù)為aI(t)S(t)/N.

Fanelli等人[7]研究得到隔離措施對(duì)疫情控制有效，即實(shí)行隔離措施可以降低感染率a.定義基礎(chǔ)感染率γ為COVID -19的自然傳播速度.隔離率β為所研究地區(qū)的感染者未被有效隔離的概率，當(dāng)β=1時(shí)，說明沒有任何隔離效果，此時(shí)a=γ.感染率a與隔離率β的關(guān)系可以記為a=βγ[8].

(4)定義隱性輸出率b1為感染者I轉(zhuǎn)變?yōu)殡[性移除者R的概率，則t時(shí)刻單位時(shí)間內(nèi)感染者沒有確診隔離而輸出的人數(shù)為b1I(t).定義顯性輸出率c1為感染者I轉(zhuǎn)變?yōu)榇_診輸出者R3的概率，則t時(shí)刻單位時(shí)間內(nèi)感染者I被確診隔離，治愈或死亡的人數(shù)為c1I(t).

(5)定義無癥狀輸出率b為無癥狀感染者I1輸出(死亡和自愈)的概率，則t時(shí)刻單位時(shí)間內(nèi)從無癥狀感染者變成移除者的人數(shù)為bI1(t).定義有癥狀非就診輸出率d為有癥狀感染者I2未被確診而直接輸出(死亡和自愈)的概率，則t時(shí)刻單位時(shí)間內(nèi)從有癥狀感染者I2轉(zhuǎn)變?yōu)橛邪Y狀隱性輸出者R2的人數(shù)為dI2(t).

(6)定義發(fā)病率c為無癥狀感染者I1轉(zhuǎn)變?yōu)橛邪Y狀感染者I2的概率，則t時(shí)刻單位時(shí)間內(nèi)出現(xiàn)癥狀的感染者人數(shù)為cI1(t).其倒數(shù)1/c表示病毒的平均潛伏期τ(單位：d).

表1 當(dāng)前環(huán)境地區(qū)人群分類

(7)定義確診率d1為有癥狀感染者I2被確診成為隔離者J的概率，則t時(shí)刻單位時(shí)間內(nèi)有癥狀感染者I2被確診隔離的人數(shù)為d1I2(t).

(8)定義治愈率e、病死率f分別為被隔離的確診者J治愈成為隔離治愈者R4的概率、死亡的概率，則t時(shí)刻單位時(shí)間內(nèi)被隔離的確診者J治愈的人數(shù)為eJ(t)，死亡的人數(shù)為fJ(t).

由于b1和c1均為I的輸出率，而I類人群每天的輸出的數(shù)量不可能超過它本身的數(shù)量，所以有b1+c1≤1;同理有b+c≤1，d+d1≤1，e+f≤1.

基于SIR模型及假設(shè)(1)-(8)，給出SIJRD -系統(tǒng)模型的流程圖及微分方程組(圖1).

圖1 SIJRD -系統(tǒng)模型流程圖和微分方程

2.2.2 SIJRD -系統(tǒng)模型參數(shù)和初值的確定

確定參數(shù)最常見的方法是輸入假定參數(shù)值用所建模型對(duì)實(shí)際數(shù)據(jù)進(jìn)行擬合，擬合度最高的一組參數(shù)即為最優(yōu)參數(shù).本研究采用最小二乘法進(jìn)行擬合，以可決系數(shù)S2來評(píng)價(jià)擬合效果.這一方法被Batista[6]和Peng等人[9]所采用并證實(shí)可行.確定初值的方法同參數(shù)確定的方法一樣.

2.2.3 基本再生數(shù)

對(duì)傳染病流行的研究可以用基本再生數(shù)R0來評(píng)價(jià)疫情的規(guī)模.基本再生數(shù)R0定義為在自然條件下，一個(gè)傳染者在其整個(gè)傳染期內(nèi)平均傳染的易感者的數(shù)量.R0越大，傳播速率越大，如果R0小于1，則意味著該疾病將逐漸消失[3].根據(jù)SIJRD -系統(tǒng)模型的特點(diǎn)，得到SIJRD -系統(tǒng)模型的基本再生數(shù)的公式：

(1)

當(dāng)沒有任何隔離效果也就是疫情剛開始爆發(fā)時(shí)，即隔離率β=1時(shí)，a=γ.此時(shí)的R0為：

(2)

由(9)式，從R00到R0的變化趨勢，可以判斷疫情的規(guī)模和疫情的發(fā)展情況.

3 SIJRD -系統(tǒng)模型對(duì)國內(nèi)COVID -19疫情的研究

本研究的全部數(shù)據(jù)均來源于中國國家衛(wèi)生健康委員會(huì)和各省市區(qū)衛(wèi)生健康委員會(huì)的官網(wǎng).將得到的相關(guān)數(shù)據(jù)(從2020 年1月27日至2020年3月20日的國內(nèi)各地每日確診、死亡、治愈人數(shù)以及2019年總?cè)丝跀?shù))進(jìn)行整理.考慮到國內(nèi)的行政區(qū)劃較多，不同行政區(qū)經(jīng)濟(jì)發(fā)展程度不同，則疫情傳播流行情況也各不相同，因此依據(jù)官方報(bào)道的疫情數(shù)據(jù)利用聚類分析方法對(duì)行政區(qū)進(jìn)行分類.由于港澳臺(tái)和湖北的前期數(shù)據(jù)不完整，西藏和青海的COVID -19流行不明顯，所以在進(jìn)行分類研究時(shí)選擇將這些地區(qū)排除在外，將其他28個(gè)行政區(qū)分為三類,分類結(jié)果如表2.

表2 28個(gè)行政區(qū)分類情況(除港澳臺(tái)、青海、西藏、湖北以外的行政區(qū))

根據(jù)表2的分類，對(duì)三類地區(qū)進(jìn)行參數(shù)和初值確定.圖2-4分別是A、B、C三類地區(qū)的確診、確診隔離、死亡和治愈人數(shù)自2020年1月27日至3月20日的數(shù)據(jù)擬合結(jié)果，其中“*”號(hào)表示實(shí)際數(shù)據(jù)，實(shí)線為擬合曲線.

圖2 A類地區(qū)模型擬合圖

圖3 B類地區(qū)模型擬合圖

圖4 C類地區(qū)模型擬合圖

結(jié)合圖2-4及其對(duì)應(yīng)的可決系數(shù)可以看出，SIJRD -系統(tǒng)模型能夠較好地?cái)M合每日確診人數(shù)，但對(duì)死亡人數(shù)、確診隔離人數(shù)、治愈人數(shù)的擬合度相對(duì)較差，這可能是因?yàn)榉酪叽胧┑募訌?qiáng)[10]會(huì)導(dǎo)致疫情發(fā)展出現(xiàn)分隔期，分隔期前后的疫情發(fā)展不同，使得部分參數(shù)發(fā)生改變.

表3和表4為擬合及計(jì)算得到的三類地區(qū)疫情的初值和參數(shù)值(COVID -19的基礎(chǔ)感染率γ大約為1[9]).

表3 三類地區(qū)COVID -19疫情的初值

表4 三類地區(qū)COVID -19疫情的參數(shù)值

利用得到的表3和表4中的初值和參數(shù)使用SIJRD -系統(tǒng)模型對(duì)每一類地區(qū)進(jìn)行各類人數(shù)預(yù)測，由于傳播并未涉及大部分人，由此把易感人群排除，通過圖2中微分方程組，得到三類地區(qū)的其他類人群人數(shù)—時(shí)間變化圖(圖5).

圖5 三類地區(qū)的其他類人群預(yù)測人數(shù)

由表4： C類地區(qū)的感染率a最小，病死率f最小，治愈率e最大，說明C類地區(qū)防控措施實(shí)行較好導(dǎo)致a小，醫(yī)療條件較優(yōu)導(dǎo)致其病死率f低而治愈率e高.C類地區(qū)包括廣東和浙江等經(jīng)濟(jì)發(fā)達(dá)省份，還包括防疫措施到位的河南等省份，與實(shí)際情況相符.三類地區(qū)的病毒潛伏期τ分別為τA=3.414 1 d,τB=3.752 3 d,τC=4.472 2 d，均在現(xiàn)有研究得到的該病毒的潛伏期2-14 d[2]內(nèi).B類地區(qū)的R00最大，即疫情剛開始時(shí)，B類地區(qū)疫情規(guī)模較大，傳播速率較快.三類地區(qū)的R0均小于1，說明各地的防疫措施行之有效，國內(nèi)的COVID -19在逐步消亡.且C類地區(qū)的R0減少量是最多的，進(jìn)一步說明C類地區(qū)的防疫措施實(shí)行最有效.

由圖5：三類地區(qū)的各類人群變化趨勢相似，無癥狀感染人數(shù)I1隨時(shí)間在減少，說明易感者被感染的人數(shù)減少，有癥狀感染人數(shù)I2的人數(shù)在2月初就達(dá)到峰值.確診隔離人數(shù)J的變化趨勢圖描繪的是確診人數(shù)輸入與輸出的情況，A類和C類地區(qū)的確診隔離者J在2月8日達(dá)到峰值，B類在2月10日達(dá)到峰值，說明這三類地區(qū)COVID -19疫情的輸出均早已大于輸入，疫情發(fā)展已逐漸良化.3月初確診人數(shù)(R3)和治愈人數(shù)(R4)增加變緩，即每日新增確診和新增治愈均在減少，說明此時(shí)三類地區(qū)疫情基本得到控制.A、B和C三類地區(qū)的治愈人數(shù)的變化曲線與確診隔離人數(shù)的變化曲線分別相交于2月15日、2月17日和2月15日.同官方所得結(jié)果：2月12日新增確診病例達(dá)到峰值，2月18日新增出院病例數(shù)開始超過新增確診病例數(shù)[10]，相差不大.說明SIJRD -系統(tǒng)模型能夠較好地描述確診隔離人數(shù)的整體變化趨勢.C類地區(qū)的I1、I2、R1和R2在2月25日左右趨于平和，變化微小，A類和B類地區(qū)的I1、I2、R1和R2則在2月底3月初趨于平和，變化微小，與實(shí)際數(shù)據(jù)走向吻合.

SIJRD -系統(tǒng)模型預(yù)測出截至2020年3月20日與隱性患者有關(guān)的如下比例：

(1)A類、B類和C類地區(qū)隱性患者(R1+R2)占總病例的50.00%、50.29%和50.32%.

(2)A類、B類和C類地區(qū)的無癥狀輸出病例(R1)占總病例的21.38%、25.59%和30,32%；有癥狀非就診輸出病例(R2)占總病例的15.25%、12.8%和8.21%；確診輸出病例(R3)占總病例的63.37%、61、61%和61.47%.

可見在國內(nèi)雖然確診輸出是感染者主要輸出路徑，但非就診輸出也占較大的比例.

4 全球疫情發(fā)展預(yù)測

基于對(duì)國內(nèi)COVID -19疫情的研究，利用SIJRD -系統(tǒng)模型擬合全球相關(guān)數(shù)據(jù)獲得相應(yīng)的初值和參數(shù)，預(yù)測全球的疫情發(fā)展趨勢(2020年5月23日—2021年6月26日，共400 d).預(yù)測圖如圖6所示.分別用虛線和實(shí)線代表模型預(yù)測的R3和J.

圖6 全球新冠肺炎疫情發(fā)展預(yù)測圖

由圖6可知，全球的確診隔離人數(shù)J將在2020年10月5日(圖中預(yù)測第136天)達(dá)到峰值，即在此之后各地的輸入病例將小于輸出病例，意味著開始疫情良化，2021年5月以后(圖中預(yù)測350天后)J增長極其緩慢，表明疫情得到基本控制.模型預(yù)測在2020年10月20日時(shí)，全球累計(jì)確診人數(shù)達(dá)到39 453 869人，而實(shí)際在10月20日時(shí)，全球累計(jì)確診人數(shù)為40 675 337人，預(yù)測與實(shí)際誤差約為0.03，SIJRD -系統(tǒng)模型對(duì)全球COVID -19疫情趨勢預(yù)測較為準(zhǔn)確.

5 小結(jié)

由于SIJRD -系統(tǒng)模型是以SIR模型為基礎(chǔ)將其R進(jìn)行細(xì)化至SIRR模型，進(jìn)而不斷細(xì)化至SIIJRD模型，所以SIJRD -系統(tǒng)模型描述的實(shí)際還是SIR模型描述的各類人群數(shù)據(jù)走向；較SIR模型而言，細(xì)化后的SIJRD -系統(tǒng)模型能夠找到導(dǎo)致這些人群數(shù)據(jù)走向的原因；由得到的隱性患者變化趨勢及相關(guān)比例推知國內(nèi)確診輸出是感染者主要輸出的途徑，但非就診輸出也占較大比例，即隱性患者在疫情傳播中起的作用不可忽視，其中未檢測出的病例占總病例的50%左右，與Ivorra等人[3]研究的結(jié)果接近.

SIJRD -系統(tǒng)模型預(yù)測COVID -19在全球的疫情將繼續(xù)發(fā)展.2021年5月以后(圖中預(yù)測350天后)得到基本控制.

綜上，SIJRD -系統(tǒng)模型可以對(duì)各地實(shí)施措施的有效性進(jìn)行探討和評(píng)估，模擬隱性人群的數(shù)據(jù)走向，探究無癥狀感染者在疫情傳播中起到的作用，同時(shí)也能預(yù)測COVID -19的后期發(fā)展趨勢，對(duì)今后可能出現(xiàn)的類似傳染病的疫情防控有指導(dǎo)作用.