一類Wolbachia氏菌在蚊群傳播數學模型的全局動力學*

武丹，郭志明

廣州大學數學與信息科學學院，廣東 廣州 510006

為了有效減少登革熱的傳播，目前應對的唯一途徑是防蚊滅蚊，比如，通過大量噴灑殺蟲劑減少蚊子的數量，可以切斷蚊媒傳染病的傳播途徑。但是這種方法會帶來極大的生態問題。另外，也可能導致蚊子會對殺蟲劑逐漸產生抗藥性，使得滅蚊效果降低?？紤]到這些情況，可見防蚊滅蚊的措施亟待改進。應該用更長遠的、更環保的方法防止蚊媒傳染病的傳播。目前，有一類新興的控制蚊媒傳染病的方法，就是使自然界中的伊蚊種群感染內共生菌Wolbachia（沃爾巴克氏）。一系列的實驗證明這樣可以阻斷蚊子傳播登革熱等疾病［1-3］。Wolbachia可以通過CI機制成功入侵自然界中的伊蚊種群，使得在伊蚊種群中，由于有CI機制的作用，野外雌蚊與感染了Wolbachia的雄蚊交配時，產下的卵不能正常發育，但是對于感染了Wolbachia的雌蚊沒有影響，并且其后代大多數都會感染Wolbachia，因此Wolbachia在入侵自然界的蚊群時擁有傳播的優勢?；赪olbachia的蚊媒控制方法是安全的。風險分析表明，釋放感染Wolbachia的蚊子對人類、動物以及環境造成的風險可以忽略不計。因此，利用Wolbachia來阻斷登革熱的傳播具有重要的應用價值。早在2005 年，生物學家奚志勇就在實驗室成功地使得Wolbachia在埃及伊蚊中穩定傳播［3］。隨著Wolbachia蚊媒控制技術的發展，建立數學模型研究Wolbachia傳播動力學具有了重要的意義。通過數學模型分析研究Wolbachia傳播動力學，可以為釋放感染Wolbachia的蚊子提供很好的策略和方案。

1990 年以來，Hoffmann 和Turelli 建立了世代不重疊的離散模型并做了解釋與分析［4-5］。2014 年，Zheng等［6］為研究Wolbachia傳播動力學建立了時滯微分方程

顯然，當t→∞時，RF(t) →0且RM(t) →0. 令RF=RM= 0，忽略（1）中的時滯，得到常微分方程組

注意到，模型（2）中出生函數為線性函數。一般說來，隨著種群規模的增加，其成年個體的出生率和未成年個體的幸存率將下降，這是由于種群內部的競爭引起的，而蚊群的內部競爭主要發生在未成年階段，所以單位時間內成年個體不再是線性增長。另外，基于成年個體對資源的競爭假設，模型（2）的死亡函數為二次函數。但是大量的事實表明，成蚊之間的競爭幾乎可以忽略?；谝陨峡紤]，本文建立一個新的常微分方程模型，出生和幸存函數為Ricker型函數，而死亡函數用線性函數。

具體來說，我們建立以下的模型的規模。

關于Wolbachia入侵動力學及其他傳染病模型可參見文獻［7-16］以及其中的參考文獻。例如Keeling等［8］建立的常微分方程模型，Huang 等［9-10］建立了反應擴散方程及其衍生的常微分方程模型，Hu 等［11］建立的隨機微分方程模型，Li等［12］建立的離散競爭模型，Huang等［13］建立的時滯微分方程。

1 平衡點的存在性和穩定性分析

下面分析各平衡點的穩定性。

定理2對于模型（3），以下結論成立。

縱觀西方國家財務管理的相關經驗而言，西方國家尤為重視專業隊伍和組織的構建，發揮隊伍的專業性帶頭作用，以提升財務管理的綜合水平。因此，對于我國而言，在大數據新型時代背景下，有必要借鑒國外西方國家的經驗，構建專業隊伍和組織，邀請行業的專業帶頭人和企業內的經驗豐富人員帶隊組建專業隊伍，為其他財務人員營造良好的學習氛圍和環境，同時能夠提供專業知識和技能的指導，幫助更多的財務人員增加知識面，豐富大數據下財務管理的相關經驗，其他人員能夠通過提高自身能力和技能以更快地融入該團隊和隊伍之中。

（i）當p1≤δ1，p2≤δ2時，E0是全局漸近穩定的；

（ii）當p1≤δ1，p2＞δ2時，E0不穩定，E2是全局漸近穩定的；

（iii）當p1＞δ1，p2≤δ2時，E0不穩定，E1是全局漸近穩定的。

由圖1可知，在情形A，軌線都趨于E1；而在情形B，取不同的初值，軌線分別趨于E1和E2. 對于這兩種情況，合理的生物學解釋是當yˉ≤xˉ時，環境對于感染了Wolbachia的蚊子更有利，此時Wolbachia具有適合度優勢；而當yˉ＞xˉ時，環境對野外蚊子更加有利，此時Wolbachia具有適合度劣勢。下面具體分析這兩種情況的平衡點穩定性。

圖1 L1和L2的兩種情況Fig.1 The two possible cases of L1 and L2

定理3在p1＞δ1，p2＞δ2的情況下，有

