劉永艷，聞敏，李瑤，Ibegbu Nnamdi Julian，郭浩

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

抑郁癥(major depression disorder，MDD)是一種嚴重的精神疾病，主要癥狀表現(xiàn)為情緒低落，悲觀，思維遲緩，缺乏主動性，飲食、睡眠差，嚴重者可出現(xiàn)自殺念頭和行為。影響著世界上超過6%的人口。同時，MDD在所有已報告的腦疾病中占比高達98%。盡管MDD在醫(yī)學領域內(nèi)多有研究，但對其臨床診斷和發(fā)病機制的研究還沒有取得較大的突破。

幸運的是，神經(jīng)影像技術的出現(xiàn)幫助人們進一步了解了人類大腦，為腦功能網(wǎng)絡研究提供了堅實的基礎。功能磁共振成像(functional magnetic resonance imaging,fMRI)[1]是一種新型的神經(jīng)影像技術，已成功應用在不同的腦疾病診斷中[2]?；趂MRI獲得的數(shù)據(jù)，人們提出了功能性大腦連接模型的分析方法[3]。傳統(tǒng)的腦網(wǎng)絡構建方法，如Guo等[4]提出了圖形化模型，Wang等[5]提出了基于相關的方法，Akin等[6]提出了基于偏相關的方法，Yu等[7]提出了稀疏表示方法，這些方法只能獲取兩兩相關的數(shù)據(jù)，不能完全反映大腦中多個區(qū)域之間的相互作用[1]。最近的研究表明，在神經(jīng)元同位素示蹤、皮層活動和局部場電位中存在更高階的相互作用[8]。神經(jīng)科學研究表明，在同一時間內(nèi)，神經(jīng)元活動過程中，一個大腦區(qū)域與其他多個大腦區(qū)域之間存在直接相互作用。這意味著大腦區(qū)域之間的高階關系被忽略了，而這些關系可能對大腦和相關疾病的研究是重要的。且腦網(wǎng)絡在本質上是動態(tài)的，相互作用的神經(jīng)元之間的動態(tài)變化可能正向影響相關功能連通性的拓撲結構和關聯(lián)強度，而這些微妙的變化可能是腦疾病發(fā)作的誘因[9]。

為了解決上述問題，Jie等[10]利用稀疏表示從靜態(tài)時間序列構建了超網(wǎng)絡。超網(wǎng)絡可以有效表現(xiàn)出人腦這一復雜系統(tǒng)中信息交互的空間層次中的多元高階關系。在現(xiàn)有的超網(wǎng)絡構建方法中，重點僅在于腦區(qū)之間的相互作用，從而存在一定的局限性，因為這種相互作用在本質上是靜態(tài)的，而功能連接的時變特性可能包含非常多的信息，從而可能無法從圖中獲取有用信息[11]。劉磊等[12]提出的高階功能連接網(wǎng)絡可以解決腦網(wǎng)絡的動態(tài)變化問題，其可體現(xiàn)出時間層次的動態(tài)變化。通過劃分時間窗的方法來構造動態(tài)低階功能連接網(wǎng)絡，堆疊所有時間窗下的動態(tài)低階功能連接網(wǎng)絡，進而計算兩兩動態(tài)時間低階功能連接網(wǎng)絡間的皮爾遜相關系數(shù)，得到高階功能連接網(wǎng)絡。此外，一些研究者嘗試將兩種網(wǎng)絡進行融合，Guo等[13]提出了一個基于感興趣區(qū)域的高階網(wǎng)絡(空間層次)，并將高階有網(wǎng)絡與低階網(wǎng)絡進行融合并分類。但是這些方法都沒有同時體現(xiàn)腦網(wǎng)絡的時間層次的動態(tài)性和空間層次的高階關系。

因此，現(xiàn)提出一種具有時變特性的超網(wǎng)絡構建方法，即融合高階功能連接網(wǎng)絡和超網(wǎng)絡的特性。使所構建的網(wǎng)絡不僅能有效反映空間中具有多重關系的多個腦區(qū)之間的相互作用，還能體現(xiàn)出時變特性對網(wǎng)絡的動態(tài)影響。

