一類Wolbachia氏菌在蚊群傳播數(shù)學(xué)模型的全局動(dòng)力學(xué)*

武丹，郭志明

廣州大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，廣東 廣州 510006

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

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

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

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

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

關(guān)于Wolbachia入侵動(dòng)力學(xué)及其他傳染病模型可參見文獻(xiàn)［7-16］以及其中的參考文獻(xiàn)。例如Keeling等［8］建立的常微分方程模型，Huang 等［9-10］建立了反應(yīng)擴(kuò)散方程及其衍生的常微分方程模型，Hu 等［11］建立的隨機(jī)微分方程模型，Li等［12］建立的離散競爭模型，Huang等［13］建立的時(shí)滯微分方程。

1 平衡點(diǎn)的存在性和穩(wěn)定性分析

下面分析各平衡點(diǎn)的穩(wěn)定性。

定理2對(duì)于模型（3），以下結(jié)論成立。

縱觀西方國家財(cái)務(wù)管理的相關(guān)經(jīng)驗(yàn)而言，西方國家尤為重視專業(yè)隊(duì)伍和組織的構(gòu)建，發(fā)揮隊(duì)伍的專業(yè)性帶頭作用，以提升財(cái)務(wù)管理的綜合水平。因此，對(duì)于我國而言，在大數(shù)據(jù)新型時(shí)代背景下，有必要借鑒國外西方國家的經(jīng)驗(yàn)，構(gòu)建專業(yè)隊(duì)伍和組織，邀請(qǐng)行業(yè)的專業(yè)帶頭人和企業(yè)內(nèi)的經(jīng)驗(yàn)豐富人員帶隊(duì)組建專業(yè)隊(duì)伍，為其他財(cái)務(wù)人員營造良好的學(xué)習(xí)氛圍和環(huán)境，同時(shí)能夠提供專業(yè)知識(shí)和技能的指導(dǎo)，幫助更多的財(cái)務(wù)人員增加知識(shí)面，豐富大數(shù)據(jù)下財(cái)務(wù)管理的相關(guān)經(jīng)驗(yàn)，其他人員能夠通過提高自身能力和技能以更快地融入該團(tuán)隊(duì)和隊(duì)伍之中。

（i）當(dāng)p1≤δ1，p2≤δ2時(shí)，E0是全局漸近穩(wěn)定的；

（ii）當(dāng)p1≤δ1，p2＞δ2時(shí)，E0不穩(wěn)定，E2是全局漸近穩(wěn)定的；

（iii）當(dāng)p1＞δ1，p2≤δ2時(shí)，E0不穩(wěn)定，E1是全局漸近穩(wěn)定的。

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

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

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

（i）當(dāng)yˉ≤xˉ時(shí)，E2不穩(wěn)定，E1是全局漸近穩(wěn)定的；

（ii）當(dāng)yˉ＞xˉ時(shí)，E1和E2都是局部漸近穩(wěn)定的，E*是鞍點(diǎn)。

證明

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

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

（ii）根據(jù)前面的計(jì)算，可以得到此時(shí)平衡點(diǎn)E1和E2都是局部漸近穩(wěn)定的。下面討論E*的穩(wěn)定性。E*處的Jacobi矩陣為

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

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

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

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

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

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

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

2 數(shù)值模擬與小結(jié)

為了更好地展示模型（3）的動(dòng)力學(xué)穩(wěn)定性，下面用MATLAB 來進(jìn)行數(shù)值模擬。在情形A（見圖4），選取的參數(shù)是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，無論初值數(shù)量取多少，E1是全局漸近穩(wěn)定的，感染蚊群會(huì)持續(xù)存在，未感染蚊群會(huì)滅絕，Wolbachia能夠成功入侵到整個(gè)蚊群當(dāng)中。

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

在情形B（見圖5），選取參數(shù)值為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都是局部漸近穩(wěn)定的，兩個(gè)平衡點(diǎn)有各自的吸引域，不同的初值數(shù)量會(huì)使得解趨于其中一個(gè)平衡點(diǎn)。當(dāng)解趨于E1時(shí)，感染蚊群會(huì)持續(xù)存在，未感染蚊群會(huì)滅絕，Wolbachia能夠成功入侵到整個(gè)蚊群當(dāng)中。而當(dāng)解趨于E2時(shí)，未感染蚊群持續(xù)存在，感染蚊群滅絕，Wolbachia不能成功入侵。

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

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