一种基于GPU的心脏电生理模拟方法与流程

本发明涉及心脏电生理领域，尤其涉及一种基于gpu的心脏电生理模拟方法。

背景技术：

近几十年，越来越多的研究者开始转向心脏电生理建模，以获取虚拟心脏进行研究。现如今对心脏建模有两种方法：一是对心脏整体进行建模，将完整心脏视为整体。这种建模方法优势就是简单且计算量少，但不关心模拟精确度只展示大致的变化情况。二是从细胞出发，建立心脏细胞模型，再由细胞模型建立心脏模型。这种方法能准确反应心脏的各个时刻的变化状态，具有极高的精度，但计算量十分巨大，模拟时间长。由于心脏研究对精度具有高要求，因此大部分的模拟都是采用第二种方法。由于从细胞建模需要大量的计算时间，因此对模拟进行加速是必不可少的。

目前主要的加速方法有两种：一是利用自适应时间步长法进行加速；二是使用图形处理器(gpu)进行加速。

自适应时间步长方法的原理是利用细胞间动作电位变化率来调节时间步长，进而加快模拟速度。心脏细胞不同时刻的动作电位的变化率是不同的。自适应时间步长方法就是通过利用电位变化率，对电位变化率大的给予较小时间步长，而电位变化率小的给与较大时间步长，通过放大部分时间步长来加快模拟速度。下一时刻的细胞电位通过常微分方程(ode)进行求解。现有的自适应时间步长主要有两种：一是victorri在1985年提出的thm方法，二是min-hungchen等在2016年提出的ccl方法。这两种方法在时间步长跨度为100(例如0.001ms-0.1ms)时，在百万规模的模型中分别能有20和80倍的加速效果。虽然使用自适应时间步长方法有了几十倍的加速，但对于心脏模拟依旧需要大量的时间。人类心脏细胞数量是亿级，这需要更长的计算时间，不能满足实际需求。

而gpu具有大量的计算单元(alu)，因此常被用于并行计算。而心脏包含大量的心脏细胞，用并行计算能有效加快模拟。现有的单个gpu加速一般在100-200倍，随着gpu数量的增多，所需的计算时间也会减少，理论上可以通过堆叠gpu达到一个临床接受的模拟时间，但这需要大量的gpu资源。

而将上述两个方案进行结合，即在gpu上运行自适应步长的方案，加速效果会出现减弱的现象，主要是由于gpu中的执行是以warp为单位同时执行(一个warp通常为32个线程)，也就是这32个线程同时开始执行也同时结束，而先执行完的线程也必须等待该warp中尚未执行完的线程执行完毕才会结束，进而开始执行下一warp，这就会导致线程等待，而自适应时间步长会给予不同的线程不同的时间步长，这就进一步加剧了线程等待问题，导致自适应时间步长方法在gpu上加速效果减弱。

技术实现要素：

本发明实施例提供一种基于gpu的心脏电生理模拟方法，能解决自适应时间步长方法在gpu上执行时，心脏电生理模拟加速效果减弱的问题，进一步提高心脏电生理模拟时的加速效果。

本发明一实施例提供一种基于gpu的心脏电生理模拟方法，包括：

gpu获取若干心脏细胞的细胞数据，并将所有所述细胞数据排列为二维数组；其中，各所述心脏细胞的细胞数据相同，且所述二维数组每一行中细胞数据的个数为32的整数倍；

所述gpu将位于同一行的细胞数据设置相同电流刺激值，继而按行为主序根据每一细胞数据以及所述电流刺激值，通过预设的自适应时间步长方法对每一心脏细胞进行心脏电生理模拟计算；其中，在模拟计算时每一初始细胞数据由所述gpu通过一个线程进行处理。

进一步的，每一所述细胞数据包括：初始电位、离子浓度、离子通道门变量、温度、最大时间步长、最小时间步长、最大变化电位和最小变化电位。

进一步的，所述预设的自适应时间步方法包括：thm或ccl。

可选的，通过以下步骤根据thm对每一心脏细胞进行心脏电生理模拟计算：

步骤a:根据所述离子浓度计算心脏细胞的总离子电流；

步骤b:通过以下公式计算心脏细胞当前时刻的时间步长：

；

