一個具有飽和增殖率的抗體免疫HIV模型的動力學與最優控制

陳沖 周英告 周凱

(中南大學數學與統計學院,長沙，410083）

1 引言

研制能夠預防艾滋病毒的疫苗仍然是全球公共衛生的一個重要優先事項.一種理想的疫苗可能需要產生能夠中和廣泛傳播的初級HIV-1分離株的抗體.這樣的抗體被稱為廣譜或廣泛中和抗體(broadly neutralizing antibody,bnAb)).廣泛中和抗體是治療和預防HIV-1感染的有希望的候選藥物.這種抗體可以暫時抑制感染個體的病毒載量.如果在接觸HIV病毒之前,一個個體體內就有抗體的話,則這個個體就能夠極大地降低被HIV感染的風險[27].迄今為止，在人類中評估的方法都沒有產生廣泛中和抗體，而開發一種能夠實現這一目標的疫苗策略的努力是目前艾滋病疫苗設計的前沿.bnAb的發現和鑒定一直是HIV疫苗科學中最具啟發性和鼓舞人心的方面，為以往疫苗失敗提供了見解，并提出了前進的方向[18].

在2019年,來自美國斯克里普斯研究所和非營利性疫苗研究組織國際艾滋病疫苗計劃（IAVI）的研究人員報道他們開發出一種實驗性HIV疫苗,它可以觸發中和多種HIV毒株的抗體–廣泛中和抗體(bnAb).他們對12只兔子進行疫苗接種,并將獲得的結果與僅接種攜帶這種聚糖屏蔽層的Env版本的對照組進行了比較.他們發現他們的疫苗策略具有更好的反應,其中的5只兔子產生能夠中和多種HIV分離株的抗體.這標志著HIV疫苗已達到了一個重要的里程碑[5].

在2021年,來自美國賓夕法尼亞大學佩雷爾曼醫學院和杜克大學等研究機構的研究人員從3個產生了bnAb的HIV-1感染者身上構建了攜帶主要傳播Env的猿-人類免疫缺陷病毒（simian-human immunodeficiency virus,SHIV），并使用這些SHIV感染了22只恒河猴.研究發現有7只恒河猴產生的bNAb表現出相當大的中和廣度和效力.這為開發高效的HIV疫苗奠定了基礎[17,19].

在以往的研究中,由于抗AIDS的疫苗一直未進入臨床實踐,因而對抗體免疫的動力學研究僅停留于抗體因病毒刺激而自然誘發的階段.但隨著疫苗研究好消息的不斷傳來,討論HIV/AIDS的接種問題不僅必要而且具有很重要的理論和實際意義.

為了研究抗體在病毒刺激下誘發免疫的動力學行為,一些數學模型被構造出來以揭示抗體免疫在不同情況下的免疫行為.在這些模型中,通常用T,L,I,V,E,和B分別表示健康的T細胞,潛在被感染的T細胞,被感染T細胞,病毒,CTL免疫細胞(效用細胞)和體液免疫B細胞(記憶細胞).

在2007年,文獻[12]僅在考慮病毒和免疫響應的情況下,構建了如下簡化的免疫病毒模型

這個模型假定病毒具有常數復制率sV,并被效用T細胞和記憶T細胞所清除.模型假定受感染細胞的數量與病毒的數量成正比,效用細胞探測和破壞受感染細胞的概率PE與受感染細胞的數量成正比,推出PE=pV,其中p是一個常數.根據文獻[2]的證據,記憶細胞積極參與免疫反應,因此假設記憶細胞也殺死感染細胞.類似于效用細胞的作用,記憶細胞探測和破壞受感染細胞的概率為qV,其中q是一個常數.效用細胞和記憶細胞受病毒刺激由一個恒定的天然細胞庫N而產生,其速率分別為sEV和sBV.效用細胞按常數率μE衰減,而記憶細胞維持內穩態受邏輯規律機制的影響.

