馬感染EHV-1病毒傳播動(dòng)力學(xué)模型分析

呂建平，羅曉峰

(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病毒主要通過空氣傳播，染病后一般有2～10天的潛伏期[3]。該病毒第一個(gè)病例發(fā)生在20世紀(jì)30年代的美國，之后在中國新疆、馬來西亞、印度等都出現(xiàn)過類似的病例。EHV-1在大型馬場(chǎng)、馬展場(chǎng)、獸醫(yī)診所、養(yǎng)馬場(chǎng)等場(chǎng)所中曾多次爆發(fā)，如2011年在美國猶他州奧格登舉行的賽馬比賽后引起的EHV-1感染。隨后，在2018年法國大型馬場(chǎng)出現(xiàn)了EHV-1感染，此次事件中，至少有200多匹馬受到影響。由于這些場(chǎng)所馬的數(shù)量比較多，疾病爆發(fā)不僅對(duì)動(dòng)物的健康造成了影響，而且?guī)砹瞬恍〉慕?jīng)濟(jì)損失。

本文基于EHV-1傳播的基本特點(diǎn)，以2011年美國猶他州奧格登舉行的賽馬比賽后EHV-1爆發(fā)為背景，建立EHV-1傳播動(dòng)力學(xué)模型，從理論和數(shù)值上分析傳染率和恢復(fù)率對(duì)有癥狀馬基本再生數(shù)、流行最終規(guī)模以及疫苗接種策略的影響，由此提出有效的預(yù)防與控制措施。

1 模型的建立

為了更好地利用模型說明實(shí)際問題，對(duì)模型作如下假設(shè)：

1) 由于賽后疾病爆發(fā)時(shí)間短且研究的參賽馬年齡集中且偏小，因此忽略馬的出生、死亡、遷入、遷出。

2) 假設(shè)康復(fù)且攜帶病毒的馬在加強(qiáng)飼養(yǎng)管理、合理休息的條件下不會(huì)受到刺激，這類馬記為恢復(fù)者。

3) 假設(shè)潛伏期的馬無感染性。

馬個(gè)體的狀態(tài)分為5類：易感S、潛伏E、無臨床癥狀A(yù)和有臨床癥狀I(lǐng)但都有傳染性、恢復(fù)R，在t時(shí)刻，各狀態(tài)的數(shù)量分別為S(t)、E(t)、A(t)、I(t)、R(t)。本文忽略其他傳播方式的影響，只考慮通過呼吸接觸感染EHV-1的馬。建立如下模型：

(1)

其中：β表示易感馬與感染EHV-1病毒且有傳染性的馬的有效接觸系數(shù)；ε表示具有傳染性且有臨床癥狀的馬未被隔離的比例；1/ρ表示感染病毒的馬的平均潛伏期；1/γ1表示具有傳染性且無臨床癥狀馬的病程；1/γ2表示具有傳染性而有臨床癥狀的馬的病程；α表示潛伏期的馬轉(zhuǎn)化為無癥狀馬的比例。

系統(tǒng)的初始條件為：S(0)=S0，E(0)≈0，A(0)≈0，I(0)≈0，R(0)≈0。易證系統(tǒng)(1)的正向不變集為

其中N0為馬的總數(shù)且保持不變，N0∈R+。

1.1 平衡點(diǎn)和基本再生數(shù)

對(duì)于系統(tǒng)(1)，顯然存在無病線平衡點(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感染馬，染病馬分為無癥狀的患病馬和有癥狀的患病馬，因此R0表達(dá)式包含兩項(xiàng)：βαN0/γ1和βε(1-α)N0/γ2。第1項(xiàng)表示一個(gè)無癥狀的馬在病程1/γ1內(nèi)傳染馬的總數(shù)；第2項(xiàng)表示有癥狀且未被隔離的馬在病程1/γ2傳染的馬的總數(shù)。

1.2 流行病最終規(guī)模

流行病最終規(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,∞)上求積分，并通過化簡(jiǎn)可得

lnS(∞)-lnS(0)=

進(jìn)而解得S(∞)，即疾病傳播過程中有幸逃避感染的易感者數(shù)。由R(∞)=S(0)-S(∞)可得感染EHV-1馬的最終規(guī)模。

2 控制措施

研究傳染病的目的是為了預(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è)體歸入倉室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)控制問題：

(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*就是控制問題(2)和(3)的一個(gè)解。下面使用Filippov-Cesari存在性定理[7]證明最優(yōu)控制對(duì)(X*,u*)的存在性。

命題1 最優(yōu)控制問題(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)就說明存在ξ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)控制問題(2)和(3)的解滿足的最優(yōu)性條件可通過Pontryagin’s極值原理得到[8]。令λ(t)=(λS(t),λE(t),λA(t),λI(t))為伴隨變量，則控制問題的哈密頓泛函定義為

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)控制問題(2)和(3)的最優(yōu)系統(tǒng)。接下來通過數(shù)值模擬來研究控制問題。

3 數(shù)據(jù)擬合以及敏感性分析

利用美國農(nóng)業(yè)部發(fā)布的2011.5.8—5.27猶他州奧格登舉行賽馬比賽新增病例數(shù)據(jù)，結(jié)合模型(1)進(jìn)行數(shù)據(jù)擬合，估計(jì)參數(shù)，提出相應(yīng)的防控措施，并預(yù)測(cè)賽事中有癥狀染病馬和無癥狀染病馬的變化趨勢(shì)。

基于模型(1)結(jié)合新增病例數(shù)據(jù)使用最小二乘估計(jì)β的值得β=0.017，根據(jù)美國農(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，這說明如果賽前在1個(gè)馬廄中放少于10匹馬，并減少它們的接觸，可有效控制疾病傳播以及爆發(fā)。

模型與數(shù)據(jù)的擬合見圖1，可見擬合結(jié)果很好，驗(yàn)證了EHV-1病毒傳播模型(1)的合理性。圖2給出了無癥狀馬和有癥狀馬隨時(shí)間的變化趨勢(shì)?？梢钥闯鲇邪Y狀的染病馬比無癥狀馬的峰值大且流行時(shí)間長(zhǎng)，因此對(duì)有癥狀馬的觀測(cè)、隔離、治療尤為重要。

圖1 模型與數(shù)據(jù)的擬合

圖2 無癥狀染病馬隨時(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)說明了恢復(fù)率γ2越小，染病馬的數(shù)量越多。在醫(yī)療資源充足的情況下，可通過加強(qiáng)對(duì)有癥狀馬的治療和減少患病病程來控制疾病。

圖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í)間的變化

4 結(jié)束語

馬皰疹病毒EHV-1會(huì)導(dǎo)致馬患有馬腦脊髓炎，尤其對(duì)賽馬會(huì)造成了很大的經(jīng)濟(jì)損失。目前治療方法少，多采用支持療法，效果也不顯著，病情嚴(yán)重者只能采用安樂死?；诖?，本文建立EHV-1疾病傳播模型，計(jì)算疾病的基本再生數(shù)和流行最終規(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ì)有癥狀染病馬的治療可以有效控制疾病。