Nature Communications：致痫脑源时空分布的无创电磁源成像

    脑网络是随时间动态变化的时空现象。功能成像方法致力于无创地评估这些潜在的过程。来自卡耐基梅隆大学大学的Sohrabpour等人提出了一种新的源成像方法：使用高密度的脑电图记录来绘制脑网络。

这种方法客观地解决了传统源成像技术的长期局限性，即难以客观地估计潜在源的空间范围以及时间演化。研究者们通过直接比较36例局灶性癫痫患者的颅内脑电图（iEEG）结果和手术切除结果来验证该方法。本研究共分析了1027次放电和86次发作，结果展示了此方法在无创性电生理测量中成像脑网络的位置和空间范围的能力，特别是对于发作期和发作期内的脑网络。本方法为无创性研究大规模动态脑网络提供新方法。本研究发表在Nature Communications杂志。

方法

   参与者信息

   本研究共纳入36例患者。这些患者由医生根据ILAE(international league against epilepsy,国际抗癫痫联盟)系统评分，平均监测18个月(随访时间)。在19个月的随访期内，21名患者被评为ILAE 1（完全无癫痫发作），3名患者被评为ILAE 2（无癫痫发作，只有先兆），12名患者在15个月的随访期内被评为ILAE 3-6（非癫痫发作）。 

       算法概要

    本研究使用ESI(electrophysiological source imaging, 电生理学源成像)优化算法，使FAST-IRES(fast spatiotemporal iteratively reweighted edge sparsity,快速时空迭代加权边缘稀疏性)算法适用于动态数据分析。该算法的基本思想是，EEG/MEG记录的源不是极度聚焦的，而是聚焦扩展的，因为EEG/MEG信号是许多同步神经元集合的突触后电位的叠加。此外，随着时间的推移，这些放射源相互叠加，在头皮上可以检测到这些放射源的活动过程。基于这些前提，我们假设，如果通过成分分析从头皮记录中描绘出潜在脑活动的TBF(time basis function,时间基函数)，则估计TBF的每一行对应于一个空间扩展的焦点源。优化问题的类型基于这些假设的求解形式如下：

    其中，ϕ(t)是感兴趣间隔内头皮电位（或磁场）测量值（E×t矩阵，其中E是测量次数，t是给定间隔内的时间点数），K是lead field矩阵（E×N矩阵，N是源数），j是大脑区域的未知电流密度（anN×Nc矩阵，其中Nc是TBFs的数目），A是时间过程激活矩阵（Nc×T矩阵）或TBF，由A=[a1（T），a2（T），…]，β²本质上是噪声功率，由差异定理确定，∑是根据基线活动度Wd确定的噪声协方差矩阵，iL−1和WiL−1是与每个j_i相关的权重，用IRES迭代中确定的相同规则进行更新，V是离散梯度算子，α是正则化项中两项之间的超参数平衡，将使用L曲线方法进行调整。

       FAST-IRES原理、限制和参数：

   FAST-IRES的主要原则可以归纳为四个部分：

1. 利用L1范数正则化项最小化源的边界。空间梯度的L1范数，即相邻两个源的振幅之差--边缘，被最小化。这使得解决方案能够在有限的位置或边缘进行突然的更改或跳跃。这是由于L1范数的性质，在有限数量的矢量元上允许振幅突变，而L2范数则强制执行平滑。

2. 最小化L1范数会强制源具有零背景(zero background)，从而抑制恒定背景。

3. 定义参数的平衡项。L曲线法可以客观地实现这一目标。

       4. 进行迭代加权。

其中正则函数的两个项被加权。为了评估L曲线对FAST-IRES估计的影响，必须承认参数α影响源估计，因为它在正则化的两个项之间提供了平衡。但是，可以使用L曲线方法来调参。L曲线法基于Pareto最优化思想，即正则化的两个项都应尽可能小，通常选择与L曲线拐点相对应的α值。这种方法本身并不是主观的。但是，由于噪声或其他不可预测的因素，曲线可能会发生一些变化。这是一种不能排除的可能性。我们在补充图13中提供了一个例子，表明即使在L曲线弯曲处附近将α改变10倍也不会对解产生太大影响。我们认为，提出的算法对α的适度变化是鲁棒的，并且使用L曲线方法可以系统地减少需要考虑的α的范围，因此它似乎不会严重影响源估计。

