乙肝免疫球蛋白阻断乙肝病毒在母婴间传播的动力学模型

李冬梅, 刘伟华, 汪 琪, 郭美静

(1- 哈尔滨理工大学理学院，哈尔滨 150080; 2- 黑龙江省医院儿科，哈尔滨 150036)

1 引言

乙型病毒(Hepatitis B Virus, HBV)性肝炎，简称乙肝，是由HBV 感染引起的传染病.母婴传播是HBV 传播的三种方式之一.实施母婴传播阻断策略能够显著降低HBV 感染，且已经成为现阶段我国乙肝预防工作的重点[1].应用乙肝免疫球蛋白(Hepatitis B Immunoglobulin, HBIG)可有效的阻断HBV 通过胎盘对胎儿的感染，进一步提高出生人口的健康质量[2].

众多数学工作者通过研究HBV 的发病机制和治疗机理建立了数学模型，并希望通过定量研究HBV 变化规律对防治乙肝疾病带来一定的帮助[3–5].Lewin 和Ribeiro[6]通过数学模型进一步研究了抗病毒药物治疗对乙肝疾病的控制问题.Stanca 等人建立了具有抗体免疫反应的数学模型，这使得人们更加精确地刻画了宿主体内的HBV 变化规律[7].实际上，在临床医学中是通过观察母体内的体征情况来预测胎儿的患病情况.周煜等人[8]运用Meta 分析的方法对临床检验的数据分析来评价药物阻断HBV 母婴宫内免疫策略的可行性，发现HBV 携带孕妇注射HBIG 可降低新生儿的感染率，但目前对于母婴之间HBV 传播的定量模型研究结果很少.

本文综合考虑HBV 进入母体后的传播机理，以及药物对病毒的抑制作用，建立了母体怀孕期间的药物阻断HBV 在母婴间传播的动力学模型，研究了母体接受HBIG 来控制体内HBV 的数量，从而减少HBV 对胎儿的感染机率.再根据临床用药量数据，运用Matlab 模拟胎儿体内的病毒数量变化曲线，从理论上优化给药方案，为指导临床治疗提供了可参考的依据.

2 乙肝免疫球蛋白阻断乙肝病毒在母婴间传播的动力学模型的建立

研究发现，当个体感染了HBV 后，其作为抗原刺激宿主细胞的免疫系统产生抗体来阻滞感染HBV.若免疫系统不强无法抵抗HBV 复制，就会带来慢性持续性感染，从而演变成乙肝患者.乙肝免疫球蛋白(HBIG)是预防HBV 入侵复制的被动免疫制剂，当人体接受这种外源性抗体后，可以在短期内中和并清除血清中游离的HBV.临床上通常采用HBIG 阻断HBV 在母婴间的传播，有关方面的研究已取得了可观的结果[9,10].这为建立数学模型提供了可依据的理论基础.在妊娠初期，由于胎盘的屏障作用，母体与胎儿血液分开，互不干扰，除了可以进行选择性的物质交换外，胎盘还可以阻止母体内的一部分病毒进入胚胎、保障胎儿正常生长发育.但在孕妇妊娠20 周后，胎盘开始传递母体抗体给胎儿的功能，但胎儿体内的抗体不足以抵抗母体传递的病毒，为此要控制母体内的病毒数量.通过应用HBIG 来降低母体内HBV 数量，可以有效阻断HBV 感染胎儿[7].由文献[6]从数学和医学角度做出如下假设：

1) 胎儿体内的HBV 主要是来源于母体内HBV 输入；

2) 在胎儿的生长发育过程中，其免疫系统尚不完善，不足以抵御外来病毒的入侵.为了降低胎儿的感染可能性，药物可以减少母体内HBV，从而阻断HBV 在母婴间的垂直传播.因而不考虑胎儿免疫系统功能；

3) 每个个体乙肝传播机理相同；

4) 认为HBV 携带者与感染者之间存在着临界值.若个体的HBV 超过临界值，可诊断为HBV 感染者.

基于上述原因与假设，将母体与胎儿看成两个个体，胎盘连接两个个体，胎儿是否感染主要是受母体HBV 量的影响.下面分别讨论母体及胎儿体内的易感肝细胞(Ti)、感染肝细胞(Ii)、HBV(Vi)和免疫抗体(C1)的情况，其中i= 1,2，下角标为1 代表母体、下角标为2 代表胎儿.考虑到HBIG 对HBV 的作用机理，建立了如下的动力学模型

模型(1)中相关参数意义，如表1 所示.

表1 模型(1)相关参数表

3 主要研究结果

