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

關(guān)鍵詞： 廣義線性模型； 動(dòng)態(tài)因果模型； 稀疏算法； 有效連接； 計(jì)算復(fù)雜度； 模型解釋力

中圖分類號(hào)： TN919?34； TP301.6" " " " " " " " " "文獻(xiàn)標(biāo)識(shí)碼： A" " " " " " " " " " 文章編號(hào)： 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" 引" 言

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

為了提升DCM的運(yùn)算速度并降低運(yùn)算成本，不斷有新變體被提出。一類變體是對(duì)模型進(jìn)行簡(jiǎn)化[6?7]，例如，簡(jiǎn)化神經(jīng)動(dòng)力學(xué)模型，將模型中的某些生物學(xué)過(guò)程或連接考慮為常數(shù)或固定值，從而減少了需要估計(jì)的參數(shù)數(shù)量；再如簡(jiǎn)化模型本身結(jié)構(gòu)，即通過(guò)限制模型中的連接類型或拓?fù)浣Y(jié)構(gòu)來(lái)減少待估計(jì)的參數(shù)數(shù)量。盡管這些簡(jiǎn)化模型降低了計(jì)算復(fù)雜度，但也帶來(lái)了一些缺點(diǎn)。例如，簡(jiǎn)化模型可能會(huì)喪失一些生物學(xué)上的細(xì)節(jié)，導(dǎo)致模型的解釋力和預(yù)測(cè)能力降低。另一類變體是通過(guò)改進(jìn)參數(shù)估計(jì)和優(yōu)化計(jì)算方法來(lái)提高DCM計(jì)算速度[8?9]，例如，通過(guò)近似推斷和變換域方法等。盡管這些方法在提高計(jì)算速度上取得了一定進(jìn)展，但它們往往需要對(duì)模型做出一些假設(shè)，從而可能影響到模型的準(zhǔn)確性和可靠性。

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

本文的貢獻(xiàn)如下：

1） 利用傅里葉變換的對(duì)稱特性，將頻域DCM的復(fù)數(shù)計(jì)算變?yōu)閷?shí)數(shù)計(jì)算，從而降低了計(jì)算復(fù)雜度；

2） 利用GLM算法和頻域?yàn)V波器，消除了觀測(cè)信號(hào)的部分干擾，提高了模型參數(shù)估計(jì)的準(zhǔn)確度；

3） 利用定義的代價(jià)函數(shù)，確定了變分推斷參數(shù)及濾波器參數(shù)，從而進(jìn)一步提高參數(shù)估計(jì)準(zhǔn)確度；

4） 使用兩組公開數(shù)據(jù)[12?14]對(duì)GSD算法進(jìn)行驗(yàn)證。

1" 算" 法

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

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

2） 為了進(jìn)一步降低觀測(cè)信號(hào)中的噪聲干擾，引入了兩種濾波技術(shù)。第一種濾波技術(shù)是通過(guò)使用GLM估計(jì)得到的擬合系數(shù)來(lái)重構(gòu)GLM信號(hào)，并將該信號(hào)替代原始BOLD信號(hào)于線性回歸方程中，實(shí)現(xiàn)第一次濾波。緊接著，采用頻域?yàn)V波器對(duì)方程進(jìn)行第二次濾波處理。經(jīng)過(guò)兩輪濾波處理的線性回歸方程盡管仍然呈現(xiàn)復(fù)數(shù)形式，但可以通過(guò)初等變換的方法將其轉(zhuǎn)換為實(shí)數(shù)方程，這一步驟進(jìn)一步簡(jiǎn)化了計(jì)算過(guò)程。

3） 通過(guò)對(duì)代價(jià)函數(shù)的尋優(yōu)確定出最佳的超參數(shù)，并通過(guò)變分推斷技術(shù)得到最終的DCM估計(jì)參數(shù)。

1.1" DCM的線性回歸

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

令[z（r）t]為第[r]個(gè)腦區(qū)在[t]時(shí)刻的神經(jīng)元活動(dòng)狀態(tài)，則根據(jù)動(dòng)態(tài)理論及DCM[2]，該腦區(qū)神經(jīng)元狀態(tài)的微分[zt]（為簡(jiǎn)化，略去[r]）可表示為：

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

式中：[ujt]代表[t]時(shí)刻的第[j]個(gè)刺激；[Zt=z（r）t∈RR×1]代表由[R]個(gè)腦區(qū)在[t]時(shí)刻的神經(jīng)元活動(dòng)構(gòu)成的列矢量；[Ut=ukt∈RK×1]代表[t]時(shí)刻的[K]個(gè)刺激構(gòu)成的列矢量；[A∈R1×R]代表與刺激無(wú)關(guān)的內(nèi)生網(wǎng)絡(luò)連接系數(shù)行矢量；[Bj∈R1×R]代表受第[j]個(gè)刺激調(diào)節(jié)的內(nèi)生網(wǎng)絡(luò)連接系數(shù)行矢量；[C∈R1×K]代表受刺激直接影響的神經(jīng)元活動(dòng)的系數(shù)行矢量。

