楊艷麗，李 瑤，谷金曄，李欣蕓，陳俊杰

(太原理工大學(xué) a.信息與計算機學(xué)院，b.軟件學(xué)院，太原 030024)

利用功能性磁共振成像技術(shù)，從復(fù)雜網(wǎng)絡(luò)分析角度研究腦疾病的發(fā)病機理，已成為腦科學(xué)領(lǐng)域的研究熱點之一[1]。該方法主要包括：節(jié)點的定義，連接的定義以及網(wǎng)絡(luò)的分析。節(jié)點定義最常用的方法是利用腦圖譜劃分節(jié)點；然而這種節(jié)點的定義方法過于依賴先驗的腦模板，不同模板建立起來的網(wǎng)絡(luò)可能會存在顯著差異，這將嚴重影響功能連接網(wǎng)絡(luò)的分析結(jié)果[2]。ALLEN在研究中使用組獨立成分分析(group independent component analysis，group ICA)提取到的獨立成分作為功能的網(wǎng)絡(luò)節(jié)點[3]。

傳統(tǒng)的靜息態(tài)功能連接假定功能連接在時間上是靜止的，功能連接網(wǎng)絡(luò)在掃描時間內(nèi)穩(wěn)定不變。這種方法的主要問題是忽略了掃描時間內(nèi)可能發(fā)生的神經(jīng)活動或相互作用[4]。然而，最近的研究已經(jīng)表明大腦的功能連接承載著豐富的時間信息，無論是任務(wù)態(tài)還是靜息態(tài)，功能連接都會隨著神經(jīng)活動的變化而變化[5]。CHEN使用滑動窗口的方法來劃分整個RS-fMRI時間序列，從而構(gòu)建出了高階功能連接網(wǎng)絡(luò)[6]。這種動態(tài)功能連接方法，在傳統(tǒng)的功能連接網(wǎng)絡(luò)上添加了時變因素，充分考慮了時間特性。

目前腦網(wǎng)絡(luò)的分析方法有很多，其中最常用的是使用網(wǎng)絡(luò)拓撲指標(biāo)來描述腦網(wǎng)絡(luò)特征[7]。然而這種方法受網(wǎng)絡(luò)規(guī)模和稀疏度的影響較大，且容易忽視網(wǎng)絡(luò)中一些有用的拓撲結(jié)構(gòu)。DU先對腦網(wǎng)絡(luò)進行閾值化處理，然后在無權(quán)圖上使用頻繁子圖挖掘和判別性特征選擇方法[8]。但閾值化處理時，無論是單閾值還是多閾值，都會造成網(wǎng)絡(luò)信息的丟失。KONG et al通過構(gòu)建腦網(wǎng)絡(luò)的加權(quán)圖(邊以一定概率存在)，進行加權(quán)圖的頻繁子圖挖掘，并將無權(quán)圖上的判別性特征選擇方法推廣到了加權(quán)圖上[9]。

1 實驗材料及方法

1.1 數(shù)據(jù)采集及預(yù)處理

本文的所有數(shù)據(jù)都來自ADNI(Alzheimer's disease neuroimaging initiative)數(shù)據(jù)庫。ADNI數(shù)據(jù)庫自建立以來，為阿爾茲海默癥(Alzheimer's disease，AD)相關(guān)疾病及早期診斷方法的研究提供了豐富的數(shù)據(jù)支持[11]，得到了廣泛的使用。本實驗數(shù)據(jù)包括30名正常被試和29名阿爾茲海默癥患者，被試基本信息如表1所示。表中數(shù)值是平均值±標(biāo)準差；NC表示正常對照組；NMSE表示簡易精神狀況檢查量表；a為雙樣本T檢驗；b為皮爾遜卡方檢驗。

表1 被試基本信息Table 1 Basic information for participants