在临床实践中，对母体应用HBIG 阻断HBV 在母婴之间传播，药物可以有效控制母亲体内病毒的数量，使母体内的病毒不传播给胎儿.因此，先考虑母体在药物作用下病毒数量变化趋势，再分析胎儿体内HBV感染可能性.可将模型(1)分解成如下的母体模型、胎儿模型

将母体模型(2)的第一、二、四方程相加，令m=min{dT1,dI1,dC1}，则由

可知，存在t1，当t ≥t1时，有

再由模型(2)的第三个方程可求得，存在t2>t1，当t ≥t2时，有

且Γ 是模型(2)的不变集.

3.1 母体模型(2)的稳定性

模型(2)在平衡点处满足如下方程

由式(4)求得

令

将式(5)代入式(4)的第三个方程，得到

其中

由式(5)、式(6)易求得无病毒平衡点

由二次函数性质知，式(6)确定的函数为f(V1) = 0.当R0> 1 时，有唯一正根V ∗1> 0，使得f(V ∗1) = 0.再将V ∗1代入式(5)，可计算出模型(2)有病毒平衡点Q∗1=(T1∗,I1∗,V1∗,C1∗).

定理1当R0< 1 时，模型(2)存在局部渐近稳定的无病毒平衡点Q∗0；当R0>1 时，模型(2)无病毒平衡点Q∗0是不稳定.

证明 模型(2)的Jacobi 矩阵为

模型(2)在Q∗0处线性近似方程对应的特征方程为

当R0< 1 时，由Routh-Hurwitz 定理[11]知，式(8)的所有特征根均为负特征根，无病毒平衡点是局部渐近稳定点；当R0>1 时，无病毒平衡点是不稳定点.

考虑系统

其中X1∈1, X2∈, X=(X1,X2), X∗=(,0)为平衡点.假设下列条件：

H2:·X1=A1(X)(X1−)系统(9)的子系统存在唯一的全局渐近稳定点；

H3: 对任意的X ∈Γ，A2(X)是Metzler 矩阵且不可约；

H4: 定义在Γ 上的A2(X)存在一个上界矩阵，即存在∈Γ，使得

H5: 矩阵的谱半径α()<0.

引理1[12]如果系统(9)满足条件H1–H5，那么平衡点X∗是系统(9)的全局渐近稳定点.

为了研究模型(2)无病毒平衡点的全局稳定性，将模型(2)写成具有形式(9)的矩阵方程形式.令

定理2当R0<1 时，模型(2)的无病毒平衡点是全局渐近稳定点.

证明 下面验证引理1 的条件.

H1: 由于Γ 是模型(2)的不变集，因而模型(2)在Γ 内是耗散的.

H2: 对于子系统

·X1=A1(X)(X1−X∗1)，即有

显然

是子系统(10)的全局渐近稳定点.

H3: 对于矩阵A2(X)，非主对角线上元素均为正，故A2(X)是Metzler 矩阵且不可约[13].

H4: 由A2(X)的定义及式(10)的稳定性可知，A2(X)在Γ 内存在上界矩阵，取=(X∗1,0)，可得

H5: ¯A2是常数矩阵，要证α()<0，等价于α(D −CA−1B)<0[11].

由引理1 知，当R0<1 时，模型(2)的无病毒平衡点是全局渐近稳定点.

定理3当R0>1 时，模型(2)在Γ 内是一致持久的.

证明 由定理1 知，当R0>1 时，在区域Γ 内模型(2)存在着两个平衡点Q∗0和Q∗1，且Q∗0是不稳定点.根据式(8)计算出模型(2)的特征值知，只有∂Γ 上初值(T10,0,0,C10)出发的轨线最终趋于唯一的平衡点Q∗0，其它轨线都将远离Q∗0.因此Q∗0是∂Γ 上孤立的最大紧的不变集.因而，存在正的常数a，对于从在Γ 的内部出发的模型(2)解有

由此模型(2)的解不会跑出∂Γ，从而，有Ws({})∩◦Γ=ϕ，其中

综上，当R0>1 时，模型(2)在空间区域Γ 内是一致持久的[13].

现考虑系统

其中F为Rn的连续函数，D ⊂Rn是连续的开集，X ∈Rn.令X(t,X0)是系统(11)的解，初始条件为X(0)=X0，假设：

1) 系统(11)存在吸引子集K ⊂D；

2) 系统(11)在D内有唯一的平衡点X∗.

引理2[14]如果假设(1)、(2)成立，且存在Lyapunov 函数L(X)、有界可导函数G(X)、正常数a1、g和a2满足：

(i)a1|X|≤L ≤a2|X|;

(ii)L′(X)≤(G′(X)−g)L(X);

则系统(11)的平衡点X∗是全局渐近稳定点.

