一类具有饱和发生率的时滞SIR模型
A Delayed SIR Model with Saturated Incidence Rate
DOI: 10.12677/pm.2026.162029, PDF, HTML, XML,   
作者: 房俊瀛, 马永峰:大连交通大学基础部,辽宁 大连
关键词: SIR模型时滞稳定性Hopf分支SIR Model Time Delay Stability Hopf Bifurcation
摘要: 本文研究了一类具有饱和发生率的时滞SIR模型,该模型通过引入饱和发生率,符合在实际生活中感染人数增长到一定程度后传播速度变缓的规律,使模型更加贴合实际情况。首先,以时滞为参数,利用二代再生矩阵的方法计算求得模型的基本再生数R0,并通过线性化、特征理论等方法对无病平衡点与地方病平衡点处的特征方程进行分析,得出时滞对无病平衡点与地方病平衡点的稳定性的影响:当R0 < 1时,无病平衡点在一定时滞范围内是渐近稳定的;当R0 > 1时,地方病平衡点随时滞的变化而发生变化,并给出了地方病平衡点处Hopf分支的存在条件,最后利用MATLAB进行数值模拟,验证理论结果的正确性,得到时滞会使系统的稳定性发生变化。
Abstract: This paper studies a class of delayed SIR models with a saturated incidence rate. By introducing the saturated incidence rate, the model conforms to the real-world law that the transmission speed slows down when the number of infections grows to a certain level, thereby making the model more consistent with practical scenarios. Firstly, taking time delay as a parameter, the basic reproduction number R0 of the model is calculated using the next-generation matrix method. Subsequently, the characteristic equations at the disease-free equilibrium point and the endemic equilibrium point are analyzed via linearization, characteristic theory, and other methods, leading to the derivation of the influence of time delay on the stability of these two equilibrium points: when R0 < 1, the disease-free equilibrium point is asymptotically stable within a certain range of time delay; when R0 > 1, the endemic equilibrium point varies with changes in time delay, and the existence condition for Hopf bifurcation at the endemic equilibrium point is provided. Finally, numerical simulations are performed using MATLAB to verify the correctness of the theoretical results, confirming that time delay can alter the stability of the system.
文章引用:房俊瀛, 马永峰. 一类具有饱和发生率的时滞SIR模型[J]. 理论数学, 2026, 16(2): 12-20. https://doi.org/10.12677/pm.2026.162029

1. 引言

随着全球人口增长和城市化进程的加快,传染病的发生和传播越来越受到人们的关注[1]。在传染病模型的研究中,SIR模型(易感–感染–恢复模型)是研究传染病传播的经典模型之一。然而,现实中的传染病传播往往存在时滞现象,如狂犬病、艾滋病等,在某个人被感染后,相应的症状并没有立刻表现出来,而是通过一段时间后爆发相应的症状。

1760年,Bernoulli通过分析分析欧洲各国天花感染数据构建了一个关于天花传染的数学模型[2],预测了天花的蔓延趋势,1927年,Kermack和McKendrick通过分析1906年孟买瘟疫数据,将人群划分为易感人群S (Susceptible)、感染人群I (Infected)和移出人群R (Removed),提出了著名的SIR模型[3]

传统的SIR模型没有将时滞这一因素考虑进来。在文中,徐瑞等引入了饱和发生率 βSI 1+αI ,探究了带有饱和发生率的传染病模型[4],考虑到时滞在SIR模型中的重要作用([5][6]等),本文综合了时滞因素、传染病传播会进入饱和状态与恢复人群会出现再次感染的情况,构建一个带有饱和发生率的时滞SIR传染病模型,以此来讨论无时滞情况与存在时滞时系统平衡点的稳定性以及在某些条件成立的情况下,时滞会在地方病平衡点 E * 产生Hopf分支。

2. 模型建立

