肖 丹， 楊智鵬， 吳 錫， 周激流

(1.成都信息工程大學計算機學院， 成都 610225； 2.成都信息工程大學電子工程學院， 成都 610225；3.四川大學計算機學院， 成都 610065)

1 引 言

彌散加權成像(Diffusion Weighted MRI，DWI)是一種能夠在大腦白質內檢測出水分子彌散運動的無創(chuàng)方法.通過估計體素中水分子的彌散方向分布函數(shù)(diffusion Orientation Distribution Functions， dODFs)來間接計算白質纖維的分布方向[1].DWI纖維束追蹤成像就是將dODFs轉換成纖維方向分布函數(shù)(fiber Orientation Distribution Functions，fODFs)，并通過其在體素間的連通性來構建腦白質連接的解剖結構.研究腦白質神經纖維的組織特性,能夠了解到大腦的發(fā)育、衰老和患病情況[2].

現(xiàn)有的DWI纖維重建算法可分為局部纖維重建方法和全局纖維重建方法.局部纖維重建方法是從初始點開始，沿著纖維走向逐步前進，最終獲得整條纖維路徑；全局纖維重建方法則是在互相連接的纖維路徑上建立代價函數(shù)，利用優(yōu)化技術尋找最佳纖維路徑.全局纖維重建方法可以消除累積噪聲及局部隨機噪聲，提高長距離成像的可靠性.

全局纖維重建方法中，目前常用的是在全局概率追蹤的貝葉斯算法中加入先驗信息，從而在兩個區(qū)域之間找到最優(yōu)纖維束的方法[3]，用類似大腦幾何形狀的模擬DWI作為可靠性估計進行定量評估[4].這種方法的缺陷在于可供使用的先驗知識只包含了兩區(qū)域間是否存在連接的信息，并不包括關于纖維束位置或功能的先驗信息.此外，由于問題過于復雜，很難求出最優(yōu)解，此方法只能通過從后驗分布的啟發(fā)式采樣來估量纖維束.2012年，Wu提出了一種將全局纖維追蹤與分層纖維聚類相結合的劃分纖維路徑的方法.這種方法采用了K均值聚類和改進的休伯特統(tǒng)計，在每個纖維束上進行迭代采樣和聚類來逼近最優(yōu)解，極大地促進了纖維束成像在人類復雜神經網絡的臨床研究[5].Eleftherios等人在進行白質纖維虛擬解剖研究時，采用了基于纖維流線的纖維束標記及聚類技術，將纖維束模型作為先驗知識，檢測出纖維束線圖中更為精準的相似流線和束[6].由此可見，評價方程的計算作為全局優(yōu)化類方法的核心，決定了優(yōu)化后的結構連接是否符合真實生理結構.

重建具有功能意義的結構連接是神經科學研究中基礎性的問題，以往的研究將DWI的結構連接與基于灰質中功能磁共振成像相結合，重建出連通多個灰質功能區(qū)域的白質結構連接.但目前的這些融合研究只是一種基本的聯(lián)合，結果只能說明在特定的灰質功能區(qū)有白質纖維連接，白質結構本身并未證明具有功能特性.最新的研究表明，白質中的功能磁共振成像(fMRI)能通過測量白質神經元功能活動中的血氧依賴水平(Blood Oxygen Level Dependent, BOLD)來分析神經纖維的功能特性，并已成功應用于病理學研究，該研究為重建具有功能特性的纖維束提供了可能[7-8].

綜上所述，本研究提出將白質fMRI融合到DWI全局優(yōu)化纖維重建中，加入功能先驗信息纖維重建出最優(yōu)功能路徑.該方法基于全局優(yōu)化類的貝葉斯最優(yōu)路徑算法，從全局纖維中找到連接特定功能區(qū)域的最優(yōu)路徑，對數(shù)據(jù)進行初始化.其可有效抑制局部噪聲，得到執(zhí)行特定功能的最優(yōu)連接路徑，避免得到局部最優(yōu)解.較之現(xiàn)有方法，打破了僅通過空間位置形成最優(yōu)路徑的框架，重建出在執(zhí)行特定腦活動時，大腦信息傳遞的最優(yōu)路徑.

2 方 法

2.1 DWI優(yōu)化方法

大腦的DWI數(shù)據(jù)定義為連接圖G=(V,E,wE)，其中，V是除去腦脊液(CSF)以外的所有體素節(jié)點集，E是邊集，wE是邊的權重.在三維圖像中每個節(jié)點都能被邊e∈E連接到其3×3×3鄰域中，并給每條邊e賦予一個權重wE(e)∈[0,1]，用于表示纖維束連接其兩個端節(jié)點的概率.路徑的似然值是路徑上所有的邊權重wE(e)的乘積，即：

