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)的实际数据进行模型拟合,并根据拟合结果对两类患者之间的差异进行分析,找出导致急症的相关原因。然后,从拟合优度、可识别性这两个方面对模型进行分析,以进一步验证模型的准确度与可信度。该研究可为类似的微观动力学建模以及相关治疗措施的制定提供一定的参考。

本章先对已有的靶细胞限制模型[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.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.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.

基于靶细胞限制模型和免疫应答模型,本文构建了一个TCL-IR模型,以定量分析COVID-19的宿主内动力学。通过结合COVID-19患者的免疫细胞数据,本文使用AM算法对模型的未知参数进行了估计,并根据实验结果分析了急症、重症患者之间的内在差异。结果显示,急症患者体内的T细胞增殖系数和免疫清除速率明显高于重症患者,但是其病毒载量反而较低。根据结果,本文认为部分急症患者的病情加重并不是因为体内过高的病毒载量,而是因为自身过激的免疫反应,该推测得到了已有研究的印证。此外,考虑到当前数据的局限性,使用不同地区、不同时期的数据来进一步增强研究的实用性是我们未来的重要工作。