CN113303815A - 一种多导经颅时间干涉电刺激电流参数的优化方法和系统 - Google Patents

Abstract

本发明是一种多导经颅时间干涉电刺激电流参数的优化方法和系统，用GPU并行化求解经颅时间干涉电刺激在每种电刺激参数下、在每个脑网格节点处的耦合电场强度，包括以下步骤：步骤1、在个性化大脑网格中提取目标脑区的网格特征，将其他部分作为无关脑区进行标记；步骤2、构建电刺激参数选择池，为经颅时间干涉电刺激构建电极位置和参数的组合矩阵，作为暴力求解算法的电刺激参数选择池；步骤3、利用GPU在每种电刺激参数下、在目标脑区每个网格格点上求解“空间两向量投影最小值问题”；步骤4、再次利用GPU求解低频耦合电场强度模长，得到多个优势电极参数；步骤5、通过增加电极数目，对预定比例的总电流进行再分配，优化出最终的参数组合。

Description

一种多导经颅时间干涉电刺激电流参数的优化方法和系统

技术领域

本发明属于认知神经科学技术领域，涉及脑部电刺激的仿真计算，尤其涉及一种个体化的多导经颅时间干涉电刺激电流参数的优化方法和系统。

背景技术

传统的经颅电刺激按照电流的特性可以分为直流、交流和随机电刺激。

本发明涉及的经颅时间干涉电刺激是自2017年以来提出的具有深部刺激特异性的新兴方法。在生理学作用原理方面，经颅电刺激一方面可以改变神经元的静息膜电位，可以促进和抑制动作电位的发放。另一方面可以影响大脑的内源性节律震荡，影响与之相关的神经活动。经颅时间干涉电刺激的作用原理与此两种情况类似。经颅电刺激效果和受试者个体化的大脑结构和目标脑区的位置非常相关。因此需要根据目标脑区个体化的最大化外加电场强度。

经颅时间干涉电刺激是利用两个相差固定频率的高频交流电刺激叠加耦合而成(例如2010和2000Hz)。两者的频率差是预估起效的低频电场频率(10Hz)。由于神经元的结构特点，其对高频的外加电场的反应是微弱的。而耦合而成的电场包络线是低频的，此低频的外加电场被认为是可以影响神经元的功能结构的。其相对于传统经颅电刺激的主要优势在于低频电场的起效位置。传统经颅电刺激由于电场由输入电极形成，因此不可避免地在靠近电极的脑部位置出现更强的外加电场强度。这使得传统经颅电刺激仅仅在刺激“表面”脑区的时候有比较清晰的效果和机制。如果是“深层”目标脑区则不可避免地混淆其电流经过的浅层脑组织。经颅时间干涉电刺激的低频电场由于是耦合而成，其最大强度可以出现在大脑深部，因此克服了传统经颅电刺激的上述弊端。但是如何设计电流使得低频电场能够最适当的出现在目标脑区。这需要进行个体化的仿真优化。这也是本发明系统所解决的最主要问题。

经颅时间干涉电刺激的仿真优化，有一些现有的解决方式。但是其核心问题是求解低频电场的包络振幅。而此包络振幅的求解是一个“空间两向量投影最小值问题”，见公式(1)。其解决的是指物理问题是在某个点上，两个高频电场如何求包络线的场强。此两种电场抽象成数学问题可以转化为两个3维向量的叠加问题。求包络问题可等价于求这两个向量在空间中任意方向内最小值的最大值。这个问题包含了最小值的求解，是数学上非光滑的，因此不能使用迭代法，只可以使用暴力求解。而这个解在脑部的每个网格节点和每种电极参数下都要求解。因此具有相当大的求解运算量。因此本系统发明首个目标就是使用GPU编程框架加速计算过程。另外之前提出的计算优化指标比较单一，仅包含峰比值，并不能很好的合理涵盖所有脑部节点。因此本系统还提出了其他合理优化指标。

关于多导电极参数的设计。原理上经颅时间干涉电刺激应使用四个电极。每两个电极通过一种高频电流。但是通过电阻分离信号源，可以比较容易的在硬件上实现电压浮动，从而在每一种高频电流内部，使用多个电极以改变其整个大脑内部的电场分布。同时耦合而成低频电场分布也会得到改变。但是重新暴力计算多个电极的最优化问题会带来难以完成的计算量。

发明内容

本发明是一种多导经颅时间干涉电刺激电流参数的优化系统和方法，利用GPU高并行化暴力求解经颅时间干涉电刺激在每种电刺激参数下、在每个脑网格节点处的耦合电场强度，利用多种合理的优化指标进行设计多导经颅时间干涉电刺激参数，本发明给出一种依次至多分配30％电流给予新增电极的方案。由于耦合特性，每对新增电极各包含一种高频成分。

