一、前言
1.学习背景
最近在学习脑网络分析方法时，笔者偶然读到了一篇发表在Physical Review Letters上的文章，文章介绍了一种名为嵌套谱分割(Nested-Spectral Partition, NSP)的方法，用于研究大脑功能网络的分离和整合特性。
传统的脑网络分析往往依赖于图论方法，而NSP方法提供了一个全新的视角，不仅能够揭示网络的层次化模块结构，还能定量衡量网络的整体平衡性。笔者认为这种方法很有意思～
在这篇博客中，笔者将尝试用自己的理解，结合手上有的数据，尝试复现这种方法的实现步骤，包括如何构建功能连接网络、如何进行模块划分，以及如何计算关键网络指标。由于笔者水平有限，如有不当之处，还请各位大佬多多指教。
2.参考文献
3.分层嵌套谱分割(Nested Spectral Partition, NSP)模型
这个方法就像是在玩俄罗斯套娃，通过特征值分解的方法，将复杂的大脑连接矩阵层层拆解。具体来说，NSP首先将整个大脑网络视为一个整体，然后基于功能连接矩阵进行特征分解，利用特征向量的符号特征，像剥洋葱一样逐层将网络划分为更细致的功能模块。是数据驱动的，可以稳定的自然展示组织方式。
目录
3.分层嵌套谱分割(Nested Spectral Partition, NSP)模型
二、连接矩阵（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个功能区域。
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方法的初步尝试，主要聚焦于功能连接网络的分析。原文中还包含了结构连接网络的分析，以及功能-结构整合的深入探讨，这些都为理解大脑的组织原理提供了更全面的视角。感兴趣的小伙伴强烈建议去阅读原文～
笔者也准备明天有时间再尝试计算一下局部指标。