用 ECoG 与 Granger 因果分析绘制脑区间频率特异性信息传递图谱

摘要

高级认知功能不是由孤立脑区完成的，而是由多个解剖上相距很远的脑区组成的大规模网络协同实现。谱 Granger 因果分析（spectral Granger causality, GC）的经典发现之一是「频率通道假说」——前馈信息主要由 gamma 节律携带，反馈信息主要由 beta 节律携带。

本 Protocol 提供了一条端到端、可复现的完整分析管线，把高密度皮层脑电图（ECoG）数据转化为电极级、频率特异性的信息流空间图谱，包含五个核心模块：

1. 数据预处理：分段、带通滤波、降采样、ICA 去伪迹、阵列内共平均重参考、正交源推导、减去 ERP 提取诱发振荡；
2. 时频分析 + 成对 GC 计算（P1 脚本）：连续小波变换（CWT）+ 非参数谱分解（Wilson 算法），对两脑区所有电极对计算双向 GC；
3. 可视化（P2–P3 脚本）：电极对 GC 时间动态、通道×通道连接矩阵、单电极信息流入/流出皮层表面热图（GC map）；
4. 差异 GC 统计图谱（P4 脚本）：两个认知任务相减（记忆任务 − 被动观看），符号翻转置换检验 + FDR 校正，叠加显著性等值线；
5. 层次聚类（P5 脚本）：对差异 GC 模式做 Spearman 相关 + 凝聚式层次聚类，客观鉴定功能电极簇，揭示脑区内功能枢纽的空间组织。

创新点：本方案不引入新的因果推断估计器，而是把从数据预处理、谱分解到地形图可视化与非参数统计映射的整个流程标准化，弥合抽象连接算法与解剖落地网络分析之间的鸿沟。示范数据为猕猴前额叶皮层（PFC）与下颞叶皮层（ITC）在联想记忆任务中的 ECoG 记录。

实验前须知（Before you begin）

• 数据要求：方法原则上适用于任何含神经振荡的多通道电生理时序数据，但需满足两条——① 信号须来自两个不同的通道组（如两个脑区）；② 每组电极须有已知的二维空间排布（如 ECoG 网格）。
• 认知状态对照：要得到「差异 GC 图」，至少需要两种认知状态的数据：
  图-颜色联想任务（PCA，记忆任务）——图片线索呈现 → 延迟期 → 在学会关联的特定颜色选择刺激出现时释放压杆，需要长期记忆提取 + 工作记忆维持；
  图片被动观看任务（PPV，对照任务）——同一批图片，延迟期后注视点变暗即释放压杆获奖励，不涉及记忆提取与维持。
• 动物伦理：示范数据来自日本猕猴（Macaca fuscata）实验，须符合本国及所在机构的动物实验与伦理规范。

关键资源清单（Key resources table）

本 Protocol 为纯计算流程，不涉及湿实验试剂，核心「材料」是软件与数据：

硬件提示：样本数据分析用了 64 GB 内存，建议至少 8 GB 起步，数据集越大需求越高。

分步方案（Step-by-step method details）

第 0 步：ECoG 数据预处理（Timing: 数小时，视数据量而定）

1. 分段：把连续原始 ECoG 信号按任务事件、以线索刺激呈现时刻对齐切分成 trial，时间窗要留足边缘余量，防止时频分析产生边缘伪影；
2. 滤波：带通滤波覆盖感兴趣的频率范围（本例 0.7 Hz 高通等）；
3. 降采样：按奈奎斯特定理降到较低采样率（如 500 Hz），显著降低计算量；
4. ICA 去伪迹：独立成分分析后人工筛查并剔除工频干扰、眼动、局部电极故障等非生理成分，重构信号；
5. 重参考：在每个电极阵列内部计算共平均信号并从各电极中减去，消除全局噪声与阵列特异性波动；
6. CRITICAL正交源推导：每个电极减去其上下左右空间相邻电极的平均电位，提升空间分辨率、抑制容积传导；
7. CRITICAL提取诱发振荡：对所有成功 trial 求平均得 ERP，再从每个 trial 中减去该 ERP，只保留诱发（induced）振荡成分。

