fMRI血液動(dòng)力學(xué)響應(yīng)模型的參數(shù)敏感性分析及優(yōu)化*

楊增宇，蔣文濤，陳宇

（四川大學(xué)工程力學(xué)系，成都610065）

1 引 言

近年來(lái)，功能磁共振成像（fMRI）已經(jīng)成為研究大腦動(dòng)作方式的首要選擇方法［1-3］。其原理主要是，BOLD-FMRI圖像以體素為單位把神經(jīng)元活動(dòng)引起的生理活動(dòng)記錄下來(lái)，針對(duì)每個(gè)體素，我們觀(guān)測(cè)得到的FMRI數(shù)據(jù)相當(dāng)于一個(gè)時(shí)間序列，如果將神經(jīng)元的活動(dòng)引起的生理活動(dòng)的變化統(tǒng)稱(chēng)為神經(jīng)活動(dòng)的血液動(dòng)力學(xué)響應(yīng)，那么可以采用如下的方式來(lái)描述：觀(guān)測(cè)到的體素時(shí)間序列＝血液動(dòng)力學(xué)響應(yīng)＋漂移＋噪聲［4］。因此，要探究fMRI數(shù)據(jù)與神經(jīng)活動(dòng)的關(guān)系，研究血液動(dòng)力學(xué)響應(yīng)是關(guān)鍵。如果用X（t）表示血流動(dòng)力學(xué)響應(yīng)，u（t）表示神經(jīng)活動(dòng)，h（t）表示相應(yīng)的血液動(dòng)力學(xué)響應(yīng)函數(shù)，那么在線(xiàn)性時(shí)不變系統(tǒng)模型［5-6］中 X（t），u（t），h（t）之間的關(guān)系可以表示為：X（t）＝u（t）×h（t），即血流動(dòng)力學(xué)響應(yīng)等于神經(jīng)活動(dòng)與血液動(dòng)力學(xué)響應(yīng)函數(shù)的卷積。因此一旦對(duì)血液動(dòng)力學(xué)響應(yīng)函數(shù)有了相當(dāng)?shù)恼J(rèn)識(shí)就可以預(yù)測(cè)由某種神經(jīng)活動(dòng)模式所導(dǎo)致的血流動(dòng)力學(xué)響應(yīng)進(jìn)而得到體素時(shí)間序列，實(shí)現(xiàn)fMRI成像功能。

2 標(biāo)準(zhǔn)血液動(dòng)力學(xué)響應(yīng)模型

由于血液動(dòng)力學(xué)響應(yīng)函數(shù)直接影響觀(guān)測(cè)的體素時(shí)間序列，所以對(duì)于血液動(dòng)力學(xué)響應(yīng)函數(shù)的研究具有重要意義。目前三參數(shù)gamma函數(shù)形式是fMRI數(shù)據(jù)處理的研究中應(yīng)用最為廣泛的血液動(dòng)力學(xué)響應(yīng)函數(shù)之一。1994年Friston等［7］提出可以用gamma函數(shù)的形式來(lái)表征血液動(dòng)力學(xué)響應(yīng)函數(shù)，在此基礎(chǔ)上Cohen［8］于1996年提出了三參數(shù)形式的血液動(dòng)力學(xué)響應(yīng)函數(shù)模型：

眾多研究［9］已經(jīng)表明上述標(biāo)準(zhǔn)的血液動(dòng)力學(xué)響應(yīng)函數(shù)模型在描述和預(yù)測(cè)血液動(dòng)力學(xué)響應(yīng)時(shí)具有一定的合理性和有效性，然而由于在實(shí)驗(yàn)數(shù)據(jù)擬合時(shí)未考慮參數(shù)間的耦合作用，導(dǎo)致得到的函數(shù)模型在求解精度和擬合效果上都并非最優(yōu)，迄今國(guó)內(nèi)外針對(duì)血液動(dòng)力學(xué)響應(yīng)函數(shù)的研究幾乎沒(méi)有對(duì)該三參數(shù)模型進(jìn)行相應(yīng)的優(yōu)化處理工作。本研究將以此模型為基本框架，通過(guò)擬合實(shí)驗(yàn)數(shù)據(jù)值對(duì)其進(jìn)行優(yōu)化。

