基于Kolmogorov前向方程評估甲型H1N1流感疫情的動態(tài)變化*

閆琴玲，唐三一

（1.長安大學(xué) 理學(xué)院，西安 710064；2.陜西師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，西安 710119）

引 言

數(shù)學(xué)模型在突發(fā)性傳染病疫情防控中發(fā)揮著至關(guān)重要的作用.它們可以用來刻畫傳染病傳播規(guī)律、評估控制措施的有效性以及預(yù)測潛在傳染性等.而針對突發(fā)性傳染病，大多數(shù)研究者主要是建立基于種群水平的確定性模型.然而，越來越多的研究表明，由于人口事件的隨機性，人口隨機效應(yīng)可能導(dǎo)致確定性模型的結(jié)果與實際情形顯著偏離[1-5].因此，一般情形下，人們基于事件驅(qū)動的模型模仿真實系統(tǒng)的現(xiàn)實行為[6-7]，即將一系列復(fù)雜的生物現(xiàn)實行為嵌入基礎(chǔ)模型，并為人們提供直觀的建?？蚣?然而，對于單次模擬，不能確定模擬的動力學(xué)是代表平均行為還是僅僅是由于偶然事件產(chǎn)生的異常值.因此，為得到對應(yīng)變量(或參數(shù))的置信區(qū)間，研究者們需要進行大量的重復(fù)模擬.

許多現(xiàn)存的近似方法(如擴散近似和矩閉合技術(shù))通過給出近似解析解以避免大量模擬種群行為及其變化[5,8].然而，只有當(dāng)種群規(guī)模足夠大時，這些方法才是比較準(zhǔn)確的.而通常研究的問題是種群規(guī)模較小的隨機模型，且隨機性對動力學(xué)行為具有相對較大的影響，需要探索適用于種群規(guī)模較小的方法.因此，我們應(yīng)用一種Markov 過程——Kolmogorov 前向方程(KFE)，來研究種群規(guī)模較小時的隨機模型.

KFE 是一個關(guān)于處于每種狀態(tài)概率的動態(tài)微分方程系統(tǒng)，其動態(tài)變化大小由狀態(tài)之間的轉(zhuǎn)換速率決定[9-10].因此，通過求解微分方程系統(tǒng)，可以得到隨機系統(tǒng)所有可能行為的完整描述.并且隨著計算機計算能力的不斷提高，這些技術(shù)越來越適用于實際問題.用KFE 方法生成大量微分方程的解，在現(xiàn)代計算機上運行相對簡單且快速[5,11-14].此外，就處于每種狀態(tài)的概率而言，它是線性的.它通常把與種群動力學(xué)相關(guān)的復(fù)雜非線性簡單地包含到矩陣項中.通過此方法，可以用更簡單的矩陣和向量運算來構(gòu)建方程，這大幅度地加快了計算速度.

另一方面，已有文獻研究表明，在突發(fā)性傳染病爆發(fā)過程中，人們的行為可能發(fā)生改變(如戴口罩、勤洗手、避免到公共場所等).而由于個體行為改變，感染風(fēng)險會降低[15-16].

為此，本文基于2009年西安市第八醫(yī)院甲型H1N1 流感數(shù)據(jù)(包括新增病例數(shù)及累計病例數(shù))，如圖1所示，建立個體決策心理模型；把行為改變率嵌入基于個體水平的SIR 模型，并基于此推導(dǎo)得到KFE，探究此傳染病在整個傳播過程中每種狀態(tài)的概率分布，從而高效地進行甲型H1N1 疫情的防控.

圖1 2009年9月3日至10月10日西安市第八醫(yī)院甲型H1N1 病例數(shù)：(a)新增病例數(shù)； (b)累計病例數(shù)Fig.1 The hospital notifications for A/H1N1 from 3rd September to 10th October 2009 in Xi’an 8th Hospital: (a)the number of new hospital notifications;(b)the accumulated number of hospital notifications

1 基于甲型H1N1 數(shù)據(jù)的個體決策心理模型

