程 超，黨偉超，白尚旺，潘理虎，2，劉春霞

(1.太原科技大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，太原 030024；2.中國科學(xué)院 地理科學(xué)與技術(shù)學(xué)院，北京 100101)

靜息態(tài)功能磁共振成像檢測大腦自發(fā)低頻神經(jīng)活動，揭示相關(guān)神經(jīng)活動的網(wǎng)絡(luò)。基于功能性磁共振成像獲得的大腦數(shù)據(jù)，已經(jīng)提出相當(dāng)多的腦功能連接建模方法，包括基于相關(guān)的方法[1]、基于偏相關(guān)的方法[2]和圖形建模方法[3]。然而，基于相關(guān)的方法僅僅能夠捕獲成對的信息，因此不能全面反映多個腦區(qū)之間的交互。此外，基于相關(guān)的網(wǎng)絡(luò)由于任意選取閾值有許多虛假的連接[4]。偏相關(guān)估計通常是通過使用逆協(xié)方差矩陣的最大似然估計(MLE)實現(xiàn)，而可靠的估計需要的數(shù)據(jù)樣本規(guī)模要比建模的大腦區(qū)域數(shù)量大得多[5]。圖形化的模型被用來研究大腦連接時缺少先驗知識。

高階信息對于疾病診斷可能是重要的，因為最近的神經(jīng)科學(xué)研究認定在神經(jīng)元同位素示蹤、局部場電位和皮層活動中有重要的高階交互[6]。為了將高階信息應(yīng)用于腦功能網(wǎng)絡(luò)研究，提出了超網(wǎng)絡(luò)構(gòu)建方法[7]。超網(wǎng)絡(luò)表示一種網(wǎng)絡(luò)，它的每一條邊代表多個腦區(qū)之間的交互作用。在已有文獻中，腦功能超網(wǎng)絡(luò)是使用稀疏線性回歸方法構(gòu)建；其中，求解稀疏線性回歸采用的是LASSO(least absolute shrinkage and selection operator)方法。即使LASSO方法已成功應(yīng)用于許多研究，它也存在局限：在超邊構(gòu)建時，選定一個腦區(qū)后，如果其它腦區(qū)之間存在較強的相關(guān)性，那么選擇與選定有關(guān)的腦區(qū)時往往只隨意選擇存在組效應(yīng)的一組腦區(qū)中的其中一個[8]，可能還有一些相關(guān)的腦區(qū)無法選擇出來，缺少解釋分組效應(yīng)信息的能力。本文考慮到腦區(qū)之間的組效應(yīng)，提出將Group Lasso方法引入到超網(wǎng)絡(luò)構(gòu)建中，對超網(wǎng)絡(luò)構(gòu)建方法進行改善，然后根據(jù)構(gòu)建的超網(wǎng)絡(luò)提取3種關(guān)于大腦區(qū)域特定的特征(3種不同定義的聚類系數(shù))。此外，利用非參數(shù)置換檢驗從3組聚類系數(shù)中選擇最具差異的特征。最后，使用多核支持向量機在正常組和自閉癥患者間進行分類模型構(gòu)建。

1 材料和方法

基于超網(wǎng)絡(luò)的腦網(wǎng)絡(luò)分類方法包括幾個主要步驟：數(shù)據(jù)預(yù)處理、超網(wǎng)絡(luò)構(gòu)建、特征提取和選擇以及分類。

1.1 被試

