一类具有病毒饱和增殖的病毒-免疫动力学模型的分析_马香香.pdf

高校应用数学学报2023,38(1):53-63一类具有病毒饱和增殖的病毒-免疫动力学模型的分析马香香,李建全(陕西科技大学 数学与数据科学学院,陕西西安 710021)摘要:在假设病毒增殖率为Michaelis-Menten函数的基础上,提出了一类病毒增殖具有饱和性的病毒与特异性免疫细胞相互作用的模型.分析发现该模型至多有两个正平衡点并会发生鞍结点分支;借助中心流形定理讨论了平衡点的局部稳定性;运用Bendixson-Dulac定理排除了周期解的存在性,进而得到模型的全局动力学性态.数值模拟显示了病毒与免疫系统相互作用的结果对初始状态的依赖性,以及在作用过程中会出现病毒载量和免疫细胞种群数量的持续振荡.关键词:病毒免疫系统;平衡点;稳定性;鞍结点分支中图分类号:O175.26文献标识码:A文章编号:1000-4424(2023)01-0053-111引言临床数据表明,长期存在于宿主体内的病毒(如HIV,HBV和HCV)在感染宿主,并形成慢性感染的过程中,病原体会抑制免疫反应1-7.这包括已经证实了的HBV 和HCV会损伤特异性辅助细胞应答8-9和HIV会破坏特定的免疫细胞4.但目前病毒感染和抗病毒免疫应答的复杂过程仍未被完全了解.为了从理论上分析病毒慢性感染形成的过程,很多的数学模型被提出,其中一些并未考虑此过程中免疫系统的作用10-11.随着现代医疗水平的提高和科学技术的发展,更多的病毒和免疫反应之间的相互作用机制被发现.相应地,一些更多考虑病毒与免疫反应相互作用机理的数学模型相继被提出.Komarova等12在既考虑免疫系统抑制病毒增长,又考虑病毒对免疫系统的抑制和刺激两方面作用的情形下提出了一类病毒免疫反应动力学模型y0=ry(1 yK)ay pyz,z0=cyz1+dy qyz bz.(1)该模型只考虑了病毒(y)和免疫细胞(z)两个种群的作用,其中假设病毒增殖率函数为ry(1 yK),病毒以pyz的速率被免疫细胞消灭.还假设病毒对免疫的刺激和抑制都依赖于病毒载量y,其中收稿日期:2022-01-10修回日期:2022-08-24基金项目:国家自然科学基金(11971281;12071268)DOI:10.13299/ki.amjcu.00225054高 校 应 用 数 学 学 报第38卷第1期免疫被病毒的抑制率用qyz表示,而免疫细胞的增殖率用具有半饱和形式的Michaelis-Menten函数cyz1+dy来表示,这里c为免疫细胞的最大增殖率系数,1/d为增殖率的半饱和常数.在模型(1)中,ay和bz分别表示在没有免疫应答作用时病毒的失活率和免疫细胞的自然死亡率.对于模型(1),Wang等13得到了决定其动力学性态的两个阈值,进一步解释了这两个阈值所对应的生物学意义,并讨论了当抗病毒治疗不能完全清除病毒时,如何引入免疫治疗来有效控制病毒的感染与增殖.同时,基于模型(1),一些反映不同机理的病毒免疫动力学模型也被不断提出.2014年,Shu等14将免疫细胞增殖的滞后效应纳入病毒与免疫应答的模型.2018年,Li等15提出了特异性免疫细胞具有钟形增殖率的慢性病毒感染模型.一些临床数据16-18显示,在感染HIV的最初几周内免疫CD4+T细胞被大量消耗,导致未被感染的靶细胞数量减少,进而使得病毒复制速度减慢.此外,先天免疫应答也会减缓病毒复制.考虑这些因素,Boer19建立了一类病原体增殖受到自身密度制约的模型,其中用具有饱和形式的Michaelis-Menten函数rL1+L来表示病原体的增殖率.这里,L代表病原体,r为最大增殖率,1/为半饱和常数.基于Boer19提出的病原体动力学模型,本文假设在没有免疫细胞作用时病毒的变化率服从模型y0=ry+y ay,(2)再根据模型(1)所反映的病毒-免疫反应机理,便可得到如下的病毒与免疫细胞相互作用的模型y0=ry+y ay pyz,z0=cyz1+dy qyz bz,(3)其中y和z含义与模型(1)相同,所有参数均为正数,且其生物学意义与模型(1)和(2)中的相同.本文将在理论分析模型(3)的动力学性态的基础上,进一步了解病毒与免疫系统作用中所能呈现的生物学现象,为控制病毒感染提供一定的理论依据.2 解的有界性显然y=0和z=0是系统(3)的两条解轨线,于是根据自治系统轨线的唯一性知,模型(3)具有正初始条件(即y(0)=y0 0,z(0)=z0 0)的解也是正的.因此模型(3)在R2+=(y,z):y 0,z 0上是不变的.对于(y,z)R2+,由系统(3)的第一个方程有y0=ry+y ay pyz y(r+y a).(4)对于u0=ru+uau,当ra 时,limt+u(t)=0;当ra 时,limt+u(t)=ra,K.因此,由(4)知,对于模型(3),当ra 时,limt+y(t)=0;当ra 时,limsupt+y(t)K.又当y=0时,模型(3)的第二个方程变为z0=bz,即有limt+z(t)=0.因此,当ra 时,对于模型(3)有limt+y(t)=limt+z(t)=0.这意味着当ra 时病毒和免疫细胞最终都消亡.于是以下仅考虑ra,即K 0的情形.马香香等:一类具有病毒饱和增殖的病毒-免疫动力学模型的分析55对于(y,z)R2+,由系统(3)有(y+pcz)0=ry+y ay(1 11+dy+qc)pyz pbczry+y ay pbcz r (y+pcz),其中=mina,b.因此limsupt+y(t)+pcz(t)r.所以,当ra 时,模型(3)是最终有界的,并且集合D=(y,z)0 y K,0 y+pcz r为模型(3)的一个不变集.于是,下面将研究ra 的情形下模型(3)在D上的动力学性态.3平衡点的存在性令系统(3)中的y0=0,z0=0,得y(r+y a pz)=0,z(cy1+dy qy b)=0.(5)易知方程组(5)有解(y,z)=(0,0)和(K,0),其中K=ra.它们对应系统(3)的边界平衡点O(0,0)和E0(K,0).由于点O的坐标均为零,所以它表示病毒和特异性免疫细胞均灭绝,而平衡点E0的y坐标为正的,z 坐标为零,这反映了免疫细胞灭绝而病毒存在的情形,故称E0为病毒控制平衡点.若方程组(5)有正解存在,则其对应模型(3)的正平衡点.其反映病毒存在但被免疫系统所控制的情形,故称正平衡点为免疫控制平衡点.由方程组(5)知,免疫控制平衡点由方程组r+y a pz=0,cy1+dy qy b=0(6)来决定.方程组(6)的第二个方程仅含有y,并且可改写为cy (1+dy)(qy+b)=0,即f(y),qdy2(c q bd)y+b=0.注意到方程f(y)=0若有实根,则必同号.记它的判别式为=(cqbd)24qbd,则f(y)=0存在正根的条件为c q+bd+2bdq,c0(q).直接求解可知,当c c0(q),即 0时,f(y)=0有两个正根y1=c q bd(c q bd)2 4qbd2qd,y2=c q bd+(c q bd)2 4qbd2qd.当c=c0(q),即=0时,f(y)=0有一个正根y3=bqd.由(6)的第一方程得z=1p(r+y a).(7)56高 校 应 用 数 学 学 报第38卷第1期于是将y=yi(i=1,2,3)代入上式,得zi=1p(r+yi a)=ap(+yi)(ra)yi=ap(+yi)(K yi).因此,当yi 0.此时方程组(6)的解(yi,zi)就对应模型(3)的免疫控制平衡点.下面讨论模型(3)免疫控制平衡点的存在性.当c c0(q)时,y1 K y2|(c q bd)2qdK|(1+dK)q+bd+bK,g1(q).注意到c0(q)g1(q)=2bdqdqKbK 0,且c0(q)g1(q)=0当且仅当q=bdK2.因此,当c g1(q)时,模型(3)有唯一的免疫控制平衡点E1(y1,z1).当c c0(q)时,y2(c q bd)2 4qbd.(8)这时,c (1+2dK)q+bd,g0(q)是必须的.进一步,(8)等价于c g1(q).由于c0(q)bdK2,且此时有g0(q)g1(q),所以当q bdK2且c0(q)c g1(q)时,模型(3)有两个不同的免疫控制平衡点E1(y1,z1)和E2(y2,z2).当c=c0(q)时,如果y3bdK2时,模型(3)有唯一的免疫控制平衡点E3(y3,z3);如果y3 K,即q bdK2且c=g1(q)时,模型(3)除病毒控制平衡点外,还存在一个免疫控制平衡点?E1(?y1,?z1),其中?y1=bqdK,?z1=1p(a(K+)+?y1 a).注意,此情形下的?y1可用y1的表达式来表示,故仍记?E1(?y1,?z1)为E1(y1,z1).综上分析,关于模型(3)的平衡点存在性有如下结论.定定定理理理3.1模型(3)总有病毒和免疫细胞均灭绝平衡点O(0,0)和病毒控制平衡点E0(K,0).当下列条件之一成立时,模型(3)有唯一的免疫控制平衡点E1(y1,z1).(i)q bdK2且c g1(q);(ii)q bdK2且c g1(q).当q bdK2且c0(q)c bdK2且c=c0(q)时,模型(3)有唯一的免疫控制平衡点E3(y3,z3).根据定理3.1的条件,模型(3)平衡点的存在条件可对应于平面(q,c)中的不同区域,其中D0=(q,c):q bdK2,c g1(q)(q,c):q bdK2,c g1(q)(q,c):q bdK2,c g1(q),D2=(q,c):q bdK2,c0(q)c bdK2,c=c0(q).于是依定理3.1,当(q,c)D0时,模型(3)无免疫控制平衡点;当(q,c)D1时,模型(3)有唯一免疫控制平衡点E1(y1,z1);当(q,c)D2时,系统(3)有两个不同的免疫控制平衡点E1(y1,z1),E2(y2,z2);当(q,c)D3时,系统(3)有一个免疫控制平衡点E3(y3,z3).图1在平面(q,c)上表示了区域Di(i=0,1,2,3).对于函数f(y)=qdy2(c q bd)y+b,根据定理3.1,当模型(3)存在唯一的免疫控制平衡点E1时,有f(K)0;当模型(3)有两个不同的免疫控制平衡点E1和E2时,有f(K)0;当马香香等:一类具有病毒饱和增殖的病毒-免疫动力学模型的分析57qcD1D0D2D3bdK2图 1系统(3)免疫控制平衡点存在性,粗线代表函数c0(q),细线代表函数g1(q)模型(3)有唯一的免疫控制平衡点E3时,有f(K)0;此外,当模型(3)没有免疫控制平衡点时,有f(K)0.此结论将用于判断病毒控制平衡点E0的稳定性.又因为f0(y)=2qdy (c q bd)=2qd(y c q bd2qd),以及y1cqbd2qd和y3=cqbd2qd,所以有f0(y1)0,f0(y3)=0.(9)(9)在讨论免疫控制平衡点的局部稳定性时将用到.4平衡点的稳定性系统(3)在O(0,0)处的Jacobi矩阵为J(O)=(K00b).于是J(O)的两个特征值分别为1=K 0,2=b 0,因此平衡点O(0,0)总是不稳定的.因此平衡点O在D上的局部稳定性如下.定定定理理理4.1对于任意情况,平衡点O在D上是不稳定的.为了方便讨论平衡点E0的局部稳定性,将模型(3)中的r用a(+K)表示,则系统(3)在E0处的Jacobi矩阵为J(E0)=(aK+KpK0cK1+dK qK b).因此系统(3)在E0处的特征值为1=aK+K 0和2=cK1+dK qK b=f(K)1+dK,58高 校 应 用 数 学 学 报第38卷第1期所以当f(K)g1(q)时,平衡点E0是不稳定的;当f(K)0,即c g1(q)时,平衡点E0是局部渐近稳定的;当f(K)=0,即c=g1(q)时,E0是高阶平衡点.下面来分析f(K)=0时系统(3)在E0附近的动力学性态.对系统(3)做平移变换u=y K,v=z,将平衡点E0平移至原点得u0=a(K+)(u+K)+(u+K)au pKv puv aK