本文的目標是處理并分析使用深度電極在難治性癲癇患者癲癇發作期間其大腦皮層中記錄到的癲癇腦電信號間的大腦效應連通性。維納-格蘭杰因果索引算法是一種眾所周知的檢測腦電信號間大腦效應連通性的有效方法。它是一種基于線性自回歸模型的方法，而模型參數估計問題在其用于腦電因果效應連通性研究中的計算準確性與魯棒性方面起著至關重要的作用。本文針對這一問題，使用了我們提出的改進的赤池信息量準則來估計算法中自回歸模型的模型階數，以提高維納-格蘭杰因果索引算法檢測大腦效應連通性的性能。實驗仿真結果表明：不管是在線性隨機系統中還是在能生成模擬癲癇信號的生理模型中，該改進的維納-格蘭杰因果索引算法在檢測腦效應連通性上都表現出良好的魯棒性。
引用本文： 楊淳沨, 向文濤, 伍家松, 孔佑勇, 姜龍玉, Le BouquinJèannes Régine, 舒華忠. 基于通用赤池信息量準則改進維納-格蘭杰因果索引算法的顱內腦電效應連通性研究. 生物醫學工程學雜志, 2018, 35(5): 665-671. doi: 10.7507/1001-5515.201709032 復制
引言
癲癇是嚴重影響人類日常生活的腦部疾病之一。其中難治性癲癇綜合征的藥物治療效果差，嚴重影響了患者的生活質量，解決手段之一就是通過手術來切除致癇區（epileptogenic zone，EZ），從而減少及消除癲癇發作[1]。在癲癇手術中，其目標是精確定位致癇區和相應的腦網絡，并由此來判斷切除手術是否會影響到正常的腦功能。因此，問題的關鍵就是確定致癇區的組成及手術要移除的區域，在抑制或者減少癲癇發作的同時有效減少認知、感覺及運動損傷。利用深度電極記錄顱內腦電信號（intracerebral electroencephalogram，iEEG）是一種術前診斷方法。其中的困難之處在于通過分析這些腦電信號之間的關系來確定哪些腦的部位屬于致癇區。關于這點，采用信號處理的方法可以提供僅靠單純的視覺分析無法獲取的一些定量信息，如已被證實有助于評估遠距離腦區間的功能耦合的相關性方法[2]。1956 年，Wiener[3]首先提出了時間序列之間因果關系的概念，認為如果憑借第一個時間序列的知識可以預測第二個時間序列，那么第一個序列對第二個序列有因果影響。這個概念由 Granger[4]在 1969 年進行了補充說明，本文中稱為 Wiener-Granger 因果索引（Wiener-Granger Causality Index，WGCI）。自 20 世紀 60 年代以來，從這個概念引申出的 Wiener-Granger 因果索引已在經濟學中廣泛使用，近年來在神經科學中也得到了廣泛應用[5-7]。WGCI 的基本思想如下：如果一個過程s1“導致”另一過程s2 的發生，則s1 過去的值包含了幫助預測s2 的信息，并且此信息大于單獨從s2 過去的值所包含的信息。WGCI 是一種基于線性自回歸模型的檢測時間序列之間因果關系的有效方法，故在該算法中，其所用的模型階數估計方法起著非常關鍵的作用，其中一種經典方法是赤池信息準則（Akaike’s information criterion，AIC）[8]。
在經典的 AIC 算法中，會假設不同信號具有相同階數，然而在實際情況中，不同信號各自的階數往往是不相同的，所以 AIC 算法會使得模型的項數產生過估計，導致后續的系數估計運算量加大且精準度不高。針對這些問題，本文采用了一種我們提出的改進的 AIC 準則，稱為通用 AIC（generalized AIC，gAIC）算法[9]，用于估計模型中的階數并將之用于 WGCI 算法。首先，我們使用 gAIC 作為在 WGCI 算法估計中起重要作用的模型階數估計方法，先在模擬信號上驗證了結合 gAIC 的 WGCI 算法的魯棒性并與經典的 AIC 準則進行了比較；然后將其應用于一個生理模型，該生理模型可模擬兩個神經元集群及其相互之間的因果聯系。這種基于生理學的模型有如下重要特征：一方面，由該生理模型生成的模擬 EEG 信號，具有癲癇患者發作時快速發作活動階段（fast onset activity，FOA）記錄所得真實癲癇 EEG 信號相似的特征；另一方面，相對于真實癲癇 EEG 信號，這種模型能設定兩個神經元集群（模擬癲癇患者腦部的兩個不同區域）之間相互聯系的強弱及傳播的方向，因此，本研究采用該模型來驗證算法在真實信號上的魯棒性。
1 方法
Wiener-Granger 因果索引算法是一種基于自回歸模型的檢測腦電信號間效應連通性的有效方法。而在該算法中，自回歸模型的模型階數估計起著至關重要的作用，在本文中，采用了我們提出的 gAIC 算法來估計自回歸模型的模型階數，減少了自回歸模型系數的估計量并提高了 Wiener-Granger 因果索引算法的魯棒性。
1.1 Wiener-Granger 因果索引算法
給出兩個均值為零的信號
和
，相應的時間觀測序列變量為
和
，
，
是序列長度。假設每個觀測序列由如下模型階數為
的一維自回歸模型分別建模，則有：
 |
