一種用于fMRI的快速動態因果模型反演方法

摘" 要： 針對動態因果模型（DCM）在分析大腦有效連接時面臨的高計算成本問題，提出一種結合廣義線性模型（GLM）和稀疏DCM的算法，即廣義稀疏DCM（GSD）算法。該算法在以下三個方面進行優化：首先，利用傅里葉變換的對稱性將頻域DCM的復數計算轉換為實數計算，降低計算復雜度；其次，應用GLM和濾波技術減少觀測信號的干擾，提高參數估計的準確度；最后，定義新的代價函數來優化變分推斷參數和濾波器參數，進一步提升參數估計的精度。該研究采用兩組公開的功能核磁共振成像（fMRI）數據對GSD算法進行驗證，包括仿真的史密斯小世界網絡數據和運動與注意力實測數據。實驗結果表明，GSD算法在保持與傳統方法相近的參數估計性能的同時，能將計算時間降低50%以上。該研究成果為平衡模型的解釋力和計算效率提供了新的視角，有望推動DCM在更廣泛領域的應用。

關鍵詞： 廣義線性模型； 動態因果模型； 稀疏算法； 有效連接； 計算復雜度； 模型解釋力

中圖分類號： TN919?34； TP301.6" " " " " " " " " "文獻標識碼： A" " " " " " " " " " 文章編號： 1004?373X（2025）07?0146?09

Fast dynamic causal modeling regression method for fMRI

HU Xinhang1， WU Haifeng1， 2， 3， ZENG Yu1， 2， 3

（1. School of Electrical and Information Engineering， Yunnan Minzu University， Kunming 650504， China;

2. Yunnan Provincial Key Laboratory of Unmanned Autonomous Systems， Kunming 650504， China;

3. Yunnan Provincial Colleges and Universities Intelligent Sensor Network and Information System Technology Innovation Team，

Yunnan Minzu University， Kunming 650504， China）

Abstract： In view of the high computational cost faced by dynamic causal modeling （DCM） when analyzing brain effective connectivity， an algorithm that combines the generalized linear model （GLM） and sparse DCM， namely the combined generalize?linear?model and sparse DCM （GSD） algorithm， is proposed. The algorithm is optimized in the following three aspects： firstly， the symmetry of Fourier transform is used to convert the complex number calculation of frequency domain DCM into real number calculation， reducing the computational complexity; secondly， GLM and filtering techniques are applied to reduce the interference of observed signals and improve the accuracy of parameter estimation; finally， a new cost function is defined to optimize the variational inference parameters and filter parameters to further improve the accuracy of parameter estimation. In this study， two sets of public functional magnetic resonance imaging （fMRI） data are used to verify the GSD algorithm， including the simulated Smith small?world network data and the motion and attention measured data. The experimental results show that the GSD algorithm can reduce the calculation time by more than 50% while maintaining parameter estimation performance similar to that of the traditional methods. The research results of this paper provide a new perspective for balancing the explanatory power and computational efficiency of the model， and are expected to promote the application of DCM in a wider range of fields.

Keywords： GLM; DCM; sparse algorithm; effective connectivity; computational complexity; model explanatory power

0" 引" 言

隨著腦成像技術的不斷發展，出現了許多結合功能核磁共振成像（Functional Magnetic Resonance Imaging， fMRI）技術來研究有效連接的方法，例如結構方程模型、多變量回歸模型、格蘭杰因果模型（Granger Causality Modeling， GCM）和動態因果模型（Dynamic Causal Modeling， DCM）等[1]，其中GCM和DCM最具代表性。DCM基于動態系統理論和狀態空間模型[2]，不僅能夠更全面地考慮大腦網絡的動態調節過程，也可以根據神經科學先驗知識對腦活動進行解釋，這使得DCM成為了研究有效連接更主流的方法，并涌現出了許多與之相結合的新方法[3?5]。盡管DCM有著上述諸多優點，但其計算成本還是相對較高，因此，DCM多用于分析腦區較少的小型網絡。

