具有年齡階段的離散SCIRS模型在我國流腦中的應用?

馬 霞,周義倉

(1-太原工業學院理學系，太原 030008;2-西安交通大學數學與統計學院，西安 710049)

1 引言

流行性腦脊髓膜炎(流腦)是世界性流行病，它是一個重要的公共衛生問題．據估計，全世界大約每年有120萬人被感染，135000人死亡[1]．因其發病率及病死率高被人們所關注．全世界每年發生30萬至35萬流腦病例，各大洲發病率在1/10萬至10/10萬，總病死率在5%～10%[2]．流腦被發現并確認已超過100多年，盡管對此病現在有些有效的預防與治療措施，但它仍是一種常見的急性傳染病，而且病死率高，繼續威脅著人類，尤其是兒童的身體健康[3]．大約10%的人在鼻咽部攜帶腦膜炎球菌，腦膜炎球菌在青少年和年輕的成年人中比較流行，嬰幼兒很少攜帶，然而侵入性腦膜炎在嬰幼兒中比較常見．

由于我國在推廣腦膜炎球菌多糖疫苗(MPV)接種前，每3～5年發生一次小流行，發病率為10/10萬～50/10萬；每8～10年發生一次大流行，發病率可高達100/10萬～500/10萬[2]．現行的A群、C群MPV對大年齡兒童和成人的免疫持久性不長，對嬰幼兒的保護效果差，且覆蓋率不高，所以常出現局部爆發[4]；目前按照國家擴大免疫規劃疫苗接種程序的規定，6月齡起即可接種A群流腦疫苗，如果提高A群流腦疫苗的及時接種率，使更多的低月齡兒童及時獲得免疫，可以進一步降低6月齡以上的低月齡兒童的流腦發病率[5]．近期的監測數據表明，疫苗應對攜帶者的保護期可持續3～10年，群體免疫的影響可持續更長的時間[6]．因此，我們的模型中將考慮疫苗因素對流腦傳播的影響．

由于接種疫苗不能阻止攜帶者攜帶腦膜炎奈瑟菌，但是它能刺激產生群體免疫進而減少疾病的爆發，而MCC疫苗被證明能顯著減少鼻咽部腦膜炎球菌的攜帶[7]，疫苗刺激產生的群體免疫對腦膜炎疾病的爆發率有很大的影響．隨著年齡的增加，攜帶者產生的免疫能減少發病的風險[8]．因此，我們的模型考慮對青少年期進行疫苗補種來提高群體的免疫力．根據流腦的傳播機理和統計數據可知，嬰幼兒期感染腦膜炎奈瑟菌發病和死亡的機率比較大．不同年齡階段的個體，感染腦膜炎奈瑟菌的機率是不一樣的，并且個體感染腦膜炎奈瑟菌后，感染時間的長短也會導致個體的發病情況不同，流腦的攜帶時間也是因人而異的，僅僅一小部分攜帶者發展成為侵襲性的染病者、不同年齡段、不同感染年齡段的個體、從攜帶潛伏期發展為染病者的情況也是不同的．因此，我們會在模型中引入年齡結構這個因素，希望利用具有年齡結構的傳染病模型來更好的探討疾病的傳播機制．

目前關于離散動力系統方面的知識還不完備，使得對一般自治的離散傳染病模型的動力學研究相對比較困難，對年齡結構的離散模型的研究更加困難．文獻[9]中討論了一類具有年齡結構和移民因素的離散SIS傳染病模型的全局穩定性態．也有學者針對水痘和帶狀皰疹建立具有年齡結構的離散模型，研究不同控制策略對于疾病傳播的影響．文中將研究具有年齡階段、免疫接種的離散SCIRS模型的動力學性態，并用模型來預測分析我國流腦的流行傳播規律．

2 模型的描述

