基于多時(shí)段變參數(shù)SEIQR模型的流感病毒傳播研究

何洋文，周大勇，方 銘，馬永峰

(大連交通大學(xué) a.土木工程學(xué)院；b.理學(xué)院；c.機(jī)械工程學(xué)院，遼寧 大連 116028)

公共衛(wèi)生事業(yè)是公共事業(yè)的重要組成部分，與人們的健康、學(xué)習(xí)、工作及生活有著密切的關(guān)系[1]。流感是一種由流感病毒引起的呼吸道傳染病，傳染性強(qiáng)，傳播面廣，極易在青少年和老年人等免疫力低的人群中傳播，特別是在學(xué)校這種相對(duì)密集和封閉的環(huán)境中，傳染病的高發(fā)期極易導(dǎo)致病毒的大范圍傳播[2]，因此對(duì)于傳染病的研究尤為重要。

利用信息技術(shù)對(duì)傳染病進(jìn)行預(yù)防和預(yù)測(cè)是公共衛(wèi)生事業(yè)非常重要的部分，可通過(guò)建立多時(shí)段變參數(shù)微分方程模型來(lái)研究突發(fā)流感的演化規(guī)律。

1 模型建立與求解分析

1.1 模型假設(shè)

對(duì)病毒傳播過(guò)程中的部分因素進(jìn)行假設(shè)(包括人數(shù)變化和采取的措施等)，以達(dá)到簡(jiǎn)化模型的效果，為模擬疫情傳播提供基礎(chǔ)條件。

假設(shè)在疫情傳播階段，學(xué)校學(xué)生總?cè)藬?shù)不變，分為易感者、潛伏者、感染者、康復(fù)者、隔離者。假設(shè)：傳播過(guò)程中每個(gè)學(xué)生被感染后的發(fā)病時(shí)間相同。在校外不發(fā)生傳播、不存在超級(jí)傳播者、每個(gè)學(xué)生自愈能力相同。

表1 符號(hào)說(shuō)明Tab.1 Symbol description

1.2 多時(shí)段變參數(shù)SEIQR模型

1.2.1 建立多時(shí)段變參數(shù)SEIQR模型

采用傳染病動(dòng)力學(xué)模型，結(jié)合疫情傳播情況，甲型H1N1流感病毒在傳播初期一般不易被察覺(jué)，此時(shí)的傳播特征符合傳統(tǒng)SEIR模型的傳播特點(diǎn)[3]，但是傳統(tǒng)的SEIR模型建立在外力無(wú)干預(yù)的情形下不能處理采取一定防控措施帶來(lái)的變化[4]。但疫情發(fā)展到一定階段時(shí)，人們會(huì)采取防控措施，并會(huì)時(shí)刻根據(jù)疫情的發(fā)展不斷調(diào)整防控措施。針對(duì)疫情傳播特點(diǎn)，建立多時(shí)段變參數(shù)的SEIQR傳染病模型，根據(jù)各時(shí)期的傳播特點(diǎn)改變相關(guān)參數(shù)，以獲得更加真實(shí)合理的模擬效果，得出更加準(zhǔn)確的結(jié)果。關(guān)于SEIQR模型分析如下：根據(jù)流感病毒的流行病學(xué)特征，將人群分為S、E、I、Q、R，建立(S-E-I-Q-R，SEIQR)模型[3]。具體傳播流程如圖1、圖2所示。

其中，β為傳染率，即易感者被感染為潛伏者的概率。

圖1 流感病毒傳播初期的SEIQR模型流程圖Fig.1 Flow chart of the SEIQR model of early transmission of influenza virus

圖2 流感病毒傳播后期的SEIQR模型流程圖Fig.2 Flow chart of the SEIQR model of late transmission of influenza virus

不同概率計(jì)算過(guò)程如下：β：與病毒傳播性質(zhì)、傳播時(shí)間有關(guān)。β1和β2參考具體實(shí)例。a：假設(shè)平均潛伏期為A，a=1/A。隔離率t：假設(shè)感染者平均隔離時(shí)間為B，g=1/B?？祻?fù)率p：假設(shè)平均發(fā)病時(shí)間為C，p=1/C。

SEIQR模型動(dòng)態(tài)全過(guò)程分析如下[4]：已知總N=S+E+I+Q+R人群保持不變。某一時(shí)刻t，根據(jù)不同人群發(fā)生變化建立模型：

已知SEIQR模型。根據(jù)t時(shí)刻各類(lèi)人群人數(shù)，可預(yù)測(cè)(t+1)時(shí)刻各類(lèi)人群人數(shù)：

