CN107273655B - 一种人体心室肌组织的心力衰竭模型的设计方法 - Google Patents

Abstract

一种人体心室肌组织的心力衰竭模型的设计方法，属于虚拟人体心脏的技术领域。本发明以心力衰竭为研究对象，对心脏兴奋‑收缩耦合下心肌电生理、力学特征和螺旋波动力学机制进行探讨，旨在发现心衰与心律失常之间相互作用的关系。TP06模型仿真人体单个心室肌细胞，再将TP06模型与NHS08模型耦联，建立ECC模型；基于改进的TP06模型和现有数据建立人心室心衰肌细胞模型，比较正常和心衰情况下人心室肌细胞电生理和力学的特性、动作电位时程恢复曲线以及早期后除极恢复曲线；最终通过在细胞间加入耦联，成功构建二维单层心肌组织动作电位传导模型，并采用“切波法”诱发螺旋波；然后，比较了正常和心衰情况下的螺旋波动力学特征。

Description

一种人体心室肌组织的心力衰竭模型的设计方法

技术领域

本发明属于虚拟人体心脏的技术领域，具体为一种人体心室肌组织的心力衰竭模型的设计方法。

背景技术

2010年流行病学资料表明，全球慢性心力衰竭(Chronic Heart Failure,CHF,简称心衰)患者的数量已高达2250万，并且仍以每年200万的速度递增，我国成年人心衰的患病率为0.9％，并呈逐年上升趋势。

当心脏发生心衰时，其生理和结构上发生巨大的改变，改变因素之间相互影响，因此，传统上从临床和解剖学角度分析心衰对心脏整体电生理和力学功能特性的影响是有局限性的。随着计算机技术的发展，越来越多的研究人员开始利用计算机虚拟心脏模型,并在细胞层次和组织层次上进行定量分析。

心脏是集电生理学、动力学、血液流体力学以及神经、生化控制等于一身的极其复杂的综合系统，虽然在目前阶段不可能将所有的因素都考虑进心脏的建模仿真过程，但心脏收缩运动源于心室肌细胞电兴奋这种最基本耦合关系应该在建模仿真中体现出来，否则仿真结果可能会有很大的片面性。目前仿真人体心室肌电生理模型比较常用的是2006年Ten Tusscher K H W J和Panfilov A V在《American Journal of Physiology-Heart andCirculatory Physiology》上发表文章《Alternans and spiral breakup in a humanventricular tissue model》提出的TP06模型。

心室肌细胞力学模型较电生理模型起步晚，但是其对心力衰竭等心脏性疾病的影响是至关重要的，过去几十年,众多研究者致力于这方面的研究并构造了很多力学心脏模型,但基本上是简化的左心室模型,这是因为心脏的动力主要来自左心室，比较经典的左心室模型有：1990年Beyar R和Sideman S在会议《Computers in Cardiology 1989，Proceedings》上发表文章《Amodel of left ventricular contractionand t ransmuralventricul ar deformation》提出的数学解析模型；1995年Gucci one JM和Cos ta KD等人在《Journal of Biomechanics》上发表文章《Finite element stress analysis of leftventricular mechanics in the beating dog heart》提出的有限元模型。虽然这些模型在心脏方面上的研究取得了不少成果,但是这些模型的很大局限性在于没有考虑心室肌电兴奋与机械收缩之间的耦合关系。近几年，越来越多的研究者将兴奋-收缩耦合这种体现了心室肌电兴奋对心脏收缩运动的调控机制加入到心脏建模仿真过程中。2011年，LouisD.Weise和Alexander V.Panfilov在《PLOS ONE》上发表文章《New Mechanism of SpiralWave Initiation in a Reaction-Diffusion-Mechanics System》提出了在连续的心脏质加入“兴奋-收缩耦联(Excitation-Contraction Coupling，简称ECC)机制”的观点。通常建立兴奋-收缩耦合模型是通过一种细胞电生理模型和一种力学模型以强耦合或者弱耦合的方式相结合，弱耦合方法中隐含一个假设，即变形和电生理之间微小的耦合效应可忽略不计，但是1977年Lab和M.J.在《Progress in Biophysics&Molecular Biology》上发表文章《Mechanosensitivity as an integrative system in heart:an audit》提出细胞的电生理和机械特征在光照射下表现出强烈的耦合效应，所以弱耦合方法具有很大的局限性。而强耦合是利用数值计算方法将描述电兴奋、主动张力、变形和它们之间的反馈机制的方程在同一时间被计算，但是只能解决计算量大的问题，适合于简化和少维模型。因此，2008年Niederer S和Smith N在《Progress in Biophysics&Molecular Biology》中发表文章《Animproved numerical method for strong coupling of excitation and contractionmodels in the heart》提出“更新方法”来建立心室肌力学模型，简称NHS08模型，并解决了部分强耦合问题中出现的数值不稳定性等现象。2013年，Louis D.Weise和AlexanderV.Panfilov在《PLOS ONE》上发表文章《ADiscrete Electromechanical Model for HumanCardiac Tissue:Effects of Stretch-Activated Currents and Stretch Conditionson Restitution Properties and Spiral Wave Dynamics》提出为了提高计算效率采用强耦合的方式将心室肌电兴奋与收缩力相结合构建人体离散化的心脏模型。

