基于改進SEIR模型的COVID-19疫情傳播建模

李 靜， 高 媛， 黃家妹

(1. 忻州師范學(xué)院 計算機系， 山西 忻州 034000； 2. 中北大學(xué) 大數(shù)據(jù)學(xué)院， 山西 太原 030051；3. 田納西大學(xué) 電子工程與計算機科學(xué)系， 田納西州 諾克斯維爾 37996)

0 引言

經(jīng)典的傳染病動力學(xué)模型SIR， SEIR等常被應(yīng)用于傳染病的傳播速度、 空間范圍、 傳播途徑和動力學(xué)機理等問題的研究， 其對染病人數(shù)的估計總體上符合傳染病的演變情況. 但是新型病毒的不斷出現(xiàn)， 人口流動性的增強， 各國醫(yī)學(xué)發(fā)展水平、 政府干預(yù)力度及疫苗研發(fā)進度等的不同， 都給疫情模擬增加了更多不確定的因素， 特別是新型冠狀病毒肺炎在全世界范圍的傳播， 使人們發(fā)現(xiàn)經(jīng)典的動力學(xué)模型已無法完成對新時期流行病傳播的模擬， 難以滿足人們對疫情發(fā)展預(yù)測的要求， 因而探索新的流行病傳播模型已成為生物信息學(xué)領(lǐng)域的研究熱點. 新冠疫情暴發(fā)以來， 研究者們對傳染病動力學(xué)模型進行了一些研究. 目前， 研究應(yīng)用較多的是SIR和SEIR倉室模型， 這些模型將傳染病流行范圍內(nèi)的人群劃分為幾個倉室， 通過求解微分方程計算各倉室人群可能的數(shù)量[1]. 對倉室模型的改進， 大多是在實際數(shù)據(jù)的驅(qū)動下對倉室人群的細分和新參數(shù)的引入這兩個方面[2-8]. 在考慮時變因素的研究方面， 朱仁杰等[9]以SIR模型為基礎(chǔ)， 通過添加新的參數(shù)表征感染系數(shù)隨時間的變化， 利用回歸分析進行參數(shù)估計， 實現(xiàn)了對7個國家疫情發(fā)展情況的模擬和預(yù)測； 喻孜等[10]采用時變的易感再生數(shù)、 當(dāng)日再生數(shù)、 潛伏再生數(shù)修正SIR模型， 研究了感染人數(shù)的變化趨勢； 王昕煒等[11]將時變的感染率和成功免疫率引入SEIR模型， 使用時滯參數(shù)描述潛伏期， 構(gòu)造了含有時滯的時變受控SEIR模型. 盡管這些模型都可以對新冠疫情的傳播進行模擬， 但由于實驗多完成于疫情暴發(fā)初期， 受數(shù)據(jù)采集等因素影響， 模型擬合效果受到一定限制.

綜上所述， 目前幾乎還沒有同時考慮潛伏期和時變參數(shù)的較為完善的流行病傳播模型. 因此， 本文在上述研究的基礎(chǔ)上， 以經(jīng)典的SEIR動力學(xué)模型為研究對象， 在引入潛伏者倉室的同時， 基于多項式回歸和多層感知機模型研究影響模型擬合效果的時變參數(shù)， 分析如何能更加真實地描述流行病傳播規(guī)律， 實現(xiàn)疫情發(fā)展趨勢的有效預(yù)測.

1 流行病傳播模型總體框架

針對傳統(tǒng)SEIR模型中參數(shù)恒定引發(fā)的瓶頸問題， 本文提出了基于時變參數(shù)的SEIR流行病傳播模型， 模型總體框架如圖 1 所示.

圖 1 流行病傳播模型總體框架

首先， 爬取疫情數(shù)據(jù)， 通過與權(quán)威發(fā)布的信息進行比對， 經(jīng)過必要的預(yù)處理操作， 確定最終實驗數(shù)據(jù)； 其次， 通過求解SEIR模型的微分方程得到模型參數(shù)的解析解， 對SEIR模型的時變參數(shù)進行假設(shè)優(yōu)化， 進而使用多項式回歸和多層感知機MLP模型進行參數(shù)擬合； 接著， 選取擬合效果較好的參數(shù)， 在SEIR模型中引入時變的感染系數(shù)和治愈系數(shù)， 形成改進的時變SEIR模型. 最后， 利用更新后的時變模型分別對各國新冠疫情數(shù)據(jù)進行建模.

