COVID-19急癥與重癥宿主內動力學建模研究

沈俊強,唐旭清

(江南大學理學院,中國江蘇 無錫 214122)

截至2022年年初,新型冠狀病毒肺炎(corona virus disease 2019,COVID-19)累計確診人數已達3億、累計死亡人數已達540萬,給世界各國人民的生命安全帶來了嚴重的威脅[1]。嚴重急性呼吸綜合征冠狀病毒2(severe acute respiratory syndrome coronavirus 2,SARS-CoV-2)是 COVID-19 的致病病毒,主要通過呼吸道飛沫和接觸的方式進行傳播[2]。在被SARS-CoV-2病毒感染的人群中,20%的患者表現出急癥或是重癥,并伴有急性呼吸窘迫綜合征和休克等癥狀,其中部分患者在被治愈后依舊會留有后遺癥[3～5]。對COVID-19的宿主內動力學進行建模研究,有助于定量分析COVID-19患者體內的相關變化與差異,從而為COVID-19的治療提供一定的參考與幫助。

當前,COVID-19的建模工作大部分是從宏觀角度出發,根據確診人數等相關數據進行模型擬合,以研究COVID-19的傳播動力學、預測疫情走勢以及評估防控措施的影響。He等[6]提出了一個包含隔離與外部輸入的傳染病模型,以模擬湖北省COVID-19疫情的演變;Sarkar等[7]使用傳染病模型來預測印度COVID-19疫情的拐點和持續時間;Moore等[8]使用傳染病模型和英國疫情數據來研究疫苗接種的影響以及COVID-19疫情的長期動態。這些宏觀水平的研究為COVID-19疫情的防控提供了重要的參考與理論依據,也為降低疫情的傳播風險和感染規模做出了巨大的貢獻。但是,這些宏觀層面的研究無法描述COVID-19患者體內的相關變化,也無法解釋不同患者之間的差異,難以為COVID-19的治療方案提供幫助。因此,開展宿主內動力學的建模研究十分必要。這類微觀水平的建模研究相對較少。Kern等[9]提出了一個改進的靶細胞限制模型(target cell limited model),以研究不同藥物方案和治療時機對SARS-CoV-2 病毒的影響;Blanco-Rodríguez等[10]使用免疫應答模型(immune response model)來定量研究COVID-19患者體內T細胞和病毒載量的變化,指出細胞的病毒清除率過低會導致疾病的惡化;Ghosh[11]根據兩名COVID-19患者的病毒載量數據擬合數學模型,研究了抗病毒藥物和疫苗接種對SARS-CoV-2病毒的影響。然而,經典的靶細胞限制模型和免疫應答模型均只能對感染與免疫中的一個方面進行描述,無法完整地刻畫患者體內的動力學變化。此外,大部分研究都側重于分析相關藥物的治療效果,而非分析不同癥狀患者之間的內在差異,難以找出導致急癥的原因。因此,構建一個新的宿主內動力學模型來分析COVID-19患者之間的內在差異是一個重要且有意義的研究課題。

為了研究SARS-CoV-2病毒及相關細胞在患者體內的動力學變化,分析急癥患者與重癥患者之間存在的差異,本文基于標準的靶細胞限制模型和免疫應答模型提出了一個TCL-IR(target cell limited-immune response)模型,并推導出模型的宿主內再生數表達式。首先,結合急癥、重癥患者體內T細胞與自然殺傷細胞(natural killer cell,NK cell)的實際數據進行模型擬合,并根據擬合結果對兩類患者之間的差異進行分析,找出導致急癥的相關原因。然后,從擬合優度、可識別性這兩個方面對模型進行分析,以進一步驗證模型的準確度與可信度。該研究可為類似的微觀動力學建模以及相關治療措施的制定提供一定的參考。

1 模型構建

本章先對已有的靶細胞限制模型[12]和免疫應答模型[13]進行介紹,然后提出改進的TCL-IR模型。

1.1 靶細胞限制模型

靶細胞限制模型與一般傳染病模型較為類似,將宿主內的相關細胞分為易感靶細胞和感染細胞兩類。不同的是,靶細胞限制模型的感染過程由病毒驅動,而病毒不會發生狀態轉移。模型的表達式如下:

靶細胞限制模型描述的過程為:病毒V侵入人體后以感染率β對靶細胞T進行感染并將其轉化為感染細胞I。被感染的細胞I以速率α進行病毒的增殖與釋放,同時因病毒效應和宿主免疫反應以速率δ進行凋亡。此外,病毒顆粒的清除速率為c。

1.2 免疫應答模型

