经颅电刺激是一种极具潜力的无创神经调控技术。基于个体经颅电刺激电场仿真模型的电极优化方法可依据个体头部组织结构差异给出个性化电刺激参数，显著提升经颅电刺激的干预靶向性。然而，目前所提出的电极优化方法存在计算耗时较长（通常大于1 d）、未考虑刺激器输出通道数受限等局限，进一步阻碍了电极优化方法在临床中的应用。因此，本文提出了一种高效可操作的电极优化方法，可以结合目标脑区强度和聚焦性两个目标进行优化并约束所使用的电极数量，且计算快速。与目前常用的电极优化方法相比，本文所提方法可以在保证优化效果的基础上，极大地节省计算时间（耗时减小85.9%），并且可依据刺激器实际输出通道数设置优化的电极使用个数，将电流分散到多个电极，提升经颅电刺激的使用耐受性。本文提出的电极优化方法具有高效率和易操作的特点，可为个体精准经颅电刺激的实施提供一定的数据支撑和参考。
引用本文： 谢旭, 王敏敏, 张韶岷. 一种高效可操作的经颅电刺激电极优化方法. 生物医学工程学杂志, 2024, 41(4): 724-731, 741. doi: 10.7507/1001-5515.202308016 复制
0 引言
经颅电刺激（transcranial electrical stimulation，TES）技术是一种非侵入式的神经调控技术，主要通过放置于头皮表面的电极施加特定的微弱电流刺激，进而调控大脑的功能活动[1]。经颅电刺激已在多种神经精神类脑疾病（比如：抑郁症[2]、癫痫[3]、帕金森病[4]和记忆障碍[5-6]）中得到应用，针对健康人群也展示出可提高运动表现[7-8]和认知功能[9-10]的潜力。
传统的经颅电刺激使用两个比较大的矩形海绵电极来进行干预，其覆盖面积一般为20～35 cm2，因而在相关脑区产生的治疗场较弥散、聚焦性差，无法进行精准的靶向刺激[11-12]。为了更加精确地控制经颅电刺激在大脑中产生的电场分布，以达到更好的刺激效果，可使用多个较小的电极（面积为1～2 cm2）代替传统的大贴片电极对脑部进行刺激[13]。高精度经颅电刺激（high definition transcranial electrical stimulation，HD-TES）的电极使用方案是一种典型的使用小电极进行刺激的应用，其包括一个阳极电极放置在靶区的正上方，四个阴极施加均匀的电流大小放在阳极四周的位置[14]，但这种方案依旧无法针对个体差异进行精准优化。为了实现个体精准靶向刺激，现已提出一些相关的经颅电刺激电极优化方法，主要是在个体经颅电刺激电场建模的基础上，结合个体脑部组织结构差异，通过电极优化算法进行计算，以获得更好的电极位置组合和对应的电流大小，使得刺激的聚焦性得到明显的改善。
如今，已经提出多种应用于经颅电刺激的电极优化算法。最早出现的是基于最小二乘法的电极优化，这种方法可以产生相对集中的刺激，即刺激的聚焦性较好，但往往在靶区产生的电场强度较低[15]。而基于线性约束最小方差法的电极优化方法通过强制设定靶区的电场强度为一个固定值，然后最小化非靶区的总能量，可以一定程度上弥补强度的不足，但容易在非目标区域形成不期望的“旁瓣”[15]。此外，有研究提出了一种最大强度法，其目标是在靶区获得尽量大的电场强度，但完全忽略了电场聚焦性的问题，导致最终优化结果在大脑产生的电场非常弥散，且电流会集中在少数的电极上，造成刺痛感[16]。还有研究提出分布约束式最大强度法，此方法可以将注入的电流分布在多个电极上，从而降低每个电极注入的电流大小[17]。另有研究使用交替方向乘数法进行电极优化，可改善目标处电场分布的聚焦性，同时降低非目标区域的强度[18]。但上述方法均为单目标优化方法，每次只能针对靶区强度或聚焦性单个目标进行优化。而平衡聚焦性和强度的刺激法（stimulation with balanced focality and intensity，SBFI）能对两个目标同时进行优化，可以得到平衡电场强度和聚焦性两个目标的优化结果，但是该方法将所有电极都置于解空间，非靶区的电场能量的计算量非常大，整个优化过程非常耗时[19]。上述方法都无法控制使用电极的数量，优化结果所使用的电极数量可多达20个，但常规的电刺激设备的通道数通常在8个以内[20]，因此这些优化方法难以投入实际应用。针对电极数量较多的问题，Zhu等[21]使用传统的多目标优化算法——非支配排序遗传算法（nondominated sorting genetic algorithm，NSGAII）进行电极优化，可以获得帕累托前沿，用于从很多电场强度和聚焦性结果相平衡的解中进行选择，此方法可以控制电极使用数量。Wang等[22]先使用遗传算法对电场强度进行单目标搜索，然后结合遗传算法与粒子群算法进行多目标优化，以获得权衡目标的帕累托前沿。Ruffini等[23]首先使用遗传算法进行计算，确定电极的位置和数量，然后使用最小二乘法分配每个电极的电流。这些经典多目标优化算法和遗传算法都具有一个很大的缺陷：优化时间过长；以遗传算法为例，通常需要耗费几个小时。综上所述，目前所提出的电极优化算法都需要通过电场建模获取候选电极的电场转移矩阵，耗时长久（> 1 d），且大部分方法都未考虑刺激器输出电流通道数，难以用于临床实践。
针对目前已有的电极优化方法存在的问题，本文提出了一种高效可操作的经颅电刺激电极优化方法，可针对目标电场强度和聚焦性两个目标进行优化，在考虑使用电极数目限制的同时，大幅提高计算效率，以期提升电极优化方法的实用性。
1 材料和方法
1.1 个体磁共振数据获取
本研究招募了10名健康成人受试者进行T1加权（T1 weighted，T1w）和T2加权（T2 weighted，T2w）磁共振成像（magnetic resonance imaging，MRI）结构扫描，以获取个体的头部结构影像学数据用于电场建模与电极优化。图像采集地点位于浙江大学华家池校区中心大楼，受试者年龄范围为20～30岁，所有受试者均签署了知情同意书。本文研究内容通过了浙江大学伦理委员会批准，并获得授权使用相关数据，伦理审查编号为：浙大生仪会审[2022] 42号。
本文所用的结构性MRI采集设备为： MAGNETOM 3T扫描仪（Siemens Inc., 德国）和32通道头部线圈。MRI的扫描参数类型包括重复时间（time of repetition，TR）、回波时间（time of echo，TE）、翻转角（flip angle）、视野范围（fields of view，FoV）和体素大小（voxel size）。T1w的具体扫描参数为：TR = 2 300 ms，TE = 2.32 ms，flip angle = 8°，FoV = 240 mm × 240 mm，voxel size = 0.9 mm3；T2w的扫描参数为：TR = 3 200 ms，TE = 408 ms，FoV = 230 mm × 230 mm，voxel size = 0.4 × 0.4 × 0.9 mm3。
1.2 个体经颅电刺激电场仿真模型构建
在进行电极优化之前，首先需要对个体人脑进行有限元建模，构建头部体积导体模型，此模型也可以被称为电场转移矩阵，即A矩阵。A矩阵表示从电极在头皮注入的电流强度到在大脑中产生的电场强度的映射关系，如式（1）所示：
 |