其中，
和
為高斯白噪聲（
和
是相應的白過程），由模型可知每個信號只依賴于它自身過去的值。若利用以下的二維自回歸模型對兩個信號同時建模，則有：
 |
此時，每個信號不僅依賴于它本身過去的值，而且還依賴于另一個信號的過去的值。變量
和
也是高斯白噪聲實現（
和
是相應的白過程），回歸模型系數
（
）由最小二乘法[10]可得，模型的階數
由 AIC 決定。
為了研究信號
和
之間的因果關系，下面以從
到
的因果關系為例：在上述式（1）中的單變量模型中，
表示的好壞由預測誤差變量
來估計，這里的
代表
的過去的值，在上述式（2）中的二變量模型中，
則代表了相應的預測誤差變量。如果
“Wiener-Granger”導致
的發生，則
大于
。因此，從
到
的 WGCI 定義如下：
 |
相應地，從
到
的 WGCI 表示為：
 |
1.2 回歸模型階數估計方法
考慮一個零均值m維的向量自回歸（vectorial autoregressive，VAR）過程
：
 |
其中，
為
的系數矩陣，
是向量自回歸模型的階數。獨立同分布過程
是協方差矩陣為Σ的零均值向量過程，對于每個變量
，
和
都是獨立的。
利用式（5）中的 VAR 過程來產生 N 個時間觀測序列
，這些觀測序列可利用最小二乘法通過
階向量自回歸模型來擬合，相應的估計模型可由以下等式給出：
 |
對應式（6）中的系數矩陣
可通過如下方程的解來求得：
 |
其中，
，
,
是
的第
行。式（7）中最小二乘法的解為：
 |
然后，向量自回歸模型的協方差矩陣計算如下：
 |
其中
。
為了選出階參數
，所用 AIC 準則的定義如下：
 |
將式（5）用二維（
）通用向量自回歸過程
來重寫，則有：
 |
其中，
為向量自回歸過程系數。通常情況下，對于式（11）中的四個階數 AIC 僅能夠給出相同的值，即
。然而，對于大多數的向量自回歸過程來說，四個階數往往是不相同的，經典的 AIC 準則僅僅使用它們中間的最大值，這就使得在后續的計算中產生過估計問題。針對這個問題，我們提出如下的針對四個不同階數
的改進準則，稱為 gAIC：
 |
