揭示癲癇過程中腦內信息交互變化對臨床診療具有重要意義。早期研究表明，相位-幅度耦合（PAC）可作為定位致癇區域和表征癲癇狀態轉換的潛在生物標志物。然而，目前癲癇的相位-幅度耦合研究多集中于θ-γ耦合，對β-γ耦合在癲癇過程中的變化及潛在應用知之甚少。對此，本文基于癲癇頭皮腦電（EEG）數據，采用調節指數（MI）的PAC方法，探究癲癇不同狀態下β-γ耦合的變化。結果發現，癲癇過程中各腦區內β-γ耦合隨狀態發生改變，多數腦區內的β-γ耦合在發作期減小，在發作后期增加，其差異具有統計學意義。癲癇過程中不同腦區間也存在β-γ耦合，且β-γ耦合強度在發作后期增加，其差異具有統計學意義。這些研究結果既有助于理解癲癇過程中腦內信息交互的改變，也為癲癇的臨床診療提供了新參考。
引用本文： 李凱杰, 盧俊峰, 余仁萍, 張銳, 陳明明. 癲癇過程中頭皮腦電信號的β-γ耦合變化. 生物醫學工程學雜志, 2023, 40(4): 700-708, 717. doi: 10.7507/1001-5515.202212024 復制
0 引言
癲癇（epilepsy，EP）是一種由神經元異常放電引起短暫大腦功能障礙的神經系統疾病[1]。反復發作是癲癇主要的臨床表現，嚴重影響患者身心健康，甚至造成患者中樞神經系統損傷[2]。當前，盡管癲癇臨床診療水平較高，多數癲癇發作能夠得到及時有效的控制，但全球范圍內仍有約30%的癲癇患者對抗癲癇藥物等不敏感，發展成為難治性癲癇[3]。而且，致癇因素的多樣性和癲癇過程中復雜的腦內信息交互嚴重阻礙了臨床準確評估和有效治療難治性癲癇。因此，揭示癲癇過程中腦內信息交互隨癲癇狀態的改變，對臨床認識和診療癲癇具有重要的指導和參考意義。
腦電圖（electroencephalography，EEG）監測技術是通過置于頭皮上的電極獲取大腦內同步神經元活動產生的EEG信號，在臨床癲癇診療中應用廣泛[4]。癲癇EEG信號含有豐富的生理和病理信息，在癲癇發作全程常伴有明顯的特征變化，而解析癲癇過程EEG信號的動態變化特征有助于預測癲癇發作、探究癲癇播散過程和揭示致癇機制等。諸如，Das等[5]基于癲癇EEG信號，采用特征提取、模式匹配等算法組合而成的模型，能夠實現對癲癇發作的準確預測。Lopes等[6]應用鎖相因子構建功能網絡，提出發作生成傳播（ictogenic spread，IS）概念，并基于全面性癲癇患者和內側顳葉癲癇患者的EEG信號，揭示癲癇發作病理性活動在整個大腦網絡中的傳播特性。Mormann等[7]應用平均相位相干性統計度量相位同步差異，表征大腦動力學狀態變化，指出癲癇患者EEG信號的相位同步性變化可能有助于深入了解癲癇發作全過程。事實上，癲癇EEG信號除相位改變之外，幅度也會發生不同程度的改變。這些相位和幅度的改變與癲癇過程中腦內信息交互的轉變可能有密切關系。
相位—幅度耦合（phase-amplitude coupling，PAC）通過計算EEG信號的低頻相位與高頻幅度間的耦合關系，反映了腦區間的信息交互，在理解神經元網絡通信和連接機制等方面具有重要作用[8-9]。諸如在癲癇研究中，PAC可呈現不同腦區信息交互上的差異性，展現了腦內信息交互的特異性改變，是輔助診療癲癇的潛在生物標記物。近期， Ma等[10]在額葉癲癇的顱內EEG信號中發現，發作前期和間期的PAC可作為定位致癇灶的潛在標記物，且耦合的首選相位隨癲癇發作規律性改變 。類似的，基于小鼠的實驗研究發現，小鼠癲癇發作期皮質耦合受損，弱化的θ-γ耦合可能是癲癇發作風險預測的早期指標[11]。另外，在帕金森病的研究中，Salimpour等[12]發現皮層的相位靶刺激調節運動皮層中異常PAC，指出PAC在神經疾病治療上具有參考作用。
然而，需要注意的是，癲癇研究中PAC的分析多集中于θ-γ耦合，對β-γ耦合在癲癇過程中的變化知之甚少[13-14]。為此，本文將基于難治性癲癇患者的EEG信號，采用調節指數（modulation index，MI）的PAC分析方法，探究癲癇過程中β-γ耦合的存在性及隨癲癇狀態的變化。通過統計比較癲癇過程通道內和通道間的β-γ耦合強度變化，探究癲癇過程中腦內信息交互改變。通過本文研究，期望可為理解癲癇不同階段腦區間的信息交互提供新思路，為臨床診療癲癇提供理論借鑒。
1 材料和方法
1.1 數據采集及預處理
本文使用的EEG數據來源于波士頓兒童醫院（Children’s Hospital Boston，CHB）和麻省理工學院（Massachusetts Institute of Technology，MIT）創建的公開數據庫CHB-MIT頭皮EEG數據庫（scalp EEG Database）[15]。
由于數據集中每個案例的發作次數、單次記錄時長等均存在差異，為使數據在分析中保持一致性，本文在選取數據時，要求選取需滿足單次發作時長大于5 s，且兩次發作間隔超過1 h的數據，故選取了22個滿足上述條件的案例。針對每個案例，本文選取時長為6 h的數據，共計65次發作。同時，為避免不同案例電極數量、電極分布位置和壞導等因素對結果的影響，本文選取前18個通道進行分析。另外，為避免噪聲、偽跡（眼電、心電和肌電信號）等干擾，本文對上述選取的原始EEG數據進行了預處理。具體流程有：①采用平均參考方法對原始EEG數據進行重參考處理。②使用陷波濾波器濾除58～62 Hz之間的信號，消除60 Hz工頻干擾。③應用獨立成分分析（independent component analysis，ICA）去除EEG信號中的噪聲偽跡[16-17]。
最后，對預處理后的EEG數據，本文以數據集標記的發作時間為準，進行了分段截取，如圖1所示。假設發作期時長為L s，本文中定義發作前期為發作前3 × L s（3L），發作后期為發作后3 × L s（3L）。
圖1
數據選擇示意圖
Figure1.
Data selection diagram
1.2 方法
1.2.1 計算方法
本文參考文獻[18]的方法，采用MI 計算PAC。MI度量的是信號高頻幅度在低頻相位上的分布相對于均勻分布的偏差大小，是通過巴克-萊布勒（Kullback-Leibler）距離實現[19]。MI方法具有靈敏度高，對數據長度、信噪比和采樣率等影響因素魯棒性強[20]的特點，其具體計算過程如下：
（1）假設有原始信號x(t)，對其進行帶通濾波，可得到低頻段和高頻段信號，分別記為xfP(t)和xfA(t)；
（2）采用希爾伯特變換分別計算得到低頻段信號xfP(t)的相位時間序列φfP(t)和高頻段信號xfA(t)的幅度包絡AfA(t)；
（3）對相位時間序列φfP(t)進行均勻區間劃分，計算每個區間相位對應幅度包絡AfA(t)的平均值AfAmean；
（4）歸一化AfAmean，如式(1)所示：
 |