免疫應答模型包含了生物學中邏輯斯諦增長的思想,主要考慮免疫T細胞和病毒之間的相互影響及數量變化。為了與上述靶細胞限制模型中的靶細胞數量T進行符號區分,這里將免疫T細胞的數量用符號D進行表示。模型的表達式如下:

免疫應答模型描述的過程為:

1)在病毒V尚未侵入人體時,體內的免疫T細胞服從穩態方程sD=δDD,其中sD表示免疫T細胞的補充速率,δDD則表示免疫T細胞的正常凋亡速率。

1.3 TCL-IR模型

盡管靶細胞限制模型和免疫應答模型已被應用于多種疾病的研究[15～17],但是依舊存在一些不足。一方面,靶細胞限制模型沒有考慮到免疫反應的影響,即病毒的增加會刺激免疫細胞進行快速增殖,從而導致病毒清除速率逐漸增大;病毒的減少會使免疫細胞增殖速率降低,從而導致清除速率逐漸減小。另一方面,免疫應答模型雖然考慮了免疫T細胞和病毒載量之間的相互影響,但是沒有考慮到靶細胞、感染細胞以及先天性免疫的NK細胞,難以對整個感染過程進行一個全面的描述。此外,相關研究表明,T細胞和NK細胞均是通過作用于感染細胞并使其凋亡來實現免疫反應的,而非直接作用于病毒本身[18～19]。由于沒有考慮感染細胞,免疫應答模型也不能對這一過程進行很好地刻畫。

因此,本文將免疫反應融入靶細胞限制模型,提出了一種改進的TCL-IR模型。在靶細胞限制模型的基礎上,增加了免疫應答模型中的T細胞數量變化表達式以及靶細胞的穩態方程sT=δTT。對于NK細胞,則參照了文獻[10]中的做法,即將NK細胞數據作為感染細胞表達式中的已知數據。此外,假設感染細胞和病毒都是均勻分布的,即被清除的感染細胞數量所占的比例等于被清除的病毒載量所占的比例。這種通過比例來間接計算病毒清除量的方式,能夠對T細胞和NK細胞的真實免疫過程進行更好地描述。相關研究指出,病毒在患者體內具有一定的壽命[15],因此模型引入了病毒的自然衰亡。最后,假設病毒的復制速率服從免疫應答模型中的邏輯斯諦函數,而非靶細胞限制模型中的線性函數,以體現病毒增長速率會受到宿主體內有限資源的制約。

綜上所述,TCL-IR模型描述的宿主內動力學過程如下:

1)在人體未感染病毒時,T細胞D和靶細胞T分別服從各自的穩態方程sD=δDD、sT=δTT。

3)在人體內的病毒被清零后,T細胞D和靶細胞T會逐漸恢復回穩態sD=δDD、sT=δTT。

根據以上過程,可得圖1所示的TCL-IR模型示意圖,圖中實線箭頭表示數量轉移,虛線箭頭表示影響。

圖1 TCL-IR模型Fig.1 The TCL-IR model

TCL-IR模型的表達式如下:

模型符號說明:

1)變量 D、T、I、V、N 分別表示免疫 T 細胞(CD8+T細胞和CD4+T細胞)、靶細胞、感染細胞、病毒、NK細胞的數量;

3)sT表示靶細胞的補充速率,δT表示靶細胞、感染細胞的自然凋亡率,β是病毒的感染率;

4)δ是感染細胞因免疫反應而凋亡的平均速率;

5)p是病毒的復制系數,K是病毒的最大環境容納量,c是病毒的自然衰亡速率。

1.4 TCL-IR模型的宿主內再生數

基本再生數(basic reproduction number)是傳染病建模研究中的一項關鍵參數,表示一個感染者在易感人群中能夠產生的下一代感染者的平均數量。基本再生數大于1意味著傳染病的蔓延,反之則意味著傳染病的消亡。類似地,可以定義出宿主內再生數(within-host reproduction number),以衡量病毒在細胞之間的傳播風險?；谙乱淮仃嚪椒╗20],本節對TCL-IR模型的宿主內再生數進行推導。

TCL-IR模型在形式上可視為媒介-宿主模型(vector-host model)的一種變體。同時,考慮到宿主內再生數衡量的是病毒在感染初期的傳播風險,所以T細胞和NK細胞的數量可視作常數D(1)、N(1)。基于此,可將模型進行簡化,具體如圖2和公式(4)所示。

圖2 模型簡化示意圖Fig.2 The model simplification diagram

根據各個倉室的新增感染數量、移出數量、非感染移入數量,得到矩陣F和W:

