帶時(shí)間約束的Louvain算法在動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化中的應(yīng)用研究*

淡楊超，王 彬,2，薛 潔，盛景業(yè)，劉 暢，詹威威

(1.昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院，云南 昆明 650500；2.昆明理工大學(xué)云南省人工智能重點(diǎn)實(shí)驗(yàn)室，云南 昆明 650500；3.云南省公安廳禁毒局，云南昆明 650228;4.提升政府治理能力大數(shù)據(jù)應(yīng)用技術(shù)國家工程實(shí)驗(yàn)室,貴州 貴陽 550022;5.中電科大數(shù)據(jù)研究院有限公司總體技術(shù)研究中心,貴州 貴陽 550022)

1 引言

功能性磁共振成像技術(shù)fMRI(functional Magnetic Resonance Imaging)是目前一種常用的醫(yī)學(xué)成像技術(shù)[1]，其原理是利用磁振造影測(cè)量神經(jīng)元活動(dòng)引發(fā)的血氧水平變化BOLD(Blood Oxygenation Level Dependent)信號(hào)對(duì)大腦活動(dòng)進(jìn)行檢測(cè)，由于其非侵入性和無輻射暴露等優(yōu)點(diǎn)而在醫(yī)學(xué)方面得到了廣泛應(yīng)用。目前，fMRI主要分為任務(wù)態(tài)ts-fMRI(task state fMRI)和靜息態(tài)rs-fMRI(resting state fMRI)[2]，由于rs-fMRI技術(shù)能夠反映大腦固有的自發(fā)活動(dòng)規(guī)律和連接模式，并且實(shí)驗(yàn)便于實(shí)施，受實(shí)驗(yàn)方案的干擾小，因此很多研究者采用rs-fMRI技術(shù)進(jìn)行動(dòng)態(tài)腦功能網(wǎng)絡(luò)的重構(gòu)[3]。目前大多數(shù)研究者基于rs-fMRI的腦網(wǎng)絡(luò)研究都假定大腦的狀態(tài)在數(shù)據(jù)采集時(shí)段內(nèi)不隨時(shí)間發(fā)生改變，但是隨著研究的深入，很多實(shí)驗(yàn)結(jié)果顯示，隨著時(shí)間的變化，大腦各部分的BOLD信號(hào)會(huì)隨著時(shí)間的變化而發(fā)生強(qiáng)弱變化，反映出大腦的狀態(tài)是隨時(shí)間改變的[4]。

