fMRI脑影像特征提取——Pearson相关与低阶功能连接LOFC（dpabi+nilearn）

一个连接由节点和连接强度（边）组成。功能连接中的节点可以是单个体素，也可以是由脑图谱定义的脑区ROI。通常的做法是将脑图谱中多个或者全部ROI作为节点，或者是将脑图谱中单个ROI与全脑体素作为节点。定义连接强度（边），最常用的方式是Pearson相关，两个节点之间可能没有物理连接，但如果他们的BOLD信号同步性很强，则它们之间可能存在很强的功能上的协作，即很强的功能连接，这是一个很好的表征两个节点

意疏

19405人浏览 · 2021-03-30 10:00:50

意疏 · 2021-03-30 10:00:50 发布

fMRI功能连接矩阵构建流程
| 图源

由于暂未发现MRI扫描对人体的伤害，加上MRI影像的超高空间分辨率，所以MRI脑影像分析理所当然的成为了脑科学研究最重要的手段之一。MRI脑影像主要包括两种，结构MRI脑影像（此后简称sMRI脑影像）和功能MRI脑影像（此后简称fMRI脑影像），前者反映大脑的组织结构，后者反映大脑的活动代谢。MRI影像是3D的，一幅图被称为一个volume，图中一个点被称为一个体素，体素以一定空间分辨率（如 $3mm^3$ ）近似的一一对应大脑的某个解剖位置。

fMRI脑影像通过BOLD信号来成像。BOLD（blood oxygenation level dependent，血氧依赖水平）信号，可以间接反映大脑的活动代谢情况，某个区域BOLD信号高表示该区域正忙于处理信息或者说被激活，是大脑功能性的体现。fMRI脑影像由多个volume组成，是一个4D序列，如上图左上角所示。一个volume代表一个时间点的影像，影像中的一个体素的值，表示空间位置与之对应的大脑某个小体积内的BOLD信号值。如果把序列中每个volume的相同大脑位置（如 $ROI_i$ ）的值都取出来，按照先后先后顺序排在一起，就形成了一个该大脑位置（如 $ROI_i$ ）的随时间变化的BOLD信号或者时间序列，如上图右上角所示。

某个大脑位置（如 $ROI_i$ ）的BOLD信号，反映了在fMRI脑影像采集过程中该大脑位置（如 $ROI_i$ ）活动代谢的变化。我们有理由相信，如果两个大脑位置BOLD信号变化的同步性很高，即大脑活动代谢的同步性很高，那么这两个大脑位置在功能上应该存在很强的联系，即存在很强的功能连接（Functional Collectivity，FC）。

常用Pearson相关系数来评估两个信号的同步性，两个大脑位置（如 $ROI_i$ 与 $ROI_j$ ）的BOLD信号的Pearson相关系数值，就称作它们之间的功能连接强度。如果将所有大脑位置BOLD信号之间两两做Pearson相关，就会形成一个对称的功能连接矩阵（FC Matrix）。功能连接矩阵可以看做是大脑功能协作性的一种体现，是基于4DfMRI脑影像提取的一种有效的描述大脑功能的特征。通过这个特征，我们可以做进一步的处理与分析。功能连接一般是做静息态fMRI分析，但是也可以用于任务态。

copyright © 意疏：https://blog.csdn.net/sinat_35907936/article/details/115253078。初学小白，若有错误，还望不吝指出。

皮尔逊Pearson相关

一个连接由节点和连接强度（边）组成。功能连接中的节点可以是单个体素，也可以是由脑图谱（atlas，parcellation）定义的脑区ROI，如下图中一个彩色块就代表一个ROI（脑区）。通常的做法是将脑图谱中多个或者全部ROI作为节点来两两建立连接，或是以多个ROI为种子点，以另外多个ROI作为目标点来建立连接，又或者是以脑图谱中单个ROI为种子点，以除此ROI外的全脑体素作为目标点来建立连接。第一中将形成一个对称的功能连接方阵，第二中将形成一个功能连接矩阵，第三个将形成一个功能连接谱（map）。暂且把第一种称基于ROI的功能连接，把后两种称作基于种子点的功能连接。

