帶時間約束的Louvain算法在動態腦功能網絡模塊化中的應用研究*

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

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

1 引言

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

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

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

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

Figure 1 Processing of dynamic brain function network construction and module division圖1 動態腦功能網絡構建及模塊劃分流程

2 帶時間約束的Louvain算法動態腦功能網絡模塊化分析

2.1 動態腦功能網絡構建及模塊化流程

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

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

Louvain算法[21]是一種基于多層次優化模塊度的算法，它能夠發現層次性的模塊，并且模塊辨識準確性高、計算復雜度低、效率高，因此被廣泛應用于模塊化研究。

圖2所示的是使用Louvain算法對動態腦功能網絡進行模塊劃分得出的不同時間點的模塊度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值的合理范圍應該在0.3～0.7。規模過大的腦功能模塊可能會使得一些重要的亞功能腦模塊被忽略，不利于深入研究腦功能網絡的屬性。

Figure 2 Dynamic brain function network modularity (Q) sequence by Louvain algorithm圖2 動態腦功能網絡Louvain算法模塊化Q值序列

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

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

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

2.2.1 帶時間約束的Louvain算法

首先，針對動態腦功能網絡給出如下定義：

定義1根據Louvain算法，定義動態腦功能網絡中模塊化劃分優劣的判定標準，即模塊度Q為：

(1)

(2)

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

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

(3)

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

定義3改進的迭代終止系數η：

(4)

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

2.2.2 動態腦功能網絡模塊化實現步驟

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

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

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

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

3 實驗與結果分析

3.1 數據來源和預處理

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

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

3.2 動態腦功能網絡的構建

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

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

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

3.3 帶時間約束的Louvain算法的動態腦功能網絡模塊化實驗及結果分析

按照2.2.2節中帶時間約束的Louvain算法對經過數據預處理的動態腦功能網絡進行模塊劃分，并分別對不同實驗的步驟和實驗結果進行分析。

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

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

首先，對預處理數據分別使用原Louvain算法和帶時間約束的Louvain算法進行模塊劃分，并計算每一個時間點的模塊度，2種算法的動態腦功能網絡模塊化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種算法進行模塊劃分后的結果，圖6中T1為使用原Louvain算法對動態腦功能網絡進行模塊劃分得到的模塊度Q值，T2為使用帶時間約束的Louvain算法得到的動態腦功能網絡模塊度Q值。由圖6可以看出，原Louvain算法對動態腦功能網絡進行劃分后得到的模塊度正常值分布上限均超過了0.7，并且有較多的異常數據點，而帶時間約束的Louvain算法其模塊度值分布位于0.35～0.56，其Q值均處在合理范圍之內，表明帶時間約束的Louvain算法能夠對動態腦功能網絡進行更好的模塊劃分。2種算法的Q值具體統計如表1所示。

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

3.3.2 模塊劃分對比實驗及結果分析

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

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

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

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

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

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

帶時間約束的Louvain算法對預處理數據進行模塊劃分，對于每一個樣本，得到一組時間序列上的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組出現了Q值降低，這顯示出由于疾病對動態腦功能網絡的影響，病人的動態腦功能網絡模塊屬性受到了破壞，因此ASD組的功能模塊性能有所降低。圖9b為ASD組和TC組的Q值分布圖，圖中TC組的Q值分布相較于ASD組更加集中，說明TC組的動態腦功能網絡模塊化程度變化相較于ASD組更小，TC組的動態腦功能網絡與ASD組相比具有更好的穩定性。圖9c中，x軸為時間點，y軸為Q值，圖9c中可以看出，在合理的Q值范圍內，ASD組的Q值時間序列普遍低于TC組。表3為上述2個樣本的Q值序列統計。

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

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

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

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

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

綜上所有實驗結果所述，對本文提出的算法和原Louvain算法在動態腦功能網絡模塊化過程中的實驗效果分析如下：

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

Table 4 Independent t test between TC&ASD

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

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

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

4 結束語

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

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

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