在2005年,文獻[14]在文獻[15]構造的一個基本TIV模型的基礎上,僅考慮體液免疫的情況下,構建了下列的病毒免疫模型

這個模型假定未感染T細胞以恒定的速率sT從體內得到補充,其自然預期壽命為1/μT天.細胞因接觸病原體(病毒)而感染,并以k1V的速度轉化為受感染細胞.感染細胞以常數率μI自然死亡.細胞的死亡導致每一個被感染的細胞釋放N個病原體,這些病原體在血液中有一個預期壽命1/μV.病原體要么死亡，要么成功感染新細胞.病原體被免疫系統以k2B的速度清除.病原體特異性淋巴細胞與病原體接觸后以sBV的速度增殖,并以μB速率自然死亡.這個模型隨后又被文獻[22]進一步擴充成一個時滯系統.

在2010年,一個類似于模型(1.2)的TIVB模型被文獻[10]所考慮.二者的差異是文獻[10]在考慮免疫細胞(B細胞)的增殖時,認為B細胞的增殖速率就是sBV,與B細胞無關.從而,將上述模型(1.2)中的第四個方程改為如下方程

模型(1.2)后來被文獻[6]和文獻[23]擴展成具有一般感染率的如下模型：

其中,新的參數sI表示受感染細胞產生游離病毒的速率常數,ψ(T,V)表示感染的發生率,ψ是一個滿足一定條件的二元函數.

在2015年,在考慮了潛在感染的前提下,文獻[7]將上述TIVB模型擴充為如下具有Beddington-DeAngelis功能響應的TLIVB模型：

其中,新的參數a和b都是正常數,1?α和α分別表示一個未感染的細胞在感染后變成潛伏感染和感染細胞的概率(0<α<1),潛在被感染細胞L以常數率μL自然死亡,并以常數率δ轉化成感染細胞.

為了同時研究CTL免疫與抗體免疫的影響,文獻[21]和[1]基于一些已有的模型構建了如下具有Beddington-DeAngelis功能響應的TIVEB模型

與此同時,文獻[9]構建了具有一般形式的TIVEB免疫病毒模型

其中,新的函數f(T,I,V)描述了易感T細胞被感染的發生率,是一個滿足一定條件的功能函數,包含了感染率的各種形式.

綜觀以上模型發現,所有的模型在構建抗體免疫動力學時,都假定B細胞以線性性自然增殖,但這似乎是不太合理的.因為我們知道,由于個體的差異,每個人的免疫能力是不一樣的,這意味著在受到病毒刺激的情況下,病毒誘發的B細胞增殖不可能是線性增長的.

基于以上考慮,本文將以飽和增長的形式構建B細胞的動力學.同時,考慮到疫苗的誘導作用和藥物的治療作用,將建立如下的HIV病毒控制模型：

其中,新的參數u1,u2分別表示疫苗的誘導力和治療艾滋病藥物PI對病毒的抑制力,0≤ui≤1,i=1,2.接種疫苗不僅會人為激活抗體的形成,而且會由此降低健康T細胞被感染的風險,因此,疫苗的作用會同時影響T細胞,I細胞及B細胞動力學.ui=1意味著接種或治療完全有效,ui=0意味著沒有接種或沒有施加治療(i=1,2).

當不施加疫苗和藥物治療時,上述控制系統簡化為如下動力系統：

以上模型中的所有參數都是正數,且a<1,初始狀態T(0),I(0),V(0)和B(0)皆為非負數.此外,基于生物學意義,我們總假定

2 動力系統的穩定性分析

進一步,在R0>1時,可得該系統唯一的無免疫地方病平衡點其中,

在HB>0且R1>1時,可得系統(1.9)唯一的免疫地方病平衡點其中,

我們知道,R0越小,病毒越不容易傳播.從以上R0的結果看,若能降低感染率,則可以降低病毒傳播的風險.

2.1 無病平衡點的穩定性

