精準定位致癇灶是手術切除致癇灶的前提，而傳統基于三維球模型或標準頭模型模板的溯源結果與真實致癇灶誤差較大。本研究擬在構建真實頭模型的基礎上，利用癲癇患者睡眠期腦電中的癇性發作棘波，通過多偶極子與真實頭模型計算出大腦皮層電流密度分布，對致癇灶進行溯源定位。實驗結果表明，基于該算法的溯源結果定位準確率為89.27%，可使植入電極的數量減少（19.34 ± 7.15）%。該方法不僅能夠提高致癇灶定位準確性，還能夠減少術中額外傷害，為神經外科醫生制定手術方案提供了更加直觀有效的參考。
引用本文： 屈若為, 王召楠, 王石峰, 王瑤, 王樂, 尹紹雅, 顧軍華, 徐桂芝. 基于真實頭模型與多偶極子算法的癲癇致癇灶腦電溯源方法研究. 生物醫學工程學雜志, 2023, 40(2): 272-279. doi: 10.7507/1001-5515.202209045 復制
0 引言
癲癇是一種具有突發性和易復發性的腦功能障礙疾病，已經成為神經科僅次于頭痛的第二大疾病。目前，世界范圍內約有5 000萬癲癇患者[1]，我國約有900萬癲癇患者，平均每年約有40萬新發病例[2]。癲癇發作時，患者肢體不受自主意識控制，即使短暫的發作也會給患者帶來生理和心理上的損傷，嚴重的癲癇發作會導致腦細胞死亡，甚至危及生命。因此，癲癇的治療無論對患者個體和家庭的生活質量，還是對我國國民的整體幸福指數，都具有重大意義。
盡管已有多種抗癲癇藥物，但僅有約67%患者的癲癇發作能夠得到有效控制[3-4]，8%～10%的患者需要接受手術或神經調控技術進行治療[5]。對于藥物難治性癲癇，通過手術方式去除致癇灶（epileptogenic zone，EZ）或神經調控是根治癲癇的重要方法。由于手術和神經調控均屬于侵入式治療手段，有潛在的感染與預后風險[6]，因此，準確的致癇灶定位能夠減少顱內電極的植入數量，降低顱內感染風險[7]。
目前用于臨床的癲癇病灶定位方法主要是基于腦電與影像數據。核磁共振成像（magnetic resonance image，MRI）[8]和電子計算機斷層掃描（compted tomography，CT）[9]可以通過大腦的靜態解剖結構圖來診斷大腦中的病變位置，從而判斷致癇灶[10]。功能MRI和正電子發射計算機斷層掃描（positron emission tomography，PET）則通過大腦功能活動時的代謝活動圖像來獲得大腦的動態圖像信息[11]。在這些檢查方法中，PET和功能MRI由于技術要求高和費用昂貴而不能作為普遍的檢測技術手段，結構MRI與CT則不能判斷非結構性的病變。頭皮腦電（electroencephalograph，EEG）具有時間分辨率高、檢查價格低廉以及非侵入式檢查等優點，在越來越多的理論研究和臨床實踐中得到廣泛的關注。因此，通過EEG建立逆問題方程對大腦皮層致癇灶進行溯源定位，在癲癇的術前評估中具有廣闊的前景和重要的價值[12]。
在通過EEG對致癇灶的溯源定位問題中，選取恰當的頭模型與建立適當的腦電傳輸模型是溯源結果準確與否的關鍵。傳統的三層同心球模型中，頭被看作是三個同心球，分為頭皮、顱骨和大腦皮層三層，因此又被稱為SSB（scalp-skull-brain）頭模型[13]。癲癇腦電源定位問題是從給定的頭皮或皮層腦電分布求解出異常放電發作起始區定位的問題，通常也稱為癲癇腦電逆問題。在其他腦神經疾病如阿爾茨海默病等全面退行性疾病中，多用標準低分辨率電磁斷層成像方法（standard low-resolution electromagnetic tomography，sLORETA）[14]對大腦皮層電流密度分布進行溯源，該算法是最小范數估計算法（minimum norm estimation，MNE）的標準化版本[15]。
盡管有很多研究利用腦電圖技術進行癲癇溯源，但這些研究仍然存在一些局限性：① 溯源使用的頭模型不是真實頭模型。癲癇大腦通常伴隨著結構或功能上的改變，利用理想球模型或標準頭模型解決溯源問題，雖然能夠在一定程度上模擬大腦的結構，但并不能細致地解釋大腦的真實物理結構和物理特性，會影響結果的可信度。② 大部分研究利用分布式源將EEG溯源到大腦皮層[16]，但分布式源會擴大溯源范圍，不適于癲癇致癇灶定位問題的研究。
因此，本文根據腦組織不同的導電特性構建真實頭模型，并基于邊界元算法剖分；采用癲癇患者EEG中的發作間期棘波，基于多偶極子模型對癲癇致癇灶進行溯源定位；基于上述溯源后的大腦皮層電流密度分布，構建不同腦區之間的功能連接網絡，分析致癇區與受累區的網絡局部特性。本研究旨在提高致癇灶定位準確性，減少術前檢查和手術操作中帶來的額外傷害和潛在風險，為神經外科醫生制定手術方案提供更加直觀、有效的參考。
1 實驗數據收集與預處理
1.1 臨床數據收集
本實驗回顧性且不間斷地收集了天津市環湖醫院2022年1月至8月共12例癲癇患者。所有患者術前均進行了EEG、MRI的T1序列和PET檢查，隨后進行了電極植入手術。排除標準為：① 嚴重腦萎縮導致無法重建大腦皮層；② 全面性癲癇，非局灶性癲癇；③ 電極植入期間未捕捉到癲癇樣放電，無法獲得病理金標準。最終共納入8例癲癇患者，其中男性患者5例，女性患者3例，年齡（29.1 ± 13.6）歲。致癇灶的診斷由神經外科醫生和影像科醫生根據國際癲癇聯盟2017年提出的診斷標準進行定位標注[17]。所有參與者均簽署了知情同意書。
1.2 真實頭模型的構建與剖分
1.2.1 核磁共振T1加權圖像的收集
設備及參數：癲癇患者腦部核磁共振T1加權影像數據使用3.0-Ingenia MR掃描儀（Philips Healthcare，荷蘭）進行11 min掃描，turbo-pin回波序列的重復時間、反轉時間和回波時間分別為6 600 ms、2 000 ms和330 ms。采集參數如下：翻轉角度：180°；渦旋回波因子：182；帶寬：1 187.1 Hz；場：220 mm×180 mm×150 mm；矩陣：220 mm×180 mm×150 mm；圖像厚度：1 mm；激發次數：2次；切片數：150片；CS-SENSE：3。
構建真實頭模型：采集到的原始數據為DICOM格式的腦軸位圖，通過Matlab中的SPM工具箱轉換為NIFTI格式。借助Matlab中的CAT工具箱分割頭皮、顱骨、皮質、白質和腦脊液，并根據解剖結構分區重建。將得到的被試圖像線性配準到標準模板空間上，統一所有圖像的坐標空間，然后基于標準模板對圖像進行坐標變換和仿射變換。在完成上述處理后，對圖像進行平滑處理，以抑制噪聲，提高信噪比。
1.2.2 基于對稱邊界元法的真實頭模型剖分與構建
常用的基于醫學成像技術構建頭部模型的方法主要包括：邊界元法、有限元法和有限差分法。由于邊界元法是基于給定的電流源計算導體表面或者邊界電位，只需對邊界離散，可以降低邊界元的自由度，降低解的空間維數，因此本研究正向模型采用對稱邊界元方法剖分真實頭模型，并設計三個表面：大腦組織和顱骨之間的接觸面、顱骨和頭皮的接觸面以及頭皮和空氣的接觸面。這些接觸面光滑且通向傳導，每一個表面被網格化為小的三角形元素。均勻邊界元導體V上任意一點r上對應的電位V(r)定義如下：
 '/> |