其中，e代表不同位置的电场强度大小，s代表不同电极的电流大小。e和A的分量如式(2)所示：
 |
其中，r1，···，rN代表在头部模型上的体素位置，N是头部模型中体素的数量，M为电极的总数量，e(r1)，···，e(rN)和a(r1) ，···，a(rN)分别代表r1，···，rN位置上的电场强度大小和A矩阵分量。
电场建模的目的是为了模拟外部电刺激对脑区的影响，不考虑大脑内本身存在的电信号，因此可以将模型视为无源，当通过电极将电流施加到体积的边界时，体积中产生的电势分布V可以作为拉普拉斯方程的解[24]，如式(3)所示：
 |
其中，J是电流密度分布，
是组织的电导率，V是电势分布。符号“
”和“
∙”分别表示梯度和散度。
电场建模总体流程如图1所示。首先通过扫描受试者的头部获得T1w和 T2w结构影像MRI数据，将MRI分割成几部分，包括头皮、头骨、脑脊液、灰质和白质等，对应组织电导率参数分别设置：头皮为0.465 S/m、头骨为0.010 S/m、脑脊液为1.650 S/m、白质为0.126 S/m、灰质为0.276 S/m、空气腔为2.5 × 10-14 S/m、电极为5.9 × 107 S/m、导电膏为0.300 S/m [25-27]。然后，将电极按照国际标准导联10-10系统的规定放置在头部模型上。利用有限元对头部模型进行网格化，并获得A矩阵。使用计算得到的A矩阵进行电极优化，获得可在靶区产生良好电场强度和聚焦性的电极与电流分配方案。最后，对电极优化结果进行电场仿真，观察电场在头部模型上的分布情况。
图1
电场建模总体流程
Figure1.
Overall process of electric field modeling
1.3 电极优化方法
在个体头部模型上对电场强度和聚焦性两个目标进行优化。在最大化目标靶区的电场强度的同时，最小化非目标区域的电场强度，以提升聚焦性。为减少计算量，提高计算效率，本文将最小化非目标区域的电场强度表征为：首先选取头部其他位置离靶点的最远间距dmax，获得该最远间距的四分之一距离上的坐标点合集，而后最小化这些坐标点的电场强度。总的电极优化公式如式(4)所示：
 |
其中，s表示电极的电流大小，e0TCs表示目标处所需电场强度的分布，E0.25dmax表示与靶点距离为dmax的四分之一的坐标点的总电场强度。使用权重参数λ来平衡聚焦性和电场强度两个目标，λ越大，聚焦性越高；λ越小，靶区的电场强度越大。
为了保证刺激的耐受性，需要对电极优化方法施加相关安全限制：
第一，电极的输入电流大小要等于输出电流大小，即电流大小总和为0。第i个电极的电流用Si表示，如式（5）所示：
 |
第二，电流的绝对值总和不大于特定值Itotal，如式（6）所示：
 |
第三，单个电极的电流的绝对值不大于特定值Imax，如式（7）所示：
 |
另外，为了更好地在临床中应用，本算法可以添加约束电极数量的条件，如约束电极的使用数量为n个。
在一般情况下，对一个靶区进行刺激时，使用的电极的位置基本都位于靶区附近。而在建模流程中进行放置电极时，如果按照国际标准导联10-10系统将72个电极全部放置到头部模型上，即计算所有电极的A矩阵向量，会耗费大量的时间（大于1 d）。因此，本文为了提高电极优化的效率，选择只提取靶点附近N个电极作为电极优化的解空间SN，如式(8)所示：
 |
