基于反應擴散方程研究新型冠狀病毒肺炎在武漢早期傳播特征

郭尊光, 李明濤, 常利利, 張 娟, 梁 娟, 邢國榮, 張 偉, 孫桂全3,,?

(1- 中北大學大數(shù)據(jù)學院，太原 030051; 2- 太原工業(yè)學院理學系，太原 030008;3- 中北大學理學院，太原 030051; 4- 太原理工大學數(shù)學學院，太原 030024;5- 山西大學復雜系統(tǒng)研究所，太原 030006;6- 疾病防控的數(shù)學技術與大數(shù)據(jù)分析山西省重點實驗室，太原 030006)

1 引言

傳染病自古以來就是危害人類健康的主要敵人之一，預防和控制傳染病仍然是當今世界的重要課題.冠狀病毒是一類具有包膜、基因組為線性單股正鏈的RNA 病毒，這種病毒會引起疾病，患者表現(xiàn)為從普通感冒到重度肺部感染等不同臨床癥狀.根據(jù)武漢市衛(wèi)健委通報的武漢肺炎情況，新型冠狀病毒是一類β 屬的冠狀病毒，有包膜、呈圓形或橢圓形顆粒，直徑60–140 nm.感染新型冠狀病毒主要表現(xiàn)發(fā)熱、干咳、乏力.潛伏期為1–14 天，多為一周左右.重癥患者出現(xiàn)呼吸困難，嚴重者表現(xiàn)為急性呼吸窘迫綜合征、膿毒癥休克及多器官功能衰竭.新型冠狀病毒肺炎已被納入《中華人民共和國傳染病防治法》規(guī)定的乙類傳染病.為防止新型冠狀病毒肺炎疫情傳播和擴散，最大限度的保障人民群眾的健康，為加強全國面上的疫情防控工作.先后有30 個省啟動了突發(fā)公共衛(wèi)生事件的一級響應.世衛(wèi)組織總干事2020 年3 月13 日在2019 冠狀病毒病(COVID-19)疫情媒體通報會上的講話中提及，全球確診病例已超過13 萬人，死亡病例達5000 多例，波及120 多個國家，已形成全球大流行的事實.

為研究新型冠狀病毒(COVID-19)的傳播規(guī)律，全世界的醫(yī)務工作者和生物數(shù)學的科研工作者已投入其中.在武漢發(fā)生新型冠狀病毒肺炎初期，Lin 等人[1]分析了最初的425 例確診病例的臨床資料，確定了新冠肺炎的一些流行病學特征.Chinazzi 等人[2]研究了旅游限制對新型冠狀病毒爆發(fā)傳播的影響，使用全球集合種群疾病傳播模型來預測旅行限制對該流行病的國內(nèi)和國際傳播的影響.Lin 等人[3]考慮了個人的行為反應和政府的行動，提出了一種SEIRNDC 概念模型，研究了疫情的發(fā)展趨勢.黃森忠等人[4]運用SEIR 模型對新冠肺炎疫情進行了預測并對控制策略的效率進行了評估.唐三一等人[5]提出了新型的七倉室的離散隨機模型，分析了不同的參數(shù)對疾病二次爆發(fā)風險的影響并給出了相應策略.嚴閱等人[6]提出了基于時滯微分方程的傳染病模型并模擬了疫情的趨勢.王霞等人[7]構(gòu)建了復雜網(wǎng)絡模型對武漢周邊地區(qū)何時復工進行了研究.Yang 等人[8]提出了一種修正的SEIR 模型，得到中國的疫情將在2 月底達到頂峰，4 月底呈現(xiàn)逐漸下降的趨勢，如果推遲5 天實施封城中國大陸的疫情規(guī)模將會擴大3 倍.還有一些醫(yī)務工作者進行了臨床研究[9–14].然而，現(xiàn)有的研究都是基于常微分方程、時滯常微分方程或是統(tǒng)計方法和臨床病例研究，沒有考慮到個體的空間擴散.事實上，空間擴散模型更能反應新型冠狀病毒在時間和空間上的傳播動力學規(guī)律.因此，本文主要是基于空間擴散的新型冠狀病毒肺炎模型研究武漢的早期傳播情況.