其中，dt为心脏细胞当前时刻的时间步长，dv为电位变化率，且dv＝iion；iion为所述总离子电流；dvmax为所述最大变化电位；dvmin为所述最小变化电位；dtmax为所述最大时间步长；

步骤c:判断心脏细胞各时刻的时间步长之和是否大于所述最大时间步长，若是则执行步骤d，若否则执行步骤e；

步骤d:计算扩散电位，继而根据所述扩散电位更新所述心脏细胞当前时刻的细胞电位，并执行步骤f；

步骤e:重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a、b和c；

步骤f:判断是否达到预设模拟时长，若是则结束模拟计算，若否则将心脏细胞各时刻的时间步长之和置零，重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a、b和c。

可选的，通过以下步骤根据ccl对每一心脏细胞进行心脏电生理模拟计算：

步骤a1:根据所述离子浓度计算心脏细胞的总离子电流；

步骤b1:通过以下公式计算心脏细胞当前时刻的时间步长：

其中，dt为心脏细胞当前时刻的时间步长，dv为电位变化率，且dv＝iion；iion为所述总离子电流；

步骤c1:判断心脏细胞各时刻的时间步长之和是否大于所述最大时间步长，若是则执行步骤d1，若否则执行步骤e1；

步骤d1:计算扩散电位，继而根据所述扩散电位更新所述心脏细胞当前时刻的细胞电位，并执行步骤f1；

步骤e1:重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a1、b1和c1；

步骤f1:判断是否达到预设模拟时长，若是则结束模拟计算，若否则将心脏细胞各时刻的时间步长之和置零，重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a1、b1和c1。

通过实施本发明的实施例具有如下有益效果：

本发明实施例提供了一种基于gpu的心脏电生理模拟方法，首先gpu获取若干心脏细胞的细胞数据，然后排列成一个二维数组，并且二维数组中每一行的细胞数据的个数为32的整数倍，即二维数组的列数为32的整数倍，而此时不同心脏细胞的细胞数据都是相同的，然后位于同一行的细胞数据设置相同电流刺激值。这样使得二维数值中同一行的细胞数据对应的细胞在模拟时处与同一电位状态，并且每一个细胞数据由gpu分配一个线程进行处理，这样由于每一行细胞的个数是32的整数倍，而gpu在执行时一个warp中有32个线程，这就使得gpu按行为主序进行处理时能将每一行的细胞数据分入整数个warp进行执行，且每个warp内的各细胞数据对应的心脏细胞都处于同一个状态，这样在通过自适应时间步长方法进行模拟计算时，能够保证每个warp中每个线程都拥有同样的时间步长，从而解决了由于自适应时间步长给不同的线程以不同的时间步长，进一步加剧了线程等待，导致自适应时间步长方法在gpu上加速效果减弱的问题，提高了心脏电生理模拟时的加速效果。

附图说明

图1是本发明一实施例提供的一种基于gpu的心脏电生理模拟方法的流程示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明一实施例提供了一种基于gpu的心脏电生理模拟加速方法，包括：

步骤s101：gpu获取若干心脏细胞的细胞数据，并将所有所述细胞数据排列为二维数组；其中，各所述心脏细胞的细胞数据相同，且所述二维数组每一行中细胞数据的个数为32的整数倍；

步骤s102：所述gpu将位于同一行的细胞数据设置相同电流刺激值，继而按行为主序根据每一细胞数据以及所述电流刺激值，通过预设的自适应时间步长方法对每一心脏细胞进行心脏电生理模拟计算；其中，在模拟计算时每一初始细胞数据由所述gpu通过一个线程进行处理。

对于步骤s101：在进行心脏电生理模拟时，每个模拟对象即每个心脏细胞都对应有一细胞数据，且一细胞数据内包含了若干个数据项；优选的上述细胞数据内包含的数据项包括：初始电位、离子浓度、离子通道门变量、温度、最大时间步长、最小时间步长、最大变化电位和最小变化电位。

初始时由cpu对每个细胞数据的各个数据项进行初始化；这里需要说明的是，在初始化时，将不同心脏细胞的细胞数据都设为相同，即各细胞数据之间对应的数据项的值都是相同的，例如细胞数据a中的初始电位，与细胞数据b中的初始电位是相同的，细胞数据a中的离子浓度与细胞数据b中的离子浓度是相同的，细胞数据a中的离子通道门变量与细胞数据b中的离子通道门变量是相同的……