本實驗共收集53名被試，其中有25名自閉癥患者，28名年齡性別匹配的健康志愿者作為對照組。自閉癥患者數(shù)據(jù)來自公開的自閉癥患者腦影像數(shù)據(jù)交換數(shù)據(jù)庫(http://fcon_1000.proj ects.nitrc.org/indi/abide/).正常被試的數(shù)據(jù)采集工作是由山西醫(yī)科大學(xué)第一醫(yī)院進行的，所有的掃描工作由熟悉磁共振操作的放射科醫(yī)生來完成。在掃描的過程中，要求被試閉眼、放松、不去想特定的事物但要保持清醒不能睡著。掃描參數(shù)設(shè)置如下：33 axial slices，repetition time (tR)=2 000 ms，echo time (tE)=30 ms，thickness/skip=4/0 mm，field of view (FOV)=192 mm×192 mm，matrix=64 mm×64 mm，flip angle=90°，248 volumes.

所有被試均被確診為自閉癥患者，通過自閉癥診斷觀察量表[9](autism diagnostic observation schedule，ADOS)及自閉癥診斷訪談量表[10](autism diagnostic interview，ADI)來診斷。同時采用嚴格的排除標(biāo)準(zhǔn)：罹患嚴重的軀體疾病或神經(jīng)系統(tǒng)疾??；妊娠或哺乳期婦女；體格檢查發(fā)現(xiàn)有異常生化指標(biāo)或腦電圖、心電圖異常者。所有被試的基本信息如表1所示。統(tǒng)計分析來評估組間差異，表中a表示雙樣本t檢驗，b表示皮爾遜卡方檢驗。

表1 被試基本情況統(tǒng)計Table 1 Demographics and clinical characteristics of the subjects

1.2 數(shù)據(jù)預(yù)處理

功能數(shù)據(jù)預(yù)處理是利用統(tǒng)計參數(shù)映射軟件包(statistical parametric mapping software package,SPM8)進行。為了保證磁化均衡，每個被試前10個獲得的功能磁共振圖像被丟棄。剩余的圖像進行時間片校正和頭動校正，根據(jù)平移參數(shù)超過3 mm的標(biāo)準(zhǔn)，2例自閉組及2例對照組數(shù)據(jù)被排除。然后，圖像進行12維度的優(yōu)化仿射變換，將其標(biāo)準(zhǔn)化得到MNI(montreal neurological institute)空間標(biāo)準(zhǔn)回波平面成像EPI(echo-planar imaging)模板，重新采樣成3 mm立方體素。由此產(chǎn)生的圖像進行空間平滑并且消除線性趨勢，最后進行低頻濾波(0.01～0.10 Hz)以降低低頻漂移及高頻的生物噪音。

1.3 超圖

圖論是描述計算機科學(xué)中許多問題和結(jié)構(gòu)的有力工具，已被廣泛應(yīng)用于腦網(wǎng)絡(luò)分析[11]。大腦連接可以被簡化為一個圖的節(jié)點以及它們相互連接的邊。節(jié)點和邊的關(guān)系定義了網(wǎng)絡(luò)的拓撲結(jié)構(gòu)，可以通過描述網(wǎng)絡(luò)結(jié)構(gòu)的局部和全局屬性來進行分析。然而，圖只能描述一些二元關(guān)系，并不足以模擬一些復(fù)雜的問題或數(shù)據(jù)。事實上，除了兩兩之間的關(guān)系，在很多應(yīng)用中，可能存在高階的關(guān)系，不能用傳統(tǒng)的圖表示。為了克服這個局限，超圖被提出用來描述節(jié)點之間的高階關(guān)系。超圖是圖的擴展，它的超邊是關(guān)于頂點的任意子集。這個概念模型比圖論更符合一般的關(guān)系類型，并且已被應(yīng)用于化學(xué)、工程和圖像處理等許多領(lǐng)域。

一個超圖H=(V，E)，節(jié)點集合V，超邊集合E，可以使用|V|×|E|維鄰接矩陣表示H：

(1)

式中：v∈V是一個節(jié)點；e?E是H的一條超邊。

基于H，每個頂點v∈V的節(jié)點度是：

(2)

超邊e的邊度是：

(3)

1.4 超網(wǎng)絡(luò)構(gòu)建

每個被試的大腦空間通過自動解剖標(biāo)記[12]模板進一步分割成90個感興趣區(qū)域(region of interest,ROI).每個區(qū)域作為一個網(wǎng)絡(luò)節(jié)點，同時該區(qū)域內(nèi)所有體素的時間序列的平均值作為該節(jié)點的時間序列。每個腦區(qū)的平均時間序列進行了回歸分析，以排除平均腦脊髓液和白質(zhì)信號以及頭動校正對信號的影響。以ROI作為網(wǎng)絡(luò)中的節(jié)點，根據(jù)R-fMRI時間序列使用稀疏線性回歸方法構(gòu)建超網(wǎng)絡(luò)，X=[x1,…,xm,…,xM]T∈RM×d表示一個訓(xùn)練被試，有M個ROIs，xm表示第m個感興趣區(qū)域的平均時間序列，d是時間序列的長度。每個感興趣區(qū)域的時間序列被視為一個響應(yīng)向量，可以利用其他M-1個感興趣區(qū)域的時間序列的線性組合估計：

xm=Amαm+τm.

(4)

式中：Am=[x1,…,xm-1,0,xm+1,…,xM]表示包含除了第m個ROI之外的其它ROIs的時間序列的數(shù)據(jù)矩陣；αm表示權(quán)重向量，衡量其它ROIs對第m個ROI的影響程度，非零元素表示相應(yīng)ROIs與第m個ROI相互作用；τm表示噪聲項。

1.4.1 傳統(tǒng)基于Lasso方法的超網(wǎng)絡(luò)構(gòu)建

腦功能超網(wǎng)絡(luò)是使用Lasso方法求解稀疏線性回歸模型進行構(gòu)建，這是l0范數(shù)問題，其優(yōu)化目標(biāo)函數(shù)為：

(5)

l0范數(shù)問題是一個NP問題，可以轉(zhuǎn)化為l1范數(shù)問題求解。

(6)

式中：λ是控制模型稀疏的正則化參數(shù)。不同的λ值對應(yīng)于不同的稀疏性的解決方案，較大的λ值表明更為稀疏的模型，即在αm中有較多的零。

在實驗中，研究不同腦區(qū)之間的相互作用，對每一個被試構(gòu)建超網(wǎng)絡(luò)，以ROI作為節(jié)點，超邊包括第m個ROI和其他在αm權(quán)重向量中非零元素對應(yīng)的ROIs.為了反映大腦區(qū)域之間信息多層次的相互作用，對每一個ROI，通過在一個特定范圍內(nèi)變化λ值產(chǎn)生一組超邊。在這里，多層次是指不同的λ值確定不同腦區(qū)之間的相互作用關(guān)系。也就是說，λ值較大的目標(biāo)函數(shù)產(chǎn)生一個更為稀疏的解，因此超邊包含更少的節(jié)點。具體地，在實驗中，為了簡單起見，改變λ值從0.1到0.9，增量為0.1.

1.4.2 基于Group Lasso方法的超網(wǎng)絡(luò)構(gòu)建

即使Lasso方法已成功應(yīng)用于許多情況，它仍然存在一些局限。在超邊構(gòu)建時，選定一個腦區(qū)后，如果其它腦區(qū)之間存在較強的相關(guān)性，那么選擇與選定有關(guān)的腦區(qū)時往往只隨意選擇存在組效應(yīng)的一組腦區(qū)中的其中一個，而不關(guān)心是哪一個；還有一些相關(guān)的腦區(qū)無法選擇出來，缺少解釋分組效應(yīng)信息的能力。

使用聚類方法將相關(guān)性強的腦區(qū)分為一組，再使用Group Lasso方法進行超邊的構(gòu)建可以幫助解決腦區(qū)之間的組效應(yīng)問題。Lasso是用來選擇單個變量[13]，Group Lasso可以用來選擇組變量，是在預(yù)先定義的變量組的基礎(chǔ)上進行變量選擇[14]。在進行超網(wǎng)絡(luò)構(gòu)建時首先要根據(jù)ROIs的平均時間序列進行聚類獲得90個腦區(qū)的分組關(guān)系。在這里，采用了k中心點聚類法[15]，首先計算腦區(qū)之間的兩兩相似度值，值越大表明兩個樣本越相似，并在此基礎(chǔ)上進行聚類。對90個腦區(qū)進行聚類時，將其劃分為k組，每個組表示一類對象，對象與組之間的關(guān)系必須滿足：1) 每個組至少包含一個對象；2) 每個對象必須屬于一個組。為了盡可能地保證聚類的穩(wěn)定性，在選擇k個初始化聚類中心時沿用k-means++[16]的思想，隨機選擇一個點作為第一個初始聚類中心，隨后的每一個初始聚類中心是從剩余的數(shù)據(jù)點中以正比于數(shù)據(jù)點與存在的最近聚類中心點的距離的概率隨機選擇。重復(fù)聚類10次選取聚類效果最好的一組作為最終的聚類結(jié)果。在實驗中，k的設(shè)置會影響到網(wǎng)絡(luò)結(jié)構(gòu)以及分類性能。經(jīng)過研究發(fā)現(xiàn)，當(dāng)k等于48時，會得到最高的分類準(zhǔn)確率(詳細分析在討論中會提到)。然后使用Group Lasso選擇腦區(qū)進行超邊的構(gòu)建，以下是優(yōu)化目標(biāo)函數(shù)：

