基于擴(kuò)展的低階多元廣義線性模型的腦節(jié)點識別方法

楊雅倩，唐紹婷

(1.數(shù)學(xué)、信息與行為教育部重點實驗室(北京航空航天大學(xué))，北京 100191；2.北京航空航天大學(xué) 大數(shù)據(jù)科學(xué)與腦機(jī)智能高精尖創(chuàng)新中心，北京 100191；3.北京航空航天大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，北京 100191)(*通信作者電子郵箱tangshaoting@buaa.edu.cn)

0 引言

功能性磁共振成像(functional Magnetic Resonance Imaging, fMRI)技術(shù)是一種通過檢測血管中血氧水平變化來測量大腦活動的技術(shù)，它因其無創(chuàng)性和高空間分辨率而被廣泛運(yùn)用于大腦活動的研究當(dāng)中，現(xiàn)已有多篇相關(guān)文章，其中一個重要研究方向是比較不同條件下的大腦反應(yīng)。具體而言，在許多心理fMRI實驗中，研究人員對兩次實驗中的大腦反應(yīng)進(jìn)行比較，旨在找出不同環(huán)境下對某一刺激具有不同反應(yīng)的大腦區(qū)域。

一種常見的fMRI實驗數(shù)據(jù)分析方法是廣義線性模型(General Linear Model, GLM)[1]，它通過血液動力學(xué)響應(yīng)函數(shù)(Hemodynamic Response Function, HRF)反映不同大腦節(jié)點對不同類型刺激的響應(yīng)情況。在此基礎(chǔ)上又發(fā)展了參數(shù)模型、非參模型、半?yún)⒛Ｐ偷榷喾NHRF估測方法，其中參數(shù)方法主要利用有限模型參數(shù)來刻畫不同HRF之間的差異特征，具體包括泊松模型[2]和經(jīng)典血液動力學(xué)響應(yīng)函數(shù)模型(Canonical HRF, Canonical)[3]等；非參數(shù)方法通常將HRF表示為一系列基函數(shù)線性組合的形式，因此在衡量不同大腦區(qū)域HRF函數(shù)特征時更加靈活，其主要包括平滑有限脈沖響應(yīng)(Smooth Finite Impulse Response, SFIR)[4]和正則化和廣義交叉驗證(Tikhonov-regularization and Generalized-Cross-Validation, Tik-GCV)[5]；此外，基于大腦活動的群體一般性和個體特異性，Zhang等[6-7]建立了更為靈活的半?yún)?shù)模型。

上述基于GLM框架的HRF計算模型都是單節(jié)點方法，即一次只對一個大腦節(jié)點的fMRI時間序列進(jìn)行分析。由于空間上相鄰的節(jié)點往往具有相似的fMRI數(shù)據(jù)，因此將大腦的空間信息并入HRF計算中將提升模型效果。基于同一區(qū)塊大腦節(jié)點享有共同的函數(shù)形狀而僅僅在振幅上有所不同這一假設(shè)，Vincent等[8]在大腦區(qū)塊貝葉斯模型[9-10]的基礎(chǔ)上提出HRF振幅空間先驗；Chaari等[11]進(jìn)一步提出了同一大腦區(qū)塊中隨空間變化的HRF振幅先驗，并開發(fā)了聯(lián)合分割檢測估計程序。同時，fMRI研究向多個體、多維度數(shù)據(jù)發(fā)展，例如Degras等[12]提出了多個體功能磁共振成像HRF估計的貝葉斯模型，Zhang等[13-14]建立了適用于fMRI時間序列中復(fù)雜時間空間相關(guān)性的貝葉斯方法。

現(xiàn)有的基于GLM框架的模型在進(jìn)行大腦反應(yīng)比較實驗時，通常需要提取HRF的某些低維特征(如高度等)，利用假設(shè)檢驗對這些低維特征進(jìn)行比較，進(jìn)而獲取兩次實驗中反應(yīng)不同的大腦節(jié)點。然而這種方法只能比較HRF某一特征(如高度)的差異而忽略了其他特征(如函數(shù)形狀)的差異，導(dǎo)致估測結(jié)果具有很大的不確定性。而假設(shè)檢驗本質(zhì)上是對每一個節(jié)點進(jìn)行分析，因此忽略了fMRI數(shù)據(jù)的空間特性，導(dǎo)致被識別節(jié)點的準(zhǔn)確度較低。針對這一缺陷，Zhang等[15]提出了綜合所有腦節(jié)點空間信息的低階多元廣義線性模型(Low-Rank Multivariate General Linear Model, LRMGLM)，該模型利用空間矩陣靈活描述了HRF的變化，識別節(jié)點準(zhǔn)確率更高，但fMRI數(shù)據(jù)的低信噪比和高變異性導(dǎo)致其計算復(fù)雜、計算效率低下。此外，由于LRMGLM只能對同一實驗的大腦反應(yīng)進(jìn)行比較，因此在比較不同實驗的大腦反應(yīng)時，需要額外設(shè)置一組相同的刺激作為參考項，適用范圍有限。