每个细胞数据的各数据项在cpu中初始化后，由cpu传输至gpu中，具体的cpu通过cudamalloc()申请显存空间，而后通过cudamemcpy()传递数据到gpu中。gpu接收到所有心脏细胞的细胞数据之后，将所有细胞数据排列为一个二维数组，数组中的一个元素代表一细胞数据，同时数组中每一行的元素个数(即细胞数据个数)为32的整数倍。

对于步骤s102:gpu将位于同一行的细胞数据设置相同电流刺激值，这样首先在数据初始化时，所有心脏细胞数据都是相同，如果给位于同一行的细胞数据设置相同的电流刺激值后，位于同一行的细胞数据所对应的心脏细胞都会在后续的模拟中，保持一个相同或大体相同的电位状态；对于不同行的细胞数据可以设置不同的电流刺激值；由于每个细胞数据由gpu分配一个线程进行处理，且每一行的细胞数据的个数是32的整数倍，这样就可以使得gpu按行为主序进行处理时，保证每个warp中各细胞数据对应的心脏细胞处于同一状态；

以下列举一实例对上述进行进一步的说明：

例如假设一个二维数组有5行32列，即说明每一行有32个细胞数据；开始时对每一行的细胞数据设置相同的电流刺激值，优选的不同行的设置不同电流刺激值；例如对于第一行所有的细胞数据设置的电流刺激值为a；第二行设置的电流刺激值为b；第三行设置的电流刺激值为c；第四行设置的电流刺激值为d；第五行设置的电流刺激值为e；在不同电流刺激值下各心脏细胞后续在同一时段下的电位状态会不一致。由于同一行的细胞数据设置相同的电流刺激值，且各细胞数据在初始化时都是相同的，那么就会使得同一行的细胞数据对应的心脏细胞的电位状态在后续的各个时刻都保持一致(或大体一致，有时会因为扩散电位的原因导致电位有细微的不同，但不会有太大的影响)；这样当cpu按行为主序依次对，每一行的细胞数据进行处理，那么此时上述二维数据的每一行的细胞数据会被分入一个wrap中进行处理(一个wrap有32个线程，每个线程对应处理一个细胞数据，而上述示例的二维数组每行有32个细胞数据)；上述5行32列的二维数组中的细胞数据就会被分为5个wrap中进行处理，而每个wrap中的每个线程所处理的细胞数据所对应的心脏细胞都处于相同的电位状态；因此当cpu在根据预设的自适应时间步长方法对每一心脏细胞进行心脏电生理模拟计算时，对于每个wrap中的每个线程都具有相同的时间步长，从而解决了由于每个wrap中的每个线程所处理的细胞数据对应的心脏细胞的电位状态不同，自适应时间步长会给予不同的线程不同的时间步长，加剧了线程等待问题，导致加速效果减弱的问题。

在一个优选的实施例中，自适应时间步长方法包括：thm或ccl。

在一个优选的实施例中，通过以下步骤根据thm对每一心脏细胞进行心脏电生理模拟计算：

步骤a:根据所述离子浓度计算心脏细胞的总离子电流；

步骤b:通过以下公式计算心脏细胞当前时刻的时间步长：

；

其中，dt为心脏细胞当前时刻的时间步长，dv为电位变化率，且dv＝iion；iion为所述总离子电流；dvmax为所述最大变化电位；dvmin为所述最小变化电位；dtmax为所述最大时间步长；

步骤c:判断心脏细胞各时刻的时间步长之和是否大于所述最大时间步长，若是则执行步骤d，若否则执行步骤e；

步骤d:计算扩散电位，继而根据所述扩散电位更新所述心脏细胞当前时刻的细胞电位，并执行步骤f；

步骤e:重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a、b和c；

步骤f:判断是否达到预设模拟时长，若是则结束模拟计算，若否则将心脏细胞各时刻的时间步长之和置零，重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a、b和c。

对于步骤a，在模拟计算时，根据例离子浓度使用rush-larsen方法计算出所有离子通道的离子电流，然后将所有离子通道的离子电流进行加总，从而获得心脏细胞的总离子电流，在这里需要说明的是，在第一次计算总离子电流时，还要将额外将设置的电流刺激值加上。后续迭代计算时就只需要计算各离子通道的离子电流然后进行加总即可。