本发明的技术方案为：一种多导经颅时间干涉电刺激电流参数的优化方法，包括以下步骤：

步骤1、个性化的从大脑中提取目标脑区的网格特征，将其他部分作为无关脑区进行标记。个性化大脑网格来源于受试者个体头部磁共振影像。将脑电的10-10系统中的76个电极作为备选电极。为每个电极和参考电极间之间利用有限元方法计算单位电流下的引导场电场分布。

步骤2、构建电刺激参数选择池。以上述76个电极，总电流设为2mA(实际应用不超过4mA)，电流分配最小精度为0.05mA为条件，为经颅时间干涉电刺激构建电极位置和参数的组合矩阵，作为暴力求解算法的电刺激参数选择池。

步骤3、利用GPU在每种电刺激参数下、在目标脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长。并依据网格提及进行加权平均。筛除电场强度小于预期阈限的参数组合。

步骤4、再次利用GPU在上一步筛除后的余下电刺激参数下、在无关脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长。并依据网格提及进行加权平均。依靠最大化目的脑区相对无关脑区的指数比例优化出最好的电极参数组合。

步骤5、从经过优化的一定数目的“优势电极参数”中分配一定比例(30％)的电流给与多导成对电极，逐步进行多导电刺激参数的优化计算。实现从4导、6导到8导的电极参数优化求解。

根据本发明的另一方面，还提出一种多导经颅时间干涉电刺激电流参数的优化系统，用于实现上述方法，包括：

磁共振设备，用于对受试者扫描个体头部磁共振影像；

脑电10-10系统，包括64个电极外加外周12个电极，共76个电极作为备选电极；

GPU，用于并行化求解经颅时间干涉电刺激在每种电刺激参数下、在每个脑网格节点处的耦合电场强度，为每个电极和参考电极间之间利用有限元方法计算单位电流下的引导场电场分布。

本发明的有益效果：

本发明基于受试者个体化的大脑三维结构和待干预认知活动对应的目标脑区，为经颅时间干涉电刺激提供了最优化的电刺激参数方案，包括电极的位置选取和各电极的电流强度。具有较快的计算速度和合理化的计算优化指标。并且还可以提供包含多导(6或8)的电刺激参数方案，进一步强化优化结果。

并且为了增加数据的有效利用率，减少内存读取时间。还利用了GPU中特有的并行规约和常量内存技术。极大的提高了计算速度，是CPU的100倍左右；虽然整体上都遵循目标脑区更大、无关脑区更小的电刺激原则。但是考虑到外加电场的均匀性、聚焦性、能量等维度，本发明给予了多种指标的计算窗口。

附图说明

图1为本发明一种多导经颅时间干涉电刺激电流参数的优化方法流程图；

图2为经颅时间干涉电刺激原理示意图；

图3为基于依据脑电10-10系统的76个备选电极位置分布图；

图4为本发明中以脑岛为目标区域，不同电极数目的无关脑区强度优化结果对比图；

图5为本发明中以脑岛为目标区域，不同优化指标的无关脑区强度优化结果对比图。

具体实施方式

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

具体结果可以请参阅图1所示一种多导经颅时间干涉电刺激电流参数的优化方法流程、图4所示为不同电极数目的无关脑区强度优化结果对比，图5所示为不同优化指标的优化结果对比图。

具体原理请参阅图2为经颅时间干涉电刺激原理示意图、图3为基于依据脑电10-10系统的76个备选电极位置分布图，本发明中高并行通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长。用于求解空间两向量投影最小值问题的低频电场耦合的计算公式(1)如下：

其中，E₁，E₂为在该点处两个高频电场的电场强度，为3维矢量，单位为V/m，

表示该点的三维坐标，α表示E₁和E₂的夹角，等式左边

为低频电场的耦合强度模长，字母上方的箭头表示矢量。

根据本发明的一个具体应用的实施例，采用一种多导经颅时间干涉电刺激电流参数的优化方法，其步骤如下：

