廖 書， 楊煒明

(重慶工商大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，重慶 400067)

越來越多的研究者利用數(shù)學(xué)模型來構(gòu)造傳染病模型以對其進(jìn)行定量和定性的研究，參見文獻(xiàn)[1-6]。為了進(jìn)一步控制水源性傳染病的流行，以霍亂傳染病為代表，首先考慮霍亂弧菌在不潔水源中會存活一段較長的時間，即增加時滯的考量，其次預(yù)防接種和利用消毒劑消毒不潔水源都是控制霍亂傳播的有效手段，建立并研究一個同時含有預(yù)防接種和消毒不潔水源雙重控制策略的霍亂時滯模型，最后采用2008年津巴布韋霍亂的數(shù)據(jù)進(jìn)行數(shù)值模擬以研究這兩種控制方法的有效性以及時滯大小對模型穩(wěn)定性的影響。

1 霍亂時滯模型

設(shè)總?cè)藬?shù)N=S+I+V+R，且總?cè)丝跀?shù)與飲水的供給和成正比，S，I，V和R分別表示易感染者、染病者、接種疫苗者和移出者，W為霍亂病菌濃度，T為消毒劑在水源中的濃度。消毒劑可以很有效地殺死水源中的霍亂弧菌并且控制其傳播，假設(shè)消毒劑的濃度與失效率成正比，也與水源中霍亂病菌濃度成正比。再假設(shè)消毒劑的自然喪失率與其濃度成比例，且消毒劑攝取了霍亂弧菌后，其吸收率是與病毒的密度以及消毒劑的濃度成比例的。模型中其他的參數(shù)βI和βW分別表示環(huán)境與人之間傳播和人與人之間傳播的傳染率系數(shù)，μ1和μ分別表示感染者和非感染者不同的死亡率，φ為疫苗接種率，α為消毒劑的失效率，π為霍亂病菌的增長率，ξ為霍亂病菌的自然喪失率，d為作為安全用水提供給住戶的水源中霍亂病菌的喪失率，χ為當(dāng)使用了消毒劑后霍亂病菌濃度的喪失率，θ為使用消毒劑的濃度，η為消毒劑在水源中的有效吸收率，γ為染病者的復(fù)原率，τ表示病菌在不潔水源中的潛伏期。σ表示疫苗的有效率，當(dāng)σ=0為該疫苗完全有效，σ=1意味著疫苗沒有效果。所有的參數(shù)都為正數(shù)。

(1)

(2)

(3)

(4)

(5)

(6)

模型的初始條件如下:

S≥0，I≥0，V≥0，W≥0，T≥0，R≥0

注意到系統(tǒng)方程中R的獨立性(在式(1)—式(5)中均不含有R)，所以為了簡化計算，在后面的分析中只考慮式(1)—式(5)即可。

2 地方病平衡點

本節(jié)主要研究地方病平衡點X*=(S*，I*，V*，W*，T*)的穩(wěn)定性。首先令s=S-S*，v=V-V*，i=I-I*，w=W-W*以及f=T-T*，其中s，i，v，w和f是圍繞X*的微小擾動。令

在地方病平衡點的特征方程可計算得：

(7)

其中

以及

由Routh-Hurwize判別法，可以證明當(dāng)τ=0時，地方病平衡點是局部漸近穩(wěn)定的。因此處計算比較冗長，故省去計算， 直接得出以下定理:

定理1 當(dāng)R0>0，τ=0時，系統(tǒng)模型式(1)—式(5)的地方病平衡點是局部漸近穩(wěn)定的。

當(dāng)τ>0時，設(shè)λ=iω是方程式(7)的一個根， 將λ代入方程式(7)并分離實部和虛部后，可得到下面的兩個方程：

(8)

(9)

將方程式(8)和方程式(9)平方相加：

(10)

再令ω2=x得到:

(11)

如果系數(shù)Ci滿足Routh-Hurwitz條件，則方程式(11)不會有任何正根，因此不會得到滿足方程式(8)和方程式(9)的正ω。在這種情況下有以下定理：

定理2 當(dāng)R0>0，τ>0時，如果Routh-Hurwitz條件滿足，系統(tǒng)式(1)—式(5) 的地方病平衡點是局部漸近穩(wěn)定的。

(12)

對方程式(7)左右兩邊同時求λ關(guān)于τ的導(dǎo)數(shù)得到：

(13)

化簡可得:

(14)

因此可以計算求得:

(15)

又由前面的假設(shè)條件C0<0可知，F(xiàn)'(x)>0，則上述方程一定大于0。即意味著當(dāng)τ>τ0時，至少存在一個根有一個正實部并且從左向右穿過虛軸。因此當(dāng)τ=τ0時，Hopf分支產(chǎn)生，并在τ=τ0附近產(chǎn)生一簇周期解。由文獻(xiàn)[7，8]的Hopf分支定理，可以得到下面的定理：

定理3 當(dāng)R0>0，時滯τ∈[0，τ0]時，如果Routh-Hurwitz條件滿足，系統(tǒng)式(1)—式(5)的地方病平衡點是局部漸近穩(wěn)定的；當(dāng)時滯τ>τ0時，系統(tǒng)式(1)—式(5)不穩(wěn)定。當(dāng)τ=τ0時，系統(tǒng)在地方病平衡點產(chǎn)生Hopf分支，并在τ=τ0附近產(chǎn)生一簇周期解。

3 模型模擬

采用世界衛(wèi)生組織發(fā)布的2008年—2009年津巴布韋霍亂的數(shù)據(jù)進(jìn)行模型的數(shù)值模擬。津巴布韋的總?cè)丝诖蠹s12 347 240人， 為了計算方便，本節(jié)按比例減少1 200倍系數(shù)使得其總?cè)丝跒? 000，模型中所用到的其他參數(shù)值為

