融合白質(zhì)功能信號(hào)的DWI纖維優(yōu)化重建

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

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

1 引 言

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

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

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

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

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

2 方 法

2.1 DWI優(yōu)化方法

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

(1)

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

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

(2)

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

(3)

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

2.2 大腦白質(zhì)fMRI信號(hào)建模

用DWI中fMRI相關(guān)張量來(lái)重建人腦中的功能結(jié)構(gòu)，通過(guò)BOLD信號(hào)的時(shí)間波動(dòng)來(lái)反映自發(fā)神經(jīng)活動(dòng)以及功能刺激下的誘發(fā)反應(yīng).對(duì)于BOLD數(shù)據(jù)集中的每個(gè)體素，可以構(gòu)造時(shí)空相關(guān)張量以表征體素與其鄰域之間的時(shí)間相關(guān)性的局部分布[10].

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

(4)

其中，t表示轉(zhuǎn)置操作.

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

D=M·FD

(5)

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

(6)

其中，-1表示逆矩陣.

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

(7)

其中，ZF是標(biāo)準(zhǔn)化常數(shù).

(8)

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

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

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

(9)

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

(10)

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

(11)

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

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

(12)

(13)

(14)

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

(15)

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

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

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

(a)

(b)

(c)

圖1 單個(gè)像素功能路徑方向的后驗(yàn)ODF

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

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

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

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

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

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

采集參數(shù)：T2*-weighted (T2* w) gradient echo (GE), echo planar imaging (EPI) 序列采集了三組BOLD信號(hào)：TR=3 s、TE=45 ms、matrix size=80×80、FOV=240×240 mm2、34層和3 mm 層厚、145 volumes、435 s.同時(shí)利用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層厚.為提供解剖學(xué)依據(jù)，所有例均采集3D高分辨T1-weighted (T1w) 解剖結(jié)構(gòu)圖像，利用multi-shot 3D GE序列采集，像素大小1×1×1 mm3.

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

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

3.3 優(yōu)化重建實(shí)驗(yàn)

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

圖2展示了從丘腦至機(jī)體感覺(jué)區(qū)域的跟蹤結(jié)果，圖中紅色標(biāo)記部分為大腦中央后回區(qū)域，從生理學(xué)分析，大腦中央后回接受背側(cè)丘腦腹后核傳來(lái)的對(duì)側(cè)軀干四肢的痛、溫、觸壓覺(jué)及位置和運(yùn)動(dòng)覺(jué)，在刺激實(shí)驗(yàn)者手掌時(shí)處于激活狀態(tài).

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

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

(16)

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

表1 丘腦至機(jī)體感覺(jué)區(qū)域真陽(yáng)性平均值和標(biāo)準(zhǔn)差

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

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

(a) (b)

(c) (d)

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

(a)～(d)分別為4個(gè)例子，其中每一例的第一排是冠狀面視角，第二排為矢狀面視角；灰色區(qū)域?yàn)榉N子點(diǎn)區(qū)域，黃色區(qū)域?yàn)閭鹘y(tǒng)DWI得到的大腦皮層區(qū)域，紅色區(qū)域?yàn)橛|覺(jué)刺激激活的皮層區(qū)域(為方便觀察流線，隨機(jī)選擇顏色進(jìn)行展示)

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 丘腦至機(jī)體感覺(jué)區(qū)域的概率密度圖

(a)～(d)分別為4個(gè)例子，其中每例第一排是為冠狀面視角，第二排為矢狀面視角；其中，黃色為高密度區(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.

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

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

表2 丘腦至島葉區(qū)域真陽(yáng)性平均值和標(biāo)準(zhǔn)差

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

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

(a) (b)

(c) (d)

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

(a)～(d)分別為4個(gè)例子的剖面視角；灰色區(qū)域?yàn)榉N子點(diǎn)區(qū)域，紅色區(qū)域?yàn)槟繕?biāo)ROI區(qū)域(為方便觀察流線，隨機(jī)選擇顏色進(jìn)行展示)

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個(gè)例子的剖面視角；其中，黃色為高密度區(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 結(jié) 論

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

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

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

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