(7)

式中：β是l2，1范數(shù)正則化參數(shù)，不同的β值對應(yīng)不同的稀疏性，β值越大表明模型越稀疏，選擇的組越少；αm通過聚類被分成了k個非重疊的組，αmGi表示第i個組。同樣，為每個被試構(gòu)建超網(wǎng)絡(luò)，以ROI為節(jié)點，根據(jù)αm中非零元素對應(yīng)的ROI構(gòu)建超邊。對每一個ROI，在一定范圍內(nèi)變化β值產(chǎn)生一組超邊，改變β值從0.1到0.9，增量為0.1.

1.5 特征提取和選擇

特征提取和特征選擇是機器學(xué)習(xí)中的關(guān)鍵預(yù)處理步驟。在實踐中，不能知道當(dāng)前的特征是否與問題有關(guān)，不相關(guān)的特征增加了預(yù)測模型的復(fù)雜性,降維方法的應(yīng)用有助于模型的建立和模型預(yù)測性能的改善。在這里，特征提取主要包含3個指標(biāo)的計算，即三個不同定義的聚類系數(shù)[17]，這3種聚類系數(shù)從不同角度衡量了超網(wǎng)絡(luò)的局部屬性。

給定一個超網(wǎng)絡(luò)H=(V，E)，u，t，v表示節(jié)點，e表示超邊，vS={ei∈E∶v∈ei}表示包含節(jié)點v的一系列超邊，vN表示包含節(jié)點v的超邊含有的其它節(jié)點的集合。然后，3種不同類型的聚類系數(shù)在節(jié)點v上可以被定義，分別如下：

