具有季節(jié)性和接種疫苗的元胞自動機(jī)傳染病模型

王燕 劉茂省 李有文

摘要：為了研究基于二維元胞自動機(jī)的一類具有季節(jié)性和疫苗接種的SIR傳染病模型對控制傳染病的影響，利用平均場近似方法建立了非線性離散模型。對模型進(jìn)行了數(shù)學(xué)分析，通過計算無病平衡點(diǎn)的Jacobian矩陣的譜半徑，得到了無病平衡點(diǎn)的局部穩(wěn)定性，通過數(shù)值模擬，找出了在正平衡點(diǎn)處染病者與元胞鄰居結(jié)構(gòu)δ的關(guān)系，并且考慮了不同的初始患者設(shè)置對疾病傳播速度的影響，不同的接種疫苗率對疾病傳播速度的影響。通過研究可知：當(dāng)疾病的初始患者為中心分布時，傳染病的傳播速度小于初始患者為隨機(jī)分布時的速度;隨著接種疫苗率的增大，傳染病的傳播速度會減小，對初始患者為中心分布的元胞自動機(jī)影響較大。研究模型在控制傳染病傳播的動力學(xué)研究方面具有一定的參考價值。

關(guān)鍵詞：微分動力學(xué)系統(tǒng);元胞自動機(jī);SIR模型;季節(jié)性;接種疫苗;平均場近似

中圖分類號：O29 文獻(xiàn)標(biāo)志碼：A doi： 10.7535/hbgykj.2018yx02002

WANG Yan，LIU Maoxing，LI Youwen.A cellular automata model for epidemics with seasonal and vaccination[J].Hebei Journal of Industrial Science and Technology，2018，35（2）：84-90.A cellular automata model for epidemics with seasonal and vaccination

WANG Yan， LIU Maoxing， LI Youwen

（School of Science， North University of China， Taiyuan， Shanxi? 030051， China）

Abstract：In order to study the influence on the infectious desease control of a class of SIR epidemic models with seasonal and vaccination based on two-dimensional cellular automaton， a nonlinear discrete model is established by means of mean-field approximation. The model is analyzed mathematically. By calculating the spectral radius of the Jacobian matrix with the disease-free equilibrium， the local stability of the disease-free equilibrium is obtained. The relationship between the number of infected persons at the positive equilibrium point and the structure of the cell neighborhood is obtained by numerical simulating. And the effects of different initial patient settings on the disease transmission speed and the different vaccination rates on disease spread are considered. Through numerical simulation， the conclusion is gotten that when the initial patient is distributed at a central distribution， the speed of transmission is less than the velocity of the infectious disease when the initial patient is randomly distributed; with the increase of vaccination rate， the infectious disease transmission speed will decrease， which has greater impact on the initial distribution of patients as the center of cellular automaton. The study model has reference value in the dynamic research of controlling infectious disease transmission.

Keywords：differential dynamical systems;cellular automata;SIR model;seasonality;vaccination;mean-field approximation

傳染病是由各種病原體引起的能在人與人、動物與動物或人與動物之間相互傳播的一類疾病，是一種可以從一個人或其他物種，經(jīng)過各種途徑傳染給另一個人或物種的傳染病。自20世紀(jì)40年代開始出現(xiàn)以微分方程為主的傳染病模型到現(xiàn)在，有大量的傳染病模型建立在常微分方程的基礎(chǔ)上，其中很有影響的是SIS[1]和SIR[2]模型，但是這種方法忽略了傳播過程的局部特征，而元胞自動機(jī)可以克服這個缺點(diǎn)，已經(jīng)成為研究流行病學(xué)的一種新的研究方法[3]。