其中，
指所需估計階數的集合，即
。階數估計的具體步驟如下：
（1）使用 AIC 方法得到自回歸模型的一個階數估計值，記為
，選取這個值確定 gAIC 方法所要估計的四個模型階數集合的選擇范圍，即
,
；
（2）估計第一個信號的兩個模型階數
和
，在估計過程中，將第二個信號的兩個模型階數
和
的值保持不變，設為
。然后，執行從式（7）到式（9）中的計算過程并通過式（12）的定義，令
遍歷步驟（1）中階數選取范圍中的所有可能值，當
取最小值時對應的模型階數即可得到第一個信號的最佳模型階數估計值
和
；
（3）估計第二個信號的兩個模型階數
和
，在估計過程中，將第一個信號的兩個模型階數
和
的值設為步驟（2）中取得的最佳模型階數估計值
和
。然后，執行從式（7）到式（9）中的計算過程并通過式（12）的定義，令
遍歷步驟（1）中階數選取范圍中的所有可能值，當gAIC（
）取最小值時對應的模型階數即可得到第二個信號的最佳模型階數估計值
和
。
1.3 顱內腦電信號生成模型
本文采用一個基于生理學的模型來模擬大腦海馬區域不同神經元集群（存在相互聯系）間的場電位活動。每個記錄到的神經元集群產生的局部活動可以認為是一個iEEG 信號。在這個模型中，每個神經元集群包含三個子神經元集群，這些子神經元集群分別由兩類局部抑制中間神經元和主錐體細胞組成，如圖 1 所示，其可用一個隨機偏微分方程組[見公式（13）]來表示并模擬iEEG 信號。讀者可以閱讀文獻[11-12]進一步了解更多細節。因為主錐體細胞是一種將其軸突投射到其他腦區域的興奮性神經元細胞，基于此，該生理模型將一個神經元集群（m）中主錐體細胞的動作場電位平均脈沖密度作為另一神經元集群（n）的興奮性輸入。此外，用另一個參數Kmn來表征神經元集群m到神經元集群n的關聯性耦合強度。然后，可以通過設置合適的Kmn參數值來構建不同神經元集群間雙向或者單向的聯系。該生理模型中其他主要參數還包含反饋回路中的抑制性和興奮性增益以及子神經元集群中突觸的平均數等，通過調整這些參數來模擬每個神經元集群的不同活動（癲癇發作時和正常情況下的腦部活動）。
圖1
大腦海馬區域不同神經元集群間相互關系。
Figure1.
Interactions between different neuronal subpopulations of the hippocampus
該模型可用下列隨機偏分方程組表示為:
 |
其中，連通常量
，
，突觸時間常量
、
和
以及輸入高斯隨機過程的參數
和
認為是已知的，其值可見參考文獻[11-12]。通過改變模型中參數
、
和
的值來模擬輸出的iEEG 信號為正常腦電信號或者顱內癲癇腦電信號。
2 實驗結果
首先在一個線性隨機系統中驗證和比較 AIC 和 gAIC 兩個準則在 WGCI 算法中的性能，然后將它們用于生理模型所模擬生成的癲癇發作時的iEEG 信號，比較和驗證其性能。
2.1 線性模型
考慮如下的線性隨機系統模型，該模型模擬生成從
到
單向聯系的兩個信號：
 |