（i）當yˉ≤xˉ時，E2不穩定，E1是全局漸近穩定的；

（ii）當yˉ＞xˉ時，E1和E2都是局部漸近穩定的，E*是鞍點。

證明

（i）根據前面的計算和分析可以得到，此時平衡點E0，E2是不穩定的，E1是局部漸近穩定的。如果系統內部有周期軌線，軌線內部必有平衡點，這是不可能的。由于平衡點E0，E2不穩定，其附近的軌線只能遠離奇點，而平衡點E1是穩定的，其附近的軌線只能靠近奇點，因此也不可能出現連接奇點的奇異閉軌線。由于系統的正半軌LP+有界，根據文獻［17］中的定理有，軌線的ω極限集Ωp只能是奇點構成的連通閉集。由于奇點是孤立的，因此任何第一象限的軌線最終都趨于奇點。

由于從E0出發的軌線必遠離E0，E0不在Ωp中；又因E2的穩定流形是y軸，所以從第一象限（除y軸外）出發的軌線都遠離E2，這說明E2?Ωp. 由于E1局部漸近穩定，因此Ωp={E1}，E1全局漸近穩定。

（ii）根據前面的計算，可以得到此時平衡點E1和E2都是局部漸近穩定的。下面討論E*的穩定性。E*處的Jacobi矩陣為

由于-q1δ1δ2x*＜0，則平衡點E*是鞍點，其穩定流形將E1和E2的吸引域分開，系統內部沒有周期軌線。由于平衡點E0，E*不穩定，其附近的軌線遠離奇點，而平衡點E1，E2穩定，其附近的軌線靠近奇點，不可能出現連接奇點的奇異閉軌線。因此任何第一象限的軌線最終趨于奇點，不存在極限環。后面數值模擬的結果也說明極限環不存在。

定理3說明如果環境對于受到Wolbachia感染的蚊群更加有利，那么給定任意初值，Wolbachia都能夠成功入侵。相反，如果環境不利于受到Wolbachia感染的蚊群，由于鞍點E*的穩定流形將E1和E2的吸引域分開，則對于不同的初值，Wolbachia可能成功入侵，也可能入侵失敗。

圖2 垂直等傾線L1和水平等傾線L2對第一象限的分解Fig.2 The decomposition of the first quadrant by vertical isocline L1 and horizontal isocline L2

通過分析第一象限各點(f，g)的方向，可知D2和D4都是正不變的。假如有一個初值(x0，y0) ∈D2，那么對于任意的時間t，有f＞0，g＜0，當t→∞時，有(x(t)，y(t)) →E1，因此，D2在E1的吸引域中。同理，D4在E2的吸引域中，因此，E*的穩定流形位于D1和D3中。根據D1和D3的定義，可知在D1和D3中的任何解(x(t)，y(t))，都有

因此穩定曲線是一條光滑且嚴格增加的函數，定義其為y=h(x). 在D1中，穩定曲線的α極限集一定是E0，因此曲線是連接E0和E*的異宿軌。在D3中，f＜0，g＜0，假設當t→-∞時，x和y的極限存在，表示為

圖3 所選取的參數值為p1= 0.5，p2= 0.9，δ1= 0.17，δ2= 0.2，q1= 0.1，q2= 0.1. 圖3 表明：E*的穩定流形的位置可以用光滑且嚴格增加的函數圖像來表示，這個圖像從E0出發，正向嚴格遞增到無窮。

圖3 h(x)的大致位置Fig.3 The approximate position of h(x)

2 數值模擬與小結

為了更好地展示模型（3）的動力學穩定性，下面用MATLAB 來進行數值模擬。在情形A（見圖4），選取的參數是p1= 9 × 10-4，p2= 4 × 10-4，δ1= 2 × 10-4，δ2= 3 × 10-4，q1= 1 × 10-9，q2= 1 × 10-9. 分別取初值x0= 5 × 106，y0= 1 × 104（圖4（a））；x0= 1 × 103，y0= 1 × 108（圖4（b）），得到圖中的解曲線。圖4表明：在情形A，無論初值數量取多少，E1是全局漸近穩定的，感染蚊群會持續存在，未感染蚊群會滅絕，Wolbachia能夠成功入侵到整個蚊群當中。

圖4 情形A的解曲線Fig.4 The solution curves of case A

在情形B（見圖5），選取參數值為p1= 4 × 10-4，p2= 9 × 10-4，δ1= 2 × 10-4，δ2= 3 × 10-4，q1= 1 ×10-9，q2= 1 × 10-9. 分別取初值x0= 5 × 106，y0= 1 × 104（圖5（a））；x0= 1 × 103，y0= 1 × 108（圖5（b）），得到圖中的解曲線。圖5表明：在情形B，E1和E2都是局部漸近穩定的，兩個平衡點有各自的吸引域，不同的初值數量會使得解趨于其中一個平衡點。當解趨于E1時，感染蚊群會持續存在，未感染蚊群會滅絕，Wolbachia能夠成功入侵到整個蚊群當中。而當解趨于E2時，未感染蚊群持續存在，感染蚊群滅絕，Wolbachia不能成功入侵。

圖5 情形B的解曲線Fig.5 The solution curves of case B

圖6 情況I時的解曲線Fig.6 The solution curves of case I