何雪晴,韦煜明
(广西师范大学 数学与统计学院,广西 桂林 541006)
世界正面临着新型冠状病毒肺炎(COVID-19)带来的严重威胁[1],人类感染病毒后常见的体征有发烧、咳嗽、气促和呼吸困难等,在较严重的病例中,感染可以导致肺炎、严重急性呼吸综合征、肾衰竭,甚至死亡。这引起了许多学者的注意,他们通过数学建模对其进行研究[2-7]。例如,Wang等[5]建立了一个模型来估计武汉市的流行趋势,对假设预防和控制措施足以或不足以控制疫情进行了讨论。Hellewell等[6]构建了一个传播模型,发现在大多数情况下,高效地对接触者进行追踪和对染病者隔离,足以在三个月内控制新的新冠肺炎疫情爆发。Mandal等[7]考虑了新冠肺炎的数学模型,其中的人类群体被细分为了5个关于时间t的类别,即易感类S(t)、暴露类E(t)、隔离类Q(t)、感染类I(t)和恢复类R(t),假设当一个易受感染的人与一个暴露的人接触时,新冠肺炎传染病就会传播。该模型由5个一阶常微分方程组成,如下所示:
(1)
其中所有参数都是正数。A为易感人群的招募率,β为疾病的传播率,ρ1(0<ρ1<1)为易感人群中保持适当预防措施的比例,ρ2(0<ρ2<1)为暴露人群中对疾病传播采取适当预防措施的比例,α和b2分别为暴露类中被感染和隔离的比例,b1和c分别为隔离类转移到易感类和感染类的比例,η和v分别为感染类和暴露类的恢复率,μ为自然死亡率,δ为因新冠肺炎诱发的死亡率。
但是,考虑到传染病的传播具有随机性,生活中的环境白噪声无所不在,所以确定性的COVID-19传染病模型难以避免地会受到环境噪声的影响,这使得模型中涉及的参数往往会随着周围环境的波动而随机地在一些平均值附近波动,因此,有必要在确定性的COVID-19传染病模型(1)中引入随机波动。研究者通常用Gauss白噪声来刻画外部环境对系统的随机扰动[8-9],考虑到除了白噪声之外,传染病模型也可能受到彩色噪声的干扰,导致系统从一种环境状态切换到另一种环境状态[10],而且在一般情况下,环境状态之间的切换是无记忆的,下一次切换的等待时间遵循指数分布[11],所以这种切换可以用连续时间马尔可夫链{r(t),t≥0}来描述,然后根据环境中的温度、湿度、光照、降雨和风度等不同条件,将其分为N种不同状态,其状态空间记为={1,2,…,N}。因此,基于以上考虑,为了让模型更符合实际情况,提出了以下具有Markov切换的随机COVID-19传染病模型
(2)
其中σi(i=1,2,3,4,5)用于刻画白噪声的强度,Bi(t)(i=1,2,3,4,5)是标准的Brown运动。设在有限状态空间={1,2,…,N}中,Markov链 {r(t),t≥0}由转移率矩阵Γ=(γij)N×N生成,即
dx(t)=b(x(t),r(t))dt+σ(x(t),r(t))dB(t),x(0)=x0,r(0)=r,
其中B(·)是d维Brown运动,r(·)是右连续的Markov 链,b(·,·):n×→n,σ(·,·):n×→n×d且满足σ(x,k)σT(x,k)=(dij(x,k)),k∈。设V(·,k)是任意二次连续可微函数,线性算子L被定义为
由广义的Ito公式[12]可知,如果x(t)∈n,那么有
dV(x(t),r(t))=LV(x(t),r(t))dt+Vx(x(t),r(t))σ(x(t),r(t))dB(t)。
引理1[13]设M={Mt}t≥0是一个实值的连续局部鞅,那么有
引理2[14]对几乎所有的ω∈Ω,存在正常数t0=t0(ω),使得当t≥t0时,有
由引理2可知,集合X是模型(2)解的一个吸引域,下面要讨论的内容仅在此区域内研究。
定理1对任意给定的初始值(S(0),E(0),Q(0),I(0),R(0),r(0))∈X,模型(2)在t≥0时,存在唯一的全局解(S(t),E(t),Q(t),I(t),R(t)),并且该解以概率1位于Γ中。
证明由于模型(2)的系数满足局部Lipschitz条件,所以对任意给定的初始值(S(0),E(0),Q(0),I(0),R(0),r(0))∈Χ,模型(2)存在唯一的局部解(S(t),E(t),Q(t),I(t),R(t),r(t)),t∈[0,τe),其中τe表示爆炸时间。要证随机模型(2)存在唯一的全局解,只需证τe=+∞,a.s.。
P{τk≤T}≥ε,k≥k1。
(3)
V(S(t),E(t),Q(t),I(t),R(t),r(t))=(S(t)-1-lnS(t))+(E(t)-1-lnE(t))+(Q(t)-1-lnQ(t))+
(I(t)-1-lnI(t))+(R(t)-1-lnR(t))。
因为对任意的x>0,有x-1-lnx≥0,所以V正定。由Ito公式可知
dV(S,E,Q,I,R)=LV(S,E,Q,I,R)dt+σ1r(t)(S-1)dB1(t)+σ2r(t)(E-1)dB2(t)+
σ3r(t)(Q-1)dB3(t)+σ4r(t)(I-1)dB4(t)+σ5r(t)(R-1)dB5(t)。
(4)
≤Ar(t)+βr(t)(1-ρ1r(t))(1-ρ2r(t))E+5μr(t)+b1r(t)+b2r(t)+αr(t)+vr(t)+cr(t)+ηr(t)+δr(t)+
这里的K1是一个正常数,所以由(4)式整理可得
dV(S,E,Q,I,R)≤K1dt+σ1r(t)(S-1)dB1(t)+σ2r(t)(E-1)dB2(t)+σ3r(t)(Q-1)dB3(t)+
σ4r(t)(I-1)dB4(t)+σ5r(t)(R-1)dB5(t)。
(5)
对(5)式从0到τk∧T积分并取期望得
(6)
V(S(τk,ω),E(τk,ω),Q(τk,ω),I(τk,ω),R(τk,ω))≥k-1-lnk,
则
V(S(0),E(0),Q(0),I(0),R(0))+K1T≥[IΩk(ω)V(S(τk,ω),E(τk,ω),Q(τk,ω),I(τk,ω),R(τk,ω))]
其中IΩk表示Ωk的示性函数。令k→+∞,则有+∞>V(S(0),E(0),Q(0),I(0),R(0))+K1T=+∞,矛盾,故τe=+∞a.s.,定理1得证。
定理2设(S(t),E(t),Q(t),I(t),R(t))是初始值为
证明由模型(2)的第二个等式及Ito公式可知
对任意的k∈,根据引理2,得到
(7)
对(7)式两边从0到t积分并同时除以t,得到
(8)
由r(t)的遍历性可知
(9)
(10)
这一部分将利用Milstein方法[16]以及来自印度的马哈拉施特拉邦(印度西部邦)和德里的真实新冠肺炎数据,说明关于模型(2)的一个主要理论结果。
假设连续时间Markov链{r(t),t≥0}只有两个环境状态,即={1,2},平稳分布π=(0.5,0.5),选取初始值为(S(0),E(0),Q(0),I(0),R(0),r(0))=(10,8,5,3,1,1),结合表1中Mandal等[17]给出的参数值,对模型(2)进行两种情况模拟。
表1 来自印度的各参数取值
1)如果固定环境状态e∈,则模型(2)的阈值参数变为那么环境在状态1的情况下,取β1=3.5,A1=0.12,σ1=0.15,σ2=0.25,σ3=0.01,σ4=0.02,σ5=0.02,则根据定理2可知疾病将会灭绝,如图1所示。