孙得宁,张存华
(兰州交通大学 数理学院,甘肃 兰州 730070)
人体主要依靠两大免疫系统灭活病原体,即先天免疫系统和适应性免疫系统.当病原体进入肺部后,首先先天免疫系统中的巨噬细胞对入侵的病原体立即作出反应,当巨噬细胞检测到自己无法遏制的病原体时,它们会召集中性粒细胞(PMNs),它是人体内比巨噬细胞更具有杀伤力的免疫细胞[1].然而通过一系列的临床研究发现,当中性粒细胞与淋巴细胞的比率(NLR)超过一定数值时,人体内会发生炎症反应.文献[2]中Wang D等发现,新冠患者在重症期间中性粒细胞计数增加且淋巴细胞计数减少(即NLR增加),这一现象表明患者的病情可能处于危重症状态,所以对这一比率(NLR)的研究对于患者的早期治疗具有十分重要的意义.基于以上医学现象与理论,Pugliese等[3]提出了一个描述病原体与宿主免疫反应之间相互作用的系统,但比较有局限性;文献[4]又对该系统进行了推广.事实上,作为先天免疫系统中对病原体最具有杀伤力的免疫细胞,PMNs的数量在正常水平下较低,但它随着病原体入侵肺部后迅速增加.因此,Young等[5]构建了如下系统:
(1)
α (2) α/β>δ/γ. (3) Young等[5]分析了在条件(3)下系统(1)的平衡点的数目与类型,并得到系统(1)的边界平衡点是全局渐近稳定的.然而参数的实际值在不同病原体间会有所差异,所以较大病原体可能不会战胜先天免疫系统,所以文献[5]又考虑了相反的情况: α/β<δ/γ. (4) Young等[5]研究了系统(1)在条件(4)下平衡点的类型,并发现在该条件下系统(1)存在更复杂的动力学现象,但没有给出具体的理论分析.Shi等[6]进一步研究了在条件(4)下系统(1)的动力学行为,并得到随着参数值的变化,系统(1)展现了复杂的动力学与分支现象:包括Hopf分支与余维二Bogdanov-Takens分支,最后给出了相应的相图与分支图,对所得到的结论进行了理论解释.然而,对于某些免疫功能较弱的老人与小孩,PMNs的自然死亡率将会增加,因此在系统(1)的基础上构建如下系统: (5) 本节中,考虑了系统(5)边界平衡点的存在条件及其全局稳定性.首先做如下尺度变换,令 则系统(5)化为(仍用x,z,t表示x1,z1,t1) (6) 其中: p1=m/α,p2=(α2b)/(βδ),p3=(αγ)/(βδ), p4=(βη)/(α2),p5=μ/β. 由(2)知p1>1,p2,p3,p4,p5均为正参数,又因为条件(4)等价于p3<1,所以在下列条件下研究系统(6), 0 (7) 且系统(6)的正不变区域为 Ω={(x,z)|x≥0,z≥0}. 证明为了讨论系统(6)的边界平衡点,设 则 令 f(x)=(p22-p22p3)x3+ 则 f′(x)=3x2(p22-p22p3)+ Δ=b2-4ac, (8) 作Poincare′变换,令x=1/ω,z=μ/ω,则系统(8)化为 (9) (10) 当ω=0时,系统(10)有两个临界点A′(0,0),B′(p3/(1-p5),0).临界点B′(p3/(1-p5),0),有两个特征值λ1,2>0且为重根,得到B′是不稳定结点.对于临界点A′(0,0),它的一个特征值为零,另一个特征值非零.接下来令μ=p2μ1-p2p3ω1,ω=p2p3μ1,dτ=-p2p3t,则系统(6)变为 (11) 根据文献[7]中的中心流形理论,得到限制在中心流形上的方程为 (12) 由于p3<1,则[p2(p3-1)]/p3<0,根据文献[8]中的定理7.1得到A′(0,0)是稳定的鞍结点,且通过做变换dτ=-p2p3t得到A′(0,0)是包含一个稳定抛物扇形的鞍结点,且ω<0的一部分是稳定的,而ω>0的一部分进入第一象限后分为两部分双曲扇形. 作Poincare′变换,令x=μ/ω,z=1/ω,则系统(8)化为 (13) (14) 则得到系统(14)临界点C′(0,0),其特征值均为零且关于临界点C′(0,0)的雅可比矩阵为零矩阵.再做极坐标变换,令v=rcos(θ),ω=rsin(θ),则系统(14)化为 (15) 其中 R1(θ)=cos(θ)(p2cos2(θ)-sin(θ)cos(θ)+ 令dτ=rdt,则系统(15)化为 (16) (a)系统(6)在C′附近的相图 图2 系统(6)边界平衡点的全局相图 本文主要研究了一类描述病原体与宿主免疫反应之间相互作用的先天免疫疫应答模型,讨论了该模型平衡点的存在条件,并通过Poincare′变换与中心流形理论证明了该系统的边界平衡点是全局渐近稳定的,最后利用MATLAB软件包对得到的理论结果进行了适当的数值模拟.1 边界平衡点的存在条件及其全局稳定性
(2p2-2p2p3+p1p2p3+
p22p4-p22p5)x2+
(1-p3+p1p3+2p2p4-2p2p5+
2p1p2p5)x+p4-p5+2p1p5-p12p5,
2x(2p2-2p2p3+p1p2p3+
p22p4-p22p5)+(1-p3+p1p3+
2p2p4-2p2p5+2p1p2p5),
sin2(θ)cos(θ)+p2p5sin2(θ)+p5sin3(θ)),
R2(θ)=-p1sin2(θ)cos2(θ)+(sin(θ)+
p2cos(θ))(p3cos(θ)+p5sin(θ)+sin(θ)),
U(θ)=cos(θ)sin(θ)(sin(θ)+p2cos(θ)).2 数值模拟
3 结语