步骤1、利用simNIBS软件从受试者头部的结构核磁共振影像中提取头部各组织3D模型。然后利用simNIBS软件建立个性化的头部立体四面体网格。以统一模板的脑岛部分为目标脑区，将其在受试者头部空间中进行对齐操作。在该大脑网格中标记脑岛作为目标脑区，将其他大脑灰质区域部分作为无关脑区进行标记。同时在头皮表面按照脑电的10-10系统中选取64个电极，并且在外围再增加12个电极，共76个电极作为备选电极位置。每个电极按照直径1.5cm，厚度3mm进行设计。利用Matlab软件依次对每个备选电极和参考电极(头顶处Cz电极)间之间利用有限元方法计算单位电流下的引导场电场分布。在有限元计算中头部各组分电导率被定义为头皮0.465S/m,头骨0.008S/m，脑脊液1.654S/m，大脑灰质0.275S/m，大脑白质0.126S/m。有限元网格使用的为1阶四面体形函数。边界条件为第一类边界条件。求解方法为基于GPU的无预条件器的共轭梯度下降法。

步骤2、构建电刺激参数选择池。以上一步头皮表面76个电极，总电流设为2mA(实际应用不超过4mA)，电流分配最小精度为0.05mA为条件，为经颅时间干涉电刺激构建电极位置和参数的组合矩阵，作为暴力求解算法的电刺激参数选择池。在测试数据中电极位置共有11546775

种情况，大脑网格节点包括25～30万个。

步骤3、利用GPU在每种电刺激参数下、在目标脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长。并依据网格提及进行加权平均。筛除电场强度小于预期阈限的参数组合。

在GPU计算中为增加计算速度，将电极的电流强度比例作为常量内存进行存储。为减少数据的吞吐量，增加可并行的数据规模，在每个线程块内部使用了共享内存和并行规约进行计算。同时为了减少合并访问冲突，大脑网格节点按照128的倍数进行零补充，使得GPU中每个线程束warp单元可以访问连续的内存数据。

筛选参数池时使用目标脑区强度最大值不低于0.2V/m或者平均强度不低于0.1V/m。这可以根据目标脑区的面积大小和电刺激要求自行选择。如果使用平均强度指标，那么其强度被每个节点代表的四面体体积加权平均。

步骤4、再次利用GPU在上一步筛除后的余下电刺激参数下、在无关脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长。并依据网格提及进行加权平均。为最大化目标脑区和最小化无关脑区强度。这里提供目的脑区相对无关脑区的最大值比值或者平均值指数比值作为优化指标。包含了电场均值和电场能量值等电学常见指标。一般来讲，越高的幂指数选择会带来更均匀的无关电场分布，符合一般实验的预期。

在图5中我们可以观察到无关脑区在幂指数(power)为4时，引起最小的激活。(此图中所有目标脑区激活强度最大值均保持不小于0.2V/m。)

在GPU计算过程中为减少数据的吞吐量，增加可并行的数据规模，在每个线程块内部使用了共享内存和并行规约进行计算。同时为了减少合并访问冲突。大脑网格节点按照128的倍数进行零补充，使得GPU中每个线程束warp单元可以访问连续的内存数据。

步骤5、经过前述优化可以得到4导(两对)电极的电机参数。在每个高频电流内部，可以增加电极数目以优化电场分布。计划每次增加一对电极，其中各有一个电极属于一个高频成分。为减少求解规模，将上一步的4导的各电机参数计算结果按照优化比值排序。从中选取前100种的“优势电极参数”进行优化。分配一定比例(30％)的电流给与新增的成对电极。电流最小分配精度与前相同可以提高10％左右的比值效果。降低和均匀无关脑区的电场分布。逐步地，可以实现从4导、6导到8导的电极参数优化求解。在图4中我们可以观察到无关脑区随着电极导数的增多(4→6→8)无关脑区激活逐渐降低。(此图中所有目标脑区激活强度最大值均保持不小于0.2V/m。)

在此步的电刺激参数池设计中，依照基尔霍夫定律，每种高频频率的电极中电流总和为0。故分配电流要考虑方向性，即是流入还是流出所属电极。每种高频频率在大脑中的电场这些电极的引导场电场依各通过的电流线性叠加计算而成。

在GPU计算中为增加计算速度，由于电流系数沿电极个数线性增长，因此不能一次性的将电极的电流强度作为常量内存进行存储。因此采取循环更新常量内存的方式。为减少数据的吞吐量，增加可并行的数据规模，在每个线程块内部使用了共享内存和并行规约进行计算。为了照顾各种GPU硬件，灵活地设置了可输入的双层循环。一层为新增电极位置。另一层为新增电流分配系数。用户可以自行选择合适的嵌套容量。每次计算同时为了减少合并访问冲突。大脑网格节点按照128的倍数进行零补充，使得GPU中每个线程束warp单元可以访问连续的内存数据。