使用ECC机制来建立心力衰竭的病态模型，可加深人类对该类疾病病理现象的认识，使用这样的模型可以帮助临床医生确定和治疗此类疾病，具有非常大的应用价值。目前国内外许多研究已经将ECC耦合模型用于建立心衰模型中，并取得了丰硕的研究成果。1999年，Winslows等人在《Circulation Research》上发表文章《Mechanisms of AlteredExcitation-Contraction Coupling in Canine Tachycardia-Induced Heart Failure,ⅡModel Studies》提出了狗的心室肌M细胞心衰模型，简称Winslow模型。2004年,在《Circulation》上发表文章《Rate dependence and regulation of action potentialand calcium transient in a canine cardiac ventricular cell model》提出了基于Winslow模型改进的HRd模型，此模型仿真了心衰的狗心室肌细胞模型。2007年，张宇在《中国生物医学工程学报》上发表文章《心室肌细胞电生理-力学复合模型在心力衰竭中的仿真研究》提出利用已有的实验数据修正兴奋-收缩耦合模型建立心力衰竭模型。2013年，臧蕴亮等人在《Journal of Molecular&Cellular Cardiology》发表文章《Theoreticalinvestigation of the mechanism of heart failure using a canine ventricularcell model:especially the role of up-regulated CaMKII and SR Ca(2+)leak.》提出基于HRd模型的心衰的狗心室肌细胞模型。

通过对目前的研究现状进行分析，发现大多数是以动物心脏作为研究对象，缺少人类心衰模型的建立，更未考虑计算的稳定性和高效性。

发明内容

本发明的目的在于克服现有心衰模型的不足，提出一种人体心室肌组织的心力衰竭模型，使得心衰模型可以仿真人体心衰生理病理特征，并提高计算的稳定性和高效性。

采用人类心室肌细胞电生理模型-TP06模型仿真人类心室肌细胞的电兴奋过程。然后将TP06模型与NHS08模型以强耦合的方式结合，形成人类心室肌细胞ECC模型，再结合已有的实验数据将ECC模型调整为人体的心力衰竭模型，最后利用特殊边界处理法实现二维心衰组织模型。

本发明的技术方案：

一种人体心室肌组织的心力衰竭模型的设计方法，步骤如下：

第一步，确定细胞刺激方案

给人体心室肌组织的单个细胞施加连续的n个刺激，前n-1个刺激每两个刺激的时间间隔不断的缩小，时间间隔为280～1000ms，最后一个刺激与前一个刺激的时间间隔为800～1000ms；刺激强度为-38～52μA/cm²，持续时间为1～2ms，刺激频率为1～2Hz；