(8)

(9)

(10)

如果?ei∈E，例如u,t∈ei，但是v?ei，則I(u,t,v1)=1，否則等于0.如果?ei∈E，例如u,t,v∈ei，則I'(u,t,v)=1，否則等于0.

對于每一種聚類系數(shù)定義，分別從超網(wǎng)絡(luò)中提取作為特征，從而每個被試產(chǎn)生三組特征。特征選擇的目的是從原始特征集合選擇最具代表性的最優(yōu)特征子集。為了選擇與自閉癥病理有關(guān)的關(guān)鍵的特征，采用統(tǒng)計分析方法非參數(shù)置換檢驗來評估自閉癥患者和正常對照之間的差異性，將具有顯著組間差異(p<0.05)的局部屬性作為分類特征進行分類模型構(gòu)建。

1.6 分類

選擇RBF(radial basis function)核函數(shù)以及LOO(leave-one-out)交叉驗證來進行分類以及評估所提出方法的性能。具體而言，一個被試用于測試，其余的用于訓(xùn)練以建立模型。對于每個被試，整個過程重復(fù)進行并且選取分類準(zhǔn)確率的算術(shù)平均值作為最終的分類結(jié)果。多核分類方法的參數(shù)是基于訓(xùn)練被試的網(wǎng)格搜索(范圍從-8到8，步長為1)，通過k折交叉驗證方法得出訓(xùn)練集驗證準(zhǔn)確率并進行比較確定。此外，對于每一種類型的聚類系數(shù)，共有90個特征從所構(gòu)建的超網(wǎng)絡(luò)中提取。對于每個提取的特征，從訓(xùn)練被試中計算平均值和標(biāo)準(zhǔn)偏差進行標(biāo)準(zhǔn)化。

2 實驗和結(jié)果

2.1 差異性的腦區(qū)

基于Group Lasso方法，對腦區(qū)的特征進行非參數(shù)置換檢驗所有被試評估自閉癥患者和正常志愿者之間的差異。統(tǒng)計分析結(jié)果表明，與對照組相比，在自閉組中出現(xiàn)顯著異常的腦區(qū)主要在部分邊緣系統(tǒng)區(qū)域(雙側(cè)內(nèi)側(cè)旁扣帶腦回，右側(cè)海馬，右側(cè)海馬旁回，右側(cè)后扣帶回)，部分額葉區(qū)域(左側(cè)三角部額下回，左側(cè)眶部額下回，雙側(cè)中央旁小葉)，以及部分頂葉區(qū)域(右側(cè)角回，左側(cè)楔前葉)等(表3)。選定的特征p值小于0.05(已校驗)，表明患者和正常人之間的差異較大。表2列出了這些異常的腦區(qū)。圖1展示了在模板空間的這些大腦區(qū)域。

表2 基于Group Lasso方法得到的異常腦區(qū)及其顯著性Table 2 Different bode attributes based on group lasso

2.2 分類準(zhǔn)確率

本文提出的基于Group Lasso的超網(wǎng)絡(luò)的腦網(wǎng)絡(luò)分類方法與原有的基于Lasso的超網(wǎng)絡(luò)的腦網(wǎng)絡(luò)分類方法以及傳統(tǒng)的連接網(wǎng)絡(luò)分類方法進行了比較。傳統(tǒng)連接網(wǎng)絡(luò)方法采用皮爾遜相關(guān)，在稀疏度為5%～40%下構(gòu)建所有被試的功能腦網(wǎng)絡(luò)，進行度、中間中心度、節(jié)點效率3個局部指標(biāo)的計算，為表征指標(biāo)在所選閾值空間內(nèi)的整體特性，計算每個指標(biāo)的AUC(area under curve)值，選擇非參數(shù)置換檢驗后具有顯著組間差異的局部屬性的AUC值作為分類特征。同時，也與其它利用腦網(wǎng)絡(luò)特征對自閉癥進行分類的研究進行了比較。分類結(jié)果如表3所示。