首先，使用獨立成分分析(independent component analysis,ICA)[14]模板構建低階功能連接(functional connection,FC)網(wǎng)絡[15]。獨立成分分析模板是一種多元數(shù)據(jù)驅動的分析方法，不需要預設模板。為每個被試構建功能連接網(wǎng)絡并將其堆疊在一起。每個被試的功能連接網(wǎng)絡的堆疊揭示網(wǎng)絡隨時間的動態(tài)特性。然后，在每個被試的堆疊的低階功能連接網(wǎng)絡基礎上，構建超網(wǎng)絡。在這項研究中，使用稀疏線性回歸模型來構建超網(wǎng)絡，并由最小絕對值收斂和選擇算子(least absolute shrinkage and selection operator，LASSO)方法進行求解。最后，根據(jù)超網(wǎng)絡的聚類系數(shù)，提取特征，通過非參數(shù)置換檢驗方法選擇具有顯著差異的特征。使用支持向量機(support vector machine,SVM)[16]對差異性特征進行分類。

1 材料與方法

1.1 方法步驟

基于時變特性的MDD患者的超網(wǎng)絡構建及分類方法包括以下5個步驟。

步驟1數(shù)據(jù)采集及預處理。采集MDD被試和正常被試的功能磁共振成像數(shù)據(jù)進行預處理。

步驟2組獨立成分分析(group independent analysis,GICA)。主要包括數(shù)據(jù)降維、獨立成分估計、數(shù)據(jù)重構以及去除噪聲。

步驟3低階功能連接網(wǎng)絡的構建。根據(jù)模板劃分成分，提取每個成分的平均時間序列。選擇固定的滑動窗口分割每個成分的平均時間序列，計算兩兩時間序列間的皮爾遜相關系數(shù)，構建低階功能連接網(wǎng)絡。將所有的低階功能連接網(wǎng)絡進行堆疊，即對每個時間窗口下低階功能矩陣中的相應元素的值進行提取，形成時變的低階功能連接網(wǎng)絡。

步驟4超網(wǎng)絡的構建。采用LASSO方法求解稀疏線性回歸模型，根據(jù)上述產(chǎn)生的時變的低階功能連接網(wǎng)絡，使用該模型構建每個被試所對應的超網(wǎng)絡。

步驟5特征提取、選擇和分類。構建超網(wǎng)絡后，引入3種不同定義的聚類系數(shù)作為超網(wǎng)絡的局部特征。以統(tǒng)計差異分析為基礎，選取最具判別性的特征作為分類的關鍵特征。即采用Kolmogorov &Smimov(KS)[17]非參數(shù)置換檢驗，并從提取的特征中選擇顯著性特征。使用支持向量機進行分類。在分類過程中，采用徑向基函數(shù)(radical basis function，RBF)和留一交叉驗證方法[18]進行有效性評估。

1.2 數(shù)據(jù)采集和預處理

嚴格按照山西醫(yī)學倫理委員會的要求(會議號：2012013)，在該項研究實施之前，與每位參與者均達成了書面協(xié)議。根據(jù)《赫爾辛基宣言》，書面知情同意書由實驗中的每一個被試簽署?？偣舱心剂?6名被試，其中包括38名首發(fā)，無用藥重度MDD患者(15名男性)和28名健康右利手志愿者(13名男性)。靜息狀態(tài)下，應用3T磁共振掃描儀(Siemens Trio 3-Tesla scanner,Siemens,Erlangen,Germany)對他們進行功能磁共振成像(fMRI)掃描。被試的具體情況如表1所示。

表1 被試的具體信息Table 1 Specific information of participants

山西醫(yī)科大學第一醫(yī)院完成數(shù)據(jù)采集任務，并且由精通磁共振技術的放射科醫(yī)師完成掃描任務。被試在進行掃描時也有相應的要求，即被試需要在放松的狀態(tài)下閉上眼睛，但不能入睡，需要保持清醒，也不進行特定的思考。每個掃描的結果是248個連續(xù)的EPI功能圖像，其中的掃描參數(shù)有如下設置：軸向切片33,回波時間TE=30 ms,重復時間TR=2 000 ms,厚度/跳過=4/0 mm,偏轉角=90°，矩陣=64 mm×64 mm,視野FOV=192 mm×192 mm。前10個功能圖像的時間序列由于被試對環(huán)境的自適應性以及初始磁共振信號的不穩(wěn)定性而被丟棄。

