基于复杂网络可视化的癫痫患者大脑状态研究
何艳1,2，于云莉3，杨帆1
1.贵州医科大学 生物与工程学院，贵州 贵阳 550004；2.西安交通大学 生命科学与技术学院，陕西 西安710049；3.贵州医科大学附属医院 神经内科，贵州 贵阳 550004
[摘 要]癫痫是一种伴随短暂认知损伤的大脑疾病，由于病理机制的异质性，导致其缺乏临床诊断和治疗评价的客观度量。多通道脑电记录是检测大脑状态的重要手段。本文采用复杂网络可视化的方法对癫痫患者大脑功能状态进行研究，利用相位锁定值方法从癫痫患者的静息态皮层脑电图（EEG）信号中提取出加权频率依赖的复杂网络，使用Pajek软件和最小生成树算法对癫痫患者大脑的复杂网络进行描述，可为量化大脑特征提供系统化、全局化的可视化思路，未来可以为进一步辅助诊断提供技术支持。
[关键词]复杂网络；脑电；可视化；癫痫；大脑网络
0 引言
癫痫是一种伴随短暂认知损伤的大脑疾病，严重威胁患者生命健康，是临床常见神经系统发作性疾病之一。由于癫痫病因复杂，早期的癫痫疾病分类主要依据发作类型可分为单纯部分性发作和全面性发作，脑电记录是其临床诊断主要的辅助手段。根据发作脑区的起源部位不同，可分为局灶性（如额叶、颞叶、枕叶、顶叶癫痫）和全身性癫痫。同一发作症状也可能对应不同的病因，准确的癫痫疾病分类诊断需要经验积累，而且大多是描述性分类，缺乏客观有效诊断指标[1]。
人脑连接组（Human Connectome）这一概念在2005年由Olaf Sporns与Patric Hagmann提出，旨在详细地刻画人脑神经结构网络，进而研究大脑的工作原理。2010年国家重点基础研究发展（973）计划中设立了“基于影像的脑网络研究及其临床应用”项目，通过多模态影像技术，搭建脑结构与脑功能的对应联系。由于大脑包含至少1010个神经元，其神经突触连接有1014个之多。脑网络数据量巨大，每个复杂网络可能包含巨量网络节点，节点间拓扑结构复杂。如何建立直观有效的可视化方法，显示脑网络特征信息，辅助脑网络观测与分析，对于研究大脑功能、揭示大脑工作机制具有重要意义。
大脑是一个高度自组织连接的复杂系统，神经元放电活动是多个脑区协同工作的基础。作为一个大规模结构和功能性整合网络，大脑包含大量弱相互作用部分。自1998年Watts和Strogatz提出小世界网络概念后，复杂网络被应用于多个领域。由于多通道脑电信号可以客观记录各个脑区信息流和连接特征，诸如核磁共振成像（Magnetic Resonance Imaging，MRI）等神经影像技术可以记录多个脑区的解剖和功能特性，基于复杂网络的脑科学研究可在图论框架里被有效研究。已有研究证实健康志愿者的大脑具备小世界特性，从而保证其快速处理外界刺激实现认知功能。基于复杂网络的脑功能研究内容大致包括静息态结构、瞬时信息处理能力、大脑活动动态演化机制等；通过脑功能神经影像数据（包括脑电）分别从连接特性（如结构连接、功能连接和有效连接）和节点特性及网络演化动力学机制等多角度揭示大脑活动规律，结合统计学分析，探索病理状态下大脑的异常指标[2]。
人类脑电的量化分析是客观、准确评价大脑状态的途径之一。大脑中神经振荡的自组织特性为大脑功能完成涉及的局部处理和全局整合提供精确信息。大脑是一个大规模结构和功能性整合网络，有大量弱相互作用部分，可以在一个数学图论框架中被有效研究。脑电功能性连接分析关注的是大脑多通道信号之间的同步或耦合情况，由于多通道脑电信号可以看做是大脑多个信号源叠加产生，从多通道数据序列中还可以提取各个脑区信息流和连接性特征。将大脑模拟为一个复杂网络，由节点集和连接边集合描述网络特性。对大脑连接性分析的核心问题是如何计算各节点之间的网络连接，即矩阵中各元素的数值如何得到。大脑系统功能连接性度量方法有许多，比如非线性相互依赖，相位同步，互信息，互相关，相干函数、互相关系数等技术可以揭示不同脑区的耦合强度，交叉互信息函数和合适替代数据可以衡量脑电信号之间的线性和非线性耦合作用，相位锁定因子可以更好地衡量窄波信号的相位同步，在脑电中也得到了极大的应用[3]。Stam等[4]利用同步似然法很好地提取了脑电中的功能连接性。相位动力学特征值分解方法的定向统计特性可以用来定义显著同步振子集群[5]。这些方法对于脑电信号连接性度量的结果大致相同，提供了有用量化信息。基于复杂网络的脑功能研究结果表明癫痫患者大脑解剖结构和网络特性与健康对照组存在显著差异[6-9]。基于图论的癫痫发作网络特征研究发现发作时脑区关联矩阵特征向量随时间变化且节点连接特性改变[10]。由于癫痫患者大脑结构和功能的网络变化（包括其连接性和拓扑结构），网络性能的破坏也与癫痫患者认知和行为损伤关联，使得当前癫痫研究从“局灶”向“网络”过渡[11]。
1 材料与方法
1.1 复杂网络提取
大脑多通道脑电记录中存在同一信号源诱发的容积效应，从而影响连接性作用判断。相位分析适合于窄波信号，比如脑电这种多频率成分的非平稳叠加信号，因此研究人员开始注重对脑电信号的相位分析。有研究证实，多通道信号的相位作用并不会受容积效应影响[3]。通过相位同步分析可发现大规模网络形成相关信息，但由于其对称性无方向度量，只适合功能连接性分析。相位锁定值（Phase Locking Value，PLV）和相位延迟指数（Phase Lag Index，PLI）是大脑功能连接性度量的一种相位方法。PLI估计的是瞬时相位延迟的对称性，若不对称则结果不为0，从而可以提取不可忽略的连接，代价是排除了有零相位延迟的相互作用。通过电极记录信号对相应脑区计算PLI均值，得到权重度或节点强度[4]，单个值反映节点在网络中的重要性，区域间的均值反映连接成本。加权PLI能较好地提取事件相关相位同步，但对于静息状态下的EEG信号则容易忽略零延迟相位的相互作用，从而不能准确提取频率依赖相位相关关系。经过对比PLI和PLV对相位同步的检测能力，本研究采用PLV作为大脑功能性连接的度量指标。PLV较好地克服了脑电产生过程中的容积效应，对频率成分有较好的敏感性。为了减少计算消耗，采用基于循环统计和错误发现率进行连接显著性分析。时间窗的选择对于分析结果很关键，之前时间窗的选择大多设定为100 ms左右，Dimitriadis等[3]提出了循环标准，即时间窗与对应频率有关。时间窗需要包含至少2.25倍频率循环周期才可以保证相位耦合度量鲁棒性。相位同步性的度量至少要通过二又四分之一个振荡循环来计算，若采样率为1 kHz，时间窗为300 ms，则最低同步频率周期需要比150 ms短，对应频率约6.6 Hz。通过滑动窗来进行计算，PLV是以Morlet小波与脑电信号Xi(n)某频率波段卷积所得到的瞬时相位为基础计算得来的[12-13]。瞬时相位φxi(n,f)的的计算公式为：
其中n代表电极位置，s代表小波函数尺度阶数，Xi(n)代表某段脑电信号时间序列。原始相位值范围为[-π,π]。经过一个相位展开去除不连续性，从而得到修正相位，随后对不同位点的信号xk(n)，xl(n)对瞬时相位差取平均得PLV指标:
N是所有时间点（样本）个数，s1/s2代表尺度限制，Δs代表尺度范围，f代表所在频率范围。若PLV为0，表示两个信号间不存在相位耦合。若PLV为1则表示两个相位完全同步。基于瑞利的PLV值检验对每个数值计算显著性[14]，显著性计算公式为：
为了提取显著性连接，引入错误发现率，当假阳性比例期望q≤0.01时设定对应的数值为显著性阈值。通过与每个通道随机打乱样本序列的替代数据的PLV对比，确定显著性PLV值；显著水平为原始数据和替代数据差值的Z分数；显著概率由假阳性率修正，从而对多个对比进行修正。
1.2 复杂网络可视化分析
复杂网络结构庞杂、节点众多、节点之间相互作用数量大且强度不一，导致其难以通过文字和表格形式一一呈现。如何能直观地还原复杂系统中各节点的相互作用关系和强度，如何提供观察网络结构、揭示有价值信息的方法，这些问题已经成为复杂科学研究中的重要部分。力导算法、可视化框架和分析工具及系统[15]的相继涌现使复杂网络能够应用于多个领域。目前有多个软件（如UCINET、Pajek、Cytoscape、Gephi等）已经实现了包含大量节点网络的可视化分析[16]，能提供基本的复杂网络特征度量。
本文使用Pajek软件，在癫痫患者大脑复杂网络邻接矩阵提取的基础上，实现其大脑网络的可视化呈现。
2 复杂网络可视化结果
本研究的脑电数据采用国际通用的10-20电极安置系统，电极名称分别为Fp1、Fp2、F3、F4、C3、C4、P3、P4、O1、O2、F7、F8、T3、T4、T5、T6、Fz、Cz、Pz，采样率为250 Hz，采样时间为长程信号记录，从中提取出10 min的脑电数据进行实验。数据来源于1例贵州医科大学附属医院门诊癫痫患者（男，30岁），排除因外伤或短期精神刺激而诱发的癫痫，其未接受过药物治疗或干预，核磁共振成像未观察到明显的脑区损伤或肿瘤，临床诊断为原发性癫痫。原始脑电信号分被为delta（0.5～4 Hz）、theta（4～8 Hz）、alpha（8～13 Hz）、beta（13～30 Hz）、gamma（30～45 Hz）等5个频率波段。多通道脑电记录各脑区对应的记录电极定义为网络节点，各电极所记录的脑电信号的相互作用定义为连边，在本研究中为相位锁定值。由于多通道脑电信号之间存在容积效应，信号源之间会相互干扰，引入相位成分计算各个脑区脑电信号之间的相位锁定值，该方法目前被广泛应用于大脑网络的功能连接分析，能通过循环统计和错误发现率确定显著性连接。在PLV的作用下，邻接矩阵单元数值等价于PLV数值，即各节点之间作用强度由PLV数值表征。各个频率波段组成的复杂网络可视化可以通过Pajek软件实现。gamma波段组成的功能连接网络，通过改变网络布局算法，可以调整网络可视化结果，见图1。
基于Pajek软件的可视化能实现网络直观表达并能区分出节点连接强度最高的节点属性，但由于不同频率波段之间大脑节点作用强度不同，当节点连接过多时，无法直观分辨出关键节点，因此引入最小生成树算法对复杂网络进行细致描述。最小生成树（Minimum Spanning Tree，MST）是指在一个具有几个顶点的连通图G中，若存在一个子图G'包含G中所有顶点和部分边，且不形成回路，则称G'为G的生成树，代价（权重）之和最小者则称为MST。依据MST确认通道网络，可以确定整个全局通道关系，提供皮层脑电数据的唯一网络结构。由于MST依赖于权重网络，需要构造距离矩阵。本研究采用Dijkstra算法对不同癫痫患者大脑的距离矩阵提取MST，该算法是典型的单源最短路径算法，用于计算一个节点到其他所有节点的最短路径，主要特点是以起始点为中心向外扩展，直到终点为止，目前已经用于癫痫疾病分析[17]。脑电信号不同频率成分（如delta和gamma）在不同节点之间的PLV作为复杂网络的关联矩阵，利用matlab中的函数“graphminspantree”实现对应复杂网络最小生成树的可视化。delta波段（图2）和gamma波段（图3）所在复杂网络的最小生成树，可以看出不同频率波段中关键节点的部位和连接强度均有所差异。
图1 基于Pajek软件的 gamma波段复杂网络可视化示意图
图2 基于MST的delta波段复杂网络可视化示意图
图3 基于MST的gamma波段复杂网络可视化示意图
3 讨论
脑电数据是大脑功能状态的客观反映，多通道脑电记录是临床常用的辅助诊断手段。由于脑电记录时间长，数据量大，无法被直观理解。如何寻找合适的数据表达方式并使得所有数据发挥作用、体现价值，是脑科学研究的重要组成部分。数据可视化的基本过程包括数据理解、数据表达及可视化。通过本文提出的复杂网络可视化方法可将脑电反映的信息直观表达。PLV方法能有效刻画脑电记录节点之间的相互作用强度。基本的节点和连接形式使得从整体上分析大脑功能状态成为可能。基于最小生成树的可视化表达则能有效提取出有效最短路径及唯一拓扑结果，这就使大量无法直观分析的临床数据得以表达。
认知神经科学研究表明，脑电信号中的不同频率成分对应着不同的功能状态。如清醒状态下脑电信号以alpha波段为主，在神经系统疾病中慢波振荡（如delta波段）占据更多能量成分，gamma波段则更多地参与学习记忆等高级认知功能。delta频率成分之间的耦合强度高于gamma波段频率成分之间的耦合强度，与癫痫疾病中的慢波振荡增强导致认知损伤加重等研究结论一致，见图2和图3。此外在不同频率波段组成的复杂脑网络中，delta波段中C3和F3表现出枢纽特性，而gamma波段中则是F3、F4、P3表现出枢纽特性，不同频率成分的枢纽节点既有不同又有重合，可能预示着记录电极F3对应的脑区参与该患者的癫痫病理，但具体作用机制还有待进一步研究。
4 结论
早期研究人员通过脑电地形图能观察到不同类型癫痫患者的EEG中，在整个大脑区域中各频率波段能量分布有所不同，但难以从脑电地形图中提出大脑网络性能对应指标或度量参数。图论是一种描述并度量网络的数学语言，在引入随机，小世界和无标度网络模型后，图论快速发展并大量运用于大脑神经科学研究。连接网络分析能够反映出大脑功能连接结构。基于PLV的大脑功能连接性分析能较好地反映大脑各脑区、各窄波频率波段相位同步状态。Pajek软件实现了多通道脑电信号组成的大脑复杂网络的有效可视化，通过改变网络布局算法能有效观测节点在网络中的连接状态和权重大小。MST有效提取出复杂网络最小生成树对应的唯一拓扑网络结构，结果发现不同频率波段组成的复杂网络其唯一拓扑网络结构与关键节点均有所差异。如何从庞杂的医疗数据中提取关键信息，并融合信息使其可视化，一直是神经科学研究中的热点和难点问题。本文设计的癫痫患者大脑网络可视化方法为深入研究癫痫作用机制提供了新途径。通过可视化手段将长程记录的多通道脑电信号表达为节点和连边的组件，并提供直观量化结果，为临床辅助诊断提供了技术支持。
[参考文献]
[1]Shorvon SD,Perucca E,Fish DR,et al[著].肖波,刘献增,龙小艳,等[译].癫痫治疗学[M].北京:人民卫生出版社,2010.
[2]张方风,郑志刚.复杂脑网络研究现状及挑战[J].上海理工大学学报,2012,34(2):138-153.
[3]Dimitriadis SI,Laskaris NA,Tsirka V,et al.Tracking brain dynamics via time-dependent network analysis[J].J Neurosci Methods,2010,193(1):145-155.
[4]Stam CJ.Nonlinear dynamical analysis of EEG and MEG:review of an emerging feld[J].Clin Neurophysiol,2005,116(10):2266-2301.
[5]Quian Quiroga R,Kraskov A,Kreuz T,et al.Performance of different synchronization measures in real data:A case study on electroencephalographic signals[J].Phys Rev E,2002,65(041903):1-14.
[6]van Diessen E,Diederen SJ,Braun KP,et al.Functional and structural brain networks in epilepsy:what have we learned?[J].Epilepsia,2013,54(11):1855-1865.
[7]Prinz AA.Understanding epilepsy through network modeling[J].Proc Natl Acad Sci USA,2008,105(16):5953-5954.
[8]Onias H,Viol A,Palhano-Fontes F,et al.Brain complex network analysis by means of resting state fMRI and graph analysis:will it be helpful in clinical epilepsy?[J].Epilepsy Behav,2014,38:71-80.
[9]Dyhrfjeld-Johnsen J,Santhakumar V,Morgan RJ,et al. Topological determinants of epileptogenesis in large-scale structural and functional models of the dentate gyrus derived from experimental data[J].J Neurophysiol,2007,97(2):1566-1587.
[10]Burns SP,Santaniello S,Yaffe RB,et al.Network dynamics of the brain and infuence of the epileptic seizure onset zone[J].Proc Natl Acad Sci USA,2014,111(49):E5321-E5330.
[11]Jirsa VK,Stacey WC,Quilichini PP,et al.On the nature of seizure dynamics[J].Brain,2014,137(Pt 8):2210-2230.
[12]闫佳运,宾光宇,梁栗炎,等.基于单路脑电的睡眠分期判别方法研究[J].中国医疗设备,2015,30(12):34-37.
[13]Freeman WJ,Rogers LJ.Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs[J].J Neurophysiol,2002,87(2):937-945.
[14]Achard S,Salvador R,Whitcher B,et al.A resilient,lowfrequency,small-world human brain functional network with highly connected association cortical hubs[J].J Neurosci,2006,26(1):63-72.
[15]吴彬彬,王京,宋海涛.基于Citespace的复杂网络可视化研究图谱[J].计算机系统应用,2014,23(2):7-16.
[16]李建勋,李维乾,郭莲丽,等.企业研发团队复杂网络的可视化与描述[J].计算机系统应用,2013,22(11):209-212.
[17]Lee U,Kim S,Jung KY.Classifcation of epilepsy types through global network analysis of scalp electroencephalograms[J].Phys Rev E Stat Nonlin Soft Matter Phys,2006,73(4 Pt 1):041920.
Investigation of Epilepsy Patients’ Brains Based on Complex Network Visualization
Abstract：Epilepsy is a kind of brain disease accompanied by transient cognitive impairment. The objective metrics for clinical diagnosis and treatment evaluation of the disease is defcient because of the heterogeneity of the pathological mechanism. Multi-channel EEG (Electroencephalogram) recording is an important technique for monitoring of brain functions. In this paper,complex network visualization is adopted to investigate the overall state of brain functions in patients with epilepsy. Transient phase locking values are brought in to extract weighted frequency-dependent complex networks from resting state cortical EEG signals of epilepsy patients,then Pajek software and minimum spanning tree algorithm are introduced for depiction of these complex brain networks,which provides a systematic and global visualization framework for quantitation of brain characteristics of epilepsy patients,and lays a solid technical foundation for development of intuitive as well as effective auxiliary diagnostic tools for the future.
Key words：complex network;electroencephalogram;visualization;epilepsy;brain network
HE Yan1,2,YU Yun-li3,YANG Fan1
1. School of Biology &Engineering,Guizhou Medical University,Guiyang Guizhou 550004,China;2. School of Life Science and Technology,Xi’an Jiaotong Univeristy,Xi’an Shaanxi 710049,China;3. Department of Neurology,The Affliated Hospital of Guizhou Medical University,Guiyang Guizhou 550004,China
[中图分类号]R742；R318
[文献标志码]A
doi：10.3969/j.issn.1674-1633.2016.09.009
[文章编号]1674-1633(2016)09-0039-04
收稿日期：2016-04-12
修回日期：2016-05-05
基金项目：国家自然科学基金项目（81460206）；贵州医科大学博士启动基金项目（院博合J2014[003]）。
通讯作者：何艳，博士，副教授，主要研究方向为复杂网络与复杂系统，神经信息工程。