一類具有Allee效應的傳染病模型時空斑圖分析

李桂花，李高峰

（1.中北大學 理學院，山西 太原 030051；2.西南大學 數(shù)學與統(tǒng)計學院，重慶 400715；3.新疆農(nóng)二師庫爾勒醫(yī)院，新疆 庫爾勒 841000）

0 引言

寄生蟲依靠宿主而生活，它對宿主的影響是多方面的，它或者攝取宿主的營養(yǎng)，對附著組織產(chǎn)生損害，從而影響宿主的生長；或者改變宿主的行為，增加發(fā)病率，使宿主易于被捕食者捕食；或者使宿主不能獲得充分的生存資源.在種群水平上就是對宿主種群的內(nèi)稟增長率產(chǎn)生影響，這種影響與寄生蟲的感染豐度和頻率分布密切相關(guān).

寄生蟲從一個宿主傳播到另外一個宿主通常是通過自由生活的感染期而取得的.在這種情況下，潛在的宿主和感染期的寄生蟲遭遇的機會無疑會受到宿主和寄生蟲密度以及它們的空間分布的影 響［1-2］.R.M.Anderson，P.J.Whitefied 和A.P.Dobson［3］研究了寄生在體表的復殖吸蟲感染期密度對該寄生蟲的傳播動態(tài)的影響.R.M.Anderson［4］研究了穿透表皮的寄生蟲感染期的密度對傳播動態(tài)的影響.A.E.Keymer，R.M.Anderson［5］和A.E.Keymer［6］研究了通過消化道而獲得感染的寄生蟲的傳播過程，其研究表明感染期的寄生蟲的密度和空間分布、宿主密度、宿主和寄生蟲接觸的長短等因素都影響著寄生蟲的傳播過程.

A.E.Keymer［7］在實驗條件下證實了Hymenolepis Diminuta 的囊尾蚴有調(diào)節(jié)Tribolium Confusum 種群密度的作用.R.M.Anderson和J.Crombie［8］證實了曼氏血吸蟲Schistosoma Mansoni的幼蟲階段具有調(diào)節(jié)中間宿主種群的作用.R.M.Anderson，P.J.Whitefield和C.A.Mills［9］在實驗室條件下研究了外寄生蟲吸蟲的種群動態(tài)，對于各種種群參數(shù)都進行了詳細的研究，其中大量的種群過程是密度制約的.

R.M.Anderson在文獻［10］中提出一個數(shù)學模型，考慮宿主染病是與自由寄生蟲v的接觸有關(guān)，這里易感者宿主的疾病發(fā)生率是雙線性的，但是許多實驗表明這個疾病的發(fā)生率是與寄生蟲的劑量有關(guān)的，且常常是具有S形的非線性函數(shù)，于是R.R.Regoes等［11］建立了下面的模型：

式中：A0表示易感者宿主的輸入量；d表示宿主的自然死亡率；β為染病的比例系數(shù)；ε為因寄生蟲引起的死亡率；c為染病宿主體內(nèi)自由病毒的釋放率；u為寄生蟲的死亡率；g（v）是寄生蟲濃度v的S形函數(shù)，且

這里β＝α＝1／mk，m表示染病的劑量；k為S形曲線在點m處的傾角.表明系統(tǒng)（2）在閾值條件下存在兩個正平衡點，且一個是不穩(wěn)定的，另一個是穩(wěn)定的.當k連續(xù)變化時，Li Guihua等［12］發(fā)現(xiàn)會有更有趣的更復雜的性態(tài)發(fā)生（見文獻［12］）.假設寄生蟲的動力學行為比染病者宿主的動力學行為快的多，即u?δ.令V＝uv，則模型（2）為

本文將考慮k＝1和k＝2時空間效應對模型（3）性態(tài)的影響，模型如下：

式中：D1，D2分別為易感者宿主S與染病者宿主I的擴散系數(shù).

1 線性穩(wěn)定性分析

為了簡便，將系統(tǒng)（4）作無量綱變換，令

并用t代替τ，得到系統(tǒng)（4）的等價系統(tǒng)

這里：

首先考慮最簡單的形式，即k＝1時系統(tǒng)（5）的動力學性態(tài).當d1＝d2＝0時，系統(tǒng)（5）的正平衡點存在性及穩(wěn)定性如下：