3 參數(shù)敏感性分析

由3參數(shù)控制的血液動(dòng)力學(xué)響應(yīng)模型，由于各參數(shù)對(duì)最終信號(hào)響應(yīng)值的影響程度是不同的：影響大的參數(shù)即使有微小攝動(dòng)，都會(huì)引起流動(dòng)應(yīng)力的較大波動(dòng)，甚至吞沒(méi)了其他相對(duì)不敏感參數(shù)的作用。因此，有必要對(duì)由實(shí)驗(yàn)數(shù)據(jù)確定的血液動(dòng)力學(xué)響應(yīng)模型的參數(shù)進(jìn)行敏感度分析。

3.1 參數(shù)敏感度分析方法

拉丁超立方抽樣（LHS）是一種均勻空間抽樣方法，最早于1979年由simpson等提出，具有抽樣次數(shù)較少、可有效避免重復(fù)抽樣等特點(diǎn)，適用于對(duì)復(fù)雜的多維參數(shù)空間抽樣［10］。抽樣過(guò)程如下：

（1）設(shè)變量X的取值范圍為［a，b］，所需抽取的樣本個(gè)數(shù)為n。

（2）將X的取值范圍［a，b］等概率的分為 n個(gè)區(qū)間，并在每個(gè)區(qū)間中隨機(jī)的抽取一個(gè)數(shù)作為Xi的一個(gè)樣本值，將所取的值進(jìn)行隨機(jī)排列，由此可得到關(guān)于Xi的樣本容量為n的樣本序列：｛x1，x2，x3……xn｝。

（3）對(duì)于m維參數(shù)空間（即含有m個(gè)變量）的抽樣，則是對(duì) m個(gè)變量（X1，X2，X3……Xm）均進(jìn)行步驟1和步驟2，得到m個(gè)樣本序列。

（4）對(duì)m個(gè)樣本序列各取一個(gè)值即可獲得一個(gè)m維的樣本參數(shù)集，共可得到n個(gè)m維的樣本參數(shù)集，見(jiàn)表1，此即完成了拉丁超立方抽樣的全過(guò)程。

表1 多參數(shù)LHS抽樣樣本參數(shù)集Table 1 Sample parameter set of Multi parameter LHS sampling

Spearman秩相關(guān)系數(shù)是一個(gè)非參數(shù)性質(zhì)（與分布無(wú)關(guān)）的秩統(tǒng)計(jì)參數(shù)［11］，由 Spearman在1904年提出，用于度量?jī)蓚€(gè)隨機(jī)變量之間的相關(guān)方向與密切程度，在統(tǒng)計(jì)學(xué)中，Spearman秩相關(guān)系數(shù)一般由希臘字母ρs表示。設(shè)變量X的樣本集為｛x1，x2，x3……xn｝，將其按照從大到小的順序進(jìn)行排列，得到一個(gè)有序序列｛x1′，x2′，x3′……xn′｝。定義原始 xi在新的有序序列｛x1′，x2′，x3′……xn′｝中所占的位置Ki為 xi在樣本集｛x1，x2，x3……xn｝中的秩。設(shè)與 X對(duì)應(yīng)的變量Y的樣本參數(shù)集為｛y1，y2，y3……yn｝。設(shè) xk在｛x1，x2，x3……xn｝中的秩為 αk，yk在｛y1，y2，y3……yn｝中的秩為 βk，則相應(yīng)的 spearman秩相關(guān)系數(shù)的計(jì)算公式為：

Spearman秩相關(guān)系數(shù)即反映了該參數(shù)的敏感度［8］。

