自闭症谱系障碍儿童静息状态下脑电微状态研究

张锁良 万灵燕 张志明 康健楠 李小俚 庞 姣*

1(河北大学电子信息工程学院,河北保定 071002)

2(河北省机器视觉技术创新中心,河北保定 071000)

3(北京师范大学认知神经科学与学习国家重点实验室,北京 100875)

引言

自闭症谱系障碍(autism spectrum disorder,ASD)是一种发展性障碍,具有语言表达能力受损、交往能力异常、兴趣范围狭窄和固定行为重复等特征,严重影响患者的社会交流、情感表达以及认知功能[1]。据美国疾控与预防中心统计,截止到2018年,每59 名儿童中约有1 名为ASD 患者[2],这一统计数据现已更新为每45 名儿童中约有1 名自闭症患者儿童。ASD 现已成为全世界范围内的亟待解决的医学难题。到目前为止,ASD 的发病机制还较为模糊,临床诊断也主要依靠医生的经验判断以及量表评测,具有一定的主观性和因此导致的误诊率[3]。儿童时期的个体是大脑神经发育和语言发展的关键时期[4],自闭症儿童在顶叶、枕叶、额叶、颞叶以及楔前叶等脑区表现出过度生长模式[5]。因此,了解ASD 儿童的发病机制和异常的大脑发育,为ASD 的临床诊断提供有效的客观指标,具有非常重大的研究意义和科学价值。

脑电微状态的概念最早由Lehmann 于1987年提出[6],被称为“思想的原子”[7],是一种成熟的基于电极阵列上的电位地形图。微状态的地形拓扑结构并不是随着时间而发生变化,是在80 ～120 ms内总体上保持相对稳定,然后又迅速变换为另一种在一个时间段内保持相对稳定的地形图拓扑结构,但在这个过程中,地形图的极性和强度可能会发生变化,而地形拓扑图极性的反转通常在微状态分析中被忽略[8]。瑞士科学家Lehmann 认为,在保持稳定的这段时间反映了大脑信息加工的基本过程,是人们思想意识组成的基本单位[6]。之前的研究显示,4～8 个微状态类别的脑电地形结构可以解释实验数据64%～83%的全局方差[9-10]。经过几十年的研究与探索,微状态分析技术已经逐渐走向成熟,Britz 等[11]确定了微状态A、B、C 和D,分别对应于语言处理、视觉网络、显著性网络和注意力相关的静息态网络。Custo[12]等使用TESS 方法,估算了微状态网络的源分布。现微状态分析方法已经被广泛应用于探究各类脑部神经疾病的病理机制。然而,微状态分析方法在ASD 中的探索少之又少,到目前为止,还没有利用这一分析方法针对ASD 儿童的研究报告。

本研究利用脑电微状态分析方法,探究ASD 儿童与正常儿童在静息状态下微状态在时间参数上的差异,根据这些参数差异表征功能意义,为ASD患者的临床诊断和治疗提供客观依据。

1 材料和方法

1.1 被试者的数据信息

本实验共招募28 名儿童,其中包括13 名ASD儿童(男孩10 名、女孩3 名,右利手,平均年龄(4.58±2.22)岁)和15 名正常儿童(男孩11 名、女孩4 名,右利手,平均年龄(4.47±1.23)岁)。所招募的ASD 儿童均符合以下标准:一是符合《精神病诊断与统计手册(第5 版)》和《孤独症诊断访谈量表修订版》标准；二是无头部损伤或者大脑疾病；三是无精神病史或服用精神药物；所招募的两组被试儿童在性别和年龄上均无统计学差异(P＞0.05)。本实验经过河北大学附属医院伦理委员会批准,所有被试者均由其监护人签署知情同意书。

1.2 数据采集与预处理

采用128 通道的EGI HydroCel Geodesic System(Eugen,OR,美国)脑电采集系统采集脑电数据。在实验过程中,保证每位儿童的采集时长为8 ～10 min,被试儿童只要静坐在无噪声的室内保持放松的状态,尽量减少面部表情和肢体动作。在采集过程中,所有通道的阻抗设置在50 kΩ 以下,采样频率为1 000 Hz,参考电极为Cz。

预处理过程首先采用0.5 ～40 Hz 的带通滤波器处理,将采集到的数据分为4 s 一段,手动地去除含有眼电、肌电、眨眼和工频等伪迹和噪声的时间段。根据10-20 国际标准导联系统选取62 个通道,然后对经上述处理过程仍存在较大伪迹的通道进行坏道差值处理,最后以对数据重参考作为平均参考,将得到的数据用于接下来的微状态分析。