本文主要研究了武漢市2020 年1 月5 日之前的新型冠狀病毒傳染病的傳播動力學，本文的傳染病數(shù)據(jù)是2019 年12 月8 日至2020 年1 月5 日武漢市首批105 例確診病例(數(shù)據(jù)來源于文獻[1]).在傳播早期，國家和武漢市還沒有采取嚴格措施，可以認為人在空間上是自由擴散的.基于此我們建立了含有空間擴散的傳染病模型，結(jié)合現(xiàn)有的文獻及各種新聞報道，基于早期染病者數(shù)據(jù)和最小二乘法，我們對反應擴散方程進行了參數(shù)估計，找到了模型的最優(yōu)參數(shù)，分析了染病者數(shù)量和潛伏者的擴散系數(shù)的敏感性.

本文的結(jié)構(gòu)如下：第2 部分基于空間擴散建立了反應擴散傳染病模型，并對模型的參數(shù)進行了解釋，第3 部分針對模型我們給出了主要結(jié)果，包含參數(shù)估計和參數(shù)的敏感性分析，第4 部分給出了結(jié)論.

2 模型

在新型冠狀病毒傳染病傳播的初期，不考慮人口動力學，即沒有自然出生和死亡.目前所見傳染源主要是新型冠狀病毒感染的患者，且潛伏期具有傳染性，傳染途徑主要是接觸傳播.根據(jù)流行病學特征，我們將研究人群分為：易感者S、潛伏者E、染病者I、移出者R 四類.這四類人群的轉(zhuǎn)移關系可用倉室圖1 表示.

圖1: 易感者、潛伏者、染病者及移出者倉室轉(zhuǎn)移圖

考慮到空間各倉室個體的空間擴散，我們建立如下傳染病模型

模型(1)各符號解釋，見表1.

表1: 模型符號說明

在新冠肺炎傳播的過程中，易感者與潛伏者或染病者接觸會有一定概率的被傳染.被傳染的初期會有一定的潛伏期，經(jīng)過潛伏期后轉(zhuǎn)化為染病者，經(jīng)過治療或自我免疫會有一定比例的移出.接下來我們結(jié)合數(shù)據(jù)對模型進行研究.

3 主要結(jié)果

3.1 參數(shù)估計

本節(jié)將對模型(1)進行參數(shù)估計，方法基于最小二乘估計.具體過程是先將模型離散化，對時間的導數(shù)使用向前差分，空間二階導數(shù)使用中心差分格式，得到模型離散格式

模擬時，我們將武漢市假設成一個方形區(qū)域[0,102]×[0,102]，選取網(wǎng)格剖分

xi=ih, i=0,1,2,··· ,M, Mh=102,

yj=jh, j =0,1,2,··· ,M, Mh=102,

tn=100nτ, n=0,1,2,··· ,N.

研究武漢市早期新冠肺炎情況，可以考慮齊次紐曼邊界條件(即零流邊界)，邊界的離散格式為：左邊界

S0,j=S2,j, E0,j=E2,j, I0,j=I2,j, R0,j=R2,j, 0 ≤j ≤M.

右邊界

SM,j=SM?2,j, EM,j=EM?2,j, IM,j=IM?2,j, RM,j=RM?2,j, 0 ≤j ≤M.

上邊界

Si,M=Si,M?2, Ei,M=Ei,M?2, Ii,M=Ii,M?2, Ri,M=Ri,M?2, 0 ≤i ≤M.

下邊界

Si,0=Si,2, Ei,0=Ei,2, Ii,0=Ii,2, Ri,0=Ri,2, 0 ≤i ≤M.

