具有非線性恢復(fù)率的隨機(jī)SIQR模型的閾值動力學(xué)

孫乾乾，趙 宇，譚德君，張樹文

(1.廈門工學(xué)院數(shù)據(jù)科學(xué)與智能工程學(xué)院，福建 廈門 361021；2.集美大學(xué)理學(xué)院，福建 廈門 361021)

0 引言

疫情的暴發(fā)對人們的生活產(chǎn)生了深遠(yuǎn)的影響[1]?？v觀歷史，也有許多其他大型流行疾病對人類的生命健康構(gòu)成嚴(yán)重威脅，如1918年西班牙的大流感[2]、歐洲20世紀(jì)中世紀(jì)的大瘟疫[3]等。由于無法及時生產(chǎn)有效的疫苗以及病毒的快速變異，隔離被認(rèn)為是控制傳染病傳播的有效方法[4-7]。此外，有效治療是保障人類生命安全和減少疾病傳播的重要方法[8]。但在疾病暴發(fā)期間，醫(yī)護(hù)人員的數(shù)量、醫(yī)療設(shè)備、醫(yī)院床位和藥品的數(shù)量都是非常有限的[9-11]。因此，為了刻畫醫(yī)療資源對疾病的影響，文獻(xiàn)[11]提出非線性恢復(fù)率γ(b，I)=γ0+(γ1-γ0)b/(I+b)，其中：γ1、γ0分別是最大和最小的恢復(fù)率；I是染病者的數(shù)量；b是醫(yī)療資源。本文將在文獻(xiàn)[5,11]研究的基礎(chǔ)上，提出SIQR模型為

(1)

其中：S、Q、R分別是易感者、隔離者、恢復(fù)者的數(shù)量；N=S+I+Q+R是總?cè)丝诘臄?shù)量；A是新的易感者的輸入率；β是疾病的傳播速率；μ是自然死亡率；δ是隔離率；α1和α2分別是染病者和隔離者的因病死亡率；ε是隔離者的恢復(fù)率。

此外，大量實驗和流行病學(xué)研究表明，病毒的活性、感染力與周圍空氣的濕度、溫度等密切相關(guān)[12-15]。因此，在流行病的傳播過程中，疾病的傳播速率不可避免會受到環(huán)境隨機(jī)波動的影響[16-19]。本文假設(shè)模型(1)的疾病傳播速率β受到白噪聲的干擾，即βdt→βdt+σdB(t)，然后得到的模型為

(2)

其中：σ是白噪聲的波動強(qiáng)度；B是一維標(biāo)準(zhǔn)布朗運(yùn)動。

1 主要結(jié)果及證明

在流行病模型的研究中，基本再生數(shù)R0是一個重要的量，它表示發(fā)病初期在一個全部是易感人群的環(huán)境中，一個染病者在平均患病期內(nèi)所傳染的人數(shù)。利用下一代矩陣可得到模型(1)的基本再生數(shù)為

R0=β/(γ1+δ+μ+α1)。

(3)

由此可知定理1。

定理1 模型(1)總是存在無病平衡點E0=(A/μ，0，0，0)。如果R0<1，則模型(1)的無病平衡點E0是局部漸近穩(wěn)定的；如果R0>1，則E0是不穩(wěn)定的。

定理2 如果R0>1，則模型(1)存在唯一的地方病平衡點E*(S*，I*，Q*，R*)。

證明由于模型(1)的地方病平衡點E*(S*，I*，Q*，R*)滿足方程

(4)

容易驗證I*滿足方程

α1-α2I*2-α3I*=0，

(5)

其中：α1=Ab(β-(δ+μ+α1+γ1))；α2={β-[1-(μ+ε)/(ε+μ+α2)]δ-α1}(δ+μ+α1+γ0)；α3={β-[1-(μ+ε)/(ε+μ+α2)]δ-α1}(δ+μ+α1+γ1)b-A[β-(δ+μ+α1+γ0)]。當(dāng)R0>1時，有α1>0且α2>0。由二次方程求根公式可知，方程(5)存在唯一的正解，即當(dāng)R0>1時，模型(1)存在唯一的地方病平衡點，定理2證畢。

顯然，當(dāng)R0<1時，方程(5)可能存在兩個正解，即模型(1)可能存在兩個正平衡點。因此，模型(1)可能存在后向分岔。接下來應(yīng)用中心流形定理來證明模型(1)在R0=1時發(fā)生后向分岔。

通過計算可得

(6)

且Λ2=(ε+μ+α2)/δ>0。然后，根據(jù)文獻(xiàn)[20]中的定理4.1可得定理3。