补充图13  L曲线和选取不同α对估计的影响示例

 

   数据记录与处理

   患者在术前评估阶段接受76导EEG记录。通道排布为10-10系统，使用CPz电极进行重参考。原始数据采样率为：500 Hz，高通滤波1 Hz，以去除数据中的虚假慢活动和可能的直流漂移。对于发作性信号，进行20–30 Hz低通滤波器，而对发作间期信号进行50 Hz低通滤波。

    对于每一位患者，使用受试者的个人磁共振成像（MRI）建立个性化的头模型，形成由K表示的lead-field矩阵。采用三层边界元方法（BEM，boundary element method）模型求正向解，得到lead-field矩阵。该模型由三层组成，分别代表大脑、颅骨和头皮，其相应的电导率分别为0.33、0.0165和0.33 s/m。采用皮层电流密度（CCD，cortical current density）模型开发FAST-IRES方法。

       ICA及成分选择

   为了去除信号的噪声，提取TBF，我们对EEG信号进行了独立成分分析(ICA)。棘波数据和发作数据的处理方式不同。

    对于发作数据，在每次发作中，选择发作开始前10s至到发作结束的数据（癫痫发作时间由训练有素的癫痫学家标记），对每次发作进行分析（在两名癫痫发作持续时间超过2-3分钟的患者中，只选择前2分钟）。使用EEGLAB工具箱（版本14.1.1b）中的logistic infomax ICA算法进行ICA。ICA将发作数据分解为空间固定源和时间独立源。

   在ICA运算结束后，人工识别并去除眨眼、眼球运动等常见伪迹成分。此外，为了识别与癫痫发作相关的成分，我们选择了与癫痫发作时间演变高度相关的独成分（ICs），即发作前无症状或活动较少，发作后活动性较高。使用这些成分的时间进程作为TBF，随后将其输入FAST-IRES算法。

    对于棘波数据，我们在每个患者中选择发作棘波周围2s的是数据并进行拼接，对拼接好的信号进行ICA分析。在存在多个棘波类型的患者中，我们对每种棘波类型独立执行该过程。数据分析与之前类似，主要区别有两个：

    1.对发作期间数据进行1-50Hz的带通滤波。

    2.剔除明显的伪迹成分后，对棘波数据的时间进程做平均。

    然后，我们检查了这些平均的时间进程，当他们在棘波点（与基线相比）出现明显的振幅增加时，将其作为与棘波相关的成分。将平均的时间进程输入FAST-IRES，以解决源成像问题。

补充图14  ICA去噪效果示例

 

   频谱分析

   在对发作数据进行去噪并形成TBF后，将FAST-IRES算法应用于癫痫发作后3-5s的数据中。以这种方式，获得了在癫痫发作初始阶段的时空分布。与我们对癫痫发作传播和节律性的预期一致，这种分布不断波动并迅速传播，迅速涉及不在SOZ(seizure onset zone,癫痫发病区)内甚至空间上与之相邻的区域。对于颞叶癫痫患者，活动可能在几毫秒内传播到对侧叶，这使得从源空间分布确定SOZ的任务变得困难。事实上，在我们研究的中，有13名患者在发作后（最早，在发作后几秒钟后）出现继发性全身性癫痫发作或向对侧扩散。在IC选择过程中，观察到对侧成分，并且这些成分通常包括在TBF（时间基函数）中，因为我们预期并观察到癫痫向对侧的传播（补充图15）。

