武丹,郭志明
廣州大學(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)證明這樣可以阻斷蚊子傳播登革熱等疾病[1-3]。Wolbachia可以通過CI機(jī)制成功入侵自然界中的伊蚊種群,使得在伊蚊種群中,由于有CI機(jī)制的作用,野外雌蚊與感染了Wolbachia的雄蚊交配時(shí),產(chǎn)下的卵不能正常發(fā)育,但是對(duì)于感染了Wolbachia的雌蚊沒有影響,并且其后代大多數(shù)都會(huì)感染W(wǎng)olbachia,因此Wolbachia在入侵自然界的蚊群時(shí)擁有傳播的優(yōu)勢(shì)?;赪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)部的競(jìng)爭(zhēng)引起的,而蚊群的內(nèi)部競(jìng)爭(zhēng)主要發(fā)生在未成年階段,所以單位時(shí)間內(nèi)成年個(gè)體不再是線性增長。另外,基于成年個(gè)體對(duì)資源的競(jìng)爭(zhēng)假設(shè),模型(2)的死亡函數(shù)為二次函數(shù)。但是大量的事實(shí)表明,成蚊之間的競(jìng)爭(zhēng)幾乎可以忽略?;谝陨峡紤],本文建立一個(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]建立的離散競(jìng)爭(zhēng)模型,Huang等[13]建立的時(shí)滯微分方程。
下面分析各平衡點(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)勢(shì);而當(dāng)yˉ>xˉ時(shí),環(huán)境對(duì)野外蚊子更加有利,此時(shí)Wolbachia具有適合度劣勢(shì)。下面具體分析這兩種情況的平衡點(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)
為了更好地展示模型(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