針對流腦的流行特征，為區分不同年齡段攜帶者和染病者的不同比例，將人群劃分為三個年齡階段，第一階段(0～4歲)，第二階段(5～14歲)和第三階段(15歲以上)，用N1,N2,N3分別表示三個階段的總人口數量，在每個階段中，將人口劃分為4個倉室，易感者Si，攜帶者Ci，染病者Ii和免疫者Ri，其中i=1,2,3．假設攜帶者Ci和染病者Ii都具有傳染性，用ε(0<ε<1)表示Ci相對于Ii的傳染力，不考慮交叉感染，當一個易感者被攜帶者和染病者感染后，隨之進入到攜帶者倉室，在攜帶期有兩種情況出現，若宿主自身免疫系統不能清除體內病毒，以比例ai(i=1,2,3)發展為染病者．若宿主自身免疫系統能清除體內病毒，攜帶者將失去攜帶以比例αi(i=1,2,3)轉變為免疫者．不考慮母嬰垂直傳染，假設新生兒的免疫失敗率相同都為p，k1為對5～14歲易感者突擊接種的有效率，k2對5～14歲攜帶者突擊接種的有效率，染病者的恢復率為γ，隨著免疫的消退，免疫者會以一定比例φi(i=1,2,3)失去免疫回到易感者類．μi(i=1,2,3)為自然死亡率，δi(i=1,2,3)為因病死亡率．Λ表示人口的常數輸入，u1,u2分別表示從第一階段到第二階段和從第二階段到第三階段的轉移率．圖1給出了流腦傳播流行的框圖，利用以上的記號和假設，建立具有三個年齡階段的流腦傳播模型如下：

圖1:具有三個年齡階段的流腦流行傳播轉換圖

根據實際問題，我們給出下面的條件

把模型(1)中的方程分別相加，可得t時刻總人口Ni(t)滿足

考慮輔助系統

由上述方程組可以得到

由文獻[10]中的比較定理可知，Ni(t)(i=1,2,3)是最終有界的，從而系統(1)的解在區間上全局存在．易驗證區域

是系統(1)的完備的正不變集，從?中出發的解始終停留在?中．

利用文獻[11,12]中再生矩陣的方法，可得系統(1)的基本再生數f R0= ρ(F(I?T)?1)，即ρ(F(I?T)?1)表示矩陣F(I?T)?1的譜半徑，其中I是6階單位矩陣

這里

3 疾病的消亡與持續

疾病的消亡或持續是由模型(1)的無病平衡點的穩定性和地方病平衡點的存在性決定的．下面的定理給出了無病平衡點的全局穩定性和地方病平衡點的存在性．

定理1 當<1時，模型(1)僅存在一個無病平衡點P0且是全局漸近穩定的，疾病將會消亡．當>1時，無病平衡點P0不穩定，模型存在一個地方病平衡點P?．

證明 模型(1)在無病平衡點P0處的雅克比矩陣為

我們考慮下面的方程

方程(3)的系數矩陣是

如果R01<1，通過計算可得ρ(A1)<1．由文獻[13]中的定理6.1.3可知，方程(3)的平衡點是全局漸近穩定的，即

同理，可得

4 SCIRS模型在我國流腦傳播中的應用