在生物學(xué)中，元胞自動機(jī)的設(shè)計思想來源于生物學(xué)自繁殖的思想，因而它在生物學(xué)上的應(yīng)用更為自然而廣泛。例如元胞自動機(jī)用于腫瘤細(xì)胞的增長機(jī)理和過程模擬，人類大腦的機(jī)理探索[4]，艾滋病病毒HIV的感染過程[5]，自組織、自繁殖等生命現(xiàn)象的研究。元胞自動機(jī)是一個在時間和空間上都離散，通過局部元胞的相互作用而引起全局變化的動力系統(tǒng)。散布在規(guī)則格網(wǎng)格中的每一個元胞取有限離散狀態(tài)，遵循同樣的作用規(guī)則，依據(jù)確定的局部規(guī)則作同步更新。大量元胞通過簡單的相互作用而構(gòu)成動態(tài)系統(tǒng)的演化。AHMED等[6-7]定義了一類基于元胞自動機(jī)的傳染病模型，討論了傳播強(qiáng)度的影響及傳播結(jié)果的分類。 SIRAKOULIS等[8]、FUENTES等[9]利用元胞自動機(jī)分別建立了SIS和SIR模型，研究了人群移動、接種疫苗以及人口初始分布等因素對疾病傳播的影響。LANGLAIS等[10]通過元胞自動機(jī)研究了SEIR傳染病模型。譚欣欣等[11]考察了元胞自動機(jī)上SEIR模型的種群移動比例和種群移動最大距離對傳染病傳播的影響。暢春玲等[12]考慮了具有垂直傳播的傳染病模型。楊嬋等[13]研究了異質(zhì)網(wǎng)絡(luò)下含有相關(guān)系數(shù)的SIR傳染病模型。但是有一類受季節(jié)性影響的傳染病卻少有人研究，如手足口病、肺結(jié)核病等。本文將研究基于二維元胞自動機(jī)的一類具有季節(jié)性和接種疫苗的SIR傳染病模型。

1元胞自動機(jī)理論

元胞自動機(jī)最基本的組成如圖1所示，可分為元胞、元胞空間、鄰居及演化規(guī)則4部分[14]，可用一個四元組表示：CA={Ld，S，N，F(xiàn)}，其中Ld代表元胞空間，d代表空間維數(shù)，S為元胞的有限狀態(tài)集，N表示一個所有鄰居內(nèi)元胞的組合，包括中心元胞在內(nèi)，是包含n個不同元胞的空間矢量，n為鄰居元胞個數(shù)，F(xiàn)表示局部轉(zhuǎn)換函數(shù)，也就是演化規(guī)則。元胞狀態(tài)的變化依賴于自身的狀態(tài)和鄰居的狀態(tài)，且某個元胞下一時刻的狀態(tài)只決定于鄰居的狀態(tài)以及自身的初始狀態(tài)。

以常用的規(guī)則四邊形網(wǎng)格劃分為例，二維元胞自動機(jī)的常用鄰居結(jié)構(gòu)可以分為Von Neumann型和Moore型2種。Von Neumann型的元胞鄰居數(shù)目為2d，Moore型的鄰居數(shù)目為3d-1，如圖3所示。本研究采用半徑為1的Moore型鄰居結(jié)構(gòu)。關(guān)于元胞自動機(jī)的具體基本知識可參考文獻(xiàn)[15]。

2模型

筆者將人群分為3類：易感者（S）、染病者（I）和恢復(fù)者（R），假設(shè)疾病的傳播是通過染病者和易感者直接接觸傳染的，在模型中考慮個體出生、自然死亡、因病死亡以及易感者接種疫苗。傳染病的傳播過程如圖4所示。

2.1元胞自動機(jī)模型

把這些個體看作是分布在具有J個點(diǎn)的二維環(huán)空間上，或者是在一個平面上彼此有聯(lián)系的J個點(diǎn)。筆者將環(huán)和復(fù)雜網(wǎng)絡(luò)或隨機(jī)網(wǎng)絡(luò)拓?fù)溆成錇橐粋€包含J節(jié)點(diǎn)的有限平面L2，這樣的結(jié)構(gòu)可以用來描述涉及個體的系統(tǒng)[16]。因此，在這個系統(tǒng)中L2上的每個節(jié)點(diǎn)最多被一個個體占據(jù)或空著。