基于文獻[17]的觀點，健康信念模型(health belief model, HBM)結(jié)構(gòu)可用于揭示突發(fā)性傳染病爆發(fā)期間的行為決策與情景信息之間的關(guān)系.因此，本節(jié)主要借助于HBM 結(jié)構(gòu)與情景信息之間的關(guān)系(包括自覺嚴(yán)重性、自覺易感性、自覺障礙、自覺效益)，建立具有以上四種HBM 結(jié)構(gòu)和行為改變率的logistic 模型(記為LHBM)，采用非線性最小二乘法對模型中的參數(shù)進行估計.LHBM 結(jié)構(gòu)的詳細描述在附錄以及圖A1 給出.

為了估計方程(A6)中的未知參數(shù) (O0，O1，O2，O3，O4，δ，λ，ξ，τ，η )，基于四種社交網(wǎng)絡(luò)、LHBM 流程圖(圖A1)以及文獻[18]，我們模擬從9月3日到10月10日，2 004 個個體(每個個體每天基于LHBM 結(jié)構(gòu)獨立作出決定)的行為改變率，并對行為改變率數(shù)據(jù)進行擬合.2 004 個個體(假設(shè)易感人數(shù)2 000，感染人數(shù)4，恢復(fù)人數(shù)0)中每個個體的行為改變率的估計值pe(i,t)，所有個體行為改變率的平均值(t)和其95% 置信區(qū)間如圖A2所示，其中黑色區(qū)域、白色曲線分別代表2004 個個體的95%的置信區(qū)間及均值((t))，圓圈表示p(t)的真實值.估計得到的參數(shù)的平均值及標(biāo)準(zhǔn)差如表 A1 所示.由于規(guī)則網(wǎng)絡(luò)殘差絕對值的和最小，因此，由規(guī)則網(wǎng)絡(luò)得到的平均行為改變率pˉ(t)將會用于以下模型的研究.

2 基于甲型H1N1 數(shù)據(jù)的KFE

2.1 KFE 的建立

基于文獻[18]中的模型選擇結(jié)果，這里，我們考慮行為改變率(t)，并將平均的行為改變率pˉ(t)以指數(shù)函數(shù)的形式嵌入經(jīng)典的SIR 模型，得到如下系統(tǒng)：

其中S(t)，I(t)，R(t)分別表示易感人數(shù)、感染人數(shù)和恢復(fù)人數(shù).

根據(jù)KFE 的有關(guān)知識[9]，模型(1)的反應(yīng)式可表示為

記X(t)=[S(t),I(t),R(t)]T(T 為轉(zhuǎn)置，t≥0)為系統(tǒng)(1)中各倉室的人數(shù)，已知初值X(0)=x(0)，則易知{X(t),t≥0}是 一個多元隨機過程.設(shè)Z(t)=[Z1(t),Z2(t)]T為模型(2)中兩種反應(yīng)在 [0,t)內(nèi)發(fā)生的次數(shù)，則{Z(t),t≥0}是一個計數(shù)過程，稱為狀態(tài)更新度(DA 過程)[19].依據(jù)DA 過程與群體過程之間的關(guān)系，有如下方程[9]：

其中

對于給定的x(0)，X(t)由Z(t)唯一確定.

基于系統(tǒng)(2)，對于任意的m(這里m∈{1,2})，在 [t,t+dt] 內(nèi)，反應(yīng)m發(fā)生的概率是 πm(x)dt+o(dt)[20].其中，π1(x)=βe?αpˉ(t)S(t)I(t)/N和 π2(x)=γI(t)是 系統(tǒng)(2)在時刻t的傾向函數(shù).因此，由于 {Zm(t),t≥0} 是Markov 過程，強度為 πm(X(t))，經(jīng)過推導(dǎo)可以得到以下KFE：

其中pZ(z;0)=δ(z)，δ (z)是 一個Kronecker δ 函數(shù)；e1=[1,0]T,e2=[0,1]T；對于m∈{1,2}，有

依據(jù)方程(3)中Z(t)和X(t)的關(guān)系，可得pX(x;t)的概率函數(shù)滿足

其中 B (X)={z|x=x(0)+Sz}.因此，通過求解方程(5)可以得到方程(3)的動態(tài)變化特性.

此外，對于 T =[0,tmax]，有限樣本空間 Z (元素zk=[z1k(t),z2k(t)]T,k=1,2,···,K)，基于概率的性質(zhì)，則pZ(z;t)滿足

令K×1 向量 φ(t)的 元素 φk(t)=pZ(zk;t)(k=1,2,···,K).則基于方程(5)，有如下K維微分方程：