前面給出了SCIRS模型的理論結果，下面我們把模型應用到我國流腦流行傳播中．基于我國法定報告的流腦流行的數據，分析預測流腦在我國人群中的發展趨勢．進行數值模擬時，所有參數的選擇以天為單位，2006～2012年我國人口的平均出生率為12/1000，平均死亡率為7.06/1000[14,15]，取常數輸入Λ=44054，0～4歲人群的平均死亡率為3.56/1000，5～14歲人群的平均死亡率為0.406/1000，15歲以上人群的平均死亡率為6.85/1000[14,15]．因此，μ1=0.00000975,μ2=0.00000111,μ3=0.00001877．假定新生兒的免疫失敗率為p=0.44，Ci相對于Ii的傳染力ε=0.7．根據法定傳染病報告的數據可得，2006～2011年腦膜炎患者平均每年因病死亡率為10.48%[16]，平均每天的因病死亡率為0.0002874，根據文獻[16,17]，我們分別取δ1=0.0004704,δ2=0.0001222,δ3=0.0002175．由文獻[18]可知，腦膜炎的平均攜帶時間為1周至1年，假設第一階段、第二階段、第三階段的攜帶時間分別為三個月、四個月、六個月，因此，α1=0.0111,α2=0.00833,α3=0.00555．由于攜帶者中只有很少一部分會發展為染病者，攜帶者失去攜帶和發展為染病者的比率一般大于100:1，通常αi>ai[19]，考慮到0～4歲人群患病率比較大，取a1=0.000111,a2=0.0000833,a3=0.0000555，染病者的恢復期為一周[18]，即γ=0.1428．假設對第二階段的易感者每年以30%的比例突擊接種，攜帶者以20%的比例突擊接種，因此k1=0.0008219,k2=0.0005479，由文獻[20]知，第一階段人群計劃接種疫苗的保護期為5年，第二階段人群突擊接種疫苗的保護期為10年，第三階段人群的免疫保護期為15年，所以φ1=0.0005479,φ2=0.000274,φ3=0.0001826．每天從第一階段到第二階段和從第二階段到第三階段的轉移率u1=0.0005479,u2=0.000274．采用文獻[21]中參數估計的方法，我們取β1=0.01648,β2=0.01048,β3=0.00748．模擬時，以2006年1月1日的數據為初始值，根據統計年鑒[14,15]中的我國人口的數據可知，2006年我國總人口數量為N(1)=1314945975，0～4歲人群總數量為N1(1)=66765160，5～14歲人群總數量為N2(1)=176101433，15歲以上的人群總數量為N3(1)=1072079382．由法定傳染病報告的數據[16]可知，2006年全年新發腦膜炎病例數為1816例，因此，我們可直接計算2006年1月1日染病者總人數大約為I(1)=1816÷365÷γ≈35，由文獻[17]中流腦病人的年齡構成可推知I1(1)=12,I2(1)=10,I3(1)=13．第一階段攜帶者人數為C1(1)=I1(1)×γ÷a1=15428，第二階段攜帶者人數為C2(1)=I2(1)×γ÷a2=17142，第三階段攜帶者人數為C3(1)=I3(1)×γ÷a3=33428．假定第一階段流腦病人的平均年齡為1歲，第二階段流腦病人的平均年齡為8歲，第一階段恢復者人數為

當參數和初始值選定后，我們用上述三個年齡階段的SCIRS模型來模擬我國2006～2020年流腦的發病情況．模擬結果如圖2所示，圖2(a)預測的是2006～2020年每天總人群中的易感者S的數量，開始有下降趨勢，最后基本趨于平衡狀態，由于考慮到對第二階段人群每年以一定比例進行突擊接種，有部分易感者接種直接進入恢復者倉室，因此，易感者的數量開始會緩慢減少，最后趨于平衡狀態．圖2(b)預測的是2006～2020年每天總人群中的染病者I的數量，有下降的趨勢，開始下降的快，最后趨于平緩但不等于0．圖2(c)是預測的2006～2020年每天總人群中的攜帶者C的數量．圖2(d)中的曲線是模擬的2006～2020年每年新發病的流腦病例，折線上的星號代表2006～2013年每年我國法定傳染病報道的流腦病例數[16]．數值模擬的結果和法定報告的傳染病數據進行對比可知，該模型可以用來預測我國流腦的流行發病情況．

5 基本再生數關于參數的敏感性分析