筆者用Z（t）={Z1，1（t），Z1，2（t），Z2，1（t），…，Zi，j（t）}，i，j∈N，表示t時刻系統(tǒng)的狀態(tài)結(jié)構(gòu)，其中Zi，j（t）表示t時刻元胞（i，j）的狀態(tài)。系統(tǒng)一共有4種狀態(tài)，分別用-1，0，1，2來表示不同的狀態(tài)，即Zi，j（t）={-1，0，1，2}，其中-1代表元胞狀態(tài)為空，0代表元胞狀態(tài)為易感者，1代表元胞狀態(tài)為染病者，2代表元胞狀態(tài)為恢復(fù)者。例如，Zi，j（t）=1表示的是t時刻（i，j）的狀態(tài)為染病者。如果某一元胞在t時刻消亡（自然死亡或因病死亡），則該元胞在t+1時刻狀態(tài)變?yōu)榭铡Ｔ谶@個系統(tǒng)中，每個元胞都會根據(jù)其鄰居的狀態(tài)以一定的概率轉(zhuǎn)變成其他狀態(tài)，空元胞會以b轉(zhuǎn)化為易感者狀態(tài)，被占據(jù)的元胞會以ε轉(zhuǎn)變?yōu)榭赵?，染病者會以μ轉(zhuǎn)變?yōu)榛謴?fù)者，或以ε+d轉(zhuǎn)變?yōu)榭赵?。?dāng)然，每個元胞也會保持自身的狀態(tài)。筆者采用半徑為1的Moore型鄰居結(jié)構(gòu)，用ni，j（Z，t）表示t時刻元胞（i，j）狀態(tài)為Z的鄰居數(shù)量，其中Z∈{-1，0，1，2}元胞的演化規(guī)則[17]由表1給出。

2.2平均場近似理論模型

平均場近似理論是研究空間問題應(yīng)用最廣泛和最簡單的方法，它認(rèn)為格子中的每個個體的狀態(tài)都是獨(dú)立存在的。HIEBELER[18]研究了元胞自動機(jī)上平均場近似的局部擴(kuò)散和無限擴(kuò)散，并且模擬了可以用于離散空間模型的一整套方法的2個極端。JIN等[19]研究了基于元胞自動機(jī)的平均場近似，同時考慮空間異質(zhì)易感性的SIR模型。

筆者用x（t）表示t時刻易感者的全局密度，y（t）表示t時刻染病者的全局密度，其中x（t），y（t）∈[0，1]，則有x（t）=N（0，t）J，y（t）=N（1，t）J。用H（t）表示系統(tǒng)的全局密度，則H（t）=N（t）J，H（t）<1，所以有1-H（t）為空元胞的全局密度。其中N（t）為t時刻L2上全部個體的數(shù)量，N（t）=∑z=-1，0，1，2N（Z，t），α（t）為有效接觸概率，則1-α（t）為沒有被傳染和保持易感者狀態(tài)的概率，即易感者不被鄰居感染的概率為q=（1-α（t））n（1，t），n（1，t）表示t時刻易感者鄰居中為染病者的數(shù)目。用δ表示鄰居中被占據(jù)的數(shù)目。因此易感者被鄰居傳染為染病者的概率為q′=1-q=1-（1-α（t））n（1，t）。在系統(tǒng)中，t+1時刻的易感者數(shù)量是t時刻剩余的易感者數(shù)量，染病者及恢復(fù)者的數(shù)量都是如此。由此得到一個非線性離散方程：

x（t+1）=b（1-x（t）-y（t）-z（t））+x（t）（1-K）-px（t）-εx（t），y（t+1）=（1-d-μ-ε）y（t）+x（t）K，z（t+1）=z（t）+px（t）+μy（t）-εz（t），（1）

其中：

K=∑δi=0δiy（t）i（1-y（t））δ-i×

[1-（1-α（t））i]，

表示易感者被鄰居感染的概率的期望值，α（t）=β（1+ξ cos πt），其中，β是平均接觸率;ξ是接觸率的波動幅度，把α（t）代入方程（1）可得到以下模型：

x（t+1）=b（1-x（t）-y（t）-z（t））+x（t）（1-y（t）α（t））δ-px（t）-εx（t），y（t+1）=（1-d-μ-ε）y（t）+x（t）[1-（1-y（t）α（t））δ]，z（t+1）=z（t）+px（t）+μy（t）-εz（t）。（2）

令方程（2）中x（t+1）=x（t）=x，y（t+1）=y（t）=y，z（t+1）=z（t）=z，得到無病平衡點(diǎn)滿足的方程：

x=b（1-x-y-z）+x（1-yα（t））δ-px-εx，y=（1-d-μ-ε）y+x[1-（1-yα（t））δ]，z=z+px+μy-εz，（3）

即可得到無病平衡點(diǎn)：

E0=（bε（b+ε）（p+ε），0，pb（b+ε）（p+ε））。

對方程（3）進(jìn)行泰勒展開，只保留前2項(xiàng)，則有：