由線性卷積模型，BOLD信號(hào)可表示為：

[yt=zt?ht]" （2）

式中[ht]為HRF。因此對(duì)式（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得到的擬合系數(shù)矩陣；[y″t=Ut?ht]。

將BOLD信號(hào)[yt]表達(dá)為離散形式，則其在第[n]個(gè)采樣點(diǎn)上的值為：

[yn=ynΔt]" （4）

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

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

將式（5）代入式（3），再對(duì)[n=1，2，…，N]個(gè)采樣點(diǎn)做數(shù)字傅里葉變換，則有：

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

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

[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）變?yōu)椋?/p>

[Yo=Xoθ] （10）

式中：向量[Yo]由BOLD信號(hào)的導(dǎo)數(shù)構(gòu)成，為已知；設(shè)計(jì)矩陣[Xo]由BOLD信號(hào)、刺激及GLM擬合系數(shù)構(gòu)成，也為已知；參數(shù)矢量[θ]是待估計(jì)的神經(jīng)連接參數(shù)，為未知。因此式（10）就可以看為一個(gè)頻域中的線性回歸方程。

1.2" GLM濾波

為提高估計(jì)準(zhǔn)確度，本節(jié)將采用GLM[10]濾波和頻域?yàn)V波兩種技術(shù)來(lái)降噪，如圖3所示。在GLM濾波中，先利用GLM重構(gòu)BOLD信號(hào)，將式（10）中的觀測(cè)和設(shè)計(jì)矩陣中的原BOLD信號(hào)替換為重構(gòu)的GLM信號(hào)，得到新的觀測(cè)和設(shè)計(jì)矩陣，如圖3中的虛線框所示。在頻域?yàn)V波中，對(duì)通過(guò)GLM濾波后得到的觀測(cè)和設(shè)計(jì)矩陣再進(jìn)行理想濾波，如圖3中實(shí)線框所示。

一個(gè)腦區(qū)的GLM輸出信號(hào)[yg=γn∈RN×1]可表示為：

[yg=xgβ] （11）

式中：[xg∈RN×R]為GLM設(shè)計(jì)矩陣；[β∈RR×1]為擬合系數(shù)矩陣。將該腦區(qū)的觀測(cè)BOLD信號(hào)[yn]用GLM輸出[γn]來(lái)替代，有：

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

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

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

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

相應(yīng)的設(shè)計(jì)矩陣[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信號(hào)[yn]的傅里葉變換替換為GLM輸出[γn]的傅里葉變換后的結(jié)果。

1.3" 頻域?yàn)V波

本節(jié)將詳細(xì)描述圖3所示的頻域?yàn)V波。

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

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

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

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

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

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

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

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

在式（17）中，濾波后又得到了新的觀測(cè)和設(shè)計(jì)矩陣。需要注意的是，濾波器設(shè)計(jì)中的參數(shù)[s]和[η]其實(shí)分別是濾波器低頻和高頻的截止頻率，可通過(guò)對(duì)代價(jià)函數(shù)尋優(yōu)得到，這將在后面的小節(jié)中描述。

1.4" 實(shí)數(shù)化

由于本文將DCM放在頻域中進(jìn)行求解，則原來(lái)的實(shí)信號(hào)就變?yōu)閺?fù)信號(hào)，而復(fù)數(shù)運(yùn)算通常比實(shí)數(shù)運(yùn)算的計(jì)算復(fù)雜度高，因此本節(jié)欲將復(fù)數(shù)方程實(shí)數(shù)化，以減少?gòu)?fù)雜度。因?yàn)閿?shù)字傅里葉變換結(jié)果的上半部分和下半部分互為共軛，并且線性方程經(jīng)過(guò)初等變換后不會(huì)改變方程的解，所以，將通過(guò)對(duì)變換結(jié)果的上半部分與下半部分進(jìn)行相加或相減的方法來(lái)完成實(shí)數(shù)化，如圖5所示。

由此，分別對(duì)式（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）中的觀測(cè)和設(shè)計(jì)矩陣替換為式（18）、式（19），有：

[Y=Xθ+Ξ] （20）

式中[Ξ]為式（17）初等變換后對(duì)應(yīng)的噪聲矢量，各矩陣和矢量的維度并沒(méi)有發(fā)生改變，但是其中的元素將不再是復(fù)數(shù)，而是實(shí)數(shù)。

1.5" 稀疏變分貝葉斯推斷