本研究基于前人已有的研究成果,将总人口 N( t )=S( t )+I( t )+R( t ) 分为易感者 S( t ) 、感染者 I( t ) 和恢复者 R( t ) ,考虑到感染过程存在时间延迟 τ (如潜伏期),即当前感染率依赖于 tτ 时刻的易感者与感染者数量,所以采用饱和函数 βS( tτ )I( tτ ) 1+αI( tτ ) 代表传染率,并且在恢复项中加入时滞可以更真实的反映感染过程中的时间延迟,即感染者在 tτ 时刻被感染后,需要经过免疫激活和病毒清除等一系列过程后,才能在 t 时刻开始恢复,加入时滞后能够使模型更加准确地描述疾病的传播动力学。因此在其基础上构建了一个带有饱和发生率的时滞SIR传染病模型。

{ dS( t ) dt =bμS( t ) βS( tτ )I( tτ ) 1+αI( tτ ) +γR( t ), dI( t ) dt = βS( tτ )I( tτ ) 1+αI( tτ ) μI( t )dI( tτ ), dR(t) dt =dI( tτ )μR( t )γR( t ). (2.1)

其中 S( t ) I( t ) R( t ) 分别表示在 t 时刻易感染人群、已感染人群和恢复人群的数量, b 代表人口迁入率, μ 代表自然死亡率, γ 代表恢复人群失去免疫力后重新成为易感染人群的比率, d 代表已感人群自然恢复率, β 代表疾病的传播系数, α 反映感染抑制强度。

3. 平衡点的稳定性分析

首先定义系统(2.1)的基本再生数 R 0 = βb μ( μ+d ) ,其中 R 0 的来源取决于人口的新增率、死亡率、治愈率和传染率。用 E=( S E , I E , R E ) 表示系统的平衡点,经过计算可得系统总有一个无病平衡点 E 0 =( b μ ,0,0 ) ,当 R 0 >1 时有一个地方病平衡点 E * =( S * , I * , R * ) ,其中

S * = ( μ+d )( 1+α I * ) β , I * = μ( μ+d )bβ μα( μ+d )β( μ+d )+ βγd μ+γ , R * = d μ+γ I *

简化系统(2.1),将系统(2.1)中的三个方程相加,并且种群大小为 N( t )=S( t )+I( t )+R( t ) ,我们有

dN( t ) dt =bμN( t ) (3.1)

显然, N( t )= b μ 是(3.1)的解,

则(3.1)的通解为

N( t )= b μ +C e μt  ( C )

因此

lim t N( t )= b μ

因此,系统(2.1)的 ω 极限集在 S+I+R= b μ 平面上,简化系统得

dI( t ) dt = βI( tτ ) 1+αI( tτ ) [ b μ I( tτ )R( tτ ) ]μI( t )dI( tτ )P( I,R ), dR( t ) dt =dI( tτ )μR( t )γR( t )Q( I,R ).

系统(2.1)的平衡点和极限环必在平面 S+I+R= b μ 平面上,因此下面考虑该简化系统的平衡点及其分支。当时滞为0时,考虑 I>0,R>0 ,取一个Dulac函数[7]

D( I,R )= 1+αI( t ) βI( t )

得到 ( DP ) I + ( DQ ) R <0

因此,我们由Bendixson-Dulac判据可以得到系统(2.1)不存在周期轨道[8]

接下来通过无量纲化处理,令 x= μ b I( t ) y= μ b R( t ) z=μt 对系统(1.1)进行简化得到

{ dx dz = βbx( zε ) μ 2 [ 1+ b μ αx( zε ) ] [ 1x( zε )y( zε ) ]xcx( zε ), dy dz =cx( zε )my. (3.2)

其中 ε=μτ c= d μ m=1+ γ μ

接下来将系统从(2.1)转化为(3.2)研究系统的平衡点及稳定性,计算得出无病平衡点 E 0 =( 0,0 ) ,地方病平衡点 E * =( x * , y * )

其中

x * = m[ βb μ 2 ( 1+c ) ] βb( c+m )+μmαb( 1+c ) ,  y * = βbc μ 2 c( 1+c ) βb( c+m )+μmαb( 1+c )

后面用二代矩阵法[9]计算系统的基本再生数 R 0

F ˜ =( βbx μ 2 +αbμx cx ),  V ˜ =( βb x 2 +βbxy μ 2 +αbμx +x+cx my )

得到

F= ( βb( μ 2 +αbμx )αβμ b 2 x ( μ 2 +αbμx ) 2 0 c 0 ) E 0 =( βb μ 2 0 c 0 )

V= ( βb( μ 2 +αbμx )( 2x+y )αβ b 2 μx( x+y ) ( μ 2 +αbμx ) 2 +1+c 0 0 m ) E 0 =( 1+c 0 0 m )

因此 F V 1 =( βb μ 2 ( 1+c ) 0 c 1+c 0 ) 得到基本再生数 R 0 = βb μ 2 ( 1+c )

X= ( x,y ) T ,对系统(3.2)线性近似后得到的线性方程为

X ( t )=AX( t )+BX( tτ )

其中

A( E 0 )=( 1 0 0 m ) B( E 0 )=( βb μ 2 c 0 c 0 )

A( E * )=( 1 0 0 m )     B( E * )=( βb μ 2 2 μ 2 bβ x * αμβ ( b x * ) 2 μ 2 βb y * ( μ 2 +αbμ x * ) 2 c βb x * μ 2 +αbμ x * c 0 )

3.1. 无病平衡点E0的稳定性

定理1 当 R 0 <1 时,无病平衡点 E 0 是局部渐近稳定的;当 R 0 >1 时,无病平衡点 E 0 是不稳定的。

证明:系统(3.2)在无病平衡点 E 0 处的特征方程为

f( λ )=| λEA e λε B |=( λ+m )[ λ+1 e λε ( βb μ 2 c ) ]

1) 当时滞 ε=0 时,我们有 λ 1 =m λ 2 =( R 0 1 )( 1+c )

① 当 R 0 <1 时, λ 1 <0 λ 2 <0 ,则无病平衡点 E 0 是局部渐近稳定的。

② 当 R 0 >1 时, λ 1 <0 λ 2 >0 ,则无病平衡点 E 0 是不稳定的。

2) 当时滞 ε>0 时,令 h( λ )=λ+1 e λε ( βb μ 2 c )