补充图15 颞叶癫痫患者的发作传播示例

 

    根据发作性源成像的文献，在主频率处进行滤波以获得更具空间相关性的信号。为此，在所有通道上，计算发作前5s和发作后5s数据的平均功率进行对比。分别将发作前和发作后的数据，按照1s时长进行分段，先计算每个时段的平均功率，再进行总平均，最终得到两个平均功率：发作前和发作后。比较两个频谱并选择存在显著差异的频段进行后续分析。选择标准是：发作后功率必须大于发作前，并且只选择频谱中最强大的一个，并按照所选频率对数据进行带通滤波。随后计算并显示癫痫发作开始时1秒间隔内该频段的平均功率，如图5a。

    定义性能指标

   使用重叠指标(overlap metric)，比较估计源的形状和相对位置与真实情况的差异情况。计算估计源和切除区域之间的重叠量，并除以切除面积或估计源的面积，得出标准化重叠比（NOR，normalized overlap ratio）。这也说明了两个分布之间的匹配程度，两个比率的理想值都是1。如果得到一个高估或低估的解，两个度量值中的一个将接近1，而另一个将显著降低。在一些统计和计算机科学文献中，这两个NORs也被称为精确度和召回率。

    本文中使用的另一个性能指标是定位误差（LE,localization error），它被定义为每个SOZ电极之间的平均最小距离（在接受iEEG植入的患者中计算）和估计的EZ(epileptogenic zone,病灶)之间的平均最小距离。 

       比较FAST-IRES时间进程与颅内痕迹。

   为了评估我们估计的时间进程与真实时间进程（即潜在源活动的时间进程）的接近程度，对一名患者进行了额外的分析。在这项分析中，FAST-IRES估计的活动时程平均分布在外侧颞前区和海马旁深部皮质，在那里放置了颅内电极，最终确定为SOZ。这些平均时间进程与这两个区域的平均颅内轨迹相关，并且两个区域的平均相关性达到0.6。对这名患者分析的所有三次发作重复这一过程，对于每一次发作，提取的时间过程与iEEG电极记录的三次随机选择的癫痫发作进行比较。

    统计检验和可重复性

   为了确定棘波和发作信号评估的性能指标之间的统计差异，我们采用了置换检验。考虑到这里所研究的性能指标不太可能来自正态分布，而且样本量有限，我们更倾向于使用置换测试。比较两组的癫痫发作分析结果和棘波分析结果。使用MATLAB进行的置换检验次数为104次，显著性水平为p<0.05。本研究还进行了其他双侧统计检验，如Welsh的双样本t检验和Wilcoxon秩和检验，并且能够得出类似的结果。

结果

   扩展源的蒙特卡罗模拟

   该方法的基本思想是同时估计潜在源的时空分布，而不是在每个时间点估计其潜在的空间分布。这使得开发更高效和准确的算法和更精确的大脑基本过程建模成为可能。假设潜在源在空间上是集中的，也就是说，局限于大脑的局部区域，不包括大脑皮层的很大区域。我们把具有这种特性的源称为聚焦扩展源。这种特性，很可能不是癫痫源特有的，在其他大规模的现象中也被提出过，在头皮水平测量时，必须同步激活延伸的皮层区域，以产生可检测的信号，如EEG和MEG。请注意，脑活动的精细微尺度组织，在EEG/MEG或任何表面记录中都无法感知，我们在这项工作中并不声称能够恢复这些来源。

图1 时空源成像方法的概念

 

    聚焦扩展源的想法是模拟大规模的脑信号和脑组织，这些信号和组织通常记录在EEG和MEG等头皮测量中。我们工作的一个关键特征是从数学上捕捉大规模脑组织的这一基本特性，并开发求解这些数学模型的计算工具。每一个聚焦扩展的空间源都有其独特的活动时间过程，对应于其振幅随时间的变化。由于所有源的活动都是线性叠加形成脑电记录的头皮电位，因此利用盲源分离（BSS，blind source separation）技术，可以从头皮测量值中描绘出源活动。此外，空间约束，如聚焦和活动与背景噪声之间存在明显的边界，可以通过迭代的方法实现稀疏性，最终促进聚焦扩展源。

图2 扩展估计的蒙特卡洛模拟结果

在四种不同信噪比条件下，我们的FAST-IRES算法估计范围的结果与模拟源范围进行了对比。四个图中的红线为参考线。

 

    为了深入研究我们的理论在可靠地确定底层源的范围方面的性能，我们进行了蒙特卡罗模拟（图2）。在这个模拟中，我们选择了大脑皮层上的随机位置，在这个位置模拟了一个扩展源，并计算了这些源产生的EEG信号。图2展示了我们的模拟结果，它是一个典型的mesio时间源，它在所有的信噪比条件下都是局部的。此源的活动估计时间过程与模拟的时间过程非常匹配。这些结果的Pearson相关值显著较高，为0.88。

   临床数据分析结果概述     

         对于所有形式的机器学习或结构化分解，识别相关的特征向量对于定量评估和结果的可解释性至关重要。在我们的例子中，为了对癫痫活动进行成像，我们首先需要从EEG记录中提取已知的癫痫特征，即发作间期的棘波和癫痫发作。我们研究了36名患有局灶性癫痫的患者，他们接受了手术治疗难以控制的癫痫发作。高密度EEG记录对于在ESI成像过程中获得准确的解决方案非常重要。我们开发了一种新的高密度（76通道）EEG技术，它允许我们在接受术前监测的病人中记录连续几天的脑电图，以捕捉发作间期事件和癫痫发作。我们分析这些新的数据来估计EZ，并定量比较我们的无创性ESI结果与有创iEEG结果以及手术切除结果。癫痫发作是由经验丰富的癫痫学家确定和标记的，同时，研究小组为每个患者提取了棘波。