基于rs-fMRI信號(hào)構(gòu)建的功能連接(Functional Connection)網(wǎng)絡(luò)研究發(fā)現(xiàn)，腦功能連接網(wǎng)絡(luò)在時(shí)間序列上是高度可變的[5]，且在靜息態(tài)會(huì)表現(xiàn)出自發(fā)的動(dòng)態(tài)變化[6]，而且對(duì)于不同的任務(wù)，腦功能網(wǎng)絡(luò)會(huì)表現(xiàn)出不同的特征[7]。由于靜態(tài)腦網(wǎng)絡(luò)難以展現(xiàn)出大腦的一些動(dòng)態(tài)任務(wù)特征，因此很多研究者對(duì)大腦網(wǎng)絡(luò)的多層表征形式進(jìn)行了探究[8,9]，針對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)，Shine等人[10]利用rs-fMRI信號(hào)分析了帕金森病(Parkinson's disease)與大腦多巴胺分泌的相關(guān)性。Rukat等人[11]利用腦電數(shù)據(jù)構(gòu)建隱馬爾科夫模型，分析大腦的狀態(tài)在時(shí)間和空間上的動(dòng)態(tài)變化。Tian等人[12]統(tǒng)計(jì)了不同年齡的人其大腦的功能連接處于不同狀態(tài)的持續(xù)時(shí)間長(zhǎng)短，研究年齡與大腦反應(yīng)時(shí)間的關(guān)系，此外還探究了功能連接的狀態(tài)波動(dòng)幅值變化與年齡的關(guān)系。當(dāng)前動(dòng)態(tài)腦功能網(wǎng)絡(luò)已經(jīng)成為了腦科學(xué)研究的主流方向，對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)的網(wǎng)絡(luò)特征進(jìn)行分析有助于了解大腦疾病的誘因，加強(qiáng)對(duì)某些心理疾病的輔助診斷，揭示大腦生理活動(dòng)的運(yùn)行機(jī)理。

上述研究致力于探究動(dòng)態(tài)腦功能網(wǎng)絡(luò)的網(wǎng)絡(luò)特征在時(shí)間和空間的變化規(guī)律，從而發(fā)現(xiàn)動(dòng)態(tài)腦功能網(wǎng)絡(luò)的網(wǎng)絡(luò)特征與大腦疾病之間的相關(guān)性。根據(jù)已有研究[13]發(fā)現(xiàn)，靜態(tài)腦功能網(wǎng)絡(luò)的模塊化程度與大腦的認(rèn)知功能存在強(qiáng)相關(guān)性，在需要不同認(rèn)知功能的任務(wù)下大腦的模塊化程度會(huì)發(fā)生改變。由于動(dòng)態(tài)腦網(wǎng)絡(luò)的模塊化程度也會(huì)隨著時(shí)間的變化而發(fā)生變化，因此很多研究者對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)的模塊化現(xiàn)象展開了分析和研究，Mucha等人[14]對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊化探究，引入了一個(gè)多時(shí)間點(diǎn)耦合參數(shù)，以連接多個(gè)時(shí)間點(diǎn)，從而將相應(yīng)的節(jié)點(diǎn)跨時(shí)間連接，得到一個(gè)新的質(zhì)量評(píng)估函數(shù)來對(duì)網(wǎng)絡(luò)進(jìn)行模塊化劃分。Damicelli等人[15]通過構(gòu)建腦網(wǎng)絡(luò)的拓?fù)湓鰪?qiáng)模型來減小腦網(wǎng)絡(luò)模塊化劃分與真實(shí)腦網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)之間的差距；Khambhati等人[16]研究了動(dòng)態(tài)腦功能網(wǎng)絡(luò)在任務(wù)態(tài)和靜息態(tài)下模塊化程度的變化，發(fā)現(xiàn)在不同任務(wù)下腦網(wǎng)絡(luò)模塊化程度有明顯差別；Baniqued等人[17]的研究揭示了動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化程度與個(gè)人的學(xué)習(xí)記憶能力呈正相關(guān)關(guān)系；Siegel等人[18]研究了中風(fēng)后恢復(fù)過程中動(dòng)態(tài)腦功能網(wǎng)絡(luò)的全局效率變化，并發(fā)現(xiàn)在中風(fēng)恢復(fù)過程中動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化程度與語言恢復(fù)程度密切相關(guān);Han等人[19]對(duì)顱腦損傷后恢復(fù)過程中動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化變化進(jìn)行了研究，發(fā)現(xiàn)在恢復(fù)過程中腦網(wǎng)絡(luò)模塊性降低，參與系數(shù)變大，可根據(jù)動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化程度評(píng)估顱腦損傷恢復(fù)程度；Hilger等人[20]發(fā)現(xiàn)多動(dòng)癥與動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊間和模塊內(nèi)的連接比率存在一定關(guān)聯(lián)，并據(jù)此建立了多動(dòng)癥的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模型，用以對(duì)嬰幼兒多動(dòng)癥進(jìn)行診斷，并使用了Louvain算法[21]對(duì)多動(dòng)癥患者動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分。

Louvain算法的核心思想是以達(dá)到模塊度的最大值為指標(biāo)來劃分模塊，而模塊度最大化面臨的一個(gè)問題是在某些情況下，可能會(huì)檢測(cè)不到低于一定規(guī)模的小模塊，這個(gè)問題在小規(guī)模網(wǎng)絡(luò)的模塊化研究中尤為突出[22]。由于基于rs-fMRI重構(gòu)的動(dòng)態(tài)腦功能網(wǎng)絡(luò)一般是90個(gè)節(jié)點(diǎn)[23]，屬于小規(guī)模網(wǎng)絡(luò)，因此模塊化分辨率對(duì)實(shí)驗(yàn)結(jié)果的影響也更加明顯。針對(duì)Louvain算法無法正確識(shí)別動(dòng)態(tài)腦功能網(wǎng)絡(luò)中小規(guī)模模塊的問題，本文提出了一種帶時(shí)間約束的Louvain算法來提高模塊的分辨率并應(yīng)用該算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊分析。

Figure 1 Processing of dynamic brain function network construction and module division圖1 動(dòng)態(tài)腦功能網(wǎng)絡(luò)構(gòu)建及模塊劃分流程

2 帶時(shí)間約束的Louvain算法動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化分析

2.1 動(dòng)態(tài)腦功能網(wǎng)絡(luò)構(gòu)建及模塊化流程

