【神经科学学习笔记】基于分层嵌套谱分割(Nested Spectral Partition)模型分析大脑网络整合与分离的学习总结

devtools/2024/11/14 20:23:01/

一、前言

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的功能连接矩阵。

我理解的FC矩阵构建过程～
俺理解的FC矩阵构建过程

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.数学原理

核心公式：C = VDV'

  接下来要把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);
endfiles = 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;elsefields = fieldnames(data);if ~isempty(fields)FC_matrix = data.(fields{1});elsecontinue;endend[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:Ncurrent_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) > 1module_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;endendendmodules{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));endideal_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_filescurrent_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
endT = 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方法的初步尝试，主要聚焦于功能连接网络的分析。原文中还包含了结构连接网络的分析，以及功能-结构整合的深入探讨，这些都为理解大脑的组织原理提供了更全面的视角。感兴趣的小伙伴强烈建议去阅读原文～

  笔者也准备明天有时间再尝试计算一下局部指标。


http://www.ppmy.cn/devtools/133297.html

相关文章

ESLint 使用教程（三）：12个ESLint 配置项功能与使用方式详解

ESLint 使用教程（三）：12个ESLint 配置项功能与使用方式详解

前言 在现代前端开发中，代码质量与一致性是至关重要的，ESLint 正是为此而生的一款强大工具，本文将带您详细了解 ESLint 的配置文件，并通过通俗易懂的方式讲解其主要配置项及其配置方法。此外，我们还将探讨一些高级配置…
阅读更多...
【前端】Svelte：动画效果

【前端】Svelte：动画效果

在现代前端开发中，动画效果可以大大提升用户体验，使应用更生动、易用。Svelte 提供了灵活的动画 API，让开发者能够快速实现从简单过渡到复杂动画的各种效果。本文将系统性地介绍 Svelte 的动画功能，并通过多个示例演示如何创建动感…
阅读更多...
前端开发中常见的ES6技术细节分享一

前端开发中常见的ES6技术细节分享一

var、let、const之间有什么区别？ var: 在ES5中，顶层对象的属性和全局变量是等价的，用var声明的变量既是全局变量，也是顶层变量​ 注意：顶层对象，在浏览器环境指的是window对象，在 Node 指的是g…
阅读更多...
【Linux】常用命令（2.6万字汇总）

【Linux】常用命令（2.6万字汇总）

文章目录 Linux常用命令汇总1. 基础知识1.1. Linux系统命令行的含义1.2. 命令的组成 2. 基础知识2.1. 关闭系统2.2. 关闭重启2.3. 帮助命令（help）2.4. 命令说明书（man）2.5. 切换用户（su）2.6.历史指令 3.目录…
阅读更多...
MongoDB笔记03-MongoDB索引

MongoDB笔记03-MongoDB索引

文章目录 一、前言1.1 概述1.2 MongoDB索引使用B-Tree还是BTree？1.3 B 树和 B 树的对比1.4 总结 二、索引的类型2.1 单字段索引2.2 复合索引2.3 其他索引 三、索引的管理操作3.1 索引的查看3.2 索引的创建3.2.1 单字段索引3.2.2 复合索引 3.3 索引的移除3.3.1 指定索…
阅读更多...
算法（第一周）

算法（第一周）

一周周五，总结一下本周的算法学习，从本周开始重新学习许久未见的算法，当然不同于大一时使用的 C 语言以及做过的简单题，现在是每天一题 C 和 JavaScript（还在学，目前只写了一题） 题单是代码随想…
阅读更多...
Chromium127编译指南 Mac篇（五）- 编译Chromium

Chromium127编译指南 Mac篇（五）- 编译Chromium

1. 前言 在Chromium127开发之旅的关键阶段，我们终于来到了编译这一激动人心的步骤。本文将详细指导您如何在macOS环境下成功编译Chromium。通过正确的配置和步骤，您将能够生成一个可运行的Chromium浏览器版本，为后续的调试、功能定制或性能优…
阅读更多...
python视频事物分析

python视频事物分析

使用 OpenCV 和 Hugging Face 的transformers库中的预视视觉模型来分析视频内容 依赖库 pip install opencv-python transformers torch代码 import cv2 import torch from transformers import pipeline# 加载预训练的目标检测模型 detector pipeline("object-detect…
阅读更多...
最新文章

Copyright @ 2022~2023