1.3 脑电的微状态分析

在确定的时间段内,脑电拓扑结构依据采集电极上获取的信息,根据EEG 数据的62 通道的电压幅值,绘制出该时间段内头皮电势分布图(见图1)。地形图的拓扑结构总是在一定的时间段内保持稳定,继而转换为另一地形拓扑状态。其中,保持一段时间稳定的地形图被称作微状态(microstates)。微状态分析,是指从时域不同的脑电地形拓扑结构中,通过合适的聚类方法分割出几类最主要的地形图,用来表示整个EEG 数据的地形拓扑结构在时间尺度上的变换,每个脑电地形拓扑结构表示一个微状态类别,每个脑电地形图保持稳定的时间大约在80～120 ms[13]。然后将分割好的几类微状态类别拟合回EEG 数据,在每个时间点标记与之相关性最大的微状态类别。通过全局解释方差(global explained variance,GEV)[13],拟合微状态地形图到EEG 数据各时间点,其值越高拟合度越高。最终,在拟合后的原始数据中提取感兴趣的时间参数[14-15]。在不同组别条件下比较时间参数的差异,进一步分析其表征在生理意义上的差异。本研究所计算的时间参数为:每类微状态的持续时间、每类微状态的发生频率、每类微状态的时间覆盖率、微状态类别间的转移概率。

图1 62 通道分布Fig.1 62 channel distribution

1.3.1 原始地形图

首先需要选取质量较好的原始地形图进行聚类,以确保研究结果清晰可观。总体电场功率(global field power,GFP)是指在某一时刻头皮上所有电极电压的标准差,是一个与参考电极选取无关、可以描述各个时间点上地形图强度的指标。之前的研究发现,GFP 较高的地形图往往具有较高的信噪比和良好的稳定性[13]。因此,在本研究中,选取GFP 局部峰值处的地形图作为聚类的原始地形图。依据EEG 数据的电极数目和每个电极上的电压值,计算每个被试的GFP 值,并绘制曲线,选取GFP 曲线局部峰值处的地形图作为接下来的聚类原子,如图2所示。GFP 的计算公式为

图2 微状态分割Fig.2 Microstates segmentation

式中,ui是某一电极上的电压,¯u是各个采样点的电极平均电压,N为电极数目。

1.3.2 分类算法

本研究中采用的是原子化与凝聚层次聚类算法[16](atomize ＆ agglomerate hierarchical clustering,AAHC),除此之外,改进的k-means聚类算法[16-17]也经常在微状态类别分割中得以应用。k-means算法对于质心的选取是随机的,对于同一数据,每次运行的结果都带有随机性和不稳定性,只有通过大量增加迭代次数,才可使其分割结果达到较为准确的值。因此,本研究选取AAHC 聚类算法。

AAHC 算法首先是找到GEV 最小的地形图从所有地形图中分离(“原子化”),然后计算它与所有地形图的相关性,将相关性较大的与之进行合并,然后再进行相关性计算,取相关性较大的再进行合并,以此类推,最终所有的地形图会被聚为一类。假设选择的聚类数目为k,则这一过程在所有剩余地形图数目为k时选择结束,得到k个微状态地形图。此算法是一次自上而下的聚类过程,不需要多次重复,仅需一次迭代,就可以得到较为稳定的聚类结果。

上述过程是在单个数据水平上的聚类,根据每个被试数据确定出一组微状态地形图,随后将组内每个被试经第一次聚类确定的微状态地形图提取出来,再次使用上述的聚类算法进行二次聚类,得到组水平上的微状态地形图,也称为模板图,然后计算其与EEG 数据各个时间点上的地形图GEV 的相似度,并对每类微状态对应的时间点进行标记,这样原来的EEG 数据段在时域上就变成了各类微状态交替变换的时间序列。在本研究中,根据Cartool 软件[18]中提供的准则,并通过计算2 ～6 个微状态类别对于EEG 数据的GEV 的解释程度,确定了在组水平上的微状态类别。GEV 越高,得出的结果就更加准确。在2 ～6 个聚类中,当数目大于4以后,GEV 几乎不再增加,因此为了便于比较,本研究确定的聚类数目为4。这4 种微状态类别分别标记为微状态A、B、C 和D,微状态分析流程如图3所示,将采集到的原始脑电信号进行预处理,去除工频、眨眼、肌电等伪迹的干扰,得到可用于后续分析的EEG 信号:首先根据需求确定合适的聚类算法,再用上述提到的Cartool 软件中的准则以及GEV 值的高低确定微状态类别的数目；然后计算并绘制每个被试EEG 数据的GFP 波形,将波形中各个局部峰值处时间点的脑电地形图提取出来；用AAHC 算法根据确定的微状态类别数目,将这些地形图进行聚类,得到单个被试EEG 的微状态模板图。依次对所有被试重复上述过程,然后使用AAHC 算法对得到的所有被试EEG 的微状态模板图进行二次聚类,得到组水平上的微状态模板图；最后根据微状态地形图与原始时间点地形图的GEV 相似度进行拟合,得到4 类微状态交替变换的时间序列,提取感兴趣的时间参数,在不同组别之间比较差异。