定理1 若A＞c，則系統(tǒng)（5）存在惟一正均勻定態(tài)解

由文獻［12］可知，正平衡點只要存在，就一定是全局漸近穩(wěn)定的.為了分析圖靈不穩(wěn)定的條件，首先計算微擾方程的系數(shù)

由于a11，a22始終是小于0，很顯然a11d2+a22d1＜0.這樣由特征方程很容易知道，系統(tǒng)的特征值始終為負，即系統(tǒng)（5）的正均勻定態(tài)解是穩(wěn)定的，不存在圖靈分岔.

接著考慮k＝2時系統(tǒng)（5）的動力學性態(tài).當d1＝d2＝0時，系統(tǒng)（5）的正平衡點存在性及穩(wěn)定性如下：

定理2（a）若A2＝4c（1+cm），則系統(tǒng)（5）有惟一正均勻定態(tài)解，E＊＝（x＊，y＊）.

（b）若A2＞4c（1+cm），則系統(tǒng)（5）存在兩個正均勻定態(tài)解，E1＝（x1，y1），E2＝（x2，y2），且y1＜y＊＜y2.

這里

很容易計算當系統(tǒng)存在惟一的均勻定態(tài)解時，此均勻定態(tài)解為非雙曲奇點，其穩(wěn)定性的計算非常復雜，在這里暫不討論，僅討論兩個均勻定態(tài)解的情形.由文獻［12］可知道E1為鞍點，E2為結(jié)點或焦點，因此只需要考慮E2附近的微擾分析即可.為了計算上的方便，用a表示a2.

定理3 假定d1＝d2＝0，且系統(tǒng)（5）存在兩個正均勻定態(tài)解，則如果下面兩個條件之一滿足，E2就為漸近穩(wěn)定的.

（b）1-2c-2c2m＞0，且A2（1-c）（1+m+cm）＞1+2cm.

系統(tǒng)（5）微擾方程的系數(shù)為

因此給出圖靈分岔的必要條件為

由微擾方程的系數(shù)可知a11的值是恒小于零的，根據(jù)上面的條件可以得出這樣一個結(jié)論：如果圖靈分岔發(fā)生，則一定有a22＞0.這意味著易感者宿主對系統(tǒng)起著阻滯子的作用，染病者宿主對系統(tǒng)起著活化子的作用.

Hopf分岔的必要條件為

下面分析圖靈分岔的穩(wěn)定性.

2 圖靈斑圖的穩(wěn)定性

首先來分析當取定一組具體參數(shù)值時，系統(tǒng)發(fā)生分岔的區(qū)域.設m為分支參數(shù)，令d2＝1，A＝1.2，c＝0.2.借助Maple軟件計算可以得到當m＜4時正均勻定態(tài)解存在.圖靈分岔發(fā)生的首要條件是系統(tǒng)對均勻微擾必須是穩(wěn)定的，由前面可以知道Δ0＞0.只需要計算tr0＜0時，m滿足的條件即可.通過計算可以知道0 ＜m＜3.884 558 064.當m＝3.884 558 064 時，tr0＝0，系統(tǒng)存在Hopf分岔.圖靈分岔滿足的另一個條件是系統(tǒng)對于某些模數(shù)的微擾是不穩(wěn)定的，會出現(xiàn)鞍結(jié)點分岔.即系統(tǒng)滿足條件

通過分析計算條件，給出了參數(shù)m與擴散系數(shù)d1的示意圖（見圖1）.參數(shù)只有在區(qū)域Ⅰ才會出現(xiàn)圖靈分岔，在區(qū)域Ⅱ與Ⅲ中，圖靈分岔條件不滿足，但系統(tǒng)穩(wěn)定與否需要進一步確定.如果固定d1或m，穩(wěn)定性如何變化，即特征值正負如何呢？若取定d1＝15，經(jīng)計算可以得到m的臨界值有兩個，分別為m1c＝3.322，m2c＝0.572 4.當0≤m≤0.572 4或3.322＜m＜3.884 6時，系統(tǒng)會出現(xiàn)圖靈不穩(wěn)定.從圖2 可以看到，當m＜0.572 4時，隨著m的增加圖靈斑圖區(qū)域逐漸減?。ㄒ妶D2的實線部分）；當m＞3.322時，隨著m的增加，圖靈斑圖區(qū)域逐漸增大（見圖2的虛線部分）.同時也可以發(fā)現(xiàn)m＞3.323 的波是m＜0.572 4時的波向左平移.