2 改進的時變SEIR模型

2.1 SEIR模型[12-13]

(1)

(2)

(3)

(4)

通過微分方程構(gòu)建的傳染病動力學(xué)模型可以刻畫疾病的動力學(xué)行為， 而基本再生數(shù)R0是決定一種病毒是否流行的標志[14].R0表示在無病平衡狀態(tài)時， 引入一個新的染病者， 該染病者在其平均染病周期內(nèi)所能感染的人數(shù)的期望.當(dāng)R0< 1時， 染病者在平均染病期內(nèi)能傳染的人數(shù)小于1， 疾病會自行逐漸消亡； 當(dāng)R0> 1時， 疾病將持久傳播[15-16]. 對于SEIR模型， 患病者數(shù)量將出現(xiàn)增長， 到達頂峰后再單調(diào)減少最終趨于零.

2.2 參數(shù)優(yōu)化

在實際建模過程中， 兩個重要參數(shù)感染系數(shù)β和治愈系數(shù)γ通常依據(jù)經(jīng)驗設(shè)置為常數(shù)值.而在實際疫情發(fā)展過程中， 隨著各國政府不斷采取有效的防控措施， 易感人群的自我隔離和保護意識逐漸增強， 會使得感染系數(shù)β隨時間發(fā)生變化； 同時， 臨床經(jīng)驗的積累及疫苗的研發(fā)等又使得治愈系數(shù)γ不斷發(fā)生著變化.因此， 本文設(shè)置感染系數(shù)β和治愈系數(shù)γ為時間t的函數(shù)， 提出改進的時變SEIR模型.

假設(shè)易感者人數(shù)S(t)與人口總數(shù)N相等， 染病成為潛伏者的人群在經(jīng)過潛伏期后， 一般都會成為感染者，E(t)與I(t)僅存在時間上的延遲， 故忽略疾病潛伏期對感染系數(shù)β造成的影響， 假設(shè)I(t)=E(t)， 將式(2)化簡為

(5)

解微分方程， 得解析解為

E(t)=Ce(β-σ)t.

(6)

當(dāng)t=0時，C=E0， 即常數(shù)C為起始時刻疾病潛伏者人數(shù)， 將C代入式(5)得感染系數(shù)β為

(7)

同理， 得治愈系數(shù)γ為

(8)

2.3 參數(shù)擬合

分別采用多項式回歸及多層感知機MLP模型對時變參數(shù)進行擬合， 從而發(fā)現(xiàn)感染系數(shù)和治愈系數(shù)隨時間的變化規(guī)律， 這樣既能優(yōu)化流行病傳播的建模過程， 又能對疫情未來的發(fā)展趨勢進行更有效的預(yù)測.

1) 多項式回歸模型

任一函數(shù)都可以用多項式逼近， 根據(jù)已知參數(shù)與采集到的數(shù)據(jù)， 可以將模型中的時變參數(shù)擬合為時間t的函數(shù).通過增加自變量t的高次項對因變量進行逼近， 以二次多項式為例， 使用多項式回歸算法實現(xiàn)參數(shù)擬合的過程如下：

算法 1 基于多項式回歸模型的參數(shù)擬合

輸入： 各國的疫情數(shù)據(jù)集DataFrame.

輸出： 3個參數(shù)擬合結(jié)果(形如二次項系數(shù)A， 一次項系數(shù)B， 常數(shù)項C).

a) 數(shù)據(jù)預(yù)處理.將數(shù)據(jù)集劃分為X和Y兩個組件，X為時間t的數(shù)組，Y為依據(jù)式(7)和式(8)計算所得參數(shù)β和γ的值.

b) 使用多項式特征擴展器PolynomialFeatures增加自變量的高次項， 構(gòu)造X的多項式特征X_Poly，X_Poly=[x0，x1，x2].

c) 使用線性回歸器LinearRegression學(xué)習(xí)X_Poly 和Y之間的映射關(guān)系.

d) 獲得各多項式特征的參數(shù)， 即待擬合的各參數(shù).

