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

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

1 引言

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

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

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

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

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

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

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

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

模塊性是復(fù)雜網(wǎng)絡(luò)系統(tǒng)的基本特征，網(wǎng)絡(luò)模塊是指一組連接緊密的節(jié)點，通常每一個模塊是專門處理信息的基礎(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時，網(wǎng)絡(luò)模塊的劃分最好[26,27]。

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

圖2所示的是使用Louvain算法對動態(tài)腦功能網(wǎng)絡(luò)進行模塊劃分得出的不同時間點的模塊度Q值時間序列，其中，圖2a和圖2c為2個健康人樣本，圖2b為自閉癥患者樣本。由圖2可以看出，3個樣本的峰值分別達到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ī)模過大的腦功能模塊可能會使得一些重要的亞功能腦模塊被忽略，不利于深入研究腦功能網(wǎng)絡(luò)的屬性。

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

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

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

2.2 帶時間約束的Louvain算法原理與步驟

2.2.1 帶時間約束的Louvain算法

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

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

(1)

(2)

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

定義2改進的模塊度增益ΔQ：

(3)

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

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

(4)

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

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

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

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

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

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

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

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

本文實驗的數(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個樣本年齡在7.5～12.9歲，其中20個自閉癥患者(ASD組)，樣本年齡在7.5～12.9歲，20個健康人(TC組)，年齡在7.8～12.4歲。重復(fù)時間TR(Repetition Time)為2 s，回波時間TE(Time of Echo)為30 ms，采集矩陣為64*64，視野FOV(Field Of View)為20 cm*20 cm，翻轉(zhuǎn)角Flip Angle為80°，共采集180個時間點數(shù)據(jù)。

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

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

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

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

Figure 5 Experimental principle and process圖5 實驗原理及流程

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

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

3.3.1 模塊度Q值分布對比實驗及結(jié)果分析

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

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

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

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

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

3.3.2 模塊劃分對比實驗及結(jié)果分析

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

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

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

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

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

3.3.3 TC組和ASD組的模塊化對比實驗及結(jié)果分析

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

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

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

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

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

3.3.4 原Louvain算法和改進的Louvain算法效率對比

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

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

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

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

Table 4 Independent t test between TC&ASD

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

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

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

4 結(jié)束語

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

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

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