為了提升DCM的運算速度并降低運算成本，不斷有新變體被提出。一類變體是對模型進行簡化[6?7]，例如，簡化神經動力學模型，將模型中的某些生物學過程或連接考慮為常數或固定值，從而減少了需要估計的參數數量；再如簡化模型本身結構，即通過限制模型中的連接類型或拓撲結構來減少待估計的參數數量。盡管這些簡化模型降低了計算復雜度，但也帶來了一些缺點。例如，簡化模型可能會喪失一些生物學上的細節，導致模型的解釋力和預測能力降低。另一類變體是通過改進參數估計和優化計算方法來提高DCM計算速度[8?9]，例如，通過近似推斷和變換域方法等。盡管這些方法在提高計算速度上取得了一定進展，但它們往往需要對模型做出一些假設，從而可能影響到模型的準確性和可靠性。

綜上所述，DCM的解釋力與復雜度是一對矛盾，解釋力強的方法，其復雜度較高；相反，復雜度低的DCM方法雖然計算速度快，但解釋力將變弱。針對以上問題，本文提出了一種結合廣義線性模型（Generalized Linear Model， GLM）[10?11]的稀疏DCM算法，即廣義稀疏DCM（Combined Generalize?Linear?Model and Sparse DCM， GSD）算法。期望在模型解釋力和計算復雜度上得到一個折中，既保證模型有一定的復雜度，同時又有較快的計算速度。

本文的貢獻如下：

1） 利用傅里葉變換的對稱特性，將頻域DCM的復數計算變為實數計算，從而降低了計算復雜度；

2） 利用GLM算法和頻域濾波器，消除了觀測信號的部分干擾，提高了模型參數估計的準確度；

3） 利用定義的代價函數，確定了變分推斷參數及濾波器參數，從而進一步提高參數估計準確度；

4） 使用兩組公開數據[12?14]對GSD算法進行驗證。

1" 算" 法

在本節中，將介紹所提算法的各個模塊，本文方法的流程圖如圖1所示。

1） 基于設計的實驗方案，獲取相應的刺激信號，并構建設計矩陣；隨后，利用GLM算法和統計參數映射（SPM）軟件，精確地識別出大腦的激活區域，并從中提取感興趣區域（Region of Interest， ROI）的血氧水平依賴性（Blood Oxygen Level Dependence， BOLD）信號。由于原始DCM中的神經元活動是不可觀測的隱藏變量，將原始模型轉化為頻域中的線性回歸方程，這一轉換可避免估計神經元隱藏量。

2） 為了進一步降低觀測信號中的噪聲干擾，引入了兩種濾波技術。第一種濾波技術是通過使用GLM估計得到的擬合系數來重構GLM信號，并將該信號替代原始BOLD信號于線性回歸方程中，實現第一次濾波。緊接著，采用頻域濾波器對方程進行第二次濾波處理。經過兩輪濾波處理的線性回歸方程盡管仍然呈現復數形式，但可以通過初等變換的方法將其轉換為實數方程，這一步驟進一步簡化了計算過程。

3） 通過對代價函數的尋優確定出最佳的超參數，并通過變分推斷技術得到最終的DCM估計參數。

1.1" DCM的線性回歸

在本文中，利用BOLD信號的線性卷積模型在原DCM方程兩邊同時卷積血液動力學響應函數（Hemodynamic Response Function， HRF），并借助GLM的擬合系數對原方程進行變換，使BOLD信號的導數[yt]被表示為外源刺激與HRF的卷積與擬合系數的乘積，如圖2所示。在新方程中，[yt]為觀測信號，外源刺激是已知的，這樣就規避了直接求解神經元隱藏量。此外，由于BOLD信號是離散信號，可以對其導數進行微分近似，然后進行傅里葉變換，從而在頻域中得到線性回歸方程。

令[z（r）t]為第[r]個腦區在[t]時刻的神經元活動狀態，則根據動態理論及DCM[2]，該腦區神經元狀態的微分[zt]（為簡化，略去[r]）可表示為：

[zt=A+j=1KujtBjZt+CUt] （1）

式中：[ujt]代表[t]時刻的第[j]個刺激；[Zt=z（r）t∈RR×1]代表由[R]個腦區在[t]時刻的神經元活動構成的列矢量；[Ut=ukt∈RK×1]代表[t]時刻的[K]個刺激構成的列矢量；[A∈R1×R]代表與刺激無關的內生網絡連接系數行矢量；[Bj∈R1×R]代表受第[j]個刺激調節的內生網絡連接系數行矢量；[C∈R1×K]代表受刺激直接影響的神經元活動的系數行矢量。

由線性卷積模型，BOLD信號可表示為：

[yt=zt?ht]" （2）