定理1當R0<1時，無病平衡點是局部漸近穩定的；當R0>1時，無病平衡點是不穩定的.

證明系統(1.9)在無病平衡點處的Jacobi矩陣為

對應的特征方程為

若記其特征根分別為λ1,λ2,λ3,λ4,則顯然有

可以看出,λ3λ4>0?R0<1,這意味著當R0<1時,所有的特征根都是負實數,故系統(1.9)的無病平衡點在R0<1時是局部漸近穩定的,而在R0>1時是不穩定的.

定理2當R0<1時,系統(1.9)的無病平衡點是全局漸近穩定的.

證明定義一個Lyapunov函數

因此,當R0<1時,這意味著系統(1.9)的無病平衡點是全局漸近穩定的.

2.2 無免疫地方病平衡點的穩定性

定理3當R0>1且HB>0時,系統(1.9)的無免疫地方病平衡點在R1<1時是局部漸近穩定的,在R1>1時是不穩定的.

證明系統(1.9)在處的Jacobi矩陣為：

它的特征多項式為

其中,

顯然,ai>0,i=1,2,3.

由特征多項式可以得到其中一個特征根為

易知,當HB>0時,

其余特征根滿足如下方程：

又

因此,由Routh–Hurwitz定理知結論成立.

定理4當R0>1且HB>0時,系統(1.9)的無免疫地方病平衡點且

時是全局漸近穩定的.

證明 定義一個Lyapunov函數

則

注意到

于是,

注意到

故,當

2.3 免疫地方病平衡點的穩定性

定理5當HB>0且R1>1時,系統(1.9)的免疫地方病平衡點是局部漸近穩定的.

證明注意到,當HB>0且R1>1時,免疫地方病平衡點是唯一的,且系統(1.9)在免疫地方病平衡點處的Jacobi矩陣為

計算可得

于是,

其特征多項式為

其中,

注意到

這導致

于是,

定理6當HB>1且R1>1時,系統(1.9)的免疫地方病平衡點是全局漸近穩定的.

證明注意到,當HB>0且R1>1時,免疫地方病平衡點是唯一的,且在免疫作用下,由此可知,當且僅當時,B取得最小值即,在免疫作用下,恒有

定義一個Lyapunov函數

則有

因為

所以

3 最優治療策略

為研究接種和藥物治療AIDS的有效性,我們將極小化感染細胞和病毒的總量,并使接種和藥物治療的代價(成本)最小.控制變量u1代表接種的有效率,u2代表藥物對HIV病毒的治療效果,ui∈[0,1],i=1,2.u1(t)=1表示完全有效的接種,u1(t)=0表示完全無效的接種,或者說無接種?u2(t)=1表示完全有效的治療,u1(t)=0表示完全無效的治療,或者說無治療.

最優控制問題描述為

并使之滿足控制系統(1.8),其中,正常數αi刻畫控制變量ui的代價(成本),i=1,2,u=(u1,u2).我們將尋找最優控制對使得

3.1 最優控制的存在性

定理7對最優控制問題(1.8)+(3.1),一定存在一個最優控制對使得J(u*)=

證明為證明存在性,按文[8]的定理III.4.1,需驗證下列5個條件:(1)容許控制和相應的狀態變量的集合非空?(2)容許控制集U是凸閉集?(3)狀態系統的右邊關于狀態變量和控制變量都是有界的?(4)目標函數中的被積函數關于容許控制U是凸的?(5)存在常數c1,c2>0和b>1,使得目標函數中的被積函數有上界

首先,由前面關于系統(1.9)的動力學討論知條件(1)和(3)成立.其次,由定義知容許控制集非空,且目標函數中的被積函數在U上是凸函數.最后,不難知道

其中,c1依賴于I和V的上界,c2可取α1與α2中的最大值的2倍,b可以取任一大于2的實數.因此,最優控制問題(1.8)+(3.1)一定有解.