图3 微状态分析流程Fig.3 Flow chart of microstates analysis

1.4 信息论分析与统计分析

基于上述微状态分析,将确定好的4 类微状态拟合回EEG 时间序列,可将其视为包含A、B、C、D 等4类微状态交替变换的离散时间序列。基于这一特征,本研究在时域上利用信息论的知识探究微状态类别间的转移是否具有独立性,以及在ASD 和TD 儿童中是否存在转移差异。通过计算0、1 和2 阶等低阶的马尔可夫模型,探究从当前的微状态转移到其他微状态类别是否与当前或者之前时间点所处的微状态类别有关。用Ii表示可能出现的微状态A、B、C、D 的类别,即Ii∈(A,B,C,D)；用Xt表示变化的时间点；P(Xt = Ii)表示t时刻微状态分布的概率。从微状态A 转移到微状态B 的转移概率表示为TAB=P(Xt+1=IB｜Xt=IA)。表达式P(Xt+1)=P(Xt+1｜Xt),表示在转移到下一个状态时依赖于当前状态,即0 阶马尔科夫模型；P(Xt+1｜Xt,Xt-1)=P(Xt+1｜Xt),表示在转移到下一状态时依赖前一个状态,即一阶马尔科夫模型；P(Xt+1｜Xt,Xt-1,Xt-2)=P(Xt+1｜Xt,Xt-1),表示转移到下一状态时依赖前两个状态,即二阶马尔科夫模型。另外,本研究通过各个微状态之间的转移概率,列出了所有被试的转移矩阵T,通过检验其对称性,探究从微状态M→N和从N→M(M,N∈(A,B,C,D))是否存在差异性[19],有

式中,PAA表示微状态A 转移到微状态A 的概率,PAD表示微状态A 转移到微状态D 的概率,以此类推。本研究不考虑转移到自身的这种情况,只讨论不同微状态类别之间的转移情况。在本研究的两组实验数据中,所有微状态类别的时间参数的统计分析均使用独立样本t检验,马尔可夫模型数据使用卡方检验,在SPSS 26.0 版本的软件中计算。

2 结果

对ASD 组以及TD 组分别进行上述微状态分析,得到两个组水平上的4 个微状态模板,如图4所示。两组微状态图的拓扑结构存在一些差异,不同颜色及颜色深浅代表采集电极所记录的电压正负和大小。

图4 组水平的4 类微状态地形图Fig.4 4 types of microstate topographic maps at group level

计算确定好4 类微状态地形图与原始EEG 数据的各个时间点的地形图GEV 相似性,并在每个时间点上标记与其相关性最大的微状态类别,原始EEG 数据在时间上变成4 类微状态在交替变换的序列。如图5所示,分别是将TD 儿童和ASD 儿童微状态类别拟合回原始数据的时间序列。可以看出,ASD 儿童在微状态类别之间转移的频率比较高,TD 儿童的微状态时间序列相对来说比较稳定,单个微状态类别的持续时间较长。

图5 两组微状态局部时间序列(每种颜色分别代表一个微状态类别，其中蓝色代表微状态A，橙色代表微状态B，黄色代表微状态C，紫色代表微状态D)。(a)ASD 组；(b)TD 组Fig.5 Local time series of two groups of microstates(Each color represents a microstate category， where blue represents microstate A， orange represents microstate B， yellow represents microstate C， and purple represents microstate D).(a)ASD group；(b)TD group

