藥物難治性癲癇可通過手術進行治療，顱內腦電圖已廣泛應用于致癇區定位。目前，多數癲癇網絡分析依賴于致癇區節點的屬性，比如中心性、出入度等，上述特征有一定片面性，且難以應用到個體化患者治療中。本文提出了一種新的致癇區定位方法，通過引入神經計算模型，重構動態化的癲癇網絡，結合虛擬切除和空間鄰域擴張方法，來實現致癇區定位。本文分析了11位難治性癲癇患者的皮層腦電圖（ECoG）數據，利用患者發作間期數據和手術切除區域和療效進行分析。結果表明，在術后療效好的患者中，定位區域與手術區域的一致率明顯高于術后療效差的患者。術后效果好的患者和效果差的患者定位區域的平均偏離距離分別為15 mm和36 mm。預后分析顯示，以本文定位區域進行手術切除，當前療效差的患者能獲得更好的治療效果。本文方法能夠為難治性癲癇患者手術治療提供個體化方案，為癲癇外科提供一種量化的分析工具。
引用本文： 黃保強, 李春勝. 基于癲癇網絡動態重構與虛擬切除的致癇區定位研究. 生物醫學工程學雜志, 2022, 39(6): 1165-1172. doi: 10.7507/1001-5515.202205048 復制
引言
癲癇是由于大腦神經元異常同步放電而引起的腦網絡功能障礙性疾病[1]，大部分癲癇患者可以使用抗癲癇藥物治愈，但是仍有約30%的患者不能藥物治愈，最終發展為藥物難治性癲癇。手術切除致癇區是治療難治性癲癇的有效手段，而準確定位致癇灶是手術切除前的重要目標[2]。皮層腦電圖（electrocorticography，ECoG）可直接記錄大腦皮層神經電活動，其信號幅值比頭皮腦電圖高，信號源位置明確，可更精準地定位致癇區[3]。腦網絡分析可以幫助研究相關腦區之間的有效連接，對于理解大腦的功能機制和治療腦疾病具有重要意義[4]。越來越多的研究表明癲癇是由腦網絡紊亂引起，患者癇性放電使癲癇網絡中的某些關鍵區域表現異常，導致整個大腦網絡進入了超同步狀態[5]。通過對腦網絡的研究，可以深入理解癲癇發作機制，實現更準確的致癇區定位。癲癇發作與患者大腦局部腦網絡異常密切相關，產生發作的過程被認為是異質性和患者特異性的。由于異質性，癲癇手術治療并不總能有效控制癲癇發作，為個體患者選擇最佳治療方案面臨很大困難。因此在致癇區定位和療效預測中，將分析限定在患者自己的數據上，可為臨床個體化治療提供更適合的方案。
構建癲癇網絡的數據主要有腦電圖、功能磁共振成像和彌散張量成像等數據，其中最為常用的是腦電圖數據[6-7]。使用腦電圖數據構建腦網絡的常見方法有格蘭杰因果[8]、有向傳遞函數（directed transfer function，DTF）[9]、自適應有向傳遞函數[10]、偏定向相干[11]、獨立有效相干[12]等，不同的網絡構建方法，分析數據的側重點不同。用于表征網絡特征的屬性也較多。網絡特征路徑長度用于估計大腦區域之間功能整合的能力，可用網絡中兩個節點的平均最短路徑長度表示。聚類系數描述網絡中節點連接的緊密程度[13]。腦網絡具有小世界網絡屬性[14]。節點在網絡中與其他節點的連接強度和數量表示節點的重要性，稱為中心性。有研究基于DTF方法對顳葉癲癇顱內腦電圖進行分析，發現了臨床診斷的病灶區和非病灶區之間的因果流差異[15]。Hao等[16]研究難治性癲癇患者的ECoG腦網絡，分析了癲癇發作期致癇區的特征向量中心性變化規律。有研究使用立體定向腦電圖（stereo-electroencephalography，SEEG）信號和獨立有效相干方法構建腦網絡，分析網絡的信息流向，并進行K-均值（K-means）聚類分析，發現網絡節點出度值對顳葉癲癇致癇區有較高的分類準確率[17]。文獻[18]研究發現發作期SEEG信號的跨頻率耦合網絡的高頻中心節點可以更準確地預測手術結果。在腦網絡分析中，較高的介數中心性（betweenness centrality，BC）表明此區域為網絡信息流的關鍵樞紐[19]。上述研究方法是通過腦電圖構建一組患者的癲癇網絡，分析癲癇網絡特征及其與病灶區之間的關系，得出致癇區節點的特征，定位致癇區。但由于癲癇發作有結構、感染、遺傳和免疫等多種病因，癲癇網絡的腦區分布不盡相同，基于一組患者網絡特征分析的定位結果不一定是某位特定患者的最佳治療方案。目前藥物難治性癲癇手術治愈率在70%左右，為了進一步提高治愈率，癲癇腦網絡外科將更關注個體化的治療策略[20]。另外，每種網絡特征都有各自的局限性，只反映網絡某一方面特征，容易形成片面結果。
本文基于癲癇網絡動態重構與虛擬切除方法，提出了一種新的個體化致癇區定位方法。在構建患者癲癇網絡基礎上，利用神經計算模型仿真癲癇網絡動態變化，進行致癇節點虛擬切除，通過搜索網絡中對癲癇發作具有重要影響的節點并進行空間鄰域擴張，形成自動定位結果。本文對11位癲癇患者ECoG數據進行了致癇區定位分析。
1 數據和方法
1.1 數據信息與預處理
本研究使用的患者數據來源于IEEG公開數據集（https://www.ieeg.org），采用了11位患者的相關數據。根據ILAE癲癇手術結果量表，手術結果由至少1年的術后隨訪確定，本文將Ⅰ類和Ⅱ類歸類為術后效果良好，Ⅲ、Ⅳ、Ⅴ類歸類為術后效果較差。根據此分類標準，將11位患者數據大致分為兩類，即術后效果良好和較差，術后效果好的患者編號為P1、P2和P3，術后效果差的患者編號為P4～P11。術后效果良好患者的切除區被認為是癲癇灶所在區域[21]。患者數據信息如表1所示。
表1
患者數據信息
Table1.
Patient information
本文分析了每位癲癇患者發作間期600 s的ECoG數據，對ECoG數據做了雙極導聯處理，突顯局部電活動。本文對Delta（1～3 Hz）、Theta（4～7 Hz）、Alpha（8～13 Hz）、Beta（14～30 Hz）、Gamma（31～80 Hz）等頻段信號進行了分析，對雙極化ECoG信號進行零相位濾波。
1.2 DTF-替代數據方法
DTF是基于格蘭杰因果提出的自回歸模型分析方法，用來分析多通道腦電信號間的相互關系。神經活動之間的因果相互作用被稱為有效連接[22]，用DTF方法描述腦網絡的定向連接，在計算中考慮兩兩通道之間的直接影響，并考慮因果流向，從而得到不同節點之間的因果關系。DTF已被用于評估發作間期活動以預測癲癇發作，以及難治性癲癇患者的致癇區定位研究[23-24]。
根據DTF算法，
通道的腦電信號可以用含未知系數的多元自回歸模型來表示，即為：
 |
 |
