節點屬性和拓撲信息相結合的腦網絡聚類模型 *

肖繼海，崔曉紅，陳俊杰

(太原理工大學信息與計算機學院，山西 晉中 030600)

1 引言

近年來，圖挖掘已經成為一個熱門的研究領域并應用到社交網絡、計算生物學及計算機網絡的分析中。此外，許多新的數據也可以表示為圖，如功能性磁共振大腦影像數據fMRI(functional Magnetic Resonance Imaging)。本文利用受試者的fMRI數據，構建腦網絡，其中每個節點代表一個腦區，每條邊代表2個腦區之間的功能連通性[1]。這些腦網絡為探索人腦的內部結構和活動提供了一個新的視角,同時也可為腦疾病臨床診斷提供有價值的輔助信息。

目前的研究大都集中于腦疾病的分類，例如，Huang等人[2]提出了一種基于樹引導的組稀疏學習多頻帶腦網絡分類模型，首先在多頻帶中獲取低頻振幅比率數據，使用基于樹的稀疏學習方法進行特征選擇，最后結合多頻段特征進行分類。Liu等人[3]從T1加權MRI(Magnetic Resonance Imaging)圖像和DTI(Diffusion Tensor Imaging)圖像中提取8個特征集，并采用兩步特征選擇方法得到最具鑒別性的特征，最后采用基于多特征集的多核SVM學習方法對精神分裂癥SZ(SchiZophrenia)患者和正常人NC(Normal Control)進行分類。楊勇哲等人[4]提取結構和靜息態功能磁共振圖像的90個感興趣區域的灰質體積、局部一致性、低頻振蕩振幅和度中心性作為特征，使用基于遞歸特征消除的支持向量機對腦疾病患者和正常人進行分類。

以上研究旨在選擇特征來訓練分類器，以實現對數據的自動分類。當訓練數據樣本量足夠大并且選擇了最優特征時，分類性能通常很好。但是，從大量的數據中提取和選擇最優的特征很費時。此外，在許多實際情況下，只能獲取到較少的標記數據樣本。因此，在有限的標記數據樣本上找到有顯著差異的特征是不可靠的。

針對以上問題，本文提出了一種節點屬性和拓撲信息相結合的腦網絡聚類模型。與分類相反，聚類的不同之處在于不需要數據的類型標簽[5]。腦網絡聚類的目的是根據腦網絡之間的相似程度自動聚類，類內腦網絡之間具有較高的相似性[6]。

該方法主要的挑戰是如何準確計算腦網絡之間的相似度。目前，文獻[7]調查了48名SZ患者和24名健康對照者的全局和節點網絡特性，發現精神分裂癥患者的大腦網絡特性與健康人的相比有顯著差異，比如SZ患者的平均介數中心性較高。另外，文獻[8]發現SZ患者功能性磁共振網絡的拓撲和空間結構被顯著破壞。以上文獻表明，SZ患者的腦網絡在節點屬性和拓撲結構上都存在顯著差異。因此，為了更全面地捕獲腦網絡的差異，本文從節點屬性和拓撲結構兩方面測量腦網絡的相似度。首先，選擇目前研究公認的SZ患者存在異常功能連接的默認模式網絡DMN(Default Mode Network)[9,10]作為感興趣區域構建腦網絡，然后分別使用余弦相似度和子網絡核來度量腦網絡的屬性相似度和結構相似度，接著將結構相似度和屬性相似度集成為一個相似度矩陣，最后利用譜聚類實現腦網絡聚類。

2 腦網絡聚類模型

腦網絡聚類模型如圖1所示，包含fMRI數據的預處理、標記DMN、構建相應腦網絡、計算腦網絡的相似度以及聚類。也就是，逐一對每個腦網絡的fMRI數據先實施預處理，找到DMN對應的腦區并構造腦網絡，隨后計算其屬性相似度和結構相似度，最后整合2種不同類型的相似度，并使用譜聚類方法對腦網絡進行聚類。

Figure 1 Brain network clustering model圖1 腦網絡聚類模型

2.1 預處理數據

使用DPARSF(Data Processing Assistant for Resting-State fMRI)工具包對fMRI圖像數據進行預處理，包括刪除前10個時間點數據、時間片的校正和頭動的校正。本文使用DARTEL(Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra)配準方法將數據配準到標準空間，期間將數據重新采樣為3×3×3 mm3體素，再進行平滑處理以去除協變量的影響，最后，進行濾波去噪，頻率范圍為[0.01,0.08] Hz，以減少高頻和極低頻生理噪聲的影響，這樣能較好地反映神經元自發活動。

2.2 標記DMN并構建功能腦網絡