對此，本文提出一種基于擴(kuò)展的LRMGLM(Extended LRMGLM, ELRMGLM)的腦節(jié)點識別方法。該方法建立了可同時處理兩次實驗所有節(jié)點數(shù)據(jù)的ELRMGLM模型，該模型通過將血液動力學(xué)響應(yīng)函數(shù)(HRF)的特征矩陣轉(zhuǎn)化為兩個低階矩陣相乘的形式，在實現(xiàn)腦節(jié)點靈活比較的同時簡化參數(shù)運(yùn)算。模型參數(shù)利用基于fMRI數(shù)據(jù)時空特性的優(yōu)化函數(shù)和迭代算法進(jìn)行估測，同時為降低模型對fMRI數(shù)據(jù)高變異性和低信噪比的敏感度，開發(fā)了基于K-means的快速選擇策略來實現(xiàn)兩次實驗中反應(yīng)不同大腦區(qū)域的快速選擇。

1 本文模型建立

(1)

其中:

以往基于GLM框架的方法一次只對一個大腦節(jié)點或一次實驗數(shù)據(jù)進(jìn)行單獨處理，本文提出一個可以同時處理兩次實驗所有節(jié)點數(shù)據(jù)的聯(lián)合模型，以更多的時空信息來減少fMRI數(shù)據(jù)中噪聲的干擾。由于不同個體大腦節(jié)點的HRF形狀不同，首先利用B-樣條插值對其進(jìn)行擬合：

(2)

(3)

(4)

相較于LRMGLM，ELRMGLM可以直接處理兩次實驗所有節(jié)點的數(shù)據(jù)信息，不僅比原模型適用范圍更廣，而且可以利用更多的時間空間信息來提高識別的準(zhǔn)確度。相較于傳統(tǒng)的單節(jié)點模型，ELRMGLM可以通過保持相應(yīng)Uk(q)不變，對Vk(q)進(jìn)行比較的方法來實現(xiàn)對HRF靈活而全面的比較。例如，在計算過程中可以令U2(1)=U2(2)，比較V2(1)和V2(2)的估計值來識別兩次實驗中對第二種刺激有不同反應(yīng)的大腦區(qū)域。

2 參數(shù)估計

2.1 優(yōu)化函數(shù)

在利用上述模型對數(shù)據(jù)進(jìn)行分析時，通常希望估計得到的Y值盡可能地接近真實值，令

q=1,2}

則有如下代價函數(shù)：

SSE(Θ)=

由于HRF在時間上連續(xù)，為避免過度擬合，得到矩陣Uk(q)上的時間平滑懲罰項：

(5)

同時由于空間上相鄰的大腦節(jié)點通常有相似的fMRI時間序列和HRF函數(shù)，因此得到矩陣Vk(q)上的空間平滑懲罰項：

(6)

τPS(Vk(q))]

(7)

2.2 迭代算法

1)給定V，找到U、β，最小化

2)給定U、β，找到V，最小化

上述迭代算法的計算效率主要取決于兩個最優(yōu)子問題的計算效率，由于步驟1)、2)中的目標(biāo)函數(shù)均為二次函數(shù)且存在最優(yōu)解析表達(dá)式，因此可對其進(jìn)行直接求解，具體推導(dǎo)過程在此不多作贅述。當(dāng)數(shù)據(jù)量較大、數(shù)據(jù)維度過高時，直接求解可能耗時較長，此時可利用最速梯度下降法加速運(yùn)算。

3 基于K-means的快速選擇策略

懲罰參數(shù)λ和τ分別控制著HRF的時間平滑性和空間平滑性。在懲罰優(yōu)化問題中，普通交叉驗證(Ordinary Cross Validation, OCV)和廣義交叉驗證(Generalized Cross Validation, GCV)[18]是選擇懲罰參數(shù)的經(jīng)典方法；在成像數(shù)據(jù)分析中，文獻(xiàn)[19-20]提出了基于GCV的選擇過程，同時文獻(xiàn)[21]對約束最大似然法(Restricted Maximum Likelihood, REML)進(jìn)行了研究。由于本文節(jié)點數(shù)量較多且有兩個懲罰參數(shù)，OCV耗時過長而REML不能直接適用，因此提出了基于K-means的快速選擇策略。

由于本文的研究重點是比較兩次實驗中反應(yīng)不同的大腦節(jié)點而非計算HRF的具體值，即只需要選擇能清楚區(qū)分反應(yīng)相同和反應(yīng)不同腦節(jié)點的懲罰參數(shù)組合，因此可利用聚類方法加快懲罰參數(shù)和大腦節(jié)點的選擇過程，將輪廓系數(shù)作為選擇標(biāo)準(zhǔn)，在提高模型對懲罰參數(shù)容忍度的同時保證所選節(jié)點的準(zhǔn)確性，具體過程如下：