第二步，建立细胞电生理模型

从数学角度来分析TP06细胞模型，将TP06细胞模型看成是多个非线性方程(微分方程)，非线性方程的数目即为TP06细胞模型中定义的变量个数；TP06细胞模型的实现变成给定初始值，求解非线性方程组的问题，非线性方程具体形式如下：

式中：τ是时间常数，即p从峰值降低到67％所需要的时间；p代表不同的门控变量，假设p有“开”、“关”两种状态，门从关到开和从开到关的速率分别为α和β，不同p的状态改变实现各离子电流，并引起膜电压V的一系列变化，即动作电位，t是计算机仿真时间；

由于心室肌细胞模型包含的变量较多，需要计算及存储的数据量很大，使用传统的一阶前向Euler算法求解非线性方程，兼顾计算精度和效率，时间步长为0.02ms～0.05ms。

第三步，建立细胞的ECC模型

细胞力学模型中的NHS08力学模型的计算方程如下：

T_ind表示肌小节的长度和速度不相关张力(主动张力水平)；T0表示肌小节在等长条件下对Ca²⁺存在依赖性(对于非最大激活条件)的张力；

T_ind＝T_ref(z/z_max)

T0＝T_inda

式中，z是可利用的肌动蛋白(TnC)位点分数，z_max是在给定肌小节长度下，可利用的肌动蛋白位点的最大的分数，a是肌小节长度和速度的尺度函数；

使用强耦合方法实现细胞电生理模型和NHS08力学模型耦联，形成ECC模型，即在相同的时间条件下，求解TP06细胞模型和NHS08力学模型中的非线性方程；使用一阶前向Euler算法求解ECC模型中所有非线性方程，时间步长为0.02～0.05ms；为了获得稳定的ECC模型，以1～2Hz的频率仿真多个电生理过程后，再记录细胞电生理模型和NHS08力学模型的人类心室肌细胞的电-力特性。

第四步，建立人类心室肌细胞的心衰模型

首先，将具有超慢失活和再激活特性的人类心室肌细胞的L-型钠电流I_NaL加入到ECC模型中，用I_NaL＝G_NaLm³h(V-E_NaL)计算，使ECC模型真实再现出人类心室肌细胞的钠电流的动态跨膜过程，式中m,h表示钠通道的激活、快速失活，G_NaL表示L-型钠通道电导，E_NaL表示L-型钠通道的平衡电势；然后对ECC模型的部分参数进行调整，如下表

第五步，建立人类心室肌组织的心衰模型

首先，将细胞之间兴奋的传导过程加入到ECC模型中，兴奋在二维心室肌组织中传导可使用扩散方程描述：

式中：

是空间梯度算子，x和y为组织的空间坐标，仿真中分别通过水平和竖直方向的二阶差商来表示，

项反映了兴奋膜电位在心室肌组织中的扩散传导过程，V(mV)是膜电压，C_m为单位细胞膜面积上的等效电容，D是扩散系数用来模拟闰盘电阻的作用；I_ion(μA/cm²)是单位细胞膜面积上受到离子浓度、膜电位和门控变量制约的总的细胞膜电流密度；I_st是单位细胞膜面积上的刺激电流强度。从数值计算的角度分析对扩散方程做近似，如下：

式中Δx和Δy称为空间步长,且可调。

然后，采用无通量(No-Flux)边界条件对二维心室肌组织进行边界处理，方程如下：

式中：x和y是组织的空间坐标，x_min和x_max分别表示x坐标的最小值和最大值，y_min和y_max分别表示y坐标的最小值和最大值，方程的求解过程如下：

沿x轴(变量为j)和y轴(变量为i)方向计算位于心室肌组织片边界处第n个细胞的膜电压(V)，兼顾计算精度和效率，细胞n位于组织边界上四个顶点时，计算公式如下

细胞n位于左右边界(即沿y轴方向变化)时,计算公式如下：

同理，细胞n位于上下边界(即沿x轴方向变化)时，计算公式为：