本文利用AAL(Anatomical Automatic Labeling)模板[23]將大腦分為90個(gè)腦區(qū)，采用大小為W，步長(zhǎng)為L(zhǎng)的滑動(dòng)窗口技術(shù)對(duì)rs-fMRI數(shù)據(jù)進(jìn)行動(dòng)態(tài)特征觀測(cè)窗口的重建，計(jì)算每個(gè)窗口內(nèi)90個(gè)腦區(qū)中每2個(gè)腦區(qū)之間的皮爾遜相關(guān)系數(shù)(Pearson Correlation Coefficient)[24]，并將腦區(qū)間的皮爾遜相關(guān)系數(shù)作為腦區(qū)間的連接強(qiáng)度，構(gòu)建動(dòng)態(tài)腦功能網(wǎng)絡(luò)。圖1所示為動(dòng)態(tài)腦功能網(wǎng)絡(luò)構(gòu)建流程。

模塊性是復(fù)雜網(wǎng)絡(luò)系統(tǒng)的基本特征，網(wǎng)絡(luò)模塊是指一組連接緊密的節(jié)點(diǎn)，通常每一個(gè)模塊是專門處理信息的基礎(chǔ)單元。功能神經(jīng)影像研究表明，人類大腦在大尺度的功能網(wǎng)絡(luò)上具有很好的模塊化結(jié)構(gòu)。一般用模塊度(Modularity)[25]Q表示模塊的劃分程度，Q值的取值為[-0.5,1)[25]，模塊度Q值越大，模塊劃分得到的單一模塊規(guī)模越大，網(wǎng)絡(luò)總的模塊數(shù)量越少。當(dāng)Q值在0.3～0.7時(shí)，網(wǎng)絡(luò)模塊的劃分最好[26,27]。

Louvain算法[21]是一種基于多層次優(yōu)化模塊度的算法，它能夠發(fā)現(xiàn)層次性的模塊，并且模塊辨識(shí)準(zhǔn)確性高、計(jì)算復(fù)雜度低、效率高，因此被廣泛應(yīng)用于模塊化研究。

圖2所示的是使用Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分得出的不同時(shí)間點(diǎn)的模塊度Q值時(shí)間序列，其中，圖2a和圖2c為2個(gè)健康人樣本，圖2b為自閉癥患者樣本。由圖2可以看出，3個(gè)樣本的峰值分別達(dá)到0.82，0.83，081，最小值分別為0.71，0.68，0.55，平均值分別為0.77，0.76,0.71，而模塊度Q值的合理范圍應(yīng)該在0.3～0.7。規(guī)模過大的腦功能模塊可能會(huì)使得一些重要的亞功能腦模塊被忽略，不利于深入研究腦功能網(wǎng)絡(luò)的屬性。

Figure 2 Dynamic brain function network modularity (Q) sequence by Louvain algorithm圖2 動(dòng)態(tài)腦功能網(wǎng)絡(luò)Louvain算法模塊化Q值序列

針對(duì)上述問題，本文提出了一種帶時(shí)間約束的Louvain算法，其原理如圖3所示，圖3a中原Louvain算法要不斷迭代至模塊度增益為0時(shí)才停止，而帶時(shí)間約束的Louvain算法，如圖3b所示，在不斷合并小模塊過程中，采用與時(shí)間相關(guān)的約束使得迭代在適當(dāng)?shù)臅r(shí)候停止，因此動(dòng)態(tài)腦功能網(wǎng)絡(luò)中的一些小規(guī)模模塊也會(huì)被識(shí)別出來。

Figure 3 Comparison of Louvain algorithm and time-constrainted Louvain algorithm圖3 原Louvain算法和帶時(shí)間約束的Louvain算法劃分結(jié)果對(duì)比圖

2.2 帶時(shí)間約束的Louvain算法原理與步驟

2.2.1 帶時(shí)間約束的Louvain算法

首先，針對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)給出如下定義：

定義1根據(jù)Louvain算法，定義動(dòng)態(tài)腦功能網(wǎng)絡(luò)中模塊化劃分優(yōu)劣的判定標(biāo)準(zhǔn)，即模塊度Q為：

(1)

(2)

其中，Ai,j表示腦區(qū)i和腦區(qū)j之間的連接強(qiáng)度之和;m表示整個(gè)動(dòng)態(tài)腦功能網(wǎng)絡(luò)中所有腦區(qū)的連接強(qiáng)度之和;ki=∑jAi,j是連接到i腦區(qū)的所有連接的連接強(qiáng)度之和;δ(ci,cj)是一個(gè)二值函數(shù)，當(dāng)腦區(qū)i和腦區(qū)j被劃分到同一個(gè)模塊時(shí),其值取1，當(dāng)腦區(qū)i和腦區(qū)j不屬于同一個(gè)模塊時(shí)，其值為0。

定義2改進(jìn)的模塊度增益ΔQ：

(3)