原始數(shù)據(jù)預(yù)處理使用SPM8(http://www.fil.ion.ucl.ac.uk/spm)進行。主要步驟包括1) 時間片校正：校正圖像在獲取時間上的差異；2) 頭動校正：刪除頭動大于3 mm或轉(zhuǎn)動大于3度的被試數(shù)據(jù)，丟棄的被試數(shù)據(jù)不包含在最終的被試中；3) 空間標(biāo)準化：對圖像進行12維度的優(yōu)化仿射變換，將其標(biāo)準化到3 mm體素的MNI(montreal neurological institute)標(biāo)準空間中；4) 平滑：消除不同被試腦結(jié)構(gòu)之間的細微差別、提高信噪比，在10 mm半高全寬高斯空間上進行平滑。

1.2 方法框架

基于組獨立成分分析的加權(quán)高階功能連接網(wǎng)絡(luò)的靜息態(tài)功能性磁共振成像數(shù)據(jù)分類方法包括以下4個步驟。

1) 數(shù)據(jù)采集及預(yù)處理。

2) 組獨立成分分析。

組獨立成分分析的主要步驟包括：數(shù)據(jù)降維、空間獨立成分估計、數(shù)據(jù)反重構(gòu)、篩選腦網(wǎng)絡(luò)成分[12]。

3) 高階功能連接網(wǎng)絡(luò)的構(gòu)建,主要包括兩個過程:

a.選擇長度固定的滑動窗口，對獨立成分的時間序列進行時間窗劃分，并計算每個時間窗下的皮爾遜相關(guān)系數(shù)矩陣，得到低階功能連接網(wǎng)絡(luò)；

b.堆棧所有的低階功能連接網(wǎng)絡(luò)，再在其上計算皮爾遜相關(guān)系數(shù)矩陣，得到高階功能連接網(wǎng)絡(luò)；

4) 加權(quán)圖的頻繁子圖挖掘及特征選擇,主要包括兩個過程：

a.對構(gòu)建的全連接加權(quán)的腦網(wǎng)絡(luò)進行頻繁子圖挖掘，得到加權(quán)的頻繁子圖模式；

b.計算基于動態(tài)規(guī)劃的判別分數(shù)函數(shù)的相關(guān)統(tǒng)計指標(biāo)。

1.3 組獨立成分分析

組獨立組分分析使用GIFT(http://mialab.mrn.org/software/gift)工具箱進行。本文采用信息極大化(Infomax)算法進行獨立成分的估計和分離，該算法的核心思想是通過最大化輸入和輸出的互信息等價地達到輸出各分量間互信息的最小化[13]。首先使用GIFT軟件估計獨立成分的個數(shù)，并結(jié)合先驗的研究信息[14]，將獨立成分的個數(shù)設(shè)為54；接著標(biāo)準化獨立成分以消除被試個體差異的影響；隨后為了確保獨立成分的穩(wěn)定性和可靠性，在ICASSO軟件上重復(fù)了20次Infomax算法；最后使用數(shù)據(jù)反重構(gòu)方法獲取被試獨立成分的空間分布和時間序列。

使用組獨立成分分析提取的獨立成分既包括本文感興趣的腦網(wǎng)絡(luò)成分，也包括其他無關(guān)成分或含噪聲較多的成分。因此需要使用先驗?zāi)０迤ヅ浞椒▽@些獨立成分進行篩選，并結(jié)合人工肉眼檢查方法對這些成分進行進一步的確認。腦網(wǎng)絡(luò)成分的篩選標(biāo)準包括：與先驗?zāi)０迤ヅ涞亩嘀鼗貧w系數(shù)較大；主要激活區(qū)域分布在灰質(zhì)區(qū)域。最終去除了32個無關(guān)成分或含噪聲較多的成分后保留了22個腦網(wǎng)絡(luò)成分，并確認了這些獨立成分分別屬于聽覺網(wǎng)絡(luò)、感覺運動網(wǎng)絡(luò)、視覺網(wǎng)絡(luò)、默認網(wǎng)絡(luò)、注意網(wǎng)絡(luò)和額葉網(wǎng)絡(luò)。

1.4 高階功能連接網(wǎng)絡(luò)

