李雅芝，任新志

1.黔南民族师范学院 数学与统计学院，贵州 都匀 558000；2.黔南州复杂系统与智能优化重点实验室，贵州 都匀 558000；3.西南大学 数学与统计学院，重庆 400715

蚊子是世界上最致命的动物之一，由其传播的疾病每年导致数百万人死亡[1]．因此，持续的蚊虫控制工作对于预防这些疾病的爆发非常重要．

Wolbachia是一种革兰氏阴性菌，是世界上分布最为广泛的共生菌，并且可以人工植入埃及伊蚊体内[2]．研究发现Wolbachia可以抑制登革热病毒在蚊子体内的复制并且可以降低病毒的传染性[3]，故使用Wolbachia控制蚊媒传染病的传播备受关注．Wolbachia对蚊子的繁殖影响主要有：垂直传播和细胞质分离(CI)，即携带Wolbachia的雌性的下一代也会携带此种菌，正常雌性与携带Wolbachia的雄性的下一代会因细胞质分离而死亡．

近年来，许多文献利用数学模型研究了Wolbachia在蚊虫种群中的传播．文献[4]提出并分析了Wolbachia感染蚊子种群与未感染蚊子种群之间的基本竞争模型，然后，利用反馈控制技术设计了Wolbachia的引入方案．文献[5-7]分别建立了具有脉冲一般出生和死亡率函数的脉冲模型、具有脉冲出生和投放的性别结构和状态依赖脉冲的综合控制模型，对使用Wolbachia控制蚊媒传染病的各种控制策略进行了研究．文献[8]考虑了环境的异质性，建立了两种机制随机切换的数学模型．研究发现：在均匀环境中维持的Wolbachia的初始状态在异质环境中会灭绝，频繁的环境转换有利于Wolbachia的传播．此外，文献[9-11]等都讨论了宿主、媒介和疾病之间的相互作用．

本文将建立具有成熟时滞的蚊子种群模型，通过理论分析和数值模拟讨论时滞对Wolbachia传播的影响．

1 模型建立

文献[12]建立了如下数学模型描述Wolbachia在蚊子种群中的传播：

(1)

其中：I(t)表示被Wolbachia感染的蚊子数量，U(t)表示未被Wolbachia感染的蚊子数量；β∈[0，1]是垂直传播的可能性；q∈[0，1]是CI影响的可能性；b>0是出生率；d>0是死亡率；D≥0是适合度损失．模型(1)考虑了Wolbachia的垂直传播和CI两种影响，由文献[3]可知，对于某种特定的Wolbachia菌株，其垂直传播率和CI影响发生的可能性都近似接近1，也即，如果雌性蚊子感染了Wolbachia，其后代也会感染；正常雌性蚊子与感染雄性蚊子的后代会因为发生CI而死亡．所以，我们假设β=1，q=1，并忽略掉适合度的影响(D=0)．蚊子的一生有4个阶段(卵、幼虫、蛹、成虫)，从卵到成虫大约需要14～21 d，雌性蚊子的寿命约3个月，故蚊子的成熟时间不能被忽略，此因素可用成熟时滞描述．而成熟时滞又受季节性影响，故考虑为周期时间依赖的时滞，记为τ(t)．根据文献[13]的方法，将τ(t)引入模型(1)中得到如下模型：

(2)

2 解的适定性与平凡解

f1(t，φ)=(1-τ′(t))e-δτ(t)bφ1(-τ(t))-d(φ1(0)+φ2(0))φ1(0)

因为τ(t)是ω-周期的，故有f(t+ω，φ)=f(t，φ)．所以，模型(2)是一个ω-周期泛函微分系统．

假设g(t)是一个连续的ω-周期函数，令

则对于系统(2)解的适定性有如下结果．

证在χ+的每个紧子集上，f(t，φ)关于φ是连续且Lipschitz的．所以，对任意φ∈χ+，系统(2)在其最大存在区间上有唯一的具有初值S0=φ的解S(t，φ)．

定义bc=(1-τ′(t))e-δτ(t)b．由系统(2)可知：

I′(t)=bcI(t-τ(t))-d(I(t)+U(t))I(t)

我们令

I′(t)=0

则有

易知

又令

V′(t)=0

并由

可得到

再令

对任意给定的ρ≥1，令

则[0，ρx*]χ是关于χ的有序区间．

容易证明：对任意ψ∈[0，ρx*]χ(ψi(0)=0～(ψi(0)=ρxi*)，i=1，2)和t∈R，有fi(t，ψ)≥0～(fi(t，ψ)≤0)．进一步地，[0，ρx*]χ对于系统(2)是正不变的．通过选择任意大的ρ即可得到解关于χ+的正性和有界性．证毕．