對 于給定的時刻t，A(t)是K×K矩陣；初始值φ (0)=[1,0,···,0]TK×1.

2.2 KFE 的數(shù)值求解

對于給定的t(t=1,2,···,38)，易知矩陣A(t)具有三角性、稀疏性、穩(wěn)定性[9].因此，對于以下方程

采用隱式Euler(IE)法計算 φ(t)在 離散時間點tj=jτ(j=1,2,···;τ=0．01)的估計值(tj).這里I是單位矩陣.因此，對于給定的(0)=φ(0)和tj(j=1,2,···)，通過方程(10)可以得到概率質(zhì)量函數(shù)pZ(z;t).

基于文獻[18]，設(shè) α =0．011 57, β =0．755 5, γ=0．545 2為參數(shù)的初始值.又因在疫情初期有4 個感染個體，因此設(shè)S(0)=2 000,I(0)=4,R(0)=0 為易感人數(shù)、感染人數(shù)和恢復(fù)人數(shù)的初始值.(t)是第1 節(jié)中N=2 004個個體行為改變率的平均值.K=(S0+1)×(S0+I0+1)=2 001×2 005=4 012 005是樣本空間的大小.因此，在給定精度下，利用IE 法計算KFE(10)的精確解，從而得到模型(5)的精確解.同時，對于模型(1)，采用Gillespie 算法進行參數(shù)估計和模型擬合(基于1 000 次Markov-chain Monte-Carlo (MCMC)模擬)，從而對結(jié)果進行進一步的比較[6,20].

考慮到R(t)=N?S(t)?I(t)及無病狀態(tài)發(fā)生的概率很高，我們只需要計算感染人數(shù)與易感人數(shù)的聯(lián)合條件概率質(zhì)量函數(shù)P(S(t),I(t)|I(t)>0).圖2(a)和(b)分別刻畫了在第38 天的概率質(zhì)量函數(shù)P(S(t),I(t)|I(t)>0)及其等高線圖.由圖2可知，假設(shè)最初有4 個感染個體，到第38 天，感染人數(shù)與易感人數(shù)的聯(lián)合條件概率達到最大(1 ×10?4)，此時新增病例數(shù)約為 5 ～15，易感人數(shù)約為 750 ～850，這與實際情況基本一致(圖1(a)).其次，圖2(c)給出了第38 天恢復(fù)人數(shù)的條件概率質(zhì)量函數(shù)P(R(t)|I(t)>0).由圖可知，到第38 天時，恢復(fù)的概率達到最大(3 ．5×10?3)，約有1 200 個個體恢復(fù).

圖2 流感流行38 天計算出的概率質(zhì)量函數(shù)的快照： (a)到第38 天時，感染人數(shù)與易感人數(shù)的聯(lián)合條件概率質(zhì)量函數(shù)P (S(t),I(t)|I(t)>0)；(b)概率質(zhì)量函數(shù)的輪廓圖；(c)第38 天恢復(fù)人數(shù)的條件概率質(zhì)量函數(shù)P(R(t)|I(t)>0)Fig.2 A snapshot of the calculated probability mass function at the end of the 38 d influenza epidemic: (a)joint conditional probability mass function P(S(t),I(t)|I(t)>0)of infected and susceptible individuals; (b)the contour of the calculated probability mass function; (c)conditional probability mass function P (R(t)|I(t)>0)of recovered individuals

再次，易感人數(shù)、感染人數(shù)和恢復(fù)人數(shù)均值(均值 ±標(biāo)準(zhǔn)差)的動態(tài)變化可以根據(jù)P(S(t),I(t)|I(t)>0)和P(R(t)|I(t)>0)直接計算得到，如圖3(a)、(b)和(c)所示.通過比較圖3(b)與圖3(d)，得到IE 法的結(jié)果比Gillespie算法的結(jié)果更好.此外，Gillespie 算法需要進行大量隨機模擬才能準(zhǔn)確估計P(S(t),I(t),R(t)).因此，與Gillespie 算法相比，IE 法更優(yōu).

3 最終規(guī)模