在構(gòu)建高階功能連接網(wǎng)絡(luò)之前，還需要對被試獨立成分的時間序列進行后處理，包括去線性漂移和低通濾波，以降低低頻漂移和高頻生理噪聲的影響。接著在每個獨立成分的時間序列上，選擇長度固定的滑動窗口，按照一定步長進行時間窗滑動。假設(shè)時間序列的長度為K，滑動窗口的長度為w，滑動窗口每次移動的步長為s，則可以將整個時間序列(維數(shù)為K×N，N為獨立成分的個數(shù))劃分為l個相互重疊的時間窗口(維數(shù)為w×H).其中l(wèi)的計算公式如下：

l=[(K-w)/s+1] .

(1)

接著在每個時間窗下，計算獨立成分兩兩之間的皮爾遜(Person)相關(guān)系數(shù)，就可以得到l個低階功能連接網(wǎng)絡(luò)(維數(shù)為N×N)，它反映了所有獨立成分之間連接強度隨時間的變化情況。其中皮爾遜相關(guān)系數(shù)的計算公式如下：

(2)

式中：ri,j表示獨立成分i與獨立成分j時間序列的皮爾遜相關(guān)系數(shù)，cov(i,j)表示兩個獨立成分時間序列的協(xié)方差，σi和σj分別表示兩個獨立成分i與j時間序列的標(biāo)準差。

最后在新生成的相關(guān)時間序列上，再次計算獨立成分兩兩之間的皮爾遜相關(guān)系數(shù)，就得到了高階功能連接網(wǎng)絡(luò)，維度為

(3)

所謂高階體現(xiàn)在矩陣中的任一元素都能反映了一條功能連接對另一條功能連接的影響，體現(xiàn)了多達4個獨立成分之間的交互作用。

1.5 加權(quán)圖的頻繁子圖挖掘

在高階功能連接網(wǎng)絡(luò)中，功能連接之間的皮爾遜相關(guān)系數(shù)可以視為功能連接之間交互的概率。因此依據(jù)概率的定義，需要將負值的皮爾遜相關(guān)系數(shù)按零值處理。本文在加權(quán)的高階功能連接網(wǎng)絡(luò)上進行頻繁子圖挖掘，頻繁子圖挖掘的目的就是在整個網(wǎng)絡(luò)中找出出現(xiàn)概率較高的連接模式(即子圖)[15]。

加權(quán)圖的頻繁子圖挖掘算法如下：

1) 遍歷加權(quán)的高階功能腦網(wǎng)絡(luò)，得到所有的1-子圖模式，按照公式(3)計算期望支持度。

2) 如果1-子圖模式的期望支持度大于給定支持度，則把它加入到1-子圖模式集合中。

3) 對每個1-子圖模式，調(diào)用模式增長算法，得到新的子圖模式集合，繼續(xù)使用公式(3)計算期望支持度。同樣地，如果新子圖模式的期望支持度大于給定支持度，就把它加入到k-子圖模式集合中(k為新子圖模式的邊數(shù))。

4) 重復(fù)上步，直到?jīng)]有滿足條件的子圖生成為止。

模式增長算法如下：

1) 將(k-1)-子圖模式所在子搜索空間標(biāo)記為i.

2) 在1-子圖模式集合中找出標(biāo)號>i的子圖，將其添加到(k-1)-子圖上，得到k-子圖。如果k-子圖模式連通，則將其加入到k-子圖模式集合中。

給定一個圖集G，一個子圖g的期望支持度定義為：

(3)

式中：Exp(g)表示子圖模式g的期望支持度，n表示加權(quán)圖的個數(shù)；Pr(g?Gj)表示子圖模式g在加權(quán)圖Gj上的概率；k表示子圖模式g的邊數(shù)；pj(ek)表示子圖模式g的邊ek在加權(quán)圖Gj上的概率。

上述加權(quán)圖的頻繁子圖挖掘中挖掘出的子圖有很多，需要從中選出具有判別性的子圖作為特征，本文采用的方法是KONG et al[9]提出的基于統(tǒng)計指標(biāo)的方法。但這種方法受統(tǒng)計指標(biāo)的影響較大且最終的分類正確率較低，因此我們提出了一些新的統(tǒng)計指標(biāo)，并結(jié)合無權(quán)圖上的判別分數(shù)函數(shù)(見表2)，通過分類正確率來衡量子圖模式的判別性。