式中[ht]為HRF。因此對式（1）兩邊均卷積上HRF，可得：

[yt=Aψt+j=1KBjy′tj+Cy″t]" （3）

式中：[ψt=Zt?ht=y（r）t∈RR×1]；[y′tj=ujtbTUt?ht]，[b]為由GLM得到的擬合系數矩陣；[y″t=Ut?ht]。

將BOLD信號[yt]表達為離散形式，則其在第[n]個采樣點上的值為：

[yn=ynΔt]" （4）

式中[Δt]為采樣間隔，即重復時間（TR）。由此[yn]可近似表示為：

[yn≈yn+1-ynΔt]" （5）

將式（5）代入式（3），再對[n=1，2，…，N]個采樣點做數字傅里葉變換，則有：

[Yn=AΨn+j=1KBjY′nj+CY″n]" （6）

式中：[Yn]、[Ψn]、[Y′nj]和[Y″n]分別是[yn]、[ψn]、[y′nj]和[y″n]的[N]點傅里葉變換結果。令：

[Xo=Ψ1Ψ2…ΨNY′11Y′21…Y′N1Y′12Y′22…Y′N2????Y′1KY′2K…Y′NKY″1Y″2…Y″NT∈CN×（R+RK+K）] （7）

[θ=AB1B2… BK CT∈R（R+RK+K）×1] （8）

[Yo=Y1Y2 … YNT∈CN×1] （9）

則式（6）變為：

[Yo=Xoθ] （10）

式中：向量[Yo]由BOLD信號的導數構成，為已知；設計矩陣[Xo]由BOLD信號、刺激及GLM擬合系數構成，也為已知；參數矢量[θ]是待估計的神經連接參數，為未知。因此式（10）就可以看為一個頻域中的線性回歸方程。

1.2" GLM濾波

為提高估計準確度，本節將采用GLM[10]濾波和頻域濾波兩種技術來降噪，如圖3所示。在GLM濾波中，先利用GLM重構BOLD信號，將式（10）中的觀測和設計矩陣中的原BOLD信號替換為重構的GLM信號，得到新的觀測和設計矩陣，如圖3中的虛線框所示。在頻域濾波中，對通過GLM濾波后得到的觀測和設計矩陣再進行理想濾波，如圖3中實線框所示。

一個腦區的GLM輸出信號[yg=γn∈RN×1]可表示為：

[yg=xgβ] （11）

式中：[xg∈RN×R]為GLM設計矩陣；[β∈RR×1]為擬合系數矩陣。將該腦區的觀測BOLD信號[yn]用GLM輸出[γn]來替代，有：

[yn=γn+ξn] （12）

式中[ξn]為BOLD信號與GLM輸出信號之間的誤差。將式（12）代入式（10）中，一個腦區的GLM輸出的傅里葉變換信號[Yg]可表示為：

[Yg=Xgθ+Ξg] （13）

式中：[Yg=Γn∈CN×1]，[Γn]為GLM輸出導數[γn=γn+1-γnΔt]的傅里葉變換；[Ξg]為將式（12）中的[ξn]代入式（10）中時產生的誤差矢量。

相應的設計矩陣[Xg]表示為：

[Xg=Ψg1Ψg2…ΨgNY′11Y′21…Y′N1Y′12Y′22…Y′N2????Y′1KY′2K…Y′NKY″1Y″2…Y″NT∈CN×R+RK+K] （14）

式中：[Ψgn]為將式（7）中[ψn]的BOLD信號[yn]的傅里葉變換替換為GLM輸出[γn]的傅里葉變換后的結果。

1.3" 頻域濾波

本節將詳細描述圖3所示的頻域濾波。

令[Hsη∈{0，1}N×1]為傅里葉變換域上的理想數字濾波器矢量，其第[n-s]個元素表示為：

[Hsηn-s=gsηn-s+g′ηsn-s] （15）

式中[gs=0ηn]是寬為[η]的門函數，即：