由于數值模擬中所采用的參數大部分是不確定的，并且這些參數直接影響到f R0的大小，因此，有必要對這些參數進行不確定性和敏感性分析，從而找出起決定作用的參數，進一步驗證模型的正確性．這里我們在LHS基礎上給出一些參數的PRCC值，有關LHS和PRCC的詳細介紹見文獻[22]，我們抽取的樣本大小為N=3000，除去Λ,γ,ε,δi,μi,ui以外的所有參數作為輸入變量，f R0的大小為輸出變量，其中各參數的初始值和數值模擬中的參數值一致，取值范圍服從[0,1]上的均勻分布，由參數和其PRCC值的散點圖可以發現，均為單調的，故PRCC方法適用．圖3給出了f R0關于參數的敏感性分析及一些參數的PRCC值，由圖可知參數p,φ1,φ2,φ3,k1不敏感，對f R0的影響均不大，β1,β2,β3,a1,a2,a3與基本再生數f R0正相關，α1,α2,α3,γ,k2與基本再生數呈負相關，這完全符合我們的期望．感染率對f R0的大小起決定性作用，感染率越大，會使疾病的蔓延越快，攜帶者的發病率越高，也會提高疾病的蔓延速度，攜帶者失去攜帶的比例越大，會降低發病的比率，減少疾病的蔓延，由于k1對f R0不敏感，因此，擴大對攜帶者補種接種的比例比擴大易感者補種接種更有效，減小基本再生數可以首先從降低和控制這些參數來入手．我們可以通過監測提高無癥狀的攜帶者的發現率，擴大免疫接種，減少染病者的高危行為，對于控制疾病有重要的意義．

6 考慮季節因素的三個年齡階段模型

由于季節性變化在很多傳染病中扮演著重要的作用，季節變化影響接觸率和出生率能最終導致復雜的動力學行為．針對我國報告的法定傳染病數據可知，流腦發病一般自每年的12月份開始明顯上升，至次年的3月份達到高峰期，4月份發病開始呈下降趨勢．流腦流行有明顯的季節性，無論流行年或散發年一般均在冬春季流行[3]．因此，在上述的年齡階段模型中可以考慮季節因素對流腦流行的影響，對于季節因素對流腦發病的影響，有一部分學者[23]認為在干燥的冬春季節，傳染病通過空氣中的飛沫更容易傳播，導致攜帶者的數量增加進而導致染病者的數量增加．也有學者[24]認為在冬春季節，天氣干燥，塵土飛揚可能損害鼻咽部位的粘膜屏障，使得無癥狀的攜帶者更容易轉變為染病者．為比較季節變化對疾病傳播速率以及疾病進展速率的影響，我們考慮疾病的流行發病周期為一年，由于流腦發病高峰一般在2～4月份，因此，考慮的參數值在干旱的季節明顯比雨季高．用下面的周期表達式替換常數速率

圖3: 基本再生數f R0關于參數的敏感性分析

假設疾病的傳播速率β受季節因素的影響，而疾病進展率a是常數時，取εβ=0.4，而其他參數的取值均不變，數值模擬結果如圖4所示，和圖2相比，可知S(t),I(t),C(t)的基本走向是一致，但是I(t),C(t)出現周期振蕩的現象，并且振幅隨時間逐漸減小，最后趨于平緩，說明當疾病的傳播速率受季節因素的影響時，對流腦的傳播和發病情況都產生很大的影響，這也驗證了季節因素影響流腦發病的第一種假設．

假設疾病的傳播速率β是常數，而疾病的進展率a受季節因素的影響，取εa=0.4，其他參數值和初始值的選取均不變，數值模擬結果如圖5所示．與圖2相比，很容易看出S(t),C(t)的變化不太明顯，I(t)出現周期振蕩的現象，開始振蕩的劇烈，振幅隨時間逐漸減小，最后趨于平緩，說明疾病的進展速率受季節變化的影響時，對流腦發病有很大的影響，和季節因素影響的第二種觀點是一致的．與圖4對比，可知I(t)振蕩的幅度比較大，而C(t)沒有出現周期振蕩的現象．說明疾病進展率a受季節因素的影響時，對攜帶者C(t)的數量影響很小，而對染病者I(t)的數量影響很大．

假設疾病的傳播速率β和疾病的進展率a同時受季節因素的影響，取εβ=0.4,εa=0.3，數值模擬結果如圖6所示，我們可以看出C(t),I(t)均出現周期振蕩的現象，而I(t)振蕩的幅度比較大，該結果和只考慮疾病的傳播速度β受季節因素的影響時模擬的結果是一致的．但是和圖4相比，I(t)振蕩的幅度較大，呈現出更明顯的季節現象．