定理3 如果Λ1>0，則模型(1)在R0=1處發(fā)生后向分岔。

為了方便理解，給出例1。

例1 首先取如下參數(shù)：A=30，μ=0.1，γ0=0.009 8，γ1=0.239 4，b=0.072，δ=0.104 6，α1=0.001 7，ε=0.239 4，α2=0.001 1，且令β在區(qū)間[0.133 71,0.467 985]內(nèi)變化，使得R0在區(qū)間[0.30,1.05]內(nèi)變化。由式(6)可得，Λ1=1.219 07>0。因此，根據(jù)定理3可知，模型(1)在R0=1處發(fā)生后向分岔(見圖1)。

此外，令β=0.231 76，此時R0=0.52。根據(jù)定理3可知，模型(1)的解要么被地方病平衡點E*(282.160 667 5，8.179 298 012，2.512 641 915，6.980 705 454)所吸引，要么被無病平衡點E*(300，0，0，0)所吸引。這意味著，R0不足以完全確定基本的動力學(xué)行為，還需要考慮初始條件等其他流行病學(xué)因素(見圖2)。

定理4的證明比較簡單，詳細(xì)的證明過程可參考文獻(xiàn)[21]。

證明首先證明，對于任意的>0，總是存在ξ>0，使得對于任意的(S(0)，I(0)，Q(0)，R(0))∈Γ∩Oξ，其中，Oξ=(A/μ-ξ，A/μ)×(0，ξ)3，有

(7)

定義：f1(S，I，Q，R)=A-βSI/N-μS，g(S，I，Q，R)=σSI/N；f2(S，I，Q，R)=βSI/N-(γ(b，I)+δ+μ+α1)I，λ=β-(γ1+δ+μ+α1)。f1(S，I，Q，R)、f2(S，I，Q，R)和g(S，I，Q，R)在無病平衡點E0=(A/μ，0，0，0)的偏導(dǎo)數(shù)分別為：(?f1)/(?S)(E0)=-μ；(?f1)/(?I)(E0)=-β； (?f1)/(?Q)(E0)=0；(?f1)/(?R)(E0)=0；(?f2)/(?S)(E0)=0；(?f2)/(?I)(E0)=λ；(?f2)/(?Q)(E0)=0；(?f2)/(?R)(E0)=0；(?g)/(?S)(E0)=0；(?g)/(?I)(E0)=σ；(?g)/(?Q)(E0)=0；(?g)/(?R)(E0)=0。因此，f1(S,I,Q,R)、f2(S，I，Q，R)和g(S，I，Q，R)在無病平衡點E0附近的泰勒展開分別為

f1=-μS-βI+o；f2=βI-(γ1+δ+μ+α1)I+o；g=σI+o，

(8)

接下來證明，對于任意的(S(0)，I(0)，Q(0)，R(0))∈Γ，模型(2)的解(S(t),I(t),Q(t),R(t))總會在某一時刻進(jìn)入?yún)^(qū)域?！蒓ξ。定義停時τξ=inf{t≥0：S(t)≥A/μ-ξ}，并定義函數(shù)V1(S，I，Q，R)=C1-(1+S)C2，其中C1和C2是兩個正常數(shù)，并且滿足如下條件：1)對于任意的(S，I，Q，R)∈Γ，都有V1>0；2)C2>2；3)對于任意的S∈(0，A/μ-ξ]，有(1+S)(A-βSI/N-μS)+0.5σ2(C2-1)S2I2/N2≥0.5ξ。

然后，根據(jù)Dynkin公式可得：E[V1(S(τξ∧t)，I(τξ∧t)，Q(τξ∧t)，R(τξ∧t))]≤V1((S(0)，I(0)，Q(0)，R(0)))-0.5ξE(τξ∧t)。令t→∞，并根據(jù)Fatou引理可得：E[V1(S(τξ)，I(τξ)，Q(τξ)，R(τξ))]≤V1((S(0)，I(0)，Q(0)，R(0)))-0.5ξE(τξ)。

注意到，對于任意的(S，I，Q，R)∈Γ，V1是有界的，因此有

E(τξ)<∞。

(9)

最后，根據(jù)強(qiáng)Markov性質(zhì)，由式(7)和式(9)可得，對任意的，對于任意的(S(0)，I(0)，Q(0)，R(0))∈Γ。定理5證畢。

