fMRI血液動力學響應模型的參數敏感性分析及優化*

楊增宇，蔣文濤，陳宇

（四川大學工程力學系，成都610065）

1 引 言

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

2 標準血液動力學響應模型

由于血液動力學響應函數直接影響觀測的體素時間序列，所以對于血液動力學響應函數的研究具有重要意義。目前三參數gamma函數形式是fMRI數據處理的研究中應用最為廣泛的血液動力學響應函數之一。1994年Friston等［7］提出可以用gamma函數的形式來表征血液動力學響應函數，在此基礎上Cohen［8］于1996年提出了三參數形式的血液動力學響應函數模型：

眾多研究［9］已經表明上述標準的血液動力學響應函數模型在描述和預測血液動力學響應時具有一定的合理性和有效性，然而由于在實驗數據擬合時未考慮參數間的耦合作用，導致得到的函數模型在求解精度和擬合效果上都并非最優，迄今國內外針對血液動力學響應函數的研究幾乎沒有對該三參數模型進行相應的優化處理工作。本研究將以此模型為基本框架，通過擬合實驗數據值對其進行優化。

3 參數敏感性分析

由3參數控制的血液動力學響應模型，由于各參數對最終信號響應值的影響程度是不同的：影響大的參數即使有微小攝動，都會引起流動應力的較大波動，甚至吞沒了其他相對不敏感參數的作用。因此，有必要對由實驗數據確定的血液動力學響應模型的參數進行敏感度分析。

3.1 參數敏感度分析方法

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

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

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

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

（4）對m個樣本序列各取一個值即可獲得一個m維的樣本參數集，共可得到n個m維的樣本參數集，見表1，此即完成了拉丁超立方抽樣的全過程。

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

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

Spearman秩相關系數即反映了該參數的敏感度［8］。

3.2 參數敏感度分析結果

對模型參數（k、n、m）進行超立方抽樣，形成 p個3維的樣本參數集，見表2，每個參數集對應一個具體血液動力學響應函數模型。

表2 血液動力學響應模型的LHS抽樣樣本參數集Table 2 LHS sampling parameter set of hemodynamic response model

將所得到的各個樣本參數集的參數帶入模型中，依據文獻［8］中實驗的時間節點值計算出相應的信號強度h*，而文獻［8］中實驗測得的信號強度為h，根據最小二乘擬合的原理，計算目標函數值W，依據spearman秩的定義和式（1），利用所得的 W值分別和相應的 k、n、m值計算spearman秩相關系數，從而得到參數敏感度，見表3。

表3 模型參數敏感度統計表Table 3 Statistics of model parameter sensitivity

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

4 優化與分析

4.1 模型優化算法—遺傳算法

遺傳算法是一種模仿生物的遺傳進化原理，通過選擇、交叉與變異等操作機制，使種群中個體的適應性不斷提高，從而找到問題的滿意解或最優解的仿生全局優化算法［12］。

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

4.2 模型優化分析

根據參數敏感性分析結果，對模型中的3個參數給定初始取值區間：k：［0，10］，m：［0，40］，n：［-20，20］，利用matlab平臺編寫遺傳算法程序，以文獻［8］中實驗值為基礎，分別取遺傳代數為50、100、500、1 000對模型參數進行最優值求解，最終，求得使目標函數W（x）最小時所對應的參數即響應模型：

圖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

將所得模型與文獻［8］中實驗值進行擬合效果比較：

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

由圖2可知隨著遺傳代數的增加，通過優化搜索所得的模型的擬合效果越來越好；算法顯示出較快的收斂性，當遺傳代數為100代時所得模型已經初步具有實驗值曲線的輪廓；當遺傳代數多于1 000代時，代與代之間差距較小，已經接近最優的模型曲線。

4.3 模型優化結果

取遺傳代數為2500代（可以認為此時目標函數已收斂到最小值），進行遺傳算法最優搜索，得到優化的模型為：

而根據文獻［8］中實驗值，由Cohen利用最小二乘法確定的三參數形式的血液動力學響應函數模型為：

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

由圖3可知Cohen模型與本研究所建立的模型都能夠描述相應的血流動力學響應，但是本研究所建立的模型能夠更加貼合實驗曲線，具有更高的求解精度與更好的擬合效果。

5 經典算例分析

為了說明優化的響應模型能夠較好的描述與預測由神經活動（u）所引起的血液動力學響應（X），本研究以一個OFF-ON任務為例，任務數據取自文獻［13］Savoy及其學生做的視覺閃光實驗，實驗采用方波刺激，有交替的45 s的刺激和休息時間，實驗記錄下fMRI信號隨時間的變化過程。利用Cohen模型以及所得到的優化后的響應模型式（4）來進行比較與驗證，可以得到圖4的血液動力學響應，圖中灰色區域表示刺激期，白色區域表示休息期，連續曲線為神經活動（此處為δ（t）階躍函數）和響應函數的卷積的形成的響應預測曲線。可以看到，（1）響應信號相對任務的開始時間有一定的延遲，在任務結束后也需要一定的時間回到基線水平，這與通常觀測到的fMRI信號特點是吻合的。（2）與傳統的Cohen模型相比較，本研究所得優化模型能夠更好描述和預測由神經活動所引起的血液動力學響應，所得結果更加接近真實值。

圖4 OFF-ON任務響應圖Fig 4 OFF-ON task response diagram

6 結論

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