王定宇, 周少波

(1.華中科技大學數(shù)學與統(tǒng)計學院, 湖北武漢 430074;2.密歇根大學統(tǒng)計學系, 密歇根安娜堡 48104)

1.引言

長期以來傳染病的傳播和流行給人類健康帶來了極大威脅, 例如曾在中世紀流行于歐洲的黑死病, 如今的新型冠狀病毒肺炎等, 都對人類的健康和經(jīng)濟的發(fā)展帶來了負面的影響.為研究1665-1666年的黑死病在英國倫敦的流行規(guī)律以及1906年瘟疫在印度孟買的流行規(guī)律,Kermack和McKendrick在1927年首先建立了SIR倉室傳染病模型[1].他們的工作為SIR模型的建立奠定了基礎(chǔ).SIRS和SEIR等傳染病模型都是在SIR模型的基礎(chǔ)上提出來的.SIR模型的具體數(shù)學表達形式如下

其中S(t)表示易感者(Susceptible), I(t)表示患病者(Infectious), R(t)表示康復者(Recovered),β表示日接觸率, μ表示日恢復率.

由于人類對某些疾病的免疫時間是有限的, 若考慮康復者會重新轉(zhuǎn)換成易感者喪失免疫力, 人們提出了如下SIRS模型

其中γ表示喪失免疫率.

由于現(xiàn)代人對疾病防控意識的增強, 心理作用在傳染病治愈過程中的作用越來越明顯.同時在絕大多數(shù)傳染病模型動力學分析的文獻中, 都只考慮了系統(tǒng)內(nèi)總?cè)丝诤愣ǖ奶厥馇闆r.因此朱晶[2]假定系統(tǒng)內(nèi)總?cè)丝诘脑鲩L符合Logistic人口增長模型, 提出了如下基于心理作用的SIR傳染病模型

其中N(t)為t時刻的人口總量; b為自然出生率; d為自然死亡率; r =b ?d為內(nèi)稟增長率; K為環(huán)境容納量; α為心理作用系數(shù)(即采取相應的預防控制措施影響發(fā)病率), 且假設(shè)1/K < α < 1;β為傳染率; μ為自然恢復率; υ為獲得終身免疫率.除內(nèi)稟增長率r以外的所有參數(shù)均為正常數(shù),參數(shù)r的大小由自然出生率b和自然死亡率d共同決定.證明了當閾值<1時, 疾病將會消亡;當閾值>1時, 疾病將會持久存在.

顯然, 在模型(1.3)中, 當心理作用系數(shù)α的數(shù)值相對較大時, 單位時間內(nèi)由易感者轉(zhuǎn)換成感染者的人數(shù)會降低, 傳染病的傳染將會得到抑制; 當α的數(shù)值較小時, 單位時間內(nèi)會有更多的易感者轉(zhuǎn)換成感染者, 則傳染病的傳染趨勢將會擴大.因此, 心理作用系數(shù)α將會對疾病的傳播趨勢產(chǎn)生深刻的影響.

基于此, 李文軒等人[3?5]提出了人口增長滿足Logistic方程的基于心理作用的SIRS模型

其中θ表示治愈率; δ表示喪失免疫率; υ表示獲得暫時免疫率.

實際上, 我們在進行參數(shù)估計(如β)時會不可避免的存在一些誤差, 即參數(shù)β的估計通常是“均值”+“誤差項”, 若假定誤差項服從正態(tài)分布, 并將其標準偏差看作是噪聲強度, 且該噪聲強度可能依賴于易感人群S(t)及染病人群I(t), 則β可由β+σ ˙ B(t)替代, σ >0代表噪聲強度,=dB(t)/dt為白噪聲(即B(t)是布朗運動).因而(1.4)相應的隨機擾動模型為

將模型(1.5)中的三個方程相加有

因此模型中的總?cè)丝谠鲩L滿足Logistic方程.

近些年來, 許多學者研究了類似的隨機擾動傳染病模型[6?7].杜金姬等[6]研究了一類具有l(wèi)ogistic增長的隨機SIRS傳染病模型, 證明了模型具有唯一的全局正解并指出模型的正解是穩(wěn)定的馬爾可夫過程, 同時也給出使模型中疾病滅絕的充分條件.朱玲等[7]建立了一類具有Logistic增長的隨機SIVR傳染病模型, 證明了隨機系統(tǒng)正解的存在唯一性, 討論了隨機模型的解在與其對應的確定性模型的無病平衡點附近的漸近行為, 并指出白噪聲強度σi對隨機系統(tǒng)的影響, 認為易感者人數(shù)變化所受到隨機擾動的強度σ1會對整個隨機系統(tǒng)產(chǎn)生決定性的作用.