第 1 步：时频分析与成对 GC 计算（P1_GC_pw_calc.m；Timing: 样本数据每个条件约 2.5 h）

把预处理后的时域信号经连续小波变换（CWT）转到时频域，用 FieldTrip 的 ft_connectivity_granger 对 Group A ↔ Group B 的所有电极对计算双向谱 GC，输出 gdata 结构供后续所有分析使用。

• ecog0path（数据目录）、analyname（条件文件夹名，如 'PCA'）、postfix（输出后缀，自由命名管理不同分析轮次）；
• trialinfo1：条件名 → trial 编码数组的映射（如 {'all', [1:72]}）；
• grpA_name/grpA_chans（如 'ITC', 1:128）与 grpB_name/grpB_chans（如 'PFC', 129:192）定义两脑区；
• freq_calc_range = [1 72]：限制 GC 计算的频率范围，大幅省内存省时间；
• toi_cwt/toi_GC（如 [-1 3] 秒）：小波变换与 GC 输出的时间窗；
• flg_full_pipeline = 1 跑完整管线，0 则复用已有小波文件重算 GC。

CRITICAL 整套方案靠「差异 GC」分离状态特异的通信动态，必须对实验条件和对照条件各跑一遍（先 analyname='PCA' 跑完，再改成 'PPV' 重跑）。

暂停点 跑完实验条件后，可在此暂停、重新配置参数再跑对照条件。

第 2 步：可视化时间动态与连接矩阵（P2_Plot_GC_TC_and_Matrix.m；Timing: 15 min）

• GC 时间动态图：选一对代表性电极（tc_chan_grpA/tc_chan_grpB），画双向 GC 全时程曲线（如 −0.5 到 3.0 s），基线校正参数 flg_abs_baseline（绝对校正）/ flg_db_baseline（dB 校正）/ flg_zscore 三选一；
• GC 连接矩阵：在 toi_matrix（如 [0.85 2.2] 秒）窗口内对时间取平均，画「A→B」与「B→A」两张通道×通道矩阵，矩阵边缘展示对侧区域平均 GC；
• freq_bands 定义频段（delta 1–4 / theta 4–8 / alpha 8–13 / beta 13–30 / gamma 30–70 Hz），脚本自动逐频段出图；
• 支持批量后台绘图（flg_plot_visible=0，macOS 强烈推荐）与 PNG/PDF/EPS 多格式保存。

第 3 步：信息流空间分布映射（P3_Plot_GC_maps.m；Timing: 15 min）

• 每个电极的值 = 它与对侧脑区所有电极的平均 GC 强度（流出 Outdegree / 流入 Indegree）；
• grpA_layout/grpB_layout 指定 FieldTrip 布局文件（如 LauraR.lay），把电极按物理排布投影到 2D 图上；
• time_windows_to_plot（如 [−0.85,−0.35; 0,0.5; 0.85,1.35; 1.7,2.2]）按基线/线索/延迟等认知分期自动出子图，左列 A→B、右列 B→A。

第 4 步：差异信息流的统计映射（P4_GC_Diff_Stat_and_Map.m；Timing: 15 min，随 nperm 线性增长）

1. 顺序加载两任务的成对 GC 数据（防内存溢出），基线校正后计算差值 ΔGC = GC(PCA) − GC(PPV)；
2. 符号翻转置换检验：随机翻转配对差值符号 nperm 次（≥1000）构建零分布，算每个电极的原始 p 值；固定 random_seed 保证 Step 5 完全可复现；
3. FDR 校正（flg_fdr）：0 = 不校正；1 = Local FDR（Benjamini-Hochberg，单张时频图内）；2 = Global FDR（强烈推荐），跨通道×时间×频段全局校正；
4. 出图：热图颜色编码 ΔGC 幅值，flg_add_sig=1 时用空间插值画显著性等值线（eq_altitude=[0.5 0.5] 控制阈值）；flg_two_sided=1 可双向着色；
5. flg_save_stat=1 导出统计结果供第 5 步复用。

第 5 步：层次聚类鉴定功能电极簇（P5_GC_Cluster.m；Timing: 30 min）