R 0 <1 时, h( λ ) 具有非实部的根,所以有

Re( h( λ ) )=Re( λ )+1 e Re(λ)ε ( βb μ 2 c )cos[ Im( λ )ε ]=0

Re( λ )=1+ e Re(λ)ε ( βb μ 2 c )cos[ Im( λ )ε ] β μ ( 1+c )=( R 0 1 )( 1+c )

由于 R 0 <1 ,则 Re( λ )<0 ,则可以得到 h( λ )=0 的所有根具有严格负实部,所以当 R 0 <1 时,无病平衡点 E 0 是局部渐近稳定的。

R 0 >1 时, h( 0 )=( 1 R 0 )( 1+c )<0

下面讨论 h( λ ) 的增减性,可以得到 h ( λ )=1+ε( βb μ 2 c ) e λε

① 当 β μ c0 时, h( λ ) ( ,+ ) 上递增,那么 Re( λ ) 单调递增且 h( 0 )<0 ,则存在一个正实数 λ 0 >0 ,使 h( λ 0 )=0 ,所以 R 0 >1 时无病平衡点 E 0 不稳定。

② 当 β μ c<0 时,由 h( 0 )=1 β μ +c<0 可以得到 β μ c>1 与前面的条件 β μ c<0 相矛盾。证毕。

3.2. 地方病平衡点E*的稳定性

计算得到系统(3.2)在地方病平衡点 E * 处的特征方程

A( E * )( a 11 0 0 a 22 ) B( E * )( b 11 b 12 b 21 0 )

f( λ )= λ 2 + p 1 λ+ p 2 + e λε ( q 1 λ+ q 2 )+s e 2λε

其中 p 1 = a 11 a 22 p 2 = a 11 a 22 q 1 = b 11 q 2 = a 22 b 11 s= b 12 b 21

1) 当时滞 ε=0 时,可以得到特征方程为

f( λ )= λ 2 +( p 1 + q 1 )λ+( p 2 + q 2 +s )=0

由Routh-Hurwitz准则知,特征方程 f( λ ) 全部解为负解的充要条件是

p 1 + q 1 >0 p 2 + q 2 +s>0

在满足上述不等式的情况下,地方病平衡点 E * 呈现局部渐近稳定状态。

2) 当时滞 ε>0 时,对于系统(3.2)的特征方程

λ 2 + p 1 λ+ p 2 + e λε ( q 1 λ+ q 2 )+s e 2λε =0

