白 迪，周武略，趙繼軍

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

0 引言

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

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

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

1 數(shù)據(jù)

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

2 方法

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

2.1 SIR-SI模型估算傳染率

圖1 SIR-SI模型示意圖

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

(1)

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

表1 參數(shù)取值

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

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

(2)

其中，r0是蚊子的平均叮咬率，ω是驅(qū)動(dòng)頻率，φ是相位，r1是蚊子季節(jié)性叮咬率。

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

(3)

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

(4)

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

為了更準(zhǔn)確地了解間日瘧傳播的動(dòng)態(tài)范圍，我們計(jì)算基本再生數(shù)R0?；驹偕鷶?shù)R0是一項(xiàng)重要的流行病學(xué)指標(biāo)，預(yù)示著疾病傳播的潛在危險(xiǎn)程度和控制的難易程度[19]。對(duì)于間日瘧來(lái)說(shuō)，R0是指一只被感染的蚊子會(huì)產(chǎn)生二次感染的數(shù)量?；驹偕鷶?shù)R0=1就是間日瘧流行與不流行的臨界水平，當(dāng)R0>1時(shí)，即1個(gè)病例至少傳播產(chǎn)生1個(gè)繼發(fā)病例時(shí)，間日瘧才可以繼續(xù)傳播下去，否則傳播將趨向終止[20]。根據(jù)式(1)的均衡狀態(tài)，可以推導(dǎo)出間日瘧的基本再生數(shù)為

(5)

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

本文選取安徽省所有地級(jí)市氣象站點(diǎn)的月氣溫和降水量?jī)蓚€(gè)因素的數(shù)值并分別求平均，與模型估計(jì)出的蚊子叮咬率進(jìn)行回歸分析，由于溫度和降水量具有不同的量綱單位，為了消除指標(biāo)之間的量綱影響，我們對(duì)兩者進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理，保證溫度和降水量處于同一數(shù)量級(jí)，適合進(jìn)行綜合對(duì)比評(píng)價(jià)。多元線性回歸函數(shù)為

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

(6)

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

3 結(jié)果

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

(7)

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

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

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

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

模型仿真結(jié)果與實(shí)際發(fā)病數(shù)的比較，-o-為實(shí)際值，——黑色為預(yù)測(cè)值。圖2 SIR-SI模型仿真結(jié)果

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

圖4 安徽省間日瘧年報(bào)告發(fā)病數(shù)

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

4 討論

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

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

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

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