目前，帕金森病發生率在逐漸上升，嚴重影響著患者生活質量，社會診療負擔不斷加重。然而，該病的早期監測手段有限，很難及時干預。為了發現其標志物，本文對服藥前后帕金森病患者及健康人的32通道靜息態腦電數據進行分頻段研究。首先利用動態加權符號互信息計算各通道腦電信號間相關性，再通過k均值聚類實現信號關聯矩陣的分類，最后得到腦電信號關聯狀態。通過統計分析發現，在Beta頻段（P = 0.034）與Gamma頻段（P = 0.010）各有一個腦電信號關聯狀態可顯著區分未服藥帕金森病患者與健康人。這表明未服藥帕金森病患者與健康人的靜息態腦電各通道信號相關性差異有統計學意義。而在服藥與未服藥帕金森病患者、服藥帕金森病患者與健康人之間，其關聯狀態差異均沒有統計學意義。這可為帕金森病的臨床診斷提供一種參考。
引用本文： 丁昊, 吳進輝, 唐旭東, 余江南, 陳軒恒, 吳占雄. 基于動態加權符號互信息與k均值聚類的帕金森病患者靜息腦電關聯狀態識別. 生物醫學工程學雜志, 2023, 40(1): 20-26. doi: 10.7507/1001-5515.202211002 復制
0 引言
帕金森病（Parkinson’s disease，PD）是一種常見的慢性行進性神經系統疾病，多見于老年人，給患者及其家庭帶來極大的痛苦與嚴重的經濟負擔。它主要是由黑質多巴胺能神經元的異常變化引起的[1-2]。帕金森病主要表現為運動緩慢、姿勢不穩、靜止性震顫，以及認知、睡眠、精神與自主神經功能障礙等非運動癥狀[3-5]。帕金森病的治療手段主要有藥物、手術、神經康復等，也可以通過增加多巴胺和神經元保護以延緩疾病進展。然而，其確切發病機制目前尚不清楚，也沒有藥物可以逆轉病程以實現根治[6]。發現簡單易用的帕金森病生物學標志，有助于及早發現、干預疾病，將為患者帶來福音。
腦電信號（electroencephalogram，EEG）利用非侵入式電極記錄腦神經元活動情況，能夠反映大腦皮層的活動規律，可為帕金森病診斷提供生物學標志[7-8]。目前已經有大量研究表明，帕金森病患者相比于健康人在腦電圖上存在異常。多巴胺神經元的變性衰亡會導致基底節-丘腦-皮層回路的同步模式明顯異常，進一步影響帕金森病患者的腦網絡信息編碼[9]。電生理學研究表明，帕金森病表現為基底節-皮層環路、皮層-皮層環路的病理振蕩，其中基底節同步活動與皮層的活動呈線性相關[10]。通過對腦電信號進行分析，可以發現基底節-丘腦-皮層回路活動改變的標記，為了解其大腦皮層異常動態變化提供新的途徑[9]。帕金森病患者的枕部Alpha與Theta頻段腦電活動信號減弱，視覺空間障礙會導致枕部右側腦電信號異常[11]。未服藥帕金森病患者的Beta頻段腦電信號波形在感覺運動電極上的振蕩比服藥后有更大的銳度不對稱和陡度不對稱性[12]。處于靜息狀態時，腦電圖顯示帕金森病患者與健康人對照組相比在Theta和Alpha頻段上功率增加，Beta和Gamma頻段上功率降低[13]。在腦電信號頻譜圖中，帕金森病患者的Beta頻段波振幅比健康人小[14]，而Alpha波段腦電信號振幅比健康人大，且其頻率減慢[15]。Zhang等[16]通過將時頻分析與深度學習相結合，應用可調Q因子小波變換深度殘余收縮網絡來對腦電數據進行分類，帕金森病預測的準確率達到99.92%。在對高分辨腦電數據源定位后，Gerard等[17]采用動態鎖相方法來評估帕金森病患者功能連接性，發現患者額區與枕顳區顯示出異常強的Theta頻帶連接性。基于腦電數據5個頻帶的相干性，Conti等[18]通過計算連接矩陣及分析網絡拓撲參數，發現帕金森病患者在Alpha、Beta頻帶連通性較低，而在Gamma頻帶連接性增強。此外，各通道腦電信號并非獨立不相關，它們之間相互影響，可能存在著某種非線性依賴關系，蘊含著可以作為檢驗帕金森病的生物學標志。為此，加權符號互信息算法（weighted symbolic mutual information，wSMI）可以用于探索腦動態功能連接性及病理條件下的腦網絡潛在變化（如阿爾茨海默病、多發性硬化癥、精神分裂癥和社交焦慮障礙）[19]。有研究使用wSMI與聚類來研究腦連接的動態過程，并對比分析了精神分裂癥患者與健康人的丘腦區域網絡效率[20]。也有人采用wSMI研究腦電信號Delta、Theta頻帶的腦皮質功能連接性[21]。
基于以上思考，本文假設帕金森病患者與正常健康人的靜息態腦電信號在各個頻段內可能存在關聯狀態差異。在計算靜息態腦電信號關聯狀態后，通過比較帕金森病患者與健康人的各通道腦電信號間的相關性來發現生物學標志。腦電信號預處理后，利用wSMI計算通道間相關性以得到動態連接矩陣，再利用k均值聚類進行關聯狀態識別，并計算出不同對象的關聯狀態發生比例，最后統計發現與帕金森病顯著相關的功能連接狀態。這里的狀態是指不同通道腦電信號之間存在的神經活動同步性。
1 方法
1.1 腦電數據
本文所采用的腦電數據為加州大學圣地亞哥分校帕金森病患者靜息狀態腦電開放數據[22]。數據集包括15名帕金森病患者［8名女性，7名男性，年齡為（63.27±8.20歲）］和15名正常健康人［9名女性，6名男性，年齡為（63.50±9.67歲）］，所有人均為右利手。在采集腦電數據時，患者繼續藥物治療方案，治療藥物主要為能影響多巴胺能的左旋多巴等。通過平衡患者用藥和停藥的時間，分別收集患者按時用藥的腦電信號（PD-on）以及患者至少停止用藥12 h（PD-off）的腦電信號。同時，收集健康人（healthy control，HC）的腦電信號作為對比分析。腦電信號采集對象信息如表1所示。腦電信號采用32通道的BioSemi ActiveTwo系統采集，并有額外的八個電極放置在眼睛的外側和下方，以監測眼睛的眨眼和運動。采集時所有人均處于舒適的靜息狀態，采樣頻率為512 Hz，靜息數據連續記錄至少3 nub。這里值得注意的是，PD-off與PD-on腦電信號采集于同一組PD患者。
表1
研究對象信息
Table1.
Demographic information of the subjects
1.2 腦電信號預處理
將腦電信號導入到EEGLAB MATLAB工具箱中進行預處理，具體步驟如下：① 坐標定位：對電極進行坐標定位，只對前32個通道采集的信號進行分析，去除后8個不反映腦電信號的電極（EXG1，EXG2，…，EXG8）。② 重定義參考：取每個電極的平均值，應用全部腦電數據的共同平均值進行重定義參考。③ 濾波：使用四階巴特沃斯數字濾波器帶通濾波，分別得到每個對象腦電信號的Delta（0.5～4 Hz）、Theta（4～8 Hz）、Alpha（8～13 Hz）、Beta（13～30 Hz）、Gamma（30～50 Hz）及高頻High-frequency（50～100 Hz）等六個頻帶區間，過濾掉其他過高頻噪聲信號和過低頻的漂移偽跡。④ 獨立成分分析：通過獨立成分分析（independent component analysis，ICA），去除眼電等生理偽影的干擾，得到純凈的腦電信號。⑤ 截段提取：考慮到每個人在儀器測量的開始和結束可能會有心理狀態的波動，未能立刻進入靜息狀態，同時各個對象所測得的腦電信號時間及數據序列長度有所不同，為了得到較為穩定的腦電信息，以腦電信號時間序列中心為基準向前后兩側各取80 s，剔除了數據最前段和最后段，只保留中間有效部分，最終得到160 s的數據，使所有對象的腦電信號時間序列長度均相等。
1.3 動態wSMI分析
本文采用wSMI來評估腦電信號隨機非線性聯合波動的程度，度量時間序列間的共享信息[21]。同時，該算法能夠發現信號增加或減少的定性“符號”模式，在不需要對交互類型進行假設的情況下能對信號的熵進行快速和穩健估計。加權符號互信息通過計算不同時間序列間同時出現的每一對符號的聯合概率來估計兩個信號之間共享符號的互信息，并通過忽略相同或相反符號信號的共同出現以解決容積傳導的問題[19]，其權值計算能夠避免來自共同來源的腦電信號間的偽相關性，使得到的結論更為嚴格。這是一種檢測非線性耦合的有效方法[23]。其計算公式如下：
 |
