呂建平,羅曉峰
(1.中北大學(xué) 理學(xué)院, 太原 030051; 2.山西大學(xué) 復(fù)雜系統(tǒng)研究所, 太原 030006)
馬皰疹病毒 (equine herpes virus,EHV) 屬于皰疹病毒科皰疹病毒α亞科,主要有EHV-1、EHV-2、EHV-3、EHV-4、EHV-5這5種類型。其中,引起馬患病的主要是EHV-1、EHV-4型,EHV-4型會(huì)引起馬患非致命的上呼吸道感染疾病,EHV-1型會(huì)導(dǎo)致馬患呼吸道疾病、母馬流產(chǎn)、新生馬死亡以及嚴(yán)重的神經(jīng)系統(tǒng)疾病馬腦脊髓炎(equine herpesvirus myeloencephalophthy,EHM)[1-2]。本文主要討論由EHV-1引起的馬感染。
EHV-1病毒主要通過(guò)空氣傳播,染病后一般有2~10天的潛伏期[3]。該病毒第一個(gè)病例發(fā)生在20世紀(jì)30年代的美國(guó),之后在中國(guó)新疆、馬來(lái)西亞、印度等都出現(xiàn)過(guò)類似的病例。EHV-1在大型馬場(chǎng)、馬展場(chǎng)、獸醫(yī)診所、養(yǎng)馬場(chǎng)等場(chǎng)所中曾多次爆發(fā),如2011年在美國(guó)猶他州奧格登舉行的賽馬比賽后引起的EHV-1感染。隨后,在2018年法國(guó)大型馬場(chǎng)出現(xiàn)了EHV-1感染,此次事件中,至少有200多匹馬受到影響。由于這些場(chǎng)所馬的數(shù)量比較多,疾病爆發(fā)不僅對(duì)動(dòng)物的健康造成了影響,而且?guī)?lái)了不小的經(jīng)濟(jì)損失。
本文基于EHV-1傳播的基本特點(diǎn),以2011年美國(guó)猶他州奧格登舉行的賽馬比賽后EHV-1爆發(fā)為背景,建立EHV-1傳播動(dòng)力學(xué)模型,從理論和數(shù)值上分析傳染率和恢復(fù)率對(duì)有癥狀馬基本再生數(shù)、流行最終規(guī)模以及疫苗接種策略的影響,由此提出有效的預(yù)防與控制措施。
為了更好地利用模型說(shuō)明實(shí)際問(wèn)題,對(duì)模型作如下假設(shè):
1) 由于賽后疾病爆發(fā)時(shí)間短且研究的參賽馬年齡集中且偏小,因此忽略馬的出生、死亡、遷入、遷出。
2) 假設(shè)康復(fù)且攜帶病毒的馬在加強(qiáng)飼養(yǎng)管理、合理休息的條件下不會(huì)受到刺激,這類馬記為恢復(fù)者。
3) 假設(shè)潛伏期的馬無(wú)感染性。
馬個(gè)體的狀態(tài)分為5類:易感S、潛伏E、無(wú)臨床癥狀A(yù)和有臨床癥狀I(lǐng)但都有傳染性、恢復(fù)R,在t時(shí)刻,各狀態(tài)的數(shù)量分別為S(t)、E(t)、A(t)、I(t)、R(t)。本文忽略其他傳播方式的影響,只考慮通過(guò)呼吸接觸感染EHV-1的馬。建立如下模型:
(1)
其中:β表示易感馬與感染EHV-1病毒且有傳染性的馬的有效接觸系數(shù);ε表示具有傳染性且有臨床癥狀的馬未被隔離的比例;1/ρ表示感染病毒的馬的平均潛伏期;1/γ1表示具有傳染性且無(wú)臨床癥狀馬的病程;1/γ2表示具有傳染性而有臨床癥狀的馬的病程;α表示潛伏期的馬轉(zhuǎn)化為無(wú)癥狀馬的比例。
系統(tǒng)的初始條件為:S(0)=S0,E(0)≈0,A(0)≈0,I(0)≈0,R(0)≈0。易證系統(tǒng)(1)的正向不變集為
其中N0為馬的總數(shù)且保持不變,N0∈R+。
對(duì)于系統(tǒng)(1),顯然存在無(wú)病線平衡點(diǎn)E0=(S0,0,0,0,0)=(N0,0,0,0,0),進(jìn)而根據(jù)下一代矩陣[6],得基本再生數(shù)為
以下給出R0的生物學(xué)解釋:R0表示在一個(gè)全部是易感馬的種群中,進(jìn)入一個(gè)染病馬,在其病程內(nèi)傳染的馬的數(shù)量[4]。對(duì)于EHV-1感染馬,染病馬分為無(wú)癥狀的患病馬和有癥狀的患病馬,因此R0表達(dá)式包含兩項(xiàng):βαN0/γ1和βε(1-α)N0/γ2。第1項(xiàng)表示一個(gè)無(wú)癥狀的馬在病程1/γ1內(nèi)傳染馬的總數(shù);第2項(xiàng)表示有癥狀且未被隔離的馬在病程1/γ2傳染的馬的總數(shù)。
流行病最終規(guī)模是反映疾病爆發(fā)嚴(yán)重程度的度量之一,在種群數(shù)量不變的情況下,流行病規(guī)模越大,疫情越嚴(yán)重,造成的損失越大[5]。以下給出系統(tǒng)(1)的最終規(guī)模。
當(dāng)時(shí)間t→∞時(shí),E(∞)=0,A(∞)=0,I(∞)=0。由系統(tǒng)(1)知
與
對(duì)上述兩式同時(shí)在區(qū)間[0,∞)上求積分,并通過(guò)化簡(jiǎn)可得
lnS(∞)-lnS(0)=
進(jìn)而解得S(∞),即疾病傳播過(guò)程中有幸逃避感染的易感者數(shù)。由R(∞)=S(0)-S(∞)可得感染EHV-1馬的最終規(guī)模。
研究傳染病的目的是為了預(yù)防和控制它,目前的措施包括疫苗接種、治療、檢疫、隔離等。疫苗接種是將微生物注射到個(gè)體體內(nèi),從而激活機(jī)體的免疫系統(tǒng)產(chǎn)生抗體,當(dāng)同種類型的微生物再次進(jìn)入體內(nèi)時(shí),它將被抗體消滅,進(jìn)而達(dá)到預(yù)防疾病的目的。
本文在系統(tǒng)的基礎(chǔ)上尋找最優(yōu)的疫苗接種措施。假設(shè)接種疫苗完全有效,且接種疫苗的個(gè)體歸入倉(cāng)室R(t),得到的模型如下
(2)
其中控制項(xiàng)u(t)是隨時(shí)間變化的有效接種率,狀態(tài)方程(2)仍在不變集D中考慮。由于馬的總數(shù)不變,僅考慮狀態(tài)方程(2)的前4個(gè)方程。為了使流行病最終規(guī)模和接種疫苗的花費(fèi)最小,應(yīng)用最優(yōu)控制理論[6],定義控制集
Ω={u(t)∈L1(0,T)|0≤u(t)≤Umax)},
Umax∈R+
建立如下具有控制受限的最優(yōu)控制問(wèn)題:
(3)
其中:ω1I(t)是目標(biāo)泛函中是使流行最小化的恒定成本;ω2u(t)S(t)是使接種數(shù)量最小化的恒定成本;u2(t)是疫苗接種的經(jīng)濟(jì)成本并假定其是接種率u(t)的非線性函數(shù)且是二次型。
對(duì)于任意一個(gè)u∈Ω,狀態(tài)方程(2)有唯一的解X=(S,E,A,I)。如果存在最優(yōu)控制u*∈Ω,相對(duì)應(yīng)的狀態(tài)變量X*=(S*,E*,A*,I*),使得目標(biāo)泛函最小化,則u*就是控制問(wèn)題(2)和(3)的一個(gè)解。下面使用Filippov-Cesari存在性定理[7]證明最優(yōu)控制對(duì)(X*,u*)的存在性。
命題1 最優(yōu)控制問(wèn)題(2)和(3)有解。
證明按照Filippov-Cesari存在性定理定義一個(gè)集合
N(t,X)={g(X,u)+ξ,f(X,u)|ξ≤0,u∈Ω}
其中:X=(S,E,A,I)表示狀態(tài)變量;g(X,u)表示目標(biāo)泛函Ψ(u)的被積函數(shù)。對(duì)任意(t,X)∈R5,f(X,u)定義為
f(X,u)=(f1(X,u),f2(X,u),f3(X,u),f4(X,u))T
其中f1(X,u)、f2(X,u)、f3(X,u)和f4(X,u)分別為狀態(tài)方程(2)前4個(gè)方程的右端項(xiàng)。
下證N(t,X)對(duì)每一個(gè)(t,X)都是凸的,即只需證對(duì)于任意y1,y2∈N(t,X),滿足對(duì)?η∈[0,1]都有
ηy1+(1-η)y2∈N(t,X)
事實(shí)上y1,y2∈N(t,X)就說(shuō)明存在ξ1,ξ2≤0及u1(t),u2(t)∈Ω,使得
yi={g(X,ui)+ξi,f(X,ui)},i=1,2
進(jìn)而,對(duì)yi的第1個(gè)分量,
令
u3(t)=ηu1(t)+(1-η)u2(t)
則u3(t)、u4(t)∈Ω。顯然,ηξ1+(1-η)ξ2=ξ3≤0,因此,y1、y2的第1個(gè)分量屬于N(t,X)。
下證yi的第2個(gè)分量屬于N(t,X)。
因此y1、y2的第2個(gè)分量屬于N(t,X)。
由以上證明知yi的2個(gè)分量都屬于N(t,X),因此N(t,X)是凸的。顯然Ω是一個(gè)緊集,且0≤S≤N0,0≤E≤N0,0≤A≤N0,0≤I≤N0有界,根據(jù)Filippov-Cesari存在性定理可證命題1成立。
最優(yōu)控制問(wèn)題(2)和(3)的解滿足的最優(yōu)性條件可通過(guò)Pontryagin’s極值原理得到[8]。令λ(t)=(λS(t),λE(t),λA(t),λI(t))為伴隨變量,則控制問(wèn)題的哈密頓泛函定義為
H(X(t),λ(t),u(t))=ω1I(t)+ω2u(t)S(t)+
u2(t)+λS(-βSA-βεSI-u(t)S)+
λE(βSA+βεSI-ρE)+λA(αρE-γ1A)+
λI((1-α)ρE-γ2I)
從Pontryagin極值原理的必要條件可以求得u*的表達(dá)式:
(4)
和伴隨系統(tǒng)方程:
(5)
狀態(tài)方程(2)、伴隨方程(5)和最優(yōu)性條件(4)構(gòu)成了最優(yōu)控制問(wèn)題(2)和(3)的最優(yōu)系統(tǒng)。接下來(lái)通過(guò)數(shù)值模擬來(lái)研究控制問(wèn)題。
利用美國(guó)農(nóng)業(yè)部發(fā)布的2011.5.8—5.27猶他州奧格登舉行賽馬比賽新增病例數(shù)據(jù),結(jié)合模型(1)進(jìn)行數(shù)據(jù)擬合,估計(jì)參數(shù),提出相應(yīng)的防控措施,并預(yù)測(cè)賽事中有癥狀染病馬和無(wú)癥狀染病馬的變化趨勢(shì)。
基于模型(1)結(jié)合新增病例數(shù)據(jù)使用最小二乘估計(jì)β的值得β=0.017,根據(jù)美國(guó)農(nóng)業(yè)部官網(wǎng)數(shù)據(jù),ρ∈[1/10,1/2],γ1∈[1/14,1/7],γ2∈[1/60,1/7]。此外,假定ε=0.12,α=0.5。將所得參數(shù)值代入R0可知當(dāng)N0<10時(shí)R0<1,這說(shuō)明如果賽前在1個(gè)馬廄中放少于10匹馬,并減少它們的接觸,可有效控制疾病傳播以及爆發(fā)。
模型與數(shù)據(jù)的擬合見(jiàn)圖1,可見(jiàn)擬合結(jié)果很好,驗(yàn)證了EHV-1病毒傳播模型(1)的合理性。圖2給出了無(wú)癥狀馬和有癥狀馬隨時(shí)間的變化趨勢(shì)??梢钥闯鲇邪Y狀的染病馬比無(wú)癥狀馬的峰值大且流行時(shí)間長(zhǎng),因此對(duì)有癥狀馬的觀測(cè)、隔離、治療尤為重要。
圖1 模型與數(shù)據(jù)的擬合
圖2 無(wú)癥狀染病馬隨時(shí)間的數(shù)量變化
圖3是感染EHV-1病毒的有癥狀馬對(duì)傳染率β、恢復(fù)率γ2的敏感性分析。圖3(a)中的β小時(shí),峰值低且到達(dá)峰值的時(shí)間晚;但當(dāng)β增加到一定程度時(shí)(圖3(a)中β>0.1),峰值和到達(dá)峰值時(shí)間相差不大。因此,減少易感的馬與有癥狀的染病馬的接觸機(jī)會(huì)可以在一定程度上控制疾病。圖3(b)說(shuō)明了恢復(fù)率γ2越小,染病馬的數(shù)量越多。在醫(yī)療資源充足的情況下,可通過(guò)加強(qiáng)對(duì)有癥狀馬的治療和減少患病病程來(lái)控制疾病。
圖4刻畫了參數(shù)β、γ1(圖4(a))和α、γ2(圖4(b))對(duì)基本再生數(shù)R0的影響,R0是β和α的增函數(shù), 是γ1和γ2的減函數(shù)。
圖3 β和γ2對(duì)I(t)的影響
圖4 β、γ1和α、γ2對(duì)R0的影響
圖5給出了在最優(yōu)控制措施下接種率u*的變化情況,表明初始階段使用最大的疫苗接種率,然后隨時(shí)間減少,有癥狀的染病馬對(duì)易感馬的影響得到明顯控制且成本最少。圖6表示接種和不接種I(t)隨時(shí)間的變化情況。采取控制措施,I(t)增長(zhǎng)相對(duì)緩慢,因此疫苗接種后,有癥狀的染病馬數(shù)量減少,易感馬受到了有效的保護(hù)。
圖5 最優(yōu)控制措施下接種率u*的變化
圖6 接種、不接種I(t)隨時(shí)間的變化
馬皰疹病毒EHV-1會(huì)導(dǎo)致馬患有馬腦脊髓炎,尤其對(duì)賽馬會(huì)造成了很大的經(jīng)濟(jì)損失。目前治療方法少,多采用支持療法,效果也不顯著,病情嚴(yán)重者只能采用安樂(lè)死。基于此,本文建立EHV-1疾病傳播模型,計(jì)算疾病的基本再生數(shù)和流行最終規(guī)模,并通過(guò)最優(yōu)控制理論給出最優(yōu)的疫苗接種策略,既節(jié)省了成本,又能有效控制疾病的爆發(fā)。利用2011年5月8日—27日猶他州數(shù)據(jù)進(jìn)行參數(shù)擬合,估計(jì)模型的參數(shù),利用獲得的參數(shù)值進(jìn)行數(shù)值模擬,表明減少易感馬與有癥狀染病馬的接觸以及加強(qiáng)對(duì)有癥狀染病馬的治療可以有效控制疾病。