文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。
导读
最近的研究表明，用于分析和解释大脑灰质(GM)中fMRI数据的数学模型不适用于检测或描述白质(WM)中的血氧水平依赖(BOLD)信号。特别是在一般线性模型中用作回归因子的血流动力学响应函数(HRF)在WM中与GM不同。研究者最近报告了WM中静息态信号时间过程的频率测量结果，显示了不同的功率谱依赖于局部结构-血管-功能之间的联系。对GM的多项研究揭示了区域之间的功能连接性如何随时间动态变化。因此，本文主要研究了静息态WM的BOLD信号是否以及如何随时间变化。研究者测量了体素频谱图，它反映了WM时间进程的时变频谱模式。结果表明，频谱模式是非平稳的，但可以分为五种模式。这些模式显示出其发生和持续时间的不同空间分布，且个体间的分布高度一致。此外，其中一种模式在个体的GM和WM之间具有很强的耦合性，并且根据模式之间转换的层次结构确定了两个WM体素群。此外，这些模式与瞬时HRF的形状相耦合。本研究结果扩展了之前的研究，揭示了BOLD信号频谱模式随时间的非平稳性质，提供了WM中静息态信号的时空频率表征。
前言
通过测量与神经元能量消耗相结合的血氧水平依赖(BOLD)信号的变化，功能性磁共振成像(fMRI)已被广泛用于映射大脑活动图。有研究者已经开发了数学模型来识别图像体素，这些体素表现为参与者在扫描期间执行的任务所触发的信号振幅增加，从而映射与任务相关的大脑激活。此外，在静息态下测量的BOLD信号波动反映了神经活动的自发变化，并揭示了根据大脑区域之间时间同步程度表征的内在功能网络。虽然灰质(GM)中BOLD信号的研究在fMRI文献中占主导地位，但关于白质(WM)BOLD效应的描述却少得多，部分原因是与GM相比，白质BOLD的血流/体素更小，因此其效应较弱。然而，越来越多的证据表明，可以在WM中可靠地检测到BOLD信号并反映神经活动，而且这些信号在神经或精神疾病的患者中有着显著变化。最近的研究揭示了WM BOLD对刺激反应的时域曲线与GM的不同之处，表现为幅值降低和峰值延迟，这反映了GM和WM之间血流动力学条件的差异。这些研究强调了了解WM时间进程特征的重要性，以便将其纳入适当的检测模型来提高灵敏度。
最近研究表明，在低频段(0.01-0.08 Hz)上，WM时间进程的功率谱与GM的功率谱不同，它们随位置的变化而变化，并取决于局部WM的神经血管和解剖结构。此外，不同位置的频谱模式预测了它们在功能整合以及特定人类行为中的变化。这些发现表明，WM信号的频率内容在空间上是异质的。
本研究数据表明，时间平均频谱实际上是五个子模式的组合，称之为频谱模式（或频率模式），其随时间反复出现。它们包括四种单峰模式（每种都在不同频率下表现出单一的频谱峰值，模式1-4）和一种均匀模式（功率在相对较低的水平上均匀分布在感兴趣的频段上，模式5）。这些模式的出现和持续时间在WM体素中呈现出不同的空间分布，并且这些分布在个体间高度一致。其中一种模式（模式1）在个体间的GM和WM之间的发生表现出强烈的耦合。更有趣的是，根据模式之间转换的层次结构确定了两个WM体素群。具体而言，在体素群1中，模式1、2和5之间的转换概率显著高于体素群2。相比之下，在体素群2中，模式3、4和5之间的转换概率显著高于体素群1。此外，瞬时频谱模式与随时间变化明显的HRF有关。本研究结果通过观察BOLD频谱的动态变化，揭示了频谱模式的非平稳性质以及它们随时间的生理耦合，从而提供了WM中静息态BOLD信号的时空频率表征。
研究方法
HCP数据
从HCP健康青年数据库中随机抽取200人（男性88人，女性112人，年龄为22-35岁）。成像协议在以前的报告中有详细描述。简而言之，所有图像均使用3T Siemens Skyra扫描仪（Siemens AG，Erlanger，Germany）获得。使用多波段梯度回波平面回波成像(EPI)获得了两个session的静息态图像。每个session由两次run（从左到右和从右到左相位编码）组成，每次run的时间为14min 33s，重复时间(TR)=720ms，回波时间(TE)=33.1ms，体素大小=2mm各向同性，体素数=1200。在整个fMRI扫描过程中记录生理数据，包括呼吸和心脏波动。作为结构参考，T1图像是使用3D磁化快速梯度回波(MPRAGE)采集获得的，TR=2400ms，TE=2.14ms，体素大小=0.7mm各向同性。
预处理
从HCP库中获得的预处理图像经过了最小预处理管道(MPP)。简而言之，T1图像使用FNIRT非线性配准到MNI空间，并通过Freesurfer管道产生体素/表面分割以及形态测量。对于fMRI，管道包括头动校正、使用反相编码方向的磁化率失真校正以及对MNI空间的非线性配准。本研究进行了额外的处理，例如对干扰变量的回归，包括12个头动参数（三个平移、三个旋转及其导数），以及通过RETROICOR方法建模的呼吸和心脏波动，然后校正线性趋势和使用带通滤波器(0.01–0.08Hz)进行时域滤波。通过平均所有被试的WM分割（源自Freesurfer），并设置阈值为0.95来重建WM掩膜，该阈值用于限制WM内的计算。由于其中较高的个体可变性，还以类似的方式使用较低的阈值来重建GM掩膜。选取0.65作为最佳阈值，首先设置阈值为0.5，然后以0.01的步长逐渐增加，直到WM和GM之间没有重叠体素。最后，为了提高信噪比并减少计算量，将预处理后的图像从2mm分辨率降采样至3mm。
频谱模式的计算
频谱模式的计算流程如图1所示。首先，使用短时傅里叶变换(STFT)计算每个体素的BOLD信号的频谱图。首先将整个时间进程分割成一组部分重叠的窗口，并计算每个窗口的变换，然后计算每个窗口的功率谱，显示不同频率下BOLD波动随时间变化的幅值。本研究选择了100.4s(140TR)的窗长以提供低至0.01Hz的频率分辨率，这被认为是基线大脑活动频带的低截止值，并且在之前的几项研究中提到过。为了提供足够的时间分辨率，同时避免过多的计算负载，本研究选择了2.88s(4TR)的步长。这就有266个窗口与每个相邻窗口重叠97.52s(136TR)，均匀地横跨BOLD时间过程的持续时间。通过获得跨体素和被试的频谱图，获得了大量跨时间的频谱模式观测结果。然后使用K-means方法将这些观察结果分成一组团簇。这些模式代表团簇的质心。然后可以根据其到质心的欧式距离将每个窗口频谱模式分配给特定模式。采用肘部法则确定聚类数的最佳值。具体来说，连续将聚类数从1增加到20，并绘制了观测到质心距离之和与聚类数的关系曲线。聚类数的最佳值由肘部在曲线上的位置确定。在实践中，为了找到肘点，研究者使用了一个函数，该函数沿着曲线一次移动一个二等分点并拟合两条线，一条到二等分点左侧的所有点，一条到二等分点右侧的所有点。肘部被判断为处于二等分点，使两次拟合的误差总和最小化(详见https://www.mathworks.com/matlabcentral/fileexchange/35094knee-point)。为了进一步验证聚类数的可靠性，研究者使用自举法随机抽样频谱特征，并应用k-means聚类和肘部法则。这个过程重复1000次，计算出聚类数的直方图。
图1.频谱模式的计算流程。
计算GM和WM模式间的耦合
在GM的模式A和WM的模式B之间耦合的情况下，对于被试i，首先提取gi，即模式A在GM体素上的平均发生率，GM体素显示出模式A的显著高发生率。对所有被试重复此过程，得到一个向量{gii=1,2,3,...,200}。然后对于WMx中的一个体素，提取模式B的发生率，wxi。对所有被试重复此过程，得到一个向量{wxii=1,2,3,...,200}。然后对WM中的每个体素执行{gi}和{wxi}之间的Pearson相关性，在所有WM体素{x}上产生r值的分布。WM体素的r值越高，表明模式A在GM中的发生率越高，预示着在同一被试的该体素中出现模式B的频率越高。
逐窗HRFs的估计
使用盲反卷积方法从每个被试的每个时间窗的静息态时间过程b(t)估计窗内HRFs。该方法不需要关于HRF的先验假设，并且基于这样的概念：相对较大的BOLD信号峰值代表可分离的、主要的、自发事件的发生。在本研究中，此类事件被检测为超出指定阈值的峰值（这里为平均值的1.5个标准差以上）。对于每个事件，使用sn(t)、事件的起始和h(t)的组合来拟合一般线性模型，h(t)表示两个伽马函数及其时间导数的线性组合。这里n表征从起始到脉冲函数s(t)的延迟时间，其中s(t)=1仅当t对应于检测到的峰值（事件）。双伽马函数与时间导数一起能够模拟具有初始倾角的一般HRF。通过搜索一个n(n∈0-12s)，最小化残差cov[b(t)–conv(sn(t)，h(t))]的协方差，可以估计出n的值和模型h(t)的几个参数，因此可以获得HRF hn(t)。
结果
频谱模式
如图1所示，使用跨时间窗、WM体素和被试的加窗频谱模式的k-means聚类得出五种频谱模式。模式1-4，称为单峰模式，在不同频率上呈现单个尖峰，而模式5是均匀模式，频谱功率相对恒定且在感兴趣的频段上较低。相同的工作流程在GM中产生了两种模式，包括与WM中模式1相似的单峰模式和类似于WM中模式5的均匀模式。
模式的发生和持续时间
图2显示了每种模式显著高发生率（p＜0.05，FWE校正）的体素分布。同时，与其他体素相比，室旁体素和颞体素表现出显著更高的模式3和4发生率。然而，模式2的效应不明显，因此未在图2中显示。注意，模式2的发生率并不明显低于图2所示的其他模式。在图2的每个板块中，箱线图表示在转换到另一种模式之前发生的次数和持续时间的组水平分布。可以观察到，模式5在图2所示的所有区域中发生频率最高且持续时间最长。此外，在额下区，模式1比模式3和4发生频率更高且持续时间更长，而在室旁和颞区则相反。总体而言，发生次数和持续时间在被试间高度一致，这反映在箱线图紧密的四分位间距上（上限和下限之间的距离）。
图2.每种模式高发生率的体素分布示意图。
图3显示，在GM中两种模式高发生率的体素分布。与GM的其余部分相比，靠近大脑表面的皮层区域显示出，模式1具有显著更高的发生率。同时，岛叶、颞区和扣带回的体素显示出模式2的发生率明显高于其他GM体素。
图3.GM中两种模式的显著高发生率区域。
GM和WM之间模式发生率的耦合
图4显示了在特定GM和WM体素中普遍存在的特定模式的发生率耦合。左图中WM体素强度显示了模式1的发生率与图3(左侧)GM区域中模式1的平均发生率之间的Pearson相关系数(r)。可以观察到，在选定的表面皮层区域中模式1的高发生率可靠地预测了同一被试的WM体素中模式1的高发生率，反之亦然。其中，GM中模式1的高发生率也预测了在同一被试的类似WM区域中的模式3和4的低发生率，这表明WM中的这些模式之间存在抑制关系。
图4.WM与GM体素中模式发生率的耦合。
图4右图中的体素强度显示了模式5的发生率与图3(右侧)GM区域中模式2的平均发生率之间的Pearson相关系数(r)。只有部分WM体素中模式5的发生率与GM区域中的模式2存在显著耦合。
模式之间的转换
图5映射了每个WM体素的五种模式之间的转换次数（用T值显示）。可以观察到，与WM的其他区域相比，额下区的转换次数明显更少。
图5.每个WM体素的转换次数。
图6显示了从模式i到j的体素转换显著高于其他WM体素。可以清楚地识别两个体素群，并用不同的颜色进行编码。左下图显示了在模式1、2和5（集群1）之间存在显著更高转换的体素，右侧显示了五种模式之间的转换概率。右下图显示了在模式3、4和5（集群2）之间表现出显著更高转换的体素，右侧显示了五种模式之间的转换概率。请注意，关于每个集群的图是通过对所有被试的模式1、2、5或模式3、4、5之间的转换图进行单样本t检验（p＜0.05，FWE校正）生成的。由此观察到，两个集群之间的转换概率模式是一致的。
图6.模式之间的转换。
瞬时HRF与频谱模式之间的关系
如图7所示，加窗后的HRFs也随时间变化，其形状与瞬时频谱模式密切相关，即初始倾角减小和到达峰值时间增加的HRF与在功率谱中识别的高频分量增加有关。这五种频谱模式对应于不同的HRFs。模式1-4的峰位移（如图1所示）反映了HRF中初始下降幅度的变化以及达到峰值的时间。由于模式5在感兴趣频段上的频谱功率相对恒定，因此其相关HRF的初始下降和峰值时间似乎是模式1-4的平均值。
图7.瞬时HRF与频谱模式之间的关系。
结论
本研究评估了静息态下WM中BOLD信号的频谱图，反映了BOLD时间进程中功率谱模式的时间变化情况。将k-means聚类应用于来自不同体素和被试的不同时间点，从而获得了五种主要模式，相当于聚类质心点。这五种模式与不同的HRFs相关。一些参数，包括发生率、持续时间和转换次数被量化，并在被试间显示出较强的一致性。除此之外，其中一种模式在GM和WM之间的发生率具有强耦合性。总而言之，本研究是白质中BOLD静息态时间进程的时空频率表征的第一份报告。本研究结果揭示了WM中BOLD波动的频谱模式的非平稳性质。BOLD频谱的测量具有高度一致性和可重复性，并揭示了功率谱的重复出现模式以及之前在WM BOLD信号中未曾报告过的HRFs，这增加了我们对人脑中时空频率生理关联的理解。
原文：Li, M. , Gao, Y. , Anderson, A. W. , Ding, Z. , & Gore, J. C. . (2022). Dynamic variations of resting-state bold signal spectra in white matter. NeuroImage, 250, 118972.
MRI图像可从HCP网站下载：
https://www. humanconnectome.org/
本文用于分析数据的工具可以从以下位置下载：
rsHRF：https://www.nitrc.org/projects/rshrf
PhysIO工具箱：https://www.nitrc.org/projects/physio/ Circos
网友评论