其中 
為電導率，
為點r上的電勢，
和 
分別表示接觸面 
內部和外部的導電率，
是與接觸面元素正交的向量，
表示接觸面元素的面積。每一個接觸面 
被分割為N個三角形。每個三角形的電勢通過上式計算可得。接觸面 
的電勢被定義為：
 '/> |
利用分割后和網格化后的數據構建頭部的真實物理模型。設置電導率：顱骨0.013 2 S/m，頭皮0.35 S/m，腦脊液1.79 S/m，皮層0.33 S/m，白質0.14 S/m。
1.3 腦電數據采集與預處理
本課題采用國際10-20系統采集癲癇患者睡眠期腦電，記錄了19個通道的腦電信號。根據受試者頭部電極的實際位置對部分電極點的位置進行了微調和校準。對采集到的EEG進行0.5～80 Hz的帶通濾波，并去除了50 Hz的工頻干擾。由于腦電源成像算法對噪聲敏感，因此采用獨立成分分析去除眼動偽跡和心電偽跡[18]。癲癇患者睡眠期EEG中與癇性發作相關的棘波由神經外科醫生手動標記，在每個被試數據中截取9～14個從棘波波峰前100 ms到棘波波峰后300 ms的EEG片段，并對這些片段進行疊加平均。
此外，在腦電誘發響應的逆問題研究中，由于作為背景噪聲的自發腦電的存在，EEG信噪比非常低。為了提高信噪比，本實驗利用特征值中位數進行噪聲協方差正規化，選取連續EEG記錄中沒有癲癇活動、快速眼動及其他偽影的50 s計算協方差矩陣，用以去除經預處理后依然存在的由儀器、傳感器和環境等帶來的噪聲。
2 原理與方法
2.1 基于多偶極子的腦電溯源算法
腦電逆問題主要分為基于空間約束的源成像和基于時空約束的源成像[19]。目前，由于分布式源成像模型具有一定的彌散性，能夠解決大部分腦彌漫性病變的溯源問題，如阿爾茨海默病[20]，但局灶性癲癇的致癇范圍相對比較局限，分布式源成像模型會帶來較高的定位假陽性。不同于其他腦電逆問題模型，基于偶極子模型的逆問題模型能夠更加精準地反演出局灶性神經元活動的位置，在致癇灶的定位問題中已有應用[21]。因此本文選擇通過計算復雜度低且定位范圍局限的偶極子模型定位大腦皮層中異常放電位置。
基于多偶極子逆問題算法的主要思想是通過調整偶極子的各個參數，擬合正向計算結果與頭皮電極測量數據。由于傳遞函數L為非線性函數，因此在未知源可用有限時，逆問題可表述為尋找一組源參數使得目標函數最小，其數學表達為 
。本文中的最小目標函數通過最小二乘法選取：
 |
