一、前言
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
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);
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代码实现
function [Hi, Mi, pi_values] = calculate_hierarchical_metrics(eigenvectors, eigenvalues, N)
Hi = zeros(N, 1);
Mi = zeros(N, 1);
pi_values = zeros(N, 1);
Mi(1) = 1;
pi_values(1) = 0;
modules = cell(N, 1);
modules{1} = ones(N, 1);
Hi(1) = (eigenvalues(1)^2 * Mi(1) * (1-pi_values(1)))/N;
for level = 2:N
current_modules = modules{level-1};
new_modules = current_modules;
current_vec = eigenvectors(:, level);
unique_mods = unique(current_modules);
next_module_id = max(unique_mods);
for m = 1:length(unique_mods)
module_mask = current_modules == unique_mods(m);
module_nodes = find(module_mask);
if length(module_nodes) > 1
module_vec = current_vec(module_nodes);
pos_nodes = module_vec > 0;
if any(pos_nodes) && any(~pos_nodes)
next_module_id = next_module_id + 1;
new_modules(module_nodes(pos_nodes)) = next_module_id;
end
end
end
modules{level} = new_modules;
Mi(level) = length(unique(new_modules));
module_sizes = zeros(Mi(level), 1);
unique_mods = unique(new_modules);
for m = 1:Mi(level)
module_sizes(m) = sum(new_modules == unique_mods(m));
end
ideal_size = N/Mi(level);
pi_values(level) = sum(abs(module_sizes - ideal_size))/N;
Hi(level) = (eigenvalues(level)^2 * Mi(level) * (1-pi_values(level)))/N;
end
end% 计算所有被试的全局指标并保存结果
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方法的初步尝试，主要聚焦于功能连接网络的分析。原文中还包含了结构连接网络的分析，以及功能-结构整合的深入探讨，这些都为理解大脑的组织原理提供了更全面的视角。感兴趣的小伙伴强烈建议去阅读原文～
笔者也准备明天有时间再尝试计算一下局部指标。