為便于得到最優控制問題(1.8)+(3.1)的最優性條件,記Hamilton函數H(T,I,V,B;u;λ1,λ2,λ3,λ4)為

其中,λi(i=1,2,···,4)是協態變量.由Pontryagin極大值原理,可得如下定理8.

定理8若是最優控制問題(1.8)+(3.1)的最優控制對,T*,I*,V*,B*是相應的最優軌道,則存在協態變量λi(i=1,2,···,4),使得,除控制系統(1.8)成立外,還滿足下列條件：

(i)協態方程：

(ii)最優性條件:

這意味著

(iii)橫截條件：

最優性系統由具有初始狀態的控制系統(1.8),結合帶有橫截條件的協態方程(3.2),并跟最優控制的表達式(3.3)共同組成.因此,最優控制問題(1.8)+(3.1)的最優性系統是

其中,

注1利用文獻[27]提供的方法,類似可以證明：對充分小的tf,最優性系統(3.4)存在唯一的解.

3.2 數值仿真

由于最優性系統是非線性的,因而很難獲得最優控制問題的解析解.本文將用梯度投影法對最優控制問題(1.8)+(3.1)進行數值模擬.數值模擬所采用的軟件是MATLAB.對控制系統(1.8)中普遍存在的參數將參照文獻[3,4,16,24,27]中的取值(范圍)進行確定,新的參數sB和a將利用文獻[20]中的數據,并運用文獻[25]提供的方法結合最小二乘法進行擬合獲得.所有參數的取值(范圍)如表1所示.

表1 參數意義及其取值表

在具體的數值模擬中,將實際選擇參數如下：μT=0.02,μI=0.25,μV=2.1,μB=0.025,k1=0.00026,k2=0.023,sT=6,sB=1,a=0.5,N=150,并取初始狀態為T(0)=200,I(0)=1,V(0)=2,B(0)=0.5,治療周期為90天.同時,假定接種與施藥的成本并不高.事實上,模擬結果顯示,在高成本下(提高成本1000倍),雙控的效果與單純的藥物治療,即u2控制相當.因此,為比較方便起見,我們只給出在較低成本下雙控制與單控制的模擬結果.

圖1 描述的是在疫苗誘導和藥物治療雙控制作用下的最優控制結果.圖2描述的是在疫苗誘導單控制作用下的最優控制結果.圖3描述的是在藥物治療單控制作用下的最優控制結果.由于我們最感興趣的是抗體免疫與藥物治療對感染細胞和病毒的影響,因此,我們僅關注最優控制下感染細胞和病毒的變化.從以上各圖可以看出,在疫苗有效的情況下,單獨的疫苗誘導控制的作用效果幾乎與雙控制下的效果一樣,說明疫苗誘導下抗體免疫具有更好的抗HIV效果,即使不配合其它抗病毒藥物治療也會極大地增加抵御AIDS風險,甚至最終戰勝AIDS.從圖2,或圖1可以看出,在抗體免疫作用下,感染細胞I和病毒V的濃度很快在抗體的作用下迅速降低,并保持在較低水平,直至消亡.但從圖3卻發現一個比較奇特的現象,即在單藥物治療的情況下,雖然在治療的后期,感染細胞和病毒的濃度最終被控制在一個較低的水平,但在用藥初期,感染細胞和病毒的濃度有明顯的繼續上升的趨勢,甚至病毒的濃度一度比不施加任何控制的情況下的濃度還要高,這是否是因為藥物的刺激作用導致了病毒短時間強力增強呢？這個問題有待得到臨床驗證.

圖1 雙控制下健康細胞、感染細胞、病毒、抗體細胞及疫苗和藥物有效性動力學性態

圖2 疫苗控制下健康細胞、感染細胞、病毒、抗體細胞及疫苗有效性動力學性態

圖3 藥物控制下健康細胞、感染細胞、病毒、抗體細胞及藥物有效性動力學性態