2) 多層感知機MLP模型[17-19]

多層感知機模型是一種前饋神經(jīng)網(wǎng)絡(luò)， 可以映射一組輸入向量到一組輸出向量. 它由輸入層、 隱藏層和輸出層組成， 不同層之間的節(jié)點間采用全連接的方式. MLP的網(wǎng)絡(luò)模型結(jié)構(gòu)如圖 2 所示.

圖 2 MLP網(wǎng)絡(luò)模型結(jié)構(gòu)

使用MLP進行參數(shù)擬合的過程如下：

算法 2 基于MLP的參數(shù)擬合

輸入： 各國的疫情數(shù)據(jù)集DataFrame.

輸出： MLP的權(quán)重矩陣及偏置項(形如mlp.coefs_和mlp.intercepts_).

a) 數(shù)據(jù)預(yù)處理. 將數(shù)據(jù)集劃分為X和Y兩個組件，X為時間t的數(shù)組，Y為依據(jù)式(7)和式(8)計算所得參數(shù)β和γ的值.

b) 隨機初始化權(quán)重矩陣wij,以X作為輸入，Y作為輸出訓(xùn)練網(wǎng)絡(luò).

c) 前向傳播過程. 當(dāng)前層神經(jīng)元與對應(yīng)的權(quán)重值相乘并累加， 再加上偏置值得到zj， 接著經(jīng)過一個激活函數(shù)f得到當(dāng)前層的輸出yj， 即下一層的輸入. 信息從輸入層向前傳播， 直到得到輸出值， 如式(9)和式(10)所示.

(9)

(10)

d) 反向傳播過程. 定義損失函數(shù)Loss， 計算模型總誤差， 利用梯度下降法迭代求解損失函數(shù)的極小值， 實現(xiàn)殘差的逐層反饋， 進而不斷更新權(quán)重值， 權(quán)值更新公式如式(11)所示.

(11)

式中：η為模型的學(xué)習(xí)率.

e) 重復(fù)以上前向傳播與反向傳播過程， 直到模型收斂， 可得到模型的權(quán)重矩陣及偏置項.

經(jīng)過以上訓(xùn)練所得到的MLP模型， 被認為可以以任意精度擬合任意復(fù)雜度的函數(shù)[20].

本文將多項式回歸與多層感知機模型運用于SEIR模型中時變參數(shù)的估計， 以優(yōu)化感染系數(shù)β及治愈系數(shù)γ， 分別得到Poly_Time_SEIR及MLP_Time_SEIR模型.

3 疫情傳播建模流程

基于兩種改進的時變SEIR模型， 根據(jù)采集到的疫情數(shù)據(jù)和已有參數(shù)， 結(jié)合改進模型所得的時變參數(shù)， 實現(xiàn)對疫情的傳播建模. 主要包含以下幾個步驟：

1) 對爬取到的疫情數(shù)據(jù)進行整理、 清洗， 存儲于CSV文件中， 用于時變參數(shù)的計算以及疫情傳播模型的建立.

2) 對文件中的疫情數(shù)據(jù)分別使用多項式回歸及MLP模型計算其感染系數(shù)β及治愈系數(shù)γ， 將得到模型參數(shù)存儲于模型文件PKL中.

3) 在疫情傳播建模階段， 首先建立傳統(tǒng)SEIR微分方程； 然后讀取步驟2)模型文件中的模型參數(shù)(多項式回歸)或權(quán)重值矩陣(MLP)， 合成時變參數(shù)β及γ； 接著用求得的時變參數(shù)替換SEIR模型原參數(shù)； 最終分別得到改進的Poly_Time_SEIR與MLP_Time_SEIR模型.

4) 分別求解改進的Poly_Time_SEIR與MLP_Time_SEIR模型建立的微分方程， 得到S，E，I，R各倉室的人群數(shù)量； 比較兩種模型， 將誤差較小且擬合效果較好的模型作為最終的疫情傳播模型.

4 實驗與結(jié)果分析

4.1 實驗環(huán)境與數(shù)據(jù)集

