安徽省間日瘧傳染率分析

白 迪，周武略，趙繼軍

(青島大學復雜性科學研究所，山東 青島 266071)

0 引言

瘧疾(malaria)是一種由于按蚊叮咬或輸入帶瘧原蟲者的血液而感染瘧原蟲所引起的蟲媒傳染病[1]，寄生于人體的瘧原蟲共有4種，即間日瘧原蟲，三日瘧原蟲，惡性瘧原蟲和卵形瘧原蟲。其中間日瘧原蟲傳播范圍最廣，惡性瘧原蟲致死率最高[2]。自20世紀50年代開始中國對瘧疾逐步開展了有計劃的全面防治，到20世紀末大部分地區的瘧疾流行已得到控制[3]。雖然一些地區感染瘧疾的病例維持在了較低的水平，甚至部分省份自2019年開始陸續通過了省級消除瘧疾評估，但是這些地區的瘧疾傳播媒介按蚊仍普遍存在，瘧疾傳播的條件尚未消除[4]，仍有再次流行的風險。江蘇省泗洪縣在20世紀90年代瘧疾發病數降低至歷史最低水平后，于2000年再次暴發了疫情，這次暴發與從鄰近省份輸入的瘧疾病例有關[5]。由于瘧疾傳播媒介的存在，輸入病例易成為傳染源，引起繼發傳播。此外對非洲高地瘧疾的研究中發現，瘧疾的傳染率受降水量周期性變化的影響[6]。當降水量的年周期降低時，瘧疾的傳播程度逐漸增大。當滿足一些特定條件的情況下，例如足夠多的輸入病例、未來某段時間氣象因素的變化都有可能導致瘧疾的重新暴發。

對瘧疾在中國流行的研究中發現，瘧疾的報告發病數與溫度、降水量等因素有關。在氣溫接近26℃與27℃時瘧疾的發病數相對較高[7-8]；瘧疾發病數的時間序列與滯后0～3個月的降水量呈現強正相關，與滯后4～9個月的降水量呈負相關[9]。這些研究是對瘧疾發病數與氣候因素之間的相關性研究。中國本地出現的瘧疾類型主要是間日瘧，其傳播媒介為中華按蚊[10]。相對于非洲高地關于岡比亞按蚊傳播惡性瘧傳染率的研究，中國通過中華按蚊傳播間日瘧傳染率的研究極其有限。傳染率是疾病傳播特性中的一個重要參數，決定了疾病流行的程度(即發病數的高低)和流行的周期。本文以安徽省2004～2009年流行的間日瘧為例，分析間日瘧的傳播動態，包括傳染率和傳染率季節性。

安徽省在地理位置上處于中國東南部，為亞熱帶和暖溫帶過渡氣候，四季分明，適宜間日瘧傳播媒介中華按蚊的繁殖，曾多次暴發間日瘧疫情[11]。2005～2015年間，安徽省間日瘧累計發病數達59 752例，居全國首位，占全國累計間日瘧發病數的51.39%[12]。本文將在研究安徽省間日瘧傳播動態的基礎上，分析氣候因素包括溫度、降水量對通過中華按蚊傳播的間日瘧傳染率的影響。本文的研究結果為進一步預測間日瘧在一個地區消除后重新暴發的可能性以及影響暴發的條件等方面的研究提供指導。

1 數據

