上QQ阅读APP看书,第一时间看更新
第六节 颅脑结构和功能成像分析方法
一、脑形态学分析方法
(一)基于体素的形态学测量
基于体素的形态学测量(voxel-based morphometry,VBM)是一种在体素水平对脑灰(白)质成分进行定量的技术,它是一种数据驱动的方法,能定量评估脑局部灰(白)质密度或体积的改变,从而精确地表征脑组织形态学变化。其基本思想是使用图像分割技术提取脑灰质和白质成分并配准到同一坐标空间,然后对每个体素的灰(白)质密度或体积进行统计学分析,从而揭示脑认知或疾病相关的局部脑结构变化。
在VBM中,灰质密度(或称浓度)是指灰质成分在某个像素内的概率大小,是统计学上的概念,白质亦然。而灰质量(或称体积)是指调制后的灰质密度,它等于灰质密度乘以空间标准化矩阵的雅可比行列式。典型分割后的脑组织成分如图1-6-1。
图1-6-1 脑组织分割成分图示例
A.原始高分辨三维T 1WI;B.灰质成分图;C.白质成分图;D.脑脊液成分图
VBM的基本处理步骤包括组织分割、空间配准、空间平滑、统计分析四个步骤。其中前两步是VBM的核心技术,后两步是脑影像统计的通用步骤。
1.脑组织分割
脑组织分割是利用MRI图像中不同脑组织的信号差异,将脑体素分成若干个独立的成分,包括脑灰质、脑白质和脑积液等。所用MRI图像通常为高分辨三维T 1WI;也有研究尝试利用两种或者两种以上的不同对比度MRI图像进行组织分割。原始T 1WI的信号对比不一定完全反映了组织T 1值的差异,其他系统因素也会对它有贡献。例如,由于扫描时不同脑区离线圈中心的距离及深度不同,导致相同脑组织成分(如脑白质)在不同位置受B1场激发的程度不一样,从而出现MRI信号差异。为了保证分割的准确性,需要对这类MRI信号进行校正,我们称之为B1场偏倚校正。
脑组织分割另一个需要考虑的问题是不同脑组织类型的信号分布是高度重叠的,其中一个主要原因是在现有空间分辨率条件下(1mm 3),一个体素可能包含多种组织成分,特别是在组织交界区(灰-白质、灰质-脑脊液、白质-脑脊液等),这种现象我们称为部分容积效应。因此,和传统硬聚类方法不同(如K-means聚类,每一个体素只能归属到某一类别中,即其值为0或者1),VBM的分割需要考虑部分容积效应,因此一般采用软聚类方法,得到的组织成分图中,每个体素的值为从0到1的连续概率分布,反映了该体素内该组织成分占有的比例,称为组织密度。在组织分割过程中,为了提高分割的准确性,可以引入额外的组织概率图谱作为组织分割的先验信息。但是使用该方法的前提是组织概率图谱与待分割的组织成分空间分布要相似,如果使用不当,反而导致分割错误。例如,成年人的组织概率图谱就不适合作为婴幼儿脑组织分割的先验模板。因此,比较合理的做法是用户自己定制适合待研究人群的组织概率图谱用于组织分割。
2.空间配准
由于扫描时被试之间头颅定位的差异和脑结构的变异,个体之间初始脑组织成分图在空间坐标上没有对应关系。为了使被试之间基于体素的脑结构比较具有可行性,需要把每个被试的脑组织配准到相同的标准空间。空间配准可以分为线性配准和非线性配准两种方法,它们在VBM中均用到。线性配准又分为刚性配准和仿射变换两种类型,其中刚性配准是指通过平移和旋转(每种操作又分别对应 x、 y、 z三个方向,故刚性配准共包含6个参数)把个体脑对齐到大体一致的空间,该方法并不改变脑的大小和形状。该方法常用于对同一被试在不同时间或用不同序列采集到的MRI图像的配准,例如,同一被试纵向采集的T 1WI(模态内配准),或者T 2WI与T 1WI之间的对齐(模态间配准)。仿射变换是在刚性配准6个参数基础上,引入了缩放和切变两类空间变换,包含12个参数。仿射变换会改变脑的大小和大体形状。仿射变换常用于不同被试间脑MRI图像的配准,或者把个体脑配准到标准模板空间。尽管线性配准能校正个体之间脑总体大小和位置的差异,但是它们不能校正脑局部结构大小和形状的差异,这导致基于体素的比较也不是特别准确。因此,VBM引入了非线性配准方法,它能校正脑局部位置、大小和形状的个体差异。非线性配准需要以线性配准的结果为基础,经历了传统一步非线性配准,联合分割配准法以及基于指数李代数的微分同胚解剖配准(diffeomorphic anatomical registration through exponential lie algebra,DARTEL)等阶段,目前已经能非常准确地对个体被试的主要沟回进行对齐。然而,如果个体之间脑结构变异非常大(例如,有些个体的局部脑回缺如),即使用DARTEL等高阶非线性配准方法,也不可能把两个脑的结构匹配得完全一样(图1-6-2)。
不管是线性配准还是非线性配准过程,均会产生一个变形场,或者称为变换矩阵,它定义了局部脑结构对齐到参考模板需要的变换(如平移、旋转、拉伸、缩放等)定量参数。可以根据变形场计算每个体素的雅可比行列式(Jacobian determinants,JD),该值代表了该体素进行空间变换后的体积变化。我们可以把雅可比行列式直接作为定量指标进行基于体素的分析,这称为基于张量的形态学测量(tensor based morphometry,TBM)或基于变形的形态学测量(deformationbased morphometry,DBM)。雅可比行列式也可以用来校正VBM分析中灰(白)质成分在配准过程中的体积变化,即用组织密度图乘以雅可比行列式,校正后的灰(白)质成分称为灰(白)质体积,这个过程称为调制。雅可比行列式既包含线性配准信息(代表了全脑体积的变化),也包含了非线性配准信息(代表了脑局部体积的变化)。在调制过程中,如果我们计算全部变形场的雅可比行列式,校正后的灰(白)质成分称为绝对灰(白)质体积;如果只计算非线性配准变形场的雅可比行列式,校正后的体积称为相对灰(白)质体积,它代表校正了个体脑体积变异的局部灰(白)质体积,是目前常用的调制方法。
图1-6-2 VBM分析流程图
3.空间平滑
空间平滑通过引入高斯平滑核对每个体素的灰(白)质含量求卷积,平滑后每个体素的值代表了它自身以及相邻体素的加权平均值。在统计分析之前需要对分割配准后的灰(白)质成分进行平滑,原因包括三方面:①大多数VBM研究都用参数检验。参数检验实施的条件是因变量残差服从高斯分布。根据中心极限定理,对图像进行空间平滑能提高数据分布的正态性,因此在使用参数统计的时候更加可靠。②由于目前的配准方法并不能把不同被试脑结构的细微差异完全对齐,因此平滑能够部分补偿这种细小的空间失配。③根据匹配过滤原理,平滑后对与平滑核大小相似的团块的探测更加敏感。实践中一般常用的平滑核大小为8mm左右,因此,由噪声导致孤立的散在体素差异不容易被探测到。
4.基于一般线性模型的统计分析
VBM等脑影像统计分析常使用一般线性模型(general linear model,GLM)。GLM假设因变量 y是多个自变量 x按照一定权重 β的线性组合,通过最小二乘法等算法求解参数 β,然后对这些参数进行统计推断,公式如下:
其中 y是因变量; x是自变量,如分组和临床认知评分; β是需求解的未知参数组合,每个 β对应一个自变量,代表该自变量的效应量; ε是误差项,即模型中不能被自变量的效应解释的 y的残差。
GLM模型的自变量 x可以是连续变量,也可以是分组变量;因变量 y必须为连续变量。因此,GLM具有普适性,能替代大多数单变量统计分析,包括各种独立测量或者重复测量的student- t检验、单因素或多因素方差分析、回归分析、协方差分析等。因此在大多数脑影像分析软件中,均用GLM进行统计建模。
在脑影像分析中(包括VBM和其他基于体素的分析方法),影像指标图[例如灰(白)质体积]即为因变量 y,通过输入自变量(例如分组信息),计算得到的 β值的空间分布图称为参数图。脑影像分析中的GLM模型如图1-6-3所示:
GLM统计分析包括统计建模、参数估计、统计假设和统计推断四个步骤。
(1)统计建模:
也叫矩阵设计或因素设计,即确定GLM模型的因变量(如VBM中的灰质体积)和自变量(如分组信息)。设计矩阵是所有自变量构成的矩阵,其中每行对应一个样本,每列对应一个自变量。在因素分析中,自变量的定义有很多方法(图1-6-4)。例如,把每个独立的组定义为一个自变量,每列自变量向量中,对应的元素属于该组就赋值为1,不属于该组赋值为0,这种因素设计称为组均值(cell mean)模型,在该模型中,每列自变量对应的 β代表了该组的均值。该设计见于SPM和FSL等统计工具中。
图1-6-3 脑影像分析与一般线性模型
图1-6-4 GLM模型常用因素设计方法
左图为组均值模型,该模型每个独立的组定义为一个自变量,每列自变量中,对应的样本属于该组就赋值为1,不属于该组赋值为0,共 N列,在该模型中,每列自变量对应的 β代表了该组的均值;右图为因子效应(factor effect)模型,该模型首先定义总均值为独立一个自变量,在余下 N-1( N为独立组数)个自变量中,需要设置一个参考组,参考组所对应的元素为固定值[其中子模型效应编码(effect coding)模型为-1,虚拟编码(dummy coding)模型为0],其余每列自变量向量元素中,属于该组就赋值为1,不属于该组赋值为0
(2)参数估计:
是在统计建模基础上,采用数学算法(如最小二乘法)求解未知参数 β和残差 ε的过程,其中每个 β元素代表了与之对应的自变量的效应量,而 ε代表了因变量 y不能被全部自变量的效应总和解释的剩余残差。
(3)统计假设:
即构建contrast,其本质是建立零假设的过程。在建立GLM模型并估计出每列自变量的参数 β后,我们可能对这些自变量本身的效应量 β不感兴趣,而是想知道这些自变量的效应量之间的关系,这就需要构建相应的 β值组合。例如,想研究第一组和第二组的灰质体积均值差是否有统计学意义,对cell mean模型设计而言,其contrast为[1-1],即用第一个 β值(代表组1的灰质体积均值)减去第二个 β值(代表组2的灰质体积均值)。需要注意的是,contrast的定义和设计矩阵是对应的。相同的统计假设,不同的设计矩阵,其contrast可能是不相同的(图1-6-5)。
(4)统计推断:
计算统计概率 P值并作出是否显著性的判断。其核心步骤为计算统计概率 P值和确定显著性阈值。根据 P值的计算方法,可以分为参数检验和非参数检验。参数检验需要确保因变量 y服从正态分布。如果这一前提条件成立,就可根据 t分布或者 F分布的特征,输入β组合、残差项以及自由度等,计算出对应的 t值或 F值,后者进一步转化成相应的 P值。然而,近年来有研究者提出脑影像指标[包括灰(白)质体积]在研究人群中并不完全服从正态分布,使用参数检验的合理性受到了质疑。因此,有学者提出了基于GLM的非参数统计模型,例如,置换检验。该方法的基本原理是随机打乱样本标签(例如,组别信息,这一步称为置换),对每次置换后的数据进行相同的GLM建模和参数估计。置换后分组标签是随机的,因此计算得到的效应量(如两组均值的差异)理论上为零;通过计算机进行上千次以上的迭代运算,便可构建出真实的零分布;最后比较真实效应相对于零分布的概率,从而获得统计概率 P值(图1-6-6)。该方法不需要确保因变量 y的残差服从正态分布,理论上更加合理。不过该统计方法耗时长,需要高性能计算机支持。
图1-6-5 GLM模型contrast定义举例
统计推断另外一个核心步骤是确定显著性阈值。基于体素的影像分析方法对每个体素都进行一次独立的统计学分析,分析整个人脑需要数万次甚至数十万次统计,如果不校正多重比较,犯Ⅰ型错误(假阳性)的概率非常大。因此,有必要进行多重比较校正,设置比较严格的显著性阈值,控制总体假阳性率。目前主流的多重比较校正方法分为体素水平和团块水平校正。体素水平常用的统计学方法为总体Ⅰ型错误率(family-wise error,FWE)校正,该方法采用随机场理论评估全脑内有效的独立比较次数(而不是实际比较次数或体素数),然后控制每次独立比较中出现假阳性错误的概率。该方法既控制了假阳性错误水平,同时又弥补了传统Bonferroni校正太过严格导致过多的假阴性结果。团块水平多重比较校正方法是先通过评估随机噪声导致的假阳性体素随机聚合成一个团块的概率(团块大小的零分布),然后比较实际结果中的团块大小相对随机分布出现的概率,从而判断该大小的团块是否由随机噪声导致。该方法需要首先确定体素水平的阈值来定义团块(体素水平 P值一般要低于0.001)。
统计推断的阳性结果可以叠加到容积图像上进行呈现,也可以投射到皮层表面呈现。
图1-6-6 置换检验基本原理及置换策略
A.置换检验基本原理(以双样本 t检验为例),置换检验的核心思想是通过打乱样本构建真实的零分布,进而直接计算 P值;B.置换策略,包括:①顺序置换,打乱样本次序,该置换适合组间比较或者相关等分析,要求组间方差齐,理论最大无重复置换次数为N!( n为样本量);②符号置换,通过改变样本的正负号,适合单样本 t检验(与0比较),或者其他对称零分布的统计,不要求方差齐,理论最大无重复置换次数为2 N;③混合置换,即顺序+符号置换,需同时满足两种置换条件,理论最大重复置换次数为N!×2 N,适合样本量比较少的统计
(二)基于皮层表面的形态学测量
VBM方法是基于脑容积的定量分析技术,其基本单元是体素。脑灰质主要由皮层和皮层下核团组成。对于脑皮层而言,VBM方法采用部分容积效应模型来刻画灰质与白质,灰质与脑脊液的关系。VBM分割得到的组织边界是渐进性的,对一些交界区的阳性结果发现往往很难确定到底是灰质还是白质起源。另外,脑皮层的沟回高度迂曲,在容积空间紧邻的体素很可能属于细胞构筑完全不同的沟回,再加上基于容积的空间平滑等因素,导致VBM等容积分析方法的定位偏倚。因此,有学者提出了脑皮层表面重建技术。其核心思想是采用硬分割算法提取皮层灰质,即确定灰质与脑脊液边界(软脑膜面)和灰质与白质边界(脑白质面),并通过三维表面重建技术对灰质皮层进行可视化。基于皮层表面的重建技术的基本单元称为顶点,它是皮层表面相邻三角形网格的汇集点。根据重建得到的软脑膜面和白质面,可以提取出皮层厚度和表面积等定量信息。利用GLM等统计学模型,既可以在脑区水平进行皮层厚度和表面积的统计学推断,也可以类似VBM,基于顶点水平进行统计学分析,这些都统称为基于皮层表面的形态学测量(surface-based morphometry,SBM)。基于皮层表面的形态学测量能有效避免VBM分析中部分容积效应导致的脑区定位偏倚问题。该方法一般包括以下几个步骤:皮层重建,空间配准,空间平滑和统计学分析(图1-6-7)。
1.皮层重建
脑皮层表面重建的方法有很多,例如,FreeSurfer软件包采用先分割确定白质面,然后在其基础上使用膨胀算法确定软脑膜面;而CAT12软件包采用的是VBM的组织分割算法先确定灰质、白质和脑脊液成分,然后再确定灰质/白质,灰质/脑脊液边界。值得注意的是,不管什么重建算法,原始采集得到的高分辨T 1WI数据都需要进行B1场偏倚校正,使脑内白质的信号均匀一致,从而确保脑白质分割的准确性。另外,对于膨胀算法而言,还需要预先进行剥脑操作,也就是去除硬脑膜、颅骨及头皮等脑外结构。剥脑的好坏决定了软脑膜面重建的准确性。如果残留了太多的脑外结构,很可能把这些成分错误地标记为灰质,导致软脑膜面错误地延伸到脑外;同样,如果把正常脑组织去掉了,那么重建出的软脑膜面将和白质面重叠,导致脑皮层表面出现镂空。为了更准确地剥脑,可以引入多种模态的MRI数据,例如,高分辨T 2WI。通过把灰(白)质边界或者灰质/脑脊液边界转化成以小三角为基本单位的网格,即构建出最初的皮层表面。由于局部分割错误,原始构建的皮层会出现拓扑结构缺损(镂空或者错连),因此还需要采用一些方法如球谐函数修复缺损。
图1-6-7 SBM分析流程
2.空间配准
与VBM一样,为了使人与人之间的脑皮层指标(如厚度和表面积)在空间上具有对应性,SBM分析中也需要用到线性配准和非线性配准。线性配准一般用于个体脑容积图像与标准脑的大体对齐(仿射变换),为非线性配准做准备;或者用于纵向研究中不同时间点的脑影像配准(刚体变换),使重复测量的数据具有可比性;或者用于不同模态脑影像的配准(刚体变换)。在非线性配准方面,由于皮层表面没有厚度的概念,因此有别于基于容积的配准方法。比较主流的方法是基于球面的配准。首先,个体脑皮层表面需要转化到球面坐标空间。然后用非线性算法(如DARTEL等)基于某些皮层特征(如脑回迂曲度,脑沟深度或者其他沟回标志等)把个体皮层向参考皮层模板对齐。也有研究把非线性容积配准和表面配准方法相结合,能显著提高皮层及皮层下核团配准的准确性。
3.空间平滑
皮层表面的空间平滑通过引入二维高斯平滑核对每类皮层形态学指标(如厚度和表面积)求卷积,平滑后每个顶点的值代表了它周围网格皮层指标的加权平均值。平滑的目的与VBM一样,包括:提高数据分布的正态性,弥补配准不良以及增加统计的敏感性等。
4.统计分析
与VBM类似,基于皮层的形态学测量也常使用GLM统计学模型,既可以在脑区水平进行皮层厚度和表面积等指标的统计学分析,也可以在顶点水平进行统计学分析。有关GLM模型的介绍参见VBM部分。
(三)脑形态学协共变网络分析
复杂网络分析能很好地揭示人脑的功能分化与整合,而图论是目前复杂网络分析最主要的数学工具。根据图论原理,复杂网络由两个基本元素组成,即节点和边。在传统人脑结构网络定义中,节点被定义为神经元或者它们构成的集群,边被定义为神经元之间的突触连接或者纤维束连接。人脑网络可以从微观(神经元)、介观(神经元柱)和宏观尺度(脑区)水平进行研究。基于弥散成像的纤维束追踪是目前活体内在宏观尺度上研究人脑结构网络的最主要方法。
除了纤维束以外,基于脑区之间形态学协共变信息也能够反映脑区之间的连接特性。形态学协共变是指在某群体中两个脑区形态学特征分布(例如,皮层厚度和灰质体积等)的统计学相关性。大量研究表明,形态学协共变与脑区之间解剖连接的强度高度相关,反映了脑区之间的营养供给、相似的遗传基础及相同经验依赖的可塑性等。根据上述研究基础,国内学者贺永教授首次提出了利用局部脑区皮层厚度的协共变信息构建脑网络。随后,脑灰质体积、皮层复杂度、皮层面积等其他形态学指标也被用于脑网络构建中,这些分析统称为脑形态学协共变网络分析。
脑形态学协共变网络分析的基本流程包括:计算形态学指标,脑区分割(定义网络节点),协共变分析(定义网络边),构建脑网络矩阵,计算拓扑属性和统计分析(图1-6-8)。
图1-6-8 脑形态学协共变网络分析流程
1.计算形态学指标
各种脑形态学指标均可以用来构建协共变网络,包括脑灰质体积、雅可比行列式、皮层厚度、皮层面积和皮层复杂度等。有关这些形态学指标的计算参见VBM分析和皮层构建部分。
2.脑区分割
在人脑网络分析中,脑区被定义为网络的节点。在宏观尺度脑网络分析中,理想的节点需满足脑区内体素之间连接特征高度相似,而脑区之间的连接特征高度不同;另外,节点在不同被试间具有较高的解剖对应性。目前最常用的节点定义是根据已发布的脑图谱划分脑区,例如,基于个体脑细胞构筑划分的Brodmann图谱、基于个体脑回特征划分的AAL图谱、基于组水平脑回特征划分的Harvard-Oxford概率图谱、基于纤维连接特征聚类的组水平Brainetome脑图谱以及根据皮层沟回信息划分的个体化皮层图谱等。也有研究组尝试进一步细化脑区分割,例如,把AAL脑区自动切割成体素大小相似的1 000个或者更多的脑区;或者把每个体素当成独立的节点。
3.协共变分析
用于定义复杂网络的边。与在个体水平构建人脑网络的边不同,形态学协共变网络是基于组水平构建网络的边,即计算任意两个脑区皮层厚度等形态学指标在人群中分布的统计相关性,例如,皮尔逊相关。对任意一组而言,只会得到该组两脑区形态学分布的协共变相关系数。每组样本量越大,得到的相关系数就越稳定,越能反映该组人群脑区间形态学指标的内在关联。
4.构建脑网络矩阵
脑网络矩阵是横纵坐标均为节点,任何两个节点之间的交点为边的二维方阵。在协共变网络中,节点为每个脑区,边是脑区之间形态学指标的协共变相关系数。根据边的定义,复杂网络可以分为二值网络和加权网络。对于基于形态学协共变分析的二值网络而言,需要首先确定一个相关系数阈值,阈值以下的矩阵元素赋值为0,阈值以上的矩阵元素赋值为1。对于加权网络而言,方阵中的值为连续变量,即协共变相关系数本身。加权网络既可以通过阈值化只保留阈值以上的那部分边(阈值以下赋值为0),也可以基于全部相关矩阵构建。
5.计算脑网络拓扑属性和统计分析
脑网络矩阵构建后,不同尺度,不同类型脑网络拓扑属性的计算都是相同的,具体参见“脑网络分析”章节。由于脑形态学协共变网络分析构建的是组水平的脑网络,拓扑属性只有一个标量,因此不能采用传统的参数统计对拓扑属性进行组水平的显著性推断。一种普遍使用的办法是置换检验。该方法通过千次以上随机打乱被试标签或者分组标签,然后用置换后的数据构建脑网络,计算拓扑属性,构建统计零分布,最后计算真实值(如拓扑属性的组间差别)在该零分布中的概率,从而做出统计学推断。
二、脑弥散分析方法
(一)弥散成像基本原理
弥散运动也称布朗运动,是分子在微观环境中相互碰撞而形成的无规则运动现象。由于弥散运动的无规则性,现有技术手段很难追踪单个分子的运动轨迹。为此,Albert Einstein采用了位移分布方程从宏观水平定量描述分子的自由弥散运动的快慢(式1-6-1)。
其中 t为检测时间, r为分子在检测时间内平均弥散位移, k是常数, D即弥散系数,它反映了分子的弥散强弱。
人体组织中水含量占全身体重的70%左右。弥散运动反映了水分子在一定时间内的空间位移信息,位移检测需要引入外源性示踪剂,而且生物组织绝大部分不透明,因此传统光学手段很难用于检测活体内水分子的弥散。磁共振现象的发展和应用为活体内检测水分子弥散运动提供了理论基础和实现工具。
DWI是指测量组织中水分子弥散运动特征的磁共振成像技术。其核心思想是把组织内水分子的位移分布转化为磁共振信号衰减。1965年,Stejskal和Tanner两位学者提出梯度脉冲序列的概念。其基本原理是在自旋回波序列重聚脉冲两侧放置一对完全相同的梯度脉冲(强度、方向和作用时间均相同),也称弥散梯度脉冲。第一个梯度脉冲使沿梯度方向上的水分子产生去相位,通过180°重聚脉冲对自旋相位进行翻转,第二个梯度脉冲使弥散位移为0的水分子产生完全的相位重聚,因此,这类水分子MR信号不受梯度脉冲的影响。对于在两梯度脉冲间隔时间内发生了位移的水分子(弥散位移不等于0)而言,由于受到的梯度场作用不再和第一次相同,因此相位重聚效应减弱甚至进一步加重散相,导致MR信号的衰减。弥散梯度权重可以用b值来度量,b值越高,弥散信号衰减就越强。同样,组织水分子弥散越快,其弥散信号衰减也越强(图1-6-9)。
图1-6-9 弥散加权成像
在自由弥散运动成立的前提下,组织内水分子的弥散系数可以用单指数衰减模型来度量(式1-6-2):
其中 S 0为没有施加弥散梯度的原始MRI信号; S为施加了弥散梯度强度为 b时测得的MRI信号,即弥散加权信号; D为弥散系数,反映了测量组织内水分子的弥散能力。因此,理论上只需要测量2次不同弥散梯度b值下的弥散加权信号,即可计算出弥散系数 D。
在生物组织中,水分子的弥散能力并不是在各个方向上都是相同的。对于均质的组织而言(例如,脑灰质),水分子在各个方向上的弥散能力大致相等,这类弥散称为各向同性弥散。对于具有高度方向性的组织而言,例如,脑白质,由于细胞膜和髓鞘的存在,细胞内的水分子在垂直于轴突方向上弥散运动明显受限,导致在该方向测量得到的 D值明显低于沿着轴突方向上的 D值,这类弥散称为各向异性弥散(图1-6-10)。通常用ADC来表示生物组织中的弥散能力。
各向异性弥散为评估白质纤维束的方向性提供了新思路,有学者在20世纪90年代初提出了DTI技术。DTI最初专指利用二阶张量模型度量生物组织水分子的弥散系数方向依赖性的弥散加权成像技术。该数学模型需要满足两个基本条件:①水分子弥散服从单指数衰减规律;②每个体素内只有一个弥散主方向。在张量模型中,各向异性组织中水分子的弥散位移分布可以被模拟成一个椭球。在特定的测量时间内(即弥散时间),每个体素内某弥散敏感梯度方向上的弥散位移分布代表了该方向上从起点到椭球表面某一个点的距离。通过采集不同弥散敏感梯度方向上的DWI,即可获得椭球表面的一系列采样点,从而构建椭球的大小、形状和方向。由于二阶张量是对称的,因此,仅需6个独立的元素即可描述张量特征,意味着只需采集6个非共线弥散梯度方向的DWI即可求解该张量模型(图1-6-10)。在实际应用中,考虑到各向异性弥散的复杂性以及噪声等因素,常需要采集更多方向上的弥散图像,然后用一定的数学算法(例如,线性最小二乘法)拟合张量的6个基本元素。通过线性变换,可以把二阶张量分解成两类基本元素,本征向量和本征值。本征向量表征张量椭球中三个正交轴的方向,代表了体素中水分子弥散的方向信息。而本征值是与弥散方向无关的标量,用来表征三个本征向量上的弥散运动强度,即这些方向上的弥散系数(图1-6-11)。式1-6-3为弥散张量的数学表达式。
图1-6-10 各向同性弥散与各向异性弥散
其中 D代表张量,包含6个基本元素, λ代表本征值, ε代表本征向量。
传统弥散成像理论假设组织内的水分子弥散需服从高斯分布,而且在空间内的弥散分布服从二阶张量分布。遗憾的是,这两个假设不是在所有脑组织中都成立。Stejskal-Tanner(S-T)梯度脉冲序列是弥散成像的基石,其核心思想是利用强梯度场放大水分子弥散运动引起的质子去相位,并检测去相位导致的MR信号衰减。在生物组织中,MR弥散信号还受微观结构、水交换、灌注等影响,导致表观弥散运动不完全服从高斯分布,因此,基于单指数弥散衰减模型和二阶张量的传统弥散成像技术在评估生物组织弥散方向特征上存在一定缺陷。例如,细胞膜等微观结构把组织大致分为细胞内和细胞外两个不同的弥散空间。由于细胞膜脂质双分子层的疏水特性,水分子不能自由通过细胞膜,导致两种成分内的水分子弥散行为是不同的。对细胞内的弥散而言,细胞膜的相对封闭性,当弥散时间(△)足够长的时候,水分子到达细胞膜后就不能自由向远处弥散,而是局限于细胞/轴突内,这种弥散称为受限弥散。显而易见,细胞内弥散不服从高斯分布,也就不能用单指数模型进行准确定量。对细胞外(细胞间隙)的弥散而言,尽管水分子弥散过程也会受到细胞膜等微观结构的阻挡而发生重定向,但是由于细胞间隙是相通的,总会找到向外弥散的出口,因此在宏观上,细胞外弥散服从高斯分布,可以用单指数模型精确定量,只是弥散轨迹相对自由弥散而言更加迂回,测得的ADC值减低,这种弥散称为受阻弥散。细胞内外的水分子尽管不能自由通过,但是可以在消耗ATP的情况下,开放细胞膜脂质双分子层上的水通道蛋白,从而实现细胞内外水分子的交换。因此,细胞膜水通道蛋白的功能状态也会对弥散信号产生显著影响。例如:在急性脑缺血时,由于葡萄糖和氧气不能及时供应,ATP消耗殆尽,导致水通道蛋白不能开放,细胞内外的水交换显著降低,此时水交换对弥散信号的贡献也会显著减低。另外,在S-T序列两脉冲梯度间隔时间内,微循环灌注也会引起显著的水分子随机位移(甚至比弥散导致的位移更大),从而产生灌注相关的信号衰减。总之,通过S-T序列测量得到的弥散信号成分比较复杂,不能简单用单指数模型进行精确定量(图1-6-12)。
图1-6-11 弥散张量椭球对脑内弥散的表征
图1-6-12 磁共振弥散加权信号的成分
另外,受到主磁场和梯度磁场等硬件技术的限制,目前针对人的磁共振弥散成像能够获得的体素大小一般都大于1mm 3(常规3T磁共振设备弥散图像体素尺寸一般在8mm 3以上)。因此,在一个白质体素内存在成千上万条轴突,如果这些轴突构成的白质纤维束不是平行走行的(例如,有交叉或者分叉),那么二阶张量就不能准确地对体素内的多个方向进行表征,导致追踪的错误,特别是对一些复杂脑白质结构,例如,半卵圆中心区的交叉纤维。
因此,自弥散成像提出三十年多来,研究者一直致力于弥散数学模型的改进。相继提出了一系列参数模型(例如,多指数模型、多张量模型、受阻受限混合模型等)和非参数模型(例如,弥散谱成像、Q球成像、弥散峰度成像等)。除此以外,近年来随着磁共振硬件性能的大幅提高,也促进了弥散相关成像技术的发展(如并行成像、多层同时成像、压缩感知等技术),在图像质量提高的同时,显著缩短了成像时间,使很多还停留在理论上的数学模型(例如,弥散谱成像)得以成为现实,促进了弥散成像技术的高速发展和应用推广。
(二)脑弥散成像数学模型及定量指标
1.经典弥散成像及定量指标
传统弥散成像假设组织内的水分子弥散需服从高斯分布,而且在每个体素内水分子弥散只有一个主弥散方向。在满足上述假设情况下,水分子弥散可以用单指数衰减模型和二阶张量表征,由此可以衍生出一系列弥散指标(图1-6-13)。
(1)表观弥散系数:
即通过单指数模型求解得到的弥散系数。在生物组织中,弥散系数不仅包含了水分子自由弥散运动的信息,还包含组织微观结构的信息(例如,细胞膜对水分子弥散的阻挡和限制),因此一般把生物组织测量得到的弥散系数称为ADC(式1-6-4)。
图1-6-13 经典弥散成像定量指标
(2)本征值(eigenvalue, λ):
通过对弥散张量进行线性变换得到的与方向无关的标量,度量张量椭球三个本征方向上的弥散系数。
(3)轴向弥散系数(axonal diffusivity, AD):
即最大本征值,反映了水分子沿轴突方向上的弥散能力,主要与轴突完整性有关。
(4)径向弥散系数(radial diffusivity, RD):
中间本征值和最小本征值的平均,反映了垂直于轴突方向上的弥散能力,主要与髓鞘及细胞间隙大小有关。
(5)平均弥散系数(mean diffusivity, MD):
即三个本征值的平均,反映了组织内水分子的平均弥散能力(式1-6-5)。
(6)分数各向异性(fractional anisotropy,FA):
用来表征组织的各向异性程度,范围为0~1,0代表各向同性弥散,值越大,各向异性程度越显著(式1-6-6)。
2.弥散受阻受限混合模型
弥散受阻受限混合模型(composite hindered and restricted model of diffusion,CHARMED)假设组织中的水分子弥散分为一个受阻成分和若干个受限成分。其中受阻弥散主要反映了轴突外弥散的贡献,而受限成分主要代表了轴突内弥散。和传统DTI使用单个参数( b值)来定义弥散梯度场不同,CHARMED模型采用两个参数来定义梯度场,其中 q值代表了弥散梯度的作用强度( q=γδG/2 π),△代表弥散时间。通过改变 q值和△值采集弥散加权像,并使用式1-6-7数学模型求解出两类不同的弥散分布。
使用CHARMED模型可以获得多种弥散指标,例如,容积分数(volume fraction, f)反映了不同弥散成分对弥散信号的贡献;以及每种成分的弥散系数。也有学者根据该模型提出了AxCaliber技术,能够估计白质纤维束内的轴突直径和密度。另外,CHARMED模型可以评估2个以上受限成分的方向分布,因此,可以区分单个体素内交叉的纤维。与传统DTI指标相比,这些指标能更特异地揭示脑白质微观结构特征,提高了白质定量和纤维束追踪的准确性。
3.弥散峰度成像模型
弥散峰度成像(diffusion kurtosis imaging,DKI)模型在传统DTI模型的基础上额外增加了一个峰度项,用于评价水分子弥散位移偏离高斯分布的程度,其数学表达式见式1-6-8:
其中 是高斯部分, 是峰度部分。峰度张量是四阶对称的矩阵形式,有15个独立元素,因此至少需扫描15个弥散敏感梯度方向。DKI线性公式中含有两个未知系数( D和 K),因此,至少需要2个非零 b值数据才可以求解,并且尽可能采集更多的高 b值数据。DKI常用指标有张量指标和峰度指标(图1-6-14)。张量指标与传统DTI一样,如 MD、 FA、 AD、 RD。常用的峰度指标包括:
(1)轴向弥散峰度(axonal kurtosis,AK):
沿张量椭球体主轴方向上(最大本征向量)的弥散峰度值。因为此方向上弥散受限相对较小,弥散偏离高斯分布的程度较低, K值也比较小。
(2)径向弥散峰度(radial kurtosis,RK):
垂直于弥散主轴方向上弥散峰度的平均值。这些方向上弥散受限最严重,因而RK值较AK高。而且白质纤维径向上弥散受限明显,白质RK值高于灰质。
(3)平均弥散峰度(mean kurtosis,MK):
组织沿各个方向弥散峰度的平均值。MK值越大表明弥散受限越严重,成分结构越复杂。
4.弥散部分容积模型(partial volume model,PVM)
也称球-棒模型(ball-stick model,BSM)。该模型假设弥散信号由一个各向同性的高斯弥散成分和数个各向异性的棒型受限弥散成分组成,后者反映了轴突内的弥散。该模型理论上能够刻画单个体素内任意方向上的弥散特征,从而区分交叉纤维。PVM模型适用于各类采集方案的弥散数据,应用比较广泛。通过PVM模型可以计算每个体素内不同纤维束成分的部分容积分数(partial volume fraction,PVF),代表了体素内各纤维对弥散信号的贡献比例(式1-6-9)。
图1-6-14 DKI峰度指标
其中 S i为弥散信号, S 0为没有施加弥散梯度的MRI信号, f为受限弥散容积分数, b为弥散强度, r为弥散方向, d为弥散系数, N为需要评估的最大纤维数目,RAR T为沿主弥散方向上的各向异性张量。
5.Q空间成像(Q-space imaging,QSI)
通过对采集的弥散信号进行直接傅里叶变换获得弥散传播函数,即概率密度函数(probability density function,PDF)。该方法能够全面地表征组织中水分子的弥散位移分布,而不需要对组织的弥散分布进行特定的假设,理论上来说是最理想的弥散模型,其数学公式见式1-6-10:
其中 q值代表了弥散梯度的作用强度,定义为( q=γδG/2 π); E代表在特定弥散时间(△)内,弥散信号随 q值变化的函数; P s代表弥散传播函数PDF; R代表净位移。 Q空间是基于 q值大小和弥散梯度方向的3D坐标空间。 Q空间成像的核心代表是弥散谱成像(DSI),它需要对3D笛卡尔空间进行密集采样填充 Q空间,然后利用快速傅里叶变换来求解弥散传播函数,这种方案常需要采集500以上方向的弥散数据才能准确求解。为了减少采样的次数,有学者提出了一种混合弥散成像(hybrid diffusion imaging,HYDI)方案,该方案通过采集Q空间的数个同心球面的数据。HYDI除了大大减少了采样数外,还可以运用于各种重建模型,例如,DSI、DKI、Q-ball等(图1-6-15)。
可以根据PDF提取出一系列弥散定量指标:
图1-6-15 Q空间梯度编码填充方案
(1)零位移概率(zero displacement probability, P 0):
描述了在弥散时间内水分子不弥散的概率,越大说明弥散能力越弱(式1-6-11)。
(2)位移均方(mean-squared displacement,MSD):
描述了体素内平均弥散系数(式1-6-12)。
同时,也可以从PDF中提取出方向信息,如弥散方向分布函数(orientation distribution function,ODF)。ODF直观描绘白质纤维在体素内各个方向上的分布情况,因此根据ODF可以表征单个体素内多条纤维束交叉的信息(图1-6-16)。
QSI不需要特定的数学假设,直接对原始信号进行定量和方向解析,因此,理论上来说是最理想的弥散模型。但是,QSI对弥散分布函数的准确评估需要满足一些前提假设:首先,弥散梯度脉冲宽度(δ)要足够小,以保证水分子在梯度作用期间弥散距离尽可能小;其次,弥散时间(△)要足够长,以便水分子能够到达细胞膜的边界。然而,实际情况下,由于梯度场等硬件技术限制,上述前提并不能完全满足,导致评估误差。另外,QSI需要采集大量不同弥散强度以及不同弥散梯度方向的数据以填充 Q空间,因此采集时间非常长,制约了其在科研和临床中的广泛应用。最近,随着多层同步采集技术(simultaneous multislice,SMS)或多频带采集技术(multi-band,MB)和压缩感知技术(compressed sensing,CS)等的提出,打破了QSI推广的瓶颈,表现出巨大的潜力。
6.高角度分辨率弥散成像(high angular resolution diffusion imaging,HARDI)
QSI需要对 Q空间密集采样,成像时间太长。为了解决这一问题,一种折中的采集方案是只需要对Q空间的最外层球面的弥散信号进行采样,即单个高 b值(2 000s/mm 2及以上)、多弥散方向(如60个方向以上),然后采用一定的数学算法直接求解ODF。Tuch等首先提出了HARDI的概念,最初用的是多张量模型(multi-tensor model,MTM)。它假设每个体素的弥散是由数个弥散张量的线性组合,每个张量代表一种纤维成分,其弥散服从高斯分布。该方法获得的ODF能够区分单个体素内的多条交叉纤维。(式1-6-13)
其中 f j代表每个弥散张量 D j的表观容积分数, E( q k)代表标准化弥散信号, q k代表弥散梯度( q k =γδg k), D代表弥散张量, τ为弥散时间。
除了MTM模型外,HARDI也衍生出了一系列其他数学模型,例如,Q球成像(QBI),该模型采用Funk-Radon变换或球谐函数直接求解ODF;球形反卷积(spherical deconvolution,SD)通过引入单条纤维弥散基函数,对每个体素的ODF进行反卷积,获得纤维ODF(fiber orientation diffusion function,fODF)。该方法显著提高了交叉纤维解析的角度分辨率,但也可能导致假阳性纤维。ODF除了能够表征纤维的方向外,也可以衍生出定量指标,例如,广义分数各向异性(generalized fraction anisotropy,GFA)、量化各向异性(quantitative anisotropy,QA)、各向同性系数(isotropic value,ISO)等。
图1-6-16 方向分布函数ODF与基于QSI的纤维束追踪
(三)弥散纤维束追踪技术
纤维束追踪(fiber tracking,FT)基本原理是假定张量或者ODF中弥散主方向与纤维走行一致,采用一定的数学算法,把相邻体素的主本征向量根据相似性连接起来,即构成了虚拟的白质纤维束。纤维束追踪大体可以分为确定性和概率性两种类型。确定性纤维束追踪假定相邻两采样点之间的连接是唯一的,而概率性追踪则用概率方式表征某采样点与周围相邻采样点之间纤维连接的可能性。
最简单也最常用的确定性纤维束追踪算法是流线法,也称连续纤维束追踪算法(fiber assignment by continuous tracking,FACT)。该方法假设三维空间是连续的,先人为确定一个种子点,然后步进很小的距离(通常1/2体素或更小),寻找主本征向量最相似所对应的采样点,并与前一个采样点连接起来,依次步进从而构建出通过该种子点的曲线,即感兴趣的白质纤维束。
由于图像噪声的存在,不管用什么数学模型评估得到的主本征向量都可能存在偏差,从而导致纤维追踪轨迹偏离真实的路径,而且这种偏差会随着追踪距离增加而不断累积,导致纤维束追踪假阳性和假阴性错误。减少纤维追踪假阳性最有效的方式之一是采用纤维编辑技术,即多感兴趣区技术。首先,在追踪一条特定的纤维束时,根据先验知识在真正纤维走行路径上定义多个种子点,只有同时通过这些种子点的纤维才保留(逻辑“AND”);如果通过上述办法还存在假纤维(也需先验知识确定),可以在假纤维走行路径上定义额外的种子点,然后删除通过这些种子点的纤维(逻辑“NOT”)。纤维编辑只能够减少纤维追踪的假阳性,要减少纤维追踪的假阴性,除了需要改进弥散评估模型外(例如,用ODF替代tensor),也有一些替代FACT的追踪算法,例如,根据基于全脑信息的纤维束追踪,或者概率性纤维束追踪技术。
概率流线法(probabilistic streamlines,PS)是最常用的概率性纤维束追踪算法。该方法首先采用Bootstrap等算法评估种子点体素中每个弥散主方向的概率分布,然后随机抽取分布内的某个方向进行流线法追踪,重复上千次从而得到通过该种子点体素纤维的概率分布图。概率性纤维束追踪表征的是通过某个体素或脑区的纤维束概率,可以有效减少图像噪声导致的追踪假阴性;但是,也会引入大量低概率的假阳性纤维轨迹,常需要根据先验知识人为设置一定的阈值来过滤这些低概率的假纤维(图1-6-17)。
(四)脑弥散成像伪影及校正技术
脑弥散成像要求超快的图像采集以减少弥散时间(△)内整体头动导致的去相位,同时也要尽量减少总的采集时间以保证受试者能耐受。单次激发自旋回波平面回波成像(single-shot spin-echo echo-planar imaging,SS-SE-EPI)序列因为满足上述条件成为脑弥散成像的首选。然而,该序列也会引入一些固有的伪影,常需要进行校正。
1.弥散成像常见伪影
基于SS-SE-EPI技术的弥散成像常见伪影主要包括以下三种类型(图1-6-18):
(1)局部磁场不均匀相关伪影:
SS-SE-EPI技术填充k空间是连续的,任何导致局部磁场不均匀的因素(例如,磁化率)均会在每次K空间相位编码线填充过程中产生相位误差,而且会不断累积,最后导致图像严重变形和(或)信号缺失。矩阵越高,伪影越严重,限制了弥散成像空间分辨率的提高。
(2)涡流伪影:
在弥散成像中,高强度弥散梯度场的快速切换会产生电涡流,引起弥散图像发生变形,而且不同弥散方向上图像变形方向也不一样,导致不同弥散梯度方向上的弥散图像对齐不良。
(3)运动伪影:
弥散成像常需要扫描很多不同方向以及不同弥散强度的图像,成像时间比较长(数分钟至数十分钟不等),期间任何头动都会影响弥散图像的质量。主要表现为三种伪影:出现在成对弥散梯度之间(弥散时间△)的剧烈头动或搏动会导致该层图像信号显著降低甚至缺失;出现在同一个TR内不同层面采集期间的头动可导致层面间对齐不良(例如,奇偶层错位);出现在不同TR间的头动导致不同脑容积对齐不良。
2.弥散成像伪影的校正或预防措施
(1)头部固定:
在现有技术条件下,固定头部是降低头动导致的各种伪影的有效方法。
(2)缩短采集时间:
能提高受试者(特别是有轻度意识障碍的患者)的耐受力。通常缩短采集时间会以损失图像信息量为代价。近年来,随着多层同时采集(MB或SMS)技术和压缩感知技术(CS)的提出和推广,可以在不损失图像信息的前提下大幅缩短采集时间。
图1-6-17 纤维束追踪基本原理
A.构建ODF/Tensor;B.提取弥散主方向;C.连接主方向相似的相邻体素;D.全脑纤维
(3)校正局部磁场不均匀所致的伪影:
针对一些体积较小结构的弥散成像(例如,视神经、脊髓等),可以选择矩形视野激发技术,相比传统的单个选层梯度和全层射频激发,这类技术通过施加两个方向上选层梯度,能够只激发层面内特定宽度的组织,因此只需要对激发区域内进行相位编码,大大减低了相位编码步数,显著降低了局部磁场不均匀导致的图像变形,而且不发生卷叠。另外,近年有学者提出了分段EPI弥散成像技术,这些方法把每层弥散数据的k空间沿着相位编码方向(multi-shot EPI序列)或者频率编码方向(RESOLVE序列)分成若干段在不同TR区间采集,显著降低了相位误差的累积,从而降低了图像变形。分段采集技术的成熟为高分辨弥散成像奠定了基础,但是其缺点是采集时间成倍增加,信噪比也有待提高。如果设备只有SS-SE-EPI弥散成像序列,可以评估扫描范围内磁场不均匀性分布,即场图(field mapping,FM);然后利用该分布对EPI图像进行校正,该方法能有效纠正EPI图像的变形,但是不能纠正磁场不均匀导致的信号丢失。
(4)电涡流和头动校正:
其基本原理是利用仿射变换,把每个弥散方向或强度采集的脑容积图像向自身参考图像(常选用第一帧脑容积图像)进行配准。该方法除了能纠正涡流所致的弥散图像变形外,同时也校正了扫描过程中不同TR间的头动。但是该方法不能纠正层面间的头动引起的错位。
图1-6-18 SS-SE-EPI弥散图像常见伪影
(五)弥散图像预处理
获得弥散加权图像后,在数学建模前,需要对数据进行预处理,以提高纤维束追踪以及弥散指标解算的可靠性。
1.图像变形校正
由于组织磁化率不同,导致单次激发EPI图像上组织交界区(例如,眶额、颞下回以及脑干等部位)的弥散图像严重变形,影响配准和对应脑区的定量统计。可以事先评估扫描范围内磁场不均匀性分布,即场图,然后利用场图对EPI图像进行校正,该方法能有效纠正EPI图像的变形。场图可以通过采集不同TE值的梯度回波序列估计,也可以通过反转弥散序列的相位编码方向估计。图像变形校正可以用FSL软件包中的FUGUE和topup工具包实现。
2.涡流和头动校正
其基本原理是把每副脑弥散图像向参考图像(常选用第一个脑容积)进行仿射线性配准。该方法除了能部分纠正涡流所致的弥散图像变形外,同时也校正了扫描过程中不同TR间的头动。但是该方法不能纠正层面间的头部运动所致的错位。
3.剥脑
去掉头皮等非脑组织结构不但可以显著减少模型拟合和追踪过程的运算量,还能提高配准的准确性。剥脑的基本原理是寻找脑组织和颅骨的边界,根据方法的不同大致可以分为膨胀法和组织分割法。膨胀法见于FSL中的BET工具包,该方法首先确定头颅的大体中心和半径,构建初始化的球面,然后向外膨胀寻找脑组织与颅骨的边界。组织分割法常用于T 1WI的剥脑,即通过组织分割获得不同组织成分,然后只把脑灰质、白质和脑脊液合并起来构成脑的蒙片,然后提取蒙片内的脑结构。
(六)脑弥散定量分析技术
除了脑梗死、脑肿瘤等少数病变可以在弥散指标图或纤维束图上肉眼观察到异常外,大部分个体变异或者疾病所致的弥散特征变化都是微弱的,需要借助计算机和统计学等工具来提取与个体变异或者疾病相关的弥散特征,由此衍生出一系列弥散定量分析技术。
1.基于感兴趣区(region of interest,ROI)的弥散定量分析
该分析方法需要根据先验知识,按照一定标准勾绘ROI,然后提取ROI内的弥散定量值进行后续的统计学分析。ROI分析是假设驱动的分析方法,其优点包括:操作简单;不需要被试间空间配准,因此不受配准不齐影响;假设明确,结果容易解释。其缺点包括:ROI的定义受预设方案和测试者影响,故有必要进行可重复性验证;只能对部分具有先验知识的脑区定量,不能反映其他脑区的特征及变化(图1-6-19)。
2.基于直方图的弥散定量分析
直方图分析是纯数据驱动的分析技术。该方法把既定脑组织范围(如全脑)自动分为若干个单元,计算每个单元平均弥散定量值,然后绘制弥散定量的直方图分布:横坐标为升序排列的定量值,纵坐标为每段取值范围对应的单元数或者比例。直方图分布特征(例如,曲线下面积、峰度系数、偏度系数等)可以表征个体之间全脑水平弥散特征的差异。该方法的优点包括:自动提取特征,避免了测试者主观和先验知识不准确的偏差;不需要被试间空间配准,因此不受配准不齐影响。缺点是:只能反映整体水平的弥散特征,不能对脑区进行定位(图1-6-20)。
图1-6-19 基于感兴趣区的弥散定量分析示例
A.人脑内囊ROI;B.猫脑内囊ROI
图1-6-20 基于直方图的弥散定量分析示例
3.基于纤维束的弥散定量分析
基于纤维束的弥散定量分析是在个体纤维束追踪的基础上,提取纤维束经过脑白质区的弥散定量信息,如纤维数目、FA值、MD值等。既可以对整条纤维束弥散特征进行定量,也可以获得纤维束路径中的每一段的定量信息。它具有ROI分析的基本特征,例如,假设驱动、不需要被试间配准、需要先验知识、局部脑组织定量等。但是相比ROI分析,纤维束提取的脑组织特征更明确(能具体到某条纤维束),边界勾绘也更加准确,因此,逐渐取代基于ROI的脑白质弥散定量分析。值得一提的是,纤维束定量的准确性取决于纤维束追踪的准确性(图1-6-21)。
4.基于纤维束图谱的弥散定量分析
上述基于纤维束的弥散定量分析尽管准确性比较高,但需要对每个个体每条纤维进行独立的纤维束追踪,工作量非常巨大,不适宜大样本分析。为了解决大数据分析问题,有学者提出了基于图谱的弥散定量分析。该方法首先获得标准空间组水平纤维束图谱,然后把每个个体的弥散指标图配准到标准空间,把纤维束图谱作为ROI投射到个体指标图上,提取每条纤维束的弥散指标并进行统计分析。纤维束图谱既可以从现有的公开数据库中获得,例如,苏黎世大学基于尸体解剖构建的组织学纤维图谱,以及约翰霍普金斯大学采用弥散纤维束追踪生成的活体脑纤维概率图谱(图1-6-22)。也可以用户自己构建纤维束图谱,其基本流程如下:
图1-6-21 基于纤维束的弥散定量分析
A.胼胝体连合纤维;B.皮质脊髓束
图1-6-22 基于纤维束图谱的弥散定量分析流程图
(1)个体水平纤维束追踪:
从总体中选择一部分被试,追踪出感兴趣的纤维。
(2)纤维束配准:
一般用b 0像向标准空间配准。根据配准方法的不同,分为线性配准和非线性配准。然后利用变形场把个体纤维束转换到标准空间。
(3)生成组水平概率模板:
由于个体变异以及配准的偏差,同一条纤维束在个体之间的走行不会完全重叠,因此,有必要对所有个体的同一条纤维束进行平均,生成组水平的纤维束概率图谱,其中每个体素的值代表了该体素中追踪出的纤维在人群中的概率。
(4)阈值化:
图谱中概率低的体素,很可能反映了配准不齐或个体变异太大,因此,不具有代表性。有必要设置一定的阈值过滤掉低概率纤维的体素分布。
基于纤维束图谱的弥散定量分析优点是效率高,适合大样本纤维束定量分析;缺点是纤维束图谱只纳入纤维束走行一致性高的体素,忽略了纤维的个体变异。另外,个体图像向标准空间配准过程的误差(配准不良)会降低弥散定量的准确性。
5.基于体素的弥散定量分析
基于体素的弥散指标定量分析和VBM类似,是对每个体素的弥散指标进行独立统计分析,获得统计分布图,从而做出统计推断(图1-6-23)。这类分析统称为基于体素的分析(voxel-based analysis,VBA),其基本流程如下:
(1)数据预处理及弥散指标提取:
参见“弥散图像预处理”和“脑弥散成像数学模型及定量指标”部分。
(2)空间配准:
为了个体间的脑弥散指标之间具有可比性,VBA分析中也需要用到空间配准。和VBM分析所用的高分辨率T 1WI相比,弥散图像的空间分辨率低,常有变形,而且灰白质对比度也比较低,给配准带来了巨大困难。最初的VBA分析是用FA图直接向标准空间的FA模板非线性配准,然而,由于FA图的灰质信号低,白质信号虽高且异质性也很高,配准后的脑结构特别是脑灰质常异位比较严重,导致定位的不准确。原始弥散图像扭曲变形是导致配准不良的主要原因之一,因此,在预处理过程中纠正扭曲变形,能显著提高配准的准确性。针对原始图像空间分辨率低的问题,可以采用两步配准法提高配准的准确性,即首先把弥散b 0像和同一被试的高分辨率T 1WI进行线性配准;然后采用DARTEL等非线性配准技术,把个体T 1WI配准到标准空间的T 1像模板;然后合并上述两步配准获得的变形场,把弥散指标写入到标准空间。配准图像与标准空间模板的对比度越相似,配准越准确,因此生成针对特定研究的自定义模板也能显著提高配准的精度。还有学者提出基于张量的配准法,与传统基于单变量信息(例如,FA值)配准法不同的是,该方法是一种多变量配准技术,即把弥散方向信息引入模型评估中,目的是使不同被试每个体素的弥散张量的方向分布保持一致,该方法也被证实能显著提高配准精度。
(3)空间平滑:
配准后弥散指标图可以用三维高斯平滑核求卷积,平滑后每个体素的值代表了它周围体素弥散指标值的加权平均。平滑的目的是提高正态性、弥补配准不良以及增加统计的敏感性等。
图1-6-23 VBA分析流程图
(4)统计分析:
与VBM类似,VBA也一般使用GLM模型进行组水平统计学分析。有关GLM模型的介绍参见VBM部分。
VBA分析是一种纯数据驱动的方法,特别适合于探测未知脑区的个体变异和疾病损伤,但是由于空间配准的问题比较严重,限制了其推广。随着弥散图像空间分辨率的提高、图像伪影的减少以及空间配准方法的改进,其运用前景会越来越广泛。
6.基于白质骨架空间统计分析(tractbased spatial statistics,TBSS)
TBSS的提出是为了解决早期VBA分析的缺陷。VBA分析面临的第一个问题是配准不良,为此,TBSS通过构建组水平的白质FA骨架,然后把每个被试垂直于骨架方向上的最高FA值投射到骨架上,部分纠正了配准不齐;其次,VBA和VBM一样,引入了空间平滑以提高参数检验的效力,然而同样引入了未知的偏差,TBSS不对数据进行平滑,而改用非参数检验,由于非参数检验不要求样本服从正态分布,因此选择非参数检验有效避免了平滑引入的偏差(图1-6-24)。在FSL软件包中,TBSS有一套流程化的处理脚本,其基本过程如下:
(1)数据预处理及弥散指标提取:
参见“弥散图像预处理”和“脑弥散成像数学模型及定量指标”部分。
(2)空间配准:
TBSS采用的默认配准方法是一步非线性配准法,即把个体FA图直接和FA图模板进行对齐。FA图模板既可以是软件自带的FA模板(FMRIB58_FA_1mm.nii),也可以是用户自定义的FA模板,还可以是被试间相互配准,找到最相似的个体FA图作为模板。FA图直接配准会导致脑灰质显著变形,也有学者对空间配准进行了改进,例如,引入两步配准法,或者基于张量的多变量配准法。
(3)构建组水平白质FA骨架:
首先对所有被试配准后的FA图进行平均,生成组水平的FA图;然后利用重力场中心理论(center of gravity,COG),在组平均FA图上寻找垂直于白质束方向上最强FA值对应的体素,并把这些体素连接起来,即构成仅有一个体素厚度的白质FA骨架。需要强调的是白质骨架与具体的白质纤维束没有直接关联,而是最强FA构成曲面或者曲线。初步生成的白质骨架包含FA比较低的体素,这些体素常不在白质区,因此需要对白质FA骨架进行阈值化(例如,FA > 0.2),仅保留位于白质区的骨架体素。最后,对白质骨架进行二值化,用于随后的个体骨架指标图生成。
(4)生成个体骨架FA图:
组水平白质骨架相当于一个容器,需要把每个个体的FA图投射到组水平骨架上,才能进行后续分析。配准的偏差,组水平骨架与个体FA图直接对应的体素FA值未必最高。因此,TBSS采用了一种优化的个体FA值填充方案:对每个骨架上的体素,寻找个体FA图上垂直于骨架方向上的最强值,作为个体FA图在骨架上的投影。通过这种投影,每个被试垂直于骨架方向的配准不齐得到了纠正。个体被试其他弥散指标,例如,MD、AD等也可以根据FA的投射轨迹填充到白质骨架上。
(5)非参数统计:
投射到白质骨架上的弥散指标不需要空间平滑,直接进行统计分析。由于骨架上的弥散指标分布未知,因此TBSS采用了一种非参数统计方案:置换检验。该方法的基本原理是随机打乱样本顺序或者符号,对每次置换后的数据进行相同的GLM建模和参数估计。置换后得到的效应量(如两组均值的差异)理论上为零,通过计算机进行成千上万的迭代运算,便可构建出真实的零分布。最后比较真实效应相对该零分布的概率,从而获得统计概率 P。有关统计的详细介绍请参阅VBM分析部分。
TBSS既纠正了VBA中的配准问题,又避免了空间平滑引入的未知误差,因此,在基于体素的脑白质统计分析中得到了广泛应用。然而,TBSS也存在固有的局限:首先,TBSS只关注了脑白质区FA值最高的部分体素,它不能揭示其他白质结构特别是皮层下脑白质的变异,因此,与VBA比较,其提供的信息明显要少得多;另外,尽管TBSS能对垂直于骨架方向上的配准不良进行纠正,但是,对于其他方向上的配准不良没有任何矫正措施;甚至有研究表明,对于相邻的脑白质束(例如,扣带束和胼胝体),TBSS很可能错误地把白质束A的体素投射到白质束B上。TBSS是为了解决VBA的配准问题而提出的,随着VBA配准问题的逐步改进,TBSS的价值可能会逐渐减低。
7.脑解剖网络分析
人脑是由大量神经元和它们之间的突触连接构成的高度复杂的系统,功能分化与整合是其两大组织原则。功能分化是指大脑中不同的脑区具有相对特异的功能;功能整合是指人脑完成某一特定的功能往往需要多个脑区相互合作。复杂网络分析能很好地揭示人脑的功能分化与整合,而图论是目前复杂网络分析最主要的数学工具。根据图论原理,复杂网络由两个基本元素组成,即节点和边。在人脑解剖网络中,节点就是神经元或者它们构成的集群,边就是神经元之间的突触连接或者它们构成的纤维束。人脑解剖网络可以从微观(神经元)、介观(神经元柱)和宏观尺度(脑区)水平研究。弥散成像能在活体内追踪脑区之间相连的纤维束,目前已经成为在宏观尺度研究脑解剖网络最重要的神经科学工具。
图1-6-24 TBSS分析流程图
图1-6-25 脑结构网络分析流程图
基于弥散成像的脑解剖网络分析主要过程包括:弥散数学建模、脑区分割(确定网络节点)、纤维束追踪(构建网络边)、构建网络矩阵、计算拓扑属性和统计分析(图1-6-25)。
(1)弥散数学建模:
选用什么类型的弥散模型是决定纤维束追踪准确性的基础。传统DTI理论的二阶张量模型假设每个体素的水分子弥散只有一个主本征向量,因此,只能表征体素内一个纤维束方向。由于成像技术的限制,目前它仍旧是最常用纤维束追踪模型,但是使用者必须清楚它很可能错误地评估纤维复杂脑区的弥散分布,导致纤维追踪的假阳性和假阴性。如果条件允许,构建人脑网络尽量采用优化的能够评估复杂纤维的弥散模型,例如,球-棒模型、HARDI、QSI等。值得注意的是,运用不同采集方案的弥散数据适用于不同的数学建模:例如传统DTI适用于低弥散强度( b≤1 000s/mm 2),梯度方向数较少(<60)的弥散数据;HARDI适用于高弥散强度( b>1 500s/mm 2),梯度方向数较多(≥60)的弥散数据;QSI适用于多种弥散强度组合以及较多梯度方向数的弥散数据等。近年来随着磁共振硬件技术的发展(如高切换率梯度线圈和多通道接收线圈),多层同时采集技术及压缩感知技术在弥散成像中的推广,HARDI/QSI等复杂模型在脑解剖网络分析中的应用越来越广泛。
(2)脑区分割:
用于定义人脑解剖网络的节点。理想的节点应满足脑区内体素之间解剖连接特征高度相似;脑区之间的解剖连接特征明显不同;节点在不同被试间的解剖对应性好。目前最常用的节点是基于个体人脑的公共模板(brodmann、AAL等),它们是基于一个被试获得的,不能反映人群的脑区变异;而且主要是基于细胞构筑或者沟回信息划分,同一个脑区内体素的连接相似性可能存在显著差异;再加上它们对脑区的划分都比较大(brodmann约48个,AAL 116个),因此不是理想的人脑网络节点。为此,中国科学院自动化所脑网络组中心蒋田仔教授团队在前人脑区划分基础上,根据结构连接特征把脑区进行了更精细的划分,而且基于大样本数据,构建出组水平节点的最大概率图谱(连接特异性脑区共246个),为大尺度人脑网络分析提供了更加理想的节点定义方法。也有研究组尝试进一步细化脑区分割,例如,把AAL脑区自动切割成体素大小相似的1 000个或者更多的脑区;或者把每个体素当成独立的节点。但是,这些节点定义方法的合理性需要验证。
(3)纤维束追踪:
对任意两个脑节点之间纤维束进行追踪,构建脑解剖网络的边。根据追踪算法可以分为确定性追踪和概率性追踪。有关追踪的具体细节参见“纤维束追踪”部分。这一步需要进行许多次纤维束追踪( N×[ N-1]/2),常采用计算机程序自动化进行(例如,PANDA工具)。
(4)构建脑解剖网络矩阵:
以每个脑区作为网络的节点,以脑区之间的纤维束定量作为网络的边构建脑解剖网络矩阵。脑解剖网络矩阵是横纵坐标均为节点,任何两个节点之间的交点为边的二维方阵。根据边的定义,复杂网络可以分为二值网络和加权网络。对于基于弥散成像的二值网络而言,如果两个脑区之间存在纤维连接,方阵中的它们的交点就赋值为1,反之为0。对于加权网络而言,方阵中的值为纤维定量,采用的纤维定量指标可以是纤维数目、纤维密度(单位体积或面积的纤维数目)或纤维的弥散指标(如平均FA、MD等)。
(5)计算脑网络拓扑属性和统计分析:
脑解剖网络矩阵构建后,不同尺度、不同类型脑网络拓扑属性的计算和分析方法都是相同的,具体参见“脑网络分析”章节。
三、动脉自旋标记脑灌注分析
(一)动脉自旋标记(arterial spin labeling,ASL)成像基本原理
脑血流量(cerebral blood flow,CBF)是度量脑组织中毛细血管床动脉血传输效率的指标,标准单位为每分钟每100g脑组织通过的动脉血毫升数。正常脑组织的平均CBF约为60ml/(100g·min),如果假设脑组织平均密度为1g/ml,脑组织平均CBF可以转化为0.01/s,代表每秒约有1%全脑组织容量的血液被新鲜动脉血替换。
CBF的测量一般是通过在动脉中团注外源性或者内源性对比剂,然后动态检测对比剂输送到脑组织以及从中清除的过程,从而计算出每个体素脑组织的CBF。脑灌注可以通过PET、CT或者MRI测量。其中PET是活体检测CBF的金标准,然而由于成本高、电离辐射等缺点,限制了其广泛使用。通过注射外源性对比剂的CT检查和MRI检查尽管成本比较低,然而也存在一定的创伤和副作用(例如,过敏风险),而且CT检查也有电离辐射。为此,20世纪90年代末有学者提出了把标记的动脉自旋作为内源性对比剂进行灌注测量的磁共振技术。
ASL采用饱和脉冲标记一段动脉血,由于饱和动脉血(弛豫不完全)产生的MRI信号明显低,当它流经特定脑组织时,也会导致该脑组织内体素的平均信号减低。通过分别采集含有或者不含有饱和动脉血(自旋标记)的MRI图像,计算它们的信号差,即可计算出脑组织的CBF(图1-6-26)。根据标记脉冲的不同,ASL可以分为以下几种采集方法:
1.连续标记ASL(continuous arterial spin labeling,CASL)
采用连续射频脉冲(2~4s)饱和感兴趣区上游一薄层组织,流过该层的动脉血将被连续标记。该方法的优点是标记充分,灌注加权图信号高。但缺点是磁化传递(magnetic transfer,MT)效应显著,而且射频能量沉积(SAR值)高。
2.脉冲式标记ASL(pulsed arterial spin labeling,pASL)
该方法使用时间非常短的脉冲标记一厚块组织内的动脉血。pASL又可以分为对称和非对称标记。对称性标记的典型代表是FAIR(flow alternating inversion recovery)技术,该技术在采集参考图时使用非选择脉冲,而在采集标记图时使用选择性脉冲。非对称性标记的典型代表是EPISTAR(echo-planar MR imaging and signal targeting with alternative radiofrequency),该方法标记感兴趣脑区动脉上游(例如,颈部)的10~15mm厚的组织。该技术衍生出一系列改进的方法,如TILT(减少MT效应)、PICORE(提高标记效率)、QUIPPS-Ⅱ和Q2TIPS(降低动脉通过时间的影响)等。pASL的优点是SAR值和MT效应均明显减小,缺点是灌注信号低,常需要数十次重复采集以提高信噪比。
图1-6-26 动脉自旋标记成像原理
3.伪连续标记ASL(pseudo-continuous arterial spin labeling,pCASL)
也称为脉冲式连续ASL。该技术通过引入持续一段时间的单脉冲串标记一薄层组织,兼顾CASL和pASL的优点,即相对pASL灌注信号明显提高的同时,相比CASL其MT效应和SAR值明显减低。
ASL技术不需要注射外源性对比剂,是无创无辐射的安全检查手段,特别适合于重复性或者纵向研究。另外,通过分析ASL的时间序列,还可以揭示脑活动过程中的灌注变化,为脑功能研究提供了新的手段。因此,基于ASL的脑灌注成像和分析在科学研究和临床运用中越来越普遍。
(二)ASL脑灌注定量
基于减影获得的ASL灌注加权信号理论上与CBF值呈正比,其定量表达式见式1-6-14:
其中△ M代表减影信号即灌注信号, M a ,0代表动脉血完全弛豫的MR信号(质子信号), c( τ)代表动脉输入函数(arterial input function,AIF), r( t-τ)代表标记自旋清除率,而 m( t-τ)代表了纵向弛豫效应。上述公式在CASL和pASL都适用。
需要注意的是,利用上述公式对CBF准确定量需要满足以下前提条件:标记动脉血到达不同脑组织的时间相等以及自旋在组织中是自由弥散的。虽然灌注信号与CBF呈正比,但是还可能受到动脉到达时间、团注时间宽度、标记质子血-组织分布系数等影响,导致CBF的评估不准确。例如,常规ASL都采用的是单个标记后延迟时间(post-labeling time,PLT),由于不同组织(特别是在脑卒中以及动脉狭窄性疾病中)动脉到达时间(arterial transit time,ATT)不尽相同,采用相同PLT引起不同组织的标记效率不一样,导致CBF的评估误差。一种行而有效的解决方法是分别采集不同PLT下的标记图像,即多反转时间ASL(multiple inversion times arterial spin echo,mTIASL)。该方法能获得动脉自旋进入以及流出脑组织的灌注信号变化,从而计算出组织特异的ATT,能显著降低单反转时间标记导致的CBF计算误差,而且还可以衍生出脑血容量(brain blood volume,CBV)等定量指标。
(三)基于体素的脑灌注定量分析
基于体素的脑灌注定量分析(voxel-based perfusion analysis,VBP)和VBM类似,是对每个体素的脑灌注指标(如CBF)进行独立统计分析,获得统计分布图,从而做出统计推断。该分析方法适用于解决未知的脑科学问题,或者揭示脑疾病的神经代谢异常。其基本过程如下(图1-6-27):
1.头动校正
图1-6-27 基于体素的脑灌注定量分析流程图
该步骤只针对采用单次激发平面回波序列(SS-EPI)的ASL成像。ASL通常要重复采集上百次SS-EPI图像,成像时间内受试者头动会导致不同TR采集的脑容积发生空间错位,因此有必要进行头动校正,把不同脑容积进行空间对齐。和经典fMRI分析类似,ASL采用刚性配准算法评估并校正脑容积间的头动。
2.灌注指标计算
CBF是最主要的ASL灌注指标。如果采用了mTI-ASL成像技术,还可以提取出ATT、CBV、MTT等灌注指标图。需要注意的是,单PLT计算得到的CBF值受ATT影响,为了校正不同被试之间脑平均ATT的差异,可以对CBF进行归一化,比较常用的办法是用每个体素的CBF值除以全脑所有体素的平均CBF值获得相对CBF。
3.空间配准
为了使得个体间的脑灌注指标在体素水平具有可比性,需要把每个个体的CBF图配准到标准空间。可以直接把原始图像用非线性算法配准到标准空间模板上,也可以把灌注加权像直接配准到标准空间模板。需要注意的是,一步非线性配准法需要保证待配准图像和标准模板的对比度一致,例如,EPI图像应该向EPI模板配准,而灌注加权像应该向标准空间的灌注加权像模板配准。灌注原始图像空间分辨率通常比较低,可以尝试用两步配准法来提高配准精度,即首先把原始像和同一被试的高分辨率T 1WI进行线性配准;然后采用DARTEL等非线性配准技术,把个体T 1像配准到标准空间的T 1模板;然后合并上述两步配准获得的变形场,把CBF等灌注定量图转换到标准空间。待配准图像与标准空间模板的对比度越相似,配准越准确,因此,生成针对特定研究的自定义模板也能显著提高配准的精度。
4.空间平滑
配准后灌注指标图需用三维高斯平滑核求卷积,平滑后每个体素的值代表了邻近体素灌注指标值的加权平均。平滑的目的与VBM一样,包括提高数据分布的正态性、弥补配准不良以及增加统计的敏感性等。
5.统计分析
和VBM类似,VBP也采用GLM模型进行组水平统计学分析。有关GLM模型的介绍参见VBM部分。
四、脑激活分析方法
(一)任务态fMRI简介
fMRI由Ogawa等人于1992年提出,是一种检测大脑神经活动的方法。目前最常用的是基于血氧水平依赖(blood-oxygen-level-dependent,BOLD)的技术,也称BOLD-fMRI。fMRI技术具有非侵入性和无电离辐射等优点,自出现以来很快在人脑功能研究中得到了大量的应用。
BOLD-fMRI技术的基本成像原理为:当脑内某局部区域的神经活动增强时,该区域的神经元需要更多包括氧在内的能量,而神经元对能量的需求通过神经血管偶联机制及局部代谢的改变引起局部血管扩张,从而导致脑血流量(CBF)和血容量(CBV)的增加。由于神经元的耗氧量小于周围血管的供氧量,从而导致该局部区域的含氧血红蛋白的浓度增加以及脱氧血红蛋白浓度的下降。脱氧血红蛋白是一种顺磁性物质,其浓度的变化会引起局部磁场强度的变化,进而引起该区域磁共振信号的变化。简单来说,局部区域神经活动强弱的变化会导致周围区域血氧浓度的变化,进而引起磁共振信号的变化。由此可见,fMRI技术测量的并不是神经活动本身,而是神经活动所引起的血流动力学变化,是对神经活动的间接测量。与神经活动的变化相比,血流动力学的变化要缓慢得多,因此,fMRI信号具有明显的延迟,在神经活动瞬态变化之后的6~8秒达到峰值,20~30秒才能完全恢复基线状态。fMRI信号的这一延迟特性可以通过血流动力学响应函数(hemodynamic response function,HRF)来刻画,如图1-6-28所示。在对fMRI数据进行分析时,必须充分考虑fMRI信号的这一特点。
图1-6-28 血流动力学响应函数(HRF)
根据实验时被试所处状态的不同,fMRI又可分为任务态(task-based fMRI)和静息态(restingstate fMRI,rs-fMRI)。简单地说,任务态fMRI的目的是为了考察人脑在执行某一特定任务时的神经活动情况,因此,要求被试在实验过程中按照要求执行某一事先设计好的任务;而静息态fMRI则不要求被试在实验过程中完成任何特定任务,只需安静地躺在MRI扫描仪中接受扫描即可。本部分主要介绍任务态fMRI数据的激活分析方法。
由上述BOLD-fMRI成像基本原理可知,当出现某种特定任务刺激时,与处理该刺激信息有关的脑区所观测到的fMRI信号会增强。基于这一现象,研究者提出检测与实验任务有关的脑激活区的思路:观察脑内哪些体素的fMRI信号变化与实验任务中刺激的出现保持同步,信号变化与刺激序列保持同步的脑区即为任务激活脑区。换句话说,将刺激出现时的fMRI信号与刺激未出现时的fMRI信号进行比较,若出现刺激时某体素的fMRI信号显著高于没有刺激时的fMRI信号(即基线状态),则说明任务刺激引起了该体素fMRI信号的变化,即激活了该体素。这一激活称为正激活,反之,若刺激条件下某体素的fMRI信号显著低于基线状态,则为负激活。
任务态fMRI实验的基本流程主要包括四个部分:实验任务设计、数据采集、数据分析、结果解释与呈现。如图1-6-29所示。下面将对这四个部分分别进行介绍。
(二)实验任务设计
根据任务刺激的呈现方式,任务态fMRI实验任务设计大致可分为三种:组块设计、事件相关设计和混合设计。
1.组块设计
图1-6-29 任务态fMRI实验基本流程
组块设计是最简单也最易操作的一种设计。在组块设计中,同一类型的刺激成组块式连续呈现,如图1-6-30所示。通常情况下每一组块的持续时间相同,具体持续时间需要根据具体实验任务确定。如果整个设计中仅有一种刺激类型,任务刺激组块需与基线(baseline)组块交替呈现(图1-6-30A)。若有两种或两种以上的刺激类型,根据实验目的的不同,可考虑包括基线组块或者不包括基线组块。若不包括基线组块,不同类型的刺激组块可交替呈现,中间没有任何间隔(图1-6-30B)。若包括基线组块,则在不同类型的刺激组块之间插入基线组块(图1-6-30C)。前面已经提到,神经活动诱发的血流动力学响应具有延迟效应(图1-6-28),因此,刺激组块内当多个刺激连续呈现时,单个刺激所诱发的血流动力学响应(即fMRI信号)会互相叠加,直至达到饱和状态后维持高水平,当刺激不再继续呈现时,fMRI信号逐渐回落至基线水平(图1-6-31)。因此,通过对比某个体素所采集到的刺激组块和基线组块分别对应的fMRI信号水平,就可判断出该体素是否与该任务刺激有关。由此可见,当组块越长(即一个组块内刺激重复次数越多或持续时间越长),fMRI信号越强,信噪比越高,但信号也会达到饱和,应同时考虑每个组块在整个实验中的重复次数,以保证足够的统计力度。根据组块设计的特点,其优势主要有两点:①组块设计的信噪比较高,所得结果稳定性强;②组块设计要求同类型刺激必须连续成块出现,其可变因素较少,设计比较简单。但同时组块设计也存在一定的局限性:①组块设计对刺激的呈现方式有较严格的限制,缺乏灵活性,无法实现一些较为复杂的实验范式;②由于多个刺激连续呈现导致每个刺激所诱发的血流动力学响应相互叠加混杂,无法提取出单个刺激所诱发的血流动力学响应,因此无法对不同刺激所诱发的fMRI信号变化的动态性进行分析;③由于组块设计中刺激的呈现非常有规律,因此被试可以较为精确地预测将要呈现的刺激的种类和时间,使被试对任务刺激有较强的预期性,从而导致脑激活结果中混杂被试对刺激的预期效应;④组块设计中同类型刺激重复出现,缺乏新鲜感,被试易疲劳,导致注意力水平降低。
图1-6-30 组块式设计
图中不同颜色的箭头代表不同类型的实验刺激。A.只有一种任务条件,与基线状态交替进行;B.两种不同的任务条件交替进行;C.两种不同的任务条件与基线状态交替进行
图1-6-31 连续刺激引起的血流动力学响应叠加效应
A.BOLD信号叠加效应,刺激越多,信号越强;B.BOLD信号叠加直至饱和,之后维持高水平,当刺激不再继续呈现时,信号逐渐回落至基线水平
2.事件相关设计
事件相关设计与组块设计最明显的差别在于其可实现对单个事件(即单个刺激)所诱发的fMRI信号变化进行分析。从设计方式上讲,事件相关设计中任务刺激不再有规律地连续呈现,相邻刺激之间必须要有一定的间隔,不同类型的刺激以随机方式呈现。根据相邻刺激呈现的间隔时间长短,事件相关设计又可分为慢速事件相关设计和快速事件相关设计。对于慢速事件相关设计,相邻刺激的间隔时间必须大于血流动力学响应的持续时间,以避免相邻刺激所诱发的血流动力学响应相互混叠(图1-6-32B)。为使任务刺激能够重复足够多的次数以保证统计力度,通常慢速事件相关设计的实验整体持续时间较长。而快速事件相关设计则可以缩短整个实验的持续时间。在快速事件相关设计中,相邻刺激呈现的时间间隔可以小于血流动力学响应的持续时间,但为保证能够提取出单个刺激所诱发的fMRI信号动态变化曲线,相邻刺激的时间间隔必须要有足够的随机化(jitter,图1-6-32C):间隔时间的随机化范围越大,越容易提取出单个刺激的fMRI信号动态变化曲线。例如,相邻刺激的时间间隔在4~16s之间的范围内进行随机要优于在8~12s之间的范围内随机。由于上述刺激呈现方式的灵活性,事件相关设计具有如下优势:①能够检测到单个刺激的血流动力学响应的瞬态变化,进而可以对fMRI信号进行时域上的分析,从而不仅可以比较不同类型刺激脑激活区的差异,而且可以比较不同类型的刺激所诱发的血流动力学响应在时间变化上的差异(例如,峰值时间的差异等);②能够进行trial水平上的分析,根据被试完成任务的情况(如,被试对每次任务刺激完成的正确性、反应时间等)在实验完成后对刺激重新进行分类,得到的脑激活结果则更能反映被试的任务完成情况;③与组块设计相比,能够最大限度地避免被试的预期效应,同时更有利于维持被试的注意力水平;④与组块设计相比,事件相关设计可以降低头动对激活检测的影响;⑤其设计方式的灵活性,可以实现更多复杂的实验范式。具有上述优势的同时,事件相关设计也有一定的局限性:①通常情况下,事件相关设计比组块设计实验耗时更长(采用快速事件相关设计可以在一定程度上减轻这一问题);②刺激诱发的fMRI信号的信噪比相对较低,导致统计结果也相对较弱;③实验设计的灵活性会使实验设计的复杂性提高,可变因素多,不易控制,对很多随机因素(如不同类型的刺激相邻呈现时,是否会产生交互作用等)所带来的影响也很难预测。
图1-6-32 不同实验设计类型示意图及对应的预期BOLD信号响应
A.组块设计;B.慢速事件相关设计;C.快速事件相关设计;D.混合设计
3.混合设计
混合设计是将组块设计和事件相关设计的特点结合在一起形成的一种设计类型。简单来讲,混合设计在实验整体上看是组块设计,但在组块内部主要采用事件相关设计。通常情况下,对于混合设计,同一组块内仍然仅有一种刺激类型,而在组块内部,刺激不再有规律地连续呈现,而是对刺激呈现的间隔进行了随机化(图1-6-32D)。组块设计和事件相关设计各有优缺点,而混合设计则将它们各自的优缺点进行了折中。
实际应用中除了需要确定刺激呈现方式(组块设计、事件相关设计、混合设计)之外,还需要根据具体情况确定实验范式的整体架构,包括每个被试需要进行几次扫描以及每次扫描过程中实验任务应如何安排和呈现等。例如,若要研究某药物对大脑功能的影响,应采用重复测量的设计,以避免被试个体间差异对结果造成的影响。即,对于每个被试,在用药前和用药后均进行两次扫描,每次扫描称为一个session。两个session唯一的差别在于被试是否用药,而两个session的实验任务设计应保持完全一致,以避免其他非药物因素对结果带来的干扰。在每个session内又可将实验任务重复几次,每次称为一段(run),每段采用完全相同的设计(组块或事件相关设计)。任务刺激每重复一次称为一个试次(trial)。如图1-6-33所示,为研究某药物对情绪处理功能的影响,可以设计如下的实验范式整体架构。其中,ISI(inter-stimulus interval)为刺激间间隔,ITI(inter-trial interval)为试次之间的间隔。
图1-6-33 实验范式整体架构举例
上述介绍了实验的整体架构模式及三种实验设计类型,每种类型又有多种参数(如:刺激重复次数,刺激间隔时间,随机化范围等)需要在实际应用中进行选择。如何选择设计类型及其参数才能实现最优的实验设计从而能够最有效地检测预期的实验效应?最优实验设计的最基本原则是对感兴趣的实验效应要有最强的统计力度,同时还要考虑实现的可行性。这两方面往往互相矛盾,例如,较强的统计力度需要同一刺激重复很多次,然而,重复很多次会导致实验时间过长而超过患者的耐受度,可行性差。因此,设计具体实验时需要在这两个方面进行折中。需要强调的是,必须依据所要研究的问题来选取最优的实验设计。同样的实验设计对不同的研究问题其有效度是不同的。而在实际应用中,往往希望通过一个实验去回答多个问题,因此,所选定的实验设计对感兴趣的多个不同问题的解答有效度可能并不相同。一般来讲,需要根据所研究的问题是否关注神经活动所引起的fMRI信号变化的时间动态性来决定采用组块设计还是事件相关设计以及采样频率(即扫描重复时间;高时间分辨率可以从信号中提取出更丰富的时域信息);根据所需要检测的脑激活区的空间范围及空间特异性等来决定其图像扫描的空间分辨率(高空间分辨率可以得到高空间特异性的脑激活)和扫描范围(全脑扫描还是仅扫描特定的感兴趣区);根据任务刺激的具体特性、参与实验的被试的具体状况等因素来决定刺激的重复次数和实验的整体持续时间。所有这些因素应整体权衡,来确定最终的实验设计。
(三)数据预处理
fMRI信号的变化除了受到神经活动所引起的血流动力学响应的影响之外,还有很多其他因素可能导致其信号的变化,即噪声。为了更加可靠地检测任务刺激所引起的脑激活,需要在进行激活分析之前对噪声进行去除以提高数据的信噪比。针对任务态数据的激活分析,通常需要进行时间校正(slice timing)、头动校正(motion correction or realign)、空间标准化(normalization)、空间平滑(smoothing)和时域高通滤波(high-pass filtering)等预处理过程。不同的fMRI数据分析软件(如:SPM、FSL、AFNI等)对数据的预处理操作过程略有不同,下面以SPM为例介绍基本的预处理步骤(图1-6-34)。
1.时间校正
fMRI的图像获取是逐层进行的,在一个TR内完成全脑所有层面的扫描,因此,每一层对应不同的获取时间。这就导致在同一个TR内获取的全脑图像其不同层面在获取时间上并不完全一致。而这种不同层面在获取时间上的错位会给进一步的激活分析带来偏差,因此,需要将所有层面的信号进行校正,使其对应同一个时间点。最常用的校正方法是首先选定某一层面(例如,每个TR的中间时刻所对应的层面)作为参考层(即:将参考层的获取时间点作为参考时间点),再对其他层面的信号通过邻近TR获取的信号值进行插值,从而估算出该层面在参考时间点上的信号值。
2.头动校正
fMRI实验均需持续几分钟甚至更长的时间,最终获得很多幅图像。而被试的头部在实验过程中难免会出现位置上的移动,这就导致从同一被试获取的不同时间点的fMRI图像在空间上并不完全对应。换句话说,同一被试不同时间点获取的图像,其同一坐标位置可能对应不同的体素。这对提取任一给定体素的fMRI信号时间曲线造成困难。为解决这一问题,所有时间点的fMRI图像需进行头部的对齐。一般采用刚体配准方法进行两步对齐:首先,所有时间点的图像与第一个时间点的图像对齐,再计算所有图像的平均图像;然后,再将第一步对齐的所有时间点的图像与第一步产生的平均图像对齐,最终实现所有时间点的图像在空间上对齐。针对每个时间点的图像,刚体配准方法会产生6个头动参数:沿 x、 y、 z三个坐标轴的平移以及围绕 x、 y、 z三个坐标轴的旋转。进而产生6条头动曲线,可以刻画实验扫描过程中被试的头部在成像空间中位置的变化。
3.空间标准化
每个个体的脑形态各不相同,而绝大多数研究需要对多个被试的数据在体素水平上进行组水平统计,这就要求将个体的脑图像进行空间标准化,使不同个体脑图像中的体素建立空间对应关系。其基本思路为:首先构建一个标准脑模板,然后将每个个体脑与标准脑模板进行配准,即对个体脑进行变形使其形态及空间位置变得与标准脑一样。目前最常用的标准脑模板为MNI坐标空间模板,在各大常用fMRI数据分析软件(如:SPM、FSL、AFNI等)中均有提供。早期研究所采用的空间标准化过程是将经过头动校正后的所有时间点的图像的平均图像配准至标准EPI模板,再将其配准参数应用于每个时间点的图像,从而实现所有图像的空间标准化。由于个体脑与标准脑模板在形态上有差异,因此配准过程需使用仿射变换和非线性变换。这种空间标准化过程仅涉及fMRI(即EPI)图像,而EPI图像所反映的脑结构形态信息远不如T 1WI图像丰富,因此,又出现了利用T 1WI作为中间桥梁进行fMRI图像的空间标准化的方法,以提高配准精度。其具体过程为:首先将个体的T 1WI与fMRI平均图像对齐,再将该个体T 1WI进行空间标准化(即:将个体T 1WI配准到标准T 1模板),最后再将所得到的标准化变形参数应用于该个体所有时间点的fMRI图像,从而实现个体fMRI图像的空间标准化。第一步中将同一个体的T 1WI与fMRI图像进行配准的过程称为coregister,是指同一个体不同模态的图像之间的配准。由于需要配准的图像来自同一个体,因此,只需使用刚体配准即可。而第二步中个体T 1WI的空间标准化则需使用仿射及非线性变换。随着配准技术的发展,目前最常用的空间标准化过程更充分地利用了不同脑组织结构(灰质、白质、脑脊液)的形态信息,在组织分割过程中同时考虑已分割好的标准脑灰质、白质及脑脊液模板,将T 1WI的组织分割过程与空间标准化过程融为一体,以实现更为精确的图像配准。
4.空间平滑
对fMRI数据进行空间平滑的主要目的是为了提高图像的信噪比,并降低配准误差及个体间功能差异对后期的组水平统计分析所带来的不利影响。平滑后的图像中每一个体素的取值均为平滑前该体素及其相邻体素信号值的加权平均。通常采用高斯核平滑,即参与平均的体素的权重服从高斯分布。高斯平滑的另一个目的是使图像中的体素更加服从高斯分布,以满足后期对统计结果进行随机场理论校正的假设条件。高斯平滑核的大小可以用半高全宽(full-widthhalf-maximum,FWHM)来表示,其决定了参与平均的体素的多少。显然,高斯平滑核越大,参与平均的体素越多,其平滑后的图像空间分辨率越低。实际应用中,平滑核的大小需要根据具体问题来选择。若与实验任务相关的脑功能区范围较大,则较大的平滑核能够更加有效地增强信噪比,提高统计力度。反之,则应选择较小的平滑核,以防止平滑后空间分辨率过低而无法检测出面积较小的任务激活区。
5.时域滤波
很多因素可能会导致信号基线随时间而产生缓慢变化,对后期的分析产生不利影响。例如,随着扫描仪的长时间运行,其温度的变化可能会引起信号基线的漂移。因此,在进行激活分析之前,通常需要对fMRI数据进行高通滤波(high-pass filter)以去除基线漂移等低频噪声。常用软件对滤波阈值均设有默认值,例如,SPM默认采用1/128Hz的滤波阈值,FSL默认采用1/100Hz的滤波阈值。在实际应用中,必须根据实验设计来确定软件设定的默认滤波阈值是否合适:研究者感兴趣的与任务相关的信号变化频率必须高于滤波阈值。若阈值设置不当,可能会将数据所包含的与任务相关的信号去掉而无法检测出应有的激活脑区。
图1-6-34 任务态fMRI数据预处理的一般过程
在上述5个预处理步骤中,时间校正与时域滤波属于时域上的处理,而其他预处理步骤则属于空间域上的处理。下图总结了一般情况下任务态fMRI数据的预处理过程。但要注意这些预处理过程并非一成不变,其顺序也可根据情况灵活调整。
(四)个体水平激活分析
前面已经提到,激活体素的检测是基于该体素的fMRI信号随时间的变化是否与实验任务刺激的呈现保持同步。这种同步性可以通过某一体素的fMRI时间曲线与任务刺激序列所对应的血流动力学响应曲线之间的相关性来度量。为更方便地考察某一体素的fMRI时间序列与多种不同类型的任务刺激之间的关系,通常采用GLM来对体素的fMRI时间序列进行建模。GLM将任一体素的fMRI信号变化看成是多种因素所引起的fMRI信号变化的线性组合,其数学形式为:
其中, Y为某选定体素的fMRI时间序列, X i为某种可能引起fMRI信号变化的因素(如某种类型的任务刺激或头动等噪声)所对应的时间序列, β i为回归系数,ε为残差。这里的回归系数 β i即衡量了 Y与 X i之间的关系。需要注意的是, β i与相关系数并不完全相同, β i值的大小还反映了 Y与 X i之间的倍数关系,即斜率。若 β i显著不等于零,则说明 Y与 X i之间存在显著关系。根据 β i符号的不同,又可分为正相关关系(即正激活)或负相关关系(即负激活)。除了考察某个感兴趣任务条件所对应的 β i(即该任务条件所对应的脑激活)之外,还可以考察不同任务条件所对应的多个β的任意线性组合的统计显著性。例如,要考察两种不同任务条件(如 X 1, X 2)所引起的脑激活的差异,可以表达为 β 1- β 2(常被称为contrast),并对其进行统计显著性检验。该contrast可以写成行向量和列向量的矩阵乘积,如式1-6-15:
根据要考察的实验效应的不同,可以定义不同的contrast。根据对contrast所采用的统计检验模型的不同,又可分为T-contrast和F-contrast。以上述公式所代表的一般线性模型为例, X 1, X 2, X 3分别代表三个实验任务条件,表1-6-1总结了常见的T-contrast和F-contrast。
表1-6-1 统计检验建模
续表
由此可见,任务态fMRI数据激活分析的核心问题即建立恰当的一般线性模型,对模型中的 β i进行估计,并检验其统计显著性。
针对任意给定的contrast,对全脑所有体素重复进行上述分析,即可得到该contrast所对应的实验效应的统计参数图(statistical parametric map,SPM),进而从该统计参数图中判断与该实验效应有关的显著激活的体素。若将全脑看成一个整体,这种针对所有体素重复进行统计分析的过程会大大提高整体犯错误概率,因此,必须进行多重比较校正以将假阳性率控制在合理范围内。关于多重比较校正及其常用方法详见本节VBM统计分析部分。
一般线性模型中的回归子 X i大致可分为两类:感兴趣的实验相关因素和不感兴趣的干扰因素。在建立一般线性模型时,应注意两点:①不同的 X i之间不应存在线性相关关系,否则, β i的最优估计存在无穷多个解,无法对 β i的估计值进行有效解释;②由于 β i的统计显著性与残差的大小成反比关系(即:在 β i的估计值相同的情况下,残差越小,显著性越强),因此,所建立的一般线性模型应尽可能包含所有可能引起fMRI信号变化的因素,以提高一般线性模型的拟合度,降低残差,从而增强对 β i的统计检验敏感性。
1 st level GLM分析的基本原理及过程如图1-6-35所示。
(五)组水平激活分析
由于噪声的影响以及脑功能的个体间差异,通常情况下,我们并不需要获得每个个体被试的脑激活区,而是需要对所有采集的被试进行组水平的统计分析,获得组水平的统计参数图,进而得到群体水平公共的激活脑区。组水平统计分析是基于每个个体得到的contrast图进行的。基于实验目的不同,需要采用不同的统计检验。例如,若只有一组被试,可采用单样本 t检验获得该组的脑激活结果;若要考察患者组和健康对照组脑激活的差异,则可采用双样本 t检验。更多关于组水平统计检验的类型以及适用的情况,详见VBM统计分析部分。
此外,组水平统计方法还分为固定效应分析(fixed-effect analysis)和随机效应分析(randomeffect analysis)。表1-6-2总结了两者各自的优缺点。从该表可以看出,在样本数允许的情况下,应尽量使用随机效应分析;但当样本数过少时,样本无法代表总体,此时可采用固定效应分析,但所得结论只可用于所采集的这些特定的样本,不可泛化至整个总体。
图1-6-35 1st level GLM分析的基本原理及过程
表1-6-2 固定效应和随机效应分析的比较
五、脑功能连接分析
(一)脑功能连接简介
脑功能连接的概念最早出现在动物的电生理研究中。20世纪90年代初,英国Friston教授等将其扩展到功能成像领域,并将其分为功能连接和效应连接两种类型。功能连接是指空间上远距离的神经生理事件之间在时间上的相关性,用于度量空间上分离的脑区的神经活动之间所存在的统计依赖关系。而效应连接是指一个神经活动单元对另一神经活动单元施加的直接或间接的影响,用于研究一个脑区的神经活动如何对另一个脑区的神经活动进行作用,因此,更强调两个脑区之间功能依赖关系的方向性。下面对功能连接和效应连接常用的分析方法进行简单介绍。
(二)功能连接分析
用于功能连接分析常用的方法有相关分析、偏相关分析、独立成分分析等。下面逐一进行介绍。
1.相关分析
相关分析是最简单也是使用最广泛的应用于fMRI数据功能连接分析的方法。相关分析是用来度量两个变量之间相互依赖关系的一种统计方法。最常见的相关分析是使用皮尔逊相关系数(pearson correlation coefficient)来度量两个变量之间的线性依赖关系。假定两个脑区 X和 Y,其对应的时间序列分别为 x( i)和 y( i), i=1, 2,…, n( n为时间点的个数)(式1-6-16)。
其中, r( x, y)代表 x和 y的皮尔逊相关系数, cov( x, y)代表 x和 y的协方差, σ 2( x)和 σ 2( y)分别代表 x和 y的方差。
根据具体科学问题,我们可以事先选定若干个感兴趣脑区,然后对这些脑区的平均fMRI时间序列两两之间进行相关分析,相关系数值即为功能连接强度,并进一步对所得到的相关系数进行统计检验,判断其显著性,以此确定这些脑区两两之间是否存在显著的功能连接。需要注意的是,由于相关系数本身不服从正态分布(-1≤ r≤1),因此在对相关系数进行组水平统计检验时,若使用单样本或双样本 t检验,需在统计检验前将相关系数进行Fisher’s r-to-z变换,使其更接近正态分布。(式1-6-17)
基于相关系数的功能连接分析的另一种常见形式是基于种子点的全脑功能连接分析,即:考察某一个特定感兴趣区与全脑所有其他体素之间的功能连接。其分析步骤通常为:
(1)根据实验目的选取特定脑区作为种子点,提取该区域内各体素的时间序列,并将所有体素的时间序列进行平均得到该种子区的平均时间序列。
(2)在全脑范围内逐个体素计算其时间序列与种子区平均序列之间的相关系数,并将相关系数值赋予该体素,代表该体素与种子点脑区的功能连接强度,进而得到种子点脑区与全脑所有体素的功能连接图。
(3)对所有体素的功能连接进行统计显著性检验,通过设定显著性阈值来确定与种子区显著相关(有显著功能连接)的脑区。全脑体素数目众多,设定的统计显著性阈值需进行多重比较校正。多重比较校正方法详见本节VBM统计分析部分。
种子点的选取方法多种多样,常见的选取方式有:①通过特定任务状态下的激活脑区来选定种子区;②通过前人研究结果确定感兴趣区中心点坐标,并以某半径画小球,将小球内包含的所有体素选定为种子区;③基于解剖标记手工勾画种子区范围;④也可以利用标准脑图谱事先定义好的脑区来选定种子区。这种基于种子点的功能连接分析方法的局限性在于所获得的脑功能连接信息对种子区域的选取具有一定依赖性,基于不同选取方法得到的种子区可能会得到不同的功能连接结果,间接增加了结果解释的难度和可信度。
2.偏相关分析(partial correlation analysis)
在功能连接分析中,多数情况下分析的不仅仅是两个脑区之间的功能连接,而是多个脑区之间的功能连接关系。而实际上,两个脑区之间的关系可能是通过第三个脑区建立起来的。例如,脑区A和脑区B之间并不存在直接关系,但两者均与脑区C有直接关系,那么,通过相关分析方法就会得到脑区A和脑区B之间也存在显著相关性。而前面提到的相关分析无法区分直接连接和间接连接。为此提出偏相关分析方法。偏相关分析在度量两个脑区之间的功能连接关系时,能够排除其他变量(即其他脑区)的影响,从而将直接关系和间接关系区分开来,因此,能够更加真实地反映脑区之间的功能连接关系。
偏相关分析方法通过控制其他变量,来计算两个脑区之间的偏相关系数。因此,偏相关系数能够度量在排除了其他变量(即控制变量)的影响后,两个脑区之间是否仍然存在显著相关性。偏相关系数的阶数定义为控制变量的个数。例如,当仅有1个控制变量时,为1阶偏相关;当有2个控制变量时,为2阶偏相关。显而易见,0阶偏相关系数与上面提到的相关系数等价。偏相关的计算方法有多种,其中一种较易理解的方法是通过线性回归的方法从感兴趣变量中排除掉控制变量的影响后得到相应的残差,然后再利用残差计算得到的相关系数即为偏相关系数。具体地,假设 X 1和 X 2为感兴趣变量, X 3, X 4,…, X p为控制变量,则通过式1-6-18可以得到残差 ε 1和 ε 2。在控制了变量 X 3, X 4,…, X p的条件下, X 1和 X 2的偏相关系数见式1-6-19。
当仅有一个控制变量时,偏相关计算公式可以简化为式1-6-20:
3.独立成分分析(independent component analysis,ICA)
ICA是一种完全数据驱动的盲源分析方法,在包括fMRI数据在内的众多生物医学信号的分析中被广泛使用。通过fMRI技术观测到的血流动力学响应信号可以看成是由很多相互独立的源所产生的信号的叠加,而ICA的目的就是将这些混合信号分解为若干个相互独立的源信号。fMRI信号可以看成是在时间维度和空间维度上变化的。根据成分分解维度的不同,ICA可分为时间独立成分分析(temporal ICA)和空间独立成分分析(spatial ICA)。需要根据数据维数的特点,来选择不同形式的ICA。通常情况下fMRI信号的空间维数(即体素的个数)远大于其时间维数(即时间采样点的个数),因此,一般采用空间ICA将fMRI数据分解为若干个空间上相互独立的成分。空间ICA在数学上可以表达为式1-6-21:
其中, X代表大小为T×V的数据矩阵,T为fMRI时间点的个数,V为全脑体素的个数,因此,每一行代表一个时间点采集的全脑图像,每一列代表一个体素的时间序列;S代表大小为K×V的独立成分矩阵,K为独立成分的个数,每一行代表一个独立成分所对应的空间图,行与行之间相互独立,每一列代表一个体素在所有独立成分上的取值;M代表大小为T×K的混合矩阵,每一列代表一个独立成分所对应的时间序列,每一行代表一个独立成分在所有时间点上的取值。当进行组水平ICA时,可以将所有人的数据矩阵在时间维度上串联起来,进行降维处理后,再进行ICA分析,这样可以针对每个空间独立成分得到每个被试对应的成分空间图,再进一步进行组水平统计分析(图1-6-36)。
对fMRI数据进行ICA分析得到的空间独立成分,属于同一成分的体素的信号随时间的变化具有相似性,因此,认为这些体素在功能上也具有相似性,而属于不同成分的体素在功能上则相对独立。利用ICA进行功能连接分析正是基于这一思想,认为成分内的脑区(或体素)具有较强的功能连接。
由于ICA是完全数据驱动的,尤其适合静息态数据的分析,因此,ICA在静息态数据的功能连接分析中得到了广泛的应用。也有研究将其应用于任务态数据,发现这一方法也能够有效提取出与任务相关的神经活动区域。此外,该方法是一种多变量分析方法,它在分析中同时考虑了所有体素的信号变化,与单变量分析方法相比,能够更加充分地利用体素间的空间关系信息。由于这些优势,ICA可以从数据中提取出传统方法不易发现的信息,并且能够有效去除与神经活动相互独立的其他噪声的影响,如:呼吸、心跳等生理性噪声以及头动等因素所引起的fMRI信号的变化。
图1-6-36 ICA分析原理示意图
图1-6-37 动态功能连接分析过程示意图
4.动态功能连接分析(dynamic functional connectivity,DFC)
上述功能连接分析方法均假定在实验过程中或在某感兴趣的实验条件下,脑区间的功能连接是静态的(即保持不变的),反映的是在特定的一段时间内,脑区之间的平均功能连接状态。然而,实际上脑区之间的功能连接并不是固定不变的,而是随时间动态变化的。因此,动态功能连接分析方法应运而生,尤其是在静息态fMRI数据分析中得到了广泛应用。研究表明,功能连接的动态变化也反映了脑网络的基本特征。目前已有多种方法可以刻画功能连接的动态变化,其中最为常用的是滑动窗方法。滑动窗方法的思路非常简单,即将整个实验过程所持续的时间分成很多个小的时间段,利用相关分析考察每个小的时间段内脑区之间的功能连接,再将所有时间段内的脑区功能连接强度按照时间顺序排列,即可得到脑区间功能连接强度随时间的动态变化曲线。其具体分析步骤为:①首先设定一个特定长度M的时间窗;②从第一个时间点开始,用该时间窗从各脑区的时间序列(长度为N)中截取出长度为M的一段序列,再利用截取出的这一段时间序列利用相关分析计算脑区间的功能连接;③将该时间窗向右滑动一个步长K(通常,1≤K≤M);④重复第二步,从各脑区的时间序列中截取出新的一段(根据步长K的不同,新截取出的一段时间序列可能与上一步截取的一段时间序列有重叠,重叠的时间点个数为M-K),再次计算脑区间的功能连接;⑤重复上述过程,直至将时间窗滑动至最后一个时间点;⑥将每一次滑动得到的功能连接值按照时间顺序排列,可以得到功能连接的动态变化曲线(图1-6-37)。上述分析步骤的示意图如下:
由上述分析步骤可以看出,所计算得到的功能连接动态性可能与所选择的时间窗长度、滑动步长等参数有关。虽然动态功能连接分析在近几年已经得到越来越多的应用,但对其结果的解释仍然面临诸多困难。例如,所观察到的功能连接动态性的根源究竟是神经活动还是其他噪声,目前仍然很难确定。即便是多个随机噪声,也可能观测到它们之间的相关性随时间变化。而且,呼吸、心跳等生理性噪声、头动或其他能够引起fMRI信号变化的因素都有可能在某些特定长度时间窗下得到时间序列相关性的动态变化。因此,对于动态功能连接分析,需特别注意数据预处理过程中噪声的去除、时间滤波等步骤对分析结果的影响。这一问题在疾病相关的研究中会一定程度地减轻,因为某些噪声(如仪器采样噪声)对两组人群通常是相似的,它们所带来的时间序列相关性的动态变化也是相似的,因此,两组间的统计比较可以在一定程度上降低这一类噪声带来的影响。
(三)效应连接分析
效应连接(effective connectivity)是指一个神经活动单元直接或间接地对另一神经活动单元的影响。它既描述了脑区之间功能连接关系的强弱,同时也能描述脑区之间信息流的传递方向。fMRI数据分析中常用的效应连接分析方法有:结构方程模型(structural equation modeling,SEM),心理-生理交互(psychophysiological interaction,PPI)模型,格兰杰因果分析(Granger causality analysis,GCA),动态因果模型(dynamic causal modeling,DCM)等。
1.结构方程模型
SEM是一般线性模型的扩展,通过建立多个脑区时间序列的协方差矩阵来提取脑区之间的神经交互信息。结构方程模型需要根据先验知识预先定义一个模型,即定义脑区间的因果关系。例如,假定三个脑区的时间序列分别为X、Y、Z,假定它们之间的效应连接关系如图1-6-38所示。
图1-6-38 结构方程模型网络连接示意图
其对应的数学模型表示如式1-6-22所示:
写成矩阵形式即为式1-6-23:
其中, β ij表示脑区 j到脑区 i的效应连接值, φ i为脑区 i所对应结构方程模型的残差。利用极大似然估计法对上述模型中的未知参数 β ij进行估计,即可得到事先假定的脑区之间的效应连接强度值,并可利用残差计算模型拟合度,在个体水平度量参数估计的可靠性和统计显著性。对于组水平统计分析,可以针对每条连接对所有个体的 β ij值进行组水平统计(例如,单样本或双样本 t检验)。
SEM既可以用于任务态数据也可用于静息态数据的分析。但这种方法是建立在协方差基础上的,它忽略了时间序列在时间上先后排列的顺序信息,即:如果将所有脑区的时间序列不同时间点的值以同样方式调换顺序后,其分析结果保持不变。此外,当脑区个数过多或连接过于复杂时,对参数的有效估计会出现困难,该方法则不适用。
2.心理-生理交互模型
PPI模型由Friston等人于1997年提出,用来表示在不同任务状态下脑区之间的连接调制关系,即脑区之间的连接关系是否会随着某一特定实验条件的变化而发生改变。理论上,PPI方法可以用于分析任意两个脑区之间的效应连接,但最常见的是基于种子区的方法,即考察从事先选定的种子区到全脑其他所有体素的效应连接关系。假定全脑任一体素的时间序列为 Y,种子区的时间序列为 X Phy(即生理项),感兴趣的实验条件为 X Psy(即心理项,通常为二值向量,每个时间点取值为1或-1),其他不感兴趣的干扰因素记为G(如头动参数等),则PPI的模型表达式如式1-6-24:
其中, X Psy× X Phy即为心理项和生理项的内积,即心理生理交互项。若模型估计出的 β 1显著不为零,则说明种子区到该体素的效应连接( )会受到该实验条件( X Psy)的调节。
一般情况下,PPI的分析步骤如下:
(1)对任务态fMRI数据执行标准的一般线性模型(GLM)分析,提取某个任务激活脑区作为种子区(根据具体情况,也可利用解剖信息定义种子区)。
(2)生成生理项、心理项及心理生理交互项所对应的时间序列,即: X Phy、 X Psy、 X Psy× X Phy。
(3)把心理生理交互项、心理项、生理项及不感兴趣的干扰因素等各项所对应的时间序列代入上述模型即可估计各项的回归系数( β 1, β 2, β 3, β 4)。我们所关心的是交互项的系数 β 1。
(4)针对每个个体被试,对全脑所有体素重复上述分析,即可得到种子区到全脑每个体素的PPI连接值,从而得到每个被试关于该种子区的全脑PPI连接图(即 β 1图)。
(5)再对所有被试的PPI连接图进行组水平统计分析,进而得到组水平的PPI连接统计参数图,再进行多重比较校正,确定统计阈值,最终得到与种子区有显著PPI连接的脑区,即种子区到这些脑区的连接会受到感兴趣实验条件的调节;根据需要,也可以进行组间比较(例如,患者组与对照组),得到受感兴趣实验条件调节有显著性组间差异的连接脑区。
上述心理生理交互模型还可进一步扩展为生理-生理交互模型,即:将上述模型中的心理项替换为第二个种子区的时间序列,并产生新的生理-生理交互项(两个种子区的时间序列的内积)。生理-生理交互模型反应的是其中某个种子区到全脑任一体素的连接受到另一个种子区神经活动的调节。依据同样道理,还可进一步扩展为更加复杂的交互模型,如两个种子区和一个实验条件的交互模型,即“生理-生理-心理”三因素交互模型,但这种交互作用解释起来比较复杂,目前的应用还较少。
3.格兰杰因果分析模型
GCA是一种多元自回归模型,通过时间先后性来推断两个脑区神经活动的因果关系,进而确定连接的方向。在格兰杰因果关系中,若脑区X过去时间点的神经活动对脑区Y当前时间点的神经活动有预测作用,则认为脑区X的神经活动为因,脑区Y的神经活动为果,即效应连接方向为 。格兰杰因果模型即可用于多脑区间的效应连接分析也可用于基于种子区的全脑效应连接分析。其具体实现形式又分为基于残差的GCA和基于系数的GCA,下面对这两种形式的分析步骤分别进行介绍。
基于系数的GCA分析步骤:
(1)根据所要研究的科学问题选取多个脑区,针对每个脑区,将脑区内所有体素的时间序列进行平均,进而得到每个脑区的平均时间序列。
(2)利用式1-6-25计算多个脑区两两之间的GCA效应连接值:
其中, Y j和 Y i分别代表 n个感兴趣脑区中的第 j个脑区和第 i个脑区的时间序列, t代表当前时间点, p为延迟时间(即多元自回归模型的阶数), t-k代表以当前时间点 t为参照的过去第 k个时间点, Z为不感兴趣的协变量(如头动参数等), ε j代表第 j个脑区所对应模型的拟合残差, 代表延迟为 k个时间点时第 i个脑区到第 j个脑区的GCA效应连接值。若 显著不为零,则说明过去第 k个时间点第 i个脑区的fMRI信号对第 j个脑区当前时间点的fMRI信号有预测作用。通常情况下,fMRI数据的采样时间为2~3秒,一阶自回归模型( P=1)最为常见,即只考察延迟为1个时间点时脑区之间的GCA效应连接。随着成像技术的发展,采用亚秒级采样率的研究也开始出现,因此,更高阶的GCA模型也将随之得到更多的应用。
当只考虑两个感兴趣脑区时,上述模型可以很方便地用于基于种子区的全脑效应连接分析(即:从种子区到全脑所有体素的效应连接分析,或从全脑任一体素到种子区的效应连接分析)(式1-6-26):
其中, Y 1为全脑任一体素的时间序列, Y 2为种子区的平均时间序列, 即为延迟为 k个时间点时从种子区到全脑任一体素的GCA效应连接。同理,若令 Y 1为种子区的时间序列, Y 2为全脑任一体素的平均时间序列,即可得到从全脑任一体素到种子区的GCA效应连接。
基于残差的GCA更多用于基于种子区的效应连接分析。假定 X为事先选定的种子区平均时间序列, Y为全脑任一体素的时间序列,则从种子区到全脑任一体素的GCA效应连接计算公式为式1-6-27:
F x → y即为从种子区 X到任一体素 Y的GCA效应连接。同理,若将 X和 Y调换位置,即: Y为事先选定的种子区平均时间序列, X为全脑任一体素的时间序列,则上述公式即为从全脑任一体素到种子区的GCA效应连接计算方法。
由于GCA分析对效应连接的方向性推断是基于时间先后性,因此,它比SEM所建立的效应连接方向具有更强的物理基础,能够更加直接地反映不同脑区神经活动之间的时间因果链。但同时也应注意其固有的局限性。首先,fMRI信号是对神经活动的间接测量,实际反映的是血流动力学响应。而不同脑区的血流动力学响应可能是不同的,因此,不同脑区fMRI信号的时间先后性未必能够真实反映这些脑区所对应的神经活动的时间先后性。为避免这一局限,有人提出改进的GCA,其核心思想在于:利用事先假定的血流动力学响应函数(HRF)对每个脑区的fMRI时间序列进行反卷积,从而估计出每个脑区的神经活动信号,再将这些估计的神经活动信号代入上述GCA计算公式中,以期获得建立在神经活动基础上(而非血流动力学信号基础上)的GCA效应连接。显然这一方法所得结果的有效性取决于所假定的HRF函数是否符合真实情况。此外,一些预处理步骤(如时间滤波)也会对信号时间先后性的推断带来显著的影响,因此,在应用中应特别注意这些预处理步骤的使用及其参数的设置。
4.动态因果模型
DCM由Friston等人于2003年提出,是建立在确定性状态方程模型基础上的一种效应连接分析方法。DCM将多个脑区的神经活动描述为一个具有因果关系的动态系统:它将大脑看成一个具有“输入-状态-输出”的网络系统,通过外界输入的刺激对该系统进行扰动,进而通过每个脑区的自连接以及脑区间的相互连接来产生网络内每个脑区的神经活动。不同脑区神经活动之间的因果性由模型所假定的效应连接方向来定义。
DCM模型由神经活动状态方程和血流动力学状态方程组成。神经活动状态方程使用双线性微分方程组表示,描述不同脑区神经活动之间的效应连接关系,如式1-6-28所示:
其中, z=[ z 1, z 2,…, z n] τ为状态变量,表示 n个脑区的神经活动信号,是不可直接观测的; 表示 z对时间 t的导数,即神经活动在单位时间内的变化; u=[ u 1, u 2,…, u m] τ表示外界输入的 m个刺激; A矩阵代表脑区之间的固有连接(intrinsic connectivity),表示实验过程中脑区之间的平均连接强度; B j代表第 j个刺激对脑区间连接的调控作用; C矩阵表示外界输入的刺激对神经活动的直接影响。 A、 B、 C均是常数矩阵。以三个脑区为例(即 n=3),假定有两个刺激输入 u 1和 u 2, u 1直接作用于第一个脑区, u 2对脑区2到脑区3之间的连接以及脑区3到脑区1之间的连接起调控作用(图1-6-39)。
图1-6-39 动态因果模型脑网络连接示意图
则式1-6-28可以写成式1-6-29:
需要注意的是,DCM对效应连接的刻画描述的是一个脑区的神经活动对另一个脑区神经活动变化(数学上表示为神经活动的导数)的影响,因此,DCM更强调神经活动的“动态性”,其效应连接值的单位为s -1,即频率(Hz),连接强度越强,表示一个脑区的神经活动所引起的另一个脑区的神经活动的变化越快。
上述效应连接模型是建立在神经活动的基础上的,而神经活动并不能被直接观测到。被直接观测到的fMRI信号反映的是神经活动所引起的血流动力学变化。因此,DCM又引入血流动力学状态方程,用以描述神经活动和血流动力学变化(即fMRI信号)之间的关系。DCM采用Balloon模型,通过血管舒张信号、血流量、脱氧血红蛋白含量等参数来刻画fMRI信号的生成。
将上述神经状态方程和血流动力学状态方程相结合,在贝叶斯理论框架下对模型中的参数进行估计,最终得到效应连接系数矩阵A、B和C。
由此可见,使用DCM进行效应连接分析,需要事先假定脑区之间的连接关系以及外界刺激输入如何对该网络系统进行调控,基于该假定的模型结构对其中的未知参数(包括连接系数矩阵)进行估计,再通过个体水平或组水平统计确定每个连接系数的显著性,进而最终确定脑区间效应连接的具体情况。但这一分析过程存在两个问题:第一,在实际研究中,哪些脑区之间存在连接以及连接的方向可能是未知的,很难事先假定其连接模型的结构;第二,模型未知参数的估计结果往往与事先假定的模型结构有关,一旦模型结构发生变化,估计得到的连接系数也可能发生变化。为此,Friston等人又进一步提出将DCM和贝叶斯模型选择(Bayesian model selection)相结合,来进行脑区间的效应连接分析,其具体思路为:首先设计出所有可能的DCM模型结构,并对所有的DCM模型进行估计;再通过贝叶斯模型对所有的模型进行评估,选择最优的模型及其所对应的参数估计结果作为最终得到的脑区间效应连接结果。其中的模型评估采用贝叶斯准则,综合考虑模型对观测数据的拟合度和模型的复杂度两个方面,拟合度高且复杂度低的模型为最优模型。
相比于其他fMRI数据效应连接分析方法,DCM的最大优势在于对脑区间效应连接的刻画直接建立在神经活动的基础上,因此更符合实际情况。而且,DCM对神经活动与fMRI信号之间关系的刻画也考虑了不同脑区甚至不同刺激条件下的血流动力学响应函数可能不同的情况。但由于DCM模型复杂,模型参数估计的计算复杂度也相对较高,因此,在实际应用中也存在一定的局限性。例如,DCM很难处理大尺度网络(网络内脑区个数过多或连接过多)的情况,难以用于全脑网络分析;一般情况下,DCM只能处理脑区个数不超过8个的网络。而且,经典DCM属于确定性状态方程模型,必须有外界刺激输入的扰动,整个系统才能处于运行状态(即:每个脑区的神经活动才能发生变化),因此,仅适用于任务态fMRI数据的分析。但DCM作为一种效应连接分析的理论框架,自出现以来一直有研究者对其经典模型不断进行改进和完善。例如,多状态变量的DCM、用于静息态数据的DCM,基于随机性状态方程的DCM等。而且,这一方法还进一步扩展到了脑电和脑磁数据的分析中。
六、脑网络分析方法
(一)脑复杂网络简介
人脑是一个庞大、复杂的系统,据估算,一个健康成年人的大脑中神经细胞数量约为10 11个,其相互之间联系的突触数量则高达10 15个,形成了一个高度复杂的系统。该系统具有精细的解剖结构和复杂的功能连接,在结构上紧密连接,在功能上既相互独立又相互影响。
通过构建脑复杂网络,可以帮助我们更好地理解组成大脑的信息处理单元以及他们之间的相互作用,不再将大脑视为大量离散的解剖单元的集合体,而是由彼此纵横交叉相互连接的神经细胞构成的复杂统一体,从而更加全面细致地了解大脑内部极为繁杂的结构组织模式和功能动态模式,从更深层次来理解大脑的信息处理和认知加工机制。
美国著名复杂脑网络分析专家Sporns教授在 Networks of the Brain一书中谈到了网络分析可能解决的若干大脑机制问题:①通过对人脑结构网络的构建和分析来揭示大脑结构的构筑原理;②可以在人脑网络的背景下研究单个组成单元的活动;③感觉通道内部以及不同感觉通道之间信息整合的网络,可以用来揭示大脑的分布式网络加工机制;④神经纤维和通路的内在网络结构可以用来研究大脑静息态的动态活动模式;⑤有感觉输入或者任务态下的特定大脑活动模式可以用来研究人脑网络的动态扰动效应;⑥在人脑网络背景下研究人脑结构损伤可以解释其造成的功能损伤后果以及可能的恢复和代偿反应;⑦个体间大脑网络的差异可以用来研究个体间的行为与认知差异;⑧随着个体发育而逐步发展的人脑网络可以用来解释个体认知能力的发展;⑨大脑-身体-环境的交互对人脑网络的形成与发展有深远的影响。因此,从网络的角度来研究人脑的功能是极为必要的。
(二)脑复杂网络的构建方法
2005年,Sporns教授指出对人脑连接组的研究可以从3个空间尺度上进行,即:微观尺度、介观尺度和宏观尺度,分别代表神经元、神经元集群和脑区这三个尺度水平。但鉴于现有的成像技术手段,目前的研究主要集中在宏观和介观尺度水平上,通过结构磁共振成像、弥散磁共振成像等成像技术来构建大脑结构网络,采用脑电图(electroencephalogram,EEG)、脑磁图(magnetoencephalography,MEG)和fMRI等技术构建大脑功能网络,然后结合基于图论的复杂网络分析方法,揭示其拓扑原理,进而理解大脑的工作机制。
图论是目前复杂网络分析领域最主要的数学工具,是应用数学的一个分支。图论是将网络描述为由一组相互连接的点组成的图形。利用图论方法可以分析脑网络组织架构的内在特性,其网络组织架构特性决定了网络内信息传递的方式和效率,进一步决定了大脑的结构和功能属性。因此,复杂网络分析方法为全面理解大脑的结构和功能开创了新的思路。
图的基本元素为节点和边,其定义根据成像技术和描述信息的不同而不同。其常见的定义如表1-6-3所示:
表1-6-3 不同脑网络节点和边的定义
数学上可以用一个方阵(行数和列数相等的矩阵,记为A)来表示一个网络。矩阵中的行或列对应网络中的节点,其行数及列数即为网络中的节点数。矩阵中的元素对应网络中的边。例如,矩阵第 i行第 j列的元素(记为 A ij)代表节点 i到节点 j的连接。
根据节点之间的边是否有方向可以将网络分为有向网络和无向网络:①无向图(即无向网络)指节点之间的边没有方向,则其对应的矩阵为对称矩阵: A ij= A ji;②有向图(即有向网络)指节点之间的边具有方向性,则其对应的矩阵为非对称矩阵: A ij≠ A ji。
根据边的取值为二值(0或1)还是连续值,还可以将网络分为二值网络和加权网络:①二值网络指对节点之间仅定义有或无连接,而不定义连接的强度。节点之间有边或无边,其对应的矩阵元素取值为1或0;②加权网络指节点之间的连接有强弱之分,其对应矩阵元素的取值为连接强度。若设定一定的阈值,当连接强度大于或等于该阈值时,保留该边,当连接强度小于该阈值时,删除该边,则可将加权网络转换为二值网络。
(三)脑复杂网络的拓扑属性分析方法
人脑网络的拓扑属性反映了网络的架构特点,拓扑属性的不同决定了网络的结构稳定性、信息传递效率等方面的不同。网络拓扑属性可以用诸多指标来描述,下面以无向图为例,对常见的几种网络拓扑属性指标进行简要介绍:
首先假定网络中节点个数为 N, i和 j为节点的标号。
1.度(degree)
用来度量网络中节点间连接的密集程度。节点 i的度 k i定义为与节点 i直接相连的边的个数(也即节点的个数),与节点 i直接相连的节点称为节点 i的邻居。网络的平均度定义为网络所有节点的度的均值(式1-6-30):
2.聚类系数(clustering coefficient)
用来度量网络中节点的聚集程度。节点 i的聚类系数 C i定义为节点 i的邻居间实际存在的边数除以邻居间最大可能存在的边数(式1-6-31):
其中, E i为节点 i的邻居节点之间实际存在的边数。
网络的平均聚类系数为所有节点的聚类系数的均值:
3.最短路径长度(shortest path length)
用来度量节点之间信息传递的效率,最短路径长度越短,信息传递效率则越高。节点 i和 j之间的最短路径长度d ij定义为从节点 i到节点 j所需通过的最少的边数。网络的特征路径长度(characteristic path length)定义为网络中任意两节点的最短路径长度的平均(式1-6-32):
上述的度、聚类系数和最短路径长度的基本概念可以用图1-6-40的具体实例来说明。
图1-6-40 图论示例
图中圆圈为节点,连接节点的线为边,节点和边构成图;灰色节点为节点 i的邻居节点,节点 i的度为5(即为邻居节点的个数);节点 i和节点 j的最短路径长度为3(加粗的边为 i和 j之间的其中一条最短路径)
4.全局效率(global efficiency)和局部效率(local efficiency)
与最短路径长度类似,也是用来度量节点之间信息传递的效率。但最短路径长度在某些特殊情况下的应用具有局限性,如:若网络中两节点 i和 j之间互不连通,会导致这两个节点间的最短路径长度无穷大,从而为计算网络特征路径长度带来困难,进而提出了全局效率和局部效率指标。全局效率定义为网络中任意两节点的最短路径长度的倒数的平均(式1-6-33):
网络节点 i的局部效率定义为节点 i的邻居所构成子图 G i的全局效率(式1-6-34):
5.紧密中心性(closeness centrality)
是用来度量节点在网络中所处地位的重要性。节点 i的紧密中心性定义为从节点 i到所有其他节点的平均最短路径长度的倒数(式1-6-35):
6.介数中心性(betweenness centrality)
是用来度量节点在网络中所处地位的重要性的另一个指标。节点 i的介数中心性定义为网络中通过该节点的最短路径的条数占网络中所有最短路径条数的比例(式1-6-36):
其中, i、 j、 k分别为三个不同的节点, σ jk为节点 j和 k之间的所有最短路径的条数, σ jk( i)为节点 j和 k之间的通过节点 i的最短路径的条数。
7.小世界属性(small-worldness)
小世界网络同时具有较高的聚类系数和较短的特征路径长度,被认为在功能分化和功能整合之间具有最佳的平衡。与随机网络相比,小世界网络具有更高的聚类系数和相似的特征路径长度,即式1-6-37:
其中下标net为真实网络,下标random为随机网络。
进一步可定义σ=γ/λ,若 σ大于1,则认为网络具有小世界属性。 σ值越大,小世界属性越强。1998年,Watts和Strogatz提出人脑网络存在小世界属性。
8.无标度(scale-free)网络
其网络节点的度分布服从幂律分布,即: P( k)~ k - α,其中, k为节点的度, P( k)表示度为 k的概率。该分布的形状在双对数坐标下通常近似一条直线(图1-6-41)。
无标度网络的特点是网络中存在少数一些拥有众多连接的核心节点(hub),而大部分节点则只有少数连接。当无标度网络受到随机性攻击时,由于核心节点受到攻击的概率很小,而非核心节点受到攻击时,不会影响网络的连通性和完整性,因此无标度网络通常具有较强的鲁棒性。
上述各种网络拓扑属性指标既相互区别,又相互联系。例如,最短路径长度与全局或局部效率成负相关,最短路径长度越短,全局或局部效率越高;聚类系数与局部效率成正相关,聚类系数越大,局部效率越高,其网络防御攻击的能力也越强;网络的聚类系数越大,特征路径长度越短,则小世界属性越强。
图1-6-41 幂律分布
(四)脑复杂网络分析应用举例
目前已有许多解剖和生理学研究揭示了人脑网络的一些基本特征:局部脑网络连接,通常位于距离几百微米的神经元之间,且距离近的神经元多具有类似的响应特性;而在较长距离上延伸的神经连接,则沟通位于不同皮层区域的神经元之间的信息。脑网络连接在局部神经元集群的特定功能中起着重要作用,还可将不同信息来源进行整合,进而决定行为及认知状态。近年来研究者们基于不同的成像方式,尝试构建正常人与脑疾病患者的脑复杂网络并分析其拓扑属性,取得了很多重要发现,举例如下。
1.基于结构磁共振成像技术构建的复杂脑网络
2004年,Sporns等人研究了猴子和猫的大脑皮层的结构网络拓扑属性,证实了这几种网络都具有“小世界”属性。2007年,贺永等人利用正常人的结构像数据,通过分析大脑皮层54个脑区皮层厚度间的相关性,构建了人脑结构网络,并验证了该网络具有“小世界”属性。随后在2008年,他们分析比较了阿尔茨海默病患者与正常被试基于皮层厚度相关性度量的脑结构网络,发现患者脑网络的聚类系数和特征路径长度都显著增大,“小世界”属性减弱,表明该类患者的脑网络拓扑结构遭到了破坏。
2.基于弥散磁共振成像技术构建的复杂脑网络
弥散磁共振成像数据可检测神经连接的白质纤维束,利用纤维的连接数目、密度、强度、概率等可以定义大脑的解剖(白质纤维)连接网络,进而实现无创地重建个体人脑的白质纤维网络。但现有的纤维束追踪方法仍然存在一定的局限性,尤其是当纤维出现交叉时,纤维跟踪容易出现错误,导致连接丢失或出现伪连接。尽管如此,前人已经试图利用该技术构建复杂脑解剖网络并对其拓扑属性进行分析,得到了很多有意义的结果。
例如,2007年,Hagmann等利用2名正常被试的弥散磁共振成像数据,构建了基于个体的大脑解剖网络,并验证了该网络具有“小世界”属性;此外,也有研究探索基于弥散磁共振成像数据所构建的复杂脑网络的核心脑区(hubs)、模块化结构分区等网络拓扑属性与性别、年龄及智商等因素的关系。一些疾病相关的研究发现精神分裂症患者基于弥散磁共振数据构建的复杂脑网络的“小世界”属性较正常被试明显降低,癫痫患者脑网络出现随机化倾向。
3.基于功能磁共振成像技术构建的复杂脑网络
2005年,Salvador等人构建了正常被试静息态fMRI大脑功能网络,采用AAL图谱定义网络节点,脑区时间序列的偏相关系数定义网络连接,进一步分析脑功能网络拓扑属性,发现其具有“小世界”属性、较高的信息传递效率、优化的连接结构以及较高的拓扑稳定性。也有研究者研究年龄增长对网络模块性的影响,不同频段大脑功能网络拓扑属性的差异等。
大量基于静息态和任务态fMRI的研究,对失连接相关神经疾病脑网络拓扑属性的异常进行了探讨。例如,2008年,刘勇等人采用静息态fMRI研究精神分裂症患者的脑功能网络,发现其聚类系数降低,特征路径长度变长,“小世界”属性显著降低,表明精神分裂症患者的脑功能网络拓扑结构恶化,信息交换出现紊乱。而且,该研究还发现脑功能网络的拓扑属性和临床指标之间是显著相关的,患者患病时间越长,信息的传递效率越低,这为进一步理解疾病的发病机制及早期诊断提供了新的线索。
4.基于脑电图和脑磁图技术构建的复杂脑网络
2004年,Stam采用脑磁图构建了健康被试静息状态下不同频段内的大脑无向功能网络,发现低频(<8Hz)和高频波段(>30Hz)的脑网络都具有“小世界”属性,而中间频段(8~30Hz)的脑网络结构则接近于规则网络,不同频段脑网络结构间的差异表明不同频段的大脑活动可能对应着不同的大脑功能和不同的信息处理功能。此外,也有研究者研究受教育程度对脑网络拓扑属性的影响、睡眠的不同阶段下的脑网络拓扑属性的差异以及正常人与精神疾病患者脑网络拓扑属性的不同。
综上所述,磁共振成像和脑电图、脑磁图等无创影像技术使我们能够在大尺度上构建人脑结构和功能网络,而图论等复杂网络分析方法则为研究人脑网络的拓扑属性提供了强有力的工具。复杂脑网络的研究促进了我们对大脑的神经生理机制,信息处理方式以及各种认知功能的工作基础的理解。此外,脑网络分析方法也可以应用于不同类型的神经精神疾病(如阿尔茨海默病、精神分裂症、儿童多动症等)的研究中,通过探索由疾病导致的脑网络拓扑结构的异常变化,在系统水平上为揭示脑疾病的病理生理机制提供新的线索,并在此基础上建立描述疾病的脑网络影像学指标,为患者的早期诊断和个体化精准医疗等提供重要的辅助工具。
七、多变量模式分析技术
(一)多变量模式分析简介
随着信息科学与人工智能的高速发展,机器学习技术开始受到神经影像领域研究者的重视,并在功能影像数据分析中得到越来越多的应用。目前在功能成像领域应用最为广泛和成功的机器学习技术之一是多变量模式分析技术(multivariate pattern analysis,MVPA)。通常情况下,变量即为医学图像中的体素,因此,也有人将该方法称为多体素模式分析(multi-voxel pattern analysis,MVPA)。与传统分析方法相比,多变量模式分析方法在信息检测方面更为敏感。图1-6-42对传统激活分析方法和多变量模式分析方法做一简单比较。
前面介绍的fMRI数据分析中传统的激活分析方法更加注重检测某种实验刺激所诱发的某个脑区整体的激活。因此,在传统的激活分析方法中,数据在预处理过程中需要进行空间平滑,即,将相邻体素的神经活动进行加权平均,以去除噪声,并突出相邻体素神经活动公共的部分。
而MVPA技术则更加注重检测脑区内部多个体素的神经信号所构成的更为精细的空间分布模式,而这种神经活动的空间分布模式可以认为在一定程度上反映了神经编码信息。研究表明,在检测两种不同实验条件下的神经活动差异时,多变量模式分析技术比传统的激活分析方法更为敏感。其原因主要有两点:①当某个脑区在两种实验条件下的激活整体水平相似,但内部神经活动空间分布模式不同时,传统的激活分析方法则无能为力,而只能用MVPA来检测;②传统的激活分析多采用大批量逐体素的分析方法,即对全脑所有的体素逐一进行分析。这就使得对每个体素都需进行一次统计分析,而全脑所有体素的分析则不可避免地面临严重的多重比较校正的问题。而MVPA方法则是同时考虑所有体素,分析其所构成的空间分布模式,不需要进行多次统计分析,因此,避免了多重比较校正问题。
(二)多变量模式分析方法原理及流程
多变量模式分析方法的基本原理及流程示意图如图1-6-43所示:
用于脑影像数据分析的多变量模式分析方法最常见的有分类和回归两种形式。两者的主要区别在于数据样本的属性是离散型还是连续型:当数据样本的属性取值为离散型时(即,所有数据样本可以被归为几类),采用分类形式的MVPA;而当数据样本的属性取值为连续型时,使用回归形式的MVPA。下面分别进行介绍。
图1-6-42 单变量激活分析与多变量模式分析方法比较
A.当两种实验条件下感兴趣内的体素信号强度的空间分布不同但平均信号强度相同时,单变量分析无法检测两种条件脑激活的差别,而MVPA则可检测出这种差别;B.MVPA的基本思想是将多体素的样本点投射到高维空间,在任意维度(即任一体素)两个条件的信号均无明显差别,但在高维空间则可很容易将两个条件的样本点分开
分类形式的多变量模式分析是通过是否能够成功预测某未知数据样本的类别来判断不同类别的数据是否存在差异。若某分类器模型能够成功预测某未知数据样本的类别(分类正确率显著高于随机分类正确率),则说明不同类别的样本之间是存在差异的。对于二分类问题,随机分类正确率通常为50%。例如,某fMRI实验有听觉刺激和视觉刺激两种实验条件,每种实验条件下均采集了某个脑区的若干fMRI数据样本。若某分类器模型能够以80%的正确率(并通过统计检验显著高于随机分类正确率)成功预测这些数据样本究竟是在听觉刺激条件还是在视觉刺激条件下获取的,那么,我们就可以做出判断:听觉刺激和视觉刺激在该脑区所诱发的脑神经活动模式(具体体现为该脑区所有体素的fMRI信号所构成的空间分布模式)必然不同。换言之,如果听觉刺激和视觉刺激在该脑区所诱发的脑神经活动模式完全相同,则任何分类器都无法正确预测这些数据样本所属的类别。
回归形式的多变量模式分析则是当数据样本的属性取值为连续型时,通过多元回归模型,基于某数据样本的信号空间分布模式来预测其对应属性的具体取值。若能够基于数据样本的信号空间分布模式准确预测其对应属性的具体取值,则可以认为,数据中包含了该属性的编码信息。例如,我们需要知道某脑区的神经活动是否编码了声音信号的频率信息,即:不同频率的声音刺激是否会诱发该脑区产生不同的神经活动模式。我们就可以通过是否能够基于该脑区内体素的fMRI信号所构成的空间分布模式来准确预测声音信号的频率,来判断该脑区的神经活动分布模式是否编码了声音的频率信息。预测的准确性可以通过预测误差来确定:若预测误差显著小于随机预测误差,则说明数据中包含属性值(如:声音的频率值)的有用编码信息。
1.预测模型的训练与测试
图1-6-43 多变量模式分析方法的基本原理及流程示意图
不论是分类形式还是回归形式的多变量模式分析,均需要首先得到一个分类模型或者回归预测模型,然后计算该模型对样本属性类别或属性取值的预测准确性。那么,如何获得分类模型或回归预测模型呢?以分类形式的MVPA为例,首先需要利用一批已知类别的样本(训练样本)来训练分类器模型或回归预测模型,这一过程称为训练过程。在训练过程中,数据样本和其对应的属性值均是已知的,每一个数据样本都可以表示为高维空间中的一个点,空间的维数即为数据样本所包含的特征个数。要训练一个分类器模型,即是要在该高维空间中寻找一个分界面,使其能够将不同类别的样本最为准确地分开,即尽可能使不同类别的样本位于分界面的两边。一旦找到这一分界面,即确定了分类器模型,下一步需要对这一训练好的分类器模型进行测试。测试的目的是为了判断该分类器是否能够准确预测未知数据样本的类别。因此,在测试过程中,我们会提供给分类器若干数据样本,这些样本所对应的真实类别对分类器是不可见的,也就是说,分类器所能够利用的信息只有数据本身,分类器将根据这些数据样本位于分界面的哪一侧来预测这些样本所对应的类别。通过对比分类器预测的类别与样本的真实类别,即可得到分类器的预测正确率。测试过程使用的数据样本称为测试样本。需要特别强调的是:测试样本必须是一组新的样本,即测试样本不能与训练样本有任何重叠。对于回归形式的MVPA,其训练和测试的基本思路与上述过程是一样的,唯一的差别在于训练过程是为了在高维空间中用一条直线(对线性模型而言)或一条曲线(对非线性模型而言)去拟合所有的训练样本,使其拟合误差达到最小。然后再将测试样本数据代入训练得到的回归模型,即可计算出该测试样本的属性预测值。
2.交叉验证
交叉验证是指将所有数据样本分成若干份,每次只拿出一份作为测试样本,其余数据作为训练样本,每次均可得到一个分类正确率或预测误差。将这一过程重复很多次,使得每份数据均作过一次测试样本。然后将得到的所有分类正确率或预测误差进行平均,作为最终的分类正确率或预测误差。这种交叉验证的方法称为“留一法”。由于通常情况下影像学实验采集到的数据非常有限,因此为充分利用数据,很多研究都采用“留一法”交叉验证。
3.预测正确率的统计显著性检验
如何判断分类模型分类正确率是否显著高于随机分类正确率呢?由于分类正确率通常表达为百分数,其取值范围为0%~100%,不服从正态分布,因此,常见的 t检验等参数统计检验方法均不太适合。目前常用的统计检验方法为置换检验。其具体过程为:首先随机打乱训练样本的类别,即将每一个训练样本随机分给某个类(原本属于A类的训练样本会被随机分配给A类或者B类,而原本属于B类的训练样本也会被随机分配给A类或者B类),但仍然保证每类内的训练样本数量保持不变。这样,不同类别的训练样本就会被混在一起而被重新分为两组,使得训练样本丧失其真正的类别信息。然后,使用这组已经打乱了类别信息的训练样本再去训练分类器,得到的分类器就不应再包含真实情况下不同类别之间的差异信息。再用该分类器对原来的测试样本进行测试,就会得到一个随机类别情况下的分类正确率(即随机分类正确率)。将上述过程重复多次(如1万次),就可以得到1万个随机分类正确率,进而构建出随机分类正确率的直方图(即:置换检验中的零分布)。通过该零分布,即可计算出真实的分类正确率所对应的 P值,从而确定其统计显著性。 P值的计算举例如下:若1万次置换检验中有5次得到的随机分类正确率大于或等于真实的分类正确率,则 P=5/10 000=0.000 05;若1万次置换检验得到的所有随机分类正确率均小于真实分类正确率,则 P<1/10 000=0.000 01。对于回归形式的MVPA置换检验的过程类似,首先将训练样本的属性值打乱,然后用已经丢失有用信息的训练样本去训练回归预测模型,再使用该预测模型对原来的测试样本的属性值进行预测,进而得到随机预测误差。重复多次后,可以得到随机预测误差的直方图。再通过检查有多少次随机置换得到的随机预测误差小于或等于真实的预测误差,进而计算出 P值。
4.“探照灯”(searchlight)MVPA
上述MVPA的分析过程是基于感兴趣区的,即需要事先选定所要研究的区域(也可将全脑作为一个感兴趣区)。当希望对某一特定感觉或认知功能相关的神经活动在全脑内进行空间定位时,可以采用“探照灯”MVPA的方法。所谓“探照灯”MVPA是针对每一个体素均进行一次MVPA,每次MVPA是基于以该体素为中心,R为半径的小球所包含的所有体素构成的一个局部区域为感兴趣区,再将得到的分类正确率赋予该中心体素。将这一分析过程对全脑所有体素重复一遍,可得到每个体素(即该体素为中心构成的小球区域)的分类正确率,从而完成对全脑所有局部区域的信息搜索。探照灯区域的半径R需要根据实际问题进行选择。半径越大,包含的体素越多,能检测到感兴趣神经活动的概率越大,但空间定位的精度也会越低。由于“探照灯”MVPA是对每个体素均进行一次分析,因此,在对分类正确率进行统计检验时,需进行多重比较校正。
(三)多变量模式分析应用举例
很多研究已经证实了多变量模式分析在检测神经活动方面比传统方法更为敏感。例如,经典神经科学理论普遍认为初级感觉皮层是单模态脑区,即仅对其对应感觉模态的刺激进行反应(例如初级视觉皮层仅对视觉刺激进行反应),而基于fMRI数据的传统脑激活研究也证实了特定模态的感觉刺激通常仅激活其对应的感觉皮层,而其他初级感觉皮层鲜有激活。然而,利用多变量模式分析方法对不同感觉模态刺激(如视觉、听觉和触觉)条件下初级感觉皮层的fMRI信号的空间分布模式进行分析,研究者发现任一初级感觉皮层能够对非对应的感觉模态刺激产生反应(例如,初级视觉皮层也会对听觉和触觉刺激产生反应),说明初级感觉皮层并不是严格的单模态脑区。利用类似的方法,也有研究者利用fMRI数据,辨识出能够编码疼痛强度信息的神经活动模式,从而可以实现对疼痛感知强度的预测。此外,基于深度学习算法构建的分类器可以根据从6~12个月大婴儿的磁共振图像中提取出的脑表面积信息成功预测婴儿在两岁时患自闭症的风险。这一研究结果不仅说明自闭症儿童的脑结构改变远早于其异常行为的出现,而且为自闭症的早期预测提供了一种可行的方案。还有研究利用脑区灰质体积信息成功区分精神分裂症、抑郁症和双相障碍患者,为精神疾病的客观诊断提供了工具。
(秦 文 梁 猛)