引言

癫痫是严重影响人类日常生活的脑部疾病之一。其中难治性癫痫综合征的药物治疗效果差，严重影响了患者的生活质量，解决手段之一就是通过手术来切除致痫区（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 的基本思想如下：如果一个过程s₁“导致”另一过程s₂ 的发生，则s₁ 过去的值包含了帮助预测s₂ 的信息，并且此信息大于单独从s₂ 过去的值所包含的信息。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）的兴奋性输入。此外，用另一个参数K^mn来表征神经元集群m到神经元集群n的关联性耦合强度。然后，可以通过设置合适的K^mn参数值来构建不同神经元集群间双向或者单向的联系。该生理模型中其他主要参数还包含反馈回路中的抑制性和兴奋性增益以及子神经元集群中突触的平均数等，通过调整这些参数来模拟每个神经元集群的不同活动（癫痫发作时和正常情况下的脑部活动）。

图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 个模型阶数的正确率

X轴表示 16 个不同长度的信号（128*n（n = 1,2,,16）），Y轴表示相应的正确率

Figure2. Percentage of correct estimation of the four model orders estimated by gAIC

The X-axis represents the 16 different lengths of the signals, i.e. 128*n (n = 1, 2, , 16), and the Y-axis displays the corresponding correct percentage of estimation

图选项

下载全尺寸图像

下载幻灯片

其次，我们使用 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（平均值和标准差）

X轴表示 16 个不同长度的信号，Y轴表示相应的 WGCI 值

Figure3. WGCI (mean value and standard deviation) using either AIC or gAIC for the linear stochastic system

The X-axis represents the 16 different lengths of the signals; the Y-axis shows the corresponding WGCI values

图选项

下载全尺寸图像

下载幻灯片

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信号例子

左半部分：两个模拟iEEG信号；右半部分：对应的功率谱密度。模型参数为：由单向连接的两个集群产生的快速准正弦波（25 Hz）活动（类似于癫痫发作时观察到的信号），方向为从集群 1（信号Y₁）到集群 2（信号Y₂）

Figure4. Example of iEEG signals simulated by the physiology-based model

left panel: two simulated iEEG signals; right panel: corresponding power spectral densities. Model parameters were such that a fast quasi-sinusoidal (25 Hz) activity (similar to that observed at seizure onset) was generated by the two populations that were unidirectionally coupled, the direction of connectivity being imposed from population 1 (signal Y₁) to population 2 (signal Y₂)

图选项

下载全尺寸图像

下载幻灯片

图5 基于生理学模型使用 AIC 和 gAIC 得到的 WGCI（平均值和标准差）

X轴表示不同的耦合参数（上半部分为K¹²，下半部分为K²¹），Y轴表示对应的 WGCI 值

Figure5. WGCI (mean value and standard deviation) using AIC and gAIC for the physiology-based model

The X-axis represents the different values of the coupling parameters (K¹² for the upper panel and K²¹ for the lower panel); the Y-axis displays the corresponding WGCI values

图选项

下载全尺寸图像

下载幻灯片

3   结论

本文中，使用了一种名为 gAIC 准则的模型阶数估计算法，并将它用于 WGCI 方法来检测信号间的因果关系。实验结果表明，相比于基于传统的 AIC 准则的 WGCI 算法，无论是通过一个线性随机系统生成的信号的实验还是用生理模型生成的模拟癫痫信号的实验，结合 gAIC 准则的 WGCI 算法都体现出了更好的鲁棒性，可以很好地识别出信号间的因果关系。在将来，计划进一步研究该算法在真实癫痫脑电信号中的性能，并同时将该算法拓展到多维信号从而检测信号间直接或者间接的因果关系。