研究已證實，在靜息態功能腦網絡中默認模式網絡(DMN)能保持相對穩定的狀態，通過它可以研究功能腦網絡連接的異常[10 - 12]。本文使用自動化解剖標記AAL(Automated Anatomical Labeling)圖譜[13]，從中找到DMN所包含的32個腦區(左右腦)，腦區名稱如表1所示。

Table 1 Brain regions included in the DMN

腦網絡的節點由表1中的32個腦區定義，節點的時間序列取腦區內所有體素時間序列的平均值，節點之間的連接權值由節點時間序列間的皮爾遜相關系數定義，加權全連接網絡由此形成。為了避免傳統二值網絡構建過程中產生的偏差，使用Kruskal算法[14]將加權全連接網絡轉換成無偏差的腦網絡，最終形成最小生成樹無偏DMN腦網絡。

2.3 測量腦網絡相似度

腦網絡不僅具有全局拓撲結構特征而且還具有局部屬性特征。為了綜合評價腦網絡的相似度，本文將上述2方面特征結合起來測量腦網絡相似度，流程圖如圖2所示。圖2中，G和H為腦網絡，S(G,H)為腦網絡G和H的相似度，Satt(G,H)為G和H的屬性相似度，Sstr(G,H)為G和H的結構相似度，δ為權重。

Figure 2 Flow chart of similarity measurement in brain network圖2 腦網絡相似度測量流程圖

(1)節點屬性相似度的測量。

介數屬性是復雜網絡的一個重要特征。節點介數越大，表示該節點是中樞節點，對整個網絡的信息傳輸起到很大的作用。研究表明[15]，SZ患者的大腦中樞受到破壞，在臨床應用中，介數常用于健康人和SZ患者腦網絡的比較。因此，本文使用節點的介數來表征腦網絡的局部屬性，腦網絡間屬性相似度由腦網絡中所有節點介數之間的相似度來評價。通過節點i的所有最短路徑的數量定義為節點i的介數bci，計算如式(1)所示：

(1)

屬性相似度Satt(G,H)使用余弦相似度計算，如式(2)所示：

(2)

其中,bcm(H)表示腦網絡H中第m個節點的介數，bcm(G)表示腦網絡G中第m個節點的介數，n表示腦網絡中節點總數。

(2)拓撲結構相似度測量。

一組對象之間的相似性可以用核函數來度量。當用核函數處理圖數據時，該核稱為圖核，圖核是把圖數據從原始的圖空間映射到特征空間，使得2幅圖像的相似性就是它們在特征空間中的點積。因此，2幅圖之間的相似性可通過圖核來度量[16]。

大腦網絡與圖的主要區別是大腦網絡中每個節點的唯一性，也就是說，大腦網絡中的每個節點都表示特定的大腦區域，是大腦網絡的一個特定特征。因此，本文采用子網絡核(Sub-network Kernel)[17]測量腦網絡間的拓撲結構相似度。與傳統意義上的圖核相比，子網絡核的一個重要特性就是能捕獲網絡節點上的多層次拓撲性質，而且還兼顧了網絡中各個節點的唯一性。子網絡核的構建過程簡要描述如下：

①在每個網絡節點上構建一組子網絡，以便反映腦網絡在多層次上的連通性。

G=(V,E)和與H=(V,E′)表示一對腦網絡，其中V代表腦網絡的節點集，E和E′分別表示G與H的邊集。各個腦網絡擁有相同的大腦區域，即它們有相同的節點。為了表達腦網絡的多級拓撲特性，需要在網絡G與網絡H的每個節點vi上分別定義2組子網絡：

(3)

(4)

由式(3)和式(4)可知，對于由n個節點組成的腦網絡，可以得到n組子網絡：

(5)

(6)

②網絡G與H中同一節點的相似度就是2個網絡中同一節點的所有子網絡的相似度，如式(7)所示：

(7)

(8)

(9)

(10)

③因此，腦網絡G與H的圖核就是腦網絡中所有節點對應子網絡對的相似度，n表示腦網絡中的節點數目，G和H的圖核如式(11)所示：

(11)

④最后得到2個腦網絡G與H的圖核即為2個腦網絡G與H之間的拓撲結構相似度，如式(12)所示：

Sstr(G,H)=k(G,H)

(12)

(3)腦網絡的相似度。

為了全面評價腦網絡之間的相似度，本文將腦網絡之間屬性方面的相似度與結構方面的相似度進行優化組合。

首先，因為屬性相似度與結構相似度的類型不同，故在組合二者之前必須歸一化，如式(13)所示：

(13)

其中，Snorm(G,H)是腦網絡G與H歸一化之后的相似度，S(G,H)是腦網絡G與H的相似度。

然后，通過權重δ來調整2種相似度的貢獻程度，使組合最優。

腦網絡G與H的相似度定義如公式(14)所示：

SG,H=δSatt_norm(G,H)+

(1-δ)Sstr_norm(G,H)

(14)