表2 判別分數(shù)函數(shù)Table 2 Discriminant fraction function

表3 統(tǒng)計指標(biāo)公式Table 3 Statistical indicator formula

表3中前3種：期望(mean)、方差(variance)和φ-probability，是KONG et al[9]在實驗中使用的方法。受該方法的啟發(fā)，本文提出了兩種新的統(tǒng)計指標(biāo)——修正的樣本均值(trimmean)和任意階中心距(moment)，來驗證這種方法的合理性。另外，表3中只給出了正樣本中的計算方法，對負樣本類似。

2 實驗結(jié)果與分析

2.1 獨立成分分析結(jié)果

從圖1中可以看出，注意網(wǎng)絡(luò)中包含的成分最多，其次是感覺運動網(wǎng)絡(luò)和視覺網(wǎng)絡(luò)，然后是默認網(wǎng)絡(luò)，額葉網(wǎng)絡(luò)和聽覺網(wǎng)絡(luò)中包含的成分最少。這說明注意網(wǎng)絡(luò)的分布較廣，而聽覺網(wǎng)絡(luò)的分布較窄，這與以往基于腦圖譜的全腦功能連接的研究結(jié)果基本一致。

圖1 靜息態(tài)腦網(wǎng)絡(luò)成分的峰值區(qū)域Fig.1 Peak area of static brain network

2.2 頻繁子圖挖掘結(jié)果

圖2中顯示的是作為特征進行分類的20個頻繁子圖模式，不同顏色的節(jié)點分別代表不同的靜息態(tài)腦網(wǎng)路。從中可以發(fā)現(xiàn)，感覺運動網(wǎng)絡(luò)、默認網(wǎng)絡(luò)和視覺網(wǎng)絡(luò)中連接的邊比較密集，這說明這些網(wǎng)絡(luò)在患者和正常被試的分類過程中起到的作用比較大。

圖2 頻繁子圖模式Fig.2 Frequent subgraph mode

2.3 分類結(jié)果與分析

由于被試人數(shù)有限，本文使用基于RBF核函數(shù)的支持向量機(support vector machine,SVM)進行分類，并使用留一交叉驗證(leave-one-out cross validation,LOOCV)方法來評估分類器的的泛化性能(具體實驗使用基于MATLAB的LIBSVM工具包進行)。分類器的性能可以量化為正確率(accuracy)，靈敏度(sensitivity)和特異度(specificity).其中，靈敏度可以衡量該方法正確識別患者的能力，特異度衡量的是正確識別正常被試的能力。進一步使用ROC曲線下的面積(the area under receiver operating characteristic curve，AUC)來評價不同方法的性能表現(xiàn)，具體結(jié)果如表4所示。

從表4中可以看出，本文提出的兩種統(tǒng)計指標(biāo)(修正的樣本均值和任意階中心距)，都可以很好地找到與分類相關(guān)的特征，獲得較高的分類正確率；但與φ-probability方法相比，在分類正確率、靈敏度、特異度以及AUC值上均有些差距。雖然在Confidence判別函數(shù)上這兩種方法的分類效果比φ-prob方法好，但分類效果最好的還是φ-prob-HSIC方法，其分類正確率達到95.24%。從判別分數(shù)函數(shù)角度來看，雖然最好的分類效果出現(xiàn)在HSIC判別函數(shù)上，但整體來看，還是Gtest判別函數(shù)的分類效果較為穩(wěn)定，分類正確率都在90%以上。因此，使用φ-prob統(tǒng)計指標(biāo)和Gtest判別函數(shù)的分類效果較好。

表4 不同特征選擇的分類結(jié)果比較Table 4 Comparison of classification results of different feature selection