其中，ΔQ即為帶時(shí)間約束的Louvain算法中的模塊度增益，用來決定動(dòng)態(tài)腦功能網(wǎng)絡(luò)中腦區(qū)模塊合并的方向；∑in表示被劃分到模塊C內(nèi)的所有腦區(qū)之間的連接強(qiáng)度之和；∑tot表示模塊C內(nèi)所有腦區(qū)與模塊C外的所有腦區(qū)之間的連接強(qiáng)度之和；ki表示腦區(qū)i與其他所有腦區(qū)之間的連接強(qiáng)度之和;ki,in表示腦區(qū)i和模塊C內(nèi)所有腦區(qū)的連接強(qiáng)度之和;m表示動(dòng)態(tài)腦功能網(wǎng)絡(luò)中所有腦區(qū)的連接強(qiáng)度之和。

定義3改進(jìn)的迭代終止系數(shù)η：

(4)

由于原Louvain算法是基于貪婪思想的算法，向著模塊度Q值最大方向合并，在合并過程中會(huì)以損失一些小規(guī)模的模塊為代價(jià)，而動(dòng)態(tài)腦功能網(wǎng)絡(luò)是與時(shí)間密切相關(guān)的，因此本文給出了一個(gè)迭代終止系數(shù)η。該系數(shù)綜合考量了所有時(shí)間點(diǎn)上Q值的平均分布情況，當(dāng)ΔQ=η時(shí)，終止迭代，并將此時(shí)的模塊劃分作為最終的模塊劃分。其中，Qmax_i為原Louvain算法計(jì)算出的n個(gè)時(shí)間點(diǎn)上除時(shí)間點(diǎn)i外所有時(shí)間點(diǎn)最大Q值，Qmin_i為原Louvain算法計(jì)算出的n個(gè)時(shí)間點(diǎn)上除時(shí)間點(diǎn)i外所有時(shí)間點(diǎn)最小Q值。

2.2.2 動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化實(shí)現(xiàn)步驟

步驟1將網(wǎng)絡(luò)中每一個(gè)節(jié)點(diǎn)視為一個(gè)獨(dú)立的模塊，合并相鄰的模塊并計(jì)算合并后的模塊度增益ΔQ，對(duì)網(wǎng)絡(luò)中所有節(jié)點(diǎn)都進(jìn)行鄰居節(jié)點(diǎn)合并，直到模塊度增益ΔQ≤η為止；

步驟2將2個(gè)合并后的模塊視為新模塊，合并新模塊和其鄰居模塊并計(jì)算模塊度增益ΔQ，對(duì)網(wǎng)絡(luò)中所有新模塊都進(jìn)行鄰居模塊合并，直到模塊度增益ΔQ≤η為止；

步驟3重復(fù)步驟2，直至整個(gè)網(wǎng)絡(luò)的模塊度增益ΔQ不大于迭代終止系數(shù)η，此時(shí)所得出的模塊結(jié)構(gòu)即為帶時(shí)間約束的Louvain算法最終所求的模塊劃分。算法流程如圖4所示。

Figure 4 Flow chart of time-constrainted Louvain algorithm圖4 帶時(shí)間約束的Louvain算法流程圖

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

3.1 數(shù)據(jù)來源和預(yù)處理

本文實(shí)驗(yàn)的數(shù)據(jù)來源于公開腦網(wǎng)絡(luò)數(shù)據(jù)庫INDI(International Neuroimaging Data-Sharing Initiative)中的斯坦福大學(xué)(Stanford University)采集的40組腦網(wǎng)絡(luò)數(shù)據(jù)樣本組成的樣本集S[28]。40個(gè)樣本年齡在7.5～12.9歲，其中20個(gè)自閉癥患者(ASD組)，樣本年齡在7.5～12.9歲，20個(gè)健康人(TC組)，年齡在7.8～12.4歲。重復(fù)時(shí)間TR(Repetition Time)為2 s，回波時(shí)間TE(Time of Echo)為30 ms，采集矩陣為64*64，視野FOV(Field Of View)為20 cm*20 cm，翻轉(zhuǎn)角Flip Angle為80°，共采集180個(gè)時(shí)間點(diǎn)數(shù)據(jù)。

首先對(duì)由rs-fMRI采集得到的原始數(shù)據(jù)進(jìn)行預(yù)處理，具體過程包括去除不穩(wěn)定時(shí)間點(diǎn)；時(shí)間層校正；頭動(dòng)校正；顱骨去除；空間標(biāo)準(zhǔn)化；帶通濾波等，根據(jù)AAL模板將全腦共劃分成90個(gè)腦區(qū)，提取得到90個(gè)腦區(qū)的BOLD信號(hào)時(shí)間序列。數(shù)據(jù)預(yù)處理采用的平臺(tái)工具是duke brain imaging analysis center的Resting State Pipeline[29]。