模擬時將1500 萬人口均勻分布在方形區(qū)域除邊界的內(nèi)部格點上，共有內(nèi)部格點數(shù)為101×101 = 10201，所以初始時刻每個格點上易感者數(shù)量約為1470，基于目前的流行病學調(diào)查，潛伏期為1–14 天，多為3–7 天(出自中國疾控中心文件《新型冠狀病毒肺炎公眾防護指南(第2 版)》[15]).我們?nèi)摲谥形粩?shù)(5 天)得到δ = 0.2，染病者的移出率[16]γ = 0.05，由于新型冠狀病毒肺炎初期人們對疾病了解甚少，潛伏者可能表現(xiàn)輕微癥狀使其活動能力較正常人低，會引起潛伏者的擴散速度低于易感者，染病者以發(fā)熱、干咳、乏力為主要表現(xiàn)，一部分患者同時伴有鼻塞、流涕、咽痛和腹瀉等癥狀，會導致染病者的擴散速度小于潛伏者的擴散速度，參考移出者的實際情況，我們假設這四類人的擴散速率有關系dS> dE> dR> dI.模擬時分別取dS= 10, dE= 9, dR= 1.5, dI= 1，假設潛伏者的傳染率與染病者的傳染率相等，在首個患者被確診時假設已有潛伏者數(shù)量E =5.模擬時時間步長和空間步長分別取為τ =0.01 和h=1.

記

圖2: 單點爆發(fā)下局部模型的染病者和真實統(tǒng)計數(shù)據(jù)隨時間的變化情況

圖3: 潛伏者隨時間在空間中的分布情況

圖4: 染病者隨時間在空間中的分布情況

3.2 參數(shù)敏感性分析

為了更好的控制傳染病的傳播，需要研究模型中各個參數(shù)對染病者數(shù)量影響的大小.若其余參數(shù)固定的情況下改變某一個參數(shù)的值引起染病者數(shù)量的相對改變量較大，此時稱染病者數(shù)量對該參數(shù)敏感性較大，否則，稱敏感性較低.本節(jié)我們將研究染病者累計數(shù)量對對潛伏者擴散系數(shù)dE、潛伏者的傳染率β1、染病者的傳染率β2的敏感性.由于這三個系數(shù)耦合在模型中，所以我們可以通過改變參數(shù)的數(shù)值，對模型進行數(shù)值求解，通過得到的整個空間上的染病者累計數(shù)量來研究敏感性的大小.敏感性指數(shù)的大小可以對疾病的防控措施提供一等的理論支持，對于敏感性較大的參數(shù)，疾病的防控要特別重視，敏感性的大小可以通過圖形反應出來，下面我們將對潛伏者的擴散率、潛伏者和染病者的傳染率進行敏感性分析，并根據(jù)敏感性的大小提出相應的預防和控制措施.

3.2.1 針對潛伏者的擴散率dE 的敏感性分析

任何一種人際傳播的傳染病，人口的自由擴散速度會影響傳染病的傳播.傳染病防控重點是管理好傳染源，切斷傳染途徑，針對新型冠狀病毒肺炎傳染病，由于潛伏者也可以傳播，所以傳染源包括染病者和潛伏者.由于潛伏者前期無癥狀或是癥狀不明顯，所以潛伏者不容易被發(fā)現(xiàn)，但潛伏者的擴散會引起傳染病的傳播，防控潛伏者的傳播對疫情的防空非常重要.

我們固定參數(shù)值dS=10, dI=1, dR=1.5, β1=0.1316, β2=0.1316, δ =0.2, γ =0.05，潛伏者的擴散率dE分別取dE= 0, dE= 1 和dE= 10，得到2019 年12 月8 日到2010 年1 月5 日的敏感性.我們將結(jié)果展示在圖5 中，為了更易看清楚染病者數(shù)量隨潛伏者擴散系數(shù)的敏感性，將圖5 中(a)的部分局部放大得(b)圖.這幅圖中清楚的看到，隨著潛伏者擴散系數(shù)的增加染病者數(shù)量也在增加，再次說明控制潛伏者的空間移動對傳染病的防控至關重要，若傳染病的早期及時控制潛伏者的空間擴散可以有效降低染病者的數(shù)量.為居家隔離等措施提供了理論依據(jù).