综上，本发明求解过程利用高并发的图形处理器(GPU)作为计算系统，相对于传统中央处理器(CPU)求解显著提高计算速度，为暴力求解此问题提供了可靠的实现方式。本发明能够个性化的依靠受试者的头部影像结构和待电刺激靶脑区，设计给与最优化的经颅时间干涉电刺激电极位置和各电极的电流强度。本发明中提出的高幂指数的电场优化指标合理计算涵盖了大脑中的大量无关脑区的节点。并且通过多导电极的逐步修正，进一步明显减少无关脑区的电流消耗，提高了靶脑区电场强度的聚焦性。

本发明为多导经颅时间干涉电刺激的个体化实施提供了电流参数方案。解决了其计算资源需求庞大、计算结果容易出现奇异值等问题。其仿真刺激效果较传统的较少电极设计有明显提高。

以上内容仅仅是对本发明的构思所作的举例和说明，所属本技术领域的技术人员对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离发明的构思或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。

Claims (8)

1.一种多导经颅时间干涉电刺激电流参数的优化方法，其特征在于：利用GPU并行化求解经颅时间干涉电刺激在每种电刺激参数下、在每个脑网格节点处的低频耦合电场强度，包括以下步骤：

步骤1、在个性化大脑网格中提取目标脑区的网格特征，将其他部分作为无关脑区进行标记；所述个性化大脑网格来源于受试者个体头部磁共振影像；将脑电的10-10系统中的64个电极外加外周12个电极，共76个电极作为备选电极，为每个电极和参考电极间之间利用有限元方法计算单位电流下的引导场电场分布；

步骤2、构建电刺激参数选择池，以上述76个电极，总电流设为2mA(实际应用不超过4mA)，电流分配最小精度为0.05mA为条件，为经颅时间干涉电刺激构建电极位置和参数的组合矩阵，作为暴力求解算法的电刺激参数选择池；

步骤3、利用GPU在每种电刺激参数下、在目标脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场，将每个网格点的低频耦合电场场强和其网格体积进行加权平均，筛除电场强度小于预期阈限的参数组合；

步骤4、再次利用GPU在上一步筛除后的余下电刺激参数下、在无关脑区每个网格格点上通过求解“空间两向量投影最小值问题”，求解低频耦合电场强度模长，并依据网格体积进行加权平均，依靠最大化目标脑区相对无关脑区的指数比例优化出最好的电极参数组合，得到多个优势电极参数；

步骤5、通过增加电极数目，对预定比例的总电流进行再分配，进一步优化多导经颅时间干涉刺激的脑部电场分布，在多导电极带来的多个参数组合数目中，优化出最终的参数组合。

2.根据权利要求1所述的方法，其特征在于：依靠目标脑区在统一脑模板中的定位，在个体化脑部网格中选取其所属对应网格，同时将大脑灰质网格中其他部分标记为无关脑区。

3.根据权利要求1所述的方法，其特征在于：将脑电的10-10系统中的76个电极作为备选电极位置，将总电流设为2mA(实际应用不超过4mA)，以0.05mA为单位精度分配给任意4、6或8个电极中，将此组合矩阵作为暴力求解的参数池。

4.根据权利要求1所述的方法，其特征在于：对每个备选电极和参考电极，在给定电流的情况下计算其引导场电场分布，相同频率的多导电极电流构成的电场，由这些电极的引导场电场依电流线性叠加计算而成。

5.根据权利要求1所述的方法，其特征在于：对于每个网格节点，对于两种频率的电场强度都要处理其空间投影的最小值最大的优化问题；其计算过程在GPU的核函数中进行实现。

6.根据权利要求1所述的方法，其特征在于：提供了依据于幂指数的多种比例型优化目标，整体上都遵循目标脑区更大、无关脑区更小的电刺激原则，考虑到外加电场的均匀性、聚焦性、能量维度，给予多种指标的计算窗口。

7.根据权利要求1所述的方法，其特征在于：从经过优化的多个优势电极参数中分配预定比例的电流给予多导成对电极，逐步进行多导电刺激参数的优化计算，实现从4导、6导到8导的电极参数优化求解。

8.一种用于实现权利要求1-7之一的方法的多导经颅时间干涉电刺激电流参数的优化系统，其特征在于，包括：

磁共振设备，用于对受试者扫描个体头部磁共振影像；

脑电10-10系统，包括64个电极外加外周12个电极，共76个电极作为备选电极；

GPU，用于并行化求解经颅时间干涉电刺激在每种电刺激参数下、在每个脑网格节点处的耦合电场强度，为每个电极和参考电极间之间利用有限元方法计算单位电流下的引导场电场分布。

Images