x=b（1-x-y-z）+x（1-δyα（t））-px-εx，y=（1-d-μ-ε）y+xδyα（t），z=z+px+μy-εz，（4）

由方程（4）解得：

x*=γδα（t），y*=bεδα（t）-γ（p+ε）（b+ε）δα（t）[b（μ+ε）+εγ]，

z*=pγδεα（t）+μy*ε，其中γ=ε+μ+d，當(dāng)δ>γ（p+ε）（b+ε）bεα（t）時，模型存在正平衡點(diǎn)E*=（x*，y*，z*）。

3模型的數(shù)學(xué)分析

為了分析無病平衡點(diǎn)在平均場近似模型中的穩(wěn)定性，將方程（2）線性化，則有：

x（t+1）=b（1-x（t）-y（t）-z（t））+x（t）（1-δy（t）α（t））-px（t）-εx（t），y（t+1）=（1-d-μ-ε）y（t）+x（t）δy（t）α（t），z（t+1）=z（t）+px（t）+μy（t）-εz（t），（5）

得到無病平衡點(diǎn)E0處的Jacobian矩陣為

JE0=1-b-ε-p-b（1+δεα（t）（p+ε）（b+ε）-b01-γ-bδεα（t）（p+ε）（b+ε）0pμ1-ε。

計算3階JE0在給定參數(shù)值時的譜半徑。如表2所示，當(dāng)δ<4時，ρ（JE0）<1，即此時無病平衡點(diǎn)是局部漸近穩(wěn)定的;當(dāng)δ>4時，ρ（JE0）>1， 則無病平衡點(diǎn)是不穩(wěn)定的。其中b=0.04，ε=0.001，μ=0.12，p=0.004，β=0.2，ξ=0.28。

4數(shù)值模擬與結(jié)果

筆者利用Python軟件對系統(tǒng)進(jìn)行了數(shù)值模擬，設(shè)置了不同的初始患者狀態(tài)，考察其對傳染病傳播速度的影響。

從圖5可以看出，在正平衡點(diǎn)處染病者隨著元胞鄰居結(jié)構(gòu)δ的增大而增加。圖6表示方程（2）在δ=4，ρ（JE0）=0.998 2時的模擬結(jié)果。圖7至圖13顯示了傳染病在時間和空間都離散的元胞自動機(jī)上的模擬。假設(shè)在200×200的規(guī)則方格中，采取半徑為1的Moore型鄰居結(jié)構(gòu)，在t=0時，即初始化狀態(tài)，取1/10的格子令狀態(tài)為空，100個格子狀態(tài)為染病者，分為2種情況：1）初始患者在元胞自動機(jī)上為中心分布。2）初始患者在元胞自動機(jī)上為隨機(jī)分布。其中β=0.2，ξ=0.28。圖7和圖8分別為初始患者為中心分布和隨機(jī)分布時的元胞自動機(jī)模擬結(jié)果，可以看到初始患者為中心分布的傳播速度要小于初始患者為隨機(jī)分布的傳播速度。圖9表示不同的2種初始患者狀態(tài)在元胞自動機(jī)上每個時刻易感者、染病者、恢復(fù)者的數(shù)量變化。圖10和圖11表示當(dāng)p=0.006時，不同的初始患者分布的元胞自動機(jī)模擬結(jié)果，與圖7和圖8比較，可以看出不同的p值對傳染病傳播速度的影響，當(dāng)p值增大時，元胞自動機(jī)上疾病的傳播速度會變緩。圖12和圖13反映了不同的p值對初始患者為中心分布的傳播速度的影響要大于隨機(jī)分布的影響。接種疫苗對季節(jié)性傳染病可以起到一定的減緩作用，這和一般的常微分方程模型的結(jié)果一樣，所以說用元胞自動機(jī)模擬傳染病的傳播是可行的，可以為研究傳染病提供更多的方法。

本研究基于二維元胞自動機(jī)，考慮了具有季節(jié)性傳染和接種疫苗的SIR模型，考察了影響傳染病傳播速度的原因，影響因素很多，模型有一定的局限性。后期可以研究當(dāng)染病者分為有癥狀和無癥狀的兩種情況，對有癥狀的染病者采取隔離措施，進(jìn)一步分析傳染病的傳播情況。

參考文獻(xiàn)/References：

[1]BAILEY N T J. The mathematical theory of infectious diseases and its applications[J]. Immmunology， 1977，28（2）：479-480.