圖1 基于Group Lasso方法得到的異常腦區(qū)Fig.1 Discriminative brain regions based on Group Lasso

表3 關(guān)于自閉癥分類方法的分類性能比較Table 3 Comparison of classification potential evaluation of classification methods

3 討論

在基于超網(wǎng)絡(luò)的腦網(wǎng)絡(luò)分類中，網(wǎng)絡(luò)構(gòu)建是非常關(guān)鍵的。在現(xiàn)有研究中已經(jīng)提出了許多功能網(wǎng)絡(luò)模型，但大多數(shù)是基于簡單圖，只反映成對的大腦區(qū)域之間的相互作用關(guān)系。本文基于超圖理論構(gòu)建超網(wǎng)絡(luò)模型，用于描述多個腦區(qū)之間的高階相互作用。這種高階相互關(guān)系可能包含有用的信息用于識別患者和正常人，并且在此基礎(chǔ)上提出了新的超網(wǎng)絡(luò)構(gòu)建方法。研究結(jié)果表明，基于超網(wǎng)絡(luò)的腦網(wǎng)絡(luò)分類方法可以提高分類性能，并且與原有超網(wǎng)絡(luò)構(gòu)建方法相比，基于Group Lasso的超網(wǎng)絡(luò)構(gòu)建方法可以實現(xiàn)更好的分類結(jié)果。同時，與其它文獻中的分類結(jié)果相比也可以得到一定的改善。

基于Group Lasso的超網(wǎng)絡(luò)構(gòu)建方法，經(jīng)過統(tǒng)計分析可以得到一些具有組間差異的腦區(qū)。這些區(qū)域主要集中在邊緣系統(tǒng)、額葉以及頂葉，其中，海馬、后扣帶回、角回等區(qū)域均為默認網(wǎng)絡(luò)關(guān)鍵區(qū)域。在自閉癥的病理研究中，默認網(wǎng)絡(luò)被廣泛認為是自閉癥的主要病理環(huán)路[22]。

本文提出的方法會受到一些參數(shù)的影響，其中，參數(shù)k是Group Lasso方法中進行聚類的組數(shù)，選取不同的k會得到不同的網(wǎng)絡(luò)結(jié)構(gòu)以及分類結(jié)果。為了比較k對于分類性能的影響，設(shè)置k的變化范圍為[6，90]，步長大小為6.由于第一個初始種子點的隨機選擇會造成結(jié)果的差異，分別在每一個k值下進行50次實驗，選取正確率的算術(shù)平均值作為最后的分類結(jié)果。圖2展示了實驗結(jié)果，結(jié)果顯示當(dāng)k=48時，最高正確率達到87.84%.

圖2 不同k值對應(yīng)的分類準(zhǔn)確率Fig.2 Classification accuracy of different k values

如圖2所示，當(dāng)k值較小或較大時，即網(wǎng)絡(luò)構(gòu)建約束過于緊張或過于寬松時，分類正確率都比較低。這一結(jié)果表明，適中的連接構(gòu)建約束可以得到更為有效的分類結(jié)果，而過于嚴格或?qū)捤傻臉?gòu)建策略，均無法達到滿意的效果。

4 結(jié)束語

本文將稀疏線性回歸模型用于構(gòu)建超網(wǎng)絡(luò)，計算相關(guān)指標(biāo)并選擇具有顯著組間差異的指標(biāo)作為分類特征，利用腦區(qū)之間的高階關(guān)系來進行自閉癥患者與正常人的分類。實驗證明，基于超圖的腦網(wǎng)絡(luò)分類方法不僅可以改善大腦的疾病分類，也便于與疾病有關(guān)的結(jié)構(gòu)的檢測。此外，基于Group Lasso的超網(wǎng)絡(luò)構(gòu)建方法可以較好地解決分組效應(yīng)，選取相關(guān)程度比較大的腦區(qū)，改善分類性能。

在目前的研究中，基于Group Lasso的方法由于第一個聚類初始種子點的隨機選取以及聚類數(shù)k的不同會造成網(wǎng)絡(luò)結(jié)構(gòu)以及分類結(jié)果的不唯一。接下來，如何建立更加穩(wěn)定的超邊則是后續(xù)工作的重點。