N=10 000，μ=0.000 442，π=70
ξ=0.233 3，d=0.000 014，χ=0.99
α=1，θ=0.001，γ=1.4，σ=0.2
φ=0.003，βI=0.00 012，ψ=0.05
βW=0.000 000 42，η=0.000 4

同時初始條件在相應(yīng)降低1 200倍之后為S0=9 890，I0=10，V0=100，T0=10和W0=0。

3.1 當(dāng)τ=0時

把所有的參數(shù)代入R0的表達(dá)式，可求得R0=1.7，即在這種情況下，霍亂會在津巴布韋流行開來，并且最終地方病平衡點的值也可求得為(S*，I*，V*，R*) = (4 938.65， 0.442， 2 960.5， 2 100.2)，這與津巴布韋的霍亂傳播實際情況相符合。圖1和圖2表示S*，I*，V*和R*隨時間的變化趨勢。

從圖1可以觀察到染病者人數(shù)在27周的時候達(dá)到最高峰72(除以1 200倍系數(shù)后的值)，然后直接下降到接近0，說明在進(jìn)行有效的防疫之后，此次疫情已消除?？紤]長期的情況，從圖2可知，在第一次疫情爆發(fā)結(jié)束之后，還會有若干次疫情再度爆發(fā)，但最高峰值會越來越小，直到約384年后S*，V*和R*最終達(dá)到自己穩(wěn)定值分別為4 938.65、2 960.5和2 100.2，此時在津巴布韋的霍亂才最終徹底消除。

圖1 當(dāng)τ=0時，隨著時間的變化I的變化趨勢

圖2 當(dāng)τ=0時，隨著時間的變化S，V，R的變化趨勢

圖3和圖4分別表示在沒有預(yù)防接種和沒有進(jìn)行水源消毒的情況下染病者人數(shù)隨時間變化的趨勢?？梢悦黠@看出，在缺失預(yù)防接種或者水源消毒后，染病者人數(shù)明顯高于圖1中的染病者人數(shù)，說明這兩種預(yù)防控制措施非常有效，都可以極大降低霍亂傳播。

圖3 缺失預(yù)防接種時，隨著時間的變化I的變化趨勢

圖4 缺失水源消毒時， 隨著時間變化I的變化趨勢

3.2 當(dāng)τ≠0時

圖5 當(dāng)τ=2時，隨著時間的變化I的變化趨勢

圖6 當(dāng)τ=4時，隨著時間的變化I的變化趨勢

在本節(jié)采用和上節(jié)同樣的參數(shù)值和初始條件對τ≠0的情況進(jìn)行數(shù)值模擬，此時時滯臨界值τ0可求得為2.88。從圖5可知，當(dāng)取τ=2時，I*的振動逐漸降低并最終趨于它們分別的穩(wěn)定值。然后當(dāng)時滯大于闕值時，系統(tǒng)的穩(wěn)定性會發(fā)生變化。

圖6中，時滯取值為4時，可以看出I*開始發(fā)生不穩(wěn)定的變化。說明對傳染病進(jìn)行模擬和預(yù)測時，在時滯變大的情況下，對未來的預(yù)測會變得更加困難。

4 結(jié) 論

建立并分析了一個同時含有預(yù)防接種和水源消毒的霍亂時滯模型。研究了模型地方病平衡點的穩(wěn)定性，并通過分析相應(yīng)特征方程根的分布，得出當(dāng)時滯大小超過一個闕值的時候，穩(wěn)定性發(fā)生變化，產(chǎn)生了Hopf分支，系統(tǒng)產(chǎn)生波動。結(jié)果說明了時滯模型比一般的ODE模型具有更大的現(xiàn)實意義，但計算上也更加困難。沒有對地方病平衡點的全局穩(wěn)定性進(jìn)行證明，這將是未來工作之一。最后要注意盡管大劑量地進(jìn)行水源殺毒可以在短時間內(nèi)有效控制霍亂疫情的流行，但在現(xiàn)實生活中，大劑量的殺毒劑使用又會對人體健康造成巨大危害，因此要同時考量對人體健康和安全的因素，合理控制殺毒劑的使用。

[1] CODECO C T. Endemic and Epidemic Dynamics of Cholera: the Role of the Aquatic Reservoir[J]. BMC Infectious Diseases， 2001(1):1

[2] HARTLEY D M， MORRIS J G，SMITH D L. Hyperinfectivity:a Critical Element in the Ability of V Cholerae to Cause Epidemics[J]. PLoS Medicine， 2006(3):63-69

[3] JOH R I， WANG H， WEISS H，et al. Dynamics of Indirectly

Transmitted Infectious Diseases with Immunological Threshold[J]. Bulletin of Mathematical Biology， 2009， 71(4):845-862

[4] SANCHES R P， FERREIRA C P，KRAENKEL R A. The Role of Immunity and Seasonality in Cholera Epidemic[J]. Bulletin of Mathematical Biology， 2011，73(12):2916-2931

[5] LIAO S，WANG J. Stability Analysis and Application of a Mathematical Cholera Model[J]. Mathematical Biosciences and Engineering， 2011(8): 733-752

[6] TIEN J H， EARN D J D. Multiple Transmission Pathways and Disease Dynamics in a Waterborne Pathogen Model[J].Bulletin of Mathematical Biology，2010，72: 1506-1533

[7] GOPALSAMY K.Stability and Oscillations in Delay Differential Equations of Population Dynamics[M].Dordrecht：Kluwer Academic，1992

[8] HALE J K.Theory of Functional Differential Equations[M].New York: Spring-Verlag，1977