采用Python 3.7作為實驗平臺， 操作系統(tǒng)為Windows 10， CPU為CoreTM i7-9750H. 實驗數(shù)據(jù)集爬取自公信度較高的網(wǎng)易官方網(wǎng)站肺炎疫情實時動態(tài)播報平臺， 并與中華人民共和國國家衛(wèi)生健康委員會官網(wǎng)及美國田納西大學(xué)電子工程與計算機科學(xué)系實驗室采集到的數(shù)據(jù)進行比對修正. 實驗共收集2020年1月21日至2021年6月31日各國的疫情數(shù)據(jù)， 分為歷史數(shù)據(jù)和實時數(shù)據(jù)兩部分. 其中， 疫情歷史數(shù)據(jù)包括累計確診人數(shù)(total_confirm)， 累計治愈人數(shù)(total_heal)和累計死亡人數(shù)(total_dead)； 實時數(shù)據(jù)包括當(dāng)日新增確診人數(shù)(today_confirm)， 新增治愈人數(shù)(today_heal)和新增死亡人數(shù)(today_dead).

4.2 參數(shù)擬合

4.2.1 基于多項式回歸模型的參數(shù)擬合

基于采集到的各國340 d左右的疫情數(shù)據(jù)， 以形如At2+Bt+C的二次多項式為例， 利用多項式回歸模型對參數(shù)A，B和C進行擬合.式(7)和式(8)中，E(t)為t時刻疾病潛伏者的人數(shù)， 新冠肺炎潛伏期分布滿足偏態(tài)分布， 鐘南山院士研究團隊在文獻[21]中表示， 潛伏期中位數(shù)為4 d， 故采用t+4時刻的確診人數(shù)表示t時刻的潛伏者人數(shù).E0為開始時刻的潛伏者人數(shù)， 潛伏系數(shù)σ取潛伏天數(shù)的倒數(shù)， 即0.25.I(t) 為t時刻感染者的人數(shù)，I0為開始時刻的感染者人數(shù).以中國和美國為例， 經(jīng)模型擬合， 得到的感染系數(shù)和治愈系數(shù)中各參數(shù)值分別如表 1 和表 2 所示.

表 1 感染系數(shù)中各參數(shù)值

表 2 治愈系數(shù)中各參數(shù)值

可以看出： 感染系數(shù)中的A1值均大于0， 表示兩國的感染系數(shù)在初期都呈下降趨勢； 而治愈系數(shù)中的A2值均小于0， 表示治愈系數(shù)在開始階段都呈上升趨勢.由決定系數(shù)R2可以看出， 多項式回歸模型對兩國感染系數(shù)和治愈系數(shù)的擬合效果較好.為了更直觀地看到時變參數(shù)隨時間的變化趨勢， 根據(jù)擬合參數(shù)計算兩國的感染系數(shù)和治愈系數(shù)， 并與真實系數(shù)進行實驗對比， 比較結(jié)果如圖 3 所示.

(a) 感染系數(shù)β隨時間的變化

從圖 3 可以看出， 隨著時間的推移， 兩國感染系數(shù)的擬合曲線整體呈下降趨勢， 中國的感染系數(shù)始終低于美國. 而兩國治愈系數(shù)的擬合曲線整體呈上升趨勢， 中國的治愈系數(shù)高于美國， 模型的擬合結(jié)果與實際情況基本相符. 但就極值數(shù)據(jù)而言， 中國的實際感染系數(shù)約在第112天左右達到最低值， 但是， 擬合曲線的最低值出現(xiàn)于第188天附近； 美國感染系數(shù)約在第29天出現(xiàn)峰值， 而擬合曲線并沒有出現(xiàn)感染系數(shù)上升的過程. 同樣， 模型對治愈系數(shù)極值的擬合也存在偏差， 且中國的實際治愈系數(shù)在第138天左右達到峰值， 之后基本保持不變， 而擬合曲線在達到峰值后， 出現(xiàn)了明顯的下降.

4.2.2 基于多層感知機MLP模型的參數(shù)擬合