其中，Satt_norm(G,H)是功能腦網絡G與H歸一化后的屬性相似度，Sstr_norm(G,H)是功能腦網絡G與H歸一化后的結構相似度，δ是權重參數。

第三步，腦網絡之間的相似矩陣S定義如下：

(15)

其中,SG,H是腦網絡G與H之間的相似度，N是腦網絡數量。

2.4 腦網絡聚類模型描述

前文已經構造出腦網絡的相似矩陣S，這樣可以方便地將腦網絡的聚類問題轉化成譜聚類[19]問題，將相似度高的腦網絡自動歸為一類。聚類模型的基本過程描述如下：

輸入：N個腦網絡。

輸出：向量Y，若向量Y的第i個元素為cj，則表示第i個被試對應的腦網絡歸屬到類cj中。

步驟1將腦網絡的相似矩陣S置為N×N零元矩陣，N表示腦網絡數目；

步驟2對N個腦網絡依次計算每個節點的介數；

步驟3由式(2)得到N個腦網絡之間的屬性相似度；

步驟4對N個腦網絡執行子網絡核算法，得到N個腦網絡之間的結構相似度；

步驟5將屬性相似度和結構相似度歸一化，由式(14)與式(15)得到相似矩陣S；

步驟6定義對角陣D，其元素D(i,i)為相似矩陣S的第i行元素之和；

步驟7定義矩陣L=D-1/2SD-1/2;

步驟8求矩陣L的m個最大特征向量x1,x2,…,xm，按列存放m個特征向量產生矩陣X=[x1,x2,…,xm]∈Rn×m；

步驟9將矩陣X的每一行標準化，得到矩陣Y；

步驟10將矩陣Y的每一行作為Rm中的一個點，執行K-means算法將這些點聚成m個類。

3 實驗數據與結果

3.1 實驗數據

