1、高校应用数学学报2023,38(1):27-36具有进化效应的SIR模型的稳定性和Hopf分支分析杨洪1,2,马秋敏1,盛江林1,仲兆满1,赵琛森1(1.江苏海洋大学 理学院,江苏连云港 222005;2.江苏海洋大学 江苏省海洋资源开发研究院,江苏连云港 222005)摘要:主要研究一类具有进化效应的SIR模型的动力学行为.首先,建立数学模型,并证明系统解的存在性,正性,有界性等基本性质.其次,分析系统正平衡点的稳定性,并给出健康个体的防御成本和已感染个体的感染能力引起健康个体和已感染个体数量振荡的充分条件.然后,利用中心流形定理,讨论Hopf 分支的性质.最后,利用合适的参数值,对系统进行
2、数值模拟.关键词:SIR模型;进化效应;稳定性;Hopf分支中图分类号:O175.13文献标识码:A文章编号:1000-4424(2023)01-0027-101引言一直以来,许多学者都利用SIR模型研究传染病的传播规律1-7.本文研究以健康个体,已感染个体和已恢复个体为研究对象的SIR模型,是基于文献8中研究的模型,进一步,考虑健康个体自身的防御成本和已感染个体感染能力的适应性变化对系统动力学行为的进化效应影响9-11.由于健康个体往往通过自身免疫系统产生防御作用来驱动健康个体适应环境,同时健康个体防御成本的降低可能导致病毒的变异,进而增加攻击能力.于是这类问题的研究更贴近现实生活,更有实际
3、意义.与现有文献中的SIR模型相比,考虑到健康个体的内禀增长率取决于自身对病毒的防御成本,同时,感染率由健康人群防御成本和已感染者感染能力共同决定.于是考虑的模型为dSdt=r(u)(1 S+IK)S(u,v)SIS+I,dIdt=(u,v)SIS+I(+)I,dRdt=I R.(1)收稿日期:2022-01-14修回日期:2022-08-15基 金 项 目:中 国 博 士 后 基 金(2020M681521);江 苏 省 创 新 创 业 训 练 计 划 一 般 项 目(202011641103Y;202111641102Y);江苏省博士后项目(2021K456C);连云港市博士后项目(LYG
4、20210012);国家自然科学基金(72174079)DOI:10.13299/ki.amjcu.00224828高 校 应 用 数 学 学 报第38卷第1期模型参数的生物学意义见表1,其中S表示健康个体的数量,I表示已感染个体的数量,R表示已恢复个体的数量.函数r(u),(u,v)的假设如下.(1)防御成本u和健康个体的内禀增长率r之间的关系假设为r(u)=r0ec1u2(参见9).因此随着人群总数接近环境承载能力,r对动力学的影响减小.当健康人群占有率低时,相同水平的防御成本较高,而当健康人群占有率高时,相同水平的防御成本较低.(2)易感人群的感染率为(u,v)=01+e(uv),这是由
5、防御成本和感染能力决定的10.感染率随着防御成本的增加而减小,随着感染能力的增强而增加.表1模型(1)符号的生物意义参数描述u健康个体的防御成本v已感染个体的感染能力r(u)内禀增长率K环境最大容纳量已感染个体的死亡率(u,v)感染率个体的恢复率已恢复个体的死亡率r0最大增长率c1防御成本0最大感染率感染率在不同人群之间差异的敏感度这里系统(1)的初始条件为S(0)=S0 0,I(0)=I0 0,R(0)=R0 0.并且为满足实际的生物意义,假设S(0)+I(0)0是存在且唯一的,并且保持非负有界.证证证在初始条件下,系统(2)关于(S,I)满足局部Lipschitz条件,说明存在某个c 0,
6、使得其解在区间0,c)内存在且唯一.为了证明正性,首先证明对于任意t 0,S和I的有界性.由于S(0)=0,意味着S(t)0.所以假设S(0)0,对于系统(2)的第一个方程来说,S0(t)=S(t)r(u)(1 S(t)+I(t)K)(u,v)I(t)S(t)+I(t),对于任意t 0得到S(t)=S(0)et0r(u)(1S(x)+I(x)K)(u,v)I(x)S(x)+I(x)dx.由此可得S的正性.同理可得I的正性.杨洪等:具有进化效应的SIR模型的稳定性和Hopf分支分析29于是(S+I)0=r(u)(1 S+IK)(+)I r0K(S+I K)S,在S(0)+I(0)0,由系统(2)
7、的解的正性可得S(t)r(u)(1 S+IK)S,从而limtS(t)K.令V(t)=S(t)+I(t),计算V(t)关于时间t的导数得V(t)=S(t)+I(t)=r(u)(1 S(t)+I(t)K)S(t)(+)I(t).对于一个非常小的正数m,不妨设m (+),有V(t)+mV(t)r(u)(1 S(t)K)S(t)+mS(t).因此存在一个正常数,使得V(t)+mV(t)L.因此可以得到V(t)(V(0)Lm)emt+Lm,则limtS(t)+I(t)Lm.因此,系统(2)中解是最终一致有界的.3 稳定性与分支分析下面对系统(2)在平衡点处进行稳定性和分支分析,并给出分支周期解存在的充
8、分条件.首先,研究系统平衡点的稳定性.显然系统(2)存在唯一的无病平衡点E1=(K,0).定定定理理理3.1当(u,v)(+)时,无病平衡点E1是局部渐近稳定的.反之E1是不稳定的.证系统(2)在无病平衡点E1=(K,0)处的Jacobian矩阵是JE1=(r(u)r(u)(u,v)0(u,v)(+),在E1处的特征方程为(r(u)(u,v)(+)=0.显然只需(u,v)(+),就有1 0,2(+)时,E1是不稳定的.定定定理理理3.2当(u,v)(+)时,无病平衡点E1是全局渐近稳定的.证构造Lyapunov函数V(t)=I(t),于是V=I=(u,v)SIS+I(+)I(u,v)S (+)
9、SSI=(u,v)(+)S.当(u,v)(+)时,V 0.当V=0 时,当且仅当S=K,I=0,故依据LaSalle不变集原理12可得,无病平衡点E1是全局渐近稳定的.下面分析正平衡点E=(S,I)的存在性和稳定性.当系统参数满足下面条件(+)(u,v)0,b11 0,b22 0.计算E处Jacobian矩阵对应的特征方程为2 Tr(JE)+Det(JE)=0,其中Tr(JE)=a11+b22,Det(JE)=a11b22 a22b11.于是特征方程的特征根为1,2=(u,v)iw(u,v)=Tr(JE)Tr(JE)2 4Det(JE)2.那么特征方程的根有负实部当且仅当a11+b22 0,(
10、u,v)0(H1).下面主要考虑(u,v)的变化所引起对系统动力学性质的影响.在分析过程中以(u,v)作为分支参数.并且,当(+)(u,v)0.定定定理理理3.3假设(H1)成立,对于系统(1)有如下结果.(1)如果a11 0且(u,v)0且(u,v)0且b22 0,当(u,v)=(+)r(u)(1 2S+IK)(S+IS I)2(+)时,则系统的平衡点E处产生Hopf 分支.接下来给出分支周期解的稳定性.为此通过平移 x=S S,y=I I,将E平移到原点.因此系统(2)变为 x=r(u)(1(x+S)+(y+I)K)(x+S)(u,v)(x+S)(y+I)(x+S)+(y+I),y=(u,
11、v)(x+S)(y+I)(x+S)+(y+I)(+)(y+I).可写成(x y)=J(x y)+(f(x,y,(u,v)h(x,y,(u,v).杨洪等:具有进化效应的SIR模型的稳定性和Hopf分支分析31其中f(x,y,(u,v)=x2(I2(+)3(u,v)2S3r(u)K)x y(r(u)K+2I(+)3(u,v)S)2)+y2(+)3(u,v)2S)x3(I2(+)4(u,v)3S4)+x2 y(2SI I2)(+)4(u,v)3S4)+x y2(2I S)(+)4(u,v)S)3)y3(+)4(u,v)3S2+O(|(x,y)|4),h(x,y,(u,v)=x2(I2(+)3(u,v
12、)2S3)+x y(2I(+)3(u,v)S)2)y2(+)3(u,v)2S)+x3(I2(+)4(u,v)3S4)+x2 y(I2 2SI)(+)4(u,v)3S4)+x y2(S 2I)(+)4(u,v)S)3)+y3(+)4(u,v)3S2+O(|(x,y)|4).记(u,v)=,则1,2()=iw,w=a11b22 a22b11时,对应于特征值1,2()的一个特征向量是=(b11,a11 iw)T,有B=(b110a11w)B1=(1b110a11wb111w).通过转换(?x?y)=B(x y),得到(?x?y)=(0ww0)(?x?y)+(f1(?x,?y)h1(?x,?y),即f
13、1(?x,?y)=1b11f(b11?x,a11?x w?y,(u,v),h1(?x,?y)=1wh(b11?x,a11?x w?y,(u,v)a11wb11f(b11?x,a11?x w?y,(u,v),其中f(b11?x,a11?x w?y,(u,v)=b211?x2(I2(+)3(u,v)2S3r(u)K)+b11?x(a11?x w?y)(r(u)K+2I(+)3(u,v)S)2)+(a11?x w?y)2(+)3(u,v)2S)+(b11?x)3(I2(+)4(u,v)3S4)+(b11?x)2(a11?x w?y)(2SI I2)(+)4(u,v)3S4)b11?x(a11?x w
14、?y)2(2I S)(+)4(u,v)S)3)(a11?x w?y)3(+)4(u,v)3S2+O(|(b11?x,a11?x w?y)|4),32高 校 应 用 数 学 学 报第38卷第1期h(b11?x,a11?x w?y,(u,v)=b211?x2(I2(+)3(u,v)2S3)b11?x(a11?x w?y)(2I(+)3(u,v)S)2)(a11?x w?y)2(+)3(u,v)2S)(b11?x)3(I2(+)4(u,v)3S4)+(b11?x)2(a11?x w?y)(I2 2SI)(+)4(u,v)3S4)b11?x(a11?x w?y)2(S 2I)(+)4(u,v)S)3)
15、+(a11?x w?y)3(+)4(u,v)3S2+O(|(b11?x,a11?x w?y)|4).将关于?x和?y的方程改写为以下极坐标形式 r=(u,v)r+(u,v)r3+O(|r|4),=w(u,v)+c(u,v)r2+O(|r|4),在(u,v)=处进行泰勒展开得 r=0()(u,v)r+()r3+O(|r|4),=w()+w0()(u,v)+c()r2+O(|r|4).周期解的稳定性是由()决定的,所以有()=116(f1?x?x?x+f1?x?y?y+h1?x?x?y+h1y?y?y)+116w()f1?x?y(f1?x?x+f1?y?y)h1?y?x(h1?x?x+h1?y?y
16、)f1?x?xh1?x?x+f1?y?yh1?y?y,其中f1?x?x?x=6b211I2(+)4(u,v)3S46a11b11(2SI I2)(+)4(u,v)3S4+6a211(2I S)(+)4(u,v)S)3+6a311(+)4b11(u,v)3S2,f1?x?y?y=2w2(2I S)(+)4(u,v)S)36a11w2(+)4b11(u,v)3S2,h1?x?x?y=2b211(I2 2SI)(+)4(u,v)3S44a11b11(S 2I)(+)4(u,v)S)3+6a211(+)4(u,v)3S2+2a11b11(2SI I2)(+)4(u,v)3S44a211(2I S)(+)4(u,v)S)3+6a311(+)4b11(u,v)3S2,h1?y?y?y=6w2(+)4(u,v)3S2+6a11w2(+)4b11(u,v)3S2,f1?x?y=wr(u)K+2I(+)3(u,v)S)2+2wa11(+)3b11(u,v)2S+2b11w?x(2SI I2)(+)4(u,v)3S4 2w(2a11?x w?y)(2I S)(+)4(u,v)S)3)+6a11w(a11?
