基于滑动窗口法的相关性能否揭示静息态fMRI中的动态功能连接?

在过去几年中,静息态功能性磁共振成像(fMRI)的研究重点已从扫描过程中平均时间序列的相关性分析转移到扫描过程中功能连接的变化分析。尽管一些研究发现动态功能连接(dFC)是存在的,但对dFC的统计评估不一定合理,在某些研究中此过程甚至被省略。
这篇文章解释了使用适当的统计检验方法检测dFC的必要性,讲述了如何进行统计检验以及如何评估dFC的性能。作者使用BOLD fMRI 数据说明了全身麻醉的猕猴被试和静息态人类被试的fMRI情况。实验结果有两方面,一方面通过模拟时间序列,作者发现在10分钟内,使用基于滑动窗口法的相关性方法几乎不可能检测到dFC这一假设在猕猴被试和人类被试中都得到了验证,在任何单次扫描实验中没有发现dFC存在的证据
其次,通过多次扫描实验或跨被试扫描求平均可以显着提高检测dFC的能力,作者最后发现大多数功能连接实际上确实是动态的。通过这项研究,作者希望人们加深使用统计检验方法对dFC进行评估的认识,以及如何使用适当的方法来避免探究FC方法中的陷阱。这篇文章发表在Neuroimage期刊。(可添加微信号siyingyxf18983979082获取原文,另思影提供免费文献下载服务,如需要也可添加此微信号入群,原文也会在群里发布)
(对文中相关术语进行解释):
1、原文中session在神经影像学中表示的是一个持续数分钟的扫描实验。
2、统计检验(test statistic):又叫“假设检验”。是根据抽样结果,在一定可靠性程度上对一个或多个总体分布的原假设做出拒绝还是不拒绝(予以接受)结论的程序。
3、零假设(null distribution):指进行统计检验时预先建立的假设。零假设成立时,有关统计量应服从已知的某种概率分布。
4、零分布(null distribution):H0的假设下样本的分布。
5、替代数据法(surrogate data):替代数据法是一种统计检验方法,它假设测得的时间序列数据是线性的,并依据该假设产生相应的一组替代数据,然后分别计算原始数据和替代数据集的检验统计量。
6、统计假设检验(statistical hypothesis test):主要是验证所选的模型和所解释的公式,在结构上、形式上、变化方向上是否能代表客观情况。一般采用最小平方法解模时必须进行统计假设检验,即应用统计推断的假设检验原理,通过特定的统计量进行各种显著性检验。
引言
基于血氧水平依赖(BOLD)的静息态的功能磁共振成像(fMRI)研究大多集中于扫描过程中静息状态的功能连接(FC)模式。最近大多数研究开始涉及时间子段中FC的波动,即动态功能连接(dFC),而前者称为静态功能连接(sFC)。Hutchison等人在关于dFC的研究中已取得相关进展,研究dFC的最常见和最直接的方法是使用滑动窗口法探究FC,它包括在连续的时间窗口段上计算FC来度量,例如计算Pearson相关系数或计算信号的锁相值,这给出了FC时间序列,可用于评估扫描过程中FC的波动情况。虽然这种分析FC的方法看起来很简单,但是在之前研究中有两个陷阱没有认识到。第一个陷阱是将检验统计的观察值与潜在值一起确定,这意味着在观察到的FC时间序列中存在波动被视为存在dFC的证据。第二个陷阱是忽略了这样一个事实,即观察到的FC值是对真实(和不可观察)值的估计,因此受统计不确定性的影响。
      因此,要确定观察到的FC时间序列的波动是由于统计不确定性引起的还是大脑FC的真实变化引起的,必须进行适当的统计检验。这通常通过计算表征FC时间序列波动的检验统计量(生物标志物)来完成,然后观察检验统计量的观察值是否落在检验统计量的零分布之外,如果相关性是静态的,它的检验统计量的观察值落在零分布之外。已经有相关研究提出了几种统计检验方法来证明dFC的存在,包括通过FC时间序列的方差证明、通过FC时间序列的傅立叶变换的统计检验方法和非线性检验统计方法等。至关重要的是,构建检验统计量分布的原假设应对应于静态FC。这可能看起来微不足道,但这种分布的构建远非微不足道,这构成了评估dFC的第二个陷阱,即使用了不适当的零假设。
由于大多数dFC度量方法无法通过数学方法推导出零分布,因此需要根据现有的数据对其进行近似。理想情况下,此类替代数据的构建方式使它们与观察到的数据共享所有统计特性,除了它们缺乏想要检验的特性。在dFC的相关文献中已经有若干种方法来近似dFC的零分布。例如通过随机打乱BOLD时间序列,或打乱傅立叶变换中的相位或通过从不同的扫描过程中随机选择BOLD时间序列。这里的陷阱是这两种方法破坏了数据中的sFC(静态功能连接),导致零假设与想要验证的假设出现了偏差。此外,先验的设定是否影响后续统计检验的结果,构建替代数据的一种更合适的方法是将时间序列模型与数据拟合,并通过从模型残差进行引导来近似零分布。另一种更容易使用的方法是对傅立叶变换的相位进行重组,以便保留sFC的分布。作者专注于基于傅立叶变换的替代数据方法。
材料与方法
动态FC的统计评估
作者首先假设已经记录了来自两个体素或感兴趣区域(ROI)的静息状态的 BOLD-fMRI时间序列,如图1A所示,所以现在任务是想确定两个时间序列之间的功能连接是否是动态的,即它在不间断扫描时是否发生变化。尽管FC的概念很广泛,但作者关注的是计算Pearson相关系数的方法,这种方法是静息态fMRI研究中衡量FC使用最广泛的方法,最直接的方法是计算重叠的时间子段中时间序列之间的Pearson相关系数,这得到如图1B中所示的FC的时间序列。请注意对于不同的窗口,相关系数是不同的,特别是相关系数有正有负,分别代表负相关和正相关。尽管在一些研究中,观察到的FC波动被视为存在dFC的证据,但大多数研究需要进行统计检验才能得出这一结论。事实上,如果相关系数的时间序列中观察到的波动可以与如果相关性是静态(即独立于时间)时观察到的波动相区别,则对dFC进行适当的统计检验就可以回答这个问题。