根據(jù)Kermack 等的定義，最終規(guī)模R∞是指，對于突發(fā)性傳染病，引入少量感染個體后，被感染人數(shù)在易感人群中(并最終恢復(fù))的比例[21].因此，根據(jù)文獻[21]，關(guān)于最終規(guī)模有以下關(guān)系式：

其中R0是基本再生數(shù).

3.1 R0的計算

對于方程(1)，根據(jù)文獻[22-23]中R0的定義(指引入一個感染個體，在感染周期內(nèi)，平均出現(xiàn)的新感染個體的數(shù)量)，設(shè)新發(fā)感染和轉(zhuǎn)移項分別用一維向量F，V表示，則

從而得到下一代一維矩陣

根據(jù)Wang 等對譜半徑的定義，R0可由感染算子的譜半徑L計算得到[23]，可表示為

對于周期的傳染病系統(tǒng)(1)，根據(jù)Becar 等的計算方法可得[22]

其中參數(shù)λ ∈(0,+∞).設(shè)系統(tǒng)(1)在 Rm上的演化算子為W(t,s,λ)，則R0為 ρ (W(ω,0,λ))=1的正根，其值為

這里最大的觀測數(shù)據(jù)時間點T1=38.

對于系統(tǒng)(1)，根據(jù)文獻[18]得到的參數(shù) α , β , γ 的值，得R0的估計值為1.314 3(95% 置信區(qū)間1.3127～1.3160).

3.2 R?∞的計算

令x=R∞(x∈[0,1))，則方程(11)可表示為

記f(x)=1?e?R0x?x，則易知x?=和x?=0是f(x)=0的兩個根.

f(x)對x求導(dǎo)，得

因此，對于方程(17)的根，根據(jù)參數(shù)R0的閾值水平，考慮以下三種情形：

(ⅰ)當(dāng)R0<1時，f′(x)<0，即f(x)在 (0,1)上 單調(diào)遞減.又因為f(0)=0，因此，當(dāng)R0<1時，x?=0是方程(17)的唯一解.

(ⅱ)當(dāng)R0=1時，f′(x)=e?x?1<0 .因此，當(dāng)R0=1時，x?=0是方程(17)的唯一解.

(ⅲ)當(dāng)R0>1時，設(shè)f′()=0，則=ln(R0)/R0>0.易驗證<1顯然成立.

又因為

(a)當(dāng)x<時，f′(x)>0； 當(dāng)x>時，f′(x)<0； 因此，f(x)在 區(qū)間(0,)上單調(diào)遞增，在區(qū)間(,1)上單調(diào)遞減.

(b)f(1)=?e?R0<0，f′()=1?(1+lnR0)/R0<1且f′>0.因此，存在唯一的解x?∈(1)使 得f(x?)=0,x?的顯式表達式如下：

綜上，方程(17)有兩個解，分別為0 和

其中

4 討論與結(jié)論

本文建立了一個基于個體水平的隨機SIR 模型，通過Markov 過程，推導(dǎo)出KFE.然而，只有對含有限個體的群體，才能得到KFE 的解析表達式[10,24].對于含有大量個體的群體，只能采用近似法或數(shù)值解方法得到其解.對于隨機的SIR 模型，基于DA 過程，Jenkinson 等[9]采用IE 法得到KFE 的精確解，并且證實了此方法優(yōu)于Krylov 子空間逼近法及Gillespie 算法.因此，本文通過IE 法得到了KFE 的數(shù)值解、傳染病傳播過程中感染人群與易感人群的聯(lián)合條件概率質(zhì)量函數(shù)、最嚴(yán)重的時間段和相應(yīng)的概率.此外，對于隨機SIR 模型的求解，IE 法比Gillespie 算法更好.

我們的研究有以下局限性.首先，考慮到計算成本和網(wǎng)絡(luò)的復(fù)雜性，在第1 節(jié)我們僅選擇對2 004 個個體進行隨機模擬.另外，由于A(t)是K×K矩陣(K=(S0+1)×(S0+I0+1))，隨著S0和I0增 大，A(t)計算量大幅度增加，考慮計算成本，在第2.2 小節(jié)數(shù)值計算仍然選擇總?cè)藬?shù)為2 004.

綜上，本研究通過基于個體水平的隨機SIR 模型，推導(dǎo)并數(shù)值求解KFE，得到了疫情爆發(fā)過程中最嚴(yán)重的時間段以及對應(yīng)的概率分布，從而能更快、更準(zhǔn)確地了解甲型H1N1 流感的傳播過程，因此有助于高效地進行甲型H1N1 的疫情防控.