事件，涉及的內(nèi)容廣泛，幾乎涵蓋了我們生活的方方面面。將事件作為一種營銷工具，有具有新聞價值的事件、大型活動、慶典、節(jié)事活動等。我們可以將事件營銷從規(guī)模和類型上進行分類。

圖1 參數(shù)m隨參數(shù)d1的變化圖Fig.1 Region graph of parameter mwith d1

圖2 參數(shù)m與特征值實部的變化圖Fig.2 Graph of parameter mand real part of eigenvalues

3 時空斑圖分析

當d1＝d2＝0時，即常微分（ODE）系統(tǒng)在一定條件下存在兩個平衡點，其中染病者宿主密度較小的為鞍點，較大的為結(jié)點或焦點.當存在兩個正平衡點時，ODE 系統(tǒng)性態(tài)可能有下面幾種：平衡點穩(wěn)定，但不存在極限環(huán)；平衡點穩(wěn)定，存在不穩(wěn)定的極限環(huán)；平衡點不穩(wěn)定，但不存在極限環(huán)；平衡點不穩(wěn)定，存在穩(wěn)定的極限環(huán).當d1，d2≠0時，即偏微分（PDE）系統(tǒng)性態(tài)比較復雜，只考慮穩(wěn)定與不穩(wěn)定兩種情形.下面取一組固定的值來分析ODE 系統(tǒng)與PDE 系統(tǒng)的可能組合.取d1＝15，d2＝1，A＝1.2，c＝0.2，具體組合見表1.在此的穩(wěn)定性指正均勻穩(wěn)態(tài)解的穩(wěn)定性.

作者發(fā)現(xiàn)，當取上面這組值時，有些組合還沒有，因此取另一組值d2，A，c不變，d1＝3，具體組合見表2.由表1 和 表2 發(fā)現(xiàn)，系統(tǒng)存在穩(wěn)定的極限環(huán)的幾種形式還不包含，在此不再列出.原因是當ODE系統(tǒng)存在穩(wěn)定的極限環(huán)時，相應的PDE系統(tǒng)是不穩(wěn)定的，通過數(shù)值模擬發(fā)現(xiàn)系統(tǒng)出現(xiàn)一般的圖靈斑圖，沒有什么新的現(xiàn)象發(fā)生.

下面分別對上面幾種情形進行模擬，看m處在不同區(qū)域時，系統(tǒng)的動力學性態(tài)如何，性態(tài)是否有規(guī)律可循.從表1 直觀分析系統(tǒng)的性態(tài)，可以得到情形1 和情形3 都屬于ODE 系統(tǒng)平衡點穩(wěn)定且不存在極限環(huán)，而相應的PDE 系統(tǒng)是不穩(wěn)定的，因此可以知道系統(tǒng)滿足圖靈不穩(wěn)定的條件.在這幾種情形下圖靈斑圖發(fā)生，圖靈斑圖究竟是什么形狀呢？下面利用第2部分的有限差分方程的方法來進行數(shù)值模擬.對于表1 中的情形1，分別取m＝0，0.1，0.3，0.5進行模擬.

表1 當d1＝15，d2＝1，A＝1.2，c＝0.2時，ODE與PDE系統(tǒng)動力學性態(tài)分類Tab.1 The different dynamic behaviors of ODE and PDE systems when d1＝15，d2＝1，A＝1.2，c＝0.2

表2 當d1＝3，d2＝1，A＝1.2，c＝0.2時，ODE與PDE系統(tǒng)動力學性態(tài)分類Tab.2 The different dynamic behaviors of ODE and PDE systems when d1＝3，d2＝1，A＝1.2，c＝0.2