其中，V為對應給定(r, q)的正向計算電位值，
為電位實測值。偶極子位置及擬合優度如圖1所示。偶極子可以根據時間著色（見圖1a），也可以根據擬合優度著色（見圖1b），值越大，代表擬合程度越高。將多組偶極子顯示在同一個MRI查看器上，可以將偶極子擬合到一系列單獨的癲癇尖峰上，觀察偶極子是否形成一個簇。如圖1c所示，將優度值低于所選閾值的偶極子予以隱藏，僅顯示具有最佳擬合度的偶極子。
圖1
3D MRI及真實頭模型中的多個偶極子的位置及擬合優度
a. 時間著色偶極子；b. 擬合優度著色偶極子；c. 最佳擬合度偶極子
Figure1. The position and value of multi dipoles in 3D MRI and realistic head modela. time coloring dipoles; b. fitting coloring dipoles; c. the best fitting dipole
2.2 基于相位傳遞熵的大腦皮層網絡構建
相位傳遞熵（phase transfer entropy，PTE）[22]是一種不依賴先驗模型的因果關系度量工具，適用于各類非平穩和高維隨機過程。對于給定的時間滯后
，從信號X到信號Y的相位傳遞熵定義為：
 '/> |
其中，
、
及 
分別為響應時間序列的概率分布。
和 
分別是時間點 
信號X與信號Y 的相位：
，
。本文中設定 
= 1。
在構建腦網絡的問題中，該方法可以分析任意多個通道或大腦區域之間的時變因果關系，不僅可以顯示給定周期內的信號信息，還可以識別兩個系統的同步性，具有良好的統計顯著性[23-24]。因此，本文基于相位傳遞熵構建了不同腦區之間的功能連接網絡，通過信息的因果流向關系區分致癇區與受累區，進一步提高致癇灶定位精確度。
3 結果與分析
3.1 基于致癇灶溯源的皮層電流密度分布結果
本文回顧性地分析了8例局灶性癲癇患者基于多偶極子算法與真實頭模型的癲癇致癇灶腦電溯源結果。實驗的金標準為立體定位腦電圖（stereoelectroencephalogram，SEEG）監測到癇性發作的起始位置或致癇灶去除手術的預后效果。使用三種經典的逆向解算法——動態統計參數映射（dynamic statistical parameter mapping，dSPM）、sLORETA和分布式源法與本文提出的多偶極子算法進行比較。同時，在每種算法中對比了無約束條件和有額外約束條件的逆向解。此處以01號患者為例，患者性別男，39歲，有家族史，癲癇發作史2～3年，藥物控制效果差。依據各算法得到的該患者皮層電流密度結果如圖2所示。
圖2
通過不同逆問題算法得到的01號癲癇患者致癇灶定位結果
Figure2.
EZ localization results of patient No. 01 based on various algorithms
影像診斷提示右側杏仁核異常且不排除左側異常，臨床醫生通過癥狀學及腦電圖診斷為右側顳葉癲癇。隨后，該患者接受了電極植入手術，植入電極數目為7根，每根電極上分別有10個電極點，如圖3所示。
圖3
手術植入電極位置
Figure3.
The position of implanted SEEG
經過SEEG電極檢測，癇性發作起始期的棘波由電極A上的第二個電極點A2檢測到，即右側杏仁核位置。檢測到的棘波由A、B、C、D擴散到F電極，即異常放電由右側顳葉擴散到右側前葉、島葉。從圖2與圖3的對比分析中可以看出，四種算法對致癇灶的定側結果準確，但定位結果各不相同。其中，dSPM算法與sLORETA算法定位范圍較大，無法對植入電極或手術切除位置提供指導。電流分布密度算法雖定位較為局限，但定位在顳葉外側，與真實致癇灶位置不符。本文提出的多偶極子逆問題算法的溯源結果為顳葉內側杏仁核，在定位準確性和范圍局限性方面均優于其他三種算法。
此外，由于術前腦電及影像不能排除左側顳葉異常放電的可能，因此手術植入電極A’檢測該處是否為致癇區。如圖4所示，依據本實驗得到的致癇灶定位結果，將截取的9個400 ms棘波片段疊加平均，定義癇性發作棘波波峰出現的時刻為0.0 s，將此時的電流密度分布圖調整為只顯示大于12 μA的皮層電流密度區域，可以看到杏仁核與海馬放電異常。若將電流密度分布的顯示閾值調整為10 μA，可以看到大部分右側顳葉及小部分左側顳葉及杏仁核位置有異常電流。植入的電極A’覆蓋了左顳區域。接著，將時間軸移動到棘波波峰出現后的100 ms時刻，此時電流密度分布圖顯示左側顳葉、額葉及枕葉均有相對較高的異常信號電流，電流密度值在7～8 μA之間，低于異常放電起始時的電流密度值。從不同時刻大腦皮層的電流密度分布差異可以推斷，異常放電的起始區位于右側杏仁核及海馬，而右顳外側及左側顳葉、額葉、枕葉均為功能受累區。SEEG結果也證實了這一推斷。
圖4
01號患者癇性發作棘波出現時刻的大腦皮層電流密度分布圖
Figure4.
The current density distribution map on cortex during seizure spike of patient No. 01
3.2 癲癇患者腦功能網絡分析
通過以上實驗得到癲癇患者大腦皮層的電流密度分布，根據Desikan-Killiany地圖集將大腦皮層分為68個感興趣區域。在構建腦網絡的問題中，相位傳遞熵可以分析任意多個通道或大腦區域之間的時變因果關系，具有良好的統計顯著性和效力。因此，本文基于相位傳遞熵構建了不同腦區之間的功能連接網絡，如圖5所示。根據腦網絡圖，大腦皮層信息由右側顳橫回流向左側額下島蓋部，由右側頂下小葉流向左側顳葉內嗅皮層，并整體呈現出皮層信息由大腦右側半球流向左側的趨勢。在左側大腦半球功能區之間，信息由左側額極流向顳葉梭形回。因此可以推斷：致癇區在大腦右側，左側為受累區；致癇區在右側顳葉，左側額葉、顳葉為受累區，這與前文分析的致癇灶定側與定位結果表現出了高度的一致性。因此，若術前參考本次實驗結果，確定左側顳葉為癇性發作的功能受累區，而非起始發作區，則可不植入圖3中所示的電極A’。
圖5
基于相位傳遞熵構建的顳葉癲癇患者腦網絡
Figure5.
The brain network of temporal lobe epilepsy based on PTE
3.3 致癇灶源定位結果與分析
與正常人大腦不同，反復的癲癇發作會使癲癇患者的大腦發生結構或功能上的改變。基于癲癇患者結構核磁圖像重建的真實頭模型，不僅能夠在逆問題方程的計算中得到更加精確的結果，還能為手術方案的制定提供更加直觀的展示[25]。8位患者基于真實頭模型、多偶極子與相位傳遞熵網絡分析的致癇灶定位結果如圖6所示。
圖6
基于多偶極子的致癇灶溯源至大腦皮層的電流密度分布圖
Figure6.
The current density distribution map of EZ localization on cortex based on multi-poles
將基于4種算法得到的放電起始位置與SEEG檢測到的起始放電位置作比較，得到的定位準確性、敏感性、特異性和F1指標如表1所示。癲癇患者的溯源結果與臨床結果對比分析如表2所示。實驗結果表明，基于多偶極子溯源與相位傳遞熵網絡分析定位致癇灶結果的各項指標均優于其他三種算法。值得一提的是，多偶極子算法結果的特異性較高，即錯診率低，如果術前參考該腦電溯源結果，就能在癲癇手術中盡可能地避免損傷正常腦組織。同時，在本文的回顧性分析中，將溯源結果、臨床及影像診斷與SEEG金標準作對比，依據本實驗溯源結果制定的SEEG植入電極數目比原方案減少（19.34 ± 7.15）%，降低了因植入電極帶來的感染風險及手術費用。
表1
基于4種算法定位致癇灶的結果
Table1.
Localization result of EZ based on 4 algorithms
表2
癲癇患者的溯源結果與臨床結果對比分析
Table2.
Comparison and analysis between localization based on multi-poles and clinical results of patients
4 結論
本研究旨在通過生物電磁成像算法，實現癲癇致癇灶的無創、精準定位。本文基于癲癇患者的多模態醫學數據，通過T1 MRI構建癲癇患者的真實頭模型，引入多偶極子算法定位致癇灶，構建了大腦皮層網絡進一步區分致癇灶與其受累區。實驗結果表明，針對癲癇致癇灶定位問題，本文提出的基于真實頭模型與多偶極子算法的定位結果精確、范圍局限。該項工作不僅能夠提高致癇灶定位的準確性，還能減少術前檢查和手術操作中帶來的額外傷害和潛在風險，為神經外科醫生制定手術方案提供更加直觀、有效的參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：屈若為是本研究實驗設計、實驗研究、數據分析和論文寫作的執行人；王召楠負責實驗數據收集和預處理；王石峰負責醫學影像數據的收集及標注；王瑤負責腦網絡的構建與該部分知識的審核；王樂參與腦電實驗設計、收集癲癇患者腦電數據并標注；尹紹雅、顧軍華和徐桂芝教授是基金的負責人，對文章知識性內容作批評性審閱與指導，參與論文校對。
倫理聲明：本研究通過了河北工業大學生物醫學倫理委員會的審批（批文編號：HEBUThMEC2022003）。
0 引言
癲癇是一種具有突發性和易復發性的腦功能障礙疾病，已經成為神經科僅次于頭痛的第二大疾病。目前，世界范圍內約有5 000萬癲癇患者[1]，我國約有900萬癲癇患者，平均每年約有40萬新發病例[2]。癲癇發作時，患者肢體不受自主意識控制，即使短暫的發作也會給患者帶來生理和心理上的損傷，嚴重的癲癇發作會導致腦細胞死亡，甚至危及生命。因此，癲癇的治療無論對患者個體和家庭的生活質量，還是對我國國民的整體幸福指數，都具有重大意義。
盡管已有多種抗癲癇藥物，但僅有約67%患者的癲癇發作能夠得到有效控制[3-4]，8%～10%的患者需要接受手術或神經調控技術進行治療[5]。對于藥物難治性癲癇，通過手術方式去除致癇灶（epileptogenic zone，EZ）或神經調控是根治癲癇的重要方法。由于手術和神經調控均屬于侵入式治療手段，有潛在的感染與預后風險[6]，因此，準確的致癇灶定位能夠減少顱內電極的植入數量，降低顱內感染風險[7]。
目前用于臨床的癲癇病灶定位方法主要是基于腦電與影像數據。核磁共振成像（magnetic resonance image，MRI）[8]和電子計算機斷層掃描（compted tomography，CT）[9]可以通過大腦的靜態解剖結構圖來診斷大腦中的病變位置，從而判斷致癇灶[10]。功能MRI和正電子發射計算機斷層掃描（positron emission tomography，PET）則通過大腦功能活動時的代謝活動圖像來獲得大腦的動態圖像信息[11]。在這些檢查方法中，PET和功能MRI由于技術要求高和費用昂貴而不能作為普遍的檢測技術手段，結構MRI與CT則不能判斷非結構性的病變。頭皮腦電（electroencephalograph，EEG）具有時間分辨率高、檢查價格低廉以及非侵入式檢查等優點，在越來越多的理論研究和臨床實踐中得到廣泛的關注。因此，通過EEG建立逆問題方程對大腦皮層致癇灶進行溯源定位，在癲癇的術前評估中具有廣闊的前景和重要的價值[12]。
在通過EEG對致癇灶的溯源定位問題中，選取恰當的頭模型與建立適當的腦電傳輸模型是溯源結果準確與否的關鍵。傳統的三層同心球模型中，頭被看作是三個同心球，分為頭皮、顱骨和大腦皮層三層，因此又被稱為SSB（scalp-skull-brain）頭模型[13]。癲癇腦電源定位問題是從給定的頭皮或皮層腦電分布求解出異常放電發作起始區定位的問題，通常也稱為癲癇腦電逆問題。在其他腦神經疾病如阿爾茨海默病等全面退行性疾病中，多用標準低分辨率電磁斷層成像方法（standard low-resolution electromagnetic tomography，sLORETA）[14]對大腦皮層電流密度分布進行溯源，該算法是最小范數估計算法（minimum norm estimation，MNE）的標準化版本[15]。
盡管有很多研究利用腦電圖技術進行癲癇溯源，但這些研究仍然存在一些局限性：① 溯源使用的頭模型不是真實頭模型。癲癇大腦通常伴隨著結構或功能上的改變，利用理想球模型或標準頭模型解決溯源問題，雖然能夠在一定程度上模擬大腦的結構，但并不能細致地解釋大腦的真實物理結構和物理特性，會影響結果的可信度。② 大部分研究利用分布式源將EEG溯源到大腦皮層[16]，但分布式源會擴大溯源范圍，不適于癲癇致癇灶定位問題的研究。
因此，本文根據腦組織不同的導電特性構建真實頭模型，并基于邊界元算法剖分；采用癲癇患者EEG中的發作間期棘波，基于多偶極子模型對癲癇致癇灶進行溯源定位；基于上述溯源后的大腦皮層電流密度分布，構建不同腦區之間的功能連接網絡，分析致癇區與受累區的網絡局部特性。本研究旨在提高致癇灶定位準確性，減少術前檢查和手術操作中帶來的額外傷害和潛在風險，為神經外科醫生制定手術方案提供更加直觀、有效的參考。
1 實驗數據收集與預處理
1.1 臨床數據收集
本實驗回顧性且不間斷地收集了天津市環湖醫院2022年1月至8月共12例癲癇患者。所有患者術前均進行了EEG、MRI的T1序列和PET檢查，隨后進行了電極植入手術。排除標準為：① 嚴重腦萎縮導致無法重建大腦皮層；② 全面性癲癇，非局灶性癲癇；③ 電極植入期間未捕捉到癲癇樣放電，無法獲得病理金標準。最終共納入8例癲癇患者，其中男性患者5例，女性患者3例，年齡（29.1 ± 13.6）歲。致癇灶的診斷由神經外科醫生和影像科醫生根據國際癲癇聯盟2017年提出的診斷標準進行定位標注[17]。所有參與者均簽署了知情同意書。
1.2 真實頭模型的構建與剖分
1.2.1 核磁共振T1加權圖像的收集
設備及參數：癲癇患者腦部核磁共振T1加權影像數據使用3.0-Ingenia MR掃描儀（Philips Healthcare，荷蘭）進行11 min掃描，turbo-pin回波序列的重復時間、反轉時間和回波時間分別為6 600 ms、2 000 ms和330 ms。采集參數如下：翻轉角度：180°；渦旋回波因子：182；帶寬：1 187.1 Hz；場：220 mm×180 mm×150 mm；矩陣：220 mm×180 mm×150 mm；圖像厚度：1 mm；激發次數：2次；切片數：150片；CS-SENSE：3。
構建真實頭模型：采集到的原始數據為DICOM格式的腦軸位圖，通過Matlab中的SPM工具箱轉換為NIFTI格式。借助Matlab中的CAT工具箱分割頭皮、顱骨、皮質、白質和腦脊液，并根據解剖結構分區重建。將得到的被試圖像線性配準到標準模板空間上，統一所有圖像的坐標空間，然后基于標準模板對圖像進行坐標變換和仿射變換。在完成上述處理后，對圖像進行平滑處理，以抑制噪聲，提高信噪比。
1.2.2 基于對稱邊界元法的真實頭模型剖分與構建
常用的基于醫學成像技術構建頭部模型的方法主要包括：邊界元法、有限元法和有限差分法。由于邊界元法是基于給定的電流源計算導體表面或者邊界電位，只需對邊界離散，可以降低邊界元的自由度，降低解的空間維數，因此本研究正向模型采用對稱邊界元方法剖分真實頭模型，并設計三個表面：大腦組織和顱骨之間的接觸面、顱骨和頭皮的接觸面以及頭皮和空氣的接觸面。這些接觸面光滑且通向傳導，每一個表面被網格化為小的三角形元素。均勻邊界元導體V上任意一點r上對應的電位V(r)定義如下：
| '/> |
其中 為電導率， 為點r上的電勢， 和 分別表示接觸面 內部和外部的導電率， 是與接觸面元素正交的向量， 表示接觸面元素的面積。每一個接觸面 被分割為N個三角形。每個三角形的電勢通過上式計算可得。接觸面 的電勢被定義為：
| '/> |
利用分割后和網格化后的數據構建頭部的真實物理模型。設置電導率：顱骨0.013 2 S/m，頭皮0.35 S/m，腦脊液1.79 S/m，皮層0.33 S/m，白質0.14 S/m。
1.3 腦電數據采集與預處理
本課題采用國際10-20系統采集癲癇患者睡眠期腦電，記錄了19個通道的腦電信號。根據受試者頭部電極的實際位置對部分電極點的位置進行了微調和校準。對采集到的EEG進行0.5～80 Hz的帶通濾波，并去除了50 Hz的工頻干擾。由于腦電源成像算法對噪聲敏感，因此采用獨立成分分析去除眼動偽跡和心電偽跡[18]。癲癇患者睡眠期EEG中與癇性發作相關的棘波由神經外科醫生手動標記，在每個被試數據中截取9～14個從棘波波峰前100 ms到棘波波峰后300 ms的EEG片段，并對這些片段進行疊加平均。
此外，在腦電誘發響應的逆問題研究中，由于作為背景噪聲的自發腦電的存在，EEG信噪比非常低。為了提高信噪比，本實驗利用特征值中位數進行噪聲協方差正規化，選取連續EEG記錄中沒有癲癇活動、快速眼動及其他偽影的50 s計算協方差矩陣，用以去除經預處理后依然存在的由儀器、傳感器和環境等帶來的噪聲。
2 原理與方法
2.1 基于多偶極子的腦電溯源算法
腦電逆問題主要分為基于空間約束的源成像和基于時空約束的源成像[19]。目前，由于分布式源成像模型具有一定的彌散性，能夠解決大部分腦彌漫性病變的溯源問題，如阿爾茨海默病[20]，但局灶性癲癇的致癇范圍相對比較局限，分布式源成像模型會帶來較高的定位假陽性。不同于其他腦電逆問題模型，基于偶極子模型的逆問題模型能夠更加精準地反演出局灶性神經元活動的位置，在致癇灶的定位問題中已有應用[21]。因此本文選擇通過計算復雜度低且定位范圍局限的偶極子模型定位大腦皮層中異常放電位置。
基于多偶極子逆問題算法的主要思想是通過調整偶極子的各個參數，擬合正向計算結果與頭皮電極測量數據。由于傳遞函數L為非線性函數，因此在未知源可用有限時，逆問題可表述為尋找一組源參數使得目標函數最小，其數學表達為 。本文中的最小目標函數通過最小二乘法選取：
|
其中，V為對應給定(r, q)的正向計算電位值，為電位實測值。偶極子位置及擬合優度如圖1所示。偶極子可以根據時間著色（見圖1a），也可以根據擬合優度著色（見圖1b），值越大，代表擬合程度越高。將多組偶極子顯示在同一個MRI查看器上，可以將偶極子擬合到一系列單獨的癲癇尖峰上，觀察偶極子是否形成一個簇。如圖1c所示，將優度值低于所選閾值的偶極子予以隱藏，僅顯示具有最佳擬合度的偶極子。
圖1
3D MRI及真實頭模型中的多個偶極子的位置及擬合優度
a. 時間著色偶極子；b. 擬合優度著色偶極子；c. 最佳擬合度偶極子
Figure1. The position and value of multi dipoles in 3D MRI and realistic head modela. time coloring dipoles; b. fitting coloring dipoles; c. the best fitting dipole
2.2 基于相位傳遞熵的大腦皮層網絡構建
相位傳遞熵（phase transfer entropy，PTE）[22]是一種不依賴先驗模型的因果關系度量工具，適用于各類非平穩和高維隨機過程。對于給定的時間滯后，從信號X到信號Y的相位傳遞熵定義為：
| '/> |
其中，、 及 分別為響應時間序列的概率分布。 和 分別是時間點 信號X與信號Y 的相位：，。本文中設定 = 1。
在構建腦網絡的問題中，該方法可以分析任意多個通道或大腦區域之間的時變因果關系，不僅可以顯示給定周期內的信號信息，還可以識別兩個系統的同步性，具有良好的統計顯著性[23-24]。因此，本文基于相位傳遞熵構建了不同腦區之間的功能連接網絡，通過信息的因果流向關系區分致癇區與受累區，進一步提高致癇灶定位精確度。
3 結果與分析
3.1 基于致癇灶溯源的皮層電流密度分布結果
本文回顧性地分析了8例局灶性癲癇患者基于多偶極子算法與真實頭模型的癲癇致癇灶腦電溯源結果。實驗的金標準為立體定位腦電圖（stereoelectroencephalogram，SEEG）監測到癇性發作的起始位置或致癇灶去除手術的預后效果。使用三種經典的逆向解算法——動態統計參數映射（dynamic statistical parameter mapping，dSPM）、sLORETA和分布式源法與本文提出的多偶極子算法進行比較。同時，在每種算法中對比了無約束條件和有額外約束條件的逆向解。此處以01號患者為例，患者性別男，39歲，有家族史，癲癇發作史2～3年，藥物控制效果差。依據各算法得到的該患者皮層電流密度結果如圖2所示。
圖2
通過不同逆問題算法得到的01號癲癇患者致癇灶定位結果
Figure2.
EZ localization results of patient No. 01 based on various algorithms
影像診斷提示右側杏仁核異常且不排除左側異常，臨床醫生通過癥狀學及腦電圖診斷為右側顳葉癲癇。隨后，該患者接受了電極植入手術，植入電極數目為7根，每根電極上分別有10個電極點，如圖3所示。
圖3
手術植入電極位置
Figure3.
The position of implanted SEEG
經過SEEG電極檢測，癇性發作起始期的棘波由電極A上的第二個電極點A2檢測到，即右側杏仁核位置。檢測到的棘波由A、B、C、D擴散到F電極，即異常放電由右側顳葉擴散到右側前葉、島葉。從圖2與圖3的對比分析中可以看出，四種算法對致癇灶的定側結果準確，但定位結果各不相同。其中，dSPM算法與sLORETA算法定位范圍較大，無法對植入電極或手術切除位置提供指導。電流分布密度算法雖定位較為局限，但定位在顳葉外側，與真實致癇灶位置不符。本文提出的多偶極子逆問題算法的溯源結果為顳葉內側杏仁核，在定位準確性和范圍局限性方面均優于其他三種算法。
此外，由于術前腦電及影像不能排除左側顳葉異常放電的可能，因此手術植入電極A’檢測該處是否為致癇區。如圖4所示，依據本實驗得到的致癇灶定位結果，將截取的9個400 ms棘波片段疊加平均，定義癇性發作棘波波峰出現的時刻為0.0 s，將此時的電流密度分布圖調整為只顯示大于12 μA的皮層電流密度區域，可以看到杏仁核與海馬放電異常。若將電流密度分布的顯示閾值調整為10 μA，可以看到大部分右側顳葉及小部分左側顳葉及杏仁核位置有異常電流。植入的電極A’覆蓋了左顳區域。接著，將時間軸移動到棘波波峰出現后的100 ms時刻，此時電流密度分布圖顯示左側顳葉、額葉及枕葉均有相對較高的異常信號電流，電流密度值在7～8 μA之間，低于異常放電起始時的電流密度值。從不同時刻大腦皮層的電流密度分布差異可以推斷，異常放電的起始區位于右側杏仁核及海馬，而右顳外側及左側顳葉、額葉、枕葉均為功能受累區。SEEG結果也證實了這一推斷。
圖4
01號患者癇性發作棘波出現時刻的大腦皮層電流密度分布圖
Figure4.
The current density distribution map on cortex during seizure spike of patient No. 01
3.2 癲癇患者腦功能網絡分析
通過以上實驗得到癲癇患者大腦皮層的電流密度分布，根據Desikan-Killiany地圖集將大腦皮層分為68個感興趣區域。在構建腦網絡的問題中，相位傳遞熵可以分析任意多個通道或大腦區域之間的時變因果關系，具有良好的統計顯著性和效力。因此，本文基于相位傳遞熵構建了不同腦區之間的功能連接網絡，如圖5所示。根據腦網絡圖，大腦皮層信息由右側顳橫回流向左側額下島蓋部，由右側頂下小葉流向左側顳葉內嗅皮層，并整體呈現出皮層信息由大腦右側半球流向左側的趨勢。在左側大腦半球功能區之間，信息由左側額極流向顳葉梭形回。因此可以推斷：致癇區在大腦右側，左側為受累區；致癇區在右側顳葉，左側額葉、顳葉為受累區，這與前文分析的致癇灶定側與定位結果表現出了高度的一致性。因此，若術前參考本次實驗結果，確定左側顳葉為癇性發作的功能受累區，而非起始發作區，則可不植入圖3中所示的電極A’。
圖5
基于相位傳遞熵構建的顳葉癲癇患者腦網絡
Figure5.
The brain network of temporal lobe epilepsy based on PTE
3.3 致癇灶源定位結果與分析
與正常人大腦不同，反復的癲癇發作會使癲癇患者的大腦發生結構或功能上的改變。基于癲癇患者結構核磁圖像重建的真實頭模型，不僅能夠在逆問題方程的計算中得到更加精確的結果，還能為手術方案的制定提供更加直觀的展示[25]。8位患者基于真實頭模型、多偶極子與相位傳遞熵網絡分析的致癇灶定位結果如圖6所示。
圖6
基于多偶極子的致癇灶溯源至大腦皮層的電流密度分布圖
Figure6.
The current density distribution map of EZ localization on cortex based on multi-poles
將基于4種算法得到的放電起始位置與SEEG檢測到的起始放電位置作比較，得到的定位準確性、敏感性、特異性和F1指標如表1所示。癲癇患者的溯源結果與臨床結果對比分析如表2所示。實驗結果表明，基于多偶極子溯源與相位傳遞熵網絡分析定位致癇灶結果的各項指標均優于其他三種算法。值得一提的是，多偶極子算法結果的特異性較高，即錯診率低，如果術前參考該腦電溯源結果，就能在癲癇手術中盡可能地避免損傷正常腦組織。同時，在本文的回顧性分析中，將溯源結果、臨床及影像診斷與SEEG金標準作對比，依據本實驗溯源結果制定的SEEG植入電極數目比原方案減少（19.34 ± 7.15）%，降低了因植入電極帶來的感染風險及手術費用。
表1
基于4種算法定位致癇灶的結果
Table1.
Localization result of EZ based on 4 algorithms
表2
癲癇患者的溯源結果與臨床結果對比分析
Table2.
Comparison and analysis between localization based on multi-poles and clinical results of patients
4 結論
本研究旨在通過生物電磁成像算法，實現癲癇致癇灶的無創、精準定位。本文基于癲癇患者的多模態醫學數據，通過T1 MRI構建癲癇患者的真實頭模型，引入多偶極子算法定位致癇灶，構建了大腦皮層網絡進一步區分致癇灶與其受累區。實驗結果表明，針對癲癇致癇灶定位問題，本文提出的基于真實頭模型與多偶極子算法的定位結果精確、范圍局限。該項工作不僅能夠提高致癇灶定位的準確性，還能減少術前檢查和手術操作中帶來的額外傷害和潛在風險，為神經外科醫生制定手術方案提供更加直觀、有效的參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：屈若為是本研究實驗設計、實驗研究、數據分析和論文寫作的執行人；王召楠負責實驗數據收集和預處理；王石峰負責醫學影像數據的收集及標注；王瑤負責腦網絡的構建與該部分知識的審核；王樂參與腦電實驗設計、收集癲癇患者腦電數據并標注；尹紹雅、顧軍華和徐桂芝教授是基金的負責人，對文章知識性內容作批評性審閱與指導，參與論文校對。
倫理聲明：本研究通過了河北工業大學生物醫學倫理委員會的審批（批文編號：HEBUThMEC2022003）。
