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

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

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

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

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

1 實(shí)驗(yàn)材料及方法

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í)驗(yàn)數(shù)據(jù)包括30名正常被試和29名阿爾茲海默癥患者，被試基本信息如表1所示。表中數(shù)值是平均值±標(biāo)準(zhǔn)差；NC表示正常對(duì)照組；NMSE表示簡易精神狀況檢查量表；a為雙樣本T檢驗(yàn)；b為皮爾遜卡方檢驗(yàn)。

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

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

1.2 方法框架

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

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

2) 組獨(dú)立成分分析。

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

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

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

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

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

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

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

1.3 組獨(dú)立成分分析

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

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

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

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

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

(1)

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

(2)

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

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

(3)

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

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

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

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

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

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

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

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

模式增長算法如下：

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

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

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

(3)

式中：Exp(g)表示子圖模式g的期望支持度，n表示加權(quán)圖的個(gè)數(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)計(jì)指標(biāo)的方法。但這種方法受統(tǒng)計(jì)指標(biāo)的影響較大且最終的分類正確率較低，因此我們提出了一些新的統(tǒng)計(jì)指標(biāo)，并結(jié)合無權(quán)圖上的判別分?jǐn)?shù)函數(shù)(見表2)，通過分類正確率來衡量子圖模式的判別性。

表2 判別分?jǐn)?shù)函數(shù)Table 2 Discriminant fraction function

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

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

2 實(shí)驗(yàn)結(jié)果與分析

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

從圖1中可以看出，注意網(wǎng)絡(luò)中包含的成分最多，其次是感覺運(yùn)動(dòng)網(wǎng)絡(luò)和視覺網(wǎng)絡(luò)，然后是默認(rèn)網(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中顯示的是作為特征進(jìn)行分類的20個(gè)頻繁子圖模式，不同顏色的節(jié)點(diǎn)分別代表不同的靜息態(tài)腦網(wǎng)路。從中可以發(fā)現(xiàn)，感覺運(yùn)動(dòng)網(wǎng)絡(luò)、默認(rèn)網(wǎng)絡(luò)和視覺網(wǎng)絡(luò)中連接的邊比較密集，這說明這些網(wǎng)絡(luò)在患者和正常被試的分類過程中起到的作用比較大。

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

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

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

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

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

此外，為了驗(yàn)證基于獨(dú)立成分分析的高階功能連接網(wǎng)絡(luò)的優(yōu)越性，本實(shí)驗(yàn)使用相同的被試數(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特征選擇方法，進(jìn)行SVM分類，結(jié)果如表5所示。

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

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

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

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

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

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

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

3 總結(jié)

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