從圖3（a）可以看到，當m＝0時，系統(tǒng)的圖靈斑圖完全是點狀的；當m＝0.1時，系統(tǒng)基本上是條狀斑圖.因此可以很自然地想到，m在0與0.1之間連續(xù)變化時系統(tǒng)的斑圖是逐漸由點狀過渡到條狀斑圖的；當m＝0.3時，發(fā)現(xiàn)系統(tǒng)的圖靈斑圖是點條共存的，但是這種斑圖與圖3（b）的點條共存正好是相反的（見圖3（c）；當m＝0.5時，系統(tǒng)的圖靈斑圖變?yōu)辄c狀的，這種點狀斑圖與圖3（a）的點狀也是正好相反的（見圖3（d）），被稱為缺口狀斑圖.因此有m在0.1與0.572 4之間連續(xù)變化時，系統(tǒng)由前一種的條狀圖靈斑圖逐步變?yōu)楹笠环N的條狀圖靈斑圖，然后由條狀逐漸變?yōu)楹笠环N點狀圖靈斑圖.換句話說，隨著染病者劑量m∈（0，0.572 4）由小到大變化，染病者宿主的密度分布由疏到密.

表1 的情形2，對于ODE 與PDE 系統(tǒng)的正均勻平衡解均穩(wěn)定，這種情形不是本文考慮的重點.對于表1的情形3和情形4，通過數(shù)值模擬發(fā)現(xiàn)，m在這兩個區(qū)域內(nèi)的圖靈斑圖類似（盡管情形4 中ODE系統(tǒng)存在極限環(huán)）.分別給出m＝3.6與m＝3.8 時，系統(tǒng)的斑圖（見圖4），可以看到圖4（b）是缺口狀斑圖到迷宮斑圖的過渡.

圖3 d1＝15，d2＝1，A＝1.2，c＝0，2時，系統(tǒng)（5）的圖靈斑圖Fig.3 Turing pattern of system（5）when d1＝15，d2＝1，A＝1.2，c＝0.2

圖4 d1＝15，d2＝1，A＝1.2，c＝0.2時，系統(tǒng)（5）的斑圖Fig.4 Spatial pattern of system（5）when d1＝15，d2＝1，A＝1.2，c＝0.2

對于表1 的最后一種情形，取m在這個區(qū)間的一個值3.95，其它參數(shù)不變.用同樣的方法進行數(shù)值模擬發(fā)現(xiàn)，反應擴散方程出現(xiàn)迷宮斑圖，如圖5 所示.由文獻［13］知道迷宮斑圖的形成需要兩個必要條件：①系統(tǒng)經(jīng)歷橫向失穩(wěn)；②當兩個波峰互相靠近時，它們之間會產(chǎn)生相互排斥作用.這種排斥作用能夠使波峰靠近的速度放慢以至最后停止.也就是說，兩個波峰不會因合并而湮滅.在此，我們想是不是由于極限環(huán)的存在，引起圖靈斑圖失穩(wěn)而引起波峰的變化，從而產(chǎn)生迷宮斑圖.

圖5 d1＝15，d2＝1，A＝1.2，c＝0.2，m＝3.95時，系統(tǒng)（5）的迷宮斑圖Fig.5 Labyrinth patterns of system（5）when d1＝15，d2＝1，A＝1.2，c＝0.2，m＝3.95

進一步，對表2中的幾種情形分別進行了討論.第1種情形不做討論，首先從第2種情形來討論.對于情形2，與情形1的區(qū)別是存在一個不穩(wěn)定的極限環(huán)，就由于這點區(qū)別，系統(tǒng)的動力學性態(tài)發(fā)生了很大的改變.取m＝3.8，利用有限差分的方法進行數(shù)值模擬發(fā)現(xiàn)：在這組參數(shù)值下，系統(tǒng)（5）有螺旋波發(fā)生（如圖6 所示）.

圖6 d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.8時，系統(tǒng)（5）的螺旋波斑圖Fig.6 Spiral waves pattern of system（5）when d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.8

對于表2的情形3，發(fā)現(xiàn)盡管不存在極限環(huán)，但在這個區(qū)間上會由共振形成一種靜態(tài)波，如圖7 所示.另外，給出了變量y隨時間t的變化圖（如圖7（b）），發(fā)現(xiàn)隨著時間的增加，變量y最終趨于一穩(wěn)態(tài).