最后在仿真程序中实现多个细胞组成的正方形二维心室肌细胞组织块。

本发明的有益效果：本发明可以真实虚拟人体心室肌心衰组织，可以通过调节ECC模型中的单一或多因素的方式加深人类对心衰生理和病理情况的认识，还可以辅助临床医生确定心衰的治疗方案。同时本发明在计算可行性的前提下保证了仿真效果。

附图说明

图1组织边界示意图

图2为本发明提供的正常和心衰下人类心室肌细胞电生理特性曲线(a)动作电位(APD)曲线、(b)Ca²⁺浓度[Ca²⁺]_i、(c)Na⁺-Ca²⁺交换电流I_NaCa、(d)L-型钙电流I_CaL。

图3为本发明提供的正常和心衰下人类心室肌细胞力学特性的曲线(a)[Ca²⁺]_TRPN变化、(b)等长张力变化(T0)

图4为正常和心衰细胞的APD恢复曲线与基本循环时长(BCL)、舒张间隔(DI)关系图(a)正常模型的DI与动作APD恢复曲线，(b)CHF模型的DI与APD恢复曲线，(c)BCL与APD在50％和90％的去极化值恢复曲线

图5,6分别为心室肌组织30cm×30cm(600×600个细胞)和50cm×50cm(1001×1001个细胞)螺旋波动态图

图7给出螺旋波传导的组织片上的伪心电图P-ECG

表1 TPO6模型中变量初始值

表2 NHS08模型中变量初始值

具体实施方式

为了使本发明的目的、技术方案和优点更加清晰明了，下面结合二维心室肌组织(假设组织面积为30cm×30cm，空间步长为Δx＝Δy＝0.05cm,则仿真的心肌组织片包含601×601个细胞)螺旋波的产生和发展过程，同时获得心衰下的伪心电图P-ECG，给出具体实施方式并参照附图，对本发明作详细说明。本实施例是以本发明的技术方案为前提进行的，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

第一步，确定细胞刺激方案

给人体心室肌组织的单个细胞施加连续的10个刺激，刺激时间间隔依次为1000ms,800ms,600ms,550ms,500ms,450ms,400ms,300ms,290ms,1000ms；刺激强度为-38μA/cm²，持续时间为1ms，刺激频率为1Hz；

第二步，建立细胞电生理模型

从数学角度来分析TP06细胞模型，将TP06细胞模型看成是多个非线性方程(微分方程)，其中非线性方程的具体数目由TP06细胞模型中定义的变量个数决定；TP06细胞模型的实现变成给定初始值，求解非线性方程组的问题，方程具体形式如下：

式中：τ是时间常数，即p从峰值降低到67％所需要的时间；p代表不同的门控变量，假设p有“开”、“关”两种状态，门从关到开和从开到关的速率分别为α和β，不同p的动态改变实现各离子电流，并引起膜电压(V)的一系列变化，即动作电位，t是计算机仿真时间；

使用传统的一阶前向Euler算法求解TP06模型中的非线性方程，初始值见表1。兼顾计算精度和效率，时间步长为0.02ms，仿真时间t＝10000ms。

第三步，建立细胞的ECC模型

细胞力学模型(NHS08力学模型)的计算方程如下：

T0＝T_inda

式中,T_ind表示肌小节的长度和速度不相关张力(主动张力水平)；T0表示肌小节在等长条件下对Ca²⁺存在依赖性(对于非最大激活条件)的张力，z是可利用的肌动蛋白(TnC)位点分数，z_max是在给定肌小节长度下可利用的肌动蛋白位点的最大的分数，a是肌小节长度和速度的尺度函数；

使用传统的一阶前向Euler算法，在每个迭代的步长时间0.02ms内，求解TP06细胞模型和NHS08力学模型中的所以非线性方程，并将上一时刻获得变量迭代值作为下一时刻迭代初始值。TP06模型非线性方程中所有变量的初始值见表1，NHS08模型非线性方程中所以变量的初始见表2。为了获得稳定的ECC模型，以1Hz的频率仿真10个电生理过程后，再记录细胞电生理模型和NHS08力学模型的人类心室肌细胞的电-力特性，主要包括心室肌细胞的动作电位(Action Potential Duration,简称APD)、各离子电流和收缩力等特性。