此外，為了驗證基于獨立成分分析的高階功能連接網(wǎng)絡(luò)的優(yōu)越性，本實驗使用相同的被試數(shù)據(jù)分別構(gòu)建了基于AAL的低階功能連接網(wǎng)絡(luò)、基于AAL的高階功能連接網(wǎng)絡(luò)、基于group ICA的低階功能連接網(wǎng)絡(luò)和基于group ICA的高階功能連接網(wǎng)絡(luò)，并使用加權(quán)圖的頻繁子圖挖掘和φ-prob-Gtest特征選擇方法，進行SVM分類，結(jié)果如表5所示。

表5 不同方法的分類結(jié)果比較Table 5 Comparison of classification results of different methods

從表5中可以看出，使用本文方法時，在group ICA提取的獨立成分的時間序列上構(gòu)建高階功能連接網(wǎng)絡(luò)，并以加權(quán)圖的頻繁子圖挖掘出的子圖作為特征進行分類識別，分類正確率達到95.24%，遠高于其他方法。而且高階功能連接腦網(wǎng)絡(luò)上的分類效果要遠好于傳統(tǒng)的低階功能連接網(wǎng)絡(luò)，這表明高階功能連接腦網(wǎng)絡(luò)上包含更多的對分類識別有用的信息。另外本文方法在靈敏度和特異度等指標(biāo)上也較傳統(tǒng)方法有一定的提升，在一定程度上能夠為臨床診斷提供客觀參考。

2.4 不同參數(shù)的影響

本文所使用的高階功能連接網(wǎng)絡(luò)會受一些參數(shù)的影響，如滑動窗口的長度w和每次移動的步長s等。由公式(1)可知，選擇不同參數(shù)時，得到的時間窗數(shù)l不同，進而構(gòu)建出來的高階功能連接網(wǎng)絡(luò)也不同。因此，本實驗先在滑動窗口每次移動步長都為1 s的情況下，分別選擇滑動窗口的長度為23，43，63，83，103 s，探究滑動窗口的長度對分類效果的影響，結(jié)果如圖3(a)所示。當(dāng)滑動窗口的長度為63 s時，分類效果最好，且分類曲線呈開口向下的拋物線形狀。這可能是因為：當(dāng)選擇的滑動窗口的長度過小時，劃分出來的時間窗數(shù)很多，其中包含的特征數(shù)更多，從中選擇判別性特征的難度將增大；而選擇的滑動窗口的長度過大時，劃分出來的時間窗數(shù)很少，其中所包含的時變特征更少，這會大大影響了分類效果。此外，本實驗還在滑動窗口的長度都為63 s時，分別選擇滑動窗口每次移動的步長為1 s，

圖3 不同參數(shù)的影響Fig.3 Impact of different parameters

2 s，3 s，4 s，5 s，探究滑動窗口每次移動的步長對分類結(jié)果的影響，結(jié)果如圖3(b)所示。當(dāng)滑動窗口每

次移動的步長為1 s時，分類效果最好，且每次移動的步長越長，分類的效果越差。這可能是因為：滑動窗口每次移動的步長越長，得到的時間窗數(shù)越少，其中所包含的時變特性就更少，分類正確率也會隨之降低。因此本實驗中選擇的滑動窗口長度為63 s和每次移動步長為1 s.

3 總結(jié)

本文在對靜息態(tài)功能性磁共振成像數(shù)據(jù)分類時，先使用組獨立成分分析方法提取靜息態(tài)腦網(wǎng)絡(luò)成分，接著構(gòu)建高階功能連接網(wǎng)絡(luò)，然后使用加權(quán)圖的頻繁子圖挖掘和判別性特征選擇方法，最后進行SVM分類。結(jié)果表明，該方法具有較好的分類效果，可以為臨床醫(yī)學(xué)診斷提供參考依據(jù)。然而，該方法構(gòu)建高階功能連接網(wǎng)絡(luò)的過程較為復(fù)雜，大大增加了運算量；而且該方法在加權(quán)圖上做頻繁子圖挖掘時存在大量的冗余操作，效率較低，只適用于小樣本數(shù)據(jù)集。因此，如何降低高階功能連接網(wǎng)絡(luò)的復(fù)雜性和加權(quán)圖的頻繁子圖挖掘的冗余度，是今后研究的重點。