功能MRI時間分辨率對基于格蘭杰因果模型的腦效應連接的影響

張欣欣，楊青青，于春水,2，梁 猛*

(1.天津醫科大學醫學影像學院，天津 300203；2.天津醫科大學總醫院醫學影像科，天津 300052)

格蘭杰因果關系分析(Granger causality analysis，GCA)是揭示腦區間效應連接(effective connectivity，EC)的強有力工具，已廣泛用于功能MRI(functional MRI，fMRI)研究[1-5]。利用GCA推斷因果性的基本原理如下：若腦區X的時間序列在過去時刻的值對預測腦區Y的時間序列當前時刻的值有顯著貢獻，則認為腦區X的功能活動為因，Y則為果。因此，fMRI數據的時間分辨率可能是檢測腦區間EC的潛在制約因素[6-7]。既往研究[6，8]基于仿真數據提出較高時間分辨率可提升GCA結果的穩定性和敏感性。目前fMRI數據最常用的重復時間(repetition time，TR)多為2～3 s。多層同步采集序列[9]的出現，將全腦fMRI采集時間縮短至亞秒級，而這是否會對GCA結果產生顯著影響目前尚不得而知。本研究觀察fMRI時間分辨率對基于GCA模型的腦EC的影響。

1 資料與方法

1.1 數據來源 基于人類大腦連接組計劃(https://neuroscienceblueprint.nih.gov/human-connectome/connectome-programs)，于公開數據庫隨機抽取40名健康受試者，男20名，女20名，年齡30～37歲，平均(34.3±1.1)歲。納入標準：右利手，無神經或精神疾病史[10]。

1.2 數據預處理 采用梯度回波-回波平面成像(gradient echo-echo planar imaging，GRE-EPI)序列采集顱腦靜息態fMRI(resting-state fMRI，rs-fMRI)數據，TE 33.1 ms，TR 720 ms，FOV 208 mm×180 mm，矩陣104×90，平面分辨率2 mm×2 mm，FA 52°，層厚2 mm，層數72，帶寬2 290 Hz/Pixel，共采集1 200個時間點數據，采集時間14.4 min。采用Matlab R2013b和DPARSF軟件行圖像預處理，包括梯度失真校正、頭動校正、EPI圖像失真校正、空間配準、去除線性漂移、回歸協變量和樣條插值及帶通濾波(0.01～0.08 Hz)，見圖1。

圖1 受試者女，36歲 A～C.分別為預處理后腦部冠狀位、矢狀位及軸位fMRI數據圖

1.3 fMRI數據下采樣 為保證數據一致性，于上述預處理后，經時間下采樣生成采集不同時間分辨率的fMRI數據。最高時間分辨率(TR 0.72 s)數據為上述預處理后共1 200個時間點的fMRI數據；中等時間分辨率(TR 1.44 s)數據為上述數據經2倍下采樣獲得，取奇數時間點，共600個時間點；最低時間分辨率(TR 2.16 s)數據為最高分辨率數據經3倍下采樣獲得，共400個時間點。

1.4 GCA分析 基于人類腦網絡組圖譜(http://atlas.brainnetome.org)[11]將大腦分為246個腦區，以每個腦區所有體素的fMRI時間序列的均值作為該腦區的平均時間序列。以Matlab和REST軟件(http://restfmri.net/forum)對所有腦區平均時間序列行兩兩間GCA。采用基于殘差的GCA分析方法，假設X、Y分別為2個腦區的時間序列，以下列公式計算腦區間的GCA關系：

(1)

var(μt)=R2

(2)

(3)

(4)

var(μ′t)=T2

(5)

(6)

公式(1)為自回歸模型，代表利用Y腦區過去時刻的神經活動信息預測其現在時刻神經活動的模型，其中p為模型階數(即使用過去p個時間點的信息)，Mi為回歸系數，εt為殘差，殘差方差即為R1。公式(2)表示在公式(1)基礎上進一步增加X腦區過去時刻的神經活動信息，以共同預測Y腦區現在時刻的神經活動，模型殘差μt的方差記為R2。利用公式(3)可基于R1和R2計算X腦區對Y腦區的效應連接Fx→y。若X腦區過去時刻的神經活動信息對Y腦區現在時刻神經活動具有預測能力，即R20，表示X腦區是Y腦區的格蘭杰原因；反之，R2≈R1，Fx→y≈0，則X腦區不是Y腦區的格蘭杰原因。同理，根據公式(4)～(6)可計算Y腦區對X腦區的效應連接Fy→x。

1.5 統計學分析 采用Matlab R2013b軟件，以置換檢驗及整體錯誤率多重比較校正方法[12]，分別于個體水平及組水平確定每條EC的統計顯著性。確定個體水平統計顯著性方法：將所有受試者同一腦區的fMRI信號序列在時間順序上以相同方式進行隨機置換，確保針對同一腦區的置換方式在受試者間一致；針對置換后的fMRI數據重新計算每名受試者所有腦區間的EC；對每名受試者重復上述過程1 000次，每次取全腦最大EC值；1 000個全腦最大EC值構成該受試者的零分布，并將其真實數據所得到的每一條EC值分別與零分布進行比較(針對每條EC值，P值計算為1 000次隨機置換中所有大于等于該真實EC值的隨機置換所獲得EC值所占比例)，即得到每名受試者每條EC經校正后的統計顯著性。

確定組水平統計顯著性方法：計算每次隨機置換后每條連接所有受試者EC的平均值，并獲取全腦所有連接的最大平均EC值；1 000次隨機置換所得全腦最大EC構成組水平零分布，將其與真實數據所得的每條連接的組水平平均EC值進行比較，即得到每條EC經校正后的組水平統計顯著性。校正后P<0.05為有顯著統計學意義。

2 結果

組水平不同時間分辨率fMRI數據的GCA結果見圖2。模型階數為1、2、3對應的全腦顯著EC條數在TR為2.16 s時分別為792、128及1 504，TR為1.44 s時分別為2 932、1 920及12 500，TR為0.72 s時分別為16 320、43 177及60 204。個體水平上所獲結果與上述組水平結果趨勢基本一致，見圖3。

圖2 基于不同時間分辨率及模型階數fMRI數據的組水平GCA結果圖

圖3 基于不同時間分辨率及模型階數的fMRI數據的個體水平GCA結果圖 A.受試者男，36歲；B.受試者女，34歲；C.受試者女，33歲；D.受試者女，36歲

3 討論

基于GCA的EC取決于對過去時刻的定義，而模型中過去時刻所包含的具體數據則由單位時間區間內所含數據及過去時刻區間長度共同決定；單位時間區間內所含的數據取決于時間分辨率，而過去時刻區間長度取決于TR(即TR越短，時間分辨率越高)與模型階數的乘積。

相同時間分辨率條件下，隨模型階數提高，單位時間區間內所含數據不變，但所含過去時刻區間長度增加；雖然檢測到的顯著EC條數也呈現增多趨勢，但不同情況下結果并不完全一致，提示模型階數增多、即過去時刻區間長度增加并不一定使EC條數顯著增加。

當過去時刻區間長度相同時，即TR為0.72 s且p=2、TR為1.44 s且p=1條件下，GCA均取決于過去1.44 s的區間；TR為0.72 s且p=3、TR為2.16 s且p=1條件下，GCA均取決于過去2.16 s的區間；TR為1.44 s且p=3、TR為2.16 s且p=2條件下，GCA均取決于過去4.32 s的區間；亦即較高時間分辨率數據檢測出的顯著EC條數總是遠高于較低分辨率數據，提示當所利用的過去時刻區間相同時，時間分辨率越高，單位時間區間內數據所含信息越豐富，檢測到的EC條數越多。

此外，相同模型階數條件下，時間分辨率越高，單位時間區間內所含數據越豐富，過去時刻區間長度越短。因此，雖然過去時刻區間長度縮短并不一定總是使檢測到的EC條數顯著減少，但總體呈減少趨勢。在此趨勢影響下，隨時間分辨率增高，檢測到的顯著EC條數仍然明顯增多。綜合以上結果，不論模型階數如何，隨著時間分辨率提高，檢測到的顯著EC條數均明顯增多，提示將fMRI數據的時間分辨率提高至亞秒級可提高檢測EC的敏感度。

WITT等[13]采用仿真fMRI時間序列數據比較結構方程建模、GCA和動態因果建模3種模型受TR的影響程度，發現結構方程建模對TR變化最不敏感而GCA最為敏感，提示從原始fMRI數據中能否正確推斷出GCA因果關系在極大程度依賴于時間分辨率。LIN等[14-15]研究表明，雖然測量血氧水平依賴信號可反映神經元的血流動力學響應，但無法直接代表神經元發生的生理事件；這是由于病理條件可能導致某些大腦區域血流動力學響應發生變異，而提高采樣率可顯著提高腦區血流動力學響應與神經元活動的相關性，故高時間分辨率fMRI數據能更真實地反映神經活動[16]。數據采樣率較低會使GCA結果在較大程度上受到血流動力學響應變化的影響，并偏離真實的神經元活動，即低時間分辨率數據的GCA結果在不同血流動力學響應下并不穩定，而高時間分辨率數據的GCA結果在血流動力學響應發生變化時仍保持穩定[7]，表明基于較高時間分辨率fMRI數據的GCA結果能更真實地反映腦區間神經元活動的因果關系。行GCA分析之前可通過采集同步電生理數據估計血流動力學響應情況，再針對去卷積后數據進行分析，但同步電生理-MR技術對受試者配合度要求較高，臨床難以普遍開展；而通過提高時間分辨率以實現真實反映腦區間的因果關系更具現實及指導意義。

綜上，提高fMRI數據時間分辨率可對GCA結果產生顯著影響，高時間分辨率可顯著提高GCA的敏感度，更加真實地反映不同腦區神經活動之間的因果關系。為得到更為穩定可靠的結果，欲行GCA分析時，建議以高時間分辨率采集fMRI數據。