矩陣F和W在無病平衡點A處的雅可比矩陣分別為:

因此,通過計算矩陣FjWj-1的譜半徑ρ(FjWj-1),可得到TCL-IR模型的宿主內再生數為:

由表達式可知,宿主內再生數相當于病毒的復制速率與清除速率之比。當其大于1時意味著復制速率更快,病毒在體內得以復制和傳播;當其小于1時意味著清除速率更快,病毒會在體內逐漸消亡。需要說明的是,TCL-IR模型的宿主內再生數不含感染率β,這是因為TCL-IR模型中的感染是病毒驅動的,而非感染細胞驅動。而且,從再生數的定義來看,宿主內再生數衡量的是病毒的產生與消亡,這與推導結果一致。盡管宿主內再生數不包含感染率β,但是這并不能否定參數β的意義,因為β的取值會對靶細胞的感染規模產生影響。

2 實驗數據與參數估計方法

2.1 實驗數據

本文使用的數據來自中國武漢707名COVID-19患者的血液樣本[21]。這些早期的COVID-19患者按照癥狀嚴重程度被劃分為中癥(moderate)、重癥(severe)和急癥(critical)3組,從第3天開始每隔5 d檢測1次血液樣本直至第68天結束。數據主要包含了這些COVID-19患者的T細胞(CD8+T細胞和CD4+T細胞之和)、CD8+T細胞、CD4+T細胞、NK 細胞的數量平均值(圖3)。從圖3可以看出,中癥患者的T細胞數量變化并不明顯,難以挖掘出有用信息。因此,本文僅選用急癥、重癥患者的T細胞數據和NK細胞數據進行研究。由于每隔5 d才檢測1次血液樣本,為了得到每一天的NK細胞數量,本研究使用分段線性擬合的方法對NK細胞的數據進行估計和補全。Blanco-Rodríguez等[10]也對相同的數據進行了研究,因此在后續的內容中會進行印證與對比。

圖3 免疫細胞數量數據源自參考文獻[21],橫坐標表示血液樣本的檢測時間。Fig.3 The number of immune cellsData is from reference[21].The x-axis depicts the detection time of the blood samples.

2.2 自適應Metropolis算法

本文使用自適應Metropolis算法(adaptive Metropolis algorithm,AM algorithm)[22]對模型中未知的參數進行估計。AM算法是一種改進的MCMC(Markov chain-Monte Carlo)方法,通過將模型中需要估計的參數視作隨機變量進行抽樣來得到參數的估計結果。算法的具體步驟如下:

1)根據先驗分布隨機產生初始參數樣本x0,并任意指定初始協方差陣C0。

2)在時刻n,將正態分布N(xn-1,Cn)作為建議分布,產生候選樣本xn。方差Cn的更新公式為:

3)根據概率α來接受或拒絕候選樣本xn。

其中,P(xn|d)表示參數的后驗概率,P(d|xn)表示似然函數,P(xn)表示先驗概率。

4)重復步驟2、3直至產生足夠數量的樣本,去掉燃燒期后的剩余樣本集合即為抽樣結果,以此得到參數的估計值。

在上述算法步驟中,需要給定參數的先驗分布和似然函數?？紤]到先驗信息不足,本文選擇均勻分布作為無信息的先驗分布,并取似然函數為:

其中,m是觀測數據的數量,δ是擬合數據與實際數據的相對誤差。

3 實驗結果與模型分析

3.1 實驗結果

為了讓研究更具可信度、更符合生物學與免疫學方面的認識,本文首先根據相關研究和文獻來固定模型中的部分參數,使其具有一定的生物學意義,然后再對其余無法查證的參數使用AM算法進行估計。對于研究結果和結論,同樣使用文獻印證的手段來增強其說服力。

根據已有的研究與認識,SARS-CoV-2病毒可感染人體內的多種細胞,但是主要的靶細胞為呼吸道上皮細胞[14,23]。Baccam等[24]估計出成年人體內的呼吸道上皮細胞數量約為4×108個,即T(1)=4×108。數據是從第3天開始的,根據Blanco-Rodríguez等[10]的做法令 D(1)=D(3)。感染初期的病毒載量可能低于檢測水平,根據文獻[14,25]建議取V(1)=0.31。假設感染細胞的初值為I(1)=3。同樣,根據文獻[14],模型中病毒和T細胞的增殖參數取 K=108、kD=1.26×105。Perelson 等[26]指出病毒的平均存活時間約為8 h,因此本文假設病毒自然衰亡速率為c=2.4。T細胞的半衰期被認為是34～255 d[27～28],所以取T細胞的凋亡速率δD=0.01,并根據穩態方程sD=δDD得到補充速率sD。