将其变形为

e λε ( λ 2 + p 1 λ+ p 2 )+( q 1 λ+ q 2 )+s e λε =0 (3.3)

当时滞 ε>0 时,假设 iθ 是方程 f( λ ) 的解,代入其中并分离实虚部可得,

{ ( θ 2 + p 2 +s )cos( θε )+( p 1 θ )sin( θε )= q 2 , p 1 θcos( θε )+( θ 2 + p 2 s )sin( θε )=θ q 1 .

解得

{ cos( θε )= q 2 ( θ 2 p 2 +s ) p 1 q 1 θ 2 θ 4 + p 1 2 θ 2 2 θ 2 p 2 + p 2 2 s 2 , sin( θε )= p 1 q 2 θ+θ q 1 ( θ 2 p 2 s ) θ 4 + p 1 2 θ 2 2 θ 2 p 2 + p 2 2 s 2 . (3.4)

根据(3.4)式的第一个式子我们可以得到

ε k = 1 θ arccos[ q 2 ( θ 2 p 2 +s ) p 1 q 1 θ 2 θ 4 + p 1 2 θ 2 2 θ 2 p 2 + p 2 2 s 2 ] 2kπ θ ,k=0,1,2,

cos 2 θε+ sin 2 θε=1 可以得到

θ 8 + h 1 θ 6 + h 2 θ 4 + h 3 θ 2 + h 4 =0. (3.5)

其中

h 1 =2 p 1 2 q 1 2 4 p 2 , h 2 = p 1 4 p 1 2 q 1 2 4 p 1 2 p 2 +2 p 2 q 1 2 +2 q 1 2 s+6 p 2 2 q 2 2 2 s 2 , h 3 =2 p 1 2 p 2 2 p 1 2 q 2 2 2 p 1 2 s 2 +4 p 1 q 1 q 2 s p 2 2 q 1 2 2 p 2 q 1 2 s q 1 2 s 2 4 p 2 3 +2 p 2 q 2 2 +4 p 2 s 2 2 q 2 2 s, h 4 = p 2 4 p 2 2 q 2 2 2 p 2 2 s 2 +2 p 2 q 2 2 s q 2 2 s 2 + s 4 .

θ 2 =w ,可以得到

w 4 + h 1 w 3 + h 2 w 2 + h 3 w+ h 4 =0

假设上式至少存在一个正实根 w k ,那么(3.5)也存在正实根 θ k = w k

接下来验证横截条件 { dRe( λ ) dε } λ=iθ 0 。对式子(3.3)两边同时关于 ε 求导并且化简得

( dλ dε ) 1 = e λε ( 2λ+ p 1 )+ q 1 sλ e λε λ e λε ( λ 2 + p 1 λ+ p 2 ) ε λ

λ=iθ 代入上式,分离实虚部可得

sign{ [ d( Reλ ) dε ] λ=iθ }=sign{ Re ( dλ dε ) λ=iθ 1 }=sign{ AC+BD C 2 + D 2 }

其中

A= p 1 cos( θε )2θsin( θε )+ q 1 , B=2θcos( θε )+ p 1 sin( θε ), C= θ 2 p 1 cos( θε )+( sθ+ p 2 θ θ 3 )sin( θε ), D= θ 2 p 1 sin( θε )+( sθ p 2 θ θ 3 )cos( θε ).

显然,当 AC+BD0 时,可以得到

sign{ Re ( dλ dε ) λ=iθ 1 }0

因此,当时滞 ε>0 ,并且满足 AC+BD0 p 1 + q 1 >0 p 2 + q 2 +s>0 时,我们有

① 若 ε[ 0, ε 0 ) ,系统(3.2)的地方病平衡点是局部渐近稳定的;

② 若 ε> ε 0 ,系统(3.2)的地方病平衡点是不稳定的;

③ 若 ε= ε 0 ,系统(3.2)的地方病平衡点正经历Hopf分支。

4. 数值模拟

在前文中,我们讨论了时滞对平衡点的影响,接下来通过数值仿真展示时滞对系统解的影响。我们选取参数 β=6,μ=0.3,c=2,m=3,α=0.4,b=0.1 。计算得 ε 0 1.1, R 0 2.22>1

ε=0.75< ε 0 时,

Figure 1. Temporal dynamics of the normalized infected and recovered individuals

1. 无量纲化后感染者与恢复者动态变化

图1可得,系统的地方病平衡点 E * 受时滞影响初期震荡,最终趋于稳定,地方病平衡点 E * 是局部渐近稳定的,疾病发展为地方病。

其他参数不变的前提下,取 ε=1.2> ε 0 时,系统稳定性发生改变,出现图2所示的持续周期震荡。

Figure 2. Dynamics of infected and recovered individuals after non-dimensionalization

2. 无量纲化后感染者与恢复者动态变化

图2可得,系统的地方病平衡点 E * 受时滞影响存在Hopf分支,模型失去稳定性,出现持续性周期震荡,地方病平衡点 E * 不稳定,系统存在周期解。

Figure 3. Bifurcation diagram in the x-y phase plane with respect to parameter ε=1.2

3. ε=1.2 时,系统在 x-y 相平面上分支

图3 ε=1.2 时, x-y 相平面上分支的情形,可以看出系统存在周期解,即产生了Hopf分支。验证结果与理论分析的结论相吻合。

5. 结论

本文构建并分析了一类具有饱和发生率的时滞SIR传染病模型。研究了时滞因素在该模型下的动力学分析,通过研究分析可以看出,当时滞 ε 没有经过临界时滞点 ε 0 时,通过数值模拟得到恢复者和感染者的人群密度趋于正平衡态,此时地方病平衡点 E * 是局部渐近稳定的,可以使疾病得到控制,当等于临界时滞点时,系统在平衡点处产生了Hopf分支,当时滞 ε 经过了临界时滞点 ε 0 时,通过数值模拟得到两人群出现周期震荡,此时系统不再稳定,这种情况下疾病难以得到控制。

综上,通过对模型的动力学分析表明,隔离延迟 ε 与临界时滞点 ε 0 的相对关系是决定疾病传播的关键阈值条件。在疾病防控中,严格将隔离延迟控制在临界值之内,在稳态阶段可以通过对人群聚集的管控和传染病传播途径的阻断等方法以此来扩大容错空间,当超过临界值时,应立即启动应急方案,通过集中隔离等方法来提高对高风险人群的管控,使疾病的传播重回稳定状态。

参考文献

[1] 马知恩, 周义仓, 王稳地, 等. 传染病动力学的数学建模与研究[M]. 北京: 科学出版社, 2004.
[2] Dietz, K. and Heesterbeek, J.A.P. (2002) Daniel Bernoulli’s Epidemiological Model Revisited. Mathematical Biosciences, 180, 1-21. [Google Scholar] [CrossRef] [PubMed]
[3] Kermack, W. and Mckendrick, A. (1991) Contributions to the Mathematical Theory of Epidemics—II. The Problem of Endemicity. Bulletin of Mathematical Biology, 53, 57-87. [Google Scholar] [CrossRef
[4] Xu, R. and Ma, Z. (2009) Stability of a Delayed SIRS Epidemic Model with a Nonlinear Incidence Rate. Chaos, Solitons & Fractals, 41, 2319-2325. [Google Scholar] [CrossRef
[5] 杨文杰, 王国际, 龚诗琪, 等. 具有时滞SIR模型的Hopf分支分析[J]. 许昌学院学报, 2023, 42(2): 7-10.
[6] 孔建云, 刘茂省, 王弯弯. 一类带有时滞的SIR模型的稳定性及分支分析[J]. 河北工业科技, 2017, 34(3): 167-171.
[7] Osvaldo, O., Joel, R.C., Cruz, V.D.L., et al. (2015) On the Existence and Construction of Dulac Functions. Value in Health, 13, A429.
[8] Xiao, D. and Ruan, S. (2007) Global Analysis of an Epidemic Model with Nonmonotone Incidence Rate. Mathematical Biosciences, 208, 419-429. [Google Scholar] [CrossRef] [PubMed]
[9] 崔玉美, 陈姗姗, 傅新楚. 几类传染病模型中基本再生数的计算[J]. 复杂系统与复杂性科学, 2017, 14(4): 14-31.