(1)

式(1)中，v∈V和v′∈V是G中的兩個節(jié)點，πv,v′是連接這兩點的路徑，可表示成節(jié)點序列πv,v′=[v1,v2,...,vn]，其中，v1=v,vn=v′,(vi,vi+1)∈E,i=1,...,n-1.路徑的基數(shù)等于其節(jié)點總數(shù)|πv,v′|=n.

我們用單位球面S2上的任意方向θ的fODFf:S2→R+求出纖維在該方向的概率，從而表示DWI的彌散情況.對每個體素的26個相鄰體素方向θi,i=1,...,26進行研究.通過計算在所有方向集Ci的fODF，得到體素在方向θi∈S2上的權重w(θi)[9].權重w(θi)表示體素連接該方向的概率，可近似表示為

(2)

wE(v,v′)=1/2·(w(v→v′)+w(v′→v))

(3)

其中，v→v′表示從體素v到體素v′的方向，于是得到對稱邊權重：wE(v,v′)=wE(v′,v).

2.2 大腦白質fMRI信號建模

用DWI中fMRI相關張量來重建人腦中的功能結構，通過BOLD信號的時間波動來反映自發(fā)神經活動以及功能刺激下的誘發(fā)反應.對于BOLD數(shù)據(jù)集中的每個體素，可以構造時空相關張量以表征體素與其鄰域之間的時間相關性的局部分布[10].

假設F是要構建的空間相關張量，估計的相關系數(shù)D沿單位向量ni(xi,yi,zi)投影得到

(4)

其中，t表示轉置操作.

D=(D1,D2,...,D26)t表示沿著26個方向觀察到的時間相關性的集合，F(xiàn)D是F重新排列后形成的列向量，則D和FD之間的關系可以表示為

D=M·FD

(5)

FD=(Mt·M)-1·Mt·D

(6)

其中，-1表示逆矩陣.

相關張量F(對應于最大特征值的特征向量)的主特征向量表示時間相關性的主要方向.本文假定該方向是局部小鄰域窗口內的神經活動傳播的方向.pF是功能ODF，它由吉布斯分布建模計算得到[11].該模型假定體素X中張量F僅取決于局部主方向VF(X).pF則可用以下公式表示.

(7)

其中，ZF是標準化常數(shù).

(8)

方程中的勢函數(shù)p隨著函數(shù)方向VF(X)和最大張量特征值λ1之間的差異而減小.分母用張量范數(shù)進行歸一化[12].對于各向異性張量，勢能給出的概率分布集中在張量F的主特征向量的方向上.對于各向同性張量，勢能函數(shù)將形成更寬的概率分布.

2.3 融合fMRI的DWI纖維優(yōu)化重建

對于DWI圖像中每個節(jié)點v∈V，pF(v)∈[0,1]表示該節(jié)點位于路徑中的功能先驗概率，使其在執(zhí)行特定腦活動時，能根據(jù)功能信息形成大腦信息傳遞的最優(yōu)路徑.本實驗提出一種有效的算法，即將腦圖G=(V,E,wE)中沿邊緣連接的貝葉斯模型，與表示節(jié)點功能信息的pF:V→+相結合. 邊緣連接的貝葉斯模型可通過之前的節(jié)點和邊e∈E的轉化來構建.對于單邊e=(v,v′)∈E，路徑的功能先驗概率P(e)定義為纖維束在v點處的功能概率pF(v)與v′點處的功能概率pF(v′)乘積的平方根.

(9)

給圖像中每條邊分配一個邊緣權重wE(e)，用于表示大腦沿邊e的結構連通性. 將公式(3)中的邊緣概率wE用概率密度函數(shù)fe來表征.

(10)

P(e|wE(e))∝P(wE(e)|e)P(e)=wE(e)P(e)

(11)

對于任意邊e(v,v′)，有

-log(wE(e)P(e))=

(12)

(13)

(14)

argmaxπv,v′P(πv,v′|G)

(15)

在以往的方法中，上式中路徑概率是體素在這條路徑上先驗信息的乘積[13].這種方法的局限性在于，只能使用特定且有限的解剖先驗信息(區(qū)域標簽).

而本文提出的方法對邊緣先驗進行了重新定義，得到圖像中可能存在的所有路徑，通過修改圖中的先驗信息來直接求出纖維束成像的最優(yōu)路徑，并不需要從后驗分布來進行路徑采樣.這種方法提供了大量的最優(yōu)路徑求解衍生算法，大大簡化了計算量.

WM中的ODF的例子如圖1所示.其中,P(πv,v′|G)是后驗ODF(c)，從DWI計算的ODF(b)和相關張量的ODF(a)計算得到，用于表示體素內的功能通路方向.