定理4当R0>1 时，模型(2)的有病毒平衡点Q∗1是全局渐近稳定点.

证明 由式(7)，可得模型(2)在Q∗1处线性近似方程对应的特征方程为

则由式(12)知，模型(2)的所有特征根均为负特征根.

由Routh-Hurwitz 定理[11]知，当R0>1 时，有病毒平衡点Q∗1是局部渐近稳定点.

再由式(7)可计算出模型(2)的第三加性复合矩阵为[15]

其中

a11=−b1V1−dT1−dI1−εC1−dV1, a22=−b1V1−dT1−dI1−kV1−dC1

a33=−b1V1−dT1−εC1−dV1−kV1−dC1, a44=−εC1−dI1−dV1−kV1−dC1.

式(13)对应的线性复合系统为

构造Lyapunov 函数

L(t,X,Y,W,Z)=max{C1|X|, V1|Y|, I1|W+Z|, I1|W|, I1|Z|}.

由定理3，在Γ 内存在a1>0, a2>0，使得

下面计算L的导数.

1) 若取L=C1|X|时，有V1|Y|≤C1|X|，可计算

其中G1=lnC1, g1=dT1+dI1+dV1.

2) 若取L=V1|Y|时，有C1|X|≤V1|Y|，则可计算

其中G2=lnI1+lnV1, g2=dT1+dC1.

3) 若取L=I1|W+Z|时，有V1|Y|≤I1|W+Z|，则可计算

其中g3=dT1+dI1+dC1.

4) 若取L=I1|W|时，有V1|Y|≤I1|W|，则可计算

其中g4=dT1+dC1.

5) 若取L=I1|Z|时，则可计算

其中G3=lnI1, g5=dI1+dV1+dC1.

由式(16)–(20)，可有

D+L ≤[G′−g]L,

其中G′=max{G1,G2}, g=min{g1,g2,g3,g4,g5}.

由引理2 知，模型(2)的有病毒平衡点Q∗2是全局渐近稳定点.

3.2 胎儿模型(3)的稳定性

在假设(1)的约定下，通过药物治疗后，母体内的HBV 最终趋于V ∗1，可以认为感染的母体输入给胎儿的HBV 量为¯η=ηV ∗1，由式(3)可得如下胎儿体内HBV 模型

同母体可行域Γ 的分析，可以求得方程(21)可行域为

其中

方程(21)的平衡点满足下列方程

整理得

将式(23)代入式(22)中第二个方程，得

由二次函数性质知，f(I2)=0 存在唯一正根，使得f()=0.再将代入式(23)中，可得到方程(21)唯一的有病毒平衡点.

定理5当R0>1 时，模型(21)存在唯一的有病毒平衡点是全局渐近稳定的.

证明 模型(21)的Jacobi 矩阵为

计算

显然有a1> 0, a3> 0，且a1a2−a3> 0.由Routh-Hurwitz 定理[11]知，有病毒平衡点是局部渐近稳定的.

模型(21)对应以ω为周期的任意周期解P(t) = (T2(t),T2(t),V2(t))的二阶复合系统为[15]

做Lyapunove 函数

由区域D内只含有稳定平衡点Q∗2，故存在c1，使得L ≥c1sup{|X|,|Y|,|Z|}.

将对(24)式求右导数，有

其中α=min{dT2,dI2}.

由式(25)得

其中

再由式(21)可知

由式(26)可得

取

由式(27)知

由此可知，方程(24)零解是稳定的.

另一方面，取H= diag(−1,1,−1)可证明的非对角线元素非正，模型(21)是竞争系统，即模型(21)满足Poincare-Bendixson 性质，由文献[11]知，模型(21)有病毒平衡点是全局渐近稳定.

4 模型应用

目前，用HBIG 阻断携带HBV 的母体感染胎儿的治疗方案没有明确统一的标准，但已取得了阶段性的进展.临床实践中，通过检测治疗过程中母体的HBV-DNA 数据(其与HBV 存在正相关关系)来确定治疗效果.并认为HBV 携带者与感染者之间存在着临界值，若检测结果没有超过临界值，可以视为HBV 转阴，无传染性可能.因此，用药原则遵循着母体内HBV 不应超过临界[7].根据检测到母体HBV-DNA 数据，可以定性推测胎儿感染HBV 可能性.用数学模型可以定量分析母婴体内HBV 数量变化趋势，预测胎儿体内HBV 数量的变化情况.并且通过调整模型的参数，进一步优选用药量使得母体的HBV 控制在临界值之下.其中HBV-DNA 的临界值为277，HBV 的临界值为1000[7].模型(1)中参数见表2[16].

表2 模型相关参数数值表