定理1系统(2)有一个种群灭绝平衡态(0，0)，且是不稳定的．

I′(t1)=(1-τ′(t1))e-δτ(t1)bI(t1-τ(t1))-d(I(t1)+U(t1))I(t1)

由于|Y(t)|<ε，则有I(t1)<ε和U(t1)<ε．所以，

我们取

则有

3 成熟时滞的影响

下分析E1和E2稳定性，将系统(2)在(I*，U*)处线性化后可写为

x′(t)=Ax(t)+Bx(t-r)

(3)

其中x(t)=(I(t)，U(t))T且

这里

a11=-d(2I*+U*)a12=-dI*a21=-dU*a22=-d(I*+2U*)

且有

假设系统(2)具有形如x(t)=ceλt(c≠0)的解，则可得到系统(3)的特征方程如下：

0=det(λI-A-e-λrB)

展开得

λ2-(tr(A))λ+det(A)+e-2λrdet(B)+e-λr[det(a1|b2)+det(b1|a2)-(tr(B))λ]=0

其中：a1，a2分别为矩阵A的第一列和第二列；b1，b2分别为矩阵B的第一列和第二列．

可以计算得到

tr(A)=-3d(I*+U*)det(A)=2d2(I*+U*)2

(λ+be-δr-be-δre-λr)(λ+2be-δr-be-δre-λr)=0

(4)

由特征根的符号易知E1是不稳定的．

(λ+e-δrb)(λ+2be-δr-be-δre-λr)=0

(5)

易知E2是渐近稳定的．

综上可得如下结果：

定理2在区域{(I，U)：I，U≥0}中，当τ(t)=r时，系统(2)有两个非平凡平衡态E1和E2，且E1是不稳定的，E2是局部渐近稳定的．

由定理2可知：如果垂直传播是完全的，Wolbachia最终将完全入侵蚊子种群．图1a将系统(2)在有时滞((μ(t)，ν(t))=(20，100))和无时滞((I0，U0)=(20，100))两种情况下的动力学行为做了比较，比较结果见表1．表1说明：无时滞的ODE系统将高估蚊子总量，这对蚊子种群的控制是不利的，因此成熟时滞对于蚊子种群的动力学行为有重要影响．

表1 图1a结果对比

图1b将不同时滞对系统(2)动力学行为的影响做了对比．可以看出：时滞越小，系统(2)解趋于稳定的时间越短，也即，如果环境有利于蚊子的成熟，Wolbachia完全入侵蚊子种群的速度会加快．但从蚊子总量看，相比于大时滞，小时滞会使环境中有更多的蚊子．在实际中，我们应该寻找平衡这两种互反影响的方法．

参数取值为b=20，d=0.07，δ=0.1．图1 时滞对系统(2)动力学行为的影响

当成熟时滞是周期函数τ(t)时，我们将通过数值模拟讨论时滞对系统(2)动力学行为的影响．根据广州市的登革热数据，拟合出成熟时滞表达式为：

其中

并且

α=23.6β=13γ=196

参数取值为b=20，d=0.07，δ=0.1．图2 不同历史值对系统(2)周期解的影响

图2显示了不同历史值对系统(2)周期解的影响．可以看出：(μ(t)，ν(t))=(20，100)时，正常蚊子种群将会灭绝；(μ(t)，ν(t))=(10，200)时，正常蚊子种群与携带Wolbachia的蚊子种群共存，并呈现周期波动．从取值来看，当携带Wolbachia的蚊子种群在历史值中所占比例较大时，Wolbachia会完全入侵野生蚊子种群．

图3在图2历史值的基础上，研究了不同出生率对系统(2)周期解的影响．图3a中(μ(t)，ν(t))=(10，200)，当b=20时两种蚊子种群共存，当b=40时正常蚊子种群灭绝，图3b中(μ(t)，ν(t))=(20，100)，当b=10时两种蚊子种群共存，当b=20时正常蚊子种群灭绝，说明出生率越大，越有利于Wolbachia的传播．

参数值取为d=0.07，δ=0.1．图3 不同出生率对系统(2)周期解的影响

4 结语

本文主要研究成熟时滞对Wolbachia传播的影响．研究发现：① 对于常数时滞，时滞越大越有利于疾病的控制；② 携带Wolbachia的蚊子种群在历史值中所占比例越大，越有利于Wolbachia的传播；③ 蚊子种群出生率越大，越有利于Wolbachia的传播．