3.2 動(dòng)態(tài)腦功能網(wǎng)絡(luò)的構(gòu)建

本文實(shí)驗(yàn)數(shù)據(jù)共采集了180個(gè)時(shí)間點(diǎn)，去除前4個(gè)不穩(wěn)定時(shí)間點(diǎn)，利用滑動(dòng)窗口技術(shù)對(duì)BOLD-fMRI信號(hào)進(jìn)行窗口觀測(cè)，取滑動(dòng)窗口長(zhǎng)度為W=40，滑動(dòng)窗口步長(zhǎng)L=1，得到包含137個(gè)時(shí)間點(diǎn)的腦網(wǎng)絡(luò)觀測(cè)數(shù)據(jù)。在每一個(gè)時(shí)間點(diǎn)上計(jì)算兩兩腦區(qū)之間的BOLD信號(hào)的皮爾遜相關(guān)系數(shù)，將腦區(qū)間的皮爾遜相關(guān)系數(shù)作為腦區(qū)之間的連接強(qiáng)度系數(shù)，可以得到每一個(gè)時(shí)間點(diǎn)上的90*90的腦網(wǎng)絡(luò)相關(guān)系數(shù)矩陣。其中，皮爾遜相關(guān)系數(shù)取值為-1～1，本文將腦網(wǎng)絡(luò)相關(guān)系數(shù)矩陣取絕對(duì)值表示其相關(guān)程度，絕對(duì)值越大相關(guān)程度越高，并去除自連接，將對(duì)角線元素置0。

此外，由于本文的實(shí)驗(yàn)組是ASD組和TC組各20個(gè)樣本，為了提取出ASD組和TC組各自的特征，對(duì)ASD組和TC組樣本分別進(jìn)行了組平均的統(tǒng)計(jì)。具體地，將ASD組20個(gè)樣本數(shù)據(jù)中相同時(shí)間點(diǎn)的腦網(wǎng)絡(luò)相關(guān)系數(shù)矩陣取平均得到組平均數(shù)據(jù)，TC組同上，最終將得到的ASD組和TC組 的組平均數(shù)據(jù)進(jìn)行閾值化處理[30]。根據(jù)已有研究[31]，腦網(wǎng)絡(luò)相關(guān)系數(shù)矩陣的閾值選取在0.45以上時(shí)腦網(wǎng)絡(luò)具有小世界性和完整性。本文選取閾值為0.5將腦網(wǎng)絡(luò)相關(guān)系數(shù)矩陣進(jìn)行閾值化，整個(gè)數(shù)據(jù)預(yù)處理流程如圖5所示。

Figure 5 Experimental principle and process圖5 實(shí)驗(yàn)原理及流程

3.3 帶時(shí)間約束的Louvain算法的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化實(shí)驗(yàn)及結(jié)果分析

按照2.2.2節(jié)中帶時(shí)間約束的Louvain算法對(duì)經(jīng)過數(shù)據(jù)預(yù)處理的動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分，并分別對(duì)不同實(shí)驗(yàn)的步驟和實(shí)驗(yàn)結(jié)果進(jìn)行分析。

3.3.1 模塊度Q值分布對(duì)比實(shí)驗(yàn)及結(jié)果分析

AD9833芯片可通過VOUT引腳提供各種輸出，可完成方波輸出、正弦波輸出和三角波的輸出。控制寄存器的 OPBITEN（D5）和 Mode（D1）bits決定著AD9833芯片輸出波形的類型。其控制邏輯如表2所示。

首先，對(duì)預(yù)處理數(shù)據(jù)分別使用原Louvain算法和帶時(shí)間約束的Louvain算法進(jìn)行模塊劃分，并計(jì)算每一個(gè)時(shí)間點(diǎn)的模塊度，2種算法的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化Q值如圖6所示。

Figure 6 Comparison of modularity (Q) between Louvain algorithm and time-constrainted Louvain algorithm圖6 原Louvain算法和帶時(shí)間約束的Louvain算法Q值對(duì)比圖

圖6a、圖6b、圖6c所示為3個(gè)樣本分別使用2種算法進(jìn)行模塊劃分后的結(jié)果，圖6中T1為使用原Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分得到的模塊度Q值，T2為使用帶時(shí)間約束的Louvain算法得到的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊度Q值。由圖6可以看出，原Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行劃分后得到的模塊度正常值分布上限均超過了0.7，并且有較多的異常數(shù)據(jù)點(diǎn)，而帶時(shí)間約束的Louvain算法其模塊度值分布位于0.35～0.56，其Q值均處在合理范圍之內(nèi)，表明帶時(shí)間約束的Louvain算法能夠?qū)?dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行更好的模塊劃分。2種算法的Q值具體統(tǒng)計(jì)如表1所示。

