具有分布時(shí)滯的毒品濫用防治模型的穩(wěn)定性分析

劉 風(fēng)

(中央司法警官學(xué)院 信息管理系,河北保定 071000;中央司法警官學(xué)院 戒毒康復(fù)研究中心,河北保定 071000)

§1 引言

按照《2020年中國毒情形勢(shì)報(bào)告》公布的數(shù)據(jù),盡管我國登記在冊(cè)的吸毒人數(shù)已連續(xù)三年持續(xù)下降,但是截至該年底全國仍有吸毒人員180.1萬名,占人口總數(shù)的0.13%,我國毒品濫用問題依舊不容樂觀.

在公共衛(wèi)生領(lǐng)域，數(shù)學(xué)模型一直是研究傳染病流行規(guī)律的有效工具之一.文獻(xiàn)[1-3]通過將海洛因的吸食擴(kuò)散與傳染病的流行進(jìn)行類比,打開了應(yīng)用傳染病模型研究海洛因吸食問題的大門.針對(duì)海洛因成癮現(xiàn)象,White 和Comiskey[4]利用常微分方程構(gòu)建了一個(gè)倉室傳染病模型,計(jì)算了模型的基本再生數(shù)和平衡點(diǎn),并對(duì)基本再生數(shù)進(jìn)行了敏感性分析,奠定了毒品濫用流行病模型的研究基礎(chǔ).Mulone 和Straughan[5]對(duì)文獻(xiàn)[4]進(jìn)行了補(bǔ)充,放松了系統(tǒng)流入率為常數(shù)的假設(shè),并證明了正平衡點(diǎn)的全局穩(wěn)定性.基于文獻(xiàn)[4-5],文獻(xiàn)[6]將吸毒人員劃分為輕度和重度兩個(gè)人群,建立了甲基苯丙胺流行病模型,分析了模型的動(dòng)力學(xué)性態(tài)并進(jìn)行了模擬仿真.

現(xiàn)實(shí)中吸毒者接受治療后往往又會(huì)出現(xiàn)復(fù)吸現(xiàn)象,從治療到復(fù)吸,狀態(tài)的轉(zhuǎn)移通常存在一個(gè)時(shí)間上的滯后.為了使模型能夠更好地刻畫現(xiàn)實(shí)情況,Liu和Zhang[7]將復(fù)吸時(shí)滯因素引入White和Comiskey的模型,建立了具有分布時(shí)滯的海洛因流行病模型,證明了無毒平衡點(diǎn)的全局穩(wěn)定性和地方病平衡點(diǎn)的局部穩(wěn)定性.Huang和Liu[8]進(jìn)而證明了文獻(xiàn)[7]中地方病平衡點(diǎn)的全局穩(wěn)定性.文獻(xiàn)[9-11]則考慮首次吸毒的時(shí)滯因素建立了相應(yīng)的海洛因流行病模型,并研究了模型平衡點(diǎn)的穩(wěn)定性.

接種疫苗是一種有效控制傳染病流行的預(yù)防措施.Parsamanesh和Erfanian[12]構(gòu)建了帶有疫苗接種措施的傳染病模型,證明了模型平衡點(diǎn)的全局穩(wěn)定性.Arino等[13]分析了疫苗接種措施對(duì)傳染病模型穩(wěn)定性的影響,研究結(jié)果表明當(dāng)疫苗并非完全有效時(shí),模型出現(xiàn)了后向分支現(xiàn)象.

針對(duì)毒品易感人群,加強(qiáng)毒品危害的宣傳教育,可以起到類似于接種疫苗的預(yù)防作用.與上述海洛因流行病模型只考慮治療措施不同,本文在預(yù)防和治療雙重措施基礎(chǔ)上引入復(fù)吸時(shí)滯影響因素,建立了一個(gè)具有分布時(shí)滯的毒品濫用防治模型,通過分析模型的動(dòng)力學(xué)性態(tài),討論了影響模型基本再生數(shù)的敏感性因素,研究結(jié)果為決策部門制定相關(guān)公共政策提供了理論依據(jù).

§2 模型構(gòu)建

毒品濫用人群的演化過程如圖1所示.S(t),I(t),T(t)和R(t)分別代表t時(shí)刻易感人群,未隔離吸毒人群,隔離治療人群和免疫人群的人數(shù).Λ是系統(tǒng)補(bǔ)充率.假設(shè)免疫人群不再吸毒,隔離治療人群不具有傳染性,β是未隔離吸毒人群與易感人群的有效接觸率,采用雙線性發(fā)生率.θ為針對(duì)易感人群進(jìn)行毒品危害宣傳教育的有效率,假定0＜θ ＜1,即宣傳教育并非完全有效.μ代表自然死亡率,δ1和δ2分別為未隔離吸毒人群和隔離治療人群的吸毒致死率,ε是未隔離吸毒人群的社區(qū)治療康復(fù)率,p為隔離治療率,φ是隔離治療康復(fù)率,σ(t,s)表示經(jīng)過時(shí)滯s到達(dá)t時(shí)刻的復(fù)吸人數(shù),根據(jù)圖1中各狀態(tài)轉(zhuǎn)移關(guān)系得到時(shí)滯微分系統(tǒng)