附錄 LHBM 的建立

LHBM 模型的詳細描述如下[18].

自覺易感性自覺易感性(perceived susceptibility)是指人們對感染傳染病最大可能性的信念[26]，記作P1，定義如下[27]：

式中

其中參數(shù)λ 是低敏感性到高敏感性的閾值水平，st，ct?1，δ分別表示累計病例數(shù)、新增病例數(shù)以及折現(xiàn)率常數(shù).

自覺嚴(yán)重性自覺嚴(yán)重性(perceived severity)是指人們對一種傳染病嚴(yán)重程度與其后遺癥所持的觀點[26].這里假設(shè)自覺嚴(yán)重性是媒體報道量的函數(shù)(記作M(t))，記作P2，即

其中參數(shù) τ 表示個體自覺嚴(yán)重性從高到低的閾值水平.

自覺障礙自覺障礙(perceived barriers)用于描述人們采取一些行動的成本，包括有形成本和心理成本[26]，記作P3.我們使用四種經(jīng)典網(wǎng)絡(luò)刻畫個體之間的連接關(guān)系，包括小世界(WS)網(wǎng)絡(luò)[28]、規(guī)則網(wǎng)絡(luò)、隨機網(wǎng)絡(luò)[29-30]、Newman 和Watts(NW)網(wǎng)絡(luò)[31].從而有(詳見文獻[18])

其中，集合 Ωi={j|j≠i} 表示所有與個體i有連接的個體，ni為集合 Ωi的大??；當(dāng)個體j戴口罩時，Aij等于 1，否則等于 0；參數(shù) η 是低障礙到高障礙的閾值；Dij表示i和j之間的親密關(guān)系，即Dij=md/dij，其中，md是i和j之間的最小距離；社交網(wǎng)絡(luò)中總?cè)藬?shù)為N1，有N2個戴口罩的人,與i連接的鄰近節(jié)點數(shù)記為c.

自覺效益自覺效益(perceived benefits)是人們?yōu)榻档陀蓚魅静∫鸬娘L(fēng)險而對采取某些行動有效性的信念[26]，這里設(shè)自覺效益接觸人數(shù)的函數(shù)[29]，記作P4，定義如下(詳見文獻[18])：

其中，參數(shù)ξ 是低效益到高效益的閾值水平.

因此, 具有以上四種HBM 結(jié)構(gòu)和行為改變率的logistic 模型(記為LHBM)，定義如下[27]：

這里,k=1,2,3,4 分 別代表自覺易感性、自覺嚴(yán)重性、自覺障礙和自覺效益；p(i,t)是第i個體第t天的行為改變率；Ok表示兩種結(jié)構(gòu)的比率；xk(i,t)的取值為 1 表示LHBM 結(jié)構(gòu)為高狀態(tài)，為 0 表示低狀態(tài);O0是校準(zhǔn)常數(shù).

基于四種網(wǎng)絡(luò)，2 004 個個體從9月3日至10月10日估計的行為改變率(pe(i,t))和實際的行為改變率(p(t))之間殘差絕對值的和分別為2.314 2E+3，2.785 7E+3，2.510 8E+3，2.403 1E+3.

圖A1 基于個體的LHBM 流程圖Fig.A1 The individual-based scheme of the LHBM

圖A2 基于模型(A6)，對2 004 個個體的個體行為改變率 pe(i,t)進行估計： (a)規(guī)則網(wǎng)絡(luò); (b)小世界(WS)網(wǎng)絡(luò); (c)NW 網(wǎng)絡(luò); (d)隨機網(wǎng)絡(luò)Fig.A2 The estimation of individual behaviour change rate pe(i,t)for 2 004 persons based on model (A6): (a)the regular network; (b)the small-world (WS)network; (c)the Newman and Watts (NW)network; (d)the random network

表A1 在四種經(jīng)典的社交網(wǎng)絡(luò)結(jié)構(gòu)下, 模型(A6)所有參數(shù)均值和標(biāo)準(zhǔn)差的估計值及取值范圍Table A1 The estimated values and ranges of means and standard deviations for all parameters of model (A6)with 4 classical social networks