Table 1 Comparison of modularity (Q) between Louvainalgorithm and time-constrainted Louvain algorithm

3.3.2 模塊劃分對(duì)比實(shí)驗(yàn)及結(jié)果分析

統(tǒng)計(jì)發(fā)現(xiàn)，基于帶時(shí)間約束的Louvain算法得到的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊劃分將90個(gè)腦區(qū)分成了6或7個(gè)不等的模塊。本文對(duì)ASD組和TC組分別進(jìn)行500次實(shí)驗(yàn)，滑動(dòng)窗口選取W=40，計(jì)算137個(gè)時(shí)間點(diǎn)的模塊劃分，對(duì)每一次實(shí)驗(yàn)統(tǒng)計(jì)模塊劃分?jǐn)?shù)量和每一個(gè)模塊內(nèi)含腦區(qū)數(shù)量。結(jié)果顯示，模塊劃分為6個(gè)模塊的次數(shù)占總實(shí)驗(yàn)的4.2%，模塊劃分為7個(gè)模塊的次數(shù)占總實(shí)驗(yàn)的95.6%，模塊劃分為5個(gè)模塊的占總實(shí)驗(yàn)次數(shù)的0.2%；而使用原Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分，得到4～5個(gè)不等的模塊，其中，模塊劃分?jǐn)?shù)量為4的次數(shù)占總實(shí)驗(yàn)次數(shù)的42.6%，模塊劃分?jǐn)?shù)量為5的次數(shù)占總實(shí)驗(yàn)次數(shù)的57.4%，部分實(shí)驗(yàn)結(jié)果如表2、圖7和圖8所示。

Table 2 Partial time point modularity and moduledivision number of time constrainted

圖7所示為原Louvain算法部分時(shí)間點(diǎn)的模塊劃分頻次統(tǒng)計(jì)圖，圖7中，縱坐標(biāo)為腦區(qū)數(shù)量，橫坐標(biāo)為動(dòng)態(tài)腦功能網(wǎng)絡(luò)劃分模塊編號(hào)，圖中最小模塊包含13個(gè)腦區(qū)，最大模塊包含34個(gè)腦區(qū);圖8所示為帶時(shí)間約束的Louvain算法部分時(shí)間點(diǎn)模塊劃分頻次圖，圖8中，最小模塊包含4個(gè)腦區(qū)，最大模塊包含23個(gè)腦區(qū)。從圖8中可以看出，帶時(shí)間約束的Louvain算法沒有高于30個(gè)腦區(qū)的模塊，而且有只包含幾個(gè)腦區(qū)的小模塊，因此與原Louvain算法相比，本文算法可以識(shí)別出小規(guī)模的模塊，對(duì)腦網(wǎng)路的模塊劃分具有合理性和有效性。

Figure 7 Brain region quantity statistics of each module at partial time points by Louvain algorithm圖7 原Louvain算法部分時(shí)間點(diǎn)各模塊包含腦區(qū)數(shù)量統(tǒng)計(jì)圖

Figure 8 Brain region quatity statistics of each module at partial time points by time-constrained Louvain algorithm圖8 帶時(shí)間約束的Louvain算法部分時(shí)間點(diǎn)各模塊包含腦區(qū)數(shù)量統(tǒng)計(jì)圖

3.3.3 TC組和ASD組的模塊化對(duì)比實(shí)驗(yàn)及結(jié)果分析

帶時(shí)間約束的Louvain算法對(duì)預(yù)處理數(shù)據(jù)進(jìn)行模塊劃分，對(duì)于每一個(gè)樣本，得到一組時(shí)間序列上的Q值，將患有自閉癥(ASD)患者的Q值和對(duì)照組(TC)的Q值進(jìn)行對(duì)比，其具體分布如圖9所示。

Figure 9 Compassion of modularity(Q) between TC&ASD圖9 TC&ASD組Q值對(duì)比