第四步，建立人类心室肌细胞的心衰模型

首先，将具有超慢失活和再激活特性的人类心室肌细胞的L-型钠电流I_NaL加入到ECC模型中，使用I_NaL＝G_NaLm³h(V-E_NaL)计算方程，式中m,h表示钠通道的激活、快速失活，G_NaL表示L-型钙通道电导，E_NaL表示L-型钙通道的平衡电势；然后对ECC模型的部分参数进行调整，如下表

第五步，建立人类心室肌组织的心衰模型

首先，将细胞之间兴奋的传导过程加入到ECC模型中，兴奋在二维心室肌组织中传导可使用扩散方程描述：

式中：

是空间梯度算子，x和y为组织的空间坐标，仿真中分别通过水平和竖直方向的二阶差商来表示，

项反映了兴奋膜电位在心室肌组织中的扩散传导过程，V(mV)是膜电压，Cm＝1μF/cm²为单位面积的膜电容，D＝0.001cm²/ms是扩散系数用来模拟闰盘电阻的作用。I_ion(μA/cm²)是受到离子浓度、膜电位和门控变量制约的细胞膜电流密度，I_st是单位面积上的刺激电流强度。从数值计算的角度分析对扩散方程做近似，如下：

式中Δx和Δy称为空间步长,且可调，这里选取Δx＝Δy＝0.05cm。

然后，采用无通量(No-Flux)边界条件对二维心室肌组织进行边界处理，方程如下：

式中：x和y是组织的空间坐标，x_min＝0和x_max＝30分别表示x坐标的最小值和最大值，y_min＝0和y_max＝30分别表示y坐标的最小值和最大值，方程的求解过程如下：

沿x轴(变量为j)和y轴(变量为i)方向计算位于心室肌组织片边界处第n个细胞的膜电压(V)，兼顾计算精度和效率，细胞n位于组织边界上四个顶点时，Δx＝Δy＝0.05cm,顶点坐标分别为(0,0)、(0,600)、(600,0)和(600,600)，计算公式如下

细胞n位于左右边界(即沿y轴方向变化)时，i＝0或者600，n＝0～600,计算公式如下：

同理，细胞n位于上下边界(即沿x轴方向变化)时，j＝0或者600，n＝0～600,计算公式为：

第六步，使用切波法诱发螺旋波。

将二维组织平面按照空间步长划分网格坐标，则每个网格点代表一个细胞。将刺激施加在正常系统的左边界600×10的条形区域内，一列行波产生，在行波向前传播到组织中央处，波行至一半截断，行波的截断处由于存在电位梯度差，随着时间的演化，截断附近区域就会逐渐卷曲，最后形成螺旋波，该方法被称为“切波法”。

第七步，获得P-ECG

分别在二维组织平面坐标A(100,100)和B(500,100)处放置点电极，在螺旋波传导的过程中，将其他点的膜电压(V)值×点到A(或者B)的距离导数叠加至A(或者B)处细胞点的V，再计算出A(或者B)点处的V值，从而获得两点的P-ECG，两者P-ECG相减可求差。

应该说明的是，本发明的上述具体实施方式仅用于示例性阐述本发明的原理和流程，不构成对本发明的限制。因此，在不偏离本发明精神和范围的情况下所做的任何修改和等同替换，均应包含在本发明的保护范围内。

表1 TP06模型中参数初始值

表2 NHS08模型中参数初始值

Claims (1)

1.一种人体心室肌组织的心力衰竭模型的设计方法，其特征在于，步骤如下：

第一步，确定细胞刺激方案