3.2 參數(shù)敏感度分析結(jié)果

對(duì)模型參數(shù)（k、n、m）進(jìn)行超立方抽樣，形成 p個(gè)3維的樣本參數(shù)集，見(jiàn)表2，每個(gè)參數(shù)集對(duì)應(yīng)一個(gè)具體血液動(dòng)力學(xué)響應(yīng)函數(shù)模型。

表2 血液動(dòng)力學(xué)響應(yīng)模型的LHS抽樣樣本參數(shù)集Table 2 LHS sampling parameter set of hemodynamic response model

將所得到的各個(gè)樣本參數(shù)集的參數(shù)帶入模型中，依據(jù)文獻(xiàn)［8］中實(shí)驗(yàn)的時(shí)間節(jié)點(diǎn)值計(jì)算出相應(yīng)的信號(hào)強(qiáng)度h*，而文獻(xiàn)［8］中實(shí)驗(yàn)測(cè)得的信號(hào)強(qiáng)度為h，根據(jù)最小二乘擬合的原理，計(jì)算目標(biāo)函數(shù)值W，依據(jù)spearman秩的定義和式（1），利用所得的 W值分別和相應(yīng)的 k、n、m值計(jì)算spearman秩相關(guān)系數(shù)，從而得到參數(shù)敏感度，見(jiàn)表3。

表3 模型參數(shù)敏感度統(tǒng)計(jì)表Table 3 Statistics of model parameter sensitivity

由表3中數(shù)據(jù)可知：（1）當(dāng)各參數(shù)的取值范圍一定時(shí)，響應(yīng)函數(shù)的參數(shù)敏感性估計(jì)受抽樣次數(shù)的影響，并且抽樣次數(shù)越大，各參數(shù)的敏感度越趨于穩(wěn)定。所以為了取得較真實(shí)、較準(zhǔn)確的模型參數(shù)敏感度，應(yīng)將抽樣次數(shù)取得足夠大，以保證所得到的敏感性分析結(jié)果有效性及合理性。當(dāng)抽樣次數(shù)由500次增加至1 000次時(shí)參數(shù)敏感度波動(dòng)很小，故可認(rèn)為500次抽樣結(jié)果已經(jīng)足夠穩(wěn)定，能夠代表樣本特征，此時(shí)的數(shù)據(jù)具有分析和應(yīng)用價(jià)值。（2）當(dāng)保持參數(shù)k、m取值范圍以及抽樣次數(shù)一定（均為500次），只改變參數(shù)n的取值范圍時(shí)，各參數(shù)敏感度均存在較大的波動(dòng)，故各參數(shù)敏感度不僅與其自身的取值范圍有關(guān)，還與其他參數(shù)的取值范圍相關(guān)，即模型參數(shù)間的耦合作用對(duì)分析結(jié)果的影響很大。（3）三個(gè)參數(shù)k、m、n中，m和 n的參數(shù)敏感度相當(dāng)，而 k的參數(shù)敏感度較小，與m、n相比有很大差異。即m、n對(duì)目標(biāo)函數(shù)的影響較大，為了使各個(gè)參數(shù)都有相當(dāng)?shù)膮?shù)敏感度，應(yīng)該使各個(gè)參數(shù)對(duì)目標(biāo)函數(shù)的spearman秩相關(guān)系數(shù)都接近0.333，然而這樣又會(huì)將k值的取值區(qū)間設(shè)置得過(guò)小，不利于后續(xù)的優(yōu)化搜索工作，故我們?nèi)”碇衚、m、n參數(shù)敏感度最接近的情況，即取初始取值區(qū)間為：k：［0，10］m：［0，40］n：［-20，20］

4 優(yōu)化與分析

4.1 模型優(yōu)化算法—遺傳算法