本文從中國疾病預防控制中心(www.phsciencedata.cn)獲得安徽省2004～2009年各月間日瘧的發病數。從安徽省統計局(http://tjj.ah.gov.cn/tjjweb/web/index.jsp)獲得2004～2009年安徽省人口出生率、死亡率。從中國氣象數據網(http://data.cma.cn/user/toLogin.html)獲得安徽省各地級市氣象站點2004～2009年溫度、降水量的數據。

2 方法

建立Susceptible Infected Recovered-Susceptible Infected(SIR-SI)模型估算間日瘧傳染率矩陣并分析傳染率的季節性，進一步通過多元線性回歸模型分析影響間日瘧傳染率的因素。

2.1 SIR-SI模型估算傳染率

圖1 SIR-SI模型示意圖

為研究安徽省間日瘧的傳播特性，我們建立病毒通過蚊子叮咬在人與蚊子之間傳播間日瘧的SIR-SI倉室模型。倉室模型是由微分方程構成的確定性建模方法。對于蚊媒傳播的疾病來說，倉室模型將人群和蚊群分為幾種不同類型，每種類型作為一個“倉室”存在。本文中的模型包括兩個部分，SIR部分描述人類易感者(Sh)、感染者(Ih)、恢復者(Rh)各倉室中人群數量的動態，SI部分描述蚊子易感(Sm)、感染(Im)各倉室中蚊子數量的動態。由于蚊子的生命周期較短，在這里我們只考慮人類從感染到恢復，因此蚊子沒有恢復者倉室[13-14]。模型的傳播機制如圖1所示。

(1)

其中Xh，Yh，Zh和Nh分別表示t時刻人類易感者，感染者，恢復者和所有人類的數量，Xm，Ym和Nm分別表示t時刻易感蚊子，感染蚊子和所有蚊子的數量。r是蚊子叮咬率，Thm(≤1)是受感染的蚊子叮咬易感者后易感者被感染的概率，Tmh(≤1)是易感蚊子叮咬感染者后蚊子被感染的概率。γh是人的恢復率；νh是人的總出生率，νm是蚊子的總出生率；μh是人的死亡率，μm是蚊子的死亡率。參數值的設置見表1。

表1 參數取值

對于蟲媒傳播的傳染病(如瘧疾、登革熱)，傳染率季節性起到了非常重要的作用。在這里我們考慮將模型中的季節性表示為蚊子叮咬率r，r是關于時間t的函數[13]：

r(t)=r0(1-r1cos(ωt+φ))

(2)

其中，r0是蚊子的平均叮咬率，ω是驅動頻率，φ是相位，r1是蚊子季節性叮咬率。

模型(1)中的傳染率矩陣為

(3)

通過最小化根據模型所得的數據與實際數據差值的平方和，即最小二乘法來得到最佳參數r0和r1的值，具體公式為

(4)

其中，DIi(t)為t時刻安徽省間日瘧的實際數據，PIi(t)為t時刻根據模型仿真得到的數據，T為仿真時長。利用R語言進行數據整理并建立模型，獲取使Sr最小的參數r0和r1，即是模型所估計的最佳參數值。

為了更準確地了解間日瘧傳播的動態范圍，我們計算基本再生數R0。基本再生數R0是一項重要的流行病學指標，預示著疾病傳播的潛在危險程度和控制的難易程度[19]。對于間日瘧來說，R0是指一只被感染的蚊子會產生二次感染的數量。基本再生數R0=1就是間日瘧流行與不流行的臨界水平，當R0>1時，即1個病例至少傳播產生1個繼發病例時，間日瘧才可以繼續傳播下去，否則傳播將趨向終止[20]。根據式(1)的均衡狀態，可以推導出間日瘧的基本再生數為

(5)

2.2 多元線性回歸模型分析傳染率影響因素

本文選取安徽省所有地級市氣象站點的月氣溫和降水量兩個因素的數值并分別求平均，與模型估計出的蚊子叮咬率進行回歸分析，由于溫度和降水量具有不同的量綱單位，為了消除指標之間的量綱影響，我們對兩者進行數據標準化處理，保證溫度和降水量處于同一數量級，適合進行綜合對比評價。多元線性回歸函數為

r(t)=aTemperature(t)+bRain(t)+c+ε

(6)

其中，蚊子叮咬率r(t)為響應變量，a、b為回歸系數，c為回歸常數，溫度Temperature(t)和降水量Rain(t)為回歸因子，ε是均值為零且方差不可觀測的隨機變量，稱為誤差項。

3 結果

應用最小二乘法對SIR-SI模型的參數進行估計，得出安徽省間日瘧的蚊子平均叮咬率r0為4.27e-06，季節性叮咬率r1為42.4%，相位φ為0.000 3。將蚊子叮咬r(t)帶入到傳染率矩陣中，即：

(7)

觀察到傳染率在每年的6、7月最高，12、1月最低。

對比模型仿真結果和實際數據，我們發現模型結果與安徽省2004～2009年間日瘧發病趨勢基本符合，且具有每年7～10月發病季節高峰的特點(見圖2)，但模型仿真結果與實際數據仍存在一定偏差。

將r0代入到式(5)中，計算出基本再生數R0的值為1.22。R0>1說明間日瘧在一定程度上流行，并觀察到發病數具有年周期性的特征。

多元線性回歸模型的結果表明溫度和降水量對間日瘧的蚊子叮咬率均有顯著影響，決定系數R2=0.73。其中溫度(℃)的回歸系數為2.609e-06(p=7.23e-11)，降水量(mm)的回歸系數為9.951e-07(p=0.004)，溫度的變化對于蚊子叮咬率的影響更顯著。此外由于溫度和降水量的方差膨脹因子VIF小于5，因此不存在顯著的共線性。

模型仿真結果與實際發病數的比較，-o-為實際值，——黑色為預測值。圖2 SIR-SI模型仿真結果

圖3 安徽省間日瘧預測感染者數(a)與預測易感者數(b)

圖4 安徽省間日瘧年報告發病數

假設2009年以后對間日瘧沒有進一步的控制措施，根據估計出的傳染率值應用SIR-SI模型對發病數趨勢進行預測，觀察到在2018年后間日瘧感染者人數維持在較低的狀態，在2034年后開始增加并于2048年達到峰值后又逐漸下降。預測的感染者人數不僅呈現出每年一個周期的特點，還存在多年一次大周期的特點(見圖3a)。在年感染者人數處于低谷期間，易感者人數則不斷累計，直至達到一定數值，感染者人數的均值和變化幅度又開始逐漸增加(見圖3b)。20世紀末的安徽省間日瘧年發病數呈現過再次上升的特點，在1993～1999年發病數較低，自2000年開始，間日瘧年發病數又不斷增加(見圖4)。

4 討論

本文以安徽省間日瘧為例，估算蚊子叮咬率、間日瘧的傳染率以及它們的季節性，并進一步分析傳染率發生季節性變化的影響因素。采用人類與蚊子相耦合的SIR-SI模型對間日瘧的傳播動態建模，估算隨時間變化的蚊子叮咬率以及傳染率。并進一步建立多元線性回歸模型探究溫度和降水量對蚊子叮咬率和間日瘧傳染率的影響。根據對安徽省間日瘧傳染率的分析結果，我們發現間日瘧的傳染率在一年中存在大幅度變化(變化幅度為42.4%)，傳染率在在6、7月最高，在12、1月最低。傳染病的傳染率季節性形成了發病數的周期性[13]，因而間日瘧發病數的年周期性波動是由較高的傳染率的季節性導致的。一些兒童傳染病的傳染率峰值早于發病數峰值若干星期或若干月[21]。本文研究發現，通過中華按蚊傳播的間日瘧的傳染率峰值也有類似的現象。間日瘧的蚊子叮咬率和傳染率又受溫度和降水量變化的影響，其中溫度對間日瘧傳染率的影響更大。有研究表明，溫度對中華按蚊的影響包括多個方面，例如繁殖、生存、活動行為等[22-23]。此外溫度對中華按蚊的吸血行為也存在顯著影響，在中華按蚊傳播間日瘧的溫度范圍內，隨著溫度的升高，雌蚊的叮咬頻率逐漸增加[24]。本文通過具有年季節性變化的叮咬率這個參數考慮了中華按蚊個體的叮咬頻率和中華按蚊數量變化對整體叮咬率的影響。

本文的研究結果顯示，間日瘧易感者人數存在多年周期性的變化趨勢，從而形成間日瘧感染者數多年的大周期。雖然間日瘧的感染者數能夠在多年維持較低的水平，但在易感者人數不斷累積并達到一定數量時，感染者數又會逐年增加，這能夠解釋90年代末安徽省間日瘧在保持若干年較低發病數后又逐年大幅增加的現象。值得注意的是90年代間日瘧開始保持較低的年發病數(<2 000)，至2000年發病數開始快速增加，保持較低發病數的時間長度為6年。本文的仿真結果顯示，2016年后間日瘧被感染人數也呈現類似的結果，但不同點是仿真結果在2016年后有更長的時間(約20年)處于較低發病數的狀態。考慮可能的原因是2016年后對間日瘧的控制措施優于90年代，降低了傳染率，這使得間日瘧的傳播速度相較于90年代更慢，因此2016年后間日瘧保持較低年發病數的時間更長。

本文中模型估計出的發病數與實際報告數據有一定偏差，產生的原因可能來自以下三個方面：1)最小二乘法對未知參數值的估計具有局限性，且仿真結果受初始值的影響較大，因此估計出的參數值并不一定是最優，未來我們考慮應用馬爾科夫蒙特卡洛方法(Markov chain Monte Carlo，MCMC)進行估計，得到更為準確的參數值；2)間日瘧的傳播動態實際上更接近于隨機過程。本文所使用的是確定性模型，該模型能夠在長時間內抓住間日瘧傳播動態的主要特征與總體趨勢，但對于具體年份上缺乏較好的描述。3)為了便于估計，我們沒有考慮各年間日瘧報告率的差異，報告率的差異會影響實際的發病數的多少。

安徽省作為中國間日瘧疫情不穩定的省份[11]，自90年代以來通過系列滅蚊措施，有效控制了間日瘧疫情，但是2000年后間日瘧發病數逐年上升并于2006年達到最高。2007年8月起安徽省在沿淮北地區的間日瘧重點防治縣實施預防服藥并采取多項綜合防治措施保證了間日瘧發病數穩步下降[25]。盡管通過人工干預取得了消除間日瘧的階段性勝利，但是由于傳播媒介中華按蚊的存在以及間日瘧可在人體內長期休眠等原因[26]，本地病例反復出現的風險依舊存在，因此我們要保持對間日瘧的監測，在通過省級消除瘧疾評估后仍要嚴防輸入再傳播，進一步防止本地疫情的反彈。此外據《第三次氣候變化國家評估報告》顯示，中國氣候變暖速率高于全球平均值。未來80年中國陸地區域增溫幅度將達到1.3～5.0℃[27]。由于溫度升高顯著增加間日瘧的傳染率，輸入性病例的增加以及累計的高易感者人數在增高的傳染率情況下是否會引起間日瘧的重新暴發是我們未來將要研究的內容。