其中，AfAmean(j)是第j個區間相位對應的幅度包絡均值；N是區間個數，通常N = 18；P(j)是第j個區間相位標準化的幅度包絡均值；k為第k個區間相位。
（5）計算Kullback-Leibler距離：假設存在離散分布P和Q，其對應的Kullback-Leibler距離如式(2)所示：
 |
其中，DKL(P, Q) ≥ 0，當且僅當P = Q時，DKL(P, Q) = 0；P(j)是離散分布P第j個區間的概率密度；Q(j)是離散分布Q第j個區間的概率密度。式(2)中，P分布的香農熵H(P)定義如式(3)所示：
 |
因此，結合式(2)和式(3)可得出Kullback-Leibler距離與香農熵的關系，如式(4)所示：
 |
其中，U是均勻分布，log(N)是最大熵值。
（6）綜合上述公式，定義MI（以符號MI表示），如式(5)所示：
 |
由MI的定義可知，MI的值在0～1之間。當高頻幅值均勻分布在低頻相位上時，MI為0，即滿足j = 1, ???, N和P( j) = 1/N，在其它情況下，MI不為0。MI的大小與離散分布P和均勻分布U的偏離程度呈正相關。在MI計算中，本文設定低頻段相位頻率為13～30 Hz，高頻段幅值頻率范圍為30～100 Hz，步長為1 Hz。
1.2.2 統計分析
本文將發作期數據劃分為5 s/段，然后分別計算每段數據β、γ頻率范圍內MI值的平均值，得到一組β-γ耦合值，并以該組β-γ耦合值的平均值作為此次癲癇發作期的MI值。其中，若存在數據段對應的β、γ頻率范圍內MI值均為0，則將該數據段對應的β-γ耦合值置為0；然后將發作前期和發作后期的數據進行類似的統計分析；最終發作前期、發作期和發作后期，各自分別獲得65個平均MI值，共可獲得195個平均MI值。由于三種癲癇狀態對應的MI值序列均不符合正態分布，故本文采用非參數統計秩和檢驗進行對比分析，即三組數據兩兩進行秩和檢驗，檢驗水準為0.05。
另外，在計算MI的過程中，為避免數據隨機性引入的偽MI干擾，本文采用了200次的隨機替代數據進行統計比較，當兩者的差異應具有統計學意義時，保留原始數據計算的MI值；否則，令MI值為0。其中，檢驗水準為0.05。
2 實驗結果
2.1 癲癇過程不同階段通道內β-γ耦合
早期研究表明，癲癇患者EEG信號的低頻（θ頻段）相位和高頻（γ頻段）振幅間存在PAC（θ-γ耦合），且θ-γ耦合對臨床癲癇診療有重要參考作用[1]。相對而言，癲癇過程中β-γ耦合是否存在及如何變化仍不清楚。為此，本文隨機選取了三個標記為FP1-F7、T7-P7和P3-O1的通道信號，這些EEG信號分別位于大腦皮層不同區域。如圖2所示，在發作前期、發作期和發作后期，每個通道信號中都可觀察到β-γ耦合。雖然三個通道信號的MI值不同，但這些MI值在發作演化過程中的變化趨勢幾乎是一致的，尤其發作后期MI值明顯增加。值得注意的是，β-γ耦合可以出現在振幅的不同中心頻率上。這些觀測結果表明，β-γ耦合不僅存在于不同的癲癇階段，而且隨癲癇狀態的演化而改變，這可為癲癇發作預測提供新的視角。
圖2
癲癇過程不同階段β-γ耦合變化
Figure2.
Alterations of β-γ coupling in different epileptic phases
2.2 癲癇過程不同階段通道內β-γ耦合的統計比較
如前所述，量化癲癇β-γ耦合在癲癇過程不同階段的變化，可為理解癲癇發作的動態演化提供新視角。因此，本文計算了每個通道內β-γ耦合的MI值，并對發作前期、發作期和發作后期三個狀態的MI值進行統計分析比較。如圖3所示，左半腦通道EEG信號在三個時期的MI值變化趨勢不同。事實上，除了FP1-F3和FZ-CZ電極通道的差異不具有統計學意義外，左半腦其它通道的EEG信號MI值在癲癇發作期的變化均有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；P7-O1（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）；C3-P3（P < 0.01）；P3-O1（P < 0.01）。發作前期較發作期：P7-O1（P < 0.01）；C3-P3（MD = 6.55 × 10?5，P = 0.0494）]。尤其在FP1-F7、F7-T7、T7-P7和F3-C3電極通道中，發作后期的MI值較發作前期和發作期增加，差異具有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）]，而發作前期和發作期差異不具有統計學意義。但在P7-O1和C3-P3電極通道中，發作期MI值較發作前期和發作后期明顯下降，而在P3-O1電極通道中，僅發作期和發作后期MI值的差異具有統計學意義（P < 0.01）。
圖3
癲癇過程不同階段左半腦各通道內MI值統計變化
Figure3.
Statistical alterations of MI values of each channel signal recorded in the left hemisphere during different epileptic phases
與左半腦的統計分析類似，右半腦通道EEG信號β-γ耦合的變化趨勢也不同，如圖4所示。顯然，FP2-F4、F4-C4、P4-O2、FP2-F8和T8-P8電極通道的MI值變化趨勢較一致。盡管這些通道MI值在發作前期和發作期的差異沒有統計學意義，但可以觀察到發作后期MI值較發作前期和發作期增加，且差異具有統計學意義[（發作前期較發作后期：FP2-F4（MD = ?1.02 × 10?4，P = 0.0185）；F4-C4（P < 0.01）；P4-O2（MD = ?1.10 × 10?4，P = 0.0182）；FP2-F8（P < 0.01）；T8-P8（P < 0.01）。發作期較發作后期：FP2-F4（P < 0.01）；F4-C4（MD = ?1.09 × 10?4，P = 0.0199）； P4-O2（P < 0.01）；FP2-F8（P < 0.01）； T8-P8（P < 0.01）]。此外，與發作期相比，C4-P4、F8-T8、P8-O2和CZ-PZ電極通道的MI值在發作后期增加，差異具有統計學意義[C4-P4（P < 0.01）； F8-T8（P < 0.01）； P8-O2（P < 0.01）； CZ-PZ（P < 0.01）]。在這些通道EEG信號中，除了在F8-T8電極通道中觀察到差異具有統計學意義外[F8-T8（MD = ?9.44 × 10?5，P = 0.0115）]，C4-P4、P8-O2和CZ-PZ通道信號的MI值在發作前期和發作期的差異均無統計學意義。
圖4
癲癇過程不同階段右半腦各通道內MI值統計變化
Figure4.
Statistical alterations of MI values of each channel signal recorded in the right hemisphere during different epileptic phases
進一步，為探究癲癇過程不同階段β-γ耦合的變化，本文依據β-γ耦合統計分析的結果，將18導EEG數據對應的電極通道進行了歸類。如表1所示，所有電極通道被分成4類，對應著癲癇過程不同階段β-γ耦合的4種改變。這些結果間接指出，各腦區的信息交互在癲癇過程不同階段具有明顯的變化，意味著β-γ耦合的改變可能反映了腦區功能的變化，為定位致癇灶提供了潛在參考。
表1
癲癇過程EEG數據的β-γ耦合分類
Table1.
Statistical classification of EEG β-γ coupling during epilepsy
需要指出的是，盡管部分通道信號的β-γ耦合變化趨勢不同，但左右半腦中也存在變化趨勢一致的電極通道，表明這些電極通道對應的腦區在功能表現上可能較相似。有趣的是，雖然每個電極通道呈現出的變化趨勢略有差異，但在癲癇過程中，整個皮層的β-γ耦合在癲癇的三個狀態下表現出的差異均具有統計學意義（P < 0.01），如圖5所示。特別是，β-γ耦合的MI值在發作后期增加，其差異具有統計學意義（P < 0.01），與部分通道EEG信號MI值的統計結果一致。
圖5
癲癇過程不同階段全腦所有通道信號內MI值統計變化
Figure5.
Statistical alterations of MI values of all channel signals recorded in the whole brain
2.3 癲癇過程不同階段通道間β-γ耦合變化
上述研究結果指出，癲癇過程不同階段通道內β-γ耦合的差異具有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）；FP2-F4（MD = ?1.02 × 10?4，P = 0.0185）；F4-C4（P < 0.01）；P4-O2（MD = ?1.10 × 10?4，P = 0.0182）；FP2-F8（P < 0.01）；T8-P8（P < 0.01）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；P7-O1（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）；C3-P3（P < 0.01）；P3-O1（P < 0.01）；FP2-F4（P < 0.01）；F4-C4（MD = ?1.09 × 10?4， P = 0.0199）；C4-P4（P < 0.01）；P4-O2（P < 0.01）；FP2-F8（P < 0.01）；F8-T8（P < 0.01）；T8-P8（P < 0.01）；P8-O2（P < 0.01）；CZ-PZ（P < 0.01）。發作前期較發作期：P7-O1（P < 0.01）；C3-P3（MD = 6.55 × 10?5，P = 0.0494）；F8-T8（MD = 9.44 × 10?5，P = 0.0115）]。進一步，本文將探究不同通道EEG信號間β-γ耦合的存在性。需要說明的是，為減少重復的統計對比分析，本文根據表1的統計歸類，從發作期較發作前、后期MI值差異具有統計學意義的3類統計結果中，分別隨機選取1導EEG數據作為交叉PAC分析的EEG信號。如圖6和圖7所示，本文選取了3導電極名稱為FP1-F7、T7-P7和P3-O1的EEG信號，分別計算和統計分析了3導EEG數據間的β-γ耦合。需要指出的是，從通道FP1-F7、T7-P7和P3-O1濾波得到的相位頻率分別記為FP1-F7(P)、T7-P7(P)和 P3-O1(P)，從通道T7-P7、FP1-F7和P3-O1濾波得到的幅值頻率分別記為T7-P7(A)、FP1-F7(A)和P3-O1(A)。
圖6
癲癇過程不同通道信號間β-γ耦合變化
Figure6.
Alterations of β-γ coupling between different channel signals in epileptic phases
圖7
癲癇過程不同通道信號MI值統計比較
Figure7.
Comparisons of MI between different channel signals in epileptic phases
類似于同通道內的結果，本文發現癲癇過程不同通道EEG信號間存在β-γ耦合，且耦合強度在發作后期有增強趨勢。其中，在3導EEG信號通道間的PAC分析中，FP1-F7(A)與P3-O1(P)和T7-P7(P)間的PAC在癲癇過程中呈現一致變化，即發作前期較發作后期和發作期的MI值降低，差異具有統計學意義[發作前期較發作后期： T7-P7(P)和FP1-F7(A)（P < 0.01）; P3-O1(P)和FP1-F7(A)（P < 0.01）。發作前期較發作期：T7-P7(P)和FP1-F7(A)（P < 0.01）; P3-O1(P)和FP1-F7(A)（MD = ?2.80 × 10?4，P = 0.0248）]。類似地，P3-O1(A)分別與T7-P7(P)、FP1-F7(P)間的PAC在發作后期較發作前期MI值增大，其差異具有統計學意義[發作前期較發作后期： FP1-F7(P)和P3-O1(A) （MD = ?3.09 × 10?4，P = 0.0153）;T7-P7(P)和P3-O1(A) （MD = ?2.88 × 10?4，P = 0.0192）]。盡管各通道間PAC差異不同，但統計分析結果指出，發作后期較發作前期MI值增加，差異具有統計學意義[發作前期較發作后期：FP1-F7(P)和T7-P7(A)（P < 0.01）；T7-P7(P)和FP1-F7(A)（P < 0.01）；FP1-F7(P)和P3-O1(A) （MD = ?3.09 × 10?4，P = 0.0153）; P3-O1(P)和FP1-F7(A)（P < 0.01）;T7-P7(P)和P3-O1(A) （MD = ?2.88 × 10?4，P = 0.0192）; P3-O1(P)和T7-P7(P) （P < 0.01）]。這些結果表明腦區間的信息交互在癲癇過程中隨癲癇狀態發生改變，且在發作后期腦區間的信息交互有可能增強。
3 討論
本文針對難治性癲癇患者的EEG數據，采用MI的PAC方法統計分析發現，癲癇過程中各腦區存在β-γ耦合，且癲癇過程不同階段的β-γ耦合強度存在差異。特別是在癲癇發作后期，多數腦區的PAC較發作前期和發作期增大，差異具有統計學意義，表明癲癇發作后，各腦區之間的信息交互可能加強。這些研究發現為理解癲癇狀態下腦區間的信息通信機制提供了新視角，有助于揭示癲癇過程中大腦認知功能的變化和臨床癲癇的診療。
癲癇過程不同階段β-γ耦合強度的差異間接反映了大腦認知狀態和信息交互的改變。早期的研究表明，PAC能夠刻畫大腦內信息交互[9, 21-23]。其中，θ-γ耦合在大腦神經節律PAC研究中最普遍。諸如，基于小鼠的研究表明，海馬θ-γ耦合可反映決策中的信息處理[24]。同樣，在人體研究中發現，人類額葉與額葉-顳葉之間的θ-γ耦合變化和聯想記憶任務中視覺感知相關[25]。不同于θ-γ耦合，盡管β-γ耦合的研究較少，但有研究指出視覺工作記憶中，大腦神經節律β-γ耦合強度會隨認知任務而改變[26]，意味著β-γ耦合與腦功能認知狀態間存在關聯。特別是，基于癲癇患者的EEG數據，本文發現癲癇過程中β-γ耦合的存在，且β-γ耦合的強度會隨癲癇狀態發生改變，意味著腦區間的信息交互可能會受到癲癇影響，這為解釋癲癇患者認知功能下降提供了參考。事實上，基于難治性癲癇患者顱內EEG數據的研究指出，發作期β-γ耦合強度增加[27]。這一結果雖與本文在EEG信號部分電極通道中觀測到的不一致，如P7-O1、F8-T8和C3-P3等，但這些不同可能與患者個體差異，如患者的病程及腦認知狀態等有關。因此，β-γ耦合與大腦的認知功能有密切聯系，其變化能間接反映腦認知狀態及腦區間信息交互的改變。
β-γ耦合為臨床精準定位致癇灶提供了新的可能。在難治性癲癇手術治療過程中，準確定位致癇灶是手術成功的關鍵。早期曾有研究嘗試應用β/γ活動探究癲癇發作起源。例如，Guggisberg等[28]使用腦磁圖和自適應空間濾波技術定位β和γ波段的棘波鎖定功率增加的來源，指出通過棘波鎖定β/γ功率變化的成像可以定位致癇組織。Lundstrom等[29]使用β-γ活動作為EEG信號的生物標志物，準確定位了內側顳葉癲癇的發作起始區域。與這些致癇灶定位研究不同，本文通過β-γ耦合將全腦的EEG信號分為四類，每類分別對應癲癇狀態中β-γ耦合的統計變化特征，如表1所示。雖然這些分類無法直接對致癇灶進行準確定位，但是它們可提供重要的參考信息以確定癲癇致癇灶的范圍。這主要是由于癲癇過程不同階段，腦區間的信息交互會發生變化，而β-γ耦合的改變能夠揭示這些腦區間信息交互的改變。因此，本文研究結果表明β-γ耦合有助于癲癇的致癇灶定位。
需要指出的是，本文研究結果雖對理解癲癇過程中腦區間信息交互的機制和臨床難治性癲癇的有效治療具有參考意義，但仍存在部分局限性。第一，CHB-MIT數據集中每位患者的臨床參考信息有限，致使在應用β-γ耦合指導致癇灶定位方面缺乏臨床參考。第二，本文僅重點關注β-γ耦合在癲癇過程中的統計變化，未對其它潛在頻段間的耦合進行全面探究。這些潛在頻段間的耦合可能為癲癇研究提供多角度信息。第三，由于EEG數據的自身局限性和PAC生理機制的不明確性，癲癇過程中β-γ耦合改變的腦機制仍需要進行深入探究。未來，本課題組將與臨床緊密合作，將臨床信息與數據分析結合，深入探究癲癇過程中多頻段的交叉耦合，為癲癇診療提供服務。
4 結論
本文基于難治性癲癇患者的EEG數據，采用MI的PAC分析方法，系統探究了癲癇過程β-γ耦合的變化。研究發現癲癇過程中β-γ耦合的改變與癲癇狀態密切相關，特別在癲癇發作后期，多數腦區的β-γ耦合增強，其差異具有統計學意義，表明大腦內的神經信息交互在癲癇發作后有所加強。這些研究結果為揭示癲癇不同狀態下腦區間的信息交互提供了新視角，為臨床癲癇診療提供了理論參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李凱杰主要負責數據處理、數據分析與論文寫作；盧俊峰主要負責繪圖；余仁萍、張銳主要負責論文審閱；陳明明主要負責實驗指導，數據分析指導和論文撰寫指導。
0 引言
癲癇（epilepsy，EP）是一種由神經元異常放電引起短暫大腦功能障礙的神經系統疾病[1]。反復發作是癲癇主要的臨床表現，嚴重影響患者身心健康，甚至造成患者中樞神經系統損傷[2]。當前，盡管癲癇臨床診療水平較高，多數癲癇發作能夠得到及時有效的控制，但全球范圍內仍有約30%的癲癇患者對抗癲癇藥物等不敏感，發展成為難治性癲癇[3]。而且，致癇因素的多樣性和癲癇過程中復雜的腦內信息交互嚴重阻礙了臨床準確評估和有效治療難治性癲癇。因此，揭示癲癇過程中腦內信息交互隨癲癇狀態的改變，對臨床認識和診療癲癇具有重要的指導和參考意義。
腦電圖（electroencephalography，EEG）監測技術是通過置于頭皮上的電極獲取大腦內同步神經元活動產生的EEG信號，在臨床癲癇診療中應用廣泛[4]。癲癇EEG信號含有豐富的生理和病理信息，在癲癇發作全程常伴有明顯的特征變化，而解析癲癇過程EEG信號的動態變化特征有助于預測癲癇發作、探究癲癇播散過程和揭示致癇機制等。諸如，Das等[5]基于癲癇EEG信號，采用特征提取、模式匹配等算法組合而成的模型，能夠實現對癲癇發作的準確預測。Lopes等[6]應用鎖相因子構建功能網絡，提出發作生成傳播（ictogenic spread，IS）概念，并基于全面性癲癇患者和內側顳葉癲癇患者的EEG信號，揭示癲癇發作病理性活動在整個大腦網絡中的傳播特性。Mormann等[7]應用平均相位相干性統計度量相位同步差異，表征大腦動力學狀態變化，指出癲癇患者EEG信號的相位同步性變化可能有助于深入了解癲癇發作全過程。事實上，癲癇EEG信號除相位改變之外，幅度也會發生不同程度的改變。這些相位和幅度的改變與癲癇過程中腦內信息交互的轉變可能有密切關系。
相位—幅度耦合（phase-amplitude coupling，PAC）通過計算EEG信號的低頻相位與高頻幅度間的耦合關系，反映了腦區間的信息交互，在理解神經元網絡通信和連接機制等方面具有重要作用[8-9]。諸如在癲癇研究中，PAC可呈現不同腦區信息交互上的差異性，展現了腦內信息交互的特異性改變，是輔助診療癲癇的潛在生物標記物。近期， Ma等[10]在額葉癲癇的顱內EEG信號中發現，發作前期和間期的PAC可作為定位致癇灶的潛在標記物，且耦合的首選相位隨癲癇發作規律性改變 。類似的，基于小鼠的實驗研究發現，小鼠癲癇發作期皮質耦合受損，弱化的θ-γ耦合可能是癲癇發作風險預測的早期指標[11]。另外，在帕金森病的研究中，Salimpour等[12]發現皮層的相位靶刺激調節運動皮層中異常PAC，指出PAC在神經疾病治療上具有參考作用。
然而，需要注意的是，癲癇研究中PAC的分析多集中于θ-γ耦合，對β-γ耦合在癲癇過程中的變化知之甚少[13-14]。為此，本文將基于難治性癲癇患者的EEG信號，采用調節指數（modulation index，MI）的PAC分析方法，探究癲癇過程中β-γ耦合的存在性及隨癲癇狀態的變化。通過統計比較癲癇過程通道內和通道間的β-γ耦合強度變化，探究癲癇過程中腦內信息交互改變。通過本文研究，期望可為理解癲癇不同階段腦區間的信息交互提供新思路，為臨床診療癲癇提供理論借鑒。
1 材料和方法
1.1 數據采集及預處理
本文使用的EEG數據來源于波士頓兒童醫院（Children’s Hospital Boston，CHB）和麻省理工學院（Massachusetts Institute of Technology，MIT）創建的公開數據庫CHB-MIT頭皮EEG數據庫（scalp EEG Database）[15]。
由于數據集中每個案例的發作次數、單次記錄時長等均存在差異，為使數據在分析中保持一致性，本文在選取數據時，要求選取需滿足單次發作時長大于5 s，且兩次發作間隔超過1 h的數據，故選取了22個滿足上述條件的案例。針對每個案例，本文選取時長為6 h的數據，共計65次發作。同時，為避免不同案例電極數量、電極分布位置和壞導等因素對結果的影響，本文選取前18個通道進行分析。另外，為避免噪聲、偽跡（眼電、心電和肌電信號）等干擾，本文對上述選取的原始EEG數據進行了預處理。具體流程有：①采用平均參考方法對原始EEG數據進行重參考處理。②使用陷波濾波器濾除58～62 Hz之間的信號，消除60 Hz工頻干擾。③應用獨立成分分析（independent component analysis，ICA）去除EEG信號中的噪聲偽跡[16-17]。
最后，對預處理后的EEG數據，本文以數據集標記的發作時間為準，進行了分段截取，如圖1所示。假設發作期時長為L s，本文中定義發作前期為發作前3 × L s（3L），發作后期為發作后3 × L s（3L）。
圖1
數據選擇示意圖
Figure1.
Data selection diagram
1.2 方法
1.2.1 計算方法
本文參考文獻[18]的方法，采用MI 計算PAC。MI度量的是信號高頻幅度在低頻相位上的分布相對于均勻分布的偏差大小，是通過巴克-萊布勒（Kullback-Leibler）距離實現[19]。MI方法具有靈敏度高，對數據長度、信噪比和采樣率等影響因素魯棒性強[20]的特點，其具體計算過程如下：
（1）假設有原始信號x(t)，對其進行帶通濾波，可得到低頻段和高頻段信號，分別記為xfP(t)和xfA(t)；
（2）采用希爾伯特變換分別計算得到低頻段信號xfP(t)的相位時間序列φfP(t)和高頻段信號xfA(t)的幅度包絡AfA(t)；
（3）對相位時間序列φfP(t)進行均勻區間劃分，計算每個區間相位對應幅度包絡AfA(t)的平均值AfAmean；
（4）歸一化AfAmean，如式(1)所示：
|
其中，AfAmean(j)是第j個區間相位對應的幅度包絡均值；N是區間個數，通常N = 18；P(j)是第j個區間相位標準化的幅度包絡均值；k為第k個區間相位。
（5）計算Kullback-Leibler距離：假設存在離散分布P和Q，其對應的Kullback-Leibler距離如式(2)所示：
|
其中，DKL(P, Q) ≥ 0，當且僅當P = Q時，DKL(P, Q) = 0；P(j)是離散分布P第j個區間的概率密度；Q(j)是離散分布Q第j個區間的概率密度。式(2)中，P分布的香農熵H(P)定義如式(3)所示：
|
因此，結合式(2)和式(3)可得出Kullback-Leibler距離與香農熵的關系，如式(4)所示：
|
其中，U是均勻分布，log(N)是最大熵值。
（6）綜合上述公式，定義MI（以符號MI表示），如式(5)所示：
|
由MI的定義可知，MI的值在0～1之間。當高頻幅值均勻分布在低頻相位上時，MI為0，即滿足j = 1, ???, N和P( j) = 1/N，在其它情況下，MI不為0。MI的大小與離散分布P和均勻分布U的偏離程度呈正相關。在MI計算中，本文設定低頻段相位頻率為13～30 Hz，高頻段幅值頻率范圍為30～100 Hz，步長為1 Hz。
1.2.2 統計分析
本文將發作期數據劃分為5 s/段，然后分別計算每段數據β、γ頻率范圍內MI值的平均值，得到一組β-γ耦合值，并以該組β-γ耦合值的平均值作為此次癲癇發作期的MI值。其中，若存在數據段對應的β、γ頻率范圍內MI值均為0，則將該數據段對應的β-γ耦合值置為0；然后將發作前期和發作后期的數據進行類似的統計分析；最終發作前期、發作期和發作后期，各自分別獲得65個平均MI值，共可獲得195個平均MI值。由于三種癲癇狀態對應的MI值序列均不符合正態分布，故本文采用非參數統計秩和檢驗進行對比分析，即三組數據兩兩進行秩和檢驗，檢驗水準為0.05。
另外，在計算MI的過程中，為避免數據隨機性引入的偽MI干擾，本文采用了200次的隨機替代數據進行統計比較，當兩者的差異應具有統計學意義時，保留原始數據計算的MI值；否則，令MI值為0。其中，檢驗水準為0.05。
2 實驗結果
2.1 癲癇過程不同階段通道內β-γ耦合
早期研究表明，癲癇患者EEG信號的低頻（θ頻段）相位和高頻（γ頻段）振幅間存在PAC（θ-γ耦合），且θ-γ耦合對臨床癲癇診療有重要參考作用[1]。相對而言，癲癇過程中β-γ耦合是否存在及如何變化仍不清楚。為此，本文隨機選取了三個標記為FP1-F7、T7-P7和P3-O1的通道信號，這些EEG信號分別位于大腦皮層不同區域。如圖2所示，在發作前期、發作期和發作后期，每個通道信號中都可觀察到β-γ耦合。雖然三個通道信號的MI值不同，但這些MI值在發作演化過程中的變化趨勢幾乎是一致的，尤其發作后期MI值明顯增加。值得注意的是，β-γ耦合可以出現在振幅的不同中心頻率上。這些觀測結果表明，β-γ耦合不僅存在于不同的癲癇階段，而且隨癲癇狀態的演化而改變，這可為癲癇發作預測提供新的視角。
圖2
癲癇過程不同階段β-γ耦合變化
Figure2.
Alterations of β-γ coupling in different epileptic phases
2.2 癲癇過程不同階段通道內β-γ耦合的統計比較
如前所述，量化癲癇β-γ耦合在癲癇過程不同階段的變化，可為理解癲癇發作的動態演化提供新視角。因此，本文計算了每個通道內β-γ耦合的MI值，并對發作前期、發作期和發作后期三個狀態的MI值進行統計分析比較。如圖3所示，左半腦通道EEG信號在三個時期的MI值變化趨勢不同。事實上，除了FP1-F3和FZ-CZ電極通道的差異不具有統計學意義外，左半腦其它通道的EEG信號MI值在癲癇發作期的變化均有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；P7-O1（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）；C3-P3（P < 0.01）；P3-O1（P < 0.01）。發作前期較發作期：P7-O1（P < 0.01）；C3-P3（MD = 6.55 × 10?5，P = 0.0494）]。尤其在FP1-F7、F7-T7、T7-P7和F3-C3電極通道中，發作后期的MI值較發作前期和發作期增加，差異具有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）]，而發作前期和發作期差異不具有統計學意義。但在P7-O1和C3-P3電極通道中，發作期MI值較發作前期和發作后期明顯下降，而在P3-O1電極通道中，僅發作期和發作后期MI值的差異具有統計學意義（P < 0.01）。
圖3
癲癇過程不同階段左半腦各通道內MI值統計變化
Figure3.
Statistical alterations of MI values of each channel signal recorded in the left hemisphere during different epileptic phases
與左半腦的統計分析類似，右半腦通道EEG信號β-γ耦合的變化趨勢也不同，如圖4所示。顯然，FP2-F4、F4-C4、P4-O2、FP2-F8和T8-P8電極通道的MI值變化趨勢較一致。盡管這些通道MI值在發作前期和發作期的差異沒有統計學意義，但可以觀察到發作后期MI值較發作前期和發作期增加，且差異具有統計學意義[（發作前期較發作后期：FP2-F4（MD = ?1.02 × 10?4，P = 0.0185）；F4-C4（P < 0.01）；P4-O2（MD = ?1.10 × 10?4，P = 0.0182）；FP2-F8（P < 0.01）；T8-P8（P < 0.01）。發作期較發作后期：FP2-F4（P < 0.01）；F4-C4（MD = ?1.09 × 10?4，P = 0.0199）； P4-O2（P < 0.01）；FP2-F8（P < 0.01）； T8-P8（P < 0.01）]。此外，與發作期相比，C4-P4、F8-T8、P8-O2和CZ-PZ電極通道的MI值在發作后期增加，差異具有統計學意義[C4-P4（P < 0.01）； F8-T8（P < 0.01）； P8-O2（P < 0.01）； CZ-PZ（P < 0.01）]。在這些通道EEG信號中，除了在F8-T8電極通道中觀察到差異具有統計學意義外[F8-T8（MD = ?9.44 × 10?5，P = 0.0115）]，C4-P4、P8-O2和CZ-PZ通道信號的MI值在發作前期和發作期的差異均無統計學意義。
圖4
癲癇過程不同階段右半腦各通道內MI值統計變化
Figure4.
Statistical alterations of MI values of each channel signal recorded in the right hemisphere during different epileptic phases
進一步，為探究癲癇過程不同階段β-γ耦合的變化，本文依據β-γ耦合統計分析的結果，將18導EEG數據對應的電極通道進行了歸類。如表1所示，所有電極通道被分成4類，對應著癲癇過程不同階段β-γ耦合的4種改變。這些結果間接指出，各腦區的信息交互在癲癇過程不同階段具有明顯的變化，意味著β-γ耦合的改變可能反映了腦區功能的變化，為定位致癇灶提供了潛在參考。
表1
癲癇過程EEG數據的β-γ耦合分類
Table1.
Statistical classification of EEG β-γ coupling during epilepsy
需要指出的是，盡管部分通道信號的β-γ耦合變化趨勢不同，但左右半腦中也存在變化趨勢一致的電極通道，表明這些電極通道對應的腦區在功能表現上可能較相似。有趣的是，雖然每個電極通道呈現出的變化趨勢略有差異，但在癲癇過程中，整個皮層的β-γ耦合在癲癇的三個狀態下表現出的差異均具有統計學意義（P < 0.01），如圖5所示。特別是，β-γ耦合的MI值在發作后期增加，其差異具有統計學意義（P < 0.01），與部分通道EEG信號MI值的統計結果一致。
圖5
癲癇過程不同階段全腦所有通道信號內MI值統計變化
Figure5.
Statistical alterations of MI values of all channel signals recorded in the whole brain
2.3 癲癇過程不同階段通道間β-γ耦合變化
上述研究結果指出，癲癇過程不同階段通道內β-γ耦合的差異具有統計學意義[發作前期較發作后期：FP1-F7（MD = ?8.99 × 10?5，P = 0.0220）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；F3-C3（MD = ?1.49 × 10?4，P = 0.0453）；FP2-F4（MD = ?1.02 × 10?4，P = 0.0185）；F4-C4（P < 0.01）；P4-O2（MD = ?1.10 × 10?4，P = 0.0182）；FP2-F8（P < 0.01）；T8-P8（P < 0.01）。發作期較發作后期：FP1-F7（MD = ?1.85 × 10?4，P = 0.0173）；F7-T7（P < 0.01）；T7-P7（P < 0.01）；P7-O1（P < 0.01）；F3-C3（MD = ?1.27 × 10?4，P = 0.0236）；C3-P3（P < 0.01）；P3-O1（P < 0.01）；FP2-F4（P < 0.01）；F4-C4（MD = ?1.09 × 10?4， P = 0.0199）；C4-P4（P < 0.01）；P4-O2（P < 0.01）；FP2-F8（P < 0.01）；F8-T8（P < 0.01）；T8-P8（P < 0.01）；P8-O2（P < 0.01）；CZ-PZ（P < 0.01）。發作前期較發作期：P7-O1（P < 0.01）；C3-P3（MD = 6.55 × 10?5，P = 0.0494）；F8-T8（MD = 9.44 × 10?5，P = 0.0115）]。進一步，本文將探究不同通道EEG信號間β-γ耦合的存在性。需要說明的是，為減少重復的統計對比分析，本文根據表1的統計歸類，從發作期較發作前、后期MI值差異具有統計學意義的3類統計結果中，分別隨機選取1導EEG數據作為交叉PAC分析的EEG信號。如圖6和圖7所示，本文選取了3導電極名稱為FP1-F7、T7-P7和P3-O1的EEG信號，分別計算和統計分析了3導EEG數據間的β-γ耦合。需要指出的是，從通道FP1-F7、T7-P7和P3-O1濾波得到的相位頻率分別記為FP1-F7(P)、T7-P7(P)和 P3-O1(P)，從通道T7-P7、FP1-F7和P3-O1濾波得到的幅值頻率分別記為T7-P7(A)、FP1-F7(A)和P3-O1(A)。
圖6
癲癇過程不同通道信號間β-γ耦合變化
Figure6.
Alterations of β-γ coupling between different channel signals in epileptic phases
圖7
癲癇過程不同通道信號MI值統計比較
Figure7.
Comparisons of MI between different channel signals in epileptic phases
類似于同通道內的結果，本文發現癲癇過程不同通道EEG信號間存在β-γ耦合，且耦合強度在發作后期有增強趨勢。其中，在3導EEG信號通道間的PAC分析中，FP1-F7(A)與P3-O1(P)和T7-P7(P)間的PAC在癲癇過程中呈現一致變化，即發作前期較發作后期和發作期的MI值降低，差異具有統計學意義[發作前期較發作后期： T7-P7(P)和FP1-F7(A)（P < 0.01）; P3-O1(P)和FP1-F7(A)（P < 0.01）。發作前期較發作期：T7-P7(P)和FP1-F7(A)（P < 0.01）; P3-O1(P)和FP1-F7(A)（MD = ?2.80 × 10?4，P = 0.0248）]。類似地，P3-O1(A)分別與T7-P7(P)、FP1-F7(P)間的PAC在發作后期較發作前期MI值增大，其差異具有統計學意義[發作前期較發作后期： FP1-F7(P)和P3-O1(A) （MD = ?3.09 × 10?4，P = 0.0153）;T7-P7(P)和P3-O1(A) （MD = ?2.88 × 10?4，P = 0.0192）]。盡管各通道間PAC差異不同，但統計分析結果指出，發作后期較發作前期MI值增加，差異具有統計學意義[發作前期較發作后期：FP1-F7(P)和T7-P7(A)（P < 0.01）；T7-P7(P)和FP1-F7(A)（P < 0.01）；FP1-F7(P)和P3-O1(A) （MD = ?3.09 × 10?4，P = 0.0153）; P3-O1(P)和FP1-F7(A)（P < 0.01）;T7-P7(P)和P3-O1(A) （MD = ?2.88 × 10?4，P = 0.0192）; P3-O1(P)和T7-P7(P) （P < 0.01）]。這些結果表明腦區間的信息交互在癲癇過程中隨癲癇狀態發生改變，且在發作后期腦區間的信息交互有可能增強。
3 討論
本文針對難治性癲癇患者的EEG數據，采用MI的PAC方法統計分析發現，癲癇過程中各腦區存在β-γ耦合，且癲癇過程不同階段的β-γ耦合強度存在差異。特別是在癲癇發作后期，多數腦區的PAC較發作前期和發作期增大，差異具有統計學意義，表明癲癇發作后，各腦區之間的信息交互可能加強。這些研究發現為理解癲癇狀態下腦區間的信息通信機制提供了新視角，有助于揭示癲癇過程中大腦認知功能的變化和臨床癲癇的診療。
癲癇過程不同階段β-γ耦合強度的差異間接反映了大腦認知狀態和信息交互的改變。早期的研究表明，PAC能夠刻畫大腦內信息交互[9, 21-23]。其中，θ-γ耦合在大腦神經節律PAC研究中最普遍。諸如，基于小鼠的研究表明，海馬θ-γ耦合可反映決策中的信息處理[24]。同樣，在人體研究中發現，人類額葉與額葉-顳葉之間的θ-γ耦合變化和聯想記憶任務中視覺感知相關[25]。不同于θ-γ耦合，盡管β-γ耦合的研究較少，但有研究指出視覺工作記憶中，大腦神經節律β-γ耦合強度會隨認知任務而改變[26]，意味著β-γ耦合與腦功能認知狀態間存在關聯。特別是，基于癲癇患者的EEG數據，本文發現癲癇過程中β-γ耦合的存在，且β-γ耦合的強度會隨癲癇狀態發生改變，意味著腦區間的信息交互可能會受到癲癇影響，這為解釋癲癇患者認知功能下降提供了參考。事實上，基于難治性癲癇患者顱內EEG數據的研究指出，發作期β-γ耦合強度增加[27]。這一結果雖與本文在EEG信號部分電極通道中觀測到的不一致，如P7-O1、F8-T8和C3-P3等，但這些不同可能與患者個體差異，如患者的病程及腦認知狀態等有關。因此，β-γ耦合與大腦的認知功能有密切聯系，其變化能間接反映腦認知狀態及腦區間信息交互的改變。
β-γ耦合為臨床精準定位致癇灶提供了新的可能。在難治性癲癇手術治療過程中，準確定位致癇灶是手術成功的關鍵。早期曾有研究嘗試應用β/γ活動探究癲癇發作起源。例如，Guggisberg等[28]使用腦磁圖和自適應空間濾波技術定位β和γ波段的棘波鎖定功率增加的來源，指出通過棘波鎖定β/γ功率變化的成像可以定位致癇組織。Lundstrom等[29]使用β-γ活動作為EEG信號的生物標志物，準確定位了內側顳葉癲癇的發作起始區域。與這些致癇灶定位研究不同，本文通過β-γ耦合將全腦的EEG信號分為四類，每類分別對應癲癇狀態中β-γ耦合的統計變化特征，如表1所示。雖然這些分類無法直接對致癇灶進行準確定位，但是它們可提供重要的參考信息以確定癲癇致癇灶的范圍。這主要是由于癲癇過程不同階段，腦區間的信息交互會發生變化，而β-γ耦合的改變能夠揭示這些腦區間信息交互的改變。因此，本文研究結果表明β-γ耦合有助于癲癇的致癇灶定位。
需要指出的是，本文研究結果雖對理解癲癇過程中腦區間信息交互的機制和臨床難治性癲癇的有效治療具有參考意義，但仍存在部分局限性。第一，CHB-MIT數據集中每位患者的臨床參考信息有限，致使在應用β-γ耦合指導致癇灶定位方面缺乏臨床參考。第二，本文僅重點關注β-γ耦合在癲癇過程中的統計變化，未對其它潛在頻段間的耦合進行全面探究。這些潛在頻段間的耦合可能為癲癇研究提供多角度信息。第三，由于EEG數據的自身局限性和PAC生理機制的不明確性，癲癇過程中β-γ耦合改變的腦機制仍需要進行深入探究。未來，本課題組將與臨床緊密合作，將臨床信息與數據分析結合，深入探究癲癇過程中多頻段的交叉耦合，為癲癇診療提供服務。
4 結論
本文基于難治性癲癇患者的EEG數據，采用MI的PAC分析方法，系統探究了癲癇過程β-γ耦合的變化。研究發現癲癇過程中β-γ耦合的改變與癲癇狀態密切相關，特別在癲癇發作后期，多數腦區的β-γ耦合增強，其差異具有統計學意義，表明大腦內的神經信息交互在癲癇發作后有所加強。這些研究結果為揭示癲癇不同狀態下腦區間的信息交互提供了新視角，為臨床癲癇診療提供了理論參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：李凱杰主要負責數據處理、數據分析與論文寫作；盧俊峰主要負責繪圖；余仁萍、張銳主要負責論文審閱；陳明明主要負責實驗指導，數據分析指導和論文撰寫指導。