感染人群：I(t+1)=aE(t)-gI(t)

隔離人群：Q(t+1)=gI(t)-pQ(t)

康復(fù)人群：R(t+1)=pQ(t)

1.2.2 采取隔離措施的SEIQR模型

甲型H1N1流感病毒有潛伏期短和傳染性強(qiáng)的特點(diǎn)，主要發(fā)生于兒童、老人等抵抗力弱的群體[7]。在疫情發(fā)展的中后期，通常會(huì)采取一定的隔離防控措施，阻斷病毒傳播。由于E在疫情中難以被發(fā)現(xiàn)，所以隔離措施只針對(duì)I。模擬疫情不同時(shí)期采取的隔離措施，分析疫情變化情況。隔離期間，假設(shè)病例病情發(fā)展過(guò)程與未隔離者相同，康復(fù)率仍為p，則t時(shí)刻，從Q中康復(fù)人數(shù)為pQ，假設(shè)康復(fù)者永久免疫[8]。

1.2.3 采取隔離+疫苗接種的SEIQR模型

冬春季節(jié)是流感的高發(fā)季節(jié)，尤其需要注意兒童的防護(hù)。其中接種疫苗是預(yù)防病毒感染的關(guān)鍵，雖然疫苗不能100%預(yù)防，但不接種疫苗被感染的概率則是0或100%[9]。通過(guò)模擬在采取一定隔離措施的情況下，接種疫苗對(duì)疫情防控的作用和在不同接種率下疫苗對(duì)疫情的有效控制程度，以反映疫苗在疫情防控中的重要性和必要性。

1.2.4 甲型H1N1流感病毒模型參數(shù)估計(jì)

通過(guò)查閱文獻(xiàn)資料，根據(jù)模擬事件的數(shù)據(jù)可知，甲型H1N1流感病毒的平均潛伏期為2 d，病程為3～5 d[10]，結(jié)合事件數(shù)據(jù)，綜合選擇潛伏期為2 d，病程為4 d，β2=0.125，取其潛伏期的傳染率為感染期的86%，即a=0.5、p=0.25、k=0.86、β1=0.108。

2 模型的評(píng)估與討論

2.1 甲型H1N1流感病毒傳播事件模擬

此次模擬的為某省較偏遠(yuǎn)地區(qū)小學(xué)的一起流感病毒感染事件[10]，主要傳播途徑為人際傳播，通過(guò)分析每日新增感染者人數(shù)(如表2所示)，可將流行曲線分為3個(gè)階段。第1階段為6月13—26日，為病毒傳播初期，26日開(kāi)始調(diào)查采樣，之前未采取任何防護(hù)措施；第2階段為6月27—28日，此時(shí)實(shí)驗(yàn)室確定病毒為甲型H1N1，并落實(shí)控制措施；第3階段為6月28日—7月10日，此時(shí)控制措施生效，新增感染者逐日下降。

表2 學(xué)校每日新增數(shù)據(jù)Tab.2 New data added daily in school

基于SEIQR算法，運(yùn)用MATLAB軟件進(jìn)行擬合優(yōu)度計(jì)算，擬合效果如圖3所示。計(jì)算得曲線擬合結(jié)果R2=0.863，P=0.000，證明此次擬合曲線與實(shí)際情況擬合效果較好。

圖3 甲型H1N1流感爆發(fā)疫情與模擬擬合結(jié)果Fig.3 Influenza A (H1N1) outbreak and simulation results

2.2 采取不同防護(hù)措施下的模型分析

2.2.1 無(wú)任何防護(hù)措施下的疫情分析

在沒(méi)有任何防護(hù)措施的情況下，疫情的傳播模型為SEIR模型。分析模擬結(jié)果可知，若此次疫情未采取任何防護(hù)措施，則共有565人感染流感病毒，患病率為63.13%，疫情持續(xù)40 d，在疫情發(fā)生的第21 d，新增病例數(shù)最多為47例，如圖4所示。實(shí)際該校累計(jì)報(bào)告病例110例，患病率為12.29%。對(duì)比可知，在無(wú)防護(hù)措施的情況下，患病率要增加50%以上，疫情多持續(xù)13 d左右。

圖4 無(wú)防護(hù)措施下甲型H1N1流感疫情模擬Fig.4 Simulation of influenza A (H1N1) epidemic without protection

2.2.2 采取及時(shí)隔離措施的疫情分析