圖5: 染病者數(shù)量對潛伏者的擴散系數(shù)的敏感性

3.2.2 針對傳染率β1 的敏感性分析

我們固定參數(shù)值dS= 10, dE= 9, dI= 1, dR= 1.5, β2= 0.1316, δ = 0.2, γ =0.05，潛伏者的傳染率β1分別取0.07, 0.1 和0.13，得到2019 年12 月8 日到2010 年1 月5日的敏感性.結(jié)果展示在圖6 和圖7 中.從圖中可以看出染病者數(shù)量隨潛伏者的傳染率增大而增加，所以為防止新冠肺炎的傳播應控制潛伏者的傳染率，盡量減少與潛伏者的接觸.所以早期居家隔離，減少去公共場所，減少走親訪友都是有效降低接觸率的有效措施.

3.2.3 針對傳染率β2 的敏感性分析

我們固定參數(shù)值dS= 10, dE= 9, dI= 1, dR= 1.5, β1= 0.1316, δ =0.2, γ = 0.05，染 病 者 的 傳 染 率β2分 別 取0.12, 0.15 和0.18，得 到2019 年12 月8 日到2010 年1 月5 日的敏感性，結(jié)果展示在圖8 和圖9 中.從圖中同樣可以看出，染病者的數(shù)量隨染病者傳染率的增大而顯著增加.在這個意義上而言，如何降低染病者的傳染率是傳染病防控的重要課題.比如外出佩戴口罩，減少無必須去醫(yī)院的次數(shù)，對染病者及時進行集中隔離，對染病者接觸人員進行有效追蹤及時隔離等措施都可以有效控制染病者的傳染率.

圖6: 染病者數(shù)量對潛伏者傳染率β1 的敏感性

圖8: 染病者數(shù)量對潛伏者傳染率β2 的敏感性

圖9: 染病者數(shù)量對潛伏者傳染率β2 的敏感性(2020 年1 月3 日的空間分布)，從這三幅圖中可以看出隨β2 的增加，每個固定位置的染病者數(shù)量同樣會升高，并且會有更多的空間位置出現(xiàn)染病者，意味著新冠肺炎會隨傳染率的升高而向空間四面八方擴散

4 結(jié)論

本文研究了新型冠狀病毒(COVID-19)肺炎在武漢早期傳播情況.在武漢封城之前，所有個體在空間上都是自由的，可以認為個體是隨機游走，在空間上存在著擴散行為.基于此，我們建立了含有空間擴散的反應擴散方程傳染病模型，在假設潛伏者和染病者的傳染率相同的情況下，通過參數(shù)估計尋找到了最優(yōu)的傳染率.對模型采用有限差分格式進行了數(shù)值模擬，染病者和潛伏者隨時間在空間中的傳播情況我們通過三維圖形進行了呈現(xiàn).并對潛伏者的擴散系數(shù)、潛伏者的傳染率和染病者的傳染率進行了敏感性分析，結(jié)果顯示對擴散系數(shù)敏感性偏小，對傳染率的敏感性較大，所以在疫情防控方面要不僅要控制患者的擴散，而且需要重視如何降低潛伏者和染病者的傳染率.因此，為了有效防控新型冠狀病毒肺炎，應該盡量不去傳染病爆發(fā)的區(qū)域，避免直接與患者的接觸，做好充分的保護措施，減少外出暴露機會.對于突發(fā)傳染病，早發(fā)現(xiàn)早隔離，及時切斷傳染源，有效降低傳染率是防控防治傳染的重要舉措.

本文聚焦于研究傳染病的早期傳播特征，即沒有干預措施或者干預措施不明顯的時間階段.而對于后期防控和治療的介入，所構(gòu)建的反應擴散模型是需要修正的.比如，居家隔離使得人群在空間無法形成擴散，進而使得人群以不同的社區(qū)尺度被隔離，需要研究傳染病在社區(qū)是如何形成傳染的.在這種情況下，基于復雜網(wǎng)絡的集合種群模型更適合去描述傳染病的傳播.