一、前言
1.学习背景
最近在学习脑网络分析方法时,笔者偶然读到了一篇发表在Physical Review Letters上的文章,文章介绍了一种名为嵌套谱分割(Nested-Spectral Partition, NSP)的方法,用于研究大脑功能网络的分离和整合特性。
传统的脑网络分析往往依赖于图论方法,而NSP方法提供了一个全新的视角,不仅能够揭示网络的层次化模块结构,还能定量衡量网络的整体平衡性。笔者认为这种方法很有意思~
在这篇博客中,笔者将尝试用自己的理解,结合手上有的数据,尝试复现这种方法的实现步骤,包括如何构建功能连接网络、如何进行模块划分,以及如何计算关键网络指标。由于笔者水平有限,如有不当之处,还请各位大佬多多指教。
2.参考文献
Wang, R., Lin, P., Liu, M., Wu, Y., Zhou, T., & Zhou, C. (2019). Hierarchical Connectome Modes and Critical State Jointly Maximize Human Brain Functional Diversity. Physical Review Letters, 123(3), 038301. https://doi.org/10.1103/PhysRevLett.123.038301
Chang, Z., Wang, X., Wu, Y., Lin, P., & Wang, R. (2023). Segregation, integration and balance in resting-state brain functional networks associated with bipolar disorder symptoms. Human Brain Mapping, 44(2), 599–611. https://doi.org/10.1002/hbm.26087
3.分层嵌套谱分割(Nested Spectral Partition, NSP)模型
这个方法就像是在玩俄罗斯套娃,通过特征值分解的方法,将复杂的大脑连接矩阵层层拆解。具体来说,NSP首先将整个大脑网络视为一个整体,然后基于功能连接矩阵进行特征分解,利用特征向量的符号特征,像剥洋葱一样逐层将网络划分为更细致的功能模块。是数据驱动的,可以稳定的自然展示组织方式。
目录
一、前言
1.学习背景
2.参考文献
3.分层嵌套谱分割(Nested Spectral Partition, NSP)模型
二、连接矩阵(C)的构建
1.功能相预处理
1.1 原文的预处理流程
1.2 DPABI预处理流程
2.脑区划分
3.功能连接矩阵计算
3.1 原文计算方法
3.2 DPABI计算方法
三、特征值的分解
1.数学原理
V(特征向量矩阵):
D(特征值矩阵):
V'(V的转置):
2.MATLAB实现方法
四、模块的划分与网络指标的计算
1.层级划分的方式
2.指标计算的数学原理
2.1 核心公式
2.2 整合公式
2.3 分离公式
3.指标计算MATLAB代码实现
二、连接矩阵(C)的构建
在进行NSP分析之前,我们首先需要构建功能连接矩阵(Functional Connectivity Matrix)。如图所示,整个构建过程可以分为三个主要步骤:首先对原始的fMRI进行预处理,包括时间层校正、头动校正等操作,得到预处理后的fMRI数据;然后基于Schaefer2018模板将大脑划分为100个感兴趣区域(ROI);最后,通过计算这100个脑区之间的时间序列皮尔逊相关性,构建出100×100的功能连接矩阵。
1.功能相预处理
1.1 原文的预处理流程
原始研究基于AFNI和FSL软件执行了标准预处理流程。预处理步骤包括以中间层为参考的时间层校正,设置0.3mm为阈值的头动处理,采用MNI模板的空间标准化,使用6mm FWHM高斯核的空间平滑,以及0.01-0.1Hz的带通滤波,值得注意的是未进行全脑信号回归。
1.2 DPABI预处理流程
笔者则采用DPABI软件平台实现ARFCWS预处理流程,即依次进行自动去噪(Automated denoising)、重建(Realignment)、滤波(Filtering)、协变量回归(Covariate regression)、去除波纹(Wave removal)和平滑(Smoothing)处理。
2.脑区划分
在预处理完成后,我们需要对大脑进行功能区域划分。采用Schaefer2018模板,该模板基于Yeo等人2011年提出的7个功能网络,将大脑皮层划分为100个功能区域。
下载链接:Schaefer2018脑图谱
3.功能连接矩阵计算
3.1 原文计算方法
首先提取每个受试者100个ROI的时间序列(每个ROI的时间序列是该区域内所有体素时间序列的平均值),然后采用Pearson相关系数计算ROI之间的功能连接强度。
3.2 DPABI计算方法
俺直接使用DPABI提取ROI信号得到的相关矩阵计算的。
但对于每个受试者得到的100×100初始矩阵,进行了两步关键修正:
- 首先将矩阵对角线上的所有值统一设置为1(因为表示脑区与自身的完全相关)。
- 其次将矩阵中所有负值改为0(由于负连接的生理意义存在争议,且为简化后续网络分析)。
这样得到的修正矩阵C = [cij]100×100既保留了重要的功能连接信息,又便于后续分析使用。
三、特征值的分解
接下来要将功能连接矩阵C分解为C = VDV' (特征向量矩阵V和特征值矩阵D),才可以使用NSP方法划分模块层级。
1.数学原理
核心公式:
接下来要把100×100功能连接矩阵C拆解成三个更简单的矩阵的乘积:
V(特征向量矩阵):
可以理解为"变换方向"。每一列代表一个主要的连接模式。100×100的矩阵。
D(特征值矩阵):
可以理解为"重要程度",对角线上的数字表示每个模式的强度,其他位置都是0.,从大到小排序,越大表示该模式越重要。
V'(V的转置):
就是V矩阵的"镜像",确保我们可以重建回原始矩阵。
2.MATLAB实现方法
% 定义数据处理路径
input_path = 'F:\HCFC';
output_path = 'F:\HCFC\result';
% 检查并创建输出目录
if ~exist(output_path, 'dir')
mkdir(output_path);
end
% 获取所有.mat文件
files = dir(fullfile(input_path, '*.mat'));
% 批量处理文件
for i = 1:length(files)
current_file = files(i).name;
file_path = fullfile(input_path, current_file);
data = load(file_path);
% 提取FC矩阵
if isfield(data, 'FC')
FC_matrix = data.FC;
else
fields = fieldnames(data);
if ~isempty(fields)
FC_matrix = data.(fields{1});
else
continue;
end
end
% 特征值分解
[evals, evecs, proc_FC] = analyze_FC_matrix(FC_matrix);
% 保存结果
results = struct();
results.eigenvalues = evals;
results.eigenvectors = evecs;
results.processed_FC = proc_FC;
results.original_filename = current_file;
[~, name, ~] = fileparts(current_file);
output_filename = fullfile(output_path, [name '_results.mat']);
save(output_filename, '-struct', 'results');
end
四、模块的划分与网络指标的计算
1.层级划分的方式
NSP方法的层级划分遵循一个从整体到局部、逐步细化的过程:
在第一个模式中,根据Perron–Frobenius定理,所有脑区在特征向量中具有相同的符号,此时将整个脑网络视为单一模块,这构成了第一层级。
进入第二个模式时,根据特征向量的正负号将脑区分为两个模块,特征向量为正的脑区组成一个模块,为负的脑区组成另一个模块,形成第二层级的两模块结构。
随后在第三个模式中,基于该模式下特征向量的正负号,将第二层级中的每个模块进一步划分为两个子模块。这个递归划分过程持续进行,直到达到预设层级或每个模块仅包含单个脑区为止。
2.全局指标计算的数学原理
2.1 核心公式
核心公式Hi了反映了网络第i层级的组织特征。
其中Λi2表示特征模式的强度,Mi表示模块数量,(1-pi)是模块大小分布的均匀度惩罚项,N作为归一化因子确保不同网络间的可比性。
这个公式通过综合考虑模块结构的强度、复杂度和均匀性,提供了对网络层级特征的量化评估
2.2 整合公式
选取第一层级(H1)作为整合度的衡量标准。这是因为在第一层级中,所有脑区被视为单一模块,特征向量的符号都相同,反映了整个网络的全局连通性。
2.3 分离公式
从第二层级开始累加所有层级的网络特征值,并进行归一化。
3.指标计算MATLAB代码实现
% 设置数据路径
folder_path = 'F:\HCFC\result';
files = dir(fullfile(folder_path, '*.mat'));
n_files = length(files);
% 初始化结果结构体
results = struct();
results.filenames = cell(n_files, 1);
results.Hin = zeros(n_files, 1); % 整合指标
results.Hse = zeros(n_files, 1); % 分离指标
results.HB = zeros(n_files, 1); % 平衡指标
% 批量处理文件
for i = 1:n_files
% 加载数据
current_file = files(i).name;
data = load(fullfile(folder_path, current_file));
% 提取特征值和特征向量
if isfield(data, 'eigenvalues') && isfield(data, 'eigenvectors')
eigenvalues = data.eigenvalues;
eigenvectors = data.eigenvectors;
% 计算层级指标
N = size(eigenvectors, 1);
[Hi, ~, ~] = calculate_hierarchical_metrics(eigenvectors, eigenvalues, N);
% 计算全局指标
Hin = Hi(1)/N; % 整合指标
Hse = sum(Hi(2:end))/N; % 分离指标
HB = Hin - Hse; % 平衡指标
% 保存结果
results.filenames{i} = current_file;
results.Hin(i) = Hin;
results.Hse(i) = Hse;
results.HB(i) = HB;
end
end
% 生成结果表格并保存
T = table(results.filenames, results.Hin, results.Hse, results.HB, ...
'VariableNames', {'Filename', 'Integration', 'Segregation', 'Balance'});
writetable(T, fullfile(folder_path, 'global_metrics.csv'));
save(fullfile(folder_path, 'global_metrics.mat'), 'results');
总结
需要说明的是,本文仅是笔者对NSP方法的初步尝试,主要聚焦于功能连接网络的分析。原文中还包含了结构连接网络的分析,以及功能-结构整合的深入探讨,这些都为理解大脑的组织原理提供了更全面的视角。感兴趣的小伙伴强烈建议去阅读原文~
笔者也准备明天有时间再尝试计算一下局部指标。