图1. dFC的统计检验方法示意图。

A. 模拟两个不同体素或ROI同时记录的fMRI时间序列。

B. 通过计算fMRI时间序列(最大重叠)的连续60秒内的Pearson相关系数从而获得相关系数的时间序列,分别以相关系数和窗口长度作为自变量和因变量绘制函数图像。虚线表示每个时间窗口中fMRI时间序列得到的相关系数的95%置信区间,水平线表示fMRI时间序列之间的平均相关性。

C. 零假设条件下相关系数的时间序列标准差的概率密度,当标准差为0.24时概率密度达到峰值,并用垂直线标记。

回答这个问题的一种方法是围绕相关系数的时间序列的值来构建置信区间,例如Kang等人(2011)和Hutchison等人(2013b)的研究。如果数据是白噪声高斯过程,则可以计算置信区间,否则它们可以通过重新采样窗口时间序列来近似,这种方法称为bootstrapping方法。图1B中的置信区间是通过这种bootstrapping方法获得的,更具体地说,对于每个时间窗口作者选择(替换)未配对的样本对构建窗口数据的抽样副本并计算其样本Pearson相关系数。通过多次重复,可以构建零分布的近似值,然后计算置信区间。请注意有几个窗口的置信区间不包含平均相关值,即从整个时间序列计算的相关系数,用红线表示。在Kang等人(2011)研究中建议通过置信区间不包含静态相关值的窗口的分数来测量dFC的范围。然而,非零分数本身并不意味着存在dFC。因此如果相关性是静态的,那么度量的值会是多少。作者认为任何用于检测dFC的度量指标都需要解决这个问题,无论是计算相关系数组成的时间序列的方差还是一些非线性度量。
这个问题可以在统计假设检验中形式化,其中零假设对应于静态相关性,而另一种假设对应于动态相关性。如果数据中的相关性是静态的,那么计算度量值的问题相当于确定零假设下度量值的分布问题,并计算从该分布中提取观察值的概率。由于在大多数测量中零分布是未知的,因此使用了一定量的随机数据,例如Chang和Glover等人(2010)以及Zalesky等人(2014)研究。随机数据也称为替代数据,是非线性和非平稳时间序列分析的基础。图1C中表示了从图1A中的fMRI时间序列计算的相关系数的时间序列的标准偏差的零分布。观察到的方差由垂直线表示。因此在零假设下,度量值从该分布中随机选出。由于观测值不包含在分布的2.5%百分位数中,因此不能拒绝原假设,因此没有证据表明dFC存在。
模拟BOLD-fMRI数据
为了衡量检验统计方法在静息态BOLD-fMRI数据中检测dFC的能力,作者首先使用模拟的数据。使用模拟的数据能够更好地控制dFC的时间序列,这是在真实数据下不可以直接观察到的,并且可以研究影响检验统计量的因素,例如它们的强度和时间尺度。为了构建模拟的fMRI时间序列,作者构建方差

的零均值高斯随机过程X=(X1,⋯,XN)和Y=(Y1,⋯,YN),ρnXnYn的Pearson相关系数:

其中E表示期望值,ρn依赖于时刻n,即是动态的。作者假设XnXm对于nm是独立的,对于YnYm也是如此。这意味着XY没有自相关,这通常是预处理BOLD-fMRI数据的情况。作者指定相关系数的时间序列是动态的,从而将动态过程参数化,以便它们的参数对应于可观察的属性,然后可以系统地改变这些属性以评估统计检验的性能。此外作者希望相关系数序列是随机的,而不是对扫描过程的开始阶段锁时,这反映了静息态BOLD波动的持续性质,所以作者将ρn建模为一个非线性变换的变量sn

其中常数

表示XY之间的平均相关性,时间序列S=(S1,⋯, SN)是通过从采样周期为TR=2的连续变量S(t)中采样获得的。变量S(t)建模为以下随机微分方程:

图2. 模拟的相关系数的时间序列,三种情况分别是对模拟的相关时间序列的独立实现。颜色对应于不同的主要时间尺度:τ=180s(蓝色)、τ=90s(红色)和τ=45s(绿色)。在所有情况下,平均相关性设置为零,dFC的强度设置为0.5(η=0.5)。观察时间为300s,相关系数的时间序列下采样到0.5Hz(TR=2s)。
本研究重点是基于滑动窗口相关性的检验统计的性能,对于同时记录的时间序列x=(x1,⋯, xN)和y=(y1,⋯, yN),可求得它们之间的相关系数

。如果xy是平稳随机过程,则X=(X1,⋯,XN)和Y=(Y1,⋯,YN)中的时间序列,说明它们的期望、方差和XnYn之间的相关系数ρn不依赖于n,是ρ的渐近无偏估计量。在这种情况下,作者将XY之间的相关性称为静态相关性。如果ρn取决于n,那么将XY之间的相关性称为动态相关性,最后将ρ=(ρ1,,ρN)称为相关系数的时间序列

通过对xy中连续的时间子段计算相关系数,可以获得相关系数的时间序列的直接估计量。窗口长度为Kρk的滑动窗口估计量是通过在时间窗口的时间序列(xkh,ykh),⋯,(xk+h,yk+h)上计算获得的,其中h=(K−1)/2是窗口的中心。随后,kM个样本上重复移位,产生相关系数的时间序列

,其中L是满足L≤(NK)/M的最大整数。因此,该时间序列的生成取决于两个参数:窗口长度K和滑动步长M接下来作者将描述两个基于相关系数的时间序列的检验统计量:广泛使用的线性度量和最近提出的非线性度量。