遺傳算法是一種模仿生物的遺傳進(jìn)化原理，通過(guò)選擇、交叉與變異等操作機(jī)制，使種群中個(gè)體的適應(yīng)性不斷提高，從而找到問(wèn)題的滿(mǎn)意解或最優(yōu)解的仿生全局優(yōu)化算法［12］。

遺傳算法基本流程框見(jiàn)圖1。在參數(shù)的初始取值區(qū)間范圍內(nèi)隨機(jī)取值，作為一個(gè)個(gè)體，將所有初始個(gè)體組成的群體稱(chēng)為初始種群，對(duì)初始種群中每個(gè)個(gè)體進(jìn)行二進(jìn)制編碼，所得到的二進(jìn)制字符串通常稱(chēng)為該個(gè)體的“染色體”；設(shè)置合理的適應(yīng)度函數(shù)，計(jì)算初始種群的適應(yīng)度值，并進(jìn)行是否終止進(jìn)化的判斷；對(duì)于種群需要多次進(jìn)化才能滿(mǎn)足求解要求的問(wèn)題，對(duì)上一代種群進(jìn)行選擇（輪盤(pán)賭方法選擇新一代進(jìn)化種群）、交叉（兩兩個(gè)體之間按照交叉概率交換他們的部分染色體片段）、變異（每個(gè)個(gè)體按照變異概率改變某一個(gè)或某些基因座上的基因值為其它的等位基因）等運(yùn)算得到新一代種群，并判斷是否滿(mǎn)足終止進(jìn)化條件。

4.2 模型優(yōu)化分析

根據(jù)參數(shù)敏感性分析結(jié)果，對(duì)模型中的3個(gè)參數(shù)給定初始取值區(qū)間：k：［0，10］，m：［0，40］，n：［-20，20］，利用matlab平臺(tái)編寫(xiě)遺傳算法程序，以文獻(xiàn)［8］中實(shí)驗(yàn)值為基礎(chǔ)，分別取遺傳代數(shù)為50、100、500、1 000對(duì)模型參數(shù)進(jìn)行最優(yōu)值求解，最終，求得使目標(biāo)函數(shù)W（x）最小時(shí)所對(duì)應(yīng)的參數(shù)即響應(yīng)模型：

圖1 遺傳算法基本流程框圖Fig 1 Basic flow diagram of genetic algorithm

50代：［k，m，n］＝［21.654，1.878，-0.891］h（t）＝21.654t1.878e-0.891t

100代：［k，m，n］＝［17.017，7.393，-2.358］h（t）＝17.017t7.393e-2.358t

500代：［k，m，n］＝［2.095，6.699，-1.538］h（t）＝2.095t6.699e-1.538t

1000代：［k，m，n］＝［0.518，7.780，-1.601］h（t）＝0.518t7.780e-1.601t

2000代：［k，m，n］＝［0.457，8.460，-1.804］h（t）＝0.457t8.460e-1.804t

將所得模型與文獻(xiàn)［8］中實(shí)驗(yàn)值進(jìn)行擬合效果比較：

圖2 不同遺傳代數(shù)擬合效果比較Fig 2 Comparison of fitting effect inthe different genetic algebra

由圖2可知隨著遺傳代數(shù)的增加，通過(guò)優(yōu)化搜索所得的模型的擬合效果越來(lái)越好；算法顯示出較快的收斂性，當(dāng)遺傳代數(shù)為100代時(shí)所得模型已經(jīng)初步具有實(shí)驗(yàn)值曲線(xiàn)的輪廓；當(dāng)遺傳代數(shù)多于1 000代時(shí)，代與代之間差距較小，已經(jīng)接近最優(yōu)的模型曲線(xiàn)。

4.3 模型優(yōu)化結(jié)果

取遺傳代數(shù)為2500代（可以認(rèn)為此時(shí)目標(biāo)函數(shù)已收斂到最小值），進(jìn)行遺傳算法最優(yōu)搜索，得到優(yōu)化的模型為：