圖1 毒品濫用人群的演化過程

考慮到復(fù)吸時(shí)滯因人而異,設(shè)τ是時(shí)滯上限,時(shí)滯分布參數(shù)s ∈[0,τ],f為非負(fù)的連續(xù)函數(shù),則是時(shí)滯s的分布函數(shù).從t-s時(shí)刻開始接受治療,經(jīng)過時(shí)滯s到達(dá)時(shí)刻t未康復(fù)的幸存者為pI(t-s)e-(μ+δ2+φ)s,故t時(shí)刻的復(fù)吸人數(shù).系統(tǒng)所有參數(shù)非負(fù),初始條件為

其中,為將區(qū)間[-τ,0]映射到的所有連續(xù)函數(shù)組成的巴拿赫空間.基于系統(tǒng)解的連續(xù)性,T(0)應(yīng)滿足

由系統(tǒng)(1)中的第三個(gè)方程和式(3)有

Γ是正不變集,對(duì)于?t ＞0,滿足初始條件(2)的解及其ω極限集仍在可行域Γ內(nèi),系統(tǒng)(1)是適定的.

由于系統(tǒng)(1)的前三個(gè)方程獨(dú)立于R(t),而T(t)又完全取決于I(t),所以下面主要研究二維系統(tǒng)

§3 平衡點(diǎn)和穩(wěn)定性

3.1 無毒平衡點(diǎn)及其穩(wěn)定性

定理3.1系統(tǒng)(4)總是有無毒平衡點(diǎn)E0=S0,I0）,當(dāng)R0＞1時(shí),E0是不穩(wěn)定的,當(dāng)R0＜1時(shí),E0全局漸近穩(wěn)定.

證系統(tǒng)(4)在E0處的特征方程為[λ+(μ+θ)]g(λ)=0,其中

可見λ=-(μ+θ)＜0是一個(gè)特征值,其他特征值是方程g(λ)=0的根.

當(dāng)R0＞1時(shí),由和g(+∞)=+∞,可知方程g(λ)=0至少存在一個(gè)正根,所以E0是不穩(wěn)定的.

當(dāng)R0＜1時(shí),設(shè)g(λ)=0的根λ=a+bi,a,b ∈R.若a ≥0,將λ=a+bi代入方程g(λ)=0,可以推出

與R0＜1矛盾,所以a ＜0,E0是局部漸近穩(wěn)定的.

為了完整地證明定理3.1,考慮Lyapunov泛函

因?yàn)楹瘮?shù)v(z)=z-1-lnz,z ∈R+在z=1處取到最小值v(1)=0,所以V(t)≥0,當(dāng)且僅當(dāng)在E0處等式成立.下面證明是非正定的.首先計(jì)算V1對(duì)t的導(dǎo)數(shù)

因?yàn)棣?(μ+θ)S0,所以

其次計(jì)算V2對(duì)t的導(dǎo)數(shù)

合并式(5)和式(6)得到

因?yàn)镽0＜1,所以正定泛函V(t)關(guān)于時(shí)間t的導(dǎo)數(shù)非正定.設(shè)M1是的最大不變子集,若要使,則需滿足S(t)=S0,I(t)=I0,因此M1=E0,由E0的局部漸近穩(wěn)定性和LaSalle不變集原理可知E0全局漸近穩(wěn)定.

3.2 地方病平衡點(diǎn)及其穩(wěn)定性

R0＞1時(shí),除不穩(wěn)定的無毒平衡點(diǎn)E0外,系統(tǒng)(4)還有一個(gè)地方病平衡點(diǎn)

定理3.2當(dāng)R0＞1時(shí),系統(tǒng)(4)的地方病平衡點(diǎn)E*是局部漸近穩(wěn)定的.

證在E*處,系統(tǒng)(4)的特征方程為

由于λ的各項(xiàng)系數(shù)和常數(shù)項(xiàng)均大于零,故方程(7)的兩個(gè)根均具有負(fù)實(shí)部,E*是局部漸近穩(wěn)定的.可見E*不穩(wěn)定只能出現(xiàn)在τ＞0,且方程(7)的根從虛軸的左側(cè)穿過虛軸進(jìn)入到復(fù)平面的右半平面時(shí).不失一般性,設(shè)λ=wi,w ＞0是方程(7)的一個(gè)根,則有分離實(shí)部和虛部得到

定理3.3對(duì)于系統(tǒng)(4),當(dāng)R0＞1時(shí),地方病平衡點(diǎn)E*是全局漸近穩(wěn)定的.