證明定義函數(shù)H(S，I，Q，R)=M[1(S+I+Q+R)-lnI-2lnS]-lnS-lnQ-lnR-ln(N-A/(μ+α1+α2))-ln(A/μ-N)∶=MV1+V2，其中：V1=1(S+I+Q+R)-lnI-2lnS；V2=-[lnS+lnQ+lnR+ln (N-A/(μ+α+α))+ln (A/μ-N)]；1=(γ1+δ+μ+α1+0.5σ2)/A；2=(γ1+δ+μ+α1+0.5σ2)/μ；M是一個充分大的常數(shù)，使得

(10)

LV1≤-βS/N+(γ1+δ+μ+α1+0.5σ2)+2(-A/S+βI/N+μ+σ2I2/(2N2))+1A-1μN(yùn)

(11)

并且有

LV2≤-A/S+β+μ+0.5σ2-δI/Q+(ε+μ+α2)-εQ/R+μ-(α1I+α2Q)/(A/μ-N)-

(α1(S+Q+R)+α2(S+I+R))/(N-A/(μ+α1+α2))+μ+α1+α2+μ。

(12)

根據(jù)式(10)～式(12)可得：LV=MLV1+LV2≤-2+M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-A/S-δI/Q-εQ/R-α2I/(N-A/(μ+α1+α2))-α1I/(A/μ-N)。令ε是一個充分小的正常數(shù)，使得-2+M2β(μ+α1+α2)/A+0.5M2σ2(μ+α1+α2)22/A2≤-1，M2β(μ+α1+α2)/μ+0.5M2σ2(μ+α1+α2)2/μ2-min{A，δ，ε，α1，α2}/<-1。

然后將ΓD分成6個區(qū)域，即ΓD=D1∪D2∪D3∪D4∪D5∪D6，其中：D1={(S，I，Q，R)∈Γ：0

對于任意的(S，I，Q，R)∈D1，有

LV≤-2+M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2
≤-2+M2β(μ+α1+α2)/A+0.5M2σ2(μ+α1+α2)22/A2≤-1。

(13)

對于任意的(S，I，Q，R)∈D2，有

LV≤M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-A/S
≤M2β(μ+α1+α2)/μ+0.5M2σ2(μ+α1+α2)2/μ2-A/≤-1。

(14)

對于任意的(S，I，Q，R)∈D3，有

LV≤M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-δI/Q
≤M2β(μ+α+α)/μ+0.5M2σ2(μ+α+α)2/μ2-δ/2≤-1。

(15)

對于任意的(S，I，Q，R)∈D4，有

LV≤M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-εQ/R
≤M2β(μ+α1+α2)/μ+0.5M2σ2(μ+α1+α2)2/μ2-ε2/3≤-1。

(16)

對于任意的(S，I，Q，R)∈D5，有

LV≤M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-α2I/(N-A/(μ+α1+α2))
≤M2β(μ+α1+α2)/μ+0.5M2σ2(μ+α1+α2)2/μ2-α2/2≤-1。

(17)

對于任意的(S，I，Q，R)∈D6，有

LV≤M2β(μ+α1+α2)I/A+0.5M2σ2(μ+α1+α2)2I2/A2-α1I/(A/μ-N)
≤M2β(μ+α1+α2)/μ+0.5M2σ2(μ+α1+α2)2/μ2-α1/2≤-1。

(18)

由式(13)～式(18)可得，對于任意的(S(t)，I(t)，Q(t)，R(t))∈ΓD，有LV≤-1。因此，Markov過程(S(t)，I(t)，Q(t)，R(t))關(guān)于D是正常返的，定理6證畢。

2 數(shù)值模擬

根據(jù)C0VID-19的實際參數(shù)值來驗證并擴(kuò)展本文的理論結(jié)果。首先假設(shè)，武漢某一地區(qū)的常住人口總規(guī)模為10 000人，每個人的平均壽命為80 a，且每千人擁有5張病床，則自然死亡率μ=1/(365×80)=3.424 65×10-5，醫(yī)院病床數(shù)b=50。假定新的易感者的輸入率A=30，疾病的傳播速率β=1.5。其他的參數(shù)取自文獻(xiàn)[22]?？偟膮?shù)取值為：A=30，β=1.5，μ=3.424 65×10-5，γ0=0.009 8，γ1=1/14，b=50，δ=0.104 6，α1=0.001 7，ε=0.239 4，α2=0.001 1。

接下來，假設(shè)隨機(jī)模型(2)的初值(S(0)，I(0)，Q(0)，R(0))=(7 000，1 000，1 000，1 000)，并考慮白噪聲強(qiáng)度σ對模型(2)動力學(xué)的影響。

3 結(jié)論