呂建平，羅曉峰

(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型會引起馬患非致命的上呼吸道感染疾病，EHV-1型會導(dǎo)致馬患呼吸道疾病、母馬流產(chǎn)、新生馬死亡以及嚴重的神經(jīng)系統(tǒng)疾病馬腦脊髓炎(equine herpesvirus myeloencephalophthy，EHM)[1-2]。本文主要討論由EHV-1引起的馬感染。

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

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

1 模型的建立

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

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

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

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

馬個體的狀態(tài)分為5類：易感S、潛伏E、無臨床癥狀A(yù)和有臨床癥狀I(lǐng)但都有傳染性、恢復(fù)R，在t時刻，各狀態(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 平衡點和基本再生數(shù)

對于系統(tǒng)(1)，顯然存在無病線平衡點E0=(S0,0,0,0,0)=(N0,0,0,0,0)，進而根據(jù)下一代矩陣[6]，得基本再生數(shù)為

以下給出R0的生物學(xué)解釋：R0表示在一個全部是易感馬的種群中，進入一個染病馬，在其病程內(nèi)傳染的馬的數(shù)量[4]。對于EHV-1感染馬，染病馬分為無癥狀的患病馬和有癥狀的患病馬，因此R0表達式包含兩項：βαN0/γ1和βε(1-α)N0/γ2。第1項表示一個無癥狀的馬在病程1/γ1內(nèi)傳染馬的總數(shù)；第2項表示有癥狀且未被隔離的馬在病程1/γ2傳染的馬的總數(shù)。

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

流行病最終規(guī)模是反映疾病爆發(fā)嚴重程度的度量之一，在種群數(shù)量不變的情況下，流行病規(guī)模越大，疫情越嚴重，造成的損失越大[5]。以下給出系統(tǒng)(1)的最終規(guī)模。

當(dāng)時間t→∞時，E(∞)=0，A(∞)=0，I(∞)=0。由系統(tǒng)(1)知

與

對上述兩式同時在區(qū)間[0,∞)上求積分，并通過化簡可得

lnS(∞)-lnS(0)=

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

2 控制措施

研究傳染病的目的是為了預(yù)防和控制它，目前的措施包括疫苗接種、治療、檢疫、隔離等。疫苗接種是將微生物注射到個體體內(nèi)，從而激活機體的免疫系統(tǒng)產(chǎn)生抗體，當(dāng)同種類型的微生物再次進入體內(nèi)時，它將被抗體消滅，進而達到預(yù)防疾病的目的。

本文在系統(tǒng)的基礎(chǔ)上尋找最優(yōu)的疫苗接種措施。假設(shè)接種疫苗完全有效，且接種疫苗的個體歸入倉室R(t)，得到的模型如下

(2)

其中控制項u(t)是隨時間變化的有效接種率，狀態(tài)方程(2)仍在不變集D中考慮。由于馬的總數(shù)不變，僅考慮狀態(tài)方程(2)的前4個方程。為了使流行病最終規(guī)模和接種疫苗的花費最小，應(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)濟成本并假定其是接種率u(t)的非線性函數(shù)且是二次型。

對于任意一個u∈Ω，狀態(tài)方程(2)有唯一的解X=(S,E,A,I)。如果存在最優(yōu)控制u*∈Ω，相對應(yīng)的狀態(tài)變量X*=(S*,E*,A*,I*)，使得目標(biāo)泛函最小化，則u*就是控制問題(2)和(3)的一個解。下面使用Filippov-Cesari存在性定理[7]證明最優(yōu)控制對(X*,u*)的存在性。

命題1 最優(yōu)控制問題(2)和(3)有解。

證明按照Filippov-Cesari存在性定理定義一個集合

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ù)。對任意(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個方程的右端項。

下證N(t,X)對每一個(t,X)都是凸的，即只需證對于任意y1,y2∈N(t,X)，滿足對?η∈[0,1]都有

ηy1+(1-η)y2∈N(t,X)

事實上y1,y2∈N(t,X)就說明存在ξ1，ξ2≤0及u1(t),u2(t)∈Ω，使得

yi={g(X,ui)+ξi,f(X,ui)},i=1,2

進而，對yi的第1個分量，

令

u3(t)=ηu1(t)+(1-η)u2(t)

則u3(t)、u4(t)∈Ω。顯然，ηξ1+(1-η)ξ2=ξ3≤0，因此，y1、y2的第1個分量屬于N(t,X)。

下證yi的第2個分量屬于N(t,X)。

因此y1、y2的第2個分量屬于N(t,X)。

由以上證明知yi的2個分量都屬于N(t,X)，因此N(t,X)是凸的。顯然Ω是一個緊集，且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*的表達式：

(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)進行數(shù)據(jù)擬合，估計參數(shù)，提出相應(yīng)的防控措施，并預(yù)測賽事中有癥狀染病馬和無癥狀染病馬的變化趨勢。

基于模型(1)結(jié)合新增病例數(shù)據(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時R0<1，這說明如果賽前在1個馬廄中放少于10匹馬，并減少它們的接觸，可有效控制疾病傳播以及爆發(fā)。

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

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

圖2 無癥狀染病馬隨時間的數(shù)量變化

圖3是感染EHV-1病毒的有癥狀馬對傳染率β、恢復(fù)率γ2的敏感性分析。圖3(a)中的β小時，峰值低且到達峰值的時間晚；但當(dāng)β增加到一定程度時(圖3(a)中β>0.1)，峰值和到達峰值時間相差不大。因此，減少易感的馬與有癥狀的染病馬的接觸機會可以在一定程度上控制疾病。圖3(b)說明了恢復(fù)率γ2越小，染病馬的數(shù)量越多。在醫(yī)療資源充足的情況下，可通過加強對有癥狀馬的治療和減少患病病程來控制疾病。

圖4刻畫了參數(shù)β、γ1(圖4(a))和α、γ2(圖4(b))對基本再生數(shù)R0的影響，R0是β和α的增函數(shù)， 是γ1和γ2的減函數(shù)。

圖3 β和γ2對I(t)的影響

圖4 β、γ1和α、γ2對R0的影響

圖5給出了在最優(yōu)控制措施下接種率u*的變化情況，表明初始階段使用最大的疫苗接種率，然后隨時間減少，有癥狀的染病馬對易感馬的影響得到明顯控制且成本最少。圖6表示接種和不接種I(t)隨時間的變化情況。采取控制措施，I(t)增長相對緩慢，因此疫苗接種后，有癥狀的染病馬數(shù)量減少，易感馬受到了有效的保護。

圖5 最優(yōu)控制措施下接種率u*的變化

圖6 接種、不接種I(t)隨時間的變化

4 結(jié)束語

馬皰疹病毒EHV-1會導(dǎo)致馬患有馬腦脊髓炎，尤其對賽馬會造成了很大的經(jīng)濟損失。目前治療方法少，多采用支持療法，效果也不顯著，病情嚴重者只能采用安樂死?；诖?，本文建立EHV-1疾病傳播模型，計算疾病的基本再生數(shù)和流行最終規(guī)模，并通過最優(yōu)控制理論給出最優(yōu)的疫苗接種策略，既節(jié)省了成本，又能有效控制疾病的爆發(fā)。利用2011年5月8日—27日猶他州數(shù)據(jù)進行參數(shù)擬合，估計模型的參數(shù)，利用獲得的參數(shù)值進行數(shù)值模擬，表明減少易感馬與有癥狀染病馬的接觸以及加強對有癥狀染病馬的治療可以有效控制疾病。