圖4: 當β受季節因素的影響時，2006～2020年我國流腦流行的情況

圖5: 當a受季節因素的影響時，2006～2020年我國流腦流行的情況

為考慮εβ的靈敏度，令εa=0.3，分別取εβ=0.3,0.4,0.5，數值模擬結果如圖7所示，εβ越大，使得C(t),I(t)振蕩的幅度越大，呈現出更明顯的季節現象．取定εβ=0.4，讓εa變動，分別取εa=0.3,0.4,0.5，數值模擬結果如圖8所示，εa越大，I(t)振蕩的幅度越大，和圖7相比，εa的變動對染病者的數量影響較大．說明季節因素對疾病進展率a的影響大于疾病傳播速率β的影響，進一步驗證了季節因素對流腦的流行傳播影響觀點的第二種假設．說明我們國家在干旱的季節，天氣干燥、塵土飛揚可能損害鼻咽部位的粘膜屏障，使得無癥狀的攜帶者更容易轉變為染病者．因此，我們可以采取控制措施減少易感者轉化為攜帶者的比率進而控制攜帶者的數量，由于疫苗不能阻止攜帶者攜帶腦膜炎奈瑟菌，但它能刺激產生群體免疫進而減少疾病的爆發，可提高免疫覆蓋率、接種率，必要時在流腦高發期可考慮開展應急接種．

圖7: 當β和a同時受季節因素影響時，εβ的靈敏度

圖8: 當β和a同時受季節因素影響時，εa的靈敏度

7 結論

文中考慮疫苗因素對流腦傳播的影響，分析了一類帶有年齡結構的離散SCIRS模型的動力學性態以及在我國流腦中的應用．得到了模型的基本再生數，證明了無病平衡點的全局穩定性．對于地方病平衡點的穩定性，由于模型自身的復雜性，這里沒有給出證明．其次，我們把模型應用到我國流腦的傳播中，由于模型中很多參數是不確定的，對基本再生數作了關于參數的敏感性分析，得出有些參數對f R0的影響比較大，因此，我們可以控制這些參數的值來達到控制疾病的目的．最后，考慮了流腦發病的季節因素對模型加以改造，預測分析了我國流腦的發病情況，數值模擬的結果顯示季節因素對疾病進展率影響的程度大于對疾病傳染率的影響，這也為控制流腦在我國的流行傳播提供了建議．

參考文獻：

[1]World Health Organisation.Epidemics of meningococcal disease of African meningitis[J].Weekly Epidemiological Record,2001,37(2):281-288

[2]耿貫一.流行病學(第2版)[M].北京：人民衛生出版社，1996:618-639 Geng G Y.Epidemiology(2nd Edition)[M].Beijing:People’s Medical Publishing House,1996:618-639

[3]Hu X J.Epidemiological surveillance and prophylaxis of epidemic cerebrospinal meningitis[J].Chinese Academy of Preventive Medicine,2001,7(5):300-303

[4]胡緒敬.流行性腦脊髓膜炎流行的監測與預防[J].中國公共衛生，2004,20(5):638-640 Hu X J.Epidemiological surveillance and prophylaxis of meningococcal meningitis[J].Chinese Journal of Public Health,2004,20(5):638-640

[5]Li Y X,Luo H M,Li J H.The epidemiological characteristics of meningococcal meningitis of China in 2006-2007[J].Chinese Journal of Vaccines and Immunization,2008,14(3):234-237

[6]Thomas H L,Andrews N,Trotter C,et al.Meningococcal C booster not recommended by evidence[J].British Medical Journal,2008,337:a1139

[7]Maiden M C,Stuart J M.Carriage of serogroup C meningococci 1 year after meningococcal C conjugate polysaccharide vaccine[J].The Lancet,2002,359(9320):1829-1831

[8]Guzzetta G,Manfredi P,Gasparini R,et al.On the relationship between meningococcal transmission dynamics and disease:remarks on humoral immunity[J].Vaccine,2009,27(25):3429-3434