其中p表示多元自回歸模型的階數，
表示N個通道的腦電信號，
表示的是 
的模型系數矩陣，
表示系統中的隨機白噪聲。
為了分析信號的頻譜特性，首先對式(2)作 Z 變換得到：
 |
再做傅里葉變換，該系統的傳遞函數可以用式(4)表示：
 |
式(4)中 f 表示的是頻率，
表示采樣間隔。
為了便于比較某一通道受其他通道的影響，通過歸一化的處理，得到式（5）。
 |
式(5)中的 
即為DTF矩陣，i 和 j 表示具體的腦網絡節點，式(5)描述了通道j 對通道 i 的影響占所有通道對通道 i 影響的比率。
考慮到直接和間接的影響，即 
對 
的影響和 
對 
的影響不同，DTF矩陣是非對稱的。DTF在分析神經信號時將腦電信號在短時間內看作平穩信號[25]，本文在計算DTF矩陣時，使用DTF-替代數據方法構建患者癲癇網絡，增加了200次數據替代，替代數據的方法打亂了原始數據的相位信息，可以保留時域和幅值特征，同時排除由于隨機因素產生的無意義連接，從而篩選出有效功能連接[26]。ECoG信號被分割為1 s時間窗，每個時間窗與上一時間窗重疊0.5 s，對一段時間內所有時間窗的DTF矩陣計算均值，得到平均DTF矩陣。
1.3 神經計算模型重構癲癇網絡
Benjamin等[27]提出了一種能反映多種癲癇活動狀態的計算模型，稱為Z6模型。模型將癲癇發作看作噪聲驅動的雙穩態系統，可用于模擬癲癇發作時節點之間信息相互傳遞、相互影響的動態過程，能直觀地描述網絡中節點的狀態轉換。神經計算模型是噪聲驅動的雙穩態系統，在噪聲的驅動下，系統可由平穩狀態轉變為癲癇發作的振蕩狀態。考慮到每個患者的差別，引入患者的癲癇網絡因果連接矩陣，模擬特定患者的癲癇發作過程。Z6模型網絡節點如下：
 |