[gs=0ηn=1，" " "0lt;n≤η" 0，" " "ηlt;n≤N ] （16）

式中：[s]為平移量；[g′ηsn-s]為門函數[gsηn-s]以軸[n=（N+1）/2]做鏡像后的結果。濾波器[Hsη]如圖4所示，其截止頻率可以通過[s]和[η]來設置。

當[s=0]時，為低通濾波；當[s≥（N+1）2-η]時，為高通濾波；其他情況為帶通濾波。對[Yg]進行濾波，即對式（13）兩邊同乘濾波器[Hsη]，有：

[Yf=Xfθ+Ξf] （17）

式中：[Yf=Yg⊙Hsη*]表示為兩矢量的哈達瑪積后再將0元素去除，因此有[Yf∈CM×1]，其中[M]為哈達瑪積后非0元素的數目。同理，設計矩陣為[Xf=Xg⊙Hsη*∈CM×（R+RK+K）]，經過濾波器后的噪音矢量為[Ξf∈CM×1]。

在式（17）中，濾波后又得到了新的觀測和設計矩陣。需要注意的是，濾波器設計中的參數[s]和[η]其實分別是濾波器低頻和高頻的截止頻率，可通過對代價函數尋優得到，這將在后面的小節中描述。

1.4" 實數化

由于本文將DCM放在頻域中進行求解，則原來的實信號就變為復信號，而復數運算通常比實數運算的計算復雜度高，因此本節欲將復數方程實數化，以減少復雜度。因為數字傅里葉變換結果的上半部分和下半部分互為共軛，并且線性方程經過初等變換后不會改變方程的解，所以，將通過對變換結果的上半部分與下半部分進行相加或相減的方法來完成實數化，如圖5所示。

由此，分別對式（17）中的[Yf]和[Xf]作初等變換，得到：

[Y=Yf1：M2+YfM2+1：MjYf1：M2-YfM2+1：M∈RM×1] （18）

[X=Xf1：M2+XfM2+1：MjXf1：M2-XfM2+1：M∈RM×（R+RK+K）] （19）

式中：[j=-1]；[m：n]表示取矩陣的第[m]～[n]行的元素。將式（17）中的觀測和設計矩陣替換為式（18）、式（19），有：

[Y=Xθ+Ξ] （20）

式中[Ξ]為式（17）初等變換后對應的噪聲矢量，各矩陣和矢量的維度并沒有發生改變，但是其中的元素將不再是復數，而是實數。

1.5" 稀疏變分貝葉斯推斷

線性回歸方程確定后，采用稀疏變分貝葉斯的方法求解模型參數。假設模型參數的變分分布為[q?]，通過對其不斷優化來逼近真實的后驗分布。本文采用如下先驗信息構建模型。

1） 為使所求的模型參數滿足稀疏性，參照文獻[8]的方法引入稀疏約束，用第[i]個連接參數的二進制隨機變量[ζi]作為似然函數中的特征選擇器，其先驗滿足概率為[pi0]的Bernoulli分布。

2） 神經連接參數[θ]的先驗滿足均值為[μ0]、協方差矩陣為[Σ0]的高斯分布。

3） 噪聲[Ξ]的精度參數[τ]的先驗滿足[α0]和[β0]的Gamma分布。其中，需要額外說明的是，伯努利分布概率[pi0]為連接參數存在的先驗概率，當GSD進行模型反演時，通過設置不同的[pi0]值來編碼腦功能網絡的稀疏性。

雖然模型反演在平均場理論[15]下進行，每個腦區的模型參數可獨立求解，但是單個腦區中模型參數的后驗分布不一定有閉合解，因此利用變分分布[q?]對[θ]、[τ]和[ζi]的后驗分布做近似估計。最終，根據變分原理可以推導出模型參數的變分更新方程，最佳的后驗分布就是通過迭代這些更新方程至收斂得到的，其中收斂條件為最大化負（變分）自由能[2]。

限于篇幅，模型參數的最終反演方案可參考文獻[8]。

1.6" 超參數尋優

在變分推斷中，超參數的選擇對性能尤為重要，本節通過代價函數尋優來獲取超參數。令所求的最優超參數[Θ*=s，η， μ0， pi0]，通過求以下表達式得到：

[Θ*=argminΘyg-fXgθΘ2] （21）

式中：[θΘ]表示采用超參數[Θ]獲得的[θ]估計；[f（?）]為反傅里葉數字變換函數。

相比傳統以最大后驗概率和自由能為代價函數的方法，本文的代價函數并未涉及到觀測值[y=yn]，而是使用了GLM的輸出[yg=γn]，這樣可以避免將[yn]中的噪聲引入到代價函數中。另外，代價函數中采用了[yg=γn]，而非其導數[yg=γn]，具體如圖6所示。

根據卷積模型，式（11）可變為：

[γn=UTt?htβ] （22）

由于[Ut]為外源刺激，可視為階躍性信號，與其微分[Ut]具有相同的周期，所以兩者的頻譜在部分頻率上具有相似性，如圖6所示。由式（22）可知，[yg]與[yg]在頻譜上也具有相似性，因此式（21）中的代價函數直接用GLM輸出替代其微分。

2" 實驗設置

2.1" 仿真數據

仿真數據是基于人類大腦的小世界架構而提出的腦功能網絡結構，本文在文獻[12]的研究基礎上，分別構建了S50和S20網絡，如圖7所示。

測試算法中的參數設置如下。

1） 稀疏回歸DCM（Sparse Regression DCM， sparse rDCM）。文獻[8]提出一種可用于分析大型腦網絡的稀疏DCM，[C]矩陣（驅動輸入矩陣）的非0位置處，伯努利先驗概率[pi0]被設置為1，其余為0；[A]矩陣（固有連接矩陣）中各位置對應的[pi0]的數值，除了自連接被設置為1之外，其余均被設置為0.15；初始值[μ0]與[Σ0]使用了SPM8（版本R4290）中DCM10的標準先驗值；頻域濾波的截止頻率由該算法中濾波器自行計算得到。

2） GSD。GSD為本文所提算法，其參數[pi0]、[μ0]、[Σ0]及濾波器參數均與sparse rDCM的設置一致。

3） DCM。DCM為文獻[2]提出的傳統DCM方法，所有參數采用SPM12軟件（修訂版7771）的DCM12.5的設置，其軟件下載地址為https：//www.fil.ion.ucl.ac.uk/spm/。

在本文仿真數據中，信噪比（Signal to Noise Ratio， SNR）被設置為1、3、5、10和100。需要說明的是，由于DCM在模型反演時的計算復雜度過大，如果在實驗平臺（操作系統為Windows 11，處理器為AMD Ryzen 5 5600H with Radeon Graphics，RAM為16.0 GB，仿真軟件為Matlab R2020a）中采用DCM對S50網絡進行反演，會導致內存溢出，因此在S50網絡中只給出了GSD與sparse rDCM的反演結果。在S20網絡中，由于節點數較少，計算量相對較小，因此給出了DCM的分析結果，性能指標包括模型的反演速度、靈敏度、特異性、精確率和估計誤差。

靈敏度被定義為反演得到的真陽性（參數期望為非0，參數反演值也非0，其中若反演值的絕對值大于閾值則認為非0，實驗結果給出了多個閾值的情形，包括取值為[10-5]、[10-3]、[10-2]）的數目與參數期望為陽（參數的期望為非0）的數目之比。特異性被定義為反演得到的真陰性（參數期望為0，參數反演也為0）的數目與參數期望為陰（參數期望為0）的數目之比。精確率被定義為反演得到的真陽性數目與反演得到的陽性數目的比。估計誤差被定義為估計參數值與真實參數值之間的二范數。

值得注意的是，由于GSD與sparse rDCM在模型反演時都引入了稀疏項，用閾值判定反演值是否為非0的操作是反演過程中的一個環節，不同的閾值設置會導致反演結果的不同。

2.2" 實測數據

注意力與運動數據常被用于研究注意力對視覺通路連通性的調節問題，使用DCM可分析這些調節作用，該數據集的獲取網址為：http：//www.fil.ion.ucl.ac.uk/spm/data/attention/，詳情可參考文獻[13?14]。

本文定義了光感、運動、注意力3個回歸因子，從3個感興趣的區域中提取了BOLD信號，分別是初級視覺皮層（V1）、運動敏感區域（V5）和注意力敏感上級頂葉皮層（SPC）。

根據文獻[13]的研究結果，本文測試的連接模型如圖8a）所示，其中固有連接（[A]矩陣）包括3個腦區的自連接，SPC與V5之間的前向連接（SPC到V5）與后向連接（V5到SPC），V1與V5之間的前向連接（V1到V5）與后向連接（V5到V1）；調制輸入（[B]矩陣）包括：運動調制了V1到V5的固有連接，注意力調制了SPC到V5的固有連接；驅動輸入（[C]矩陣）方面，光感刺激驅動V1中的活動。在本文實驗中，GSD與傳統DCM分析的是圖8a）的模型，而sparse rDCM由于無法估計[B]矩陣，因此分析的是圖8b）的模型，兩種模型的唯一區別就是有無調制輸入。

測試的相關算法參數如下。

sparse rDCM：[μ0]和[Σ0]的設置使用了SPM8中DCM10的標準先驗值；濾波器的截止頻率由該算法自行計算得到；[pi0]通過最大化負自由能得到，具體為[C]矩陣的非0位置處，[pi0]被設為0.5，其余為0；[A]矩陣中各位置對應的[pi0]的數值，除了自連接被設置為1之外，其余均被設置為0.5。

GSD：[pi0]的設置由式（21）的尋優得到，具體為[C]矩陣與[B]矩陣的非0位置處的[pi0]被設為0.8，其余為0；[A]矩陣中各位置對應的[pi0]，除了自連接被設置為1之外，其余均被設置為0.8；[Σ0]的設置采用了SPM8中DCM10使用的標準先驗值；[μ0]與頻域濾波的截止頻率，有以下兩種設置。

1） 通過式（21）的尋優得到分別為0.8和60（高通截止頻率）。

2） [μ0]使用SPM8中DCM10的標準神經先驗，截止頻率由sparse rDCM的算法計算所得。

DCM：所有參數采用SPM12中DCM12.5的設置。

在該實測數據的結果中，除了展示連接模型中各連接參數的反演結果，還給出了擬合誤差的結果，即擬合BOLD信號與真實BOLD信號之間的差值。其中擬合BOLD信號被表示為神經連接參數與設計矩陣乘積的反傅里葉數字變換。

3" 結果分析

3.1" 仿真數據

本小節給出的結果均是經過20次重復獨立計算后的均值。圖9～圖11分別給出了GSD與sparse rDCM在閾值為[10-5]、[10-3]、[10-2]時，反演S50網絡的性能結果。在特異性和靈敏度上，兩者的性能較為接近，除了在SNR=1時。然而，在精確率上，GSD算法在所有SNR值下均優于sparse rDCM。特別地，在SNR=100時，sparse rDCM的精確率還有所下降，但GSD的精確率仍在上升。另外，兩種算法隨著閾值從[10-5]增大到[10-2]時，精確率均在提升，并且GSD的精確率仍然高于sparse rDCM。值得注意的是，GSD在計算時間上要遠低于sparse rDCM，約為后者的[13]。

圖12給出了SNR=3下DCM、GSD和sparse rDCM算法對S20網絡的反演結果。由圖12可看到，DCM的精確率和特異性均低于其他兩種算法，并且運算時間也大幅高于其他兩種算法。另外，雖然本文GSD算法的特異性等性能與sparse rDCM算法的性能接近，但計算時間仍然少于后者，該結果與圖9中關于S50網絡的結果一致。

3.2" 實測數據

本節給出的結果仍是經過20次重復獨立計算后的均值，閾值設置為[10-3]。首先，三種模型的反演時間如圖13所示。GSD的反演時間不足1 s，sparse rDCM的反演時間約在1～2 s之間，DCM的反演時間已超過12 s，仍然是GSD所用時間最少。

在固有連接方面，雖然三種模型中的3個腦區均顯示有自連接，但是GSD與DCM中的這些自連接都為激勵效應，而sparse rDCM是抑制效應；在GSD與DCM中發現了SPC與V5之間的前向連接（SPC到V5）與后向連接（V5到SPC），但是sparse rDCM只存在后向連接；DCM和sparse rDCM的連接模型均沒有顯示V1與SPC之間的連接，而GSD中V1到SPC的連接強度只達到[10-2]量級，也可近似地視為沒有連接；V1與V5之間的后向連接（V5到V1）在3種模型中均被發現，而它們之間的前向連接（V1到V5）只在DCM中被發現，但也僅為[10-2]量級。在輸入調制方面，GSD與DCM顯示出了相同的調制效應，sparse rDCM未顯示任何調制效應。在驅動輸入方面，GSD與DCM均有光感刺激輸入在V1區域，sparse rDCM依然未顯示任何效應。從以上結果來看，GSD與DCM中各連接參數的反演結果基本一致，而sparse rDCM的反演結果與前兩者相差較大，且幅值也較小。

最后，圖15給出了2種超參數尋優方法的擬合誤差，其中第1種為采用式（21）的尋優方法，第2種為sparse rDCM的尋優方法。在3個ROI中，第1種方法的擬合誤差均小于第2種方法。

4" 結" 論

從實驗結果看，GSD算法的計算速度有著較大的提升。在仿真數據集中，GSD的反演速度相較于sparse rDCM提高了約70%，相較于傳統DCM提高了近98%。在實測數據集中，GSD算法的速度優勢也得到了延續。

GSD在有著高反演速度的同時，還具備著較強的模型解釋力。比如在實測數據的實驗中，無論是固有連接還是輸入調制，GSD與傳統DCM的反演結果都具有較好的一致性，兩者顯示出了相同的效應（激勵效應或抑制效應）和相似的連接強度。sparse rDCM在該數據集中則完全“失靈”。

GSD算法的超參數設置更具靈活性與科學性。在實測數據中，使用了GSD的尋優函數所得到的結果有著較低的擬合誤差。

注：本文通訊作者為吳海鋒。

參考文獻

[1] 馬園園，鄭罡，周潔敏，等.基于fMRI的腦功能整合數據分析方法綜述[J].生物物理學報，2011，27（1）：18?27.

[2] FRISTON K J， HARRISON L， PENNY W. Dynamic causal modeling [J]. Neuroimage， 2003， 19（4）： 1273?1302.

[3] LIANG Q J， LI J H， ZHENG S N， et al. Dynamic causal modelling of hierarchical planning [J]. Neuroimage， 2022， 258： 119384.

[4] OU Y L， DAI P S， ZHOU X Y， et al. A strategy of model space search for dynamic causal modeling in task fMRI data exploratory analysis [J]. Physical and engineering sciences in medicine， 2022， 45（3）： 867?882.

[5] SHI Y， LI Y. The effective connectivity analysis of fMRI based on asymmetric detection of transfer brain entropy [J]. Cerebral cortex， 2024， 34（3）： bhae070.

[6] FR?SSLE S， LOMAKINA E I， RAZI A， et al. Regression DCM for fMRI [J]. Neuroimage， 2017， 155： 406?421.

[7] FRISTON K J， KAHAN J， BISWAL B B， et al. A DCM for resting state fMRI [J]. Neuroimage， 2014， 94： 396?407.

[8] FR?SSLE S， LOMAKINA E I， KASPER L， et al. A generative model of whole?brain effective connectivity [J]. Neuroimage， 2018， 179： 505?529.

[9] SAHOO D， HONNORAT N， DAVATZIKOS C. Sparse low?dimensional causal modeling for the analysis of brain function [C]// Medical Imaging 2019： Image Processing. Bellingham， WA： SPIE， 2019： 701?707.

[10] WORSLEY K J， FRISTON K J. Analysis of fMRI time?series revisited： again [J]. Neuroimage， 1995， 2（3）： 173?181.

[11] POLINE J B， BRETT M. The general linear model and fMRI： does love last forever？ [J]. Neuroimage， 2012， 62（2）： 871?880.

[12] SMITH S M， MILLER K L， SALIMI?KHORSHIDI G， et al. Network modelling methods for FMRI [J]. Neuroimage， 2011， 54（2）： 875?891.

[13] BüCHEL C， FRISTON K J. Modulation of connectivity in visual pathways by attention： cortical interactions evaluated with structural equation modelling and fMRI [J]. Cerebral cortex， 1997， 7（8）： 768?778.

[14] FRISTON K J， FRITH C D， TURNER R， et al. Characterizing evoked hemodynamics with fMRI [J]. Neuroimage， 1995， 2（2）： 157?165.

[15] H?JEN?S?RENSEN P A， WINTHER O， HANSEN L K. Mean?field approaches to independent component analysis [J]. Neural computation， 2002， 14（4）： 889?918.

作者簡介：胡新航（1998—），男，山東菏澤人，碩士研究生，主要研究方向為生物醫學信號處理。

吳海鋒（1977—），男，云南昆明人，博士研究生，教授，碩士生導師，主要研究方向為生物醫學信號處理、機器學習。

曾" 玉（1981—），女，云南昆明人，博士研究生，講師，主要研究方向為生物醫學信號處理、機器學習。

收稿日期：2024?07?16" " " " " "修回日期：2024?08?13

基金項目：國家自然科學基金資助項目（62161052）；云南省教育廳科學研究基金項目（2024Y432）