一类Wolbachia氏菌在蚊群传播数学模型的全局动力学*

武丹，郭志明

广州大学数学与信息科学学院，广东 广州 510006

为了有效减少登革热的传播，目前应对的唯一途径是防蚊灭蚊，比如，通过大量喷洒杀虫剂减少蚊子的数量，可以切断蚊媒传染病的传播途径。但是这种方法会带来极大的生态问题。另外，也可能导致蚊子会对杀虫剂逐渐产生抗药性，使得灭蚊效果降低。考虑到这些情况，可见防蚊灭蚊的措施亟待改进。应该用更长远的、更环保的方法防止蚊媒传染病的传播。目前，有一类新兴的控制蚊媒传染病的方法，就是使自然界中的伊蚊种群感染内共生菌Wolbachia（沃尔巴克氏）。一系列的实验证明这样可以阻断蚊子传播登革热等疾病［1-3］。Wolbachia可以通过CI机制成功入侵自然界中的伊蚊种群，使得在伊蚊种群中，由于有CI机制的作用，野外雌蚊与感染了Wolbachia的雄蚊交配时，产下的卵不能正常发育，但是对于感染了Wolbachia的雌蚊没有影响，并且其后代大多数都会感染Wolbachia，因此Wolbachia在入侵自然界的蚊群时拥有传播的优势。基于Wolbachia的蚊媒控制方法是安全的。风险分析表明，释放感染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