1)對每個懲罰參數(shù)，在e-1～e5范圍內(nèi)選取大量候選參數(shù)值，利用2.2節(jié)中提到的迭代算法對帶有不同懲罰參數(shù)組合的ELRMGLM進(jìn)行模型參數(shù)估計，得到相應(yīng)的時間矩陣Uk(q)和空間矩陣Vk(q)。

2)計算每組懲罰參數(shù)對應(yīng)矩陣Vk(q)在兩次實驗中的差值S=|Vk(1)-Vk(2)|，其中k為所要比較的刺激類型。

3)對每組懲罰參數(shù)，利用K-means聚類將所有節(jié)點的S值分為兩類，選擇輪廓系數(shù)最大的懲罰參數(shù)組合和聚類結(jié)果，其中S均值較大的群組為兩次實驗中對刺激k反應(yīng)不同的大腦節(jié)點集合。

上述快速選擇策略的原理在于兩方面：首先，大部分懲罰參數(shù)組合對應(yīng)的大腦節(jié)點S值可以被自然地分為兩部分，一部分在一個較小值附近上下浮動，另一部分的S值則顯著較大，它們分別代表了由生理噪聲引起的低頻漂移和兩次實驗中不同的大腦反應(yīng)，因此可通過聚類方法加以區(qū)分；其次，不同的懲罰參數(shù)組合影響了兩組節(jié)點S值的差異性以及各組節(jié)點S值的穩(wěn)定性，從而影響了聚類結(jié)果的準(zhǔn)確性，而輪廓系數(shù)是一種評價聚類效果好壞的方式，它衡量了個體相較其他群集與其所屬群集的相似程度，其范圍從-1到1，輪廓系數(shù)值越大表明聚類效果越好，因此可將輪廓系數(shù)作為懲罰參數(shù)的選擇標(biāo)準(zhǔn)。一組輪廓系數(shù)較高的懲罰參數(shù)組合會使大腦節(jié)點的S值具有高差異性和高穩(wěn)定性，聚類效果較好，從而保證了被選節(jié)點的準(zhǔn)確性和可靠性。

聚類方法提高了算法對不同懲罰參數(shù)組合的容忍度，高輪廓系數(shù)保證了所選節(jié)點的準(zhǔn)確性，通過采用基于K-means的快速選擇策略，利用聚類加快懲罰參數(shù)和大腦節(jié)點的選擇過程，可以在保證準(zhǔn)確性的同時快速找到兩次實驗中反應(yīng)不同的大腦區(qū)域。

4 實驗分析

4.1 實驗數(shù)據(jù)

本文采用與文獻(xiàn)[22-24]中真實fMRI實驗相同的實驗設(shè)計進(jìn)行分析，該實驗使用了以頻率27.9%、21.1%、50.8%和0.2%隨機(jī)出現(xiàn)的4種不同刺激，共有106個受試者，每個受試者的fMRI時間序列包含205次掃描，每次掃描時間為2 s。考慮到LRMGLM只適用于存在參考項的實驗數(shù)據(jù)，本文采用文獻(xiàn)[15]中使用的實驗數(shù)據(jù)集。該數(shù)據(jù)集共含有三組fMRI數(shù)據(jù)，分布在15×15×15的大腦網(wǎng)格上。在組1中，HRF函數(shù)服從經(jīng)典形式且前兩種刺激的HRF在振幅和延遲上均不同；在組2中，HRF函數(shù)形狀改變，其余與組1保持相同；在組3中，改變第二種刺激的振幅使得前兩種刺激的HRF僅僅在延遲時間上有所不同，其余與組2保持相同。此外，該數(shù)據(jù)集中與第二種刺激有關(guān)的HRF參數(shù)只針對中心9×9×9的節(jié)點網(wǎng)格，其余大腦節(jié)點的HRF參數(shù)與第一種刺激相同。為比較實際問題中大腦節(jié)點在兩次實驗的中不同反應(yīng)，本文增設(shè)含三組相同fMRI數(shù)據(jù)的實驗二，并令其中與第二種刺激有關(guān)的HRF參數(shù)只針對中心偏右的5×9×9大腦網(wǎng)格，其余與實驗一保持相同。此時約9.6%的大腦節(jié)點在兩次實驗中對第二個刺激反應(yīng)不同。

值得注意的是，本文采用的實驗設(shè)置并不完全遵循提出的ELRMGLM，盡管如此，實驗結(jié)果將表明ELRMGLM方法能夠超越更適合該實驗數(shù)據(jù)的單節(jié)點分析方法和LRMGLM。

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