下面根据模型(1)结果和表2 中乙肝患者的相关参数选取，针对于HBIG 治疗HBV 孕妇案例，分析模型结果的合理性，数值模拟HBIG 剂量对母体内HBV 的数量影响，预测胎儿感染HBV 数量的变化情况.

案例：孕妇怀孕第20 周，孕期检查时发现孕妇体内携带HBV，检测HBV-DNA 没有超过临界值为277，认为还没有发病.为避免将病毒传染给胎儿，在家属知情的情况下，对孕妇进行注射HBIG 的阻断治疗.治疗期间检测到母体内HBV-DNA 含量如下表3.采取的治疗方案是标准用药方案，即每月注射一针HBIG，每次剂量200 IU，直至分娩[17].

表3 孕妇血清中HBV-DNA 含量

由表3 中数据可以绘制出治疗期间孕妇血清中HBV-DNA 含量的曲线如下图1.

由图1 可以看出，在怀孕20 周时开始用药，母体内HBV-DNA 含量先增加，到了第22 周左右时达到峰值，随着时间的推移，母体内HBV-DNA 含量逐渐减少，直至40 周后趋于平稳.

图1 母体内HBV-DNA 含量

4.1 模型的合理性分析

结合案例数据，用模型(1)数值模拟母体HBV 数量变化趋势，分析模型结果的合理性.选取母体内的HBV 数量的初值V1(20)<，将表2 中数据代入模型(1)，计算出R0= 1.05> 1，由定理4 知病毒平衡点是全局渐近稳定点.应用Matlab 软件对模型(1)分别数值模拟出母体、胎儿体内的HBV 数量随时间变化的曲线图像如下图2 和图3.

图2 当µ2 =200IU 时，母体内V1(t)随时间的变化曲线

图3 当µ2 =200IU 时，胎儿体内V2(t)随时间的变化曲线

从图2 可以看出，在20 周时接受治疗，HBV 数量增加到第22 周时达到峰值，随后逐渐降低.图2 的HBV 数量曲线变化趋势与图1 反映的HBV-DNA 含量曲线变化趋势基本相同.在治疗期间内模型给出HBV 数量变化趋势与临床实际监测HBV-DNA 含量变化趋势一致，说明模型是合理的.同时，由图2 的数据可知，母体内HBV 达到的峰值数量900，低于在治疗期间(20 周–45 周)平均HBV 数量为483，HBV 数量在可控的范围内，降低了胎儿的感染HBV 的风险.由图3 数据可计算出孕期胎儿的HBV 数量是362.

4.2 药物剂量对携带HBV 孕妇病毒数量的影响

当孕妇携带HBV 的数量没有超过临界值数量时，胎儿感染HBV 风险很大.临床上应用HBIG 阻断HBV 在母婴间传播，4.1 中给出了标准用药方案下孕妇、胎儿体内HBV 数量，能够有效控制HBV 数量不超出临界数值.在不超出最小中毒剂量下，当增加用药剂量至时400IU[17]，分析数值模拟HBV数量的变化趋势.选取母体内的HBV数量的初值V1(20) = 600<，其它参数仍取表2，计算出R0= 0.64<1.由定理3 知，无病毒平衡点是全局渐近稳定点.由于孕周期时间是有限的，在此只考虑45 周之内孕妇体内HBV 数量变化情况，用Matlab 数值模拟出母体、胎儿体内HBV 数量随时间变化的曲线图像如下图4 和图5.

由图4 数值可知，当改变药物剂量为时400IU，母体内HBV 数量峰值795 低于1000，平均病毒数量367.对比图2、图4，发现在孕期内HBV 数量均下降，说明调整剂量后的用药方案能更好的控制体内的病毒数量.由图5 数值可以计算出胎儿体内的HBV 的数量是340.由此说明增加药剂量可以快速降低母体内HBV 数量，胎儿的感染风险减小，达到更好的治疗效果.

图4 当µ2 =400IU 时，母体内V1(t)随时间的变化曲线

图5 当µ2 =400IU 时，胎儿体内V2(t)随时间的变化曲线

4.3 药物剂量对感染HBV 孕妇病毒数量的影响

图6 不同剂量下，母体内V1(t)随时间的变化曲线

图7 不同剂量下，胎儿体内V2(t)随时间的变化曲线

从图6 可以看出，随着药物剂量的增大，母体内HBV 的数量下降的速度越快，且在第12 周时HBV 数量分别为1050、830、810，说明给药剂量为200IU 时，不能严格的控制母体HBV 传染给胎儿，不能达到阻断的效果，通过剂量为400IU 与500IU 对比知，500IU 的母体内的HBV 数量下降得更快，则在临床实践中，在患者经济允许情况下，可适当地加大药物的治疗剂量，做出针对性的治疗方案.