其中
和
為高斯白噪聲。我們對 16 個不同長度（長度分別為
個點，
）的信號各自都做了 200 次實驗。每次實驗中，分別用 AIC 和 gAIC 準則估計模型階數并將之用于計算上述模型中兩個信號間的 Wiener-Granger 因果關系。正如 1.2 節中所指出的，我們可以知道式（14）中線性模型的真實模型階數為
。首先，我們給出了當信號長度不同的情況下，使用 gAIC 準則得到模型階數的正確率（正確識別指估計得到的 4 個模型階數與真實模型階數完全一致），其結果如圖 2 所示。在圖中，X軸表示信號的 16 個不同長度，即從 128 個點到 2 048 個點，步長為 128 個點，Y軸表示相應的模型階數估計的正確率。從圖 2 中我們可以得到如下結論：當信號長度超過 512 個點時，gAIC 準則已經幾乎可以正確估計出各個不同的模型階數。以信號長度為 2 048 點為例，使用 gAIC 準則得到的四個模型階數分別為｛2,0,1,4｝,即用于 WGCI 估計時所需估計模型系數的項數為 7 項，而使用 AIC 時項數為 4*4 = 16 項，顯然 WGCI 方法中相比于 AIC 準則，采用 gAIC 準則時運算量降低了一半以上。
圖2
由 gAIC 算法估計 4 個模型階數的正確率
The
其次，我們使用 AIC 或者 gAIC 來估計 WGCI，結果如圖 3 所示，將 WGCI 的平均值和標準差作為考察 AIC 或者 gAIC 準則估計參數的評價指標，WGCI 的值越大，標準差越小，表明所采用的準則越有效。圖中的上半部分，對 16 個長度不同的信號，我們畫出了分別使用 AIC 和 gAIC 方法的從信號
到信號
的 Wiener-Granger 因果索引
；圖中的下半部分，則是相應的 16 個長度不同信號的從信號
到信號
的 Wiener-Granger 因果索引
。從該圖可以看出，無論是使用 AIC 準則還是 gAIC 準則，當信號長度逐漸增加時，所有
的平均值都趨近于同一極限（接近 0.713 9），其相應的標準差則趨近于極限值 0.025 8。由模型可知，該系統中兩個信號間的聯系為從信號 1 到信號 2 的單向聯系，故
的值理論上應為 0。由圖 3 中下半部分可知，使用 gAIC 準則的
的均值總是低于使用 AIC 準則得到的
的平均值，所有使用 gAIC 準則的
的均值都快速減少到一個極小值（約為 0.000 4），甚至當信號長度為 256 點時。實際上，當信號長度大于 1 664 點時，使用 gAIC 準則得到的
平均值就已經是零了。對于從
到
的 WGCI 相應的標準差也有上述相同的結論。因此，綜上所述，我們得出如下的結論：基于 gAIC 準則的 WGCI 方法要比基于 AIC 準則的具有更好的魯棒性。
圖3
在線性隨機系統中使用 AIC 和 gAIC 的 WGCI（平均值和標準差）
The
2.2 生理模型
使用 1.3 節中描述的生理模型來模擬兩個神經元集群間具有固定連通模式的癲癇iEEG 信號，每次模擬生成一組長度為 400 s 的信號，采樣頻率為 256 Hz，如圖 4 所示。兩個神經元集群間的連通模式由K值來表示，本實驗中測試了不同的
值（
），而令
，即該兩個神經元集群間的連通模式為單向聯系。對于每次不同K值（例如，每對
（
）），都進行了 250 次實驗，每次實驗的信號長度為 2 048 點。在本節中，我們也比較了 AIC 和 gAIC 兩種準則的性能，將它們分別用于模型階數的估計，圖 5 給出了相應 WGCI 的平均值和標準差。圖中的上半部分展示了 3 個不同
值（
被設為 0）的從
到
的 WGCI 值，同樣地，圖中的下半部分展示了對應的從
到
的 WGCI 值。乍看起來，不管是使用了 AIC 準則還是 gAIC 準則的算法都表明了信號
和
間的連通模式為從
到
的單向聯系。從該圖的上半部分，我們可以看到，當這兩個信號間單向的因果聯系被檢測到時，WGC 索引值隨著耦合強度（即
）的增強而增長。此外，不管是在哪個耦合強度（不同
值）下，使用 gAIC 準則得到的 WGCI 的平均值總是大于使用 AIC 準則得到的 WGCI 平均值。從圖的下半部分可以看到，無論使用哪種準則，從
到
的 WGCI 的平均值都非常低，但是對于相應的標準差，使用 gAIC 準則得到的要更小一些。因此，總體上來說，我們得出如下結論：在用生理模型生成的模擬癲癇iEEG 信號上，使用 gAIC 模型階數估計準則的 WGCI 算法也表現出比 AIC 準則更好的性能。
圖4
由基于生理模型模擬出的iEEG信號例子
左半部分：兩個模擬
left panel: two simulated
圖5
基于生理學模型使用 AIC 和 gAIC 得到的 WGCI（平均值和標準差）
The
3 結論
本文中，使用了一種名為 gAIC 準則的模型階數估計算法，并將它用于 WGCI 方法來檢測信號間的因果關系。實驗結果表明，相比于基于傳統的 AIC 準則的 WGCI 算法，無論是通過一個線性隨機系統生成的信號的實驗還是用生理模型生成的模擬癲癇信號的實驗，結合 gAIC 準則的 WGCI 算法都體現出了更好的魯棒性，可以很好地識別出信號間的因果關系。在將來，計劃進一步研究該算法在真實癲癇腦電信號中的性能，并同時將該算法拓展到多維信號從而檢測信號間直接或者間接的因果關系。
引言
癲癇是嚴重影響人類日常生活的腦部疾病之一。其中難治性癲癇綜合征的藥物治療效果差，嚴重影響了患者的生活質量，解決手段之一就是通過手術來切除致癇區（epileptogenic zone，EZ），從而減少及消除癲癇發作[1]。在癲癇手術中，其目標是精確定位致癇區和相應的腦網絡，并由此來判斷切除手術是否會影響到正常的腦功能。因此，問題的關鍵就是確定致癇區的組成及手術要移除的區域，在抑制或者減少癲癇發作的同時有效減少認知、感覺及運動損傷。利用深度電極記錄顱內腦電信號（intracerebral electroencephalogram，iEEG）是一種術前診斷方法。其中的困難之處在于通過分析這些腦電信號之間的關系來確定哪些腦的部位屬于致癇區。關于這點，采用信號處理的方法可以提供僅靠單純的視覺分析無法獲取的一些定量信息，如已被證實有助于評估遠距離腦區間的功能耦合的相關性方法[2]。1956 年，Wiener[3]首先提出了時間序列之間因果關系的概念，認為如果憑借第一個時間序列的知識可以預測第二個時間序列，那么第一個序列對第二個序列有因果影響。這個概念由 Granger[4]在 1969 年進行了補充說明，本文中稱為 Wiener-Granger 因果索引（Wiener-Granger Causality Index，WGCI）。自 20 世紀 60 年代以來，從這個概念引申出的 Wiener-Granger 因果索引已在經濟學中廣泛使用，近年來在神經科學中也得到了廣泛應用[5-7]。WGCI 的基本思想如下：如果一個過程s1“導致”另一過程s2 的發生，則s1 過去的值包含了幫助預測s2 的信息，并且此信息大于單獨從s2 過去的值所包含的信息。WGCI 是一種基于線性自回歸模型的檢測時間序列之間因果關系的有效方法，故在該算法中，其所用的模型階數估計方法起著非常關鍵的作用，其中一種經典方法是赤池信息準則（Akaike’s information criterion，AIC）[8]。
在經典的 AIC 算法中，會假設不同信號具有相同階數，然而在實際情況中，不同信號各自的階數往往是不相同的，所以 AIC 算法會使得模型的項數產生過估計，導致后續的系數估計運算量加大且精準度不高。針對這些問題，本文采用了一種我們提出的改進的 AIC 準則，稱為通用 AIC（generalized AIC，gAIC）算法[9]，用于估計模型中的階數并將之用于 WGCI 算法。首先，我們使用 gAIC 作為在 WGCI 算法估計中起重要作用的模型階數估計方法，先在模擬信號上驗證了結合 gAIC 的 WGCI 算法的魯棒性并與經典的 AIC 準則進行了比較；然后將其應用于一個生理模型，該生理模型可模擬兩個神經元集群及其相互之間的因果聯系。這種基于生理學的模型有如下重要特征：一方面，由該生理模型生成的模擬 EEG 信號，具有癲癇患者發作時快速發作活動階段（fast onset activity，FOA）記錄所得真實癲癇 EEG 信號相似的特征；另一方面，相對于真實癲癇 EEG 信號，這種模型能設定兩個神經元集群（模擬癲癇患者腦部的兩個不同區域）之間相互聯系的強弱及傳播的方向，因此，本研究采用該模型來驗證算法在真實信號上的魯棒性。
1 方法
Wiener-Granger 因果索引算法是一種基于自回歸模型的檢測腦電信號間效應連通性的有效方法。而在該算法中，自回歸模型的模型階數估計起著至關重要的作用，在本文中，采用了我們提出的 gAIC 算法來估計自回歸模型的模型階數，減少了自回歸模型系數的估計量并提高了 Wiener-Granger 因果索引算法的魯棒性。
1.1 Wiener-Granger 因果索引算法
給出兩個均值為零的信號
和
，相應的時間觀測序列變量為
和
，
，
是序列長度。假設每個觀測序列由如下模型階數為
的一維自回歸模型分別建模，則有：
|
其中，
和
為高斯白噪聲（
和
是相應的白過程），由模型可知每個信號只依賴于它自身過去的值。若利用以下的二維自回歸模型對兩個信號同時建模，則有：
|
此時，每個信號不僅依賴于它本身過去的值，而且還依賴于另一個信號的過去的值。變量
和
也是高斯白噪聲實現（
和
是相應的白過程），回歸模型系數
（
）由最小二乘法[10]可得，模型的階數
由 AIC 決定。
為了研究信號
和
之間的因果關系，下面以從
到
的因果關系為例：在上述式（1）中的單變量模型中，
表示的好壞由預測誤差變量
來估計，這里的
代表
的過去的值，在上述式（2）中的二變量模型中，
則代表了相應的預測誤差變量。如果
“Wiener-Granger”導致
的發生，則
大于
。因此，從
到
的 WGCI 定義如下：
|
相應地，從
到
的 WGCI 表示為：
|
1.2 回歸模型階數估計方法
考慮一個零均值m維的向量自回歸（vectorial autoregressive，VAR）過程
：
|
其中，
為
的系數矩陣，
是向量自回歸模型的階數。獨立同分布過程
是協方差矩陣為Σ的零均值向量過程，對于每個變量
，
和
都是獨立的。
利用式（5）中的 VAR 過程來產生 N 個時間觀測序列
，這些觀測序列可利用最小二乘法通過
階向量自回歸模型來擬合，相應的估計模型可由以下等式給出：
|
對應式（6）中的系數矩陣
可通過如下方程的解來求得：
|
其中，
，
,
是
的第
行。式（7）中最小二乘法的解為：
|
然后，向量自回歸模型的協方差矩陣計算如下：
|
其中
。
為了選出階參數
，所用 AIC 準則的定義如下：
|
將式（5）用二維（
）通用向量自回歸過程
來重寫，則有：
|
其中，
為向量自回歸過程系數。通常情況下，對于式（11）中的四個階數 AIC 僅能夠給出相同的值，即
。然而，對于大多數的向量自回歸過程來說，四個階數往往是不相同的，經典的 AIC 準則僅僅使用它們中間的最大值，這就使得在后續的計算中產生過估計問題。針對這個問題，我們提出如下的針對四個不同階數
的改進準則，稱為 gAIC：
|
其中，
指所需估計階數的集合，即
。階數估計的具體步驟如下：
（1）使用 AIC 方法得到自回歸模型的一個階數估計值，記為
，選取這個值確定 gAIC 方法所要估計的四個模型階數集合的選擇范圍，即
,
；
（2）估計第一個信號的兩個模型階數
和
，在估計過程中，將第二個信號的兩個模型階數
和
的值保持不變，設為
。然后，執行從式（7）到式（9）中的計算過程并通過式（12）的定義，令
遍歷步驟（1）中階數選取范圍中的所有可能值，當
取最小值時對應的模型階數即可得到第一個信號的最佳模型階數估計值
和
；
（3）估計第二個信號的兩個模型階數
和
，在估計過程中，將第一個信號的兩個模型階數
和
的值設為步驟（2）中取得的最佳模型階數估計值
和
。然后，執行從式（7）到式（9）中的計算過程并通過式（12）的定義，令
遍歷步驟（1）中階數選取范圍中的所有可能值，當gAIC（
）取最小值時對應的模型階數即可得到第二個信號的最佳模型階數估計值
和
。
1.3 顱內腦電信號生成模型
本文采用一個基于生理學的模型來模擬大腦海馬區域不同神經元集群（存在相互聯系）間的場電位活動。每個記錄到的神經元集群產生的局部活動可以認為是一個iEEG 信號。在這個模型中，每個神經元集群包含三個子神經元集群，這些子神經元集群分別由兩類局部抑制中間神經元和主錐體細胞組成，如圖 1 所示，其可用一個隨機偏微分方程組[見公式（13）]來表示并模擬iEEG 信號。讀者可以閱讀文獻[11-12]進一步了解更多細節。因為主錐體細胞是一種將其軸突投射到其他腦區域的興奮性神經元細胞，基于此，該生理模型將一個神經元集群（m）中主錐體細胞的動作場電位平均脈沖密度作為另一神經元集群（n）的興奮性輸入。此外，用另一個參數Kmn來表征神經元集群m到神經元集群n的關聯性耦合強度。然后，可以通過設置合適的Kmn參數值來構建不同神經元集群間雙向或者單向的聯系。該生理模型中其他主要參數還包含反饋回路中的抑制性和興奮性增益以及子神經元集群中突觸的平均數等，通過調整這些參數來模擬每個神經元集群的不同活動（癲癇發作時和正常情況下的腦部活動）。
圖1
大腦海馬區域不同神經元集群間相互關系。
Figure1.
Interactions between different neuronal subpopulations of the hippocampus
該模型可用下列隨機偏分方程組表示為:
|
其中，連通常量
，
，突觸時間常量
、
和
以及輸入高斯隨機過程的參數
和
認為是已知的，其值可見參考文獻[11-12]。通過改變模型中參數
、
和
的值來模擬輸出的iEEG 信號為正常腦電信號或者顱內癲癇腦電信號。
2 實驗結果
首先在一個線性隨機系統中驗證和比較 AIC 和 gAIC 兩個準則在 WGCI 算法中的性能，然后將它們用于生理模型所模擬生成的癲癇發作時的iEEG 信號，比較和驗證其性能。
2.1 線性模型
考慮如下的線性隨機系統模型，該模型模擬生成從
到
單向聯系的兩個信號：
|
其中
和
為高斯白噪聲。我們對 16 個不同長度（長度分別為
個點，
）的信號各自都做了 200 次實驗。每次實驗中，分別用 AIC 和 gAIC 準則估計模型階數并將之用于計算上述模型中兩個信號間的 Wiener-Granger 因果關系。正如 1.2 節中所指出的，我們可以知道式（14）中線性模型的真實模型階數為
。首先，我們給出了當信號長度不同的情況下，使用 gAIC 準則得到模型階數的正確率（正確識別指估計得到的 4 個模型階數與真實模型階數完全一致），其結果如圖 2 所示。在圖中，X軸表示信號的 16 個不同長度，即從 128 個點到 2 048 個點，步長為 128 個點，Y軸表示相應的模型階數估計的正確率。從圖 2 中我們可以得到如下結論：當信號長度超過 512 個點時，gAIC 準則已經幾乎可以正確估計出各個不同的模型階數。以信號長度為 2 048 點為例，使用 gAIC 準則得到的四個模型階數分別為｛2,0,1,4｝,即用于 WGCI 估計時所需估計模型系數的項數為 7 項，而使用 AIC 時項數為 4*4 = 16 項，顯然 WGCI 方法中相比于 AIC 準則，采用 gAIC 準則時運算量降低了一半以上。
圖2
由 gAIC 算法估計 4 個模型階數的正確率
The
其次，我們使用 AIC 或者 gAIC 來估計 WGCI，結果如圖 3 所示，將 WGCI 的平均值和標準差作為考察 AIC 或者 gAIC 準則估計參數的評價指標，WGCI 的值越大，標準差越小，表明所采用的準則越有效。圖中的上半部分，對 16 個長度不同的信號，我們畫出了分別使用 AIC 和 gAIC 方法的從信號
到信號
的 Wiener-Granger 因果索引
；圖中的下半部分，則是相應的 16 個長度不同信號的從信號
到信號
的 Wiener-Granger 因果索引
。從該圖可以看出，無論是使用 AIC 準則還是 gAIC 準則，當信號長度逐漸增加時，所有
的平均值都趨近于同一極限（接近 0.713 9），其相應的標準差則趨近于極限值 0.025 8。由模型可知，該系統中兩個信號間的聯系為從信號 1 到信號 2 的單向聯系，故
的值理論上應為 0。由圖 3 中下半部分可知，使用 gAIC 準則的
的均值總是低于使用 AIC 準則得到的
的平均值，所有使用 gAIC 準則的
的均值都快速減少到一個極小值（約為 0.000 4），甚至當信號長度為 256 點時。實際上，當信號長度大于 1 664 點時，使用 gAIC 準則得到的
平均值就已經是零了。對于從
到
的 WGCI 相應的標準差也有上述相同的結論。因此，綜上所述，我們得出如下的結論：基于 gAIC 準則的 WGCI 方法要比基于 AIC 準則的具有更好的魯棒性。
圖3
在線性隨機系統中使用 AIC 和 gAIC 的 WGCI（平均值和標準差）
The
2.2 生理模型
使用 1.3 節中描述的生理模型來模擬兩個神經元集群間具有固定連通模式的癲癇iEEG 信號，每次模擬生成一組長度為 400 s 的信號，采樣頻率為 256 Hz，如圖 4 所示。兩個神經元集群間的連通模式由K值來表示，本實驗中測試了不同的
值（
），而令
，即該兩個神經元集群間的連通模式為單向聯系。對于每次不同K值（例如，每對
（
）），都進行了 250 次實驗，每次實驗的信號長度為 2 048 點。在本節中，我們也比較了 AIC 和 gAIC 兩種準則的性能，將它們分別用于模型階數的估計，圖 5 給出了相應 WGCI 的平均值和標準差。圖中的上半部分展示了 3 個不同
值（
被設為 0）的從
到
的 WGCI 值，同樣地，圖中的下半部分展示了對應的從
到
的 WGCI 值。乍看起來，不管是使用了 AIC 準則還是 gAIC 準則的算法都表明了信號
和
間的連通模式為從
到
的單向聯系。從該圖的上半部分，我們可以看到，當這兩個信號間單向的因果聯系被檢測到時，WGC 索引值隨著耦合強度（即
）的增強而增長。此外，不管是在哪個耦合強度（不同
值）下，使用 gAIC 準則得到的 WGCI 的平均值總是大于使用 AIC 準則得到的 WGCI 平均值。從圖的下半部分可以看到，無論使用哪種準則，從
到
的 WGCI 的平均值都非常低，但是對于相應的標準差，使用 gAIC 準則得到的要更小一些。因此，總體上來說，我們得出如下結論：在用生理模型生成的模擬癲癇iEEG 信號上，使用 gAIC 模型階數估計準則的 WGCI 算法也表現出比 AIC 準則更好的性能。
圖4
由基于生理模型模擬出的iEEG信號例子
左半部分：兩個模擬
left panel: two simulated
圖5
基于生理學模型使用 AIC 和 gAIC 得到的 WGCI（平均值和標準差）
The
3 結論
本文中，使用了一種名為 gAIC 準則的模型階數估計算法，并將它用于 WGCI 方法來檢測信號間的因果關系。實驗結果表明，相比于基于傳統的 AIC 準則的 WGCI 算法，無論是通過一個線性隨機系統生成的信號的實驗還是用生理模型生成的模擬癲癇信號的實驗，結合 gAIC 準則的 WGCI 算法都體現出了更好的魯棒性，可以很好地識別出信號間的因果關系。在將來，計劃進一步研究該算法在真實癲癇腦電信號中的性能，并同時將該算法拓展到多維信號從而檢測信號間直接或者間接的因果關系。
