1. 引言
功能磁共振成像(fMRI)用来研究大脑功能已经有十余年了 [1] [2] 。虽然fMRI相对于其他的成像方法,在空间和时间分辨上都占有一定的优势,但它的数据量很大。所以fMRI的数据处理方法一直是研究的热点。一般说来,用于分析功能核磁数据的方法可以分为两大类:模型驱动和数据驱动 [3] 。基于模型驱动的方法需要事先根据实验设计假设大脑激活的时间模式,而数据驱动方法就不需要这种假设,例如,主成分分析(Principal Component Analysis PCA) [4] ,独立成分分析(Independent Component Analysis ICA) [5] [6] 和模糊簇分析方法(Fuzzy Clustering Analysis) [7] [8] 都是数据驱动方法。然后这些数据驱动的方法在实际应用中都存在一定的局限性。
TCA方法的出现是为了得到大脑反应的时间和空间信息,利用这种方法第一次得到了被试服用葡萄糖之后大脑出现的时间峰,并进一步得到了大脑激活的空间信息 [9] 。这种方法称为原始时间簇分析方法(original temporal clustering analysis, OTCA)。TCA的基本假设是,大脑图像中的所有像素,对于每个像素来讲,只认为达到最大值的时间点是激活点而其他没有达到最大值的时间点是非激活点也就是参照点。对于在每个时间点有m ´ n个空间像素(m行,n列),而对于每个像素有p个时间点的一层fMRI数据可以表示成二维的矩阵S,该矩阵有m ´ n行和p列:
(1)
在方程(1)中,矩阵元素表示的是第i个空间像素在第j时刻的像素灰度值。
然后由矩阵S通过以下的公式就可以得到矩阵W:
(2)
在方程(2)中,代表第i个空间像素在整个时间变化过程中在第j时刻达到最大值的像素。然后把矩阵W中对应的列相加,就可以得到一维的时间向量K表示如下:
(3)
结果的每一个像素代表每个时间点达到最大值的所有像素的个数。
从而,画出向量K对应于时间点的曲线,在对大脑的时间反映模式完全未知的情况下就可以通过该曲线看出激活峰的时间窗。
在方程(1)中,矩阵S的行数代表的是空间像素的个数,其中包括脑外像素和脑内像素。然而,脑外像素和脑内非激活像素的存在会影响TCA的检测灵敏度。对于脑外像素比较容易去除,因为脑外像素的灰度值比脑内像素的灰度值低几倍,利用阈值法即可基本把脑外像素去除。但对于脑内非激活像素则要复杂得多,因为脑内像素的灰度值相差不大。但脑内组织可以分为三种,即灰质、白质和脑脊液,它们所对应的像素灰度值是有差别的,利用这一特性可以将核磁共振图像中的三种成分分开。SPM软件中已经包含了这种功能,可以利用SPM软件将核磁共振图像分割成灰质、白质和脑脊液图像。因为大脑的激活区理论上出现在灰质区域,因此只需对分割后的灰质图像进行TCA分析,这样就可以去除白质和脑脊液中所对应的像素对TCA的影响,使得TCA方法也可以用于全脑数据。而且我们用针刺合谷穴的fMRI功能数据进行验证改进TCA方法的可行性。
在改进的TCA方法中,先用SPM2软件(Statistical Parametric Mapping统计参数图,SPM是在数学软件包Matlab平台上开发的,用SPM进行数据处理分析的过程主要分为两大部分:预处理过程和统计分析过程)将每个被试的第一幅核磁共振图像分割成灰质、白质、脑脊液图像(如图1(1-1)~图1(1-3)所示)。其中白质和脑脊液所包含的像素是与神经活动无关的噪声成份,所以只需考虑灰质内所包含的像素。由SPM分割后得到的灰质直接就可以把脑外像素去除,但脑内灰质外的空隙也就是白质或脑脊液所在位置的这些像素灰度值并不为零,不过这些像素的灰度值相对于灰质内的像素的灰度值要小很多,利用阈值法即可把这些像素去除。经过多次实验发现,如果取灰质内灰度值最大像素的七分之一作为阈值去除脑内位于白质和脑脊液的像素,改进TCA方法的灵敏度是最好的。其中比阈值小的像素都被认为是灰质外的像素并赋值为0,其余剩下的就是灰质所包含的像素并赋值为1。这样每个被试都有一个归一化的灰质图像。然后把每个被试的归一化灰质图像与这个被试的所有功能图像相与得到的就是只包含灰质的图像。这样就避免了脑外像素和脑内白质和脑脊液所包含的非激活像素对TCA检测灵敏度的影响。通过这样的处理,得到的是每个被试都是只包含灰质的功能图像。然后用TCA算法对这样的功能图像进行处理就可以得到相应的时间曲线。
2. 实验与数据处理
本研究有十名健康自愿者参加了实验,男6例,女4例,年龄(26.8 ± 3.6)岁,右利手,无精神神经系统疾病史。被试都没有神经或精神病史。志愿者被告知他们将接受一次针灸刺激。所有被试在实验前被告知实验内容,签署了知情同意书,并被告知有权在实验的任何时刻退出,实验的设计方案如图2所示。每个被试只经历1个block。
实验在西门子1.5特斯拉全身核磁扫描仪上进行(Sonata,西门子,德国),使用标准的头部线圈。20层平行于AC-PC线的图像覆盖全脑。功能图像使用T2*加权梯度回波EPI序列采集,平面分辨率3.44毫米(TR 3000 ms,TE 50 ms,翻转角90˚度,FOV为220 × 220 mm,矩阵64 × 64,6 mm层厚,1.2 mm层间隔)。
自愿者在平卧位、视听封闭状态下参加实验,用黑色眼罩遮盖双眼,海绵耳塞堵塞双耳,针刺右合谷穴(握拳在第二掌指关节与第一掌骨腕端连线的中点)实验。62次静息状态的扫描后,一根直径0.30 mm
Functional images 1-1 GM 1-2 WM 1-3 CSF
Figure 1. Functional image and corresponding gray matter, white matter and CSF
图1. 功能图像及相应的灰质、白质、脑脊液
Figure 2. After a 60 seconds resting state, a sterile silver needle was inserted and lasted 60 s, then followed by a rest of 280 s
图2. 静息60 s秒后,进针,刺激持续60秒后拔针,接着是280秒的静息
长25 mm的银针被插入并旋转60个扫描点,然后,针被拔出并且继续扫描直到全部402个扫描点完成。所有的针灸操作由同一位针灸专家实施。针灸过程中,采用“平补平泻”手法以1Hz的频率顺时针和逆时针捻转。
首先用统计参数图软件(SPM2, http://www.fil.ion.ucl.ac.uk/spm/) [10] 对所有被试的功能图像进行头部校正。为了避免刚开始扫描时磁场场强的波动,每个被试所扫描图像的前两幅图像在进行数据处理时被去掉,这样每个被试有400幅图像。为了校正在扫描期间图像与图像之间的头部移动,用最小二乘法把每个被试所得的功能图像都配准到所得的第一幅图像。在这个研究中,不管TCA或是改进TCA方法,取的都是每个时间点上达到最大值的像素个数。用这两种方法分别对预处理后的每个被试的功能数据单独进行处理,然后对所得的每个被试的时间曲线进行平均得出所有被试的平均时间曲线。这样就得到了每个被试的时间变化曲线和所有被试的平均时间变化曲线。
3. 结果
利用TCA和改进TCA两种方法分别求得了十个被试的时间变化曲线和十个被试的平均时间变化曲线。其中,十个被试的平均时间变化曲线如图3所示,而两个被试的时间变化曲线如图4所示。
Figure 3. The average over 10 subjects TCA signal. (a) is the detected TCA signals with modified TCA method. (b) is the time-series detected by TCA
图3. 两种方法检测的样本的时间簇分析的平均结果
Figure 4. Time series of two single-subjects detected by the two methods. (a) and (b) are temporal response of single-subject detected by modified TCA. (c) and (d) are temporal response of single-subject detected by TCA
图4. 显示的是用两种TCA方法得到的两个被试的TCA时间曲线
用TCA和改进TCA所得到的所有被试的平均时间变化曲线如图3所示,从图可以看出,由改进TCA所检测出来的大脑反应曲线和实验设计方案很接近。大脑在进针和拔针的两个时间点上都有很明显的激活峰出现,表明中枢神经系统对针灸有了反应。
对平均结果图进行时间窗分析,发现第一个峰的最大值发生在第63个采样点,第二个峰的最大值出现在第122个采样点。因为TR为3 s,所以进针7~9 s后血氧水平依赖性信号达到了最大值,拔针4~6 s后血氧水平依赖性信号达到另一个峰值。而从TCA所检测出来的时间曲线来看,就看不出有明显的激活峰存在。
改进TCA方法检测的结果与实验设计很相符,在进针和拔针时分别检测到了两个明显的峰,但是TCA方法检测出来的点很分散,没有与刺激相对应的峰。在这个研究中我们对功能图像做了一些处理并分别用TCA方法对经过处理和没有作处理的同一组针刺合谷穴的实验数据进行计算。对两种结果进行比较可以看出,在处理全脑fMRI功能数据时改进TCA比TCA更加适用。
4. 讨论
在以前的大多数研究中,TCA都局限于处理单层fMRI数据 [9] [11] 。为了解决TCA不能用于处理全脑数据的这个局限性,在改进TCA中,我们利用SPM2把功能图像分割成灰质、白质和脑脊液三部分,而我们只要灰质图像,另外通过阈值法把灰质空隙处像素灰度值不为零的像素去掉。这样剩下的像素就只有灰质内的像素,而在理论上激活也只存在于灰质当中。只让灰质参与最后的TCA计算,就能排除脑外像素和脑内非激活像素对TCA检测灵敏度的影响。
改进TCA所得的平均时间曲线(图3(a)),我们可以看出改进TCA方法能有效地检测出对应于进针和拔针两个明显的激活峰,而且所检测出的激活峰的位置也和实验设计很吻合。而由TCA所检测出的时间曲线中(图3(b)),就没有检测出明显的激活峰。
因为人和人之间有很大的差异存在,所以对单个被试进行统计分析也很有意义。Wu和Napadow在对多被试统计分析的结果和单被试统计分析的结果进行比较时发现两种结果之间有很大的差别 [12] [13] 。然而,因为TCA方法的灵敏度较低,在处理单个被试的数据时,许多情况下检测不出激活的存在,在以前的研究中,所给出的TCA曲线也都只是平均结果 [10] 。尽管Liu等人给出了单个被试的TCA时间曲线 [9] ,但他们处理的都是单层的大脑数据。在这个研究中,我们对TCA进行了改进并把这种改进的TCA方法用于针刺合谷实验的全脑数据。通过比较改进TCA和TCA处理单个被试的结果可以看出(图4(a)~图4(d)),改进TCA处理单个被试全脑数据的能力比TCA有很大程度的提高。
另外,由改进TCA检测结果可以看出(图3(a),图4(a),图4(b)),改进TCA不仅检测出进针、拔针和捻针过程中脑内有激活,而且还检测出在拔针之后大约10分钟后大脑也有明显的激活。这和传统中医理论里说的后效应相符合。但是,从图3(a)和图4(a),图4(b)也可以看出,由改进TCA检测出大脑在进针前的静息扫描期间也有明显激活。这可能是因为被试刚开始对实验中的针灸有所顾忌,心里紧张造成的。
总的来说,改进TCA在处理全脑数据是可行并有效的一种方法。本文并用一组针刺合谷的全脑数据证明了改进TCA的检测能力。
基金项目
本研究受广东省卫生厅基金资助(A2014280).