圖7 d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.9時，系統(tǒng)（5）的斑圖與變量y隨時間t變化的關(guān)系圖Fig.7 Spatial patterns of system（5）and relationship diagram of variable y and time twhen d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.9

對于情形4，發(fā)現(xiàn)在這個區(qū)間上會出現(xiàn)迷宮斑圖，給出了m＝3.95時分別取空間步長為1，時間步長為0.05，運行1萬次和10萬次時斑圖的變化情況（見圖8）.由圖8 可以發(fā)現(xiàn)，隨著時間的增加，染病者宿主的分布逐漸變得稀疏.另外，從表2 的這幾種情形分析，將宿主的擴散率固定在某一數(shù)值時，隨著染病者劑量函數(shù)m的由小變大，系統(tǒng)（5）的斑圖由行波、靜態(tài)波到迷宮斑圖.

圖8 d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.95時，系統(tǒng)（5）的迷宮斑圖Fig.8 Labyrinth patterns of system（5）when d1＝3，d2＝1，A＝1.2，c＝0.2，m＝3.95

通過上面的數(shù)值模擬可以發(fā)現(xiàn)，螺旋波的出現(xiàn)是由于擴散率的變化引起的，因此考慮擴散率的變化對系統(tǒng)斑圖的變化是非常自然的事情.固定m＝3.8來分析擴散系數(shù)d1的變化，動力學性態(tài)會有什么改變.當m＝3.8，A＝1.2，c＝0.2，d2＝1時，常微分系統(tǒng)始終存在一不穩(wěn)定的極限環(huán).分別取d1＝0.1，0.5，1，3，3.5，4，利用有限差分法進行數(shù)值模擬（如圖9），其中（a）～（e）為典型的螺旋波斑圖，（f）為混沌斑圖.

由圖9 可以發(fā)現(xiàn)，當d1＝0.1，1，3.5時，系統(tǒng)中的缺陷數(shù)目有多個；當d1＝0.5，3時，系統(tǒng)中的缺陷數(shù)目只有一個；當d1＝4時，發(fā)現(xiàn)螺旋波失穩(wěn)導致時空混沌（見圖9（f））；繼續(xù)增大d1，化為表1的第4種情形.我們知道螺旋波是由缺陷為中心自組織形成的一類特殊的行波，對于一個穩(wěn)定的螺旋波斑圖態(tài)，系統(tǒng)中的缺陷（或缺陷密度）很少，并且他們的數(shù)目不隨時間變化.但是，如果系統(tǒng)中的控制變量超過某些臨界值時，螺旋波會自發(fā)地產(chǎn)生出新的缺陷，而每個缺陷都趨向于產(chǎn)生新的螺旋波.因此，系統(tǒng)中的缺陷數(shù)目會隨著時間以指數(shù)的形式增加，直到系統(tǒng)達到一個飽和的缺陷密度.此時系統(tǒng)中被缺陷充滿，它的長期有序現(xiàn)象不復存在，系統(tǒng)進入時空混沌態(tài).

圖9 不同d1，d2＝1，A＝1.2，c＝0.2，m＝3.8 時，系統(tǒng)（5）螺旋波的不同斑圖Fig.9 Different spiral waves patterns of system（5）when different d1，d2＝1，A＝1.2，c＝0.2，m＝3.8

4 討論

討論了當k＝1 和k＝2 時系統(tǒng)（5）的時空斑圖.當k＝1時，系統(tǒng)始終是穩(wěn)定的，不會發(fā)生圖靈斑圖；當k＝2時，根據(jù)系統(tǒng)穩(wěn)定性及極限環(huán)的存在性進行了分類.發(fā)現(xiàn)一些非常有趣的現(xiàn)象，當ODE系統(tǒng)穩(wěn)定，且不存在極限環(huán)相應的PDE系統(tǒng)不穩(wěn)定時，則PDE 系統(tǒng)出現(xiàn)圖靈斑圖，且圖靈斑圖可能有兩種形式，每種形式包括點狀與條狀斑圖.若ODE系統(tǒng)存在一不穩(wěn)定的極限環(huán)，而PDE系統(tǒng)是穩(wěn)定的話，出現(xiàn)螺旋波.若ODE 系統(tǒng)盡管不存在極限環(huán)，但平衡點是不穩(wěn)定的，而PDE系統(tǒng)是穩(wěn)定的時，系統(tǒng)同樣會出現(xiàn)行波，但這種行波是靜態(tài)的.若ODE 系統(tǒng)不存在極限環(huán)，但平衡點是不穩(wěn)定的，而PDE 系統(tǒng)也是不穩(wěn)定的時，則系統(tǒng)出現(xiàn)迷宮斑圖.通過數(shù)值模擬猜測這種規(guī)律是可循的.在本文中，由于時間的關(guān)系，僅考慮了一種特殊形式的系統(tǒng)的性態(tài)，即k＝2.若對于k取更大的值，PDE系統(tǒng)的性態(tài)可能更復雜，詳細分析可以查閱文獻［14］.