分析圖5模擬結(jié)果可知，如果CDC在調(diào)查取樣當(dāng)天僅實(shí)施隔離，疫情的峰值和感染人數(shù)都會(huì)比實(shí)際要高；如果CDC在疫情傳播第10 d采取隔離措施，則疫情的峰值會(huì)更早到來(lái)且處于一個(gè)較低的狀態(tài)。越早實(shí)施隔離措施，疫情持續(xù)的時(shí)間和峰值到來(lái)的時(shí)間越短，對(duì)于疫情的控制也愈有利。

圖5 隔離措施下甲型H1N1流感疫情模擬Fig.5 Simulation of A (H1N1) influenza epidemic under quarantine measures

2.2.3 采取隔離+接種疫苗的疫情分析

模擬在第14 d采取隔離措施，分析學(xué)生的疫苗接種率分別為10%、30%和50%時(shí)的疫情傳播情況。分析圖6的模擬結(jié)果可知，疫苗接種率為30%時(shí)，疫情就基本被控制住，且隨著疫苗接種率的提高，疫情的發(fā)展得到了明顯抑制。由于模擬未考慮疫苗的防護(hù)效果不是100%，所以實(shí)際情況會(huì)比模擬結(jié)果略高，但是趨勢(shì)是一致的。研究表明，疫苗的保護(hù)率為45.15%～100%，免疫成功率達(dá)到80%以上[11]，且當(dāng)疫苗接種率達(dá)到70%時(shí)，甲型H1N1流感的發(fā)病率不足30%[12]，所以隨著疫苗接種率的增加，疫苗保護(hù)率的影響越低。

圖6 隔離+接種疫苗措施下甲型H1N1流感疫情模擬Fig.6 Simulation of influenza A (H1N1) epidemic under quarantine + vaccination measures

2.3 傳染病模型靈敏度分析

2.3.1 靈敏度相關(guān)參數(shù)分析

由于本次模型使用的參數(shù)a、p、k是通過(guò)查閱文獻(xiàn)資料和疫情數(shù)據(jù)綜合確定的，所以參數(shù)存在一定的不確定性。為了評(píng)估參數(shù)不確定性下輸入變量對(duì)模型失效概率的影響[13]，需對(duì)這些參數(shù)進(jìn)行靈敏度分析。結(jié)合文獻(xiàn)資料，參數(shù)的取值區(qū)間分別為a=0.3～1(潛伏期為1～3d)，p=0.2～0.3(病程為3～5d)，k=0.8～0.9，根據(jù)不同區(qū)間，在區(qū)間內(nèi)分別取a1=0.4、a2=0.6、a3=0.7；p1=0.2、p2=0.3；k1=0.8、k2=0.9，對(duì)數(shù)值進(jìn)行SEIQR模型的靈敏度分析。

2.3.2 靈敏度分析

運(yùn)用MATLAB軟件得出各個(gè)參數(shù)結(jié)果，通過(guò)分析結(jié)果可知，模型對(duì)參數(shù)a、p、k具有一定的敏感性，但模型的變化率均在合理范圍內(nèi)，模擬的曲線間的變化情況如圖7、8、9所示。

圖7 模型對(duì)參數(shù)a的靈敏度Fig.7 Sensitivity of the model to parameter a

圖8 模型對(duì)參數(shù)p的靈敏度Fig.8 Sensitivity of the model to parameter p

圖9 模型對(duì)參數(shù)k的靈敏度Fig.9 Sensitivity of the model to parameter k

3 結(jié)束語(yǔ)

綜合考慮在疫情實(shí)際傳播中可能遇到的情況建立多時(shí)段變參數(shù)的SEIQR模型，針對(duì)疫情傳播過(guò)程中采取的措施，將疫情傳播分為數(shù)個(gè)時(shí)段，根據(jù)每個(gè)時(shí)段內(nèi)的防控措施優(yōu)化模型參數(shù)，模擬疫情傳播中的變化情況，分析得出在接種疫苗的基礎(chǔ)上，即使沒(méi)有及時(shí)采取隔離防控措施，依然可以及時(shí)有效地控制住疫情發(fā)展，此研究結(jié)果對(duì)其他傳染病的模擬與預(yù)防具有一定的參考價(jià)值。由于模型參數(shù)的不確定性，運(yùn)用靈敏度分析方法，對(duì)不確定的參數(shù)進(jìn)行可靠性靈敏度分析，通過(guò)分析最后的模型圖像，評(píng)估輸入?yún)?shù)對(duì)模型的影響均在合理范圍內(nèi)。