本文將通過構(gòu)造Lyapunov函數(shù), 利用It?o引理, 強大數(shù)定理和停時等隨機分析理論證明模型(1.5)全局正解的存在唯一性, 并給出使疾病滅絕和持久的充分條件.其次, 本文將考慮時滯對系統(tǒng)(1.5)的影響, 證明基于心理作用的時滯隨機SIRS傳染病模型全局正解的存在唯一性.最后將應用Euler方法和Milstein方法進行數(shù)值模擬, 驗證本文建立的結(jié)論并進行分析.

2.基于心理作用的隨機SIRS傳染病模型

進而

令k →∞, 有∞> V(S(0),I(0),R(0))+AT = ∞.由此出現(xiàn)矛盾, 假設(shè)不成立.因此有τ∞=∞ a.s., 從而τe=∞ a.s.證畢.

Ⅱ 滅絕性

其中

于是

其中ηj是相互獨立且服從標準正態(tài)分布N(0,1)的隨機變量.

相應的, 系統(tǒng)(1.5)對應的Milstein高階近似解Sj,Ij,Rj滿足如下方程組:

其中ξj是相互獨立且服從標準正態(tài)分布N(0,1)的隨機變量.

若取初始時刻易感者人數(shù)250萬人, 感染者人數(shù)150萬人, 免疫者人數(shù)5萬人, 即S(0) =250,I(0) = 150,R(0) = 5; 取傳染率β = 0.05, 自然出生率b = 0.5, 環(huán)境容納量K = 500, 心理作用系數(shù)α = 0.1, 自然死亡率d = 0.2, 自然恢復率μ = 0.05, 治愈率θ = 0.2, 喪失免疫率δ =0.01, 獲得暫時免疫率υ =0.1, 隨機擾動σ =0.02.由于

?0.001=2α[β ?α(d+μ+θ+υ)]<σ2=0.004<βα=0.005,

因此在該參數(shù)下系統(tǒng)存在唯一的全局正解且滿足定理2.2的條件.又由定理2.2易得?0.07 < 0 a.s.故疾病最終將會以指數(shù)形式滅亡.同時通過計算可得= 0.873 < 1.圖1中數(shù)值模擬的結(jié)果映證了此時疾病將會滅絕.

圖1 σ =0.02,<1時的數(shù)值模擬圖

相反的, 若提高傳染率β = 0.1, 降低心理作用系數(shù)α = 0.05, 保持模型中的其他參數(shù)不變,則由計算可得

圖2展示了此時的數(shù)值模擬結(jié)果, 表明此時疾病不再滅絕而是持久存在的.

圖2 σ =0.02,>1時的數(shù)值模擬圖

更進一步, 若保持模型中參數(shù)不變, 僅將σ增大至0.1, 經(jīng)計算可得

在某些情況下, 隨機性模型中數(shù)值模擬的結(jié)果顯示疾病將會滅亡, 但在其對應的確定性模型中數(shù)值模擬的結(jié)果顯示疾病將會持續(xù).因此隨機擾動會對模型產(chǎn)生極大的影響.相比于確定性系統(tǒng)(2.2), 隨機系統(tǒng)(1.5)中疾病滅絕的條件會相對弱.因此當系統(tǒng)(1.5)中的參數(shù)可以使疾病滅絕時, 同樣的參數(shù)并不一定能保證系統(tǒng)(2.2)中的疾病滅絕.數(shù)值模擬圖3(d)刻畫了這一現(xiàn)象.

圖3 σ =0.1,<1時的數(shù)值模擬圖

值得注意的是步長?t的取值會對系統(tǒng)解的路徑產(chǎn)生影響.當步長?t = 2 × 10?5時,S(t),I(t),R(t)的波動要比步長?t = 0.02時更加顯著.這說明越小的步長越能反應出在一個細微的時間段內(nèi)系統(tǒng)解的變化情況, 而較大的步長可能會掩蓋了在該時間段內(nèi)系統(tǒng)解的某些顯著變化.因此在計算機算力允許的條件下, 我們應該盡可能地選擇較小的步長?t.這在處理白噪聲強度σ值較大的隨機微分系統(tǒng)時更加的重要, 因為較大的σ值會使系統(tǒng)的擾動更加的劇烈.最后通過比較不難發(fā)現(xiàn)歐拉方法與Milstein方法在處理這類隨機微分方程的效果上并沒有顯著的差異, 也沒有出現(xiàn)系統(tǒng)解不收斂的情況.而在確定性系統(tǒng)中, 由于白噪聲強度σ = 0, 因此Milstein方法退化成為歐拉方法, 故在這類系統(tǒng)中這兩種方法的遞推式是一樣的.