将4 类微状态拟合回EEG 数据以后,提取出感兴趣的时间参数来进行计算,本研究所选的时间参数为每类微状态的持续时间、发生频率、时间覆盖率以及转移概率。结果及统计分析如表1所示。由表可见,ASD 组的各微状态的持续时间均小于TD组的相应时间,微状态A、C、D 在两组中的持续时间更是存在显著性差异(P＜0.05)；ASD 组微状态A的时间覆盖率显著小于TD 组的相应状态,微状态B的时间覆盖率显著大于TD 组(P＜0.05)；ASD 组的4 个微状态的发生频率均显著大于TD 组的相应状态(P＜0.05)。

表1 时间参数的统计分析结果Tab.1 Statistical analysis results of time parameters

图6所示为两组中4 类微状态的持续时间、发生频率以及时间覆盖率的箱型图,可以直观地观察到这组数据的分布和离散情况,“箱子”越长表示该组数据分布越离散,“箱子”越短表示该组数据分布越集中。其中,(a)为两组中各个微状态的持续时间对比,ASD 组微状态A、B、C 的整体数值分布都远小于TD 组的状态,微状态B 的中位数值略小于TD 组的状态,但是不存在明显差异；(b)为两组中各个微状态的时间覆盖率对比,ASD 组微状态A 的数值分布整体小于TD 组的状态,微状态B 的整体数值分布大于TD 组的状态,微状态C 和微状态D在两组中的数值分布存在略微差异,但不够显著,ASD 组微状态D 的数据离散程度显著大于TD 组的状态；(c)为两组中各个微状态的发生频率对比。ASD 组的4 类微状态数值分布都显然大于TD 组的状态,并且ASD 组中的微状态A、B、C 数据分布的离散程度远大于微状态D 的状态。

图6 两组微状态时间参数差异箱型图(蓝色线条代表ASD 组，红色线条代表TD 组；箱型结构自上而下的横线位置分别为这组数据的上边缘、上四分位数、中位数、下四分位数及下边缘；+表示这组数据中的异常值)。(a)持续时间；(b)时间覆盖率；(c)发生频率Fig.6 Box diagram of the difference in time parameters between the two groups of microstates(The blue line represents the ASD group and the red line represents the TD group.The horizontal lines from top to bottom of the box structure are the upper edge， upper quartile， median， lower quartile and lower edge of the data.“+” indicates an outlier in this set of data).(a)Duration；(b)Time coverage；(c)Frequency of occurrence

4 类微状态的转移概率情况如表2所示。从统计结果来看,在ASD 组中,微状态A 时,转移到微状态B 的概率最大；微状态B 时,转移到微状态C的概率最大；微状态C 时,转移到微状态B 的概率最大；微状态D 时,转移到微状态C 的概率最大。在TD 组中,微状态A 时,转移到微状态C 的概率最大；微状态B 时,转移到微状态C 的概率最大；微状态C 时,转移到微状态A 的概率最大；微状态D 时,转移到微状态A 的概率最大。在两组之间A→B、A→C、A→D、B→A、B→C、C→A、C→B、C→D、D→A和D→C 的转移过程中,都存在较为显著的差异。另外,ASD 组中有5 名被试的转移概率矩阵显示出了对称性,这种情况TD 组中有6 名。从各个类别之间转移概率的平均值来看,M→N和N→M(M,N∈(A,B,C,D))存在的差异极小,并且在两组中都是如此。

表2 转移概率统计分析结果Tab.2 Transition probability statistical analysis results

本研究建立的零假设为时序的微状态序列具有独立性,即P(Xt+1)= P(Xt+1｜Xt)、P(Xt+1｜Xt,Xt-1)= P(Xt+1｜ Xt)、P(Xt+1｜ Xt,Xt-1,Xt-2)=P(Xt+1｜Xt,Xt-1)应用卡方检验,验证零假设与实际观测值的显著差异性,显著性水平设置在0.05,结果如表3所示。3 组结果均显示P＜0.05,即当前结果拒绝零假设,当前微状态的出现至少依赖于前两个微状态。

表3 马尔科夫模型计算结果Tab.3 Markov model calculation results

3 讨论

本研究在13 名ASD 儿童和15 名TD 儿童之间进行了微状态分析比较,确定了4 个微状态类别,用来解释这两组数据之间的差异。本研究使用IBM SPSS 26.0 版本的软件进行统计分析,显著性水平设置在P=0.05。结果清晰地表明,ASD 组和TD 组儿童在地形图拓扑结构和4 个时间参数之间存在差异。