零假设和检验统计量
作者认为dFC缺失对应于零假设,即H0:η=0,并且dFC的存在对应于另一种假设,即H1:η>0。
η2是相关系数的时间序列的方差,所以如果能从xy检测出dFC,则对应于拒绝H0。尽管这看起来微不足道,但实际上对dFC的评估是非常重要的,而且如上所述它仍然是静息态fMRI动力学中的一个混乱来源。
作者专注于两个检验统计量及其检测dFC的能力。第一个是静息态fMRI研究中使用最广泛的检验统计量,也是最直接的。它是标准偏差或等效的样本相关序列的方差,表示为κ:
其中

是估计的相关系数的时间序列,

表示的样本均值。接下来作者将用ζ表示第二个检验统计量,其定义如下:设m

的中位数,n1,n2,…,nj

m交叉的样本。因此,

m开始J−1次连续偏移,第j次偏移的长度In和高度Hn定义为In=nj+1-nj

。检验统计量现在定义为:

其中αβ控制偏移长度和高度的相对权重,作者设置α=0.9和β=1用于计算ζ。

检测概率
为了验证上文所述的FC动力学模型允许评估任意检验统计量T检测dFC存在的能力,作者模拟了两个长度为N的fMRI时间序列,其中重复时间为TR,其余参数取值固定

。在上述模型中,零假设对应于H0:η=0,而其他假设对应于H1:η>0。如果FC是静态的,拒绝H0的概率为100α%,其中α通常选择0.05。因此拒绝H0的概率,即犯第一类错误的概率为5%。对于模型参数θ的给定值,用πT(θ)表示拒绝H0的概率。这个概率由下式给出:

为了说明幂函数的用处,作者模拟了两个fMRI时间序列xy,长度为N= 301,重复时间为TR=2,θ=(0,120,0.5)。因此平均相关性为零

,并且相关性是动态的,强度为η=0.5,时间尺度τ=120s,其中强度η表示相关值真实波动的标准偏差。假设使用相关系数的时间序列的标准偏差作为dFC的度量,即T=σρ并取K=20为最大重叠的样本的窗口长度。那么,在模拟的时间序列中实际检测到dFC的概率是多少?图3B显示了H0T的概率密度

(蓝线)。它的第95个百分位数由水平线表示,等于ξα=0.27(α= 0.05)。给定η=0.5,

的概率密度以红色显示。根据πT的定义,检测xy之间的相关性是动态的概率是通过从ξα到∞积分

得到的,得到πT≈0.69。因此在这种情况下,检测到dFC的概率约为70%。因此,如果要扫描20个被试,每个被试10分钟,TR为2秒,作者希望在20个被试中的14个被试里检测到dFC。图3A显示了η=0和η=1之间dFC的幂函数。如果存在弱dFC,即η=0.1,那么检测到它的机会不到10%。

替代数据的构建
假设已经从两个同时记录的BOLD时间序列计算了一个检验统计量T的值,并获得了一个值T*。如何确定观察值是否具有统计显着性?这就意味着需要执行假设检验,其中零假设和另一个假设H0H1分别对应于不存在dFC和存在dFC。通常是否拒绝零假设是通过计算H0T分布中T*的百分位数来完成的。如果T*落在5%的最大范围内,那么则拒绝H0并得出H1为真的结论即得出相关性是动态的结论。
图3. 检测概率。
A.在H0(蓝色)和η=0.5(红色)下的密度。垂直线的位置对应于 H0κ的95%百分位,它是κ的临界值,κ*=0.2702,检测到η>0等于从κ*向上的红色密度下面积的概率等于π=0.6914。参数选择为Ns=301,η=0(蓝色密度)和η=0.5(红色密度),K=20个样本,τ=120s,Ndata=1000。
B.检测概率πκ作为dFC强度η的函数,使用10×103蒙特卡罗样本获得,黑色曲线是通过蒙特卡罗模拟获得,误差线表示每个η值的最小值和最大值。
然而,大多数检验统计量的零分布是未知的,因此作者使用手上的数据进行近似。具体地,通过生成大量适当随机化的数据副本(所谓的替代数据)来近似零分布,然后计算每个副本的检验统计值。通过适当地随机化选择保证每个副本都在H0下构建,因此结果的分布确实是T的未知零分布的近似值。目前有两种替代数据,第一称为约束方法,直接从观察值中构建。第二个称为典型方法,它是通过将数据拟合到模型然后使用该模型生成替代数据来构建的。
Prichard等人(1994)首次引入了适用于动态FC上下文的约束随机化。虽然它可以应用于多变量数据,但在本研究中作者只使用两个变量的情况,因此让x1,⋯, xNy1,⋯, yN同时记录BOLD信号。该方法分别求xy的离散傅立叶变换X1, ⋯, XNY1, ⋯, YN,并与随机相位相乘:

,其中ϕ1,⋯, ϕN是在区间[0, 2π]中均匀分布的独立随机变量的向量。至关重要的是,XY都乘以相同的相位ϕ1, ⋯, ϕN以保持(静态)相关结构。随后作者使用逆离散傅立叶变换应用于X1, ⋯, XN和Y1, ⋯, YN, ⋯,并生成x的随机副本

y的随机副本