因此,模型中需要估計的參數為:T細胞增殖系數r、靶細胞自然凋亡速率δT、感染率β、免疫凋亡速率δ、病毒復制系數p。根據急癥、重癥患者的相關數據對以上參數進行估計,結果如表1所示。根據參數估計值,可得急癥、重癥患者的宿主內再生數分別為1.094、1.140。這說明初期的免疫細胞數量不足以遏制病毒的傳播,必須進行快速增殖來提高病毒的清除速率。同時,得到的模型模擬結果如圖4所示。其中離散點表示實際數據,實線表示擬合數據。

表1 參數估計結果Table 1 Parameter estimation results

從T細胞的擬合結果(圖4A)中可以看出,急癥患者T細胞的快速增殖于30～40 d開始,在40～50 d達到峰值;重癥患者T細胞的快速增殖于20～30 d開始,在30～40 d達到峰值。在快速增殖的過程中,急癥患者與重癥患者的曲線斜率較為接近,表明兩類患者的T細胞數量增長率非常接近。但是,由于急癥患者的T細胞數量水平低于重癥患者,因此急癥患者體內要進行更劇烈的增殖活動才能保證T細胞增長率的接近。從參數估計結果中也能看出,急癥患者的T細胞增殖系數為0.652,約為重癥患者的3倍。同時,急癥患者的免疫凋亡速率為4.137×10-6,高于重癥患者的3.283×10-6。這些結果說明急癥患者體內的免疫反應比重癥患者的更為劇烈。

從病毒載量的模擬結果(圖4D)中可以看出,急癥患者的病毒載量在35 d左右達到峰值,約為7.20×104copies/mL,而重癥患者的病毒載量在28 d左右達到峰值,約為1.38×105copies/mL。這與Blanco-Rodríguez等[10]的研究結果較為接近。根據圖4D和C可知,兩類患者的病毒載量和感染細胞數量在達到峰值后都迅速下降,并在10 d左右接近清零。這種快速清零的結果與已有研究[10～11]一致。模擬結果還顯示,急癥患者的病毒載量低于重癥患者,這可能是因為急癥患者的免疫反應更為劇烈。相關研究指出,在癥狀更為嚴重的患者體內,病毒載量可能更低[10,29]。此外,從圖4C中可知,病毒載量較低導致急癥患者的感染細胞數也較少。這些結果說明更高的病毒載量雖然會增大靶細胞的感染規模,但是患者本身未必會出現更嚴重的癥狀。

圖4 模型模擬結果橫坐標表示血液樣本的檢測時間。Fig.4 The model simulation resultsThe x-axis depicts the detection time of the blood samples.

綜合以上結果,本文認為部分患者急癥癥狀的出現與自身過度的免疫反應有一定關系。在急癥患者體內,T細胞的數量水平較低,而T細胞增殖系數和免疫凋亡速率較高。這說明急癥患者體內通過更劇烈的免疫反應來彌補自身T細胞數量的不足。然而,這種過度的免疫反應可能會使患者表現出更為嚴重的癥狀甚至死亡。Hernandez-Vargas等[25]通過研究指出,老年患者病癥的嚴重程度并非直接歸因于病毒載量,而應歸因于宿主免疫反應;Wang等[4]發現,免疫系統為應對感染而釋放的細胞因子風暴可導致敗血癥并使患者死亡;Huang等[30]也通過實驗數據推測,疾病的嚴重程度與細胞因子風暴有一定關系。這些研究在一定程度上證實了過度的免疫反應確實會導致急癥的產生。這樣的結論不僅為COVID-19急癥患者的治療提供了一定的參考與建議,也肯定了接種疫苗的積極作用:

1)在治療急癥患者時,不僅需要使用抗病毒藥物,還應進行免疫調節治療,以減輕過度免疫反應導致的其他并發癥。

2)急癥患者體內劇烈的免疫反應與較低的T細胞數量水平有關。這意味著可以通過檢測T細胞數量水平來預測可能成為急癥的患者,從而提前做好相應的準備和治療方案。

3)提前接種疫苗能夠有效降低患者病情發展成為重癥、急癥的概率。因為提前接種疫苗能夠提高免疫細胞的數量水平,這不僅能使免疫系統對SARS-CoV-2病毒做出更快的反應,還能使免疫反應更為平穩地進行,從而防止過度免疫反應帶來的危害。

3.2 擬合優度分析