經(jīng)多次實驗對比， 選用含2層隱藏層的MLP來擬合感染系數(shù)， 含1層隱藏層的MLP來擬合治愈系數(shù). 估計感染系數(shù)β時， 第一隱層設(shè)置110個神經(jīng)元， 第二隱層設(shè)置100個神經(jīng)元， 經(jīng)過約 4 000次訓(xùn)練， 損失函數(shù)達到收斂. 估計治愈系數(shù)γ時， 使用含110個神經(jīng)元的單隱層神經(jīng)網(wǎng)絡(luò)， 約2 000次訓(xùn)練達到收斂. 基于MLP訓(xùn)練所得的權(quán)重矩陣來計算β和γ值， 得到兩國時變參數(shù)變化趨勢如圖 4 所示.

(a) 感染系數(shù)β隨時間的變化

從圖 4 可以看出， 相比基于多項式回歸模型的參數(shù)擬合， 基于MLP模型擬合的參數(shù)不僅整體變化趨勢與實際相符， 而且更好地捕捉到了極值數(shù)據(jù). 由于中國衛(wèi)健委于2020年2月12日公布的確診病例中包含了臨床診斷病例， 導(dǎo)致實際感染系數(shù)與之前相差較大， 出現(xiàn)了明顯的斷層， 但并沒有影響模型對之后感染系數(shù)的擬合. 擬合曲線在達到極小值后， 感染系數(shù)基本保持不變. 同樣， 中國治愈系數(shù)的擬合曲線在達到峰值后， 也沒有出現(xiàn)明顯下降. 可見， MLP模型對時變參數(shù)的擬合與實際數(shù)據(jù)更吻合.

4.3 基于最優(yōu)參數(shù)的疫情模擬

4.3.1 疫情傳播建模

在中國2020年1月21日至12月31日的疫情數(shù)據(jù)集上， 使用經(jīng)典的SIR， SEIR， 文獻[9]中提出的基于線性回歸感染系數(shù)的SIR模型(文中記作Linear_Time_SIR)及本文提出的基于多項式回歸時變參數(shù)的Poly_Time_SEIR模型及基于多層感知機MLP時變參數(shù)的MLP_Time_SEIR模型， 對各模型計算所得的感染人數(shù)I及治愈人數(shù)R與實際人數(shù)做了實驗比較， 實驗結(jié)果如圖 5 所示. 圖中， 實線代表感染人數(shù)I， 虛線代表治愈人數(shù)R.

圖 5 各模型與真實數(shù)據(jù)的對比

從圖 5 可以看出： 經(jīng)典的SIR模型由于忽略了疾病的潛伏者人群， 對感染、 治愈兩類人群的估計都偏低， 在103數(shù)量級； SEIR模型及本文提出的MLP_Time_SEIR模型估計值都在104數(shù)量級， 但由于SEIR模型采用恒定的感染系數(shù)與治愈系數(shù)， 對兩類人群的估計值均低于實際值； 雖然Linear_Time_SIR， Poly_Time_SEIR及MLP_Time_SEIR 3種模型都是基于時變參數(shù)的， 但前兩種模型在中國疫情數(shù)據(jù)集上的擬合效果并不理想， 對兩類人群的擬合誤差分別在103和104數(shù)量級上， 而本文提出的MLP_Time_SEIR模型在引入潛伏者人群的SEIR模型基礎(chǔ)上， 由于使用前饋神經(jīng)網(wǎng)絡(luò)對時變參數(shù)進行了較好的擬合， 使其對感染和治愈兩類人群的估計誤差最小， 實驗取得了較好的效果.

4.3.2 基本再生數(shù)R0的變化趨勢

傳染病模型中基本再生數(shù)R0的導(dǎo)出方法有多種， 本文采用定義法將模型的基本再生數(shù)定義為感染系數(shù)β與治愈系數(shù)γ的比值， 得到2020年1月24日至2021年12月29日共338 d兩種模型在中國數(shù)據(jù)集上的基本再生數(shù)隨時間的演化趨勢， 如圖 6 所示. 圖中， 實線為模型的時變基本再生數(shù)R0， 空心點表示基本再生數(shù)的平均值.

(a) MLP_Time_SEIR模型的R0

由圖 6 知： MLP_Time_SEIR模型R0的均值為1.69， 方差為9.63， 波動較大； Poly_Time_SEIR模型R0的均值為 1.49， 方差為1.31. 兩種模型導(dǎo)出的基本再生數(shù)均小于鐘南山院士早期對疫情基本再生數(shù)的估計(2～3之間)[21]， 這得益于疫情發(fā)展過程中政府與民眾的積極防疫. 同時，R0值的降低也意味著會延緩疫情峰值的到來， 緩解醫(yī)療系統(tǒng)所面臨的壓力.R0均值大于1說明疫情不會自行消亡， 疾病的有效預(yù)防、 控制和免疫策略的制定仍然具有重要意義.