本文將提出的ELRMGLM運(yùn)用于實驗數(shù)據(jù)，為保持運(yùn)算簡潔，令P=2來刻畫不同個體大腦節(jié)點HRF在振幅和延遲上的差異。圖1展示了兩次實驗中對第二種刺激反應(yīng)不同的大腦區(qū)域，其中深色區(qū)域表示被選中節(jié)點，白色區(qū)域表示未選節(jié)點，X、Y、Z分別表示大腦節(jié)點的三維坐標(biāo)。可以看到，ELRMGLM方法選中了大腦中心偏左的4×9×9網(wǎng)格，準(zhǔn)確度較高。

圖1 被選中節(jié)點圖像Fig. 1 Image of selected voxels

圖2 三組實驗數(shù)據(jù)的ROC曲線Fig. 2 ROC curves of three experimental datasets

下面將ELRMGLM分別與canonical方法[3]、SFIR方法[4]和Tik-GCV方法[5]進(jìn)行比較。在進(jìn)行計算時，先用上述單節(jié)點方法對HRF進(jìn)行估計，再提取得到的HRF低維特征(如高度)，通過t假設(shè)檢驗對提取的低維特征進(jìn)行比較，進(jìn)而獲得兩次實驗中反應(yīng)不同的大腦區(qū)域。通過改變t檢驗中的P值，得到了顯示不同臨界值對應(yīng)的真正類率(True Positive Rate, TPR)和負(fù)正類率(False Positive Rate, FPR)組合的受試者工作特征曲線(Receiver Operator characteristic Curve, ROC)。為進(jìn)行比較，本文還通過使用不同的懲罰參數(shù)組合畫出了ELRMGLM的ROC曲線，具體如圖2所示。

從圖2可知：Canonical方法在第一組和第三組實驗數(shù)據(jù)中表現(xiàn)較差，但在第二組實驗數(shù)據(jù)中表現(xiàn)較好，其TPR和FPR分別達(dá)到了約80%和4%。與之相反，SFIR在第一組和第三組數(shù)據(jù)中表現(xiàn)較好，其TPR和FPR分別達(dá)到了約95%、3%和99%、20%，但在第二組數(shù)據(jù)中表現(xiàn)一般。Tik-GCV在所有單節(jié)點分析方法中表現(xiàn)最為穩(wěn)定，其TPR和FPR在三組數(shù)據(jù)中均達(dá)到99%和20%左右。通過使用第3章基于K-means的快速選擇策略，本文提出的ELRMGLM在三組實驗數(shù)據(jù)中均表現(xiàn)優(yōu)異，其TPR和FPR均實現(xiàn)了99%以上和1%以下 (其TPR和FPR分別達(dá)到99.73%、0.44%，99.99%、0.09%和99.99%、0.01%)，比以上三種方法的最優(yōu)結(jié)果分別提升了約20%、8%、20%，不僅實現(xiàn)了高敏感度(sensitivity)和高特異度(specificity)，同時在不同的數(shù)據(jù)集上表現(xiàn)穩(wěn)定。

下面以第一種刺激作為參考項，對ELRMGLM和LRMGLM方法進(jìn)行比較，其結(jié)果如表1所示。可以看到在三組實驗中，ELRMGLM對兩次實驗中反應(yīng)不同的大腦節(jié)點識別準(zhǔn)確度略高于LRMGLM。在計算時間上，ELRMGLM的平均迭代次數(shù)少于100而LRMGLM需要迭代上萬次；同時ELRMGLM的單次迭代時間為8 s，遠(yuǎn)小于LRMGLM的單次迭代時間60 s。不難算出，ELRMGLM的計算時間是LRMGLM的1/750，算法效率大幅提高。

表1 ELRMGLM與LRGMLM方法比較Tab. 1 Comparison of ELRMGLM and LRGMLM

5 結(jié)語

本文提出了一種用于識別兩次實驗中反應(yīng)不同大腦區(qū)域的擴(kuò)展的低階多元廣義線性模型(ELRMGLM)。該模型同時綜合了兩次實驗的數(shù)據(jù)信息，通過帶懲罰項的優(yōu)化函數(shù)考慮了fMRI數(shù)據(jù)的時空特性，并利用K-means聚類提高了模型對參數(shù)的容忍度，進(jìn)而實現(xiàn)了對大腦節(jié)點的快速準(zhǔn)確識別。通過在三組實驗數(shù)據(jù)集上的分析，該模型在準(zhǔn)確度、計算效率和穩(wěn)定性方面均高于現(xiàn)有模型。ELRMGLM主要用于群體大腦活動的評估，對個體大腦反應(yīng)的比較還稍有欠缺，因此如何擴(kuò)展ELRMGLM使其適用于大腦活動的個體差異將成為下一步的研究重點。