圖9a中，箱型圖x軸為健康人(TC)組和自閉癥(ASD)患者組，y軸為Q值。由圖9a中可以看出，ASD組的Q值整體上小于TC組的，而且，ASD組有很多時(shí)間點(diǎn)的Q值異常低，由此可以看出，TC組相比較于ASD組，有明顯的模塊化屬性，而ASD組出現(xiàn)了Q值降低，這顯示出由于疾病對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)的影響，病人的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊屬性受到了破壞，因此ASD組的功能模塊性能有所降低。圖9b為ASD組和TC組的Q值分布圖，圖中TC組的Q值分布相較于ASD組更加集中，說明TC組的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化程度變化相較于ASD組更小，TC組的動(dòng)態(tài)腦功能網(wǎng)絡(luò)與ASD組相比具有更好的穩(wěn)定性。圖9c中，x軸為時(shí)間點(diǎn)，y軸為Q值，圖9c中可以看出，在合理的Q值范圍內(nèi)，ASD組的Q值時(shí)間序列普遍低于TC組。表3為上述2個(gè)樣本的Q值序列統(tǒng)計(jì)。

Table 3 Statistics of modularity(Q) between TC&ASD

另外，本文還使用spss軟件對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行了樣本分析，如表4所示，首先對(duì)2組Q值序列進(jìn)行萊文方差等同性檢驗(yàn)，得到P=0.199，2組數(shù)據(jù)為等方差，再對(duì)2組數(shù)據(jù)進(jìn)行獨(dú)立樣本t檢驗(yàn)，得到P<<0.05，說明ASD組和TC組的Q值序列具有統(tǒng)計(jì)學(xué)意義上的差別，具體來說，自閉癥(ASD)患者組的模塊度Q值遠(yuǎn)低于健康人(TC)組，而模塊度Q值大小表明了動(dòng)態(tài)腦功能網(wǎng)絡(luò)的模塊性高低，自閉癥患者會(huì)表現(xiàn)出語言交流障礙和智能障礙，其大腦功能相較于正常人有一定的缺失。本實(shí)驗(yàn)得出自閉癥患者動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊性低于健康人，與自閉癥對(duì)大腦的影響相符，因此可以使用動(dòng)態(tài)腦功能網(wǎng)絡(luò)的模塊化特征對(duì)自閉癥致病原理以及其病理特征進(jìn)行探究，從而為自閉癥的診斷提供一定的依據(jù)。其t檢驗(yàn)結(jié)果如表4所示。

3.3.4 原Louvain算法和改進(jìn)的Louvain算法效率對(duì)比

對(duì)同一組數(shù)據(jù)分別使用原Louvain算法和帶時(shí)間約束的Louvain算法進(jìn)行模塊劃分，統(tǒng)計(jì)其所需要時(shí)間。本次實(shí)驗(yàn)平臺(tái)為：處理器Intel i5-4200H 2.8 GHz，內(nèi)存4 GB，顯卡NIVDIA GTX860M。

對(duì)同一樣本分別采用2個(gè)不同的算法進(jìn)行模塊劃分，并計(jì)算每一次實(shí)驗(yàn)所需要的運(yùn)行時(shí)間。由于帶時(shí)間約束的Louvain算法設(shè)置了迭代終止系數(shù)η，其迭代次數(shù)會(huì)小于原Louvain算法，因此帶時(shí)間約束的Louvain算法運(yùn)行時(shí)間少于原Louvain算法的，計(jì)算的效率也要高于原Louvain算法的。對(duì)3個(gè)樣本分別使用2種算法進(jìn)行模塊劃分，原Louvain算法平均計(jì)算時(shí)間為0.131 s，帶時(shí)間約束的Louvain算法平均計(jì)算時(shí)間為0.066 s，帶時(shí)間約束的Louvain算法相較于原Louvain算法效率提高了49.6%。

綜上所有實(shí)驗(yàn)結(jié)果所述，對(duì)本文提出的算法和原Louvain算法在動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化過程中的實(shí)驗(yàn)效果分析如下：

(1)從算法的原理上看，Louvain算法是以模塊度增益最大為目標(biāo)來進(jìn)行模塊劃分，且原Louvain算法針對(duì)的是大規(guī)模網(wǎng)絡(luò)(數(shù)百萬節(jié)點(diǎn))，而動(dòng)態(tài)腦功能網(wǎng)絡(luò)從規(guī)模上說屬于小規(guī)模網(wǎng)絡(luò)(數(shù)百個(gè)節(jié)點(diǎn))，因此若用原Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行模塊劃分將導(dǎo)致嚴(yán)重的偏差，影響模塊化效果；本文提出了基于時(shí)間序列的迭代終止系數(shù)，通過提前終止模塊合并得到具有更好效果的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊劃分，得到更加合理的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊度Q值。

Table 4 Independent t test between TC&ASD