AAL脑图谱
定义连接强度（边），最常用的方式是Pearson相关，两个节点之间可能没有物理连接，但如果他们的BOLD信号同步性很强，则它们之间可能存在很强的功能上的协作，即很强的功能连接，这是一个很好的表征两个节点间关系的特征。如果一个节点包含不只一个体素（通常情况），那该节点的BOLD信号为该节点内所有体素BOLD信号的平均。

现在想想同步性在描述的是什么？很容易想到，就是线性关系，如果一个信号值变大，另一个信号值也跟着变大，一个信号值变小，另一个信号值也跟着变小，或者相反，那么以两个信号的所有取值对应作为坐标画成的散点图，一定会在一条直线附近分布。Pearson相关，是两个信号或者随机变量线性相关性的衡量，值在[-1,1]，当值为1或者-1时，是完全的线性相关，当值为0时，是完全的线性无关。其定义式如式（1），可见其是归一化了的协方差，而协方差正是在衡量两个随机变量或者信号X,Y之间的协同性。

$\tag 1 r = \cfrac {cov(X,Y)}{\sigma_X \sigma_Y}=\cfrac {E((X_t-\overline X)(Y_t- \overline Y))} {\sigma_X \sigma_Y}$

其中 $X_t,Y_t$ 表示两个随机变量的第 $t$ 个取值，或者两个信号第 $t$ 个时刻的信号取值。

Pearson相关计算式由式（3）给出，如果把信号或者随机变量X，Y的T次取值依次排列成T维向量X，Y，那么信号或者随机变量X与Y的Pearson相关就等于T维向量X，Y夹角的余弦，如式（4），我们知道余弦值可以表示两个向量的共线程度，这也是一种各方向值变化（等价于原信号的各取值变化）的协同程度的衡量。
$\tag 2 r = \cfrac {\cfrac 1 T \sum^T_{t=1}(X_t - \overline X)(Y_t- \overline Y)} {\sqrt {\cfrac 1 T\sum^T_{t=1} (X_t - \overline X)^2} \sqrt {\cfrac 1 T\sum^T_{t=1} (Y_t - \overline Y)^2}}$

$\tag 3 = \cfrac { \sum^T_{t=1}(X_t - \overline X)(Y_t - \overline Y)} {\sqrt {\sum^T_{t=1} (X_t - \overline X)^2} \sqrt {\sum^T_{t=1} (Y_t - \overline Y)^2}}$

$\tag 4 = \cfrac {(\bold X- \overline X)^T(\bold Y - \overline Y)} {|| (\bold X - \overline X)||_2||\bold Y- \overline Y||_2 } =cos \angle(X,Y)$

低阶功能连接——LOFC

LOFC

两个节点BOLD信号直接相关（默认Pearson相关）建立的功能连接，被称为低阶功能连接（LOFC），它也常被直接称作功能连接FC。LOFC作为fMRI脑影像特征的稳定性可靠性已经有相关文章证明，至于高阶功能连接HOFC，笔者后面再写。就基于ROI的LOFC而言，所有的节点的BOLD信号两两做Pearson相关，就可以得到LOFC矩阵，它是对称方阵。如下图中，第 $i$ 个节点 $ROI_i)$ 与第 $j$ 个节点 $ROI_j)$ 之间的低阶功能连接强度由式（5）给出，图源。

$\tag 5 LOFC_{ij} = \cfrac { \sum^T_{t=1}(ROI_i(t) - \overline {ROI_i})(ROI_j(t) - \overline {ROI_j})} {\sqrt {\sum^T_{t=1} (ROI_i(t) - \overline {ROI_i})^2} \sqrt {\sum^N_{i=1} (ROI_j(t) - \overline {ROI_j})^2}}$

其中， $ROI_i(t)$ 为第 $t$ 个时间点时， $ROI_i$ 上的BOLD信号值， $\overline {ROI_i}$ 为 $ROI_i$ 上T个时间点的BOLD信号值的平均。
fMRI功能连接矩阵构建流程