線性回歸方程確定后，采用稀疏變分貝葉斯的方法求解模型參數(shù)。假設(shè)模型參數(shù)的變分分布為[q?]，通過(guò)對(duì)其不斷優(yōu)化來(lái)逼近真實(shí)的后驗(yàn)分布。本文采用如下先驗(yàn)信息構(gòu)建模型。

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

2） 神經(jīng)連接參數(shù)[θ]的先驗(yàn)滿足均值為[μ0]、協(xié)方差矩陣為[Σ0]的高斯分布。

3） 噪聲[Ξ]的精度參數(shù)[τ]的先驗(yàn)滿足[α0]和[β0]的Gamma分布。其中，需要額外說(shuō)明的是，伯努利分布概率[pi0]為連接參數(shù)存在的先驗(yàn)概率，當(dāng)GSD進(jìn)行模型反演時(shí)，通過(guò)設(shè)置不同的[pi0]值來(lái)編碼腦功能網(wǎng)絡(luò)的稀疏性。

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

限于篇幅，模型參數(shù)的最終反演方案可參考文獻(xiàn)[8]。

1.6" 超參數(shù)尋優(yōu)

在變分推斷中，超參數(shù)的選擇對(duì)性能尤為重要，本節(jié)通過(guò)代價(jià)函數(shù)尋優(yōu)來(lái)獲取超參數(shù)。令所求的最優(yōu)超參數(shù)[Θ*=s，η， μ0， pi0]，通過(guò)求以下表達(dá)式得到：

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

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

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

根據(jù)卷積模型，式（11）可變?yōu)椋?/p>

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

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

2" 實(shí)驗(yàn)設(shè)置

2.1" 仿真數(shù)據(jù)

仿真數(shù)據(jù)是基于人類大腦的小世界架構(gòu)而提出的腦功能網(wǎng)絡(luò)結(jié)構(gòu)，本文在文獻(xiàn)[12]的研究基礎(chǔ)上，分別構(gòu)建了S50和S20網(wǎng)絡(luò)，如圖7所示。

測(cè)試算法中的參數(shù)設(shè)置如下。

1） 稀疏回歸DCM（Sparse Regression DCM， sparse rDCM）。文獻(xiàn)[8]提出一種可用于分析大型腦網(wǎng)絡(luò)的稀疏DCM，[C]矩陣（驅(qū)動(dòng)輸入矩陣）的非0位置處，伯努利先驗(yàn)概率[pi0]被設(shè)置為1，其余為0；[A]矩陣（固有連接矩陣）中各位置對(duì)應(yīng)的[pi0]的數(shù)值，除了自連接被設(shè)置為1之外，其余均被設(shè)置為0.15；初始值[μ0]與[Σ0]使用了SPM8（版本R4290）中DCM10的標(biāo)準(zhǔn)先驗(yàn)值；頻域?yàn)V波的截止頻率由該算法中濾波器自行計(jì)算得到。

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

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

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

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

值得注意的是，由于GSD與sparse rDCM在模型反演時(shí)都引入了稀疏項(xiàng)，用閾值判定反演值是否為非0的操作是反演過(guò)程中的一個(gè)環(huán)節(jié)，不同的閾值設(shè)置會(huì)導(dǎo)致反演結(jié)果的不同。

2.2" 實(shí)測(cè)數(shù)據(jù)

注意力與運(yùn)動(dòng)數(shù)據(jù)常被用于研究注意力對(duì)視覺(jué)通路連通性的調(diào)節(jié)問(wèn)題，使用DCM可分析這些調(diào)節(jié)作用，該數(shù)據(jù)集的獲取網(wǎng)址為：http：//www.fil.ion.ucl.ac.uk/spm/data/attention/，詳情可參考文獻(xiàn)[13?14]。

本文定義了光感、運(yùn)動(dòng)、注意力3個(gè)回歸因子，從3個(gè)感興趣的區(qū)域中提取了BOLD信號(hào)，分別是初級(jí)視覺(jué)皮層（V1）、運(yùn)動(dòng)敏感區(qū)域（V5）和注意力敏感上級(jí)頂葉皮層（SPC）。

根據(jù)文獻(xiàn)[13]的研究結(jié)果，本文測(cè)試的連接模型如圖8a）所示，其中固有連接（[A]矩陣）包括3個(gè)腦區(qū)的自連接，SPC與V5之間的前向連接（SPC到V5）與后向連接（V5到SPC），V1與V5之間的前向連接（V1到V5）與后向連接（V5到V1）；調(diào)制輸入（[B]矩陣）包括：運(yùn)動(dòng)調(diào)制了V1到V5的固有連接，注意力調(diào)制了SPC到V5的固有連接；驅(qū)動(dòng)輸入（[C]矩陣）方面，光感刺激驅(qū)動(dòng)V1中的活動(dòng)。在本文實(shí)驗(yàn)中，GSD與傳統(tǒng)DCM分析的是圖8a）的模型，而sparse rDCM由于無(wú)法估計(jì)[B]矩陣，因此分析的是圖8b）的模型，兩種模型的唯一區(qū)別就是有無(wú)調(diào)制輸入。