[9]Zhou Y,Ma Z.Global stability of a class of discrete age-structured SIS models with immigration[J].Mathematical Biosciences and Engineering,2009,6(2):409-425

[10]Elaydy S.An Introduction to Diff erence Equations[M].New York:Springer-Verlag,2004:204-255

[11]Diekmann O,Heesterbeek J A P,Metz J A J.On the defi nition and the computation of the basic reproduction ratio R0in models for infectious diseases in heterogeneous populations[J].Journal of Mathematical Biology,1990,28(4):365-382

[12]Allen L J S,van den Driessche P.The basic reproduction number in some discrete-time epidemic models[J].Journal of Diff erence Equations and Applications,2008,14(11):1127-1147

[13]Gil M I.Diff erence Equations in Normed Spaces:Stability and Oscillations[M].Amsterdam:Elsevier Science,2007

[14]National Bureau of Statistics of China.2006-2012 China Statistical Yearbook[M].Beijing:China Statistics Press,2012

[15]National Bureau of Statistics of China.China Population Census Data in 2010[M].Beijing:China Statistics Press,2010

[16]Chinese Center for Disease Control and Prevention.Yearly report on infectious disease in China[OL].http://www.chinacdc.cn/tjsj/fdcrbbg/index.html,2015-12-22

[17]Wang M,Li Y X,Li J H.The epidemiological characteristics of meningococcal meningitis of China in 2009-2010[J].Disease Surveillance,2012,27(6):435-438

[18]Irving T J,Blyuss K B,Colijn C,et al.Modelling meningococcal meningitis in the African meningitis belt[J].Epidemiology and Infection,2012,140(5):897-905

[19]Greenwood B.The changing face of meningococcal disease in West Africa[J].Epidemiology and Infection,2007,135(5):703-705

[20]Trotter C L,Edmunds W J.Reassessing the cost-eff ectiveness of meningococcal serogroup C conjugate(MCC)vaccines using a transmission dynamic model[J].Medical Decision Making,2006,26(1):38-47

[21]Zhou Y C,Khan K,Feng Z L,et al.Projection of tuberculosis incidence with increasing immigration trends[J].Journal of Theoretical Biology,2008,254(3):215-228

[22]Marino S,Hogue I B,Ray J C,et al.A methodology for performing global uncertainty and sensitivity analysis in systems biology[J].Journal of Theoretical Biology,2008,254(1):178-196

[23]Trotter C L,Greenwood B M.Meningococcal carriage in the African meningitis belt[J].Lancet Infectious Diseases,2007,7(12):797-803

[24]Moore P S.Meningococcal meningitis in sub-Saharan Africa:a model for the epidemic process[J].Clinical Infectious Diseases,1992,14(2):515-525

Received:18 Sep 2015. A ccep ted:29 July 2016.

Found ation item:The National Natural Science Foundation of China(11171267);the Technology Department Fund of Taiyuan Institute of Technology(2015LQ19);the International Development Research Center,Ottawa,Canada(104519-010).

Abstract:The dynamic characteristic of meningococcal meningitis SCIRS model with three age structures is studied.First,the basic reproduction numberis defi ned by using the regeneration matrix.It is proved that the disease-free equilibrium is globally asymptotically stable when<1.The disease-free equilibrium is unstable,there exists an endemic equilibrium and the system is uniformly persistent when>1.Second,using the data from the report of notifi able infectious diseases in China,the model is applied to describe the spread of meningococcal meningitis in China.For the uncertainty of many parameters in the model,we make sensitivity analysis about the parameters of the basic reproductive number.Finally,we consider the infl uence of seasonal factors to the incidence of meningococcal meningitis to modify the model,and predict the meningococcal meningitis in population development trend of China.The numerical simulation results show that the impact of seasonal factors on the rate of disease progression of meningococcal meningitis is greater than the eff ect on infection rate,this also provides advice to control the spread of meningococcal meningitis in our country.

K eyword s:discrete model;age structure;meningococcal meningitis;basic reproduction number;sensitivity analysis