(a)

(b)

(c)

圖1 單個像素功能路徑方向的后驗ODF

(a)功能ODF；(b)彌散ODF；(c)后驗ODF

Fig.1 Posterior ODF of functional pathway directions for a voxel

(a) Function ODF; (b) Diffusion ODF; (c) Posterior ODF

3 實驗結果

3.1 數(shù)據(jù)采集

全腦MRI數(shù)據(jù)來自健康的成年志愿者.實驗儀器使用3T Philips Achieva scanner (Philips Healthcare, Inc., Best, Netherlands)，32通路頭部線圈.實驗數(shù)據(jù)集為四位成年人在進行感覺刺激實驗時的觸覺刺激功能圖像.感覺刺激被設計為方波形式，刷子刺激手掌30 s然后無刺激30 s，周期重復.

采集參數(shù)：T2*-weighted (T2* w) gradient echo (GE), echo planar imaging (EPI) 序列采集了三組BOLD信號：TR=3 s、TE=45 ms、matrix size=80×80、FOV=240×240 mm2、34層和3 mm 層厚、145 volumes、435 s.同時利用a single-shot, spin echo EPI序列采集了diffusion weighted images (DWI) 數(shù)據(jù)：b=1000 s/mm2、32 diffusion-sensitizing directions、TR=8.5 s、 TE=65 ms、SENSE factor=3、matrix size=128×128、FOV=256×256、68 層和2 mm層厚.為提供解剖學依據(jù)，所有例均采集3D高分辨T1-weighted (T1w) 解剖結構圖像，利用multi-shot 3D GE序列采集，像素大小1×1×1 mm3.

3.2 數(shù)據(jù)預處理

采集的數(shù)據(jù)均使用SPM12工具箱進行預處理.BOLD信號依次經過時間層矯正、頭動矯正、FWHM=4 mm高斯平滑.如果頭動位移超過2 mm以下、旋轉大于2°，數(shù)據(jù)將被剔除.以b=0的DWI數(shù)據(jù)為參考， 將所有被試者的平滑后的數(shù)據(jù)配準到各自的DWI數(shù)據(jù)空間.T1w數(shù)據(jù)進行偏移矯正和分割得到白質、灰質和腦脊液，然后共同以b=0 DWI為參考配準到DWI圖像空間.

3.3 優(yōu)化重建實驗

本研究分析了人腦臨床數(shù)據(jù)的兩個白質區(qū)域：丘腦至機體感覺區(qū)域和丘腦至島葉區(qū)域.

圖2展示了從丘腦至機體感覺區(qū)域的跟蹤結果，圖中紅色標記部分為大腦中央后回區(qū)域，從生理學分析，大腦中央后回接受背側丘腦腹后核傳來的對側軀干四肢的痛、溫、觸壓覺及位置和運動覺，在刺激實驗者手掌時處于激活狀態(tài).

如圖2所示，采用傳統(tǒng)DWI算法得到了圖中黃色所示的大面積皮層區(qū)域的通路，而融合白質功能信號的優(yōu)化算法則能直接找到功能激活區(qū)域的通路.說明本實驗優(yōu)化算法能夠重建執(zhí)行特定腦活動時功能激活的白質纖維通路.如圖3所示，優(yōu)化算法相較于傳統(tǒng)DWI算法，其概率密度更集中.

本文定義真陽性值(TP)來反映高概率區(qū)域中包含多少體素，從而對兩種方法進行定量比較：

(16)

表1展示了傳統(tǒng)DWI方法和本文優(yōu)化算法的真陽性平均值和標準差，由表可知，優(yōu)化算法相較于傳統(tǒng)DWI算法，其真陽性參數(shù)平均值更大，方差更小.由此可知，本實驗優(yōu)化算法重建功能激活狀態(tài)的通路時，獲得纖維束更集中緊湊，較之現(xiàn)有方法具有更強的魯棒性.

表1 丘腦至機體感覺區(qū)域真陽性平均值和標準差

Tab.1 Mean and standard deviation of true positive scores from thalamus to sensory area

算法真陽性平均值和標準差傳統(tǒng)DWI0.073 2±0.026 優(yōu)化算法0.283 5±0.011

(a) (b)

(c) (d)

圖2 丘腦至機體感覺區(qū)域的跟蹤結果

(a)～(d)分別為4個例子，其中每一例的第一排是冠狀面視角，第二排為矢狀面視角；灰色區(qū)域為種子點區(qū)域，黃色區(qū)域為傳統(tǒng)DWI得到的大腦皮層區(qū)域，紅色區(qū)域為觸覺刺激激活的皮層區(qū)域(為方便觀察流線，隨機選擇顏色進行展示)