盡管從前文模型構建的角度來看,TCL-IR模型對宿主體內病毒和細胞的刻畫更為全面,但是依舊需要將其與已有模型進行實際比較,以說明TCL-IR模型的準確度。赤池信息量準則(Akaike information criterion,AIC)是一種常用的模型比較標準,用于描述模型的擬合優度。AIC值越低意味著模型對數據的描述越好。以AIC值為標準,將 TCL-IR 模型與 Blanco-Rodríguez 等[10]使用的模型進行比較,具體計算公式如下:

其中,R表示均方根誤差,yi表示觀測數據,表示對應的估計值,m表示未知參數的個數,n表示數據點的個數。

根據公式得到的模型比較結果如表2所示。

表2 模型的AIC值Table 2 The AIC values of the models

其中,文獻模型1是僅考慮了CD8+T細胞的標準免疫應答模型,也是文獻[10]選取的最優模型;文獻模型2是進一步考慮了CD4+T細胞和NK細胞而提出的改進版免疫應答模型,盡管考慮得更為全面,但是其AIC值反而較高,意味著對數據的描述較差。本文的TCL-IR模型不僅考慮了CD4+T細胞和NK細胞,而且AIC值比文獻中的兩個模型都低,表明TCL-IR模型對數據的描述更好。除了AIC值,由于免疫應答模型僅考慮了免疫反應,所以無法推導出宿主內再生數,而TCL-IR模型是融合了免疫反應的靶細胞限制模型,可以推導出宿主內再生數,從而能夠衡量病毒在體內的傳播風險。綜上所述,相較于已有模型,本文所提出的TCL-IR模型能夠更好地對數據進行描述。

3.3 可識別性分析

利用動力學建模的方法對實際問題進行研究,通常是以模型的參數估計值為基礎進行分析。這就意味著模型的參數估計需要從生物學、數學兩個方面的合理性進行考慮。為了考究參數估計值在生物學上的合理性,本文主要通過文獻印證和理論支撐的手段,來證實大部分參數和研究結果的可信度。為了考究參數在數學上的合理性,除了要讓模型盡可能地擬合數據外,還應該對當前數據和參數下的模型進行可識別性分析。

定義(可識別性[31]) 如果數學模型中的參數能夠根據數據唯一確定,那么這個數學模型就是可識別的。

當數據的數量不足或類別單一時,即使是理論上可識別的數學模型,在實際問題中可能依舊難以識別,這種問題在宏觀和微觀的建模研究中普遍存在[32～33]。由于本文模型僅依靠T細胞數據作為擬合標準,而NK細胞數據只是作為模型擬合的輸入數據,因此其可能會出現可識別性問題。為了判斷能否根據當前數據對模型參數進行識別,這里以急癥患者的主要參數值為例,使用輪廓似然法[34](profile likelihood method)進行分析,結果如圖5所示。由于本文使用的是最大似然估計,而非最小二乘估計,因此輪廓似然法的結果曲線呈凸形時意味著對應參數可被識別,且結果曲線的峰值對應參數的最優估計。

在圖5A、B中,參數r、δT的結果曲線先增大后減小,表明參數能夠被識別。在圖5C、D中,由于縱坐標跨度較大,曲線的最大值在圖中并不明顯。事實上,圖5C在4×10-6附近有一個最大值,圖5D在5附近有一個最大值,表明參數δ、p也能被識別,且本文估計的 δ=4.137×10-6、p=4.909均在結果曲線的最大值附近。盡管在當前情況下參數能夠被識別,但是考慮到參數敏感性和圖5C、D中存在的多個局部極值,在今后的研究中使用更多更全面的數據來進一步增強研究的可信度十分必要。

圖5 輪廓似然法結果橫坐標表示參數取值的變化。Fig.5 The results of the profile likelihood methodThe x-axis depicts the change of parameter values.

4 結論

基于靶細胞限制模型和免疫應答模型,本文構建了一個TCL-IR模型,以定量分析COVID-19的宿主內動力學。通過結合COVID-19患者的免疫細胞數據,本文使用AM算法對模型的未知參數進行了估計,并根據實驗結果分析了急癥、重癥患者之間的內在差異。結果顯示,急癥患者體內的T細胞增殖系數和免疫清除速率明顯高于重癥患者,但是其病毒載量反而較低。根據結果,本文認為部分急癥患者的病情加重并不是因為體內過高的病毒載量,而是因為自身過激的免疫反應,該推測得到了已有研究的印證。此外,考慮到當前數據的局限性,使用不同地區、不同時期的數據來進一步增強研究的實用性是我們未來的重要工作。