在本课题中,结合了国内外对于脑电微状态的研究。在前人对于4 个微状态类别的研究得出的一致结论中[20-24],这4 个微状态类别分别对应于一个从血氧水平依赖(blood oxygen level dependent,BOLD)信号获得的大尺度静息态网络。对微状态进行源定位的研究发现,微状态A 与双侧颞上回和颞中回的负BOLD 激活有关,在语音加工过程中起着关键作用[10]；微状态B 与双侧纹外视觉皮层的负BOLD 激活相关,表征视觉静息网络[11,13]；微状态C与前扣带回的后部、双侧下额叶皮层与右侧前脑岛中的正BOLD 激活相关,这些区域与自我表征、情感交流和认知控制网络有关,反映默认网络模式的一部分,在执行认知网络的过程中会减少出现的频率和时间[25]；微状态D 与额叶和顶叶皮层的右侧区域和腹侧的负BOLD 激活相关,这些是与注意力网络相关的脑区[11,13],每个微状态都代表了一种静息态网络状态。对于两组中存在的不同微状态转移偏好,反映了ASD 患者异于正常人的信息加工过程和动态脑网络。

从本研究中对时间参数分析显示的结果来看,ASD 组的4 个微状态类别的发生频率都显著高于TD 组的相应状态,在ASD 组中各个类别具有更高的切换频率。除此之外,TD 组的4 个微状态类别的持续时间均高于ASD 组的相应状态,根据Koenig 等的研究,平均持续时间可反映神经元集群的稳定性[26]。

微状态A 的持续时间在TD 组显著大于ASD组,时间覆盖率在TD 组显著大于ASD 组,即微状态A 在TD 组中具有更高的存在率。这表明,ASD 儿童患者的大脑参与语言加工和处理过程的时间过少,而语言加工能力的受损是导致ASD 患者社交障碍的原因之一[27]。

微状态B 在ASD 组中的发生频率和时间覆盖率显著高于TD 组的相应状态,ASD 患者的一个典型特征是较强的知觉加工能力。有研究者曾提出“知觉功能增强模型”[28],认为ASD 患者的感觉加工系统更加深入地参与信息的加工过程。即使是在非任务状态下,ASD 患者与视觉加工相关的脑网络出现的频率和时间仍然高于正常儿童的相应状态。本研究的结果恰好与该理论契合。

微状态C 在ASD 组中的持续时间显著低于TD组的持续时间,反映了ASD 患者在自我表征能力和情感交流方面存在缺陷,这是由于前突神经与内侧前额叶皮层/前扣带回皮层之间的功能连接性降低[29],也是导致ASD 患者社交障碍的原因之一[27]。此外,微状态C 在ASD 组和正常组中,相比其他3种微状态类别,都表现出了较高的存在率。根据Michel 和Koenig 等的研究,微状态C 的潜在源与默认模式网络(default mode network,DMN)存在部分重叠[10],且DMN 网络在静息状态下出现较多。

微状态D 在ASD 组的发生频率显著高于TD组的发生频率,这和Jia 等的研究结论[23]一致,这可能是由于在ASD 患者中微状态D 所对应的电生理活动在右侧额中回显著比TD 组的活跃[23]。右侧额中回被认为是腹侧注意网络的重要组成部分,与注意力过程相关[24],另外,在其他的研究中发现,额区的异常是导致ASD 脑发育缺陷的主要原因[30],本研究关于微状态D 的结果也证实了ASD 儿童在额区的异常发育。

最后,通过计算0、1 和2 阶马尔科夫模型,探索了微状态时间序列上的独立性特征。结果显示,无论是在ASD 组还是TD 组,微状态序列在时序上都存在依赖性,当前时间点的微状态至少与前两个微状态存在关联性,表明了微状态序列上的信息存在共享性,当前时间点的微状态携带了其他微状态的信息。

4 结论

本研究利用微状态分析方法,探究了ASD 儿童和TD 儿童的脑机制差异。研究发现,ASD 与TD 儿童微状态的时间特征表现出显著性差异,ASD 儿童患者相比于TD 儿童,在语言加工、视觉网络、自我表征、情感交流和注意力过程等相关的大脑神经网络和信息加工方面发生了显著变化,并且在微状态时间序列中,微状态之间相互依赖、信息互享。EEG微状态分析可以揭示ASD 儿童的异常脑发育,而微状态的时间参数特征可以作为研究ASD 病理机制的生理指标。