3.基于心理作用的時滯隨機SIRS傳染病模型

由于大部分傳染性疾病具有潛伏期, 因此當易感者被傳染而轉(zhuǎn)換成感染者之后并不會馬上具有傳染性.故在t時刻受到感染的易感者并不是被t時刻的感染者感染的, 而實際上是被t ?τ時刻的感染者所感染的, 只是t ?τ時刻的感染者在(t ?τ)時刻沒有表達出傳染性.這體現(xiàn)了疾病傳染的滯后性.在具有時滯的系統(tǒng)中, 時滯可能會影響到系統(tǒng)的穩(wěn)定性, 引發(fā)系統(tǒng)的周期震蕩, 因此時滯隨機微分方程的動力學行為會更加復雜.考慮系統(tǒng)(1.5)具有時滯的情況,建立模型如下

且S(t)=N(t)?I(t)?R(t), 因此可以考慮上述方程組的等價模型

本節(jié)將證明系統(tǒng)(3.3)具有唯一的全局正解, 并繪制S(t),I(t),R(t)在不同時滯τ下對應的路徑圖.

Ⅰ 正解的存在性和唯一性

可以計算

其中A是獨立于N,I,R,t的正常數(shù).因此

對上式兩端從0到τk∧T積分, 并取期望可得

Ⅱ 數(shù)值模擬

在這一部分我們將使用高階Milstein方法對系統(tǒng)的解進行模擬分析.

取初始時刻易感者人數(shù)250萬人, 感染者人數(shù)150萬人, 免疫者人數(shù)5萬人, 即S(0) =250,I(0) = 150,R(0) = 5; 取傳染率β = 0.05, 自然出生率b = 0.5, 環(huán)境容納量K = 500,心理作用系數(shù)α = 0.1, 自然死亡率d = 0.2, 自然恢復率μ = 0.05, 治愈率θ = 0.2, 喪失免疫率δ =0.01, 獲得暫時免疫率υ =0.1, 隨機擾動σ =0.03, 分別取不同的τ值.則系統(tǒng)(3.1)的模擬結(jié)果如圖4所示.

圖4 T =500,?t=0.02時的數(shù)值模擬圖

可以看到疾病最終都將滅絕.但疾病滅絕所需要的時間隨著時滯τ的增大而增大, 這說明當傳染病的潛伏期較長時, 疾病往往需要更長的時間才能夠滅絕.因此具有潛伏期的疾病將會更加的狡猾, 需要政府采取更長時間的隔離措施才能夠有效地預防疾病的傳播.COVID-19就具有較長的潛伏期, 我國政府也是根據(jù)該傳染病的這種特性制定了入境人員14+7+7天的隔離措施, 收到了較好的成效.圖4的(c)和(d)反映出感染者和易感者人數(shù)的變化具有波動性, 且振蕩頻率隨著時滯τ的增大而增大.這說明當疾病的潛伏期較長時, 疾病的變化趨勢是具有波動性的, 即疫情會出現(xiàn)反復, 這提醒決策部門在分析疫情時不能盲目樂觀, 要考慮到這樣的波動特性, 做好打持久戰(zhàn)的準備, 時刻緊繃, 防患未然.COVID-19在我國的傳染趨勢就具有這樣的波動性.復旦大學附屬華山醫(yī)院感染科主任張文宏也曾在2020年接受采訪的時候表示要做好應對疫情反復, 甚至是第二波疫情的準備.

4.結(jié)論

本文研究了一類具有心理作用的隨機SIRS傳染病模型的動力學性質(zhì), 即模型解的存在唯一性, 滅絕性以及持久性.證明了當系統(tǒng)的白噪聲強度較大或閾值≤1且白噪聲強度不大時, 疾病將會滅絕; 當閾值> 1時, 對應系統(tǒng)的解將會持續(xù).本文還研究了時滯對模型的影響.同時, 我們應用歐拉方法和Milstein高階法模擬了不同參數(shù)下S(t),I(t),R(t)的路徑圖, 證實了本文建立的結(jié)論.