證考慮Lyapunov泛函

顯然W(t)≥0,當(dāng)且僅當(dāng)在E*處等式成立.下面證明是非正定的.首先計(jì)算W1對(duì)t的導(dǎo)數(shù)

其次計(jì)算W2對(duì)t的導(dǎo)數(shù)

合并式(11)和式(12)得到

§4 數(shù)值模擬

以某市相關(guān)數(shù)據(jù)為基礎(chǔ)設(shè)定參數(shù)值

利用MatLab軟件模擬不同初始條件下未隔離吸毒人群的演化過程,驗(yàn)證平衡點(diǎn)的存在性和穩(wěn)定性.簡單起見,設(shè).R0=1時(shí),τ ≈1.3035.

首先考慮τ=15,即R0≈0.4173＜1時(shí),如圖2(a)所示,無論初始狀態(tài)如何,未隔離吸毒人群的人數(shù)最終都將趨近于0，模擬結(jié)果驗(yàn)證了無毒平衡點(diǎn)E0的存在性和全局漸近穩(wěn)定性.其次當(dāng)τ=1,β=1.3×10-7,即R0≈1.4207＞1時(shí),如圖2(b)所示,對(duì)于不同的初始值,未隔離吸毒人群的人數(shù)最終都趨近于唯一的正平衡點(diǎn),模擬結(jié)果表明系統(tǒng)唯一的地方病平衡點(diǎn)E*具有全局漸近穩(wěn)定性.

圖2 平衡點(diǎn)模擬結(jié)果

由于基本再生數(shù)R0是決定毒品濫用是否持續(xù)流行的一個(gè)關(guān)鍵閾值,因此有必要就R0對(duì)其參數(shù)進(jìn)行敏感性分析.下面分別計(jì)算R0關(guān)于參數(shù)β,θ和p的敏感系數(shù).

可以看出R0與β的變化方向是一致的,而與θ和p的變化方向相反,若要使R0＜1,可以通過降低β或者提高θ和p得以實(shí)現(xiàn).顯然相對(duì)于p,R0對(duì)β的變化更為敏感,圖3(a)中R0關(guān)于β和p的等值線也支持了這個(gè)結(jié)論.

R0對(duì)θ的相對(duì)變化率取決于各參數(shù)的具體取值,固定Λ=1×104,β=1×10-7,ε=0.005,φ=0.3,μ=δ1=δ2=0.004,設(shè)置θ ∈[0.02,0.08],p ∈[0.02,0.08],畫出R0關(guān)于θ和p的等值線圖.如圖3(b)所示，此時(shí),R0對(duì)θ的變化比對(duì)p的變化更為敏感.

圖3 R0的等值線圖

綜合上述分析,相對(duì)于p,β和θ是對(duì)R0更為敏感的影響因素,這意味著若通過降低R0減小吸毒人員規(guī)模,進(jìn)而消除毒品濫用現(xiàn)象,事前預(yù)防比事后治療是更為有效的措施.

§5 結(jié)論

通過引入復(fù)吸時(shí)滯影響因素,本文建立了一個(gè)具有分布時(shí)滯的毒品濫用防治模型,在計(jì)算模型基本再生數(shù)R0的基礎(chǔ)上,分析和驗(yàn)證了模型的平衡點(diǎn)及其全局漸近穩(wěn)定性,并就R0對(duì)其參數(shù)進(jìn)行了敏感性分析.當(dāng)R0＜1時(shí),模型具有全局漸近穩(wěn)定的無毒平衡點(diǎn),毒品濫用現(xiàn)象將最終消亡.當(dāng)R0＞1時(shí),模型具有唯一全局漸近穩(wěn)定的正平衡點(diǎn),毒品濫用現(xiàn)象將持續(xù)流行.相對(duì)于參數(shù)p,R0對(duì)參數(shù)β和θ的變化更為敏感,抑制毒品濫用的流行,事前預(yù)防比事后治療更為有效.這意味著若要根除毒品濫用現(xiàn)象,必須嚴(yán)厲打擊制毒販毒犯罪活動(dòng),深入實(shí)施毒品預(yù)防教育,尤其是青少年毒品預(yù)防教育工程,創(chuàng)新完善堵源截流機(jī)制.唯有這樣,才能不斷推動(dòng)禁毒,戒毒工作高質(zhì)量發(fā)展，為全面落實(shí)健康中國戰(zhàn)略做出貢獻(xiàn).

今后需要進(jìn)一步結(jié)合現(xiàn)實(shí)數(shù)據(jù)深入研究各參數(shù)與仿真結(jié)果的定量關(guān)系及影響規(guī)律,進(jìn)而將模型用于預(yù)測毒品濫用的流行趨勢(shì),為決策部門制定更有效的公共政策提供理論依據(jù).