赵润东, 孙梅慈, 刘启明
陆军工程大学石家庄校区 军政基础系, 石家庄 050003
传染病历来是危害人类健康的大敌, 为了遏制疾病传播, 许多学者通过建立数学模型来研究其传播过程, 其中主要使用的是“仓室”(Compartment)模型. 1927年, 文献[1]建立了著名的SIR传染病仓室模型. 其后, 经过许多学者的不断研究, 建立了适用不同疾病的传染病模型, 如SIS[2], SEIR[3], SEIRS[4]等.
遏制疾病的传播, 我们通常采用疫苗接种和隔离两种方法. 但是针对新出现的传染病, 疫苗的研发和生产往往需要很长时间, 因此在疾病传播初期最为有效的方法就是对人群进行隔离[5]. 1995年, 文献[6]首次在传染病模型中考虑隔离的影响, 建立了SIQR模型; 2002年, 文献[7]在随机网络传染病模型中加入隔离项, 建立并研究了SIQS传染病模型. 上述研究都是基于随机网络研究的, 其特点是每个个体是均匀接触的.
然而, 文献[8]发现现实中大多数网络的节点分布符合无标度性(异质性), 也就是服从幂律分布p(k)=Ck-γ(2<γ≤3), 因此基于异质复杂网络来建立模型就更加贴合实际. 2001年, 文献[9]首次在无标度网络上对一类SIS传染病模型进行了研究. 此后, 许多学者开始研究复杂网络上的传染病动力学. 另一方面, 现实中许多传染病当前的传播状态会受到过去状态的影响, 因此, 建立时滞传染病模型就更具有现实意义, 其中时滞可以用来描述病人的平均感染周期、 潜伏周期、 免疫周期和隔离周期等[10]. 近期许多学者将网络的无标度性和时滞结合在一起研究传染病模型, 取得了丰富成果. 2012年, 文献[11]建立了时滞SEIRS网络传染病模型, 其中时滞代表平均免疫周期. 2018年, 文献[12]建立并研究了时滞SEIR网络传染病模型, 其时滞代表疾病的平均潜伏周期. 2019年, 文献[13]研究了复杂网络上一类新的时滞SIS模型, 其时滞代表病人的平均感染周期. 但是鲜有人在网络上用时滞表示隔离周期来建立数学模型对传染病动力学性态进行研究.
根据以上分析, 本文基于无标度网络建立一类新的具有时滞的SIQR传染病模型, 其中时滞代表平均隔离周期. 通过泛函微分方程稳定性理论, 研究了该模型的动力学行为, 得到疾病传播的基本再生数, 分析了平衡点的全局稳定性, 并通过数值模拟验证了研究结果的正确性.
假设总人群的接触网络是一个无标度网络, 一个节点表示一个人, 网络上的连边表示人与人之间的接触. 我们作如下假设:
1) 整个网络的出生率和自然死亡率分别为A和d, 并且出生的个体都为易感染者. 依据文献[14], 添加和删除节点和边在网络中只占很小的比例, 对网络结构的改变很小, 因此可以假设网络上的总节点数N是不变的, 是静态的, 也就是A=d.
2) 网络上的人分为4类: 易感染者S(Susceptible)、 感染者I(Infected)、 隔离者Q(Quarantine)、 恢复者R(Recovered). 令Sk(t),Ik(t),Qk(t)和Rk(t)分别代表度为k的易感染者、 感染者、 隔离者和恢复者在t时刻的相对密度. 标准化后,Sk(t)+Ik(t)+Qk(t)+Rk(t)=1.
3) 每个易感染者S都有概率λ(k)(与节点的度k有关)被感染者感染, 成为感染者I. 感染者I有概率δ被隔离, 成为隔离者Q. 感染者I同时有概率γ恢复健康, 成为恢复者R. 隔离者Q经过τ时间的隔离治疗后, 成为恢复者R.
图1 时滞SIQR模型的仓室框图
基于以上假设, 可建立SIQR模型的仓室框图(图1), 对应的微分方程系统如下:
(1)
其中λ(k)为感染率, 其形式一般有如下两种[14]:λk,λc(k).Θ(t)代表一个度为k的易感染者每次接触感染者的概率[14], 其形式如下:
(2)
给出系统(1)的初始条件
(3)
其中
由泛函微分方程的基本理论[18]易知在初始条件(3)下, 系统(1)存在唯一解, 并且当t≥0时, 系统(1)存在唯一正解. 同时, 不难得出区域Ω是系统(1)的正向不变集, 本文将在Ω内讨论系统(1)的性态.
(4)
建立模型后, 需得出疾病的基本再生数R0[18], 即单位病程内一个病人所传染的人数. 当R0<1时, 一个病人在单位病程能传染的人数小于1, 疾病将自然消失, 不会流行; 当R0>1时, 一个病人在单位病程能传染的人数大于1, 疾病将持久存在, 成为流行病.
定理1令
1) 系统存在无病平衡点E0
2) 当R0>0 时, 系统(1)存在地方病平衡点E*
证由系统(1)中的方程组, 不难得出系统(1)始终存在无病平衡点E0
接下来, 假设系统存在地方病平衡点E*,
则E*满足系统(1)
(5)
其中,
(6)
解方程组(5), 得到
(7)
(8)
显然,Θ*=0是平凡解. 当Θ*≠0时对(8)式两边同除以Θ*, 研究函数
因此, 我们定义基本再生数R0如下
综合上述分析, 当R0>1时, 系统(1)存在地方病平衡点E*.
注1由R0的表达式得出,R0与出生率A, 感染率λ(k)和非线性传染系数φ(k)正相关, 与死亡率d, 恢复率γ和隔离率δ负相关.
定理2当R0<1时, 系统(1)的无病平衡点E0全局渐近稳定.
计算V(t)沿系统(1)的导数可得
定理3当R0>1时, 系统(1)的地方病平衡点E*全局渐近稳定.
证注意到, 系统(1)中的前两个方程不涉及Q(t)和R(t), 由此, 只需考虑如下系统:
(9)
其中
联合(9)式和(5)式, 得到
考虑如下函数
(10)
U(t)沿系统(9)求导, 得
(11)
定义函数
H(x)=-x+lnx,G(x)=x-1-lnx
注意到当x>0时,G(x)≥0, 当且仅当x=1时, 等号成立. 则(11)式可化为
(12)
考虑如下两个矩阵
v=(v1,v2, …,vn)=(C11,C22, …,Cnn)
即
(13)
现在定义Lyapunov函数
注3当A=0,d=0,λ(k)=λk,φ(k)=k, 并取时滞τ=ε-1, 时滞微分系统(1)简化为文献[22]中的常微分系统(1). 此时基本再生数的表达式R0=λ〈k2〉·[〈k〉(γ+δ)]-1, 结论与文献[22]一致.
取λ=0.1,A=0.1,d=0.1,τ=3,δ=0.1,γ=0.03, 则R0≈0.774 1<1.I40(t),I80(t),I120(t),I160(t)和I(t)随时间t的变化趋势见图2, 可以看出当R0<1时, 疾病逐渐消失. 取λ=0.2,A=0.1,d=0.1,τ=3,δ=0.1,γ=0.03, 则R0≈1.548 2>1.I40(t),I80(t),I120(t),I160(t)和I(t)随时间t的变化趋势见图3, 可以看出当R0>1时, 疾病持续存在, 并且逐渐趋向到一个稳定值. 图2和图3分别验证了定理2和定理3.
最后, 我们进行参数敏感性分析, 用PRCC(偏秩相关系数)检测基本再生数R0对于参数的依赖性. 取样本空间n=1 200, 计算6个影响R0参数的PRCC值. 如图4所示,λ和A对R0有正影响,d,γ,δ对R0负影响, 而R0对τ不敏感. 图4验证了注1 中对R0表达式的说明. 因此, 增大隔离率δ, 提高恢复率γ, 降低感染率λ可以控制疾病传播.
图4 R0关于不同参数的PRCC值
为了研究隔离周期对传染病的影响, 本文基于无标度网络建立了具有时滞的SIQR传染病模型, 其中时滞代表平均隔离周期. 通过微分方程定性与稳定性理论, 得到了疾病传播的基本再生数. 通过构造Lyapunov函数, 证明了当R0<1时, 无病平衡点是全局渐近稳定的; 当R0>1时, 地方病平衡点是全局渐近稳定的. 最后, 对模型进行数值模拟, 验证了结论的正确性, 并对不同参数进行了敏感性分析. 研究结果表明隔离周期的长短不影响易感染者和感染者的最终人数, 但是影响隔离者和恢复者的最终人数. 基于R0的表达式以及不同参数的PRCC值, 得出控制疾病传播的有效方法为增大隔离率, 提高恢复率和降低感染率.