4.4 模型測試

4.4.1 中美疫情測試

選取中國和美國2020年1月21日至11月30日的數(shù)據(jù)作為訓(xùn)練集， 12月1日至12月29日的數(shù)據(jù)作為測試集， 對各模型所得感染人數(shù)、 治愈人數(shù)與真實數(shù)據(jù)的均方根誤差RMSE進行對比實驗， 結(jié)果如表 3 所示.

表 3 各種模型的均方根誤差(RMSE)比較

由表 3 可以看出： 本文提出的MPL_Time_SEIR模型在中國疫情數(shù)據(jù)集上對感染人數(shù)和治愈人數(shù)的預(yù)測誤差都是最低的， 模型的回歸效果相比真實值平均只相差6.82和25.85； 在美國疫情數(shù)據(jù)集上， 由于美國疫情相比中國要嚴重得多， 感染人數(shù)和治愈人數(shù)都遠遠高于中國， 各模型的預(yù)測誤差較大， MPL_Time_SEIR對感染人數(shù)的估計誤差相對較小， 比經(jīng)典的SEIR模型均方根誤差下降了42.9%； Poly_Time_SEIR模型對治愈人數(shù)的估計誤差最小， MPL_Time_SEIR僅次于Poly_Time_SEIR模型.

經(jīng)MLP訓(xùn)練時變參數(shù)， 基于擬合得到的最優(yōu)參數(shù)值， 用改進的MLP_Time_SEIR模型對中國2020年12月1日～12月29日共29 d內(nèi)的感染人數(shù)和治愈人數(shù)進行預(yù)測， 預(yù)測結(jié)果如圖 7 所示. 圖中， 實線表示中國2020年12月感染人數(shù)和治愈人數(shù)的真實值， 標記點表示MLP_Time_SEIR模型預(yù)測到的感染人數(shù)和治愈人數(shù).

(a) 感染人數(shù)預(yù)測

由圖 7 可以看出， 中國感染人數(shù)和治愈人數(shù)均呈現(xiàn)緩慢上升趨勢， 預(yù)測標記點很好地擬合了真實曲線， 模型對中國疫情傳播的預(yù)測達到了很好的效果.

4.4.2 意大利疫情測試

在意大利數(shù)據(jù)集上采用MPL_Time_SEIR模型對2021年5月和6月的感染人數(shù)及治愈人數(shù)進行預(yù)測， 并與真實數(shù)據(jù)進行對比， 比較結(jié)果如圖 8 所示.

圖 8 意大利2021年5月～6月疫情趨勢預(yù)測

從圖 8 可以看出， 模型對感染人數(shù)的擬合效果較好， 治愈人數(shù)與真實曲線的整體趨勢一致， 但高于實際治愈人數(shù)， 均方根誤差為59 745.9， 仍然優(yōu)于其他模型.

5 結(jié) 論

本文在經(jīng)典SEIR模型的基礎(chǔ)上， 提出了基于改進模型的疫情傳播建模方法. 針對經(jīng)典模型的參數(shù)恒定問題， 進行了參數(shù)優(yōu)化， 給出了基于多項式回歸與多層感知機模型的參數(shù)擬合算法， 并選取最優(yōu)參數(shù)對疫情傳播過程進行了建模. 實驗結(jié)果表明：

1) 在參數(shù)擬合效果方面， 兩種模型Poly_Time_SEIR 和MPL_Time_SEIR對時變參數(shù)的擬合整體與實際情況相符， 但MPL_Time_SEIR模型可以更好地捕捉到極值點.

2) 在疫情趨勢預(yù)測方面， MPL_Time_SEIR模型對感染人數(shù)和治愈人數(shù)的預(yù)測相對于其他4種模型的預(yù)測誤差更低， Poly_Time_SEIR模型對美國12月治愈人數(shù)預(yù)測的效果相對更好， 體現(xiàn)了兩種模型的有效性.