记录和预处理
文中使用的猕猴fMRI数据是从三只健康雄性猕猴中收集的(Logothetis等人, 2012)。所有实验数据的使用程序均经地方当局(Regierungsprasidium, Tubingen, Germany)批准,完全符合欧洲共同体(EUVD 86/609/EEC)的实验室动物护理和使用准则。文中实验在被试全身麻醉下进行,瑞芬太尼(0.5-2 g/kg/min)联合速效麻痹剂(mivacurium chloride,5-7 mg/kg/h)。瑞芬太尼是一种超速作用的μ-阿片受体激动剂,因此对神经血管活动没有显着影响,特别是在疼痛基质以外的大脑区域(Wise等人, 2002; Pattinson等人, 2007)。此外存在已知仅轻微影响神经和血管反应的幅度和时间进程(Logothetis等人, 2001年; Goense等人, 2008年)。实验过程中持续监测动物的生理状态,并严格保持在正常范围内,并且通过含2.5%葡萄糖的乳酸林格氏液以10 ml/kg/h的速率输入以预防酸中毒。
作者使用4.7T扫描仪(BioSpec 47/40v,Bruker BioSpin,Ettlingen,Germany)对被试进行测量。定制的正交体积射频线圈用于对大脑深部结构进行成像,通常可以获得整个大脑22个轴向切片。这些切片的BOLD活动是在2秒的时间分辨率下获得的,具有两次梯度回波EPI图像(重复时间=1000毫秒,回波时间=20毫秒,带宽=150 kHz,翻转角= 60度,FOV = 96×96毫米,2毫米切片厚度)。最后作者获得具有相同FOV的T2加权RARE图像。解剖图像后来被变形以匹配EPI图像。为了避免扫描过程中的图像强度不同,每个体素通过除以具有足够强度的所有体素的平均值并乘以被试体素值的1000进行归一化。
T1加权图像被分割成灰质、白质和脑脊液,作者使用12参数仿射变换配准到猕猴MNI图谱(http://www.bic.mni.mcgill.ca/ServicesAtlases/Macaque)。fMRI数据的预处理包括大脑提取(使用结合灰质、白质和脑脊液图像的掩模)、运动校正、T1加权图像的共同配准以及使用先前的猕猴MNI图谱的空间对齐计算12参数仿射变换。在配准到MNI空间后,作者对头部运动参数以及白质和脑脊液信号进行回归,以及使用Chebychev II型滤波器和局部平均以0.05 Hz的截止频率进行低通滤波以获得每个体素的时间序列猕猴MNI空间。随后,根据Paxinos 2008分区,时间序列被分为436个感兴趣区域(ROI)。所有预处理步骤均使用spm工具箱(参见http://www.fil.ion.ucl.ac.uk/spm)中的代码和函数,预处理使用的编译软件是Matlab。
作者还分析了在24名人类被试中收集的静息状态BOLD-fMRI数据,这些数据也用于之前的许多研究(Deco等人, 2013年;Mantini等人, 2013年;Deco等人, 2014年;Ponce-Alvarez等人, 2015年)。志愿者被告知实验程序并签署了书面知情同意书,研究得到了伦理委员会的批准,使用3TMR Philips Achieva扫描仪进行MR数据采集。被试被要求连续注视背景中没有视觉刺激的点,并且在扫描过程中不要移动。功能图像是使用T2*加权回波平面图像(EPI)获得的,使用SENSE成像具有BOLD对比度。EPI包括按升序获取并覆盖整个大脑的32个轴向切片(32个切片,230×230平面内矩阵,TR/TE=2000/35,翻转角= 90°,体素大小=2.875×2.875×3.5mm3)。对于每个被试,进行了两次10分钟的连续扫描。此外,通过MP-RAGE序列(TR/TE=8.1/3.7,体素大小=0.938×0.938×1mm3)收集了用于解剖参考的3D高分辨率T1加权图像。人类fMRI数据的预处理类似于上面描述的猕猴fMRI数据的预处理,唯一的不同是作者将功能图像空间配准到MNI-ICBM人类图谱而不是MNI Macaque Atlas。此外,文中通过局部平均获得ROI时间序列的大脑分割方法是参考的Hagmann等人研究中的方法,大脑分割为66个区域。
动态FC分析
介绍完实验数据的来源和预处理方法,作者正式将前面部分中描述的检测dFC的方法应用于猕猴和人类被试BOLD-fMRI数据。如前所述,目前有3名猕猴被试每个被试各25次的扫描实验,以及24名人类被试每个被试各2次的扫描实验。记录持续时间为5分钟,每次扫描的采样周期为2秒。猕猴被试的数据被分割成436个ROI,而人类被试数据被分割成66个ROI。
对于每次扫描实验中,作者使用两种方法分别获得线性和非线性检验统计数据。在第一种方法中,作者使用两分钟的窗口长度和两秒的间隔获得每对ROI的相关系数的时间序列,然后计算每对ROI的检验统计数据。然后对于每次扫描,作者为每个ROI生成250个替代数据,这样在每一组替代数据中,每对ROI之间的平稳相关性得到了保留。最后将所有对的ROI的值合并在一起,以获得检验统计数据观察值的p值。作者还计算了不同被试观察和替代检验统计值的平均值,并获得了相应的p值。
在第二种方法中,作者计算了ROI的平均检验统计量以提升统计效果。因此对于每次扫描实验,作者使用两分钟的窗口长度和一分钟的步长获得了每对ROI的相关系数的时间序列,然后计算每对ROI的两种检验统计数据,然后取平均值得到每个ROI相对于其他ROI的检验统计数据。因此作者只有436个和66个值与替代数据进行比较,而不是分别为猕猴被试和人类被试构建94830和2145个检验统计量。在第二种方法中,作者为猕猴被试一共生成了10000个替代数据,为人类被试一共生成了1500个替代数据,每个数据都保留了每对ROI之间的平稳相关性,并计算了线性和非线性统计值的代理分布和获得了相应的p值。作者重复这个过程,求每次扫描中观察到的和替代的检验统计值的平均值,并获得相应的p值。
如果您对脑功能及脑结构数据处理感兴趣,欢迎浏览思影科技数据处理课程及服务(可添加微信号siyingyxf18983979082咨询,另思影提供免费文献下载服务,如需要也可添加此微信号入群):
上海:
第四十六届磁共振脑影像基础班(上海,9.1-6)

第十五届磁共振脑影像结构班(上海,9.7-12)

第四届DWI提高班(上海,9.23-28)

第三届磁共振脑网络数据处理提高班(上海,10.12-17)

重庆:

第四十五届磁共振脑影像基础班(重庆,8.26-31)

第十八届脑影像机器学习班(重庆,9.12-17)

北京:

第一届任务态fMRI提高班(北京,8.25-30)

第十七届脑影像机器学习班(北京,9.3-8)

第九届小动物磁共振脑影像数据处理班(9.11-16)

第二十三届磁共振脑网络数据处理班(北京,9.17-22)

第二届磁共振脑网络数据处理提高班(北京,9.23-28)

第四十七届磁共振脑影像基础班(北京,10.22-27)

南京:

第二十届DWI数据处理班(南京,10.21-26)

数据处理业务介绍:

思影科技功能磁共振(fMRI)数据处理业务
思影科技弥散加权成像(DWI/dMRI)数据处理
思影科技脑结构磁共振成像数据处理业务(T1)
思影科技啮齿类动物(大小鼠)神经影像数据处理业务
思影数据处理业务三:ASL数据处理
思影科技脑影像机器学习数据处理业务介绍
招聘及产品:
招聘:脑影像数据处理工程师(上海,北京)
BIOSEMI脑电系统介绍
目镜式功能磁共振刺激系统介绍
结果
对时间尺度的依赖性
为了评估线性检验统计量κ的检测能力是否取决于ρt的动态性和强度η,作者计算了不同相关系数的时间序列尺度τ值和dFC的强度η的检测概率πκ之间的关系。其余参数固定为

、TR=2s、N = 300个样本和M = 1,定义τ的范围从20s到600s,间隔为10s,定义η在0.1和1之间,间隔为0.1,窗口长度K在与τ相同的范围内变化。

图4A为概率

作为τη的函数可视化结果,每个图对应每个图上方显示的η值。从图中作者发现了几个有趣的现象:

首先虽然(绝对)概率随着η的增加而增加,但相对概率保持大致稳定。例如无论强度η如何变化,选择K <τ总是比选择K>τ更好地检测出dFC。
其次如果τ很小,比如τ<30s,那么dFC永远不会被检测到,因为无论强度η多强,检测概率πκα(设置为0.05)。
第三,检测到dFC的最高概率对应的K的最大值,用Kopt表示,大致是相关系数的时间序列时间尺度的三分之一:

然而,从这些图中还可以发现,相关系数的时间序列的波动很慢(即τ数值大),

,因此K的选择应仅被视为经验法则,一方面,为了减小检验统计量的方差,K应该尽可能大;另一方面,在给定采样频率的情况下,规定≈τ/3仅具有理论意义,因为在实际实验中,(静息状态)BOLD-fMRI时间序列之间的相关系数波动的时间尺度是未知的。

然而,通过在所有时间尺度上计算检验统计量的平均值仍然可以说明一些结论。在图4B中,作者绘制了一组检测概率曲线,它们表示在所有时间尺度上计算得到的平均概率,并以窗口长度为自变量绘制出函数图像。十条曲线对应于dFC强度的十个值,位置较高的曲线对应于较高的η值。这些曲线表明,检测概率与dFC的强度无关,并且最佳窗口长度约为50秒。有趣的是,这种窗口长度在实验研究中很常见的。

对扫描持续时间和扫描实验次数的依赖性

为了评估检测概率是否取决于扫描持续时间,作者将其计算为样本数N的函数,以150为间隔在150到1500之间变化,TR为2秒,扫描持续时间的间隔从5分钟到50分钟不等。模型参数设置为

和滑动窗口参数为(K, M)=(30, 1)。

实验结果如图4C所示,可以发现检测概率和扫描持续时间之间的函数近似线性增加。例如5分钟前的静息状态扫描中检测到dFC的概率约为15%,5分钟后增加到约50%。

对于人类被试,长时间的扫描实验是不切实际的,提高检测概率的另一种方法是在多个扫描实验期间测量BOLD的波动程度,然后对扫描实验的检验统计数据求平均。为了评估检测概率是否取决于扫描实验的次数,作者计算了扫描实验的次数(每次5分钟)和检测概率的函数,扫描实验的次数的变化范围是1到25。结果如图4C(右)所示,该图显示检测概率随扫描实验次数的增加而近似线性增加。

作者从理论上对这些实验结果进行理解,并举了一个简单的例子,即线性检验统计量κ的分布可以通过解析来计算。假设κ是使用非重叠窗口(M = K)计算的,并且总体相关系数的时间序列ρt在每个窗口内是恒定的,比如ρt =ρn,在第n个窗口内,pn服从方差为η2的正态分布。虽然这种形式的相关系数的时间序列不是dFC动力学模型的特例,但η2具有相同的解释,即它是总体相关系数的时间序列的方差。在这些假设下,如果KN很大并且N远大于K,则κ的期望和方差近似由下式给出:

第一个公式表明κ2的期望值近似等于η2,即dFC的强度。第二个公式表明κ2的方差收敛到零,这意味着在不考虑K的情况下,在H0 :η=0和H1:η>0下κ2密度重叠收敛到零,即检测概率πκ收敛到1。在这个例子中,dFC的存在因此总是可以被检测到,前提是观察时间足够长。然而第二个公式也表明,如果dFC很小(η2≈0),则需要更多的样本来检验它。
线性与非线性度量
上文中作者发现在持续时间为10分钟的单次扫描实验中检测到dFC的概率比较低,这种低概率可能是由于使用的检验统计量对FC的变化相当不敏感所致,或者原因可能在于相关系数的时间序列本身,即相关系数的方差具有不确定性,这意味着从相关系数的时间序列中得出的任何检验统计量的敏感性都很低。为了证明后者,作者将(线性)检验统计量κ的检测概率与最近提出的非线性检验统计量ζ进行了比较,并将其表示为ξ
图4. 基于窗口相关的检测概率。
A. 窗口长度K和时间尺度τ与检测概率之间的函数,这十个数字对应于不同的η值(dFC的强度),如上图所示。η在0.1到1的范围内以0.1为间隔递增,图中每个白色直线是相同的,表示线(K,3K)。添加它是为了说明选择满足最佳窗口长度的经验法则,粗略地说,它等于相关时间大小的三分之一。
B. 检测概率在所有相关系数的时间序列的时间尺度上取平均值并绘制为滑动窗口长度的函数,十条曲线对应于η的十个值(较高的曲线对应于较高的η值)。
C. 作为观察时间的函数的检测概率,即扫描过程的持续时间(左)和扫描实验次数的函数关系(右)。在所有情况下,模型参数都设置为

和滑动窗口的参数为(K, M)=(30,1)。此外,在计算图中右侧的检验概率时,观察时间设置为T=300 s。图中的两个误差条表示10次模拟的最小和最大概率。

D. 表示的是线性(左)和非线性(右)检验统计的检测概率。时间尺度的范围为20至600秒,间隔为1秒,dFC的强度范围为0.1至1,间隔为0.1。窗口长度设置为时间尺度的三分之一,并且窗口一次滑动通过模拟的fMRI时间序列一个样本。平均相关性设置为

,重复时间设置为TR=2s。

在以下的模拟实验中,作者将观察时间设置为10分钟,并以10秒的间隔在20到600秒之间改变时间尺度τ。dFC的强度由η衡量,以0.1的间隔在0.1和1之间变化。计算相关系数的时间序列的窗口长度设置为时间尺度的三分之一,因为该窗口长度被发现是最佳的。窗口在模拟的fMRI时间序列中每次滑动一个样本(M = 1)。此外,模拟的fMRI时间序列之间的平均相关性设置为零(

),TR =2s。

线性和非线性检验统计量的检测概率如图4D所示。作者进行如下两个观察:
首先,与非线性检验统计量相比,线性检验统计量对较低的dFC强度产生较高的检测概率。
其次,线性检验统计量对相关系数的时间序列中的时间尺度不敏感,但非线性检验统计量的高检测概率局限于较短的时间尺度,后者原因是非线性检验统计量由相关系数的时间序列中偏移的宽度和高度的乘积组成,如果时间尺度很大,可能只有一到两个偏移。
在任何情况下,由于非线性检验统计量在dFC检测方面并没有做得更好,所以作者怀疑问题出在相关系数的时间序列本身。
低通滤波的影响

静息态BOLD-fMRI时间序列的预处理通常包括低通滤波以去除伪影并选择合适的频率。然而,这种滤波会直接或通过在BOLD时间序列本身中引入自相关来改变整个相关系数的时间序列。为了深入了解这两种效应的相对强度,作者用ρt和滤波器的归一化脉冲响应ht来表示低通滤波时间序列的总体相关系数的时间序列,用ρfilt表示。

其中⊗表示卷积算子,如果一个BOLD时间序列被一个具有(归一化)脉冲响应h(t)的滤波器过滤,即与h(t)卷积,则相关系数的时间序列将被脉冲响应过滤。在频域中该方程的形式为:
其中PfiltH分别表示ρfilth的傅立叶变换后的结果,ω表示(角)频率。请注意,对于零相位滤波器,比如频率响应H0,它认为:
因此H是实数不会导致ρ的相位谱失真。由于频率响应HH具有与H相似的截止频率,如果1/τ<fc,则ρfilt(t)≈ρ(t)。也就是说如果(总体)相关系数的时间序列ρ中的主导周期小于零相位低通滤波器(具有陡峭滚降)的截止频率fc,则(总体)相关系数的时间序列的过滤后的时间序列不会因过滤而失真。当1/τ>fc时,滤波会衰减ρ的波动。
为了说明这一点,作者再一次进行了本节第一段中描述的模拟实验,其中模拟了η=0.5的BOLD时间序列,并用截止频率fc = 0.05Hz对它们进行了低通滤波。图5A显示了人类被试相关系数的时间序列(黑色轨迹)及其过滤版本(绿色轨迹),如从上述等式计算得到,显然过滤不影响时间序列的相关。图5B和图5C显示了原始和过滤的BOLD时间序列,过滤后的时间序列中可见的自相关是由低通滤波引起的,降低了dFC的检验概率。
图5. 低通滤波的影响。
A. 总体相关示例(黑色线条,Nsample=1501,η=0.5,τ=600),滤波版本(绿色线条)和截止频率(fc=0.05 Hz)。
B. 模拟fMRI时间序列。
C. 相同的时间序列,但使用截止频率为fcHz进行低通滤波。
因此,对于通常用于过滤fMRI时间序列(0.05<fc)的截止频率,dFC的波动很可能不受影响。也就是说,可以通过滑动窗口法检测到的dFC的波动,因为当fc低至0.05Hz时,时间尺度为20s或更短的波动将被衰减。作者还得出结论,滑动窗口的检测能力下降,低通滤波fMRI时间序列上的估计量κ完全是由于低通滤波在时间序列中产生自相关这一事实。值得注意的是,当以相对于窗口长度的适当截止频率执行高通滤波时,可以最大限度地减少相关系数的时间序列中的虚假波动。
全麻猕猴BOLD-fMRI的应用

首先,作者计算了每个ROI对的线性和非线性检验统计量的值(参见动态FC分析部分)。图6A显示了所有ROI对的线性(顶部)和非线性(底部)检验统计量的值,对于每次单独的扫描实验,以及多次实验的平均值。对于每次扫描,根据图6B中所示的相应的p值对ROI进行升序排列。请注意,在每次单独的扫描实验中,对检验统计值的可变性大于多次实验的平均值同时,与任何单次扫描实验相比,平均情况下有更多具有较低p值的ROI对。重要的是,多次实验平均值p值的分布会被转移到较小的值,这表明多次扫描求平均会提升统计效果。尽管多次扫描求平均提高了统计的效果,但没有一次实验平均值超过5%,即显著性阈值。在多个测试中,保守度较低的错误发现率(FDR)校正也没有产生任何ROI对具有统计学意义的动态FC(见图6C和D)

为了提升统计效果,作者计算了ROI的检验统计值的平均值(参见dFC分析部分)。图7A表示线性(左)和非线性(右)检验统计数据的ROI平均值,图7B显示了相应的FDR校正值。图7C为显著的动态连接的数量,包括扫描实验次数以及ROI检验统计数据的平均值的数量。作者补充了一些备注,首先,尽管在线性和非线性检验统计数据中显著的动态连接的比例几乎为零,但在跨扫描实验求平均时,显著的动态连接的比例要大得多。这完全符合作者在“Section dependency on scanning duration and session number”中的模拟结果,该结果预测了在整个扫描实验中统计效果的提升。在线性检验统计中,有213个ROI(左半球109个,右半球104个)具有显著的动态连接非线性检验统计中有73个ROI(左半球38个,右半球35个)具有显著的动态连接,该结果体现了对dFC进行预测时,线性检验统计量比非线性测试统计量更敏感。
图6. 全身麻醉下猕猴大脑成对ROU之间的动态功能相关性。
A. 94830对ROI的检验统计值,分别为25个扫描实验中获得,并对所有扫描实验的结果求平均值(最后一列)。
B. 对应于A中的检验统计量的p值,使用Bonferroni校正针对多重比较进行校正,未校正的p值是通过将检验统计量的观测值与250对替代数据的观测值进行比较而得到的。在A和B中,每对ROI按(Bonferroni校正)p值升序排列。
C. 每对ROI之间的检验统计值,所有25次扫描实验结果的平均值。
D. 对应于C中扫描实验的平均检验统计值的p值,使用错误发现率多次比较并进行校正。

图7. 全身麻醉下猕猴的ROI的平均动态功能相关性。

A. 436个ROI中线性(左)和非线性(右)检验统计量的平均值,以及25个扫描实验的平均值(最终列)。每次扫描实验中每对ROI的相关系数的时间序列是使用2分钟的窗口长度和1分钟的步长获得的。

B. 对应于A中的检验统计值的p值,使用错误发现率对多次比较进行校正。未经校正的p值是通过将检验统计量的观测值与10000个阶段随机替代时间序列的观测值进行比较得到的。

C. 为每次扫描实验计算的具有统计显著性(95%置信水平)的检验统计值以及多次扫描实验求得的平均ROI数目。

图8A表示线性(左图)和非线性(右图)检验统计量平均值的可视化结果。有趣的是,这些值在半球上分布相当对称,尤其是线性统计。图8B表示观察值(品红色),以及从替代数据(灰色)中获得的10000个值。黑线对应于95%的置信水平。图8C表示相应的未校正和FDR校正的p值。图8D中,具有显著的动态连接的ROI用红色表示。再次注意,在线性统计的情况下,半球的分布相当对称。另一个有趣的观察结果是,在线性统计的213个显著的ROI中,有一定数量的脑区是属于皮质下的脑区,特别是小脑似乎具有高度的动态连接性。此外,在线性统计的情况下,属于默认模式网络(DMN)猕猴的9个皮质ROI中有7个是动态的。具体来说,这些ROI是左前扣带皮层的区域24和32、左背外侧前额叶皮层的区域9和46d、弓状沟的区域8b、后扣带皮层的区域23、顶叶区域、颞顶枕叶的尾部区,以及顶叶区的中间部分17。此外几乎所有属于海马-内嗅皮层区域的区域(CA1CA2CA3CA4、旁和前突以及内嗅皮层的尾部、嗅觉、嘴部和中间部分)与其他区域均表现出显著的平均动态相关性。
图8. 全身麻醉下猕猴的扫描实验的平均动态功能相关性。
A. 在猕猴MNI模板的水平切片上显示颜色编码的线性(左)和非线性(右)ROI检验统计的平均值。
B. 在数据中观察到的线性(左)和非线性(右)检验统计量的值(品红)和来自10000个阶段随机时间序列的检验统计量(灰色)。黑线表示每个ROI的替代分布的95%百分比。
C. 线性(左)和非线性(右)检验统计量的nn值,未校正(蓝色)和校正(绿色)使用错误发现率进行多次比较。
D. 显示ROI空间分布的可视化结果,显示统计显著值ROI平均动态相关性(红色),使用线性(左)和非线性(右)检验统计量。
在人体静息态BOLD-fMRI中的应用
接下来,作者试图识别人类被试在静息态下的动态功能连接。因此,首先我们计算了在24个人类被试的2次扫描实验中的每2145个ROI对的线性和非线性检验统计量。与猕猴被试的情况类似,与单次扫描实验相比,平均情况下检验统计量的变异性较低(A)ROI对的数量更多,p值更低(B)。因此,该分析还证实了多次扫描实验的平均值提升了统计性能。然而,使用FDR校正后的平均p值均未超过5%显著性阈值,这与作者在猕猴数据中的发现一致。
图9表示平均线性和非线性检验统计量对应的ROI的可视化结果。值得指出的是它们在两个半球之间的对称性:在线性和非线性检验统计量下,左右半球值之间的(Pearson)相关性分别为0.79和0.48。这些观察结果与对猕猴数据做的检验一致。
图9. 人类静息态BOLD-fMRI的平均动态功能相关性。
该图显示了线性(左)和非线性(右)检验统计量的平均值对应的ROI的可视化结果。
讨论
本文研究的第一个目标是证明在静息态BOLD-fMRI数据中的dFC分析中进行适当的统计检验的重要性,并描述如何使用适当的替代数据进行此类检验。我们注意到这种检验方法也适用于电生理学(EEG/MEG)和与任务态相关的数据,作者以此试图解释为什么不能直接将通过滑动窗口法得到的相关系数的时间序列中存在的波动作为dFC存在的证据。第二个目标是强调选择适当的零假设的重要性,这是通过将数据随机化以产生替代数据的方式完成的。最后一个目标是评估基于滑动窗口的检验统计方法在检测静息状态BOLD-fMRI数据中的dFC方面的性能。作者的主要结论是:
在单次静息态扫描过程中,无论使用的检验统计数据是线性还是非线性,dFC的检验概率都很低,原因在于滑动窗口法本身。然而通过模拟,作者还发现,可以通过对多次扫描实验的结果求平均以提高检验概率。作者在猕猴数据上也发现和强调了这些结论。
具体来说,作者试图通过对3只全身麻醉的猕猴的25组BOLD-fMRI数据,使用线性和非线性检验统计量,识别每一对ROI的dFC。重要的是作者没有发现成对dFC存在的证据,无论是在单次的扫描实验中,还是通过多次扫描的检验统计值的平均值。然而,通过计算ROI的平均检验统计量来减少比较的次数,作者发现了分布式dFC网络存在的有力证据。所识别的网络在大脑半球是对称的,包括几个皮层以及皮层下区域。特别是几乎所有属于海马体和内嗅皮层的区域都是这个网络的一部分。这一发现与Logothetis等人在该数据集中记录的每次扫描实验期间海马区域的结果相一致(Logothetis等人,2012)。他们观察到,许多皮层和皮层下区域分别表现出短暂的激活和弱激活。因此,这些事件和反应可能促成了与皮层区域的短暂连接。
这项研究在几个方面受到限制。首先,模拟的 BOLD-fMRI 数据是使用功能相关性动态的特定模型构建的。具体来说,作者假设相关系数的时间序列表现为噪声驱动的阻尼谐波振荡。事实上,其他模型已被用于模拟相关系数的时间序列的动态性,例如向量自回归 (VAR) 模型、确定性振荡和瞬变。在某种程度上,这些选择是任意的,因为潜在的相关动态是未知的。除了一些基本要求,例如信号在模拟开始时没有进行时间锁定,仍存在一些不足。例如,作者对检测概率与动态相关性的强度和时间尺度是否有关感兴趣。其次,作者仅限于研究两个检验统计量的性能,而在早期研究中已经使用了更多检验统计量,结果证明这些检验统计量在检测dFC方面表现更好。
本研究中使用代理数据的一个缺点是,Chang等人使用的受约束的代理数据也受到影响,即生成代理数据的原假设比作者感兴趣的原假设更具体,即没有关于dFC的假设。更具体地说,替代原假设等同于数据是静止的。虽然这意味着没有dFC,但后者并不意味着数据的平稳性。因此,当使用两种替代数据中的任何一种检测dFC时,作者可以得出的唯一结论是数据是非平稳的。这是检测dFC的一个相当严重的障碍,由于dFC可能影响个体 BOLD-fMRI 时间序列的动态性,使得相关系数的时间序列的动态性和BOLD-fMRI 时间序列的动态性变得更加复杂。另一种可能是BOLD-fMRI过程的信噪比在扫描期间发生变化,这将影响观察到的相关值,尽管真正的相关性也可能是静态的。因此,作者检测出dFC的检测绝非易事。
几项研究报告了BOLD-fMRI的波动与来自电生理记录的测量值以及行为测量值之间的相关性。例如在Tagliazucchi等人发现功能连接的波动与不同频段的 EEG功率相关。在Di等人研究中,发现几个静息态网络内功能连接的波动与相应网络的活动水平相关。然而,这种相关性本身的存在不能作为功能相关性本身是动态的证据,同样适用于报告条件之间dFC差异的研究。在作者看来,这些研究的有趣之处在于,它们表明BOLD-fMRI时间序列之间的协调与电生理过程共变,从而扩展了研究电生理相关性的研究,例如EEG α功率,单个BOLD-fMRI 时间序列。因此,作者并不声称在样本dFC检验统计中观察到的波动是无意义的。然而,作者声称的是如果没有适当的统计检验,这种波动不能解释为dFC,因为dFC具有特定的统计意义。这引出了一个重要的问题,即根据(时间相关的)FC对静息态BOLD-fMRI数据的描述是否是最自然或最有意义的。
一个可能导致混淆并因此应仔细解决的重要问题如下:假设有一种强大的统计方法可供使用,当应用于特定数据集时,无法拒绝静态FC的原假设。这如何与神经相互作用的高度动态性质一致?事实上,在静止状态下,网络的动力学极其复杂性,具有非线性和动态性。为了调和两者,作者指出数据的统计描述本身并不能说明底层系统的动态特性。相反,统计描述通常是在缺乏明确的确定性描述的情况下采用的,因此通常反映了对系统内部工作原理的缺乏了解。现在的关键是,只有在数据的统计描述中,平稳性(或缺乏平稳性)的概念才有意义,因为它是一个固有的统计概念(这意味着随机过程的转移概率密度与时间无关)。事实上,正是复杂性和大量自由度本身使作者能够首先采用系统动力学的统计描述。因此,在作者的研究中并不声称没有dFC意味着系统本身的平稳性,因为这个概念在统计框架之外没有任何意义。作者研究的重点是:如果采用数据的统计描述并询问FC是否是动态的,则需要使用适当的统计工具来回答。在下一段中,作者认为对静息状态 BOLD-fMRI 数据的纯粹统计描述可能不是最具信息量的描述。
总结:
总体来说,本文研究探究了dFC是否是将静息状态 BOLD 波动进行动态概念化的最自然的方式。作者首先使用了模拟数据来进行测试,然后在24名人类被试的48次扫描实验中寻找动态FC,并于与猕猴被试的情况进行了比较。
结果发现,无论是在单次扫描实验中还是通过多次扫描实验平均,作者都没有发现任何成对dFC存在的证据。然而,当作者对所有48次扫描实验求平均时,计算ROI的平均检验统计值得出了几个区域非平稳性的证据。非线性检验统计量对大约三分之一的区域产生显著值,而线性检验统计量对所有区域都产生显著值。这表明多次扫描提高了对dFC进行检测的统计检验力,同时,作者的模拟实验也表明,对dFC的检测概率随着求平均的实验次数的增加而增加。两种统计检验数据的值显示出半球的对称性可能也一定程度上说明了dFC的存在性。
本文的结果可能反映了FC波动包含一系列瞬时相关的激活模式的假设,但从实际来看,短时的单次的FC扫描可能产生不正确的dFC的构建,但多次扫描的平均可以加强统计检验力。在人脑中的这种共激活模式可能是自发BOLD活动的组成部分,而dFC则是这些活动的反映,并且共激活模式可能是神经和认知事件的血流动力学的相关因素。
(0)

相关推荐