本文實驗的數據來自于openfMRI數據庫(https://www.openfmri.org/)，經篩選，滿足數據收集條件的有99名，其中包含50名SZ患者與49名正常對照組。掃描參數：TR(Repetition Time)重復時間為2 000 ms；TE(Echo Time)回波時間為30 ms；翻轉角度為90°；切片厚度為4 mm；切片數量為34；時間點數為152個。被試的基本信息統計情況如表2所示。

Table 2 Statistical table of basic information of subjects

3.2 實驗參數設置

實驗中，因為僅有2個可能的標簽(SZ，NC)，故將類別數m設置為2。在子網絡核算法中，參數t與d分別設置為3與3。權重δ通過網格搜索法來優化，δ初始值是0.1，步長是0.1，最大值是0.9。

3.3 無偏腦網絡

為了構建合理的腦網絡，本文腦網絡的節點由DMN包含的32個腦區定義，節點之間的連接權值由節點時間序列間的皮爾遜相關系數定義。使用Kruskal算法得到SZ與NC所有被試的無偏腦網絡。SZ(No.01)與NC(No.01)的無偏腦網絡如圖3所示。從圖3中看出，SZ患者和正常人的大腦在額葉、顳葉和頂下小葉之間的功能連接存在明顯的差異。

Figure 3 Unbiased brain network of SZ and NC圖3 SZ與NC的無偏腦網絡

3.4 相似矩陣

譜聚類的關鍵在于相似性度量。腦網絡構建后，腦網絡之間的相似度通過整合腦網絡在屬性與結構2方面的相似度來綜合評價。本文首先使用余弦相似度來計算腦網絡間的介數屬性相似度，腦網絡的屬性相似度矩陣如圖4所示；其次，使用子網絡核來度量腦網絡間的結構相似性，其相似度矩陣如圖5所示；最后，選擇最優權重參數δ，將屬性相似度與結構相似度進行優化組合，得到最終相似度矩陣，如圖6所示。

Figure 4 Attribute similarity matrix圖4 屬性相似度矩陣

3.5 聚類性能

將屬性相似度矩陣與結構相似度矩陣按一定比例優化組合，形成最終相似度矩陣。然后通過譜聚類[20]完成上述腦網絡的聚類。通過以下指標來評估聚類方法的性能： Rand指數、召回率、精確率和F-measure[21]。為準確評價上述指標，實驗中要求重復運行程序50次，將其平均值作為最終的性能評價結果。

Figure 5 Structural similarity matrix圖5 結構相似度矩陣

Figure 6 Final similarity matrix圖6 最終相似度矩陣

(1)Rand指數表示正確聚類的百分比，Rand指數越大，意味著聚類結果與真實情況越吻合。

Rand指數由式(16)計算：

(16)

其中,TP是真正類的個數，TN是真負類的個數，FP是假正類的個數，FN是假負類的個數。

(2)召回率(Recall)表示在被識別為正類別的樣本中，真實情況是正類別的比例，即查全率，由式(17)計算：

(17)

(3)精確率(Precision)表示在所有正類別樣本中，被正確識別為正類別的比例，即查準率,由式(18)計算：

(18)

(4)F-measure綜合了精確率和召回率的結果，是精確率與召回率加權調和平均值。F-measure由式(19)定義：

(19)

其中,P表示精確率，R表示召回率。在本文中，使用F1(β=1)評價聚類的性能。當F1的值較高時，則說明聚類模型的效果比較理想。

為評估本文所提方法SA-cluster(Structure Attribute-cluster)的聚類性能，在相同數據集上設計實驗,將本文的方法與2種方法進行比較：

(1)基于介數屬性的譜聚類算法(A-cluster)：構造基于介數屬性的相似矩陣，將相似矩陣輸入譜聚類模型，對腦網絡進行聚類。

(2)基于結構的譜聚類方法(S-cluster)：構造基于拓撲結構的相似矩陣，將相似矩陣輸入譜聚類模型，對腦網絡進行聚類。

相同數據集上3種方法的聚類性能如表3所示。結果表明，本文方法(SA-cluster)的聚類性能最好，其Rand指數是0.91，精確率是0.86，召回率是0.98，F1是0.92。

Table 3 Clustering performance of three methods

4 實驗結果討論

4.1 聚類性能的比較

如表3所示，本文方法在Rand指數、精確率、召回率和F1方面表現最出色。3種方法中，A-cluster方法使用腦網絡的介數屬性計算每對腦網絡之間的相似性，并用于聚類。A-cluster方法的Rand指數與召回率比較低，這說明該方法不能較準確地識別出SZ患者，存在漏診。結果表明，僅僅從局部屬性角度不能準確反映腦網絡之間的相似性。

S-cluster通過在腦網絡上執行子網絡核算法來得到相似矩陣。從表3可以看出，與A-cluster方法相比，該方法的聚類結果有所改善，查準率和查全率提高了，漏診率下降，能較準確地識別出患者。結果說明，當網絡的局部屬性發生變化時，網絡的全局連接結構也受到影響。因此，全局結構的相似性描述方法對聚類有很大貢獻。

與前2種方法相比，本文方法(SA-cluster)的Rand指數是0.91，精確率是0.86，召回率是0.98，F1是0.92，準確率和查全率最高，漏診率最低，更能準確地識別出患者，聚類的效果比較理想。這主要是因為節點屬性和拓撲結構從不同的角度描述腦網絡的相似度，將兩者有效地結合可以更好地表達腦網絡之間的相似度。這一結論正好與文獻[7,8]的研究結果一致，進一步證實了SZ患者的大腦網絡在節點屬性和拓撲結構上與正常人的腦網絡存在差異。

基于上述分析，本文方法巧妙地整合了屬性相似度與結構相似度，其中以介數屬性相似度來刻畫腦網絡的局部特征，以結構相似度來刻畫腦網絡的全局拓撲結構特征。特別地，本文使用子網絡核算法來度量腦網絡的拓撲結構相似性，該算法不僅考慮到腦網絡中每個腦區節點的唯一性，而且還能夠捕獲到腦區節點的多級拓撲特征，這一點對于腦網絡的相似性度量非常重要。可以看出，屬性特征與內部節點拓撲結構對圖聚類具有重要意義。因此，基于屬性與結構相結合的相似性度量方法可以更準確地刻畫腦網絡的相似性，從而有效提高聚類性能。

4.2 參數δ對聚類性能的影響

權重參數δ用于調整屬性相似度與結構相似度的貢獻程度。為分析權重參數δ對本文方法聚類性能的影響，本文為參數δ設置了不同的值，δ初始值是0.1，遞增步長是0.1，終值是0.9。參數δ取不同值的聚類結果如圖7所示。

Figure 7 Influence of different δ values on clustering performance圖7 參數δ不同取值對聚類性能的影響

圖7表明，無論δ取何值，Rand指數、精確率、召回率和F1值都在0.8以上，說明從屬性和結構2方面能較準確地衡量腦網絡相似度，腦網絡的拓撲結構和屬性特征對腦網絡的聚類起著重要作用。當δ取值為0.2時，聚類性能最好，其Rand指數為0.91,精確率為0.86,召回率為0.98,F1值為0.92，這說明在相似度的計算中，結構相似度貢獻較大。

5 結束語

本文提出一種節點屬性和拓撲信息相結合的腦網絡聚類模型。首先，使用余弦相似度來計算腦網絡之間的屬性相似度，然后，使用子網絡核來度量腦網絡之間的結構相似度，最后，將屬性相似度與結構相似度相結合，得到相似度矩陣，并用譜聚類算法聚類。實驗結果表明，本文提出的聚類模型在Rand指數、召回率、精確率和F1值方面都表現較優。