用SPM8(http://www.fil.ion.ucl.ac.uk/spm)來完成數(shù)據(jù)的預處理過程。首先，頭動校正和時間片校正必不可少，而在校正過程中，抑郁組和對照組中分別出現(xiàn)2例轉動大于3°或者頭動大于3 min的被試，因此棄除這些被試的掃描數(shù)據(jù)。需要注意的是，最后的66例樣本數(shù)據(jù)中不包含那些被丟棄的數(shù)據(jù)。接著，圖像經(jīng)過優(yōu)化仿射變換會被標準化到蒙特利爾神經(jīng)研究所(Montreal neurological institute，MNI)的標準空間中。最后，為降低生物高頻噪音和低頻漂移的影響，對數(shù)據(jù)進行帶通濾波(0.01～0.10 Hz)和線性降維。

1.3 組獨立成分分析

獨立成分分析(ICA)已成為處理和分析fMRI數(shù)據(jù)的重要方法之一。ICA是一種多元數(shù)據(jù)驅動的分析方法，主要優(yōu)點之一是它不需要預設模板。從fMRI數(shù)據(jù)分析中找到不同被試之間的一致模式很重要，但是被試之間的個體差異會導致從不同的被試中獲得的成分彼此不一致，并且難以組合這些成分進行分析。通過將所有被試的fMRI數(shù)據(jù)鏈接到ICA估計中，提出組獨立成分分析(GICA)方法，找出在被試之間表現(xiàn)出較高一致性的通用成分。這樣就可以將所得的獨立成分用于計算統(tǒng)計顯著性。

研究采用GIFT(http://mialab.mrn.org/software/gift)。具體來說，使用最小描述長度(minimum description length,MDL)[19]來估計兩組中最優(yōu)的分解次數(shù)。MDL最終被設置為54。為保證獨立成分的穩(wěn)定性和可靠性，對每個fMRI實例和54個空間獨立成分使用信息最大化(information maximization,Infomax)算法，在ICASSO上隨機化初始分解矩陣，重復20遍Infomax算法，得到相同的收斂閾值。最后，使用GICA3算法反轉數(shù)據(jù)，并對其進行重構，獲得被試獨立成分的空間分布和時間序列。

經(jīng)過GICA處理后，提取到很多成分，但其中一些是不需要的，而某些成分還包含了較大的噪音。為了選擇重要的成分，采用了匹配法和人工檢查法來確定所篩選出來的成分[20]。排除成分的標準包括較大的激活區(qū)域，其中的多元回歸系數(shù)與先前模板匹配；主要激活區(qū)域在灰質中的分布，以及這些區(qū)域與已知成分(例如血管和低頻空間中的頭部運動)的重疊。低頻功率在激活區(qū)域中對時間序列功率譜的控制[21]。

經(jīng)過篩選，共獲得22個成分。發(fā)現(xiàn)這些成分是聽覺網(wǎng)絡、感覺運動網(wǎng)絡、視覺網(wǎng)絡、默認模式網(wǎng)絡、注意力網(wǎng)絡或額葉網(wǎng)絡的一部分。如圖1所示，顯示了屬于上述各個網(wǎng)絡的成分。

圖1 成分及其各自網(wǎng)絡的空間地圖Fig.1 Spatial maps of components and their respective networks

1.4 網(wǎng)絡構建

1.4.1 動態(tài)低階功能連接網(wǎng)絡的構建

使用ICA模板，大腦被分為22個(左右半腦各11個)成分。計算每個成分中體素的血氧水平依賴(blood oxygen level-dependent,BOLD)信號的平均值，這代表每個節(jié)點的信號值。通過提取所有體素在不同點的BOLD信號值并求平均值，得到22個成分的平均時間序列。選取一個固定的窗口大小，對提取的平均時間序列進行分割[22]。

采用滑動窗口法將時間序列劃分為重疊子序列的片段，每個片段代表較短時間內(nèi)的序列。假設rs-fMRI的平均時間序列是針對第l個被試的第m個感興趣成分，K個重疊段的表達式為

K=[(M-N)/s+1]

(1)

對于第l個被試的R個成分，第k段的矩陣形式表示為

(2)

(3)

(4)

1.4.2 具有時變特性的超網(wǎng)絡的構建

超網(wǎng)絡的構建是將上述矩陣Y(l)中的每一個時間序列分別作為新的節(jié)點集，并利用稀疏線性回歸模型[23]來構建超網(wǎng)絡。使用該模型，可以獲得Y(l)中所有時間序列兩兩之間的交互。

稀疏線性回歸模型采用的公式為

Xi=Aiαi+τi

(5)

式(5)中：Xi為指定Y(l)的第i行；Ai=[X1,…,Xi-1,0,Xi+1,…,Xu]為包含Y(l)中除第i行以外的其他所有行的數(shù)據(jù)矩陣；αi為權重向量，表示其他行時間序列對第i行的影響程度；τi為噪聲項。

采用LASSO方法對稀疏線性回歸模型進行處理和計算，優(yōu)化函數(shù)如下。

min‖Xi-Aiαi‖2+λ‖αi‖1

(6)

式(6)中：‖·‖1為L1范式；‖·‖2為L2范式；λ為調(diào)節(jié)模型稀疏度的參數(shù)，該稀疏度代表不同稀疏解。λ的取值是一個范圍，取值越大表示αi中有更多的零值。模型具有的零值數(shù)目越高，模型將越稀疏。λ值的范圍為0.1～0.9。此外，為了解決優(yōu)化問題，使用SLEP軟件包[24]。

每個被試可以構建得到一個超網(wǎng)絡，其中一個時間序列被認為是一個節(jié)點，超網(wǎng)絡中的一條超邊ei不僅考慮到中心時間序列(第i個時間序列)，還考慮了用式(6)計算得到權重αi中其他非零權重對應的時間序列。由于時間序列之間可能存在多層次交互作用，可以通過在0.1～0.9范圍內(nèi)給出不同的λ，為每個時間序列生成一組超邊。具體來說，對于每個被試，總共有u個時間序列，其中u=[R(R-1)]/2，R表示每個被試中成分的個數(shù)。假如設置num個不同的λ，則可以生成得到的表示超網(wǎng)絡連接關系的矩陣A，其中A∈Ru×u×num，num設置為9。這個最終的矩陣就是用LASSO方法構造的超網(wǎng)絡[25]。

1.5 特征提取

構建功能連接超網(wǎng)絡之后，選擇超網(wǎng)絡中每個頂點的屬性值作為特征。對已構建的超網(wǎng)絡進行指標計算，定義了3種聚類系數(shù)，分別具有不同的定義[26]，分別記作HCC1、HCC2以及HCC3。這3種聚類系數(shù)從不同的角度反映了超網(wǎng)絡的局部聚類屬性。

(7)

(8)

(9)

式中：u、t、v為節(jié)點；N(v)為集合，指超邊中除了節(jié)點v以外，還包括其他節(jié)點的集合；|e|為超邊中包含的節(jié)點數(shù)目；S(v)為超邊集合，這些超邊中都含有節(jié)點v，且S(v)={ei∈E,v∈ei}。

式(7)計算與節(jié)點v不存在連接的相鄰節(jié)點的數(shù)目。如果E表示邊集，ei表示某一條超邊，則I(u,t,-v)=1當且僅當?ei∈E,u,t∈ei，但v?ei；否則I(u,t,-v)=0。HCC1查找不包含u的鄰居之間的連接，HCC1的優(yōu)點是，在這個集合中發(fā)現(xiàn)的任何交互都可能表示鄰居之間的真實連接。它的缺點在于可能過分關注于那些次要的共享連接，這些連接與u的交互沒有什么關系。

式(8)計算與節(jié)點v存在連接關系的相鄰節(jié)點的數(shù)目。如果?ei∈E,u,t∈ei，則I′(u,t,-v)=1。HCC2查找那些包含u的鄰居的連接，這種方式找到的邊真實地反映了u和鄰居之間的聚集。但是需要注意這種連接可能只是與u共享連接的人工數(shù)據(jù)。

式(9)計算超邊之間的重疊量，具體指節(jié)點v的相鄰超邊。|e|表示超邊中包含的節(jié)點數(shù)目；S(v)指超邊集合，這些超邊中都含有節(jié)點v，且S(v)={ei∈E,v∈ei}。其分子表示與u相關聯(lián)的超邊的頂點個數(shù)。分母表示這種重疊可能的數(shù)目。HCC3則通過鄰域超邊的重疊量來衡量鄰域的密度。與以上兩種定義都不同，它從節(jié)點的角度來定義重疊量。

對于HCC1、HCC2以及HCC3，從超網(wǎng)絡中提取一組聚類系數(shù)作為特征，因此每個被試對應3組特征。

1.6 特征選擇和分類

從超網(wǎng)絡中提取了多個特征，但是其中一些提取的特征是多余的。以統(tǒng)計差異分析為基礎，選擇最具判別性的特征作為分類的關鍵特征。計算每個被試所對應的平均聚類系數(shù)(HCC1、HCC2以及HCC3，分別對231個腦區(qū)進行平均)，使用KS非參數(shù)置換檢驗來分析兩組特征是否具有差異。用KS檢驗比較正常對照組和抑郁組的所有693個特征及其性質，選取局部屬性(P<0.05)(聚類系數(shù))。將所選擇的屬性作為分類特征，用于建立分類模型。

將N個特征中的每個特征都分別作為測試集，其他的N-1個特征作為訓練集。訓練集使用K折交叉驗證(K-fold cross validation,K-CV)[27]對懲罰因子c和核參數(shù)g進行尋優(yōu)。在測試集驗證中，c和g的最優(yōu)值取分類準確率最高的那組值，并由此建立N個不同的模型。在標準化分類特征的平均值和標準偏差后，最后的分類結果取不同模型的分類準確率的平均值。

2 結果與分析

2.1 異常成分

經(jīng)過特征選擇和特征提取，選出了27個顯著的功能連接。選擇這些連接的原因在于，和這些連接有關的成分記錄了與其他成分的最大連接數(shù)。在22個成分中，IC25、IC33和IC43的出現(xiàn)次數(shù)最高，各7次，IC39出現(xiàn)6次。這些成分共有27個功能連接，如表2所示。

表2 27個功能連接對應的數(shù)據(jù)信息Table 2 Data information corresponding to 27 functional connections

2.2 分類結果

研究的主要目的是創(chuàng)建一個與多個成分交互的超網(wǎng)絡，并反映神經(jīng)交互作用隨時間的動態(tài)變化。

對于MDD的分類結果，對比了不同的方法，包括稀疏時間動態(tài)網(wǎng)絡(sparse temporally dynamic networks，DNwee)方法，動態(tài)高階網(wǎng)絡(dynamic high-order network，DNH)方法、分層高階功能連接網(wǎng)絡(Hierarchical high-order FC networks，HHON-SFS)方法和基于時變的超網(wǎng)絡構建方法。如表3所示，可以看出，所提方法具有明顯更高的準確率，可以更好地對MDD進行分類。

此外，不同研究采用的數(shù)據(jù)集與方法都存在差異，因此，為了評估本文方法的有效性，還在同一數(shù)據(jù)集下構建了超網(wǎng)絡(hyper-network,HN)、低階功能連接網(wǎng)絡(low-order function connection networks,LON)以及高階功能連接網(wǎng)絡(high-order function connection networks,HON)，并進行了對比。如表4所示，顯示了具有時變特性的超網(wǎng)絡與上述網(wǎng)絡在分類精度上的顯著優(yōu)勢。

表4 同一數(shù)據(jù)集下不同構建方法的分類結果Table 4 Classification results of different network construction methods with the same data set

此外，經(jīng)過特征提取和選擇過程后，分類準確率為86.36%，特異度為92.10%，靈敏度為78.57%。

3 方法論

具有時變特性的超網(wǎng)絡構建和所選分類方法的性能在很大程度上依賴于一組特定參數(shù)的選擇，如滑動窗口長度(大小)、步長、超網(wǎng)絡構建模型參數(shù)(λ)和SVM分類模型參數(shù)(c和g)。這些參數(shù)值的變化會顯著影響我們得到的結果。為了得到最佳的分類結果，必須為每個參數(shù)選擇最優(yōu)值。另外，為了獲得盡可能高的分類精度，選擇了不同的參數(shù)值組合。

3.1 滑動窗口長度

滑動窗口長度的變化會改變時間窗口的數(shù)量，這也將影響為每個被試構建的低階功能連接網(wǎng)絡的數(shù)量。在其他參數(shù)保持不變，步長固定為2 s的情況下，選取窗口長度分別為40、50、60、70、80、90 s，研究各自的結果。當窗口長度設為60 s時，得到了最優(yōu)結果。值得注意的是，當滑動窗口長度過大或過小時，該方法的精度顯著降低。與使用不同窗口長度的高階網(wǎng)絡相比，本文方法的分類結果如圖2所示。

圖2 不同滑動窗口長度的分類結果Fig.2 Classification results of different sliding window lengths

3.2 正則化參數(shù)的影響

λ是超網(wǎng)絡建設中最重要的參數(shù)之一，因為它控制著網(wǎng)絡的稀疏性。網(wǎng)絡拓撲和模塊性高度依賴于稀疏性。λ太大會導致網(wǎng)絡稀疏，而太小會導致網(wǎng)絡過于粗糙甚至可能包含噪聲[30]。這些都將導致網(wǎng)絡的缺陷。因此，重要的是找到λ的最佳值以給出最佳結果。結果表明，當λ設置為0.1時，網(wǎng)絡是最可靠的。λ設置為0.1意味著網(wǎng)絡中每個節(jié)點都至少在一條超邊上[31]。先前的研究[10]表明，多級λ設置是最好的方法，因為它可以組合多個λ以提供有關網(wǎng)絡的更多信息。使用一系列值升序排列組合表示λ，即為{0.1}，{0.1,0.2}，…，{0.1,0.2,…,0.9}。較小和較大的λ被舍棄，這證明了節(jié)點是在超網(wǎng)絡中連接的，使得特征提取更加具體。

3.3 SVM分類參數(shù)的影響

SVM分類器在各領域中的應用都比較廣泛，而分類時經(jīng)常涉及到核函數(shù)的選取問題。由于RBF核函數(shù)應用廣泛，無論是小樣本還是大樣本均試用。因此，在分類中選取RBF核函數(shù)。SVM模型中有兩個參數(shù)對分類影響較大，即懲罰因子c和核參數(shù)g：c表示調(diào)節(jié)優(yōu)化方向中兩個指標(間隔大小，分類準確度)偏好的權重，即對誤差的寬容度，c過大或過小，其泛化能力都會變差；g隱含地決定了數(shù)據(jù)映射到新的特征空間后的分布，g越大，支持向量的個數(shù)越少，反之則越多，而支持向量的個數(shù)影響訓練與預測的速度。最優(yōu)的c和g能使SVM的分類性能達到最佳。對于如何找出最佳的c和g，首先選取某一組給定的c和g的值，把訓練集作為原始數(shù)據(jù)集，利用K-CV方法得到在該組c和g下的訓練集驗證的分類準確率，然后不斷更換c、g的值，最終取使測試集驗證分類準確率最高的那組c和g作為最佳參數(shù)。參數(shù)c、g可在[2-8，28]內(nèi)變化，并將1設置為步長大小。圖3表示參數(shù)c和g的尋優(yōu)結果，結果顯示，測試集驗證分類準確率最高時的c、g參數(shù)的值分別為2和0.062 5，最高準確率是91.923 1%。

圖3 c和g的參數(shù)尋優(yōu)結果Fig.3 Parameter optimization results ofc and g

4 結論

超網(wǎng)絡已被證明是腦診斷中重要的腦網(wǎng)絡構建方法之一。然而，現(xiàn)有的超網(wǎng)絡構建方法忽略神經(jīng)相互作用在一段時間內(nèi)的動態(tài)變化。為了解決這一問題，提出了一種具有時變特性的超網(wǎng)絡的MDD分類方法。提出的網(wǎng)絡不僅顯示了兩個以上成分(大腦網(wǎng)絡)在同一時間內(nèi)的內(nèi)在相互作用，而且考慮了相互作用之間的時變特性，為分類和診斷提供了大量信息。此外，本文方法具有較高的分類準確率。由此證明，具有時變特性的超網(wǎng)絡是一種更強大的網(wǎng)絡構建方法，可以為MDD和其他腦疾病的分類提供更豐富的信息。