測(cè)試的相關(guān)算法參數(shù)如下。

sparse rDCM：[μ0]和[Σ0]的設(shè)置使用了SPM8中DCM10的標(biāo)準(zhǔn)先驗(yàn)值；濾波器的截止頻率由該算法自行計(jì)算得到；[pi0]通過(guò)最大化負(fù)自由能得到，具體為[C]矩陣的非0位置處，[pi0]被設(shè)為0.5，其余為0；[A]矩陣中各位置對(duì)應(yīng)的[pi0]的數(shù)值，除了自連接被設(shè)置為1之外，其余均被設(shè)置為0.5。

GSD：[pi0]的設(shè)置由式（21）的尋優(yōu)得到，具體為[C]矩陣與[B]矩陣的非0位置處的[pi0]被設(shè)為0.8，其余為0；[A]矩陣中各位置對(duì)應(yīng)的[pi0]，除了自連接被設(shè)置為1之外，其余均被設(shè)置為0.8；[Σ0]的設(shè)置采用了SPM8中DCM10使用的標(biāo)準(zhǔn)先驗(yàn)值；[μ0]與頻域?yàn)V波的截止頻率，有以下兩種設(shè)置。

1） 通過(guò)式（21）的尋優(yōu)得到分別為0.8和60（高通截止頻率）。

2） [μ0]使用SPM8中DCM10的標(biāo)準(zhǔn)神經(jīng)先驗(yàn)，截止頻率由sparse rDCM的算法計(jì)算所得。

DCM：所有參數(shù)采用SPM12中DCM12.5的設(shè)置。

在該實(shí)測(cè)數(shù)據(jù)的結(jié)果中，除了展示連接模型中各連接參數(shù)的反演結(jié)果，還給出了擬合誤差的結(jié)果，即擬合BOLD信號(hào)與真實(shí)BOLD信號(hào)之間的差值。其中擬合BOLD信號(hào)被表示為神經(jīng)連接參數(shù)與設(shè)計(jì)矩陣乘積的反傅里葉數(shù)字變換。

3" 結(jié)果分析

3.1" 仿真數(shù)據(jù)

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

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

3.2" 實(shí)測(cè)數(shù)據(jù)

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

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

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

4" 結(jié)" 論

從實(shí)驗(yàn)結(jié)果看，GSD算法的計(jì)算速度有著較大的提升。在仿真數(shù)據(jù)集中，GSD的反演速度相較于sparse rDCM提高了約70%，相較于傳統(tǒng)DCM提高了近98%。在實(shí)測(cè)數(shù)據(jù)集中，GSD算法的速度優(yōu)勢(shì)也得到了延續(xù)。

GSD在有著高反演速度的同時(shí)，還具備著較強(qiáng)的模型解釋力。比如在實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)中，無(wú)論是固有連接還是輸入調(diào)制，GSD與傳統(tǒng)DCM的反演結(jié)果都具有較好的一致性，兩者顯示出了相同的效應(yīng)（激勵(lì)效應(yīng)或抑制效應(yīng)）和相似的連接強(qiáng)度。sparse rDCM在該數(shù)據(jù)集中則完全“失靈”。

GSD算法的超參數(shù)設(shè)置更具靈活性與科學(xué)性。在實(shí)測(cè)數(shù)據(jù)中，使用了GSD的尋優(yōu)函數(shù)所得到的結(jié)果有著較低的擬合誤差。

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

參考文獻(xiàn)

[1] 馬園園，鄭罡，周潔敏，等.基于fMRI的腦功能整合數(shù)據(jù)分析方法綜述[J].生物物理學(xué)報(bào)，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.

作者簡(jiǎn)介：胡新航（1998—），男，山東菏澤人，碩士研究生，主要研究方向?yàn)樯镝t(yī)學(xué)信號(hào)處理。

吳海鋒（1977—），男，云南昆明人，博士研究生，教授，碩士生導(dǎo)師，主要研究方向?yàn)樯镝t(yī)學(xué)信號(hào)處理、機(jī)器學(xué)習(xí)。

曾" 玉（1981—），女，云南昆明人，博士研究生，講師，主要研究方向?yàn)樯镝t(yī)學(xué)信號(hào)處理、機(jī)器學(xué)習(xí)。

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（62161052）；云南省教育廳科學(xué)研究基金項(xiàng)目（2024Y432）