而根據(jù)文獻(xiàn)［8］中實(shí)驗(yàn)值，由Cohen利用最小二乘法確定的三參數(shù)形式的血液動(dòng)力學(xué)響應(yīng)函數(shù)模型為：

圖3 Cohen模型與本研究模型的比較Fig 3 Comparison of the Cohen model with the model presented in this study

由圖3可知Cohen模型與本研究所建立的模型都能夠描述相應(yīng)的血流動(dòng)力學(xué)響應(yīng)，但是本研究所建立的模型能夠更加貼合實(shí)驗(yàn)曲線(xiàn)，具有更高的求解精度與更好的擬合效果。

5 經(jīng)典算例分析

為了說(shuō)明優(yōu)化的響應(yīng)模型能夠較好的描述與預(yù)測(cè)由神經(jīng)活動(dòng)（u）所引起的血液動(dòng)力學(xué)響應(yīng)（X），本研究以一個(gè)OFF-ON任務(wù)為例，任務(wù)數(shù)據(jù)取自文獻(xiàn)［13］Savoy及其學(xué)生做的視覺(jué)閃光實(shí)驗(yàn)，實(shí)驗(yàn)采用方波刺激，有交替的45 s的刺激和休息時(shí)間，實(shí)驗(yàn)記錄下fMRI信號(hào)隨時(shí)間的變化過(guò)程。利用Cohen模型以及所得到的優(yōu)化后的響應(yīng)模型式（4）來(lái)進(jìn)行比較與驗(yàn)證，可以得到圖4的血液動(dòng)力學(xué)響應(yīng)，圖中灰色區(qū)域表示刺激期，白色區(qū)域表示休息期，連續(xù)曲線(xiàn)為神經(jīng)活動(dòng)（此處為δ（t）階躍函數(shù)）和響應(yīng)函數(shù)的卷積的形成的響應(yīng)預(yù)測(cè)曲線(xiàn)?？梢钥吹?，（1）響應(yīng)信號(hào)相對(duì)任務(wù)的開(kāi)始時(shí)間有一定的延遲，在任務(wù)結(jié)束后也需要一定的時(shí)間回到基線(xiàn)水平，這與通常觀(guān)測(cè)到的fMRI信號(hào)特點(diǎn)是吻合的。（2）與傳統(tǒng)的Cohen模型相比較，本研究所得優(yōu)化模型能夠更好描述和預(yù)測(cè)由神經(jīng)活動(dòng)所引起的血液動(dòng)力學(xué)響應(yīng)，所得結(jié)果更加接近真實(shí)值。

圖4 OFF-ON任務(wù)響應(yīng)圖Fig 4 OFF-ON task response diagram

6 結(jié)論

血液動(dòng)力學(xué)響應(yīng)函數(shù)模型參數(shù)的敏感度受其他參數(shù)的取值范圍的影響，故不能單獨(dú)研究某一材料參數(shù)的敏感度，表明本研究的LHS和spearman秩相關(guān)結(jié)合的整體分析方法對(duì)于該問(wèn)題是適合的。本研究采用參數(shù)敏感度估計(jì)搜索區(qū)域和遺傳算法結(jié)合的方法能夠很好的實(shí)現(xiàn)血液動(dòng)力學(xué)響應(yīng)函數(shù)模型的參數(shù)識(shí)別。所用于最優(yōu)化處理的遺傳算法具有較好的收斂性與普適性。建立的優(yōu)化模型能夠很好的描述和預(yù)測(cè)由神經(jīng)活動(dòng)引起的血流動(dòng)力學(xué)響應(yīng)，所得到的預(yù)測(cè)曲線(xiàn)與通常觀(guān)測(cè)到的fMRI數(shù)據(jù)特點(diǎn)吻合。此外，本研究建立的優(yōu)化模型在求解精度和擬合效果上較傳統(tǒng)的Cohen模型有一定的提高。