經顱電刺激是一種極具潛力的無創神經調控技術。基于個體經顱電刺激電場仿真模型的電極優化方法可依據個體頭部組織結構差異給出個性化電刺激參數，顯著提升經顱電刺激的干預靶向性。然而，目前所提出的電極優化方法存在計算耗時較長（通常大于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號）