LOFC实现

上图写明了求LOFC的基本流程：根据需要选择合适的脑图谱，将被试的fMRI影像配准到脑图谱，根据脑图谱提取每个ROI上的平均BOLD信号，相关。

用dpabi可以很轻易的算LOFC：勾选操作，选择图谱。
dpabi算功能连接

Functional Connectivity： 以图谱中每个ROI为种子点，计算其BOLD信号与全脑体素的BOLD信号的相关，将形成与种子点个数相同的LOFC谱（map），如果以AAL116脑图谱做为ROI划分标准，一个被试将会形成116个LOFC谱（map），其中之一如下图所示。种子点需要根据我们的实际任务来定义。

功能连接谱

Extract ROI courses： 提取图谱中每个ROI的BOLD信号，并两两做相关，形成一个对称的LOFC矩阵，如果以AAL116脑图谱做为ROI划分标准，一个被试将会形成一个116x116的对称LOFC矩阵，如下图所示。这是最常见的求LOFC的方式。

load('E:\LAB\rMRI\test\Results\ROISignals_FunImgARWCFS\ROICorrelation_sub-10517.mat');

ROICorrelation = ROICorrelation - diag(diag(ROICorrelation));
figure(1);
imagesc(ROICorrelation);
colorbar

用nilearn也可以很轻易的算LOFC：加载数据，选择图谱，提取BOLD信号（时间序列），相关，显示。代码参考：nilearn官方例程

import numpy as np
import scipy.io as io
from nilearn import plotting
from nilearn import datasets
from nilearn.image import load_img
from matplotlib import pyplot as plt
from nilearn.input_data import NiftiLabelsMasker
from nilearn.connectome import ConnectivityMeasure

fMRIData = load_img(r'E:\LAB\rMRI\test\FunImgARWCFS\sub-10517\sFiltered_4DVolume.nii')

# atlas

atlas = datasets.fetch_atlas_aal()
atlas_filename = atlas.maps
labels = atlas.labels

# Extracting times series
masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True,
                           memory='nilearn_cache', verbose=5)
time_series = masker.fit_transform(fMRIData)

# correlation
correlation_measure = ConnectivityMeasure(kind='correlation')
correlation_matrix = correlation_measure.fit_transform([time_series])[0]
print(correlation_matrix)

# Plot the correlation matrix
plt.figure()
np.fill_diagonal(correlation_matrix, 0)
print(correlation_matrix)
plt.imshow(correlation_matrix, interpolation="nearest", cmap="RdBu_r")
plt.colorbar()
plt.show()

# save correlation result
io.savemat('correlation_matrix.mat', {'correlation_matrix': correlation_matrix})

# plot connectome with 99% edge strength threshold in the connectivity
coordinates = plotting.find_parcellation_cut_coords(labels_img=atlas_filename)
plotting.plot_connectome(correlation_matrix, coordinates,
                         edge_threshold="99%",
                         title='AAL 116')
plotting.show()

116x116的LOFC矩阵：

只保留1%的边时的功能连接体，体现了116个ROI之间较强的功能连接。

功能连接体

从可视化后的LOFC矩阵来看，忽略颜色差异，两个工具形成的图是基本一致的。但是如果仔细对比两个工具产生的LOFC矩阵，会发现一些差异，下图是两者差异的可视化。笔者看过两者的代码，都是做的皮尔逊相关，按理说就算差异存在，也应该差异非常小，与实际情况有所不同，笔者暂时也还没有发现原因。（差异应当是两个工具对AAL图谱采用了不同的重采样方法，导致ROI的分割有偏差导致的）

load('E:\LAB\rMRI\test\Results\ROISignals_FunImgARWCFS\ROICorrelation_sub-10517.mat');load('E:\PycharmProjects\NImage\correlation_matrix.mat');
ROICorrelation = ROICorrelation - diag(diag(ROICorrelation));
err = ROICorrelation - correlation_matrix;
magesc(err);
colorbar