1. 计算所有电极对的差异 GC 向量之间的 Spearman 秩相关——相关高说明两个电极的任务调制通信模式高度相似；
2. 对相关距离矩阵（1 − r）做凝聚式层次聚类，sig_lvel_clust=0.05 被换算成临界相关距离，用来剪树枝状图、划出离散功能网络；
3. 输出三类图：重排序的块对角相关矩阵 + 树状图、颜色编码的皮层表面聚类地图；flg_use_sele_chan=0 时纳入全部电极（推荐全局探索用）；
4. flg_load_stat=1 直接复用第 4 步存好的区域平均置换结果（省时且保证等值线一致）；CRITICAL 两脚本的核心参数（freq_bands、time_windows_to_plot、flg_fdr、nperm、postfix、基线标志）必须完全一致，否则触发参数校验报错。

注意事项

• 内存分配策略：脚本按「Group A 逐通道顺序加载、Group B 整块加载」设计以提速——通道数多的脑区应分给 Group A（GC 双向计算，交换分组不改变结果）；
• FieldTrip 路径：绝对禁止 addpath(genpath(...)) 整树添加，只 addpath 主目录并执行 ft_defaults，并清除旧版本与冲突工具箱路径；
• 平稳性假设：GC 的数学前提是宽平稳。小波非参数谱分解对轻微瞬态非平稳有一定稳健性，但严重漂移或过长分析时窗会产生虚假连接。GC 计算前应目检慢漂移或用 ADF 检验（注意 ADF 只查单位根非平稳，目检仍不可省），必要时多项式去趋势或加严高通滤波。本示范数据的标准预处理（0.7 Hz 高通 + 减 trial 平均 ERP）已足够；
• 统计严格性：Global FDR 极其严格，trial 数少或任务对比弱时可能全军覆没——考虑合并多个 session/被试、或先看未校正 p 值做诊断。

Troubleshooting

预期结果

跑通整条管线，你将得到多层次、频率特异性的脑区信息传递可视化——把「整个脑区当单个节点」的旧范式，推进到皮层表面功能网络的精细空间拓扑：

1. 第一层（P2）：代表性电极对的时间分辨 GC 曲线呈现 ITC ↔ PFC 双向信息流的时间动态；通道×通道连接矩阵呈现整体网络结构；
2. 第二层（P3）：地形热图逐时间窗展示每个电极位置的平均信息流入/流出，直接定位通信「枢纽」；
3. 第三层（P4）：差异 GC 图（PCA − PPV）+ 显著性等值线，锁定任务特异调制的信息流区域；
4. 第四层（P5）：树状图 + 块对角相关矩阵 + 颜色编码皮层地图，揭示任务调制区域内部的功能模块化组织——相似通信模式的电极聚成空间上连续的簇。

原文所有附图（含交互版）均开放获取，theta 频段延迟期的结果仅作示范。

局限性

1. 成对 GC 的固有局限：无法区分直接因果与经第三个未记录脑区中介的间接影响（共同驱动或中继节点）。GC 图应解读为「预测性信息流与功能枢纽」而非严格解剖通路。多变量/条件 GC 理论上可缓解，但数学复杂度与数据需求剧增；
2. 计算量大：CWT + Wilson 谱分解对每一对通道×每个时间点逐一计算、不做时间平均，通道数上千或采样率极高时会内存溢出或跑不完——要么用 freq_step 重降采样，要么上 HPC 集群；
3. 宽平稳假设：严重非平稳数据会产生虚假连接，GC 计算前必须评估（见注意事项）；
4. MATLAB 专属实现：核心算法依赖 FieldTrip，目前不提供 Python/MNE 等开源等价实现，需要 MATLAB 许可证才能复现。

资源可用性

• 通讯作者：Hisashi Tanigawa（hisashi.q@gmail.com）；技术联系人：Xiangsheng Zeng（painter1126@163.com）
• 代码与样本数据：Zenodo, DOI: 10.5281/zenodo.21380181
• 资助：深圳医学研究专项资金（B2504002）、国家自然科学基金（31872776）、JSPS KAKENHI（23H00413）、AMED（JP26wm0625205）