式（6）中，
是一個復變量。
和 
是實常數，λ 是系統自身的吸引子，決定系統進入癲癇振蕩態的可能性。λ 選取在0和1之間，取值越趨近于1，模型興奮性越強。β 是網絡節點間連接強度的系數，神經計算模型的參數β選取在0.01和1之間。ω表示信號的振蕩頻率，對癲癇發作的影響較小[27]。本文仿真中 
，
，
。
是患者癲癇網絡矩陣，表示節點 j 對節點 i 的影響。N 是患者的網絡節點數，
是利用高斯白噪聲作為模擬大腦外部因素的驅動，系數是1.5。本文對神經計算模型仿真采用了Euler-Maruyama法，時間步長為0.05[28]。圖1是三節點網絡模型及其仿真結果。仿真開始時三個節點 
、
、
都處于平穩態（發作間期）。在外加噪聲的影響下，節點 
先于節點 
和 
進入振蕩態（癲癇發作）。從仿真開始時刻到節點進入振蕩態所用的時間稱為逃逸時間，如圖1所示。
圖1
三節點模型的逃逸時間示意圖
Figure1.
Escape time of a three-node network
1.4 空間鄰域擴張方法的致癇區定位
癲癇發作的起源位置被稱為致癇區，然后傳播至大腦其他區域，從而引起癲癇發作。癲癇的傳播可以通過對癲癇網絡的干預來控制[29]。癲癇發作是從致癇性高的節點開始，沿著周邊節點傳導，而切除致癇性最高的節點，能夠使整個腦網絡致癇性下降最多，達到減少和消除癲癇發作的效果。本文提出在重構的癲癇網絡中，采用虛擬切除手段，由致癇性最高的節點開始，沿著切除后使致癇性下降最多的節點擴張。空間鄰域擴張定位方法的流程如圖2所示。
圖2
空間鄰域擴張定位方法流程
Figure2.
Flowchart of the spatial neighbor expansion method
空間擴張方法具體步驟如下：
（1）使用患者ECoG信號，計算DTF腦網絡連接矩陣。
（2）將DTF網絡矩陣納入神經計算模型，通過改變腦網絡，依次虛擬切除單個節點，計算逃逸時間。
（3）以致癇性最高的節點為起點，在切除起點的基礎上，虛擬切除與起點相鄰的每個節點，選擇使逃逸時間增大最多的節點作為新定位的致癇節點。
（4）在新定位節點和此前定位節點的基礎上，以新定位節點為起點，按照步驟（3）的方法尋找下一個致癇節點，重復此過程。
（5）當切除新的節點使逃逸時間不再明顯增大，或新的節點周圍沒有直接相鄰的未切除節點，空間擴張過程結束。
（6）對定位結果做預后分析。
1.5 手術預后分析
根據文獻[28]方法來分析手術切除定位區域后的治療效果。腦網絡的致癇性用逃逸時間表述，每個通道都有自己的逃逸時間，取所有通道逃逸時間的均值，稱為平均逃逸時間T。T越大，則癲癇發作的可能性越低。在神經計算模型中，使用隨機噪聲驅動癲癇發作，通過多次的虛擬仿真消除隨機因素的影響。首先在未進行虛擬切除的時候，通過模擬仿真計算出平均逃逸時間T，然后在虛擬切除后計算出平均逃逸時間
，如果
，則說明此切除方案使腦網絡整體的致癇性降低了，可認為此切除方案對手術康復有利。然后進行隨機切除。本文的隨機切除是保留原切除方案中的節點，在剩下的節點中隨機選擇相同數量的節點，進行隨機切除仿真，計算出平均逃逸時間
。
總體而言，需要
，說明此手術切除方案有效。在此前提下，需要
，可進一步說明該切除方案優于隨機切除，方案可行；若
，則該方案不可行。
2 結果
2.1 定位區域與手術切除區相關性
本文采用重構癲癇網絡和空間鄰域擴張方法對患者致癇區定位，在腦電圖5個不同頻段上定位區域與臨床切除區域的平均一致率如圖3所示，手術結果良好的患者在腦電圖Delta、Theta、Alpha、Beta、Gamma頻段上一致率的平均值分別為：35.6%、35.7%、47.6%、30.3%、62.1%。定位結果在Gamma頻段上的平均一致率最高，且在手術結果良好與較差患者間的差別最大。術后效果良好的患者在Gamma頻段的定位一致率均大于50%，術后效果差的患者中，P4、P5、P9和P11的定位一致率為0%，如表2所示。本文進一步對采用腦電圖Gamma頻段信號的定位結果進行了分析。
圖3
各頻段定位區域與臨床切除區域一致性
Figure3.
Consistency between the localization area of each frequency band and the clinical resection area
表2
Gamma頻段致癇區定位結果
Table2.
Localizing results of epileptogenic zone using gamma rhythm
2.2 致癇區定位結果
采用本文方法和腦電圖Gamma頻段數據的定位區域和逃逸時間變化如表2所示。圖4顯示了6位患者的致癇區定位過程和結果。對于3位術后效果好的患者P1、P2和P3，空間鄰域擴張定位方法的起點位于切除區，即定位的致癇性最高的區域位于切除區。而3位術后效果差的患者P4、P5和P6，空間鄰域擴張定位方法的起點位于切除區之外的區域。
圖4
致癇區定位區域分布
Figure4.
Distribution of the localized epileptogenic zone
進一步計算了本文定位區域與手術切除區不一致部分的偏離距離，每個患者定位區域的偏離距離如圖5a所示。平均偏離距離用來說明兩組患者的定位區域相對于切除區域的偏離情況，結果顯示手術效果好和效果差的患者平均偏離距離分別為15 mm和36 mm，前者平均偏離距離明顯小于后者，如圖5b所示。
圖5
定位區域的平均偏離距離
a. 患者個體的偏離距離；b. 兩組患者的平均偏離距離
Figure5. The average deviation distance of the localized areaa. deviation distance of individual patients; b. mean deviation distance for both groups of patients
2.3 定位區域的預后分析
進一步分析以本文定位區域為手術切除方案的可行性，對定位結果進行了預后分析。圖6顯示了患者P1～P6定位區域虛擬切除后的結果。在未進行虛擬切除的時候，通過模擬仿真計算出平均逃逸時間 T，虛擬切除本文定位區域，計算出平均逃逸時間 
，以及隨機切除的平均逃逸時間 
。每個切除方案進行3 000次計算，取平均逃逸時間。結果顯示，如以本文方法定位區域進行手術切除，所有11位患者都將取得較好的手術結果（
），見表2。
圖6
切除定位區域的預后結果
Figure6.
Surgical outcome prediction using the localized area
3 討論
本文基于DTF方法構建患者的致癇腦網絡，利用神經計算模型，結合虛擬切除和空間擴張的方法定位致癇區。在3位術后效果良好的患者中，采用腦電圖高頻率段Gamma的定位區域與臨床切除區域的一致率最高，與文獻[30]的結果相近，說明癲癇腦電圖中高頻活動與致癇灶關系密切[31]。本文方法定位區域與臨床切除區域的一致率高于50%，且定位方法的起點都位于切除區內，說明本文提出的方法能有效定位致癇性高的區域。術后效果差的患者中，定位區域與切除區域的一致率明顯低于手術效果良好患者的一致率，平均偏離距離明顯大于手術效果良好患者的偏離距離，說明采用本文方法可以更好地確定致癇區。預后分析顯示，如果按本文方法定位的區域進行切除，8位術后效果差的患者均能獲得更好的手術治療效果。
文獻[32]使用19名患者的術前彌散張量成像數據得到個體的結構網絡，以臨床假設的致癇區的特征向量中心性（eigenvector centrality，EC）為核心，使用“易感-感染-恢復”[33]（susceptible-infected-recovered，SIR）模型模擬切除EC到大腦各區域的連接，并與隨機切除相同數量的連接進行比較，發現要達到90%的效果可以通過移除顯著少于90%的連接來獲得。因此，范圍更小、優化的虛擬切除實現了與實際手術相近的效果，而成本要低得多。本文選擇DTF方法構建癲癇網絡，基于癲癇發作動力學模型和虛擬切除的方法，分析方法與文獻[32]具有一致性。
本文分析采用了顱內ECoG數據，顱內腦電圖具有測量信號準確和位置明確的優點，而其記錄區域常受限于特定腦區，空間覆蓋往往不夠廣泛。在利用顱內腦電圖進行癲癇網絡分析時需要充分考慮上述限制因素，來構建全面準確的癲癇網絡。另外，本研究采用的患者數據樣本量仍相對較少，為了實現本文方法的臨床應用，將在下一步與臨床單位合作，利用多中心大樣本數據進行分析改進。
4 結論
致癇區定位是癲癇外科治療的核心問題。本文在利用神經計算模型重構患者癲癇網絡的基礎上，引入虛擬切除方法定位網絡中高致癇性節點，可為難治性癲癇患者手術治療提供個體化方案。通過對11位癲癇患者ECoG數據、手術區域和療效的深入分析，驗證了本文方法的有效性。本文提出的致癇區定位和預后分析相結合的方法，可為臨床癲癇手術治療提供一種量化的分析工具。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李春勝和黃保強共同提出本文致癇區定位方法，黃保強負責數據處理和分析，李春勝和黃保強共同完成本文的寫作及修改。
引言
癲癇是由于大腦神經元異常同步放電而引起的腦網絡功能障礙性疾病[1]，大部分癲癇患者可以使用抗癲癇藥物治愈，但是仍有約30%的患者不能藥物治愈，最終發展為藥物難治性癲癇。手術切除致癇區是治療難治性癲癇的有效手段，而準確定位致癇灶是手術切除前的重要目標[2]。皮層腦電圖（electrocorticography，ECoG）可直接記錄大腦皮層神經電活動，其信號幅值比頭皮腦電圖高，信號源位置明確，可更精準地定位致癇區[3]。腦網絡分析可以幫助研究相關腦區之間的有效連接，對于理解大腦的功能機制和治療腦疾病具有重要意義[4]。越來越多的研究表明癲癇是由腦網絡紊亂引起，患者癇性放電使癲癇網絡中的某些關鍵區域表現異常，導致整個大腦網絡進入了超同步狀態[5]。通過對腦網絡的研究，可以深入理解癲癇發作機制，實現更準確的致癇區定位。癲癇發作與患者大腦局部腦網絡異常密切相關，產生發作的過程被認為是異質性和患者特異性的。由于異質性，癲癇手術治療并不總能有效控制癲癇發作，為個體患者選擇最佳治療方案面臨很大困難。因此在致癇區定位和療效預測中，將分析限定在患者自己的數據上，可為臨床個體化治療提供更適合的方案。
構建癲癇網絡的數據主要有腦電圖、功能磁共振成像和彌散張量成像等數據，其中最為常用的是腦電圖數據[6-7]。使用腦電圖數據構建腦網絡的常見方法有格蘭杰因果[8]、有向傳遞函數（directed transfer function，DTF）[9]、自適應有向傳遞函數[10]、偏定向相干[11]、獨立有效相干[12]等，不同的網絡構建方法，分析數據的側重點不同。用于表征網絡特征的屬性也較多。網絡特征路徑長度用于估計大腦區域之間功能整合的能力，可用網絡中兩個節點的平均最短路徑長度表示。聚類系數描述網絡中節點連接的緊密程度[13]。腦網絡具有小世界網絡屬性[14]。節點在網絡中與其他節點的連接強度和數量表示節點的重要性，稱為中心性。有研究基于DTF方法對顳葉癲癇顱內腦電圖進行分析，發現了臨床診斷的病灶區和非病灶區之間的因果流差異[15]。Hao等[16]研究難治性癲癇患者的ECoG腦網絡，分析了癲癇發作期致癇區的特征向量中心性變化規律。有研究使用立體定向腦電圖（stereo-electroencephalography，SEEG）信號和獨立有效相干方法構建腦網絡，分析網絡的信息流向，并進行K-均值（K-means）聚類分析，發現網絡節點出度值對顳葉癲癇致癇區有較高的分類準確率[17]。文獻[18]研究發現發作期SEEG信號的跨頻率耦合網絡的高頻中心節點可以更準確地預測手術結果。在腦網絡分析中，較高的介數中心性（betweenness centrality，BC）表明此區域為網絡信息流的關鍵樞紐[19]。上述研究方法是通過腦電圖構建一組患者的癲癇網絡，分析癲癇網絡特征及其與病灶區之間的關系，得出致癇區節點的特征，定位致癇區。但由于癲癇發作有結構、感染、遺傳和免疫等多種病因，癲癇網絡的腦區分布不盡相同，基于一組患者網絡特征分析的定位結果不一定是某位特定患者的最佳治療方案。目前藥物難治性癲癇手術治愈率在70%左右，為了進一步提高治愈率，癲癇腦網絡外科將更關注個體化的治療策略[20]。另外，每種網絡特征都有各自的局限性，只反映網絡某一方面特征，容易形成片面結果。
本文基于癲癇網絡動態重構與虛擬切除方法，提出了一種新的個體化致癇區定位方法。在構建患者癲癇網絡基礎上，利用神經計算模型仿真癲癇網絡動態變化，進行致癇節點虛擬切除，通過搜索網絡中對癲癇發作具有重要影響的節點并進行空間鄰域擴張，形成自動定位結果。本文對11位癲癇患者ECoG數據進行了致癇區定位分析。
1 數據和方法
1.1 數據信息與預處理
本研究使用的患者數據來源于IEEG公開數據集（https://www.ieeg.org），采用了11位患者的相關數據。根據ILAE癲癇手術結果量表，手術結果由至少1年的術后隨訪確定，本文將Ⅰ類和Ⅱ類歸類為術后效果良好，Ⅲ、Ⅳ、Ⅴ類歸類為術后效果較差。根據此分類標準，將11位患者數據大致分為兩類，即術后效果良好和較差，術后效果好的患者編號為P1、P2和P3，術后效果差的患者編號為P4～P11。術后效果良好患者的切除區被認為是癲癇灶所在區域[21]。患者數據信息如表1所示。
表1
患者數據信息
Table1.
Patient information
本文分析了每位癲癇患者發作間期600 s的ECoG數據，對ECoG數據做了雙極導聯處理，突顯局部電活動。本文對Delta（1～3 Hz）、Theta（4～7 Hz）、Alpha（8～13 Hz）、Beta（14～30 Hz）、Gamma（31～80 Hz）等頻段信號進行了分析，對雙極化ECoG信號進行零相位濾波。
1.2 DTF-替代數據方法
DTF是基于格蘭杰因果提出的自回歸模型分析方法，用來分析多通道腦電信號間的相互關系。神經活動之間的因果相互作用被稱為有效連接[22]，用DTF方法描述腦網絡的定向連接，在計算中考慮兩兩通道之間的直接影響，并考慮因果流向，從而得到不同節點之間的因果關系。DTF已被用于評估發作間期活動以預測癲癇發作，以及難治性癲癇患者的致癇區定位研究[23-24]。
根據DTF算法， 通道的腦電信號可以用含未知系數的多元自回歸模型來表示，即為：
|
|
其中p表示多元自回歸模型的階數， 表示N個通道的腦電信號， 表示的是 的模型系數矩陣， 表示系統中的隨機白噪聲。
為了分析信號的頻譜特性，首先對式(2)作 Z 變換得到：
|
再做傅里葉變換，該系統的傳遞函數可以用式(4)表示：
|
式(4)中 f 表示的是頻率， 表示采樣間隔。
為了便于比較某一通道受其他通道的影響，通過歸一化的處理，得到式（5）。
|
式(5)中的 即為DTF矩陣，i 和 j 表示具體的腦網絡節點，式(5)描述了通道j 對通道 i 的影響占所有通道對通道 i 影響的比率。
考慮到直接和間接的影響，即 對 的影響和 對 的影響不同，DTF矩陣是非對稱的。DTF在分析神經信號時將腦電信號在短時間內看作平穩信號[25]，本文在計算DTF矩陣時，使用DTF-替代數據方法構建患者癲癇網絡，增加了200次數據替代，替代數據的方法打亂了原始數據的相位信息，可以保留時域和幅值特征，同時排除由于隨機因素產生的無意義連接，從而篩選出有效功能連接[26]。ECoG信號被分割為1 s時間窗，每個時間窗與上一時間窗重疊0.5 s，對一段時間內所有時間窗的DTF矩陣計算均值，得到平均DTF矩陣。
1.3 神經計算模型重構癲癇網絡
Benjamin等[27]提出了一種能反映多種癲癇活動狀態的計算模型，稱為Z6模型。模型將癲癇發作看作噪聲驅動的雙穩態系統，可用于模擬癲癇發作時節點之間信息相互傳遞、相互影響的動態過程，能直觀地描述網絡中節點的狀態轉換。神經計算模型是噪聲驅動的雙穩態系統，在噪聲的驅動下，系統可由平穩狀態轉變為癲癇發作的振蕩狀態。考慮到每個患者的差別，引入患者的癲癇網絡因果連接矩陣，模擬特定患者的癲癇發作過程。Z6模型網絡節點如下：
|
式（6）中， 是一個復變量。 和 是實常數，λ 是系統自身的吸引子，決定系統進入癲癇振蕩態的可能性。λ 選取在0和1之間，取值越趨近于1，模型興奮性越強。β 是網絡節點間連接強度的系數，神經計算模型的參數β選取在0.01和1之間。ω表示信號的振蕩頻率，對癲癇發作的影響較小[27]。本文仿真中 ，，。 是患者癲癇網絡矩陣，表示節點 j 對節點 i 的影響。N 是患者的網絡節點數， 是利用高斯白噪聲作為模擬大腦外部因素的驅動，系數是1.5。本文對神經計算模型仿真采用了Euler-Maruyama法，時間步長為0.05[28]。圖1是三節點網絡模型及其仿真結果。仿真開始時三個節點 、、 都處于平穩態（發作間期）。在外加噪聲的影響下，節點 先于節點 和 進入振蕩態（癲癇發作）。從仿真開始時刻到節點進入振蕩態所用的時間稱為逃逸時間，如圖1所示。
圖1
三節點模型的逃逸時間示意圖
Figure1.
Escape time of a three-node network
1.4 空間鄰域擴張方法的致癇區定位
癲癇發作的起源位置被稱為致癇區，然后傳播至大腦其他區域，從而引起癲癇發作。癲癇的傳播可以通過對癲癇網絡的干預來控制[29]。癲癇發作是從致癇性高的節點開始，沿著周邊節點傳導，而切除致癇性最高的節點，能夠使整個腦網絡致癇性下降最多，達到減少和消除癲癇發作的效果。本文提出在重構的癲癇網絡中，采用虛擬切除手段，由致癇性最高的節點開始，沿著切除后使致癇性下降最多的節點擴張。空間鄰域擴張定位方法的流程如圖2所示。
圖2
空間鄰域擴張定位方法流程
Figure2.
Flowchart of the spatial neighbor expansion method
空間擴張方法具體步驟如下：
（1）使用患者ECoG信號，計算DTF腦網絡連接矩陣。
（2）將DTF網絡矩陣納入神經計算模型，通過改變腦網絡，依次虛擬切除單個節點，計算逃逸時間。
（3）以致癇性最高的節點為起點，在切除起點的基礎上，虛擬切除與起點相鄰的每個節點，選擇使逃逸時間增大最多的節點作為新定位的致癇節點。
（4）在新定位節點和此前定位節點的基礎上，以新定位節點為起點，按照步驟（3）的方法尋找下一個致癇節點，重復此過程。
（5）當切除新的節點使逃逸時間不再明顯增大，或新的節點周圍沒有直接相鄰的未切除節點，空間擴張過程結束。
（6）對定位結果做預后分析。
1.5 手術預后分析
根據文獻[28]方法來分析手術切除定位區域后的治療效果。腦網絡的致癇性用逃逸時間表述，每個通道都有自己的逃逸時間，取所有通道逃逸時間的均值，稱為平均逃逸時間T。T越大，則癲癇發作的可能性越低。在神經計算模型中，使用隨機噪聲驅動癲癇發作，通過多次的虛擬仿真消除隨機因素的影響。首先在未進行虛擬切除的時候，通過模擬仿真計算出平均逃逸時間T，然后在虛擬切除后計算出平均逃逸時間，如果，則說明此切除方案使腦網絡整體的致癇性降低了，可認為此切除方案對手術康復有利。然后進行隨機切除。本文的隨機切除是保留原切除方案中的節點，在剩下的節點中隨機選擇相同數量的節點，進行隨機切除仿真，計算出平均逃逸時間。
總體而言，需要，說明此手術切除方案有效。在此前提下，需要，可進一步說明該切除方案優于隨機切除，方案可行；若，則該方案不可行。
2 結果
2.1 定位區域與手術切除區相關性
本文采用重構癲癇網絡和空間鄰域擴張方法對患者致癇區定位，在腦電圖5個不同頻段上定位區域與臨床切除區域的平均一致率如圖3所示，手術結果良好的患者在腦電圖Delta、Theta、Alpha、Beta、Gamma頻段上一致率的平均值分別為：35.6%、35.7%、47.6%、30.3%、62.1%。定位結果在Gamma頻段上的平均一致率最高，且在手術結果良好與較差患者間的差別最大。術后效果良好的患者在Gamma頻段的定位一致率均大于50%，術后效果差的患者中，P4、P5、P9和P11的定位一致率為0%，如表2所示。本文進一步對采用腦電圖Gamma頻段信號的定位結果進行了分析。
圖3
各頻段定位區域與臨床切除區域一致性
Figure3.
Consistency between the localization area of each frequency band and the clinical resection area
表2
Gamma頻段致癇區定位結果
Table2.
Localizing results of epileptogenic zone using gamma rhythm
2.2 致癇區定位結果
采用本文方法和腦電圖Gamma頻段數據的定位區域和逃逸時間變化如表2所示。圖4顯示了6位患者的致癇區定位過程和結果。對于3位術后效果好的患者P1、P2和P3，空間鄰域擴張定位方法的起點位于切除區，即定位的致癇性最高的區域位于切除區。而3位術后效果差的患者P4、P5和P6，空間鄰域擴張定位方法的起點位于切除區之外的區域。
圖4
致癇區定位區域分布
Figure4.
Distribution of the localized epileptogenic zone
進一步計算了本文定位區域與手術切除區不一致部分的偏離距離，每個患者定位區域的偏離距離如圖5a所示。平均偏離距離用來說明兩組患者的定位區域相對于切除區域的偏離情況，結果顯示手術效果好和效果差的患者平均偏離距離分別為15 mm和36 mm，前者平均偏離距離明顯小于后者，如圖5b所示。
圖5
定位區域的平均偏離距離
a. 患者個體的偏離距離；b. 兩組患者的平均偏離距離
Figure5. The average deviation distance of the localized areaa. deviation distance of individual patients; b. mean deviation distance for both groups of patients
2.3 定位區域的預后分析
進一步分析以本文定位區域為手術切除方案的可行性，對定位結果進行了預后分析。圖6顯示了患者P1～P6定位區域虛擬切除后的結果。在未進行虛擬切除的時候，通過模擬仿真計算出平均逃逸時間 T，虛擬切除本文定位區域，計算出平均逃逸時間 ，以及隨機切除的平均逃逸時間 。每個切除方案進行3 000次計算，取平均逃逸時間。結果顯示，如以本文方法定位區域進行手術切除，所有11位患者都將取得較好的手術結果（），見表2。
圖6
切除定位區域的預后結果
Figure6.
Surgical outcome prediction using the localized area
3 討論
本文基于DTF方法構建患者的致癇腦網絡，利用神經計算模型，結合虛擬切除和空間擴張的方法定位致癇區。在3位術后效果良好的患者中，采用腦電圖高頻率段Gamma的定位區域與臨床切除區域的一致率最高，與文獻[30]的結果相近，說明癲癇腦電圖中高頻活動與致癇灶關系密切[31]。本文方法定位區域與臨床切除區域的一致率高于50%，且定位方法的起點都位于切除區內，說明本文提出的方法能有效定位致癇性高的區域。術后效果差的患者中，定位區域與切除區域的一致率明顯低于手術效果良好患者的一致率，平均偏離距離明顯大于手術效果良好患者的偏離距離，說明采用本文方法可以更好地確定致癇區。預后分析顯示，如果按本文方法定位的區域進行切除，8位術后效果差的患者均能獲得更好的手術治療效果。
文獻[32]使用19名患者的術前彌散張量成像數據得到個體的結構網絡，以臨床假設的致癇區的特征向量中心性（eigenvector centrality，EC）為核心，使用“易感-感染-恢復”[33]（susceptible-infected-recovered，SIR）模型模擬切除EC到大腦各區域的連接，并與隨機切除相同數量的連接進行比較，發現要達到90%的效果可以通過移除顯著少于90%的連接來獲得。因此，范圍更小、優化的虛擬切除實現了與實際手術相近的效果，而成本要低得多。本文選擇DTF方法構建癲癇網絡，基于癲癇發作動力學模型和虛擬切除的方法，分析方法與文獻[32]具有一致性。
本文分析采用了顱內ECoG數據，顱內腦電圖具有測量信號準確和位置明確的優點，而其記錄區域常受限于特定腦區，空間覆蓋往往不夠廣泛。在利用顱內腦電圖進行癲癇網絡分析時需要充分考慮上述限制因素，來構建全面準確的癲癇網絡。另外，本研究采用的患者數據樣本量仍相對較少，為了實現本文方法的臨床應用，將在下一步與臨床單位合作，利用多中心大樣本數據進行分析改進。
4 結論
致癇區定位是癲癇外科治療的核心問題。本文在利用神經計算模型重構患者癲癇網絡的基礎上，引入虛擬切除方法定位網絡中高致癇性節點，可為難治性癲癇患者手術治療提供個體化方案。通過對11位癲癇患者ECoG數據、手術區域和療效的深入分析，驗證了本文方法的有效性。本文提出的致癇區定位和預后分析相結合的方法，可為臨床癲癇手術治療提供一種量化的分析工具。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李春勝和黃保強共同提出本文致癇區定位方法，黃保強負責數據處理和分析，李春勝和黃保強共同完成本文的寫作及修改。