给人体心室肌组织的单个细胞施加连续的n个刺激，前n-1个刺激每两个刺激的时间间隔不断的缩小，时间间隔为280～1000ms，最后一个刺激与前一个刺激的时间间隔为800～1000ms；刺激强度为-38～52μA/cm²，持续时间为1～2ms，刺激频率为1～2Hz；

第二步，建立细胞电生理模型

从数学角度分析TP06细胞模型，TP06细胞模型即为多个非线性方程，非线性方程的数目即为TP06细胞模型中定义的变量个数；非线性方程具体形式如下：

式中：τ是时间常数，即p从峰值降低到67％所需要的时间；p代表不同的门控变量，假设p有“开”、“关”两种状态，门从关到开和从开到关的速率分别为α和β，不同p的状态改变实现各离子电流，并引起膜电压V的一系列变化，即动作电位，t是计算机仿真时间；设定时间步长为0.02ms～0.05ms；

第三步，建立细胞的ECC模型

细胞力学模型中的NHS08力学模型的计算方程如下：

T_ind＝T_ref(z/z_max)

T0＝T_inda

式中，T_ind表示肌小节的长度和速度不相关张力；T0表示肌小节的等长条件下对Ca²⁺存在依赖性的张力，z是可利用的肌动蛋白TnC位点分数，z_max是在给定肌小节长度下，可利用的肌动蛋白位点的最大的分数，a是肌小节长度和速度的尺度函数；

使用强耦合方法实现细胞电生理模型和NHS08力学模型耦联，形成ECC模型，即在相同的时间条件下，求解TP06细胞模型和NHS08力学模型中的非线性方程；使用一阶前向Euler算法求解ECC模型中所有非线性方程，时间步长为0.02～0.05ms；为了获得稳定的ECC模型，以1～2Hz的频率仿真多个电生理过程后，再记录细胞电生理模型和NHS08力学模型的人类心室肌细胞的电-力特性；

第四步，建立人类心室肌细胞的心衰模型

首先，将具有超慢失活和再激活特性的人类心室肌细胞的L-型钠电流I_NaL加入到ECC模型中，用I_NaL＝G_NaLm³h(V-E_NaL)计算，使ECC模型真实再现出人类心室肌细胞的钠电流的动态跨膜过程，式中m,h表示钠通道的激活、快速失活，G_NaL表示L-型钠通道电导，E_NaL表示L-型钠通道的平衡电势；然后对ECC模型的部分参数进行调整，如下表

第五步，建立人类心室肌组织的心衰模型

首先，将细胞之间兴奋的传导过程加入到ECC模型中，兴奋在二维心室肌组织中传导使用扩散方程描述：

式中：

是空间梯度算子，x和y为组织的空间坐标，仿真中分别通过水平方向和竖直方向的二阶差商来表示，

项反映了兴奋膜电位在心室肌组织中的扩散传导过程，V是膜电压，C_m为单位细胞膜面积上的等效电容，D是扩散系数，即用来模拟闰盘电阻的作用；I_ion是单位细胞膜面积上受到离子浓度、膜电位和门控变量制约的总的细胞膜电流密度；I_st是单位细胞膜面积上的刺激电流强度；从数值计算的角度分析对扩散方程做近似，如下：

式中：Δx和Δy称为空间步长，且可调；

然后，采用无通量边界条件对二维心室肌组织进行边界处理，方程如下：

式中：x和y是组织的空间坐标，x_min和x_max分别表示x坐标的最小值和最大值，y_min和y_max分别表示y坐标的最小值和最大值，方程的求解过程如下：

沿x轴和y轴方向计算位于心室肌组织片边界处第n个细胞的膜电压，其中：沿x轴的变量为j，沿y轴的变量为i，兼顾计算精度和效率，细胞n位于组织边界上四个顶点时，计算公式如下：

细胞n位于左右边界时，即沿y轴方向变化，计算公式如下：

同理，细胞n位于上下边界时，即沿x轴方向变化，计算公式为：

最后实现多个细胞组成的正方形二维心室肌细胞组织块。

Images