其中，Elec表示电极，dElec表示电极距靶点的距离，dElecN表示第N个距靶点最近的电极的距离。
1.4 评价指标
为验证本文所提出方法的性能，选取背外侧前额叶皮层（dorsolateral prefrontal cortex，DLPFC）、初级运动皮层（primary motor cortex，M1）和初级视觉皮层（primary visual cortex，V1）三个脑区作为候选刺激脑区。位置系统采用蒙特利尔神经研究所（montreal neurological institute，MNI）坐标。对应选择的靶点分别命名为靶点1（MNI: [55, 31, 31]）、靶点2（MNI: [−48, −15, 50]）和靶点3（MNI: [−22, −105, −3]）。
对于每个模型，通过计算两个指标来进行评价：① 靶区的平均电场强度；② 电场强度的聚焦性。聚焦性的衡量表征为以靶点为球心，内部累积电场强度占全部电场一半的球体的半径r0.5。电场强度越大、r0.5越小，代表优化结果越好。电场仿真通过经颅电刺激建模仿真软件ROAST（NYU，美国）实现[28]，以获取不同电极配置在头部模型上产生的电场分布。同时，在Windows 10，Intel Core i5-13490F，64 GB RAM计算机（Lenovo Inc.，中国）上进行计算，采用计算时间作为一个指标，比较不同情况下的优化时间，以说明本文所提出的方法的高效性。
2 结果
2.1 不限制候选电极位置下优化结果分布
首先，为探究不同受试者在同一刺激靶点下优化的电极空间分布，本文先提取全部候选电极的A矩阵对三个靶区进行电极优化，约束电极数量为5个，单个电极最大电流绝对值为1 mA，电流绝对值总和为4 mA。在进行优化之后发现，由于考虑了电场的聚焦性，优化得到的电极位置基本上都在靶点附近，如图2所示。本次对11个头部仿真模型进行优化，包括MNI的152标准头模（MNI152）和10个受试者头模（s1～s10）。对于靶点1，对11个头部模型进行优化的结果仅使用了靶点附近11个电极（阳极包括F4、F6，阴极包括F2、F8、FC2、FC6、FT8、C6、AFz、AF4、AF8）；对于靶点2和靶点3，对所有头部模型进行优化的结果为分别使用了靶点附近10个电极和9个电极。由此可见，本文所提优化方法得到的电极空间分布主要集中在靶区附近。在其他传统方法正常计算获取A矩阵时，通常需要对于每个电极都循环计算一次，以获得对应的A矩阵系数。然而，提取所有电极的A矩阵系数往往需要耗费大量的时间（24 h左右）。另外，从理论层面进行分析，由于电场在大脑上的分布是连续的，且刺激电流为从阳极流向阴极，在刺激电流流经的区域产生的电场强度较大，因此在进行经颅电刺激时，电场强度较大的区域一般都集中在以外部的电极位置为边界、靶点为圆心的圆周以内，因此可以考虑只选择靶点附近部分电极的A矩阵来进行电极优化。
图2
不限制候选电极位置下的电极优化结果
Figure2.
Electrode optimization results without limiting candidate electrode position
2.2 限制候选电极位置下的优化结果分布
为对比不同候选电极数目差异对优化结果的影响，对刺激靶点1只提取靶区附近10个和15个电极的A矩阵来进行优化，优化后的电场仿真结果如图3所示。图3中，E和F分别表示靶点的电场强度和电场的聚焦性。对于只选取附近10个电极进行优化，其结果为在靶区建立的电场强度可达0.31 V/m，电场聚焦性可达52.91 mm；对于只选取附近15个电极进行优化，其结果为在靶区建立的电场强度最大为0.39 V/m，最好的电场聚焦性为54.30 mm。不同数目候选电极的结果证明，只选取靶点附近10个和15个电极的A矩阵进行优化，可以实现足够的刺激强度（≥ 0.10 V/m）[15]，且电场的聚焦性也在合适范围内（≤ 70.00 mm）[26]。
图3
选取10个和15个候选电极进行优化的电场仿真结果
Figure3.
Electric field simulation results for the optimization of 10 and 15 candidate electrodes
然后，本文进一步比较了靶点附近选取10个和15个电极的A矩阵与选取全部电极的A矩阵的电极优化结果，如图4所示。对于靶点1，选取附近10个和15个电极的A矩阵进行优化，结果显示在靶区的电场强度为（0.22 ± 0.05）V/m和（0.25 ± 0.06）V/m，电场聚焦性为（57.60 ± 4.50）mm、（60.11 ± 6.44）mm。而对于全部电极的A矩阵进行优化，靶区的电场强度为（0.22 ± 0.04）V/m，电场平均聚焦性为（63.05 ± 5.66）mm。三组优化结果的电场强度和聚焦性大小都比较接近。对靶点2和靶点3进行优化也得出类似结果，即对不同电极的A矩阵进行优化之后的电场强度和聚焦性结果相似。
图4
对3个靶点在3种电极解空间（C1、C2、C3）中进行优化的结果及所用时间
Figure4.
Optimization results and time consuming of three targets in three electrodes solution spaces (C1, C2, C3)
最后，本文还统计了计算不同数量电极的A矩阵所耗费的时间，如图4所示。选取所有电极的A矩阵需要对于每个电极都循环计算一次，总共花费的时间平均为19.8 h，比较耗时。计算10个电极和15个电极的A矩阵所需要花费的平均时间分别为2.8 h和4.2 h，相比于计算所有电极的A矩阵分别节省了85.9%和78.8%，大大提高了电极优化的时间效率。因此，本文所提方法只合理计算靶点附近部分电极的A矩阵（如10个或15个）不会明显影响电极优化的结果，且可以大大节省时间成本。
3 讨论
本文提出了一种针对经颅电刺激的高效可操作的电极优化方法，在保证刺激效果的同时可节省大量时间成本，此方法通过调节权重参数来权衡电场强度和刺激聚焦性两个目标，选取靶区附近的部分电极进行优化以提升计算效率，且通过限制单个电极的电流和总电流大小保证刺激的安全性与耐受性，还可依据所使用刺激器的实际情况约束使用电极的数量，更适用于临床。
传统电极优化方法往往只针对一个目标，比如只对电场强度这一目标进行优化，这样可能会造成电场在大脑皮层分布的聚焦性太差，导致电场刺激到周围脑区，引起不可预估的后果。本文提出的电极优化方法在保证靶区电场强度足够的同时，兼顾到电场聚焦性，使得电场分布尽量集中在靶区，以保证刺激效果和安全性。另外，目前大部分的电极优化算法都不能控制使用的电极数量，这样会导致临床实际应用时存在困难，因为目前临床的电刺激设备通道数有限，一般不超过8个通道，本文提出的电极优化方法可以约束优化结果的电极数量，将其控制在8个以内，从而便于实际使用，增加可操作性。
理想情况下，使用电极优化方法提升刺激的聚焦性时，应最小化除靶区以外的所有脑区的电场强度。但在实际情况下，将所有非靶区位置的电场强度进行计算会导致计算量过大，使得算法运行时间过长甚至无法收敛。考虑到电场分布具有连续性，为了提升计算效率，本文将最小化非目标区域的电场强度表征为：首先选取头部其他位置离靶点的最远间距，获得该最远间距的四分之一距离上的坐标点合集，而后最小化这些坐标点的电场强度。通过该方法大大减少计算量，且对优化结果的影响可以忽略不计。对于人脑而言，脑部其他点距靶点的最远间距通常在200 mm左右，因此选取距靶点50 mm左右的坐标点进行电场的计算。在使用5个电极的情况下进行刺激时，聚焦性最好的结果一般为阳极放置于靶区正上方，其余四个阴极为距阳极电极最近的周围四个电极，依据以往的计算结果来看，这种电极位置情况下的电场聚焦性一般在40～50 mm，即使用5个电极的最优聚焦性在40～50 mm范围。因此本文选择抑制距靶点50 mm左右的坐标点处的场强具有合理性，有利于在提升靶区场强的同时，保持一个良好的电场聚焦性。
此外，传统电极优化方法在实际操作过程中费时费力，计算所有电极的A矩阵往往需要花费大量时间（24 h左右），极大限制了其临床应用。所以，本研究通过选取靶点附近部分电极的A矩阵向量进行优化，并将优化结果与包含全部电极的A矩阵的优化结果进行对比。本研究发现对于不同的靶点，不同三种电极数量的A矩阵的优化结果在头部模型上产生的靶区电场强度和聚焦性差别不大。此外，在进行电极优化的研究过程中，电场强度最大的位置往往在阳极和阴极之间，这个规律在文献[29-30]也有提及。在将单个电极电流限制不超过1 mA、电流绝对值之和不超过4 mA时，优化的结果中通常在靶点两边的位置包含一个阳极和一个阴极，对应电流大小分别为1 mA和—1 mA，这样就可以使得靶点处的电场强度为最大值。由此可知，要保证一定的电场聚焦性，使用的电极之间的距离一般不会太远，都会靠近靶点所在位置。也就是说，使用的电极距离靶点的位置在一定程度上可以反映出电场的聚焦性。这也是本文在进行电极优化时，发现使用到的电极几乎都分布在靶点周围的原因。所以，只获取靶点附近特定数量的电极的A矩阵进行优化，对优化结果几乎不造成影响，因为距离靶点过远的电极一般使用不到。这意味着在对某一个靶区进行电极优化时，可以仅计算靶点附近10个或15个电极的A矩阵并进行优化，而不需要获取全部电极的A矩阵，这样可以大大节省时间成本，提升效率。
除了针对经颅电刺激的电极优化的研究之外，Lynch等[31]也提出了应用于经颅磁刺激的线圈位置优化方法，即靶向功能网络刺激（targeted functional network stimulation，TANS）。TANS可以根据每位受试者的大脑网络功能拓扑特征和皮层折叠模式来优化得到一个头皮上放置的最佳线圈位置，以实现精准刺激指定功能网络靶区的同时最大程度地抑制非靶区的刺激。同TANS相比，本文提出的电极优化方法在选择靶区方面没有考虑刺激脑区的形状特征，仅简化地将刺激脑区的中心点作为靶点进行优化，这会降低干预脑区的精准度，进而可能影响到刺激的效果，TANS提供了一个优化思路，将会是未来研究的一个重点。
4 结论
本文提出了一种高效可操作的经颅电刺激电极优化方法，在多位受试者上进行应用，验证了所提方法的有效性。本文方法针对电场强度和聚焦性两个目标进行优化，可以约束电极的电流大小以及使用电极的数量，通过选取靶点附近一定数量电极进行优化，缩短了计算时间，提升了电极优化方法在临床使用中的实用性和便利性。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突
作者贡献声明：谢旭主要负责数据记录与分析、算法程序设计、论文编写；王敏敏主要负责项目主持、数据分析指导、协调沟通以及论文审阅修订；张韶岷主要负责项目主持、算法程序设计指导以及论文审阅修订。
伦理声明：本研究通过了浙江大学伦理委员会的审批（批文编号：浙大生仪会审[2022] 42号）
0 引言
经颅电刺激（transcranial electrical stimulation，TES）技术是一种非侵入式的神经调控技术，主要通过放置于头皮表面的电极施加特定的微弱电流刺激，进而调控大脑的功能活动[1]。经颅电刺激已在多种神经精神类脑疾病（比如：抑郁症[2]、癫痫[3]、帕金森病[4]和记忆障碍[5-6]）中得到应用，针对健康人群也展示出可提高运动表现[7-8]和认知功能[9-10]的潜力。
传统的经颅电刺激使用两个比较大的矩形海绵电极来进行干预，其覆盖面积一般为20～35 cm2，因而在相关脑区产生的治疗场较弥散、聚焦性差，无法进行精准的靶向刺激[11-12]。为了更加精确地控制经颅电刺激在大脑中产生的电场分布，以达到更好的刺激效果，可使用多个较小的电极（面积为1～2 cm2）代替传统的大贴片电极对脑部进行刺激[13]。高精度经颅电刺激（high definition transcranial electrical stimulation，HD-TES）的电极使用方案是一种典型的使用小电极进行刺激的应用，其包括一个阳极电极放置在靶区的正上方，四个阴极施加均匀的电流大小放在阳极四周的位置[14]，但这种方案依旧无法针对个体差异进行精准优化。为了实现个体精准靶向刺激，现已提出一些相关的经颅电刺激电极优化方法，主要是在个体经颅电刺激电场建模的基础上，结合个体脑部组织结构差异，通过电极优化算法进行计算，以获得更好的电极位置组合和对应的电流大小，使得刺激的聚焦性得到明显的改善。
如今，已经提出多种应用于经颅电刺激的电极优化算法。最早出现的是基于最小二乘法的电极优化，这种方法可以产生相对集中的刺激，即刺激的聚焦性较好，但往往在靶区产生的电场强度较低[15]。而基于线性约束最小方差法的电极优化方法通过强制设定靶区的电场强度为一个固定值，然后最小化非靶区的总能量，可以一定程度上弥补强度的不足，但容易在非目标区域形成不期望的“旁瓣”[15]。此外，有研究提出了一种最大强度法，其目标是在靶区获得尽量大的电场强度，但完全忽略了电场聚焦性的问题，导致最终优化结果在大脑产生的电场非常弥散，且电流会集中在少数的电极上，造成刺痛感[16]。还有研究提出分布约束式最大强度法，此方法可以将注入的电流分布在多个电极上，从而降低每个电极注入的电流大小[17]。另有研究使用交替方向乘数法进行电极优化，可改善目标处电场分布的聚焦性，同时降低非目标区域的强度[18]。但上述方法均为单目标优化方法，每次只能针对靶区强度或聚焦性单个目标进行优化。而平衡聚焦性和强度的刺激法（stimulation with balanced focality and intensity，SBFI）能对两个目标同时进行优化，可以得到平衡电场强度和聚焦性两个目标的优化结果，但是该方法将所有电极都置于解空间，非靶区的电场能量的计算量非常大，整个优化过程非常耗时[19]。上述方法都无法控制使用电极的数量，优化结果所使用的电极数量可多达20个，但常规的电刺激设备的通道数通常在8个以内[20]，因此这些优化方法难以投入实际应用。针对电极数量较多的问题，Zhu等[21]使用传统的多目标优化算法——非支配排序遗传算法（nondominated sorting genetic algorithm，NSGAII）进行电极优化，可以获得帕累托前沿，用于从很多电场强度和聚焦性结果相平衡的解中进行选择，此方法可以控制电极使用数量。Wang等[22]先使用遗传算法对电场强度进行单目标搜索，然后结合遗传算法与粒子群算法进行多目标优化，以获得权衡目标的帕累托前沿。Ruffini等[23]首先使用遗传算法进行计算，确定电极的位置和数量，然后使用最小二乘法分配每个电极的电流。这些经典多目标优化算法和遗传算法都具有一个很大的缺陷：优化时间过长；以遗传算法为例，通常需要耗费几个小时。综上所述，目前所提出的电极优化算法都需要通过电场建模获取候选电极的电场转移矩阵，耗时长久（> 1 d），且大部分方法都未考虑刺激器输出电流通道数，难以用于临床实践。
针对目前已有的电极优化方法存在的问题，本文提出了一种高效可操作的经颅电刺激电极优化方法，可针对目标电场强度和聚焦性两个目标进行优化，在考虑使用电极数目限制的同时，大幅提高计算效率，以期提升电极优化方法的实用性。
1 材料和方法
1.1 个体磁共振数据获取
本研究招募了10名健康成人受试者进行T1加权（T1 weighted，T1w）和T2加权（T2 weighted，T2w）磁共振成像（magnetic resonance imaging，MRI）结构扫描，以获取个体的头部结构影像学数据用于电场建模与电极优化。图像采集地点位于浙江大学华家池校区中心大楼，受试者年龄范围为20～30岁，所有受试者均签署了知情同意书。本文研究内容通过了浙江大学伦理委员会批准，并获得授权使用相关数据，伦理审查编号为：浙大生仪会审[2022] 42号。
本文所用的结构性MRI采集设备为： MAGNETOM 3T扫描仪（Siemens Inc., 德国）和32通道头部线圈。MRI的扫描参数类型包括重复时间（time of repetition，TR）、回波时间（time of echo，TE）、翻转角（flip angle）、视野范围（fields of view，FoV）和体素大小（voxel size）。T1w的具体扫描参数为：TR = 2 300 ms，TE = 2.32 ms，flip angle = 8°，FoV = 240 mm × 240 mm，voxel size = 0.9 mm3；T2w的扫描参数为：TR = 3 200 ms，TE = 408 ms，FoV = 230 mm × 230 mm，voxel size = 0.4 × 0.4 × 0.9 mm3。
1.2 个体经颅电刺激电场仿真模型构建
在进行电极优化之前，首先需要对个体人脑进行有限元建模，构建头部体积导体模型，此模型也可以被称为电场转移矩阵，即A矩阵。A矩阵表示从电极在头皮注入的电流强度到在大脑中产生的电场强度的映射关系，如式（1）所示：
|
其中，e代表不同位置的电场强度大小，s代表不同电极的电流大小。e和A的分量如式(2)所示：
|
其中，r1，···，rN代表在头部模型上的体素位置，N是头部模型中体素的数量，M为电极的总数量，e(r1)，···，e(rN)和a(r1) ，···，a(rN)分别代表r1，···，rN位置上的电场强度大小和A矩阵分量。
电场建模的目的是为了模拟外部电刺激对脑区的影响，不考虑大脑内本身存在的电信号，因此可以将模型视为无源，当通过电极将电流施加到体积的边界时，体积中产生的电势分布V可以作为拉普拉斯方程的解[24]，如式(3)所示：
|
其中，J是电流密度分布，是组织的电导率，V是电势分布。符号“”和“∙”分别表示梯度和散度。
电场建模总体流程如图1所示。首先通过扫描受试者的头部获得T1w和 T2w结构影像MRI数据，将MRI分割成几部分，包括头皮、头骨、脑脊液、灰质和白质等，对应组织电导率参数分别设置：头皮为0.465 S/m、头骨为0.010 S/m、脑脊液为1.650 S/m、白质为0.126 S/m、灰质为0.276 S/m、空气腔为2.5 × 10-14 S/m、电极为5.9 × 107 S/m、导电膏为0.300 S/m [25-27]。然后，将电极按照国际标准导联10-10系统的规定放置在头部模型上。利用有限元对头部模型进行网格化，并获得A矩阵。使用计算得到的A矩阵进行电极优化，获得可在靶区产生良好电场强度和聚焦性的电极与电流分配方案。最后，对电极优化结果进行电场仿真，观察电场在头部模型上的分布情况。
图1
电场建模总体流程
Figure1.
Overall process of electric field modeling
1.3 电极优化方法
在个体头部模型上对电场强度和聚焦性两个目标进行优化。在最大化目标靶区的电场强度的同时，最小化非目标区域的电场强度，以提升聚焦性。为减少计算量，提高计算效率，本文将最小化非目标区域的电场强度表征为：首先选取头部其他位置离靶点的最远间距dmax，获得该最远间距的四分之一距离上的坐标点合集，而后最小化这些坐标点的电场强度。总的电极优化公式如式(4)所示：
|
其中，s表示电极的电流大小，e0TCs表示目标处所需电场强度的分布，E0.25dmax表示与靶点距离为dmax的四分之一的坐标点的总电场强度。使用权重参数λ来平衡聚焦性和电场强度两个目标，λ越大，聚焦性越高；λ越小，靶区的电场强度越大。
为了保证刺激的耐受性，需要对电极优化方法施加相关安全限制：
第一，电极的输入电流大小要等于输出电流大小，即电流大小总和为0。第i个电极的电流用Si表示，如式（5）所示：
|
第二，电流的绝对值总和不大于特定值Itotal，如式（6）所示：
|
第三，单个电极的电流的绝对值不大于特定值Imax，如式（7）所示：
|
另外，为了更好地在临床中应用，本算法可以添加约束电极数量的条件，如约束电极的使用数量为n个。
在一般情况下，对一个靶区进行刺激时，使用的电极的位置基本都位于靶区附近。而在建模流程中进行放置电极时，如果按照国际标准导联10-10系统将72个电极全部放置到头部模型上，即计算所有电极的A矩阵向量，会耗费大量的时间（大于1 d）。因此，本文为了提高电极优化的效率，选择只提取靶点附近N个电极作为电极优化的解空间SN，如式(8)所示：
|
其中，Elec表示电极，dElec表示电极距靶点的距离，dElecN表示第N个距靶点最近的电极的距离。
1.4 评价指标
为验证本文所提出方法的性能，选取背外侧前额叶皮层（dorsolateral prefrontal cortex，DLPFC）、初级运动皮层（primary motor cortex，M1）和初级视觉皮层（primary visual cortex，V1）三个脑区作为候选刺激脑区。位置系统采用蒙特利尔神经研究所（montreal neurological institute，MNI）坐标。对应选择的靶点分别命名为靶点1（MNI: [55, 31, 31]）、靶点2（MNI: [−48, −15, 50]）和靶点3（MNI: [−22, −105, −3]）。
对于每个模型，通过计算两个指标来进行评价：① 靶区的平均电场强度；② 电场强度的聚焦性。聚焦性的衡量表征为以靶点为球心，内部累积电场强度占全部电场一半的球体的半径r0.5。电场强度越大、r0.5越小，代表优化结果越好。电场仿真通过经颅电刺激建模仿真软件ROAST（NYU，美国）实现[28]，以获取不同电极配置在头部模型上产生的电场分布。同时，在Windows 10，Intel Core i5-13490F，64 GB RAM计算机（Lenovo Inc.，中国）上进行计算，采用计算时间作为一个指标，比较不同情况下的优化时间，以说明本文所提出的方法的高效性。
2 结果
2.1 不限制候选电极位置下优化结果分布
首先，为探究不同受试者在同一刺激靶点下优化的电极空间分布，本文先提取全部候选电极的A矩阵对三个靶区进行电极优化，约束电极数量为5个，单个电极最大电流绝对值为1 mA，电流绝对值总和为4 mA。在进行优化之后发现，由于考虑了电场的聚焦性，优化得到的电极位置基本上都在靶点附近，如图2所示。本次对11个头部仿真模型进行优化，包括MNI的152标准头模（MNI152）和10个受试者头模（s1～s10）。对于靶点1，对11个头部模型进行优化的结果仅使用了靶点附近11个电极（阳极包括F4、F6，阴极包括F2、F8、FC2、FC6、FT8、C6、AFz、AF4、AF8）；对于靶点2和靶点3，对所有头部模型进行优化的结果为分别使用了靶点附近10个电极和9个电极。由此可见，本文所提优化方法得到的电极空间分布主要集中在靶区附近。在其他传统方法正常计算获取A矩阵时，通常需要对于每个电极都循环计算一次，以获得对应的A矩阵系数。然而，提取所有电极的A矩阵系数往往需要耗费大量的时间（24 h左右）。另外，从理论层面进行分析，由于电场在大脑上的分布是连续的，且刺激电流为从阳极流向阴极，在刺激电流流经的区域产生的电场强度较大，因此在进行经颅电刺激时，电场强度较大的区域一般都集中在以外部的电极位置为边界、靶点为圆心的圆周以内，因此可以考虑只选择靶点附近部分电极的A矩阵来进行电极优化。
图2
不限制候选电极位置下的电极优化结果
Figure2.
Electrode optimization results without limiting candidate electrode position
2.2 限制候选电极位置下的优化结果分布
为对比不同候选电极数目差异对优化结果的影响，对刺激靶点1只提取靶区附近10个和15个电极的A矩阵来进行优化，优化后的电场仿真结果如图3所示。图3中，E和F分别表示靶点的电场强度和电场的聚焦性。对于只选取附近10个电极进行优化，其结果为在靶区建立的电场强度可达0.31 V/m，电场聚焦性可达52.91 mm；对于只选取附近15个电极进行优化，其结果为在靶区建立的电场强度最大为0.39 V/m，最好的电场聚焦性为54.30 mm。不同数目候选电极的结果证明，只选取靶点附近10个和15个电极的A矩阵进行优化，可以实现足够的刺激强度（≥ 0.10 V/m）[15]，且电场的聚焦性也在合适范围内（≤ 70.00 mm）[26]。
图3
选取10个和15个候选电极进行优化的电场仿真结果
Figure3.
Electric field simulation results for the optimization of 10 and 15 candidate electrodes
然后，本文进一步比较了靶点附近选取10个和15个电极的A矩阵与选取全部电极的A矩阵的电极优化结果，如图4所示。对于靶点1，选取附近10个和15个电极的A矩阵进行优化，结果显示在靶区的电场强度为（0.22 ± 0.05）V/m和（0.25 ± 0.06）V/m，电场聚焦性为（57.60 ± 4.50）mm、（60.11 ± 6.44）mm。而对于全部电极的A矩阵进行优化，靶区的电场强度为（0.22 ± 0.04）V/m，电场平均聚焦性为（63.05 ± 5.66）mm。三组优化结果的电场强度和聚焦性大小都比较接近。对靶点2和靶点3进行优化也得出类似结果，即对不同电极的A矩阵进行优化之后的电场强度和聚焦性结果相似。
图4
对3个靶点在3种电极解空间（C1、C2、C3）中进行优化的结果及所用时间
Figure4.
Optimization results and time consuming of three targets in three electrodes solution spaces (C1, C2, C3)
最后，本文还统计了计算不同数量电极的A矩阵所耗费的时间，如图4所示。选取所有电极的A矩阵需要对于每个电极都循环计算一次，总共花费的时间平均为19.8 h，比较耗时。计算10个电极和15个电极的A矩阵所需要花费的平均时间分别为2.8 h和4.2 h，相比于计算所有电极的A矩阵分别节省了85.9%和78.8%，大大提高了电极优化的时间效率。因此，本文所提方法只合理计算靶点附近部分电极的A矩阵（如10个或15个）不会明显影响电极优化的结果，且可以大大节省时间成本。
3 讨论
本文提出了一种针对经颅电刺激的高效可操作的电极优化方法，在保证刺激效果的同时可节省大量时间成本，此方法通过调节权重参数来权衡电场强度和刺激聚焦性两个目标，选取靶区附近的部分电极进行优化以提升计算效率，且通过限制单个电极的电流和总电流大小保证刺激的安全性与耐受性，还可依据所使用刺激器的实际情况约束使用电极的数量，更适用于临床。
传统电极优化方法往往只针对一个目标，比如只对电场强度这一目标进行优化，这样可能会造成电场在大脑皮层分布的聚焦性太差，导致电场刺激到周围脑区，引起不可预估的后果。本文提出的电极优化方法在保证靶区电场强度足够的同时，兼顾到电场聚焦性，使得电场分布尽量集中在靶区，以保证刺激效果和安全性。另外，目前大部分的电极优化算法都不能控制使用的电极数量，这样会导致临床实际应用时存在困难，因为目前临床的电刺激设备通道数有限，一般不超过8个通道，本文提出的电极优化方法可以约束优化结果的电极数量，将其控制在8个以内，从而便于实际使用，增加可操作性。
理想情况下，使用电极优化方法提升刺激的聚焦性时，应最小化除靶区以外的所有脑区的电场强度。但在实际情况下，将所有非靶区位置的电场强度进行计算会导致计算量过大，使得算法运行时间过长甚至无法收敛。考虑到电场分布具有连续性，为了提升计算效率，本文将最小化非目标区域的电场强度表征为：首先选取头部其他位置离靶点的最远间距，获得该最远间距的四分之一距离上的坐标点合集，而后最小化这些坐标点的电场强度。通过该方法大大减少计算量，且对优化结果的影响可以忽略不计。对于人脑而言，脑部其他点距靶点的最远间距通常在200 mm左右，因此选取距靶点50 mm左右的坐标点进行电场的计算。在使用5个电极的情况下进行刺激时，聚焦性最好的结果一般为阳极放置于靶区正上方，其余四个阴极为距阳极电极最近的周围四个电极，依据以往的计算结果来看，这种电极位置情况下的电场聚焦性一般在40～50 mm，即使用5个电极的最优聚焦性在40～50 mm范围。因此本文选择抑制距靶点50 mm左右的坐标点处的场强具有合理性，有利于在提升靶区场强的同时，保持一个良好的电场聚焦性。
此外，传统电极优化方法在实际操作过程中费时费力，计算所有电极的A矩阵往往需要花费大量时间（24 h左右），极大限制了其临床应用。所以，本研究通过选取靶点附近部分电极的A矩阵向量进行优化，并将优化结果与包含全部电极的A矩阵的优化结果进行对比。本研究发现对于不同的靶点，不同三种电极数量的A矩阵的优化结果在头部模型上产生的靶区电场强度和聚焦性差别不大。此外，在进行电极优化的研究过程中，电场强度最大的位置往往在阳极和阴极之间，这个规律在文献[29-30]也有提及。在将单个电极电流限制不超过1 mA、电流绝对值之和不超过4 mA时，优化的结果中通常在靶点两边的位置包含一个阳极和一个阴极，对应电流大小分别为1 mA和—1 mA，这样就可以使得靶点处的电场强度为最大值。由此可知，要保证一定的电场聚焦性，使用的电极之间的距离一般不会太远，都会靠近靶点所在位置。也就是说，使用的电极距离靶点的位置在一定程度上可以反映出电场的聚焦性。这也是本文在进行电极优化时，发现使用到的电极几乎都分布在靶点周围的原因。所以，只获取靶点附近特定数量的电极的A矩阵进行优化，对优化结果几乎不造成影响，因为距离靶点过远的电极一般使用不到。这意味着在对某一个靶区进行电极优化时，可以仅计算靶点附近10个或15个电极的A矩阵并进行优化，而不需要获取全部电极的A矩阵，这样可以大大节省时间成本，提升效率。
除了针对经颅电刺激的电极优化的研究之外，Lynch等[31]也提出了应用于经颅磁刺激的线圈位置优化方法，即靶向功能网络刺激（targeted functional network stimulation，TANS）。TANS可以根据每位受试者的大脑网络功能拓扑特征和皮层折叠模式来优化得到一个头皮上放置的最佳线圈位置，以实现精准刺激指定功能网络靶区的同时最大程度地抑制非靶区的刺激。同TANS相比，本文提出的电极优化方法在选择靶区方面没有考虑刺激脑区的形状特征，仅简化地将刺激脑区的中心点作为靶点进行优化，这会降低干预脑区的精准度，进而可能影响到刺激的效果，TANS提供了一个优化思路，将会是未来研究的一个重点。
4 结论
本文提出了一种高效可操作的经颅电刺激电极优化方法，在多位受试者上进行应用，验证了所提方法的有效性。本文方法针对电场强度和聚焦性两个目标进行优化，可以约束电极的电流大小以及使用电极的数量，通过选取靶点附近一定数量电极进行优化，缩短了计算时间，提升了电极优化方法在临床使用中的实用性和便利性。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突
作者贡献声明：谢旭主要负责数据记录与分析、算法程序设计、论文编写；王敏敏主要负责项目主持、数据分析指导、协调沟通以及论文审阅修订；张韶岷主要负责项目主持、算法程序设计指导以及论文审阅修订。
伦理声明：本研究通过了浙江大学伦理委员会的审批（批文编号：浙大生仪会审[2022] 42号）