对于步骤b,首先计算电位变化率，由于当前时刻的细胞电位v(tn)＝dt*iion+v(tn-1),因此有：dv＝iion；dv为电位变化率，iion为步骤a求得的总离子电流；v(tn)为当前时刻的细胞电位，v(tn-1)为上一时刻的细胞电位

然后通过以下公式计算当前时刻心脏细胞的时间步长：

dt为心脏细胞当前时刻的时间步长；dvmax为所述最大变化电位；dvmin为所述最小变化电位；dtmax为所述最大时间步长。

对于步骤c:将心脏细胞在各个时刻的时间步长进行加总，然后判断是否大于最大时间步长，如果是则执行步骤d计算扩散电位，如果不是则执行步骤e。

对于步骤d:细胞间电位相互传递，即当前细胞的电位受周围细胞电位影响，同时也影响周围的细胞。

使用拉普拉斯算子计算扩散电位：

在二维平面中：其中x、y为二维平面坐标轴。

在三维空间中：其中x、y、z为三维空间坐标轴。

将计算得到的扩散电位与原当前时刻的细胞电位进行加总，作为更新后的当前时刻的细胞电位即：v1(tn)＝v(tn)+diff。v1(tn)为更新后的当前时刻的细胞电位。

对于步骤e：通过rush-larsen方法重新计算离子浓度。

对于步骤f:判断是否达到预设模拟时长，若是则结束模拟计算，将电位数据回传至cpu；

若否，则将心脏细胞各时刻的时间步长之和置零，重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a、b和c。

上述实施例是通过thm对每一心脏细胞进行心脏电生理模拟计算，实现模拟加速的方法示例；具体效果如下表所示：

上表中，“优化前”的数据项，表示直接在gpu上使用thm方法模拟600ms心脏动作电位所需计算时间，“优化后”的数据项，表示通过实施本发明的实施例在gpu上使用thm方法模拟600ms心脏动作电位所需计算时间，可以看出计算时间明显减短了，加速效果更好。

在本发明另外一个实施例中，也可以根据ccl对每一心脏细胞进行心脏电生理模拟计算,步骤如下:

步骤a1:根据所述离子浓度计算心脏细胞的总离子电流；

步骤b1:通过以下公式计算心脏细胞当前时刻的时间步长：

其中，dt为心脏细胞当前时刻的时间步长，dv为电位变化率，且dv＝iion；iion为所述总离子电流；

步骤c1:判断心脏细胞各时刻的时间步长之和是否大于所述最大时间步长，若是则执行步骤d1，若否则执行步骤e1；

步骤d1:计算扩散电位，继而根据所述扩散电位更新所述心脏细胞当前时刻的细胞电位，并执行步骤f1；

步骤e1:重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a1、b1和c1；

步骤f1:判断是否达到预设模拟时长，若是则结束模拟计算，若否则将心脏细胞各时刻的时间步长之和置零，重新计算所述离子浓度，并根据更新后的离子浓度重复执行步骤a1、b1和c1。

其中步骤a1、c1、d1、e1、f1与上一实施例中的步骤a、c、d、e、f是相同的不再进行解释说明，主要对步骤b1进行说明。

对于步骤b1：ccl主要考虑电位二阶导数，提高模拟精度。

首先计算电位变化率dv：由于当前时刻的细胞电位v(tn)＝dt*iion+v(tn-1)因此有：dv＝iion。(相同的符合定义与之前的实施例相同，不再赘述)

由二阶泰勒展开式可得：

其中

有

上式为关于dt的二次方程，求解可得

则当前时刻细胞电位：v(tn)＝v(tn-1)+dt*iion。

具体效果如下表所示：

上表中，“优化前”的数据项，表示直接在gpu上使用ccl方法模拟600ms心脏动作电位所需计算时间，“优化后”的数据项，表示通过实施本发明的实施例在gpu上使用ccl方法模拟600ms心脏动作电位所需计算时间，可以看出计算时间明显减短了，加速效果更好。

通过实施本发明实施例能解决了由于自适应时间步长给不同的线程以不同的时间步长，进一步加剧了线程等待导致自适应时间步长方法在gpu上加速效果减弱的问题，提高了心脏电生理模拟时的加速效果。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。