其中
、
為時間序列，
為
中元素，
為
中元素。
為構成符號的元素個數，
為加權函數，
為概率函數，
為自然對數函數。
首先根據預先定義的連續采樣數據的振幅趨勢將腦電信號轉變為一系列離散符號。通過估計每對符號的聯合概率來計算腦電信號序列間的加權符號互信息值。對于一對相同的符號，權值設置為零，表示這對符號可以由唯一的共同源引出[9]。為了分析腦電信號時間序列的動態特性，不是一次性對整個時間序列進行加權符號互信息計算，而是通過施加時間窗口，每次取一段腦電信號進行處理計算。隨著時間窗口的滑動，可以得到一系列加權符號互信息計算結果，構成三維矩陣。這樣可以更好地體現出腦狀態的時間波動性。本文里，wSMI采用Python語言實現，底層調用MNE庫（MEG、EEG analysis & visualization，用來處理EEG/MEG等信號的開源程序，
1.4 k均值聚類
k均值聚類，是一種經典的高效快速聚類方法，其計算過程是把多個樣本分為k個簇，按照距離中心點最“小”的原則，使得各簇內的樣本點具有較高的相似性，而劃為不同簇的樣本相似性較低。相似度是依據一個簇中樣本的平均值來進行計算的。每次迭代后重新計算聚類中心，重新劃分類別，由此不斷迭代重復，直至聚類中心值趨于穩定[26]。
將加權符號互信息得到的關聯矩陣集合按照每個頻段分別聚類。這里聚類數目取k = 6[20]。對每個矩陣樣本計算它到聚類中心的距離并將它分到距離最小的類中，然后針對每個矩陣樣本類別重新計算聚類中心。兼顧速度和效率，這里設置重復迭代上限為200次。最后得出所有頻段的聚類中心、樣本所屬類別、概率分布和轉移概率等數據。最后得到的關聯矩陣如圖1所示。
圖1
腦電信號關聯矩陣
Figure1.
Correlation matrix between each pair of EEG channels
1.5 統計分析
不考慮其他因素的影響，把腦電數據分為健康人HC、未服藥帕金森病患者PD-off、服藥帕金森病患者PD-on三類（每組15人，共三組），作為ANOVA統計分組變量。在驗證所得的樣本數據符合正態分布后，計算HC、PD-off、PD-on各組對象在不同頻段關聯狀態占比，之后對這些樣本進行ANOVA統計分析。在發現不同組間差異有統計學意義后，采用最小顯著差別方法（least significant difference，LSD）進行校驗。檢驗水準為0.05。本論文所實現的算法源代碼與分析結果數據可從
2 結果
圖2展示了不同頻段靜息態腦電信號關聯矩陣聚類結果，每個頻帶共聚類為6個關聯狀態。本論文采用狀態所占比例作為指標來區分不同分組，探究將靜息態腦電信號關聯狀態作為帕金森病生物學標志的可能性。健康人、服藥帕金森病患者、未服藥帕金森病患者三個分組的腦電信號關聯狀態占比數據如表2所示。
圖2
對三組對象采用加權符號互信息及k均值聚類獲取靜息態腦電信號關聯狀態
每個頻段聚類得到6個關聯狀態，采用關聯狀態占比作為指標來區分不同組
Figure2. Correlation matrices extracted from EEG signals using wSMI and k-means clustering in three groupssix correlation estimates were obtained for each frequency band. The correlation status proportion was used as an index to distinguish different groups
表2
三組對象靜息態腦電信號關聯狀態占比（
±s）
Table2.
Correlation status proportion of HC, PD-off, and PD-on groups (
±s)
ANOVA統計分析之后，發現HC與PD-off組間Beta頻帶State1占比差異有統計學意義（P = 0.034，見表3），同時這兩組間Gamma頻帶State5占比差異也有統計學意義（P = 0.010，見表4）。這說明在帕金森病患者與健康人之間，靜息態腦電信號關聯性在Beta和Gamma頻帶顯著不同。
表3
Beta頻帶關聯狀態1占比的顯著性水平（LSD校正）
Table3.
Significance level of the proportion of state1 in Beta band (LSD test)
表4
Gamma頻帶關聯狀態State5占比的顯著性水平（LSD校正）
Table4.
Significance level of the proportion of state5 in Gamma band (LSD test)
3 討論與結論
本論文利用動態wSMI及k均值聚類實現了帕金森病患者靜息態腦電信號關聯性分析。發現在Beta段范圍內，健康人與未服藥帕金森病患者狀態State1占比差異有統計學意義（P = 0.034），而未服藥與服藥帕金森病患者之間State1占比差異沒有統計學意義（P = 0.437），健康人與服藥患者間State1占比差異也沒有統計學意義（P = 0.169）。在Gamma頻帶，健康人與未服藥帕金森病患者間State5占比差異有統計學意義（P = 0.010），而未服藥與服藥患者間State5占比差異沒有統計學意義（P = 0.560），健康人與服藥患者間State5占比差異也沒有統計學意義（P=0.053）。
考慮到k均值聚類可能存在誤差，本文進行了多次重復聚類計算，均得到同樣的結果。同時，觀察到Delta頻帶的狀態圖與其他頻段有明顯的不同（見圖2），猜想在Delta頻帶通道相關性可能具有與其他頻帶完全不同的獨特性質。此外，還發現帕金森病患者的腦電波形在后半段有更大的波動和陡峭性，而健康人的則平坦穩定，這與之前的研究結果是一致的[12]。在k均值聚類時，Delta與High-frequency頻段具有比較大的波動性，有可能是頻段過窄或過寬的原因。
另外，本文也驗證了當聚類數目分別設為4、5、7時，并沒有發現差異具有統計學意義的腦網絡關聯狀態，而只有在聚類數為6時發現兩個有顯著差異的狀態。這也驗證了Fu等[20]的工作。
同時，分析本文研究方法可能產生誤差的原因：① 獨立成分分析去除生理偽影。采用獨立成分分析去除生理偽影時，對生理偽影的判斷可能會有一定的誤差。② 高頻帶的噪聲干擾會對互信息計算有所影響。③ 采集靜息態數據時對象可能存在情緒或認知活動波動，例如患者腦神經功能的缺陷使得對象不能長時間集中精力等。
在Beta頻帶State1和Gamma頻帶State5占比上，盡管未服藥帕金森病患者與服藥帕金森病患者之間、服藥帕金森病患者與健康人之間差異沒有統計學意義，但還是可以觀察到患者服藥后其狀態占比的統計均值介于未服藥和健康人之間，說明服藥能夠縮小這兩者狀態的差距，因此服藥對于病變的程度還是有所改善的。關于帕金森病患者服藥與未服藥的差異，還可以在此基礎上進一步深入研究。同時，在本研究采用32通道腦電數據得出了上述結論，如果基于更多采集通道（如64、128等）的腦電信號數據，腦電信號關聯狀態的了解與分析可以更為透徹與精確。
wSMI能夠評估腦電信號間非隨機聯合波動的程度，度量其共享信息值。相較于其他相關性算法，加權符號互信息主要有三個優點：① 能夠尋找信號增加或減少的定性或“符號”模式，從而快速和準確地估計信號的熵，符號的轉換取決于符號的長度和時間間隔；② 提供了一種檢測非線性耦合的有效方法，并且不需要對相互作用的類型進行假設，提高了結果的真實性；③ 權重計算舍棄了腦電信號之間的偽相關性[21, 23]。這些特點使得wSMI非常合適腦電信號關聯狀態的識別。
總之，本研究采取動態加權符號互信息方法計算腦電信號間關聯性，通過聚類評估帕金森病患者腦電信號關聯狀態的改變模式，從而發現帕金森病的生物學標志。這為帕金森病患者的病情診斷與篩查提供了一種參考方法。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：丁昊負責數據分析和論文寫作，吳進輝與唐旭東負責腦電數據預處理，余江南與陳軒恒負責算法編程，吳占雄負責試驗設計與文字修改。
0 引言
帕金森病（Parkinson’s disease，PD）是一種常見的慢性行進性神經系統疾病，多見于老年人，給患者及其家庭帶來極大的痛苦與嚴重的經濟負擔。它主要是由黑質多巴胺能神經元的異常變化引起的[1-2]。帕金森病主要表現為運動緩慢、姿勢不穩、靜止性震顫，以及認知、睡眠、精神與自主神經功能障礙等非運動癥狀[3-5]。帕金森病的治療手段主要有藥物、手術、神經康復等，也可以通過增加多巴胺和神經元保護以延緩疾病進展。然而，其確切發病機制目前尚不清楚，也沒有藥物可以逆轉病程以實現根治[6]。發現簡單易用的帕金森病生物學標志，有助于及早發現、干預疾病，將為患者帶來福音。
腦電信號（electroencephalogram，EEG）利用非侵入式電極記錄腦神經元活動情況，能夠反映大腦皮層的活動規律，可為帕金森病診斷提供生物學標志[7-8]。目前已經有大量研究表明，帕金森病患者相比于健康人在腦電圖上存在異常。多巴胺神經元的變性衰亡會導致基底節-丘腦-皮層回路的同步模式明顯異常，進一步影響帕金森病患者的腦網絡信息編碼[9]。電生理學研究表明，帕金森病表現為基底節-皮層環路、皮層-皮層環路的病理振蕩，其中基底節同步活動與皮層的活動呈線性相關[10]。通過對腦電信號進行分析，可以發現基底節-丘腦-皮層回路活動改變的標記，為了解其大腦皮層異常動態變化提供新的途徑[9]。帕金森病患者的枕部Alpha與Theta頻段腦電活動信號減弱，視覺空間障礙會導致枕部右側腦電信號異常[11]。未服藥帕金森病患者的Beta頻段腦電信號波形在感覺運動電極上的振蕩比服藥后有更大的銳度不對稱和陡度不對稱性[12]。處于靜息狀態時，腦電圖顯示帕金森病患者與健康人對照組相比在Theta和Alpha頻段上功率增加，Beta和Gamma頻段上功率降低[13]。在腦電信號頻譜圖中，帕金森病患者的Beta頻段波振幅比健康人小[14]，而Alpha波段腦電信號振幅比健康人大，且其頻率減慢[15]。Zhang等[16]通過將時頻分析與深度學習相結合，應用可調Q因子小波變換深度殘余收縮網絡來對腦電數據進行分類，帕金森病預測的準確率達到99.92%。在對高分辨腦電數據源定位后，Gerard等[17]采用動態鎖相方法來評估帕金森病患者功能連接性，發現患者額區與枕顳區顯示出異常強的Theta頻帶連接性。基于腦電數據5個頻帶的相干性，Conti等[18]通過計算連接矩陣及分析網絡拓撲參數，發現帕金森病患者在Alpha、Beta頻帶連通性較低，而在Gamma頻帶連接性增強。此外，各通道腦電信號并非獨立不相關，它們之間相互影響，可能存在著某種非線性依賴關系，蘊含著可以作為檢驗帕金森病的生物學標志。為此，加權符號互信息算法（weighted symbolic mutual information，wSMI）可以用于探索腦動態功能連接性及病理條件下的腦網絡潛在變化（如阿爾茨海默病、多發性硬化癥、精神分裂癥和社交焦慮障礙）[19]。有研究使用wSMI與聚類來研究腦連接的動態過程，并對比分析了精神分裂癥患者與健康人的丘腦區域網絡效率[20]。也有人采用wSMI研究腦電信號Delta、Theta頻帶的腦皮質功能連接性[21]。
基于以上思考，本文假設帕金森病患者與正常健康人的靜息態腦電信號在各個頻段內可能存在關聯狀態差異。在計算靜息態腦電信號關聯狀態后，通過比較帕金森病患者與健康人的各通道腦電信號間的相關性來發現生物學標志。腦電信號預處理后，利用wSMI計算通道間相關性以得到動態連接矩陣，再利用k均值聚類進行關聯狀態識別，并計算出不同對象的關聯狀態發生比例，最后統計發現與帕金森病顯著相關的功能連接狀態。這里的狀態是指不同通道腦電信號之間存在的神經活動同步性。
1 方法
1.1 腦電數據
本文所采用的腦電數據為加州大學圣地亞哥分校帕金森病患者靜息狀態腦電開放數據[22]。數據集包括15名帕金森病患者［8名女性，7名男性，年齡為（63.27±8.20歲）］和15名正常健康人［9名女性，6名男性，年齡為（63.50±9.67歲）］，所有人均為右利手。在采集腦電數據時，患者繼續藥物治療方案，治療藥物主要為能影響多巴胺能的左旋多巴等。通過平衡患者用藥和停藥的時間，分別收集患者按時用藥的腦電信號（PD-on）以及患者至少停止用藥12 h（PD-off）的腦電信號。同時，收集健康人（healthy control，HC）的腦電信號作為對比分析。腦電信號采集對象信息如表1所示。腦電信號采用32通道的BioSemi ActiveTwo系統采集，并有額外的八個電極放置在眼睛的外側和下方，以監測眼睛的眨眼和運動。采集時所有人均處于舒適的靜息狀態，采樣頻率為512 Hz，靜息數據連續記錄至少3 nub。這里值得注意的是，PD-off與PD-on腦電信號采集于同一組PD患者。
表1
研究對象信息
Table1.
Demographic information of the subjects
1.2 腦電信號預處理
將腦電信號導入到EEGLAB MATLAB工具箱中進行預處理，具體步驟如下：① 坐標定位：對電極進行坐標定位，只對前32個通道采集的信號進行分析，去除后8個不反映腦電信號的電極（EXG1，EXG2，…，EXG8）。② 重定義參考：取每個電極的平均值，應用全部腦電數據的共同平均值進行重定義參考。③ 濾波：使用四階巴特沃斯數字濾波器帶通濾波，分別得到每個對象腦電信號的Delta（0.5～4 Hz）、Theta（4～8 Hz）、Alpha（8～13 Hz）、Beta（13～30 Hz）、Gamma（30～50 Hz）及高頻High-frequency（50～100 Hz）等六個頻帶區間，過濾掉其他過高頻噪聲信號和過低頻的漂移偽跡。④ 獨立成分分析：通過獨立成分分析（independent component analysis，ICA），去除眼電等生理偽影的干擾，得到純凈的腦電信號。⑤ 截段提取：考慮到每個人在儀器測量的開始和結束可能會有心理狀態的波動，未能立刻進入靜息狀態，同時各個對象所測得的腦電信號時間及數據序列長度有所不同，為了得到較為穩定的腦電信息，以腦電信號時間序列中心為基準向前后兩側各取80 s，剔除了數據最前段和最后段，只保留中間有效部分，最終得到160 s的數據，使所有對象的腦電信號時間序列長度均相等。
1.3 動態wSMI分析
本文采用wSMI來評估腦電信號隨機非線性聯合波動的程度，度量時間序列間的共享信息[21]。同時，該算法能夠發現信號增加或減少的定性“符號”模式，在不需要對交互類型進行假設的情況下能對信號的熵進行快速和穩健估計。加權符號互信息通過計算不同時間序列間同時出現的每一對符號的聯合概率來估計兩個信號之間共享符號的互信息，并通過忽略相同或相反符號信號的共同出現以解決容積傳導的問題[19]，其權值計算能夠避免來自共同來源的腦電信號間的偽相關性，使得到的結論更為嚴格。這是一種檢測非線性耦合的有效方法[23]。其計算公式如下：
|
其中、為時間序列，為中元素，為中元素。為構成符號的元素個數，為加權函數，為概率函數，為自然對數函數。
首先根據預先定義的連續采樣數據的振幅趨勢將腦電信號轉變為一系列離散符號。通過估計每對符號的聯合概率來計算腦電信號序列間的加權符號互信息值。對于一對相同的符號，權值設置為零，表示這對符號可以由唯一的共同源引出[9]。為了分析腦電信號時間序列的動態特性，不是一次性對整個時間序列進行加權符號互信息計算，而是通過施加時間窗口，每次取一段腦電信號進行處理計算。隨著時間窗口的滑動，可以得到一系列加權符號互信息計算結果，構成三維矩陣。這樣可以更好地體現出腦狀態的時間波動性。本文里，wSMI采用Python語言實現，底層調用MNE庫（MEG、EEG analysis & visualization，用來處理EEG/MEG等信號的開源程序，
1.4 k均值聚類
k均值聚類，是一種經典的高效快速聚類方法，其計算過程是把多個樣本分為k個簇，按照距離中心點最“小”的原則，使得各簇內的樣本點具有較高的相似性，而劃為不同簇的樣本相似性較低。相似度是依據一個簇中樣本的平均值來進行計算的。每次迭代后重新計算聚類中心，重新劃分類別，由此不斷迭代重復，直至聚類中心值趨于穩定[26]。
將加權符號互信息得到的關聯矩陣集合按照每個頻段分別聚類。這里聚類數目取k = 6[20]。對每個矩陣樣本計算它到聚類中心的距離并將它分到距離最小的類中，然后針對每個矩陣樣本類別重新計算聚類中心。兼顧速度和效率，這里設置重復迭代上限為200次。最后得出所有頻段的聚類中心、樣本所屬類別、概率分布和轉移概率等數據。最后得到的關聯矩陣如圖1所示。
圖1
腦電信號關聯矩陣
Figure1.
Correlation matrix between each pair of EEG channels
1.5 統計分析
不考慮其他因素的影響，把腦電數據分為健康人HC、未服藥帕金森病患者PD-off、服藥帕金森病患者PD-on三類（每組15人，共三組），作為ANOVA統計分組變量。在驗證所得的樣本數據符合正態分布后，計算HC、PD-off、PD-on各組對象在不同頻段關聯狀態占比，之后對這些樣本進行ANOVA統計分析。在發現不同組間差異有統計學意義后，采用最小顯著差別方法（least significant difference，LSD）進行校驗。檢驗水準為0.05。本論文所實現的算法源代碼與分析結果數據可從
2 結果
圖2展示了不同頻段靜息態腦電信號關聯矩陣聚類結果，每個頻帶共聚類為6個關聯狀態。本論文采用狀態所占比例作為指標來區分不同分組，探究將靜息態腦電信號關聯狀態作為帕金森病生物學標志的可能性。健康人、服藥帕金森病患者、未服藥帕金森病患者三個分組的腦電信號關聯狀態占比數據如表2所示。
圖2
對三組對象采用加權符號互信息及k均值聚類獲取靜息態腦電信號關聯狀態
每個頻段聚類得到6個關聯狀態，采用關聯狀態占比作為指標來區分不同組
Figure2. Correlation matrices extracted from EEG signals using wSMI and k-means clustering in three groupssix correlation estimates were obtained for each frequency band. The correlation status proportion was used as an index to distinguish different groups
表2
三組對象靜息態腦電信號關聯狀態占比（±s）
Table2.
Correlation status proportion of HC, PD-off, and PD-on groups (±s)
ANOVA統計分析之后，發現HC與PD-off組間Beta頻帶State1占比差異有統計學意義（P = 0.034，見表3），同時這兩組間Gamma頻帶State5占比差異也有統計學意義（P = 0.010，見表4）。這說明在帕金森病患者與健康人之間，靜息態腦電信號關聯性在Beta和Gamma頻帶顯著不同。
表3
Beta頻帶關聯狀態1占比的顯著性水平（LSD校正）
Table3.
Significance level of the proportion of state1 in Beta band (LSD test)
表4
Gamma頻帶關聯狀態State5占比的顯著性水平（LSD校正）
Table4.
Significance level of the proportion of state5 in Gamma band (LSD test)
3 討論與結論
本論文利用動態wSMI及k均值聚類實現了帕金森病患者靜息態腦電信號關聯性分析。發現在Beta段范圍內，健康人與未服藥帕金森病患者狀態State1占比差異有統計學意義（P = 0.034），而未服藥與服藥帕金森病患者之間State1占比差異沒有統計學意義（P = 0.437），健康人與服藥患者間State1占比差異也沒有統計學意義（P = 0.169）。在Gamma頻帶，健康人與未服藥帕金森病患者間State5占比差異有統計學意義（P = 0.010），而未服藥與服藥患者間State5占比差異沒有統計學意義（P = 0.560），健康人與服藥患者間State5占比差異也沒有統計學意義（P=0.053）。
考慮到k均值聚類可能存在誤差，本文進行了多次重復聚類計算，均得到同樣的結果。同時，觀察到Delta頻帶的狀態圖與其他頻段有明顯的不同（見圖2），猜想在Delta頻帶通道相關性可能具有與其他頻帶完全不同的獨特性質。此外，還發現帕金森病患者的腦電波形在后半段有更大的波動和陡峭性，而健康人的則平坦穩定，這與之前的研究結果是一致的[12]。在k均值聚類時，Delta與High-frequency頻段具有比較大的波動性，有可能是頻段過窄或過寬的原因。
另外，本文也驗證了當聚類數目分別設為4、5、7時，并沒有發現差異具有統計學意義的腦網絡關聯狀態，而只有在聚類數為6時發現兩個有顯著差異的狀態。這也驗證了Fu等[20]的工作。
同時，分析本文研究方法可能產生誤差的原因：① 獨立成分分析去除生理偽影。采用獨立成分分析去除生理偽影時，對生理偽影的判斷可能會有一定的誤差。② 高頻帶的噪聲干擾會對互信息計算有所影響。③ 采集靜息態數據時對象可能存在情緒或認知活動波動，例如患者腦神經功能的缺陷使得對象不能長時間集中精力等。
在Beta頻帶State1和Gamma頻帶State5占比上，盡管未服藥帕金森病患者與服藥帕金森病患者之間、服藥帕金森病患者與健康人之間差異沒有統計學意義，但還是可以觀察到患者服藥后其狀態占比的統計均值介于未服藥和健康人之間，說明服藥能夠縮小這兩者狀態的差距，因此服藥對于病變的程度還是有所改善的。關于帕金森病患者服藥與未服藥的差異，還可以在此基礎上進一步深入研究。同時，在本研究采用32通道腦電數據得出了上述結論，如果基于更多采集通道（如64、128等）的腦電信號數據，腦電信號關聯狀態的了解與分析可以更為透徹與精確。
wSMI能夠評估腦電信號間非隨機聯合波動的程度，度量其共享信息值。相較于其他相關性算法，加權符號互信息主要有三個優點：① 能夠尋找信號增加或減少的定性或“符號”模式，從而快速和準確地估計信號的熵，符號的轉換取決于符號的長度和時間間隔；② 提供了一種檢測非線性耦合的有效方法，并且不需要對相互作用的類型進行假設，提高了結果的真實性；③ 權重計算舍棄了腦電信號之間的偽相關性[21, 23]。這些特點使得wSMI非常合適腦電信號關聯狀態的識別。
總之，本研究采取動態加權符號互信息方法計算腦電信號間關聯性，通過聚類評估帕金森病患者腦電信號關聯狀態的改變模式，從而發現帕金森病的生物學標志。這為帕金森病患者的病情診斷與篩查提供了一種參考方法。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：丁昊負責數據分析和論文寫作，吳進輝與唐旭東負責腦電數據預處理，余江南與陳軒恒負責算法編程，吳占雄負責試驗設計與文字修改。