图3 研究设计概览

图中描述了分析流程：

上：这项研究的主要两个方面展示了如何从EEG记录中提取棘波和发作数据，去噪，并确定它们的TBF（时间基函数），并将其输入FAST-IRES中。

中：FAST-IRES源成像方法考虑了源的空间范围和焦点。算法的输出是潜在源的时空分布，从中提取EZ并与临床结果进行比较。

下：最后，通过比较癫痫发作间期活动和FAST-IRES的发作活动来评估癫痫特征的表现，并将估计的EZ与临床结果进行比较。

 

    我们分析了癫痫发作和发作间期棘波，对比了每种特征的评估结果与临床发现的质量和一致性（图3）。独立成分分析用于抑制噪声，并确定潜在源的TBF。如前文介绍，选择TBF的标准对于棘波和癫痫发作数据是不同的。一旦计算出TBF，对棘波和癫痫进行源成像处理，并将每个特征的估计值与手术切除的体积和/或癫痫发作区（SOZ）的颅内电极进行比较（图3下）。 

       棘波成像

      一般认为，棘波来自刺激区，这与SOZ不同。虽然刺激区和SOZ的定义很简单，也很容易评估，但事实证明EZ(发病区）是难以捉摸的。临床上，EZ被定义导致癫痫发作的最小组织量。这个定义很抽象，很难评估。最接近EZ的可测量特征很可能是术后无癫痫发作的患者的手术切除量。为了便于比较，手术切除体积可以被认为是伪EZ或真EZ。虽然先前的文献明确地证明棘波成像结果与临床确定的EZs一致，但是棘波成像结果与来自发作记录的成像结果之间的关系仍然不确定。我们在本研究中调查了这种关系，以确定这些特征中的哪一个可以提供关于EZ的更多相关信息。

图4 棘波成像概览和结果

a.FAST-IRES的输出是一个时空分布，其中空间分布对应于活动的时间过程。

b.同一患者临床表现的棘波成像结果。

c.所有患者的棘波成像结果的定量结果（顶），并根据手术结果分开（底）。

 

    在每个病人术前的EEG记录中，选择棘波周围2s的对称窗口进行分析。然后对这些波形进行平均，并进行源成像（图4a）。将平均棘波及其对应的TBF输入FAST-IRES后，输出的是在2s间隔内变化的脑电活动分布。在棘波时间周围的40 ms窗口内对结果进行平均，以获得单个源分布，并与每个患者的临床结果进行比较（图4a）。图4b给出了一些棘波成像结果的例子，与每个病例的相应临床结果重叠。平均而言，棘波成像结果的精确度和召回率为0.6（图4c）。这大致相当于EZ与切除组织重叠60%的棘波成像估计值。这表明，可以通过棘波成像结果来估计EZ，特别是考虑到手术切除的面积通常足以确保所有致痫组织都被切除（因此他们高估了真实的EZ）。对于棘波成像结果，估计的EZ到SOZ电极的距离的定位误差平均约为18mm（而在癫痫成像中为6mm，图5c）。这符合这样一个事实，即癫痫发作时的棘波和癫痫发作并非来自大脑中完全相同的区域，因此传递了有关潜在癫痫回路的不同信息。

图5 发作成像概览和结果

a.FAST-IRES的输出是一个时空分布，其中空间分布对应于活动的时间过程。

b.同一患者临床表现的发作成像结果。

c.所有患者的发作成像结果的定量结果（顶），并根据手术结果分开（底）。

 

   发作成像

   虽然发作间期棘波已经过多年的分析，但癫痫发作（发作记录）通常没有分析。这在一定程度上是因为发作记录非常嘈杂，并且常常包含大量不需要的伪影，例如眼球运动、眨眼、肌肉和运动伪影。因此，对原始发作记录进行成像非常困难。目视检查原始EEG记录不能提供关于潜在癫痫发生组织的起源和范围的准确信息。这些障碍在推动发作数据定量研究的主要动机，这种定量方法可以提供大量关于癫痫发作网络的起源和程度的临床相关信息。       

       我们对癫痫发作的前几秒（3-5秒）进行成像，并估计在这个时间间隔内大脑中的电活动分布，然后在主要发作节律下对我们的解决方案进行带通滤波。为了找出每名患者的主要发作节律，我们将发作后的前5s数据的频谱与发作前5s的进行比较，并选择两个频谱差异最大的频带（图5a）。

    我们估计的EZ的精确度和召回率都很高，约为0.75，这表明我们推断的来源既没有低估也没有高估真正的潜在EZ（图5c）。这是因为精度和召回率是标准化指标，描述了源和真实情况的重叠，即手术切除的体积，相对于估计的源大小或真实性。因此，如果一个来源低估或高估了手术量的范围，那么只有一个值是高的。只有当估计的放射源和手术体积大致相同且重叠良好时，才会这两个值都高。

此外，在接受iEEG记录的患者中，定位误差约为5mm，这已经接近iEEG电极分辨率极限（定义为iEEG电极之间最小距离的一半，在我们的研究中通常为10 mm）。当绘制无癫痫发作和非癫痫发作患者组的结果时，没有观察到统计上的显著差异（图5c）。

    在接受有创iEEG记录的患者中，有87%的患者植入深部电极，以确定深层结构是否参与癫痫的发生，从而导致深部结构，如海马体被切除。然而，小的定位误差表明，使用我们提出的方法，即使是深源也可以被无创地描绘出来。此外，我们还计算了前颞叶外侧部分、海马体和海马旁皮质附近深部结构内活动时间进程的相关性，并将其与该患者的颅内记录进行了对比。

    我们的结果显示，在这两个脑区，估计的时间进程和iEEG记录之间的相关值为0.61。这种对大脑深层结构的检查是我们的无创性ESI方法的主要优点之一，它有助于研究癫痫网络。它也可以用来研究一般的大脑网络。 

       连通性成像

   FAST-IRES提供了潜在脑源的时空估计。这意味着，除了潜在的大脑活动的位置和空间范围外，这些活动的时间进程也可以估计出来。可根据FAST-IRES结果进行后续网络或连接分析。我们根据我们的FAST-IRES结果概述并执行了连通性分析。

图6 棘波成像和发作成像结果对比

    在定位误差、几何平均值和谐波平均值方面，癫痫发作成像的结果在统计学上明显优于棘波成像（图6a、b）。在进一步的调查中，我们发现这种差异（棘波和癫痫成像结果之间）只存在于术后没有癫痫发作的患者（图6a，b）。这很可能是因为有些患者有多种类型的尖棘波，其中一种类型可能是对侧（图6c），而癫痫发作始终发生在与临床结果相同的一侧。

    这表明，当出现不一致的棘波类型时，仅仅依靠棘波成像结果可能会产生误导，必须考虑癫痫发作的电生理评估。我们进一步比较一致的棘波与癫痫发作分析之间的差异。结果表明，在持续性癫痫发作和癫痫发作之间没有明显差异。此外，我们还观察了连通性和发作性成像结果之间的差异，发现两组之间没有统计上的显著差异。这表明尽管连通性成像方法提供了一个更具分析性的框架，但它并没有显著改善发作期成像结果。 

总结

   本研究展示了无创源成像方法的优点，以及它们可以用于脑网络成像的程度。换言之，此方法可以定位网络节点，确定这些节点的空间范围，估计节点活动的时间变化，计算节点间的连通性和动态性本研究展示了无创源成像方法的优点，以及它们可以用于脑网络成像的程度。换言之，此方法可以定位网络节点，确定这些节点的空间范围，估计节点活动的时间变化，计算节点间的连通性和动态性。此外，本研究证明，癫痫网络可以成功地成像，以确定个别患者的EZ（发病区），其结果与侵入性临床发现一致。特别是，当患者的脑电图记录中有多个空间不一致的棘波类型时，发作期成像优于棘波成像。 

原文：Noninvasive electromagnetic source imaging of spatiotemporally distributed epileptogenic brain sources
如需原文及补充材料请加思影科技微信：siyingyxf 或者18983979082（杨晓飞）获取,如对思影课程及服务感兴趣也可加此微信号咨询。觉得对您的研究有帮助，请给个转发，以及右下角点击一下在看，是对思影科技莫大的支持。

微信扫码或者长按选择识别关注思影

非常感谢转发支持与推荐