(2)從樣本的模塊劃分來看，原Louvain算法進(jìn)行模塊劃分得到的模塊數(shù)量少于本文算法的且其模塊數(shù)量不穩(wěn)定，具有很大的隨機(jī)性，而本文算法在多次模塊劃分實(shí)驗(yàn)中得到的模塊劃分?jǐn)?shù)量基本不變，具有較好的實(shí)驗(yàn)穩(wěn)定性。此外，原Louvain算法得到的模塊包含大量腦區(qū)，缺乏對(duì)小模塊的識(shí)別能力，而本文算法能夠準(zhǔn)確識(shí)別出小模塊。

(3)從對(duì)不同類的樣本識(shí)別上看，本文所提出的帶時(shí)間約束的Louvain算法得到的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊劃分在自閉癥患者和健康人2類樣本上能夠有特別明顯的區(qū)別，自閉癥患者由于疾病影響導(dǎo)致其動(dòng)態(tài)腦功能網(wǎng)絡(luò)的模塊性降低，二者的實(shí)驗(yàn)結(jié)果滿足統(tǒng)計(jì)意義上的區(qū)別。

(4)從計(jì)算的效率上看，相較于原Louvain算法，本文提出的帶時(shí)間約束的Louvain算法具有更高的計(jì)算效率，其平均運(yùn)行時(shí)間也小于原Louvain算法的，在針對(duì)大量樣本的計(jì)算上能夠節(jié)省大量的時(shí)間和計(jì)算資源。

4 結(jié)束語

由于原Louvain算法過分追求模塊度Q值最大化而忽略了對(duì)一些小規(guī)模模塊的辨識(shí)，因此在使用Louvain算法進(jìn)行動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊劃分時(shí)，不能夠準(zhǔn)確辨識(shí)動(dòng)態(tài)腦功能網(wǎng)絡(luò)。針對(duì)這一問題，本文提出了帶時(shí)間約束的Louvain算法，考慮到動(dòng)態(tài)腦功能網(wǎng)絡(luò)的構(gòu)建是與時(shí)間密切相關(guān)的，并且其模塊化屬性也會(huì)隨著時(shí)間的不同而發(fā)生變化，因此本文給出了一個(gè)基于時(shí)間約束的迭代終止系數(shù)。該系數(shù)綜合考量了所有時(shí)間點(diǎn)上Q值的平均分布狀況，并設(shè)定尋找最佳模塊劃分結(jié)果的過程中，當(dāng)滿足該系數(shù)要求時(shí)就終止迭代，從而在避免Q值過大的同時(shí)保證模塊劃分結(jié)果的合理性。

使用該算法對(duì)健康人和自閉癥患者的動(dòng)態(tài)腦功能網(wǎng)絡(luò)進(jìn)行了模塊化分析，分別進(jìn)行了3個(gè)對(duì)比實(shí)驗(yàn)并進(jìn)行了分析。使用原Louvain算法和本文算法進(jìn)行動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊劃分，原Louvain算法得到的模塊度Q值基本超出了正常范圍，而經(jīng)過改進(jìn)后的算法得到的模塊度Q值均處于合理范圍之內(nèi)。此外，通過多次實(shí)驗(yàn)結(jié)果可知，原算法得到的模塊劃分存在隨機(jī)性，實(shí)驗(yàn)結(jié)果差異很大，模塊數(shù)量少且其單一模塊包含腦區(qū)數(shù)量過多，無法識(shí)別出小規(guī)模的動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊；而本文算法得到的模塊劃分較原算法有更高的穩(wěn)定性和更小的模塊規(guī)模，能夠識(shí)別動(dòng)態(tài)腦功能網(wǎng)絡(luò)中的小規(guī)模模塊。最后，通過對(duì)自閉癥患者樣本和健康人樣本的對(duì)照實(shí)驗(yàn)，發(fā)現(xiàn)自閉癥患者(ASD)動(dòng)態(tài)腦功能網(wǎng)絡(luò)模塊化程度比健康人(TC)的動(dòng)態(tài)腦功能網(wǎng)絡(luò)有所下降，且二者具有統(tǒng)計(jì)學(xué)意義上的區(qū)別。

本文提出的帶時(shí)間約束的Louvain算法對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)的模塊化過程中存在的Q值高于合理范圍并且對(duì)動(dòng)態(tài)腦功能網(wǎng)絡(luò)小規(guī)模模塊辨識(shí)度不高的問題提出了一種解決方案，為自閉癥的動(dòng)態(tài)腦功能網(wǎng)絡(luò)屬性研究及自閉癥的輔助診斷提供了一種輔助方法，同時(shí)也對(duì)其它同類動(dòng)態(tài)網(wǎng)絡(luò)模塊屬性的研究具有一定的參考價(jià)值。