關(guān)于考慮空間因素的傳染病模型的研究，已有一些文獻，有考慮一般傳染病模型（見文獻［15-16］），有考慮具體疾病的模型（見文獻［17-18］）.若考慮空間因素還可以應用到各個領(lǐng)域，比如生態(tài)環(huán)境等（見文獻［19］）.另外，考慮反應擴散的傳染病模型還有很多文獻，在此不一一列出.

［1］Crofton H D.A quantitative approach to parasitism［J］.Parasitology，1971，62：179-194.

［2］Anderson R M.The regulation of host population growth by parasitic speies［J］.Parasitology，1978，76：119-157.

［3］Anderson R M，Whitefield P J，Dobson A P.Experimental studies of infection dynamics：infection of the definitive host by the cercariae of transversotrema patialense［J］.Parasitology，1978，77：189-200.

［4］Anderson R M.Population dynamics of snail infection by microcidia［J］.Parasitology，1978，77：201-224.

［5］Keymer A E，Anderson R M.The dynamics of infection of tribolium confusum by hymenolopis diminuta：the influence of infective-stage density and spatial distribution［J］.Parasitology，1979，79：195-207.

［6］Keymer A E.The dynamics of infection of tribolium confusum by hymenolopis diminuta，the influence of exposure time and host density［J］.Parasitology，1982，84：157.

［7］Keymer A E.Population dynamics of hymenolepis diminutain the intermediate host［J］.Journal of Animal Ecology，1981，50：941-950.

［8］Anderson R M.Experimental studies of age-prevalence curves of schistosoma mansoni infection in populations of biomphalaria glabrata［J］.Parasitology，1984，89：79-104.

［9］Anderson R M，Whitefield P J.An experimental study of the population dynamics of anectoparasitic digenean transversotrema patialense（soparker）：cercarial and adult stages［J］.Journal of Animal Ecology，1977，46：555-580.

［10］Anderson R M，May R M.The population dynamics of microparasites and their invertebrate hosts［J］.Phil.Trans.R.Soc.Lond.B，1981，291：451-524.

［11］Regoes R R，Ebert D.Dose-dependent infection rates of parasites produce the allee effectin epidemiology［J］.Proc.R.Soc.Lond.B，2002，269：271-279.

［12］Li Guihua.Bifurcation analysis of an epidemic model with nonlinear incidence［J］.Applied Mathematics and Computation，2009，214（2）：411-423.

［13］Britton N F.Essential Mathematical Biology［M］.London：Springer，2003.

［14］李桂花，王穩(wěn)地.傳染病動力學模型性態(tài)分析［D］.西南大學，2008.

［15］Wang Weiming，Cai Yongli.Complex dynamics of a reaction-diffusion epidemic model［J］.Nonlinear Analysis：Real World Applications，2012，13（5）：2240-2258.

［16］Sun Guiquan，Jin Zhen，Liu Quqnxing，et al.Pattern formation in a spatial S-I model with non-linear incidence rates［J］.J.Stat.Mech.，2007，11：11011.

［17］Wang K，Wang W，Song S.Dynamics of an HBV model with diffusion and delay［J］.J.Theoret.Biol.，2008（1）：36-44.

［18］Xu R，Ma Z.An HBV model with diffusion and time delay［J］.J.Theoret.Biol.，2009，257（3）：499-509.

［19］Sun Guiquan，Li Li，Wang Zike.Spatial dynamics of a vegetation model in an arid flat environment［J］.Nonlinear Dynamics，2013，73（4）：2207-2219.