Fig.2 Results of tracking from the thalamus to sensory area

There are 4 examples which the seed area are in gray color. The target ROI of traditional DWI method is in yellow color and activation area is in red. For each subject, the upper rows a coronal view and lower row is an sagittal view. Note that the streamlines are randomly colored for visual effect.

圖3 丘腦至機體感覺區(qū)域的概率密度圖

(a)～(d)分別為4個例子，其中每例第一排是為冠狀面視角，第二排為矢狀面視角；其中，黃色為高密度區(qū)域.

Fig.3 Probability density map from the thalamus to sensory area

There are 4 examples which the mean path directions with high probability in yellow color.

本實驗還重建了被實驗者受到手掌刺激時丘腦到島葉的流線.文獻[14-15]發(fā)現(xiàn)島葉前部與丘腦有神經相連，并且該路徑與觸覺表達相關.由圖4可知，在重建丘腦與島葉區(qū)域的通路時，由于該區(qū)域白質纖維流向復雜，傳統(tǒng)DWI算法用更多的追蹤次數(shù)重建出來的纖維束更為分散，可靠性降低.而融合白質功能信號的優(yōu)化算法用較少的纖維束追蹤次數(shù)直接重建出島葉前半部分的通路.

由圖5和表2也可知，優(yōu)化算法的實驗結果更集中緊湊.

表2 丘腦至島葉區(qū)域真陽性平均值和標準差

Tab.2 Mean and standard deviation of true positive scores from thalamus to insula

算法真陽性平均值和標準差傳統(tǒng)DWI0.105 4±0.057 優(yōu)化算法0.251 1±0.014

(a) (b)

(c) (d)

圖4 丘腦至島葉區(qū)域的跟蹤結果

(a)～(d)分別為4個例子的剖面視角；灰色區(qū)域為種子點區(qū)域，紅色區(qū)域為目標ROI區(qū)域(為方便觀察流線，隨機選擇顏色進行展示)

Fig.4 Results of tracking from the thalamus to insula

The seed area is in red color and the target ROI is in gray color. There are axial views. Note that the streamlines are randomly colored for visual effect.

(a) (b)

(c) (d)

圖5 丘腦至島葉區(qū)域的概率密度圖

(a)～(d)分別為4個例子的剖面視角；其中，黃色為高密度區(qū)域

Fig.5 Probability density map from the thalamus to insula

There are 4 examples which the mean path directions with high probability in yellow color.

4 結 論

本文闡述了彌散磁共振纖維跟蹤重建可以融合白質fMRI信號的功能信息，用于跟蹤特定神經活動狀態(tài)下的功能活動通路. 時空相關張量對白質纖維沿線的功能特性建模，為基于DWI跟蹤流程的跟蹤方向提供有效的補充信息.此優(yōu)化重建方法是研究人腦的結構功能相互關系的有用補充.

首先，本研究分析了連接丘腦和機體感覺區(qū)域的纖維束通路.中央后回接受背側丘腦腹后核傳來的對側軀干四肢的痛、溫、觸壓覺及位置和運動覺，在刺激實驗者手掌時處于激活狀態(tài).傳統(tǒng)DWI算法無法精準地找到功能激活狀態(tài)的通路，而本研究優(yōu)化方法能精準有效地找到了白質纖維束的功能活性通路，說明了本研究有助于解決白質功能結構跟蹤困難的問題.

接著，本研究還分析了連接丘腦和島葉皮層的纖維束通路.島葉皮層與多種腦功能相關，如感知、運動控制、自我意識，認知功能和個體經驗.腦島的前部接收來自丘腦腹側內側核基部的直接信息.這條路徑穿過復雜區(qū)域與其他功能路徑交叉，基于DWI的跟蹤重建在經過大量的跟蹤次數(shù)后才能找到其通路，實驗結果繁多且分散.相反，本文優(yōu)化方法可以聯(lián)合丘腦和腦島之間的功能信息重建這條路徑，重建結果集中且緊湊.實驗結果表明本文優(yōu)化方法具有研究島葉皮層功能結構以及與其他區(qū)域連接的潛力.

綜上所述，本研究提出了一種新的白質纖維束優(yōu)化重建方法，它基于全局優(yōu)化類的貝葉斯最優(yōu)路徑算法，利用白質中的功能信息來找到在執(zhí)行特定腦活動時，大腦信息傳遞的最優(yōu)路徑.白質中fMRI信號的各向異性被建模為時空相關張量，并調制用于跟蹤的彌散信號導出的ODF.本文融合白質功能信號的DWI纖維優(yōu)化重建具有在特定功能回路中重建纖維通路的巨大潛力.