[2]ANDERSON R M， MAY R M. Infectious Diseases of Humans： Dynamics and Control[M]. England：Cambridge University Press， 1991.

[3]CHOPARD B. Cellular automata modeling of physical systems[J]. Computational Complexity， 1998：865-892.

[4]VICTOR J D. What can automaton theory tell us about the brain？[J]. Physica D： Nonlinear Phenomena， 1990， 45（1/2/3）：205-207.

[5]SIEBURG H B， MCCUTCHAN J A， CLAY O K， et al. Simulation of HIV infection in artificial immune system[J]. Physica D：Nonlinear Phenomena，1990，45（1/2/3）：208-227.

[6]AHMED E， AGIZA H N， HASSAN S Z. On modeling hepatitis B transmission using cellular automata[J]. Journal of Statistical Physics， 1998， 92（3/4）：707-712.

[7]AHMED E， AGIZA H N. On modeling epidemics including latency， incubation and variable susceptibility[J]. Physica A Statistical Mechanics & Its Applications， 1998， 253（1/2/3/4）：347-352.

[8]SIRAKOULIS G C， KARAFYLLIDIS I， THANAILAKIS A. A cellular automaton model for the effects of population movement and vaccination on epidemic propagation[J]. Ecological Modelling， 2000， 133（3）：209-223.

[9]FUENTES M A， KUPERMAN M N. Cellular automata and epidemiological models with spatial dependence[J]. Physica A Statistical Mechanics & Its Applications， 1999， 267（3/4）：471-486.

[10]LANGLAIS M， SUPPO C. A remark on a generic SEIRS model and application to cat retroviruses and fox rabies[J].Mathematical and Computer Modeling，2000，31（4/5）：117-124.

[11]譚欣欣， 戴欽武， 史鵬燕，等. 基于元胞自動機(jī)的個體移動異質(zhì)性傳染病傳播模型[J]. 大連理工大學(xué)學(xué)報， 2013，53（6）：908-914.

TAN Xinxin，DAI Qinwu，SHI Pengyan，et.al.CA-based epidemic propagation model with inhomogeneity and mobility[J]. Journal of Dalian University of Technology， 2013，53（6）：908-914.

[12]暢春玲， 隋英， 孫艷玲. 具有垂直傳播和一般接觸的元胞自動機(jī)傳染病模型[J]. 科技廣場， 2010（6）：132-134.

CHANG Chunling，SUI Ying，SUN Yanling.Cellular automata model for epidemics with the characteristic of? vertical tarnsmission and contact[J].Science Mosaic，2010（6）：132-134.

[13]楊嬋，張菊平.異質(zhì)網(wǎng)絡(luò)下含有相關(guān)系數(shù)的SIR傳染病模型的建立和分析[J].河北工業(yè)科技，2018，35（1）：7-11.

YANG Chan，ZHANG Juping. Establishment and analysis of SIR epidemic model with correlation coefficient under heterogeneous network[J]. Hebei Journal of Industrial Science and Technology， 2018，35（1）：7-11.

[14]陜振沛， 寧寶權(quán)， 郭亞丹. 采用元胞自動機(jī)的甲型H1N1流感傳播模型研究[J]. 六盤水師范學(xué)院學(xué)報， 2013，25（2）：73-76.

SHAN Zhenpei，NING Baoquan，GUO Yadan.Transmission of influenza a （H1N1） model of cellular automata adopted for research[J].Journal of Liupanshui Normal University，2013，25（2）：73-76.

[15]YACOUBI S E， JAI A E. Cellular automata modelling and spreadability[J]. Mathematical & Computer Modelling， 2002， 36（9）：1059-1074.

[16]LI Ying， LIU Yang， SHAN Xiuming，et al. Dynamic properties of epidemic spreading on finite size complex networks[J]. Chinese Physics， 2005， 14（11）：2153-2157.

[17]LIU Quanxing， JIN Zhen. Cellular automata modelling of SEIRS[J]. Chinese Physics， 2005， 14（7）：1370-1377.

[18]HIEBELER D. Stochastic spatial models： From simulations to mean field and local structure approximations[J]. Journal of Theoretical Biology， 1997， 187（3）：307-319.

[19]JIN Zhen， LIU Quanxing. A cellular automata model of epidemics of a heterogeneous susceptibility[J]. Chinese Physics， 2006， 15（6）：1248-1256.