CN112867441A

CN112867441A - 用于监测神经信号的***及方法

Info

Publication number: CN112867441A
Application number: CN201980047997.7A
Authority: CN
Inventors: P·L·珀登; H·苏拉; A·M·贝克; E·P·斯蒂芬
Original assignee: General Hospital Corp
Current assignee: General Hospital Corp
Priority date: 2018-07-16
Filing date: 2019-07-16
Publication date: 2021-05-28
Also published as: EP3823527A4; JP2021531096A; WO2020018595A1; US20220142554A1; EP3823527A1

Abstract

本发明的一个方面提供一种用于估计一病人大脑状态的至少一个振荡分量的方法。所述方法包括：接收脑电图(EEG)数据；拟合一状态空间模型至所述脑电图(EEG)数据，所述状态空间模型包括至少一个振荡分量；根据所述状态空间模型，估计所述脑电图(EEG)数据的跨频耦合的至少一个值；根据跨频耦合的所述至少一个值，产生一报告；以及将所述报告输出到一显示器及一存储器的至少其中之一。

Description

用于监测神经信号的***及方法

相关申请的交叉引用

本申请基于提交于2018年7月16日、标题为“用于监测神经振荡的系及方法”的62/698,435号美国专利申请、对该美国专利申请要求优先权，并通过引用该美国专利申请将其整体并入本申请。

关于联邦资助研究的声明

本发明是在由美国“国家卫生研究院”(National Institutes of Health)授予的R01AG056015-01号、R01AG054081-01A1号及GM118269-01A1号项目资助的政府支持之下完成的。美国政府对本发明有一定权利。

背景技术

振荡在神经信号中无处不在，并且反映神经元群在广泛的时间和空间尺度范围内的协调活动。神经振荡在调节大脑功能的许多方面扮演重要的角色，包括注意力、认知、感觉介入处理以及意识。因此，神经认知障碍、精神障碍以及意识改变状态与脑震荡的改变或中断有关。神经振荡可以跨越单个神经元膜电位、神经回路及非侵入性头皮电磁场，大规模地记录在动物及人体模型中。因此，神经振荡可提供有价值的机理支架，将不同规模及模型的观测结果联系起来。

遗憾的是，与许多其他生理信号相比，神经振荡相当不易察觉，而且可能被无数的影响完全或部分地所遮蔽。因此，在试图获取神经信号及/或分析神经信号之前，通常先进行稳健信号处理及/或信号鉴别。

例如，脑振荡通常用频域方法来分析，比如用非参数谱分析或基于线性带通滤波的时域分析方法。典型的分析可能试图估计特定频率范围内的振荡中的功率，例如alpha频段振荡，其频率为8～12Hz之间。通常，8～12Hz之间的功率是使用功率频谱，在频域中计算，或通过估计带通滤波信号的功率或变化，在时域中计算。

这种性质的处理基于各种各样的假设，其中一些假设可能是错误的，而其他的假设可能是基于有意识的权衡而作出。例如，上述频率选择及滤波天然地区别对待某“信号”，以消除噪声。然而，与所述“噪声”一起丢失的、该“信号”的部份可能比想象的更有价值。此外，虽然许多信号处理技术(比如带通频率滤波)易于实施，但可能并非捕获该“信号”的理想信号处理技术。

因此，有必要拥有能更好地促进神经信号的采集及分析的***及方法。

发明内容

本发明提供多种用于处理与振荡状态有关的神经信号的***及方法，比如用于病人监护以及脑监测意外的应用中的信号解调。在一个非限制性方面，本发明提供用于捕获作为一状态空间模型内的线性振荡表达式组合的神经信号的***及方法。这些模型的参数可以使用一期望最大化算法或其他方法来选择适当模型，该模型用于识别数据中存在的振荡，并从而提取所需的神经信号。通过这些***及方法所提供的适当提取，本发明能够区分或将感兴趣的振荡与会混淆后续分析的潜在次级生理信号分离。因此，本发明的***及方法能够估计对监测应用有用的振荡的性质，比如相位与振幅互相关性、振荡相位、振荡振幅、振荡功率、相对振荡功率、振荡带宽、振荡谐振、振荡品质因数等等。

根据本发明的一个方面，其提供一种用于估计病人脑状态的至少一个振荡分量的方法。该方法包括：接收脑电图(EEG)数据；根据所述脑电图(EEG)数据，拟合一状态空间模型，所述状态空间模型包括至少一个振荡分量；根据所述状态空间模型，估计所述脑电图(EEG)数据的跨频耦合的至少一个值；根据所述跨频耦合的至少一个值，产生一报告；以及将所述报告输出到一显示器及一存储器的至少其中之一。

在所述方法中，所述跨频耦合的至少一个值可包括所述状态空间模型中包括的一第一波分量与一第二波分量之间的相位振幅调制的一估计值。

在所述方法中，所述第一波可以是一慢速波分量，而所述第二波可以是一alpha波分量。

在所述方法中，所述相位振幅调制的估计值可根据所述慢速波分量的一相位及所述alpha波分量的一振幅来计算。

在所述方法中，所述跨频耦合的至少一个值可包括所述状态空间模型中包括的一振荡分量与所述脑电图(EEG)数据的一频带之间的一相关值。

在所述方法中，所述拟合所述状态空间模型的步骤可包括拟合所述至少一个振荡分量中包括的每一振荡分量的谐振频率。

在所述方法中，所述至少一个振荡分量可包括一alpha分量，而所述alpha分量与一中心频率的先验分布有关。所述先验分布可以是一“冯·米塞斯”(Von Mises)先验分布。

在所述方法中，所述状态空间模型的一阻尼因数可用一先验分布来约束。所述先验分布可以是一beta分布。

根据本发明的另一个方面，其提供一种用于估计病人脑状态的至少一个振荡分量的***。该***包括：多个传感器，所述多个传感器用于接收脑电图(EEG)数据；一处理器，所述处理器用于接收所述脑电图(EEG)数据，并执行包括“根据所述脑电图(EEG)数据，拟合一状态空间模型，所述状态空间模型包括至少一个振荡分量；根据所述状态空间模型，估计所述脑电图(EEG)数据的跨频耦合的至少一个值；根据所述跨频耦合的至少一个值，产生一报告”等步骤；以及一显示器，所述显示器用于显示所述报告。

在所述***中，所述跨频耦合的至少一个值可包括所述状态空间模型中包括的一第一波分量与一第二波分量之间的相位振幅调制的一估计值。

在所述***中，所述第一波可以是一慢速波分量，而所述第二波可以是一alpha波分量。

在所述***中，所述相位振幅调制的估计值可根据所述慢速波分量的一相位及所述alpha波分量的一振幅来计算。

在所述***中，所述跨频耦合的至少一个值可包括所述状态空间模型中包括的一振荡分量与所述脑电图(EEG)数据的一频带之间的一相关值。

在所述***中，所述拟合所述状态空间模型的步骤可包括拟合所述至少一个振荡分量中包括的每一振荡分量的谐振频率。

在所述***中，所述至少一个振荡分量可包括一alpha分量，而所述alpha分量与一中心频率的先验分布有关。所述先验分布可以是一“冯·米塞斯”(Von Mises)先验分布。

在所述***中，所述状态空间模型的一阻尼因数可用一先验分布来约束。所述先验分布可以是一beta分布。

本发明的前述优点及其他优点将出现在以下说明中。在所述说明中，参阅了构成其一部份的附图，所述附图以图解方式描绘本发明的优选实施例。这样的实施例并非必定代表本发明的全部范围；因此，为解释本发明的范围，需参考权利要求及此中内容。

附图说明

图1A显示一生理检测***10的实施例。

图1B显示一示例性脑电图(EEG)传感器。

图2显示一用于病人监护的示例性***。

图3显示一用于获取数据及产生一报告的示例性流程。

图4显示另一用于获取数据及产生一报告的示例性流程。

图5A显示接受丙泊酚治疗的人的一自回归模型的脑电图(EEG)频谱。

图5B显示一自回归模型在休息状态期间的脑电图(EEG)频谱。.

图6A显示一组合振荡时间序列。

图6B显示图6A的一慢速振荡时间序列。.

图6C显示图6A的alpha振荡时间序列。

图7显示带通慢速分量及带通alpha分量的多锥频谱。

图8A显示一带通方法及一状态空间方法的总信号。

图8B显示一带通方法及一状态空间方法的慢速振荡分量。

图8C显示一带通方法及一状态空间方法的alpha振荡分量。

图8D显示一带通方法及一状态空间方法的残差。

图9A显示脑电图(EEG)数据的自回归(AR)分量的频谱。

图9B显示所述脑电图(EEG)数据的估计振荡的频谱。

图10A显示没有“冯·米塞斯”(Von Mises)先验分布的一慢速及alpha模型的自回归(AR)分量。

图10B显示没有“冯·米塞斯”(Von Mises)先验分布的所述慢速及alpha模型的估计振荡的频谱。

图10C显示一慢速无alpha模型的自回归(AR)分量的频谱。

图10D显示所述慢速无alpha模型的估计振荡的频谱。

图11显示一示例性流程1100，该流程1100以脑电图(EEG)数据拟合一状态空间模型。

图12A显示原始脑电图(EEG)数据。

图12B显示所述脑电图(EEG)数据的一个六秒时窗。

图12C显示在所述六秒时窗中拟合的一慢速振荡。

图12D显示在所述六秒时窗中拟合的一快速振荡。

图12E显示所述状态空间模型的一示例性形式。

图12F显示一具有可信区间的慢速振荡相位。

图12G显示一具有可信区间的快速振荡振幅。

图12H显示一估计K^mod及相关分布。

图12I显示一估计φ^mod及相关分布。

图12J显示已回归的alpha波振幅。

图13A显示已知数据的一参数功率谱密度数据的一第一带通拟合。.

图13B显示所述参数功率谱密度数据的一阈值。

图13C显示所述参数功率谱密度数据的一第二带通拟合。

图13D显示一矫正功率谱密度数据。.

图13E显示参数功率谱密度数据的一拟合线。

图14A显示供应予一病人脑电图(EEG)数据的丙泊酚浓度。

图14B显示一对应于所述病人脑电图(EEG)数据的响应概率。

图14C显示所述病人脑电图(EEG)数据的调制回归值。

图14D显示所述病人脑电图(EEG)数据的参数谱。

图14E显示所述病人脑电图(EEG)数据的“相位振幅调制图”(PAM)。

图14F显示所述病人脑电图(EEG)数据的调制参数。

图15A显示一响应概率图。

图15B显示一丙泊酚浓度图。

图15C显示一使用六秒时窗的标准方法的调制图。

图15D显示一与图15C相关的调制指数图。

图15E显示一使用一百二十秒时窗的标准方法的调制图。

图15F显示一与图15E相关的调制指数图。

图15G显示一用于“状态空间跨频分析”(SSCFA)方法的调制图。

图15H显示一与图15G相关的调制指数图。

图15I显示一用于“双状态空间跨频分析”(dSSCFA)方法的调制图。

图15J显示一与图15I相关的调制指数图。

图16A为一响应概率图。图16B为一丙泊酚浓度图。

图16C为一使用六秒时窗的标准方法的调制图。

图16D为一与图16C相关的调制指数图。

图16E为一使用一百二十秒时窗的标准方法的调制图。

图16F为一与图16E相关的调制指数图。

图16G为一用于“状态空间跨频分析”(SSCFA)方法的调制图。

图16H为一与图16G相关的调制指数图。

图16I为一用于“双状态空间跨频分析”(dSSCFA)方法的调制图。

图16J为一与图16I相关的调制指数图。

图17A显示使用一先前方法对振荡分量的分量进行分解。

图17B显示使用图17A的先前方法确定的功率谱密度。

图17C显示一通过图17A的先前方法进行的相位振幅互相关性估计。

图17D显示使用一“状态空间跨频分析”(SSCFA)方法对振荡分量的分量进行分解。

图17E显示使用图17D的“状态空间跨频分析”(SSCFA)方法确定的功率谱密度。

图17C显示一通过先前方法进行的相位振幅互相关性估计。

图17F显示一通过图17D的“状态空间跨频分析”(SSCFA)方法进行的相位振幅互相关性估计。

图18A显示振幅调制的动力学状态。

图18B显示相位调制的动力学状态。

图19A显示一用于额部电极的基于互相关性的跨频耦合。

图19B显示一组合调制图概要。

图20A显示跨频耦合主频众数为频率的函数。

图20B显示所述主频众数的总能量百分比。

图21显示映射到一第一主频众数的跨频耦合图案。

图22显示一在不同状态期间通过大脑不同部位处的低频活动映射到所述第一主频众数的宽带频率耦合图。

图23显示一使用拟合状态空间模型来执行跨频耦合分析的示例性流程。

具体实施方式

如下文所述，本发明提供一示例病人监护***及使用方法。作为一非限制性例子，如下文所述，所述***包括可用于提供对病人进行生理监测，包括对被施用至少一种具有麻醉特性的药物的病人进行的监测。

例如，如下文所述的一种用于所述***及方法的非限制性应用乃为监测全身麻醉、镇静状态或其他药物引起的状态。因此，根据本发明的一个方面，其提供一种用于监测被施用至少一种具有麻醉特性的药物的病人的***。所述***可包括：多个传感器，所述多个传感器在所述病人接受至少一种具有麻醉特性的药物时用于获取所述病人的生理数据；以及至少一个处理器，所述至少一个处理器用于获取来自所述多个传感器的生理数据，并根据所述生理数据，确定可用于确定所述病人在全身麻醉或镇静状态或其他药物引起的状态期间的信号标记。

所述***及方法亦可提供给药的闭环控制，以产生及调节全身麻醉、镇静状态或其他药物引起的状态所需要的大脑状态。因此，根据本发明的一个方面，其提供一种用于对被施用至少一种具有麻醉特性的药物的病人进行闭环控制的***。所述***可包括：多个传感器，所述多个传感器在所述病人接受至少一种具有麻醉特性的药物时用于获取所述病人的生理数据；以及至少一个处理器，所述至少一个处理器用于获取来自所述多个传感器的生理数据，并根据所述生理数据，确定可用于推断所述病人在全身麻醉或镇静状态或其他药物引起的状态期间的大脑状态。所述***亦可包括药物传输***，以根据所述监测器提供的大脑状态信息进行给药。

现参阅附图，图1A显示一生理监测***10的例子。在所述生理检测***10中，一内科病人12由一个或多个传感器13进行监测，每个传感器13通过一电缆15或其他通信链路或媒介，向一生理监测器17传输一信号。所述生理监测器17包括一处理器19以及(可选择地)一显示器11。所述一个或多个传感器13包括传感元件，比如脑电图(EEG)传感器或类似物。所述传感器13可通过测量所述病人12的生理参数，产生相应的信号。所述信号接着由一个或多个处理器19处理。所述一个或多个处理器19接着将已处理的信号传输到所述显示器11(如果提供一个显示器11)。在一个实施例中，所述显示器11并入所述生理监测器17中。在另一个实施例中，所述显示器11与所述生理监测器17分开。所述生理监测***10是在一个配置中的一种便携式监测***。在另一个例子中，所述生理监测***10是一个吊舱，没有显示器，而且适合向一显示器提供生理参数数据。

为了表达清楚起见，图1A中显示的一个或多个传感器13以一单块描绘。应该理解，图中所示的传感器13旨在代表一个或多个传感器。在一实施例中，所述一个或多个传感器13包括以下描述的多种类别的其中之一的单一传感器。在另一实施例中，所述一个或多个传感器13包括至少两个脑电图(EEG)传感器。在又另一实施例中，所述一个或多个传感器13包括至少两个脑电图(EEG)传感器以及一个或多个脑氧合传感器等。在前述的每一实施例种，亦可选择地包括不同类别的额外传感器。其他不同数目或类别的传感器组合亦适合用于所述生理监测***10。

在图1A中显示的***的一些配置中，用于接收及处理来自所述传感器的信号的所有硬件安装在同一外壳中。在其他实施例中，一些用于接收及处理信号的硬件安装在一独立的外壳内。此外，某些实施例的所述生理监测器17包括硬件、软件或硬件及软件(安装在一个或多个外壳中)，用于接收及处理由所述传感器13传输的信号。

如图1B中所示，所述脑电图(EEG)传感器13可包括一电缆25。所述电缆25可包括一电气屏蔽体中的三根导线。一根导线26可提供电力至一生理监测器17，一根导线28可提供一接地信号至所述生理监测器17，而一根导线28可将来自所述传感器13的信号传输至所述生理监测器17。对于多个传感器而言，可提供一根或多根额外电缆15。

在一些实施例中，所述接地信号为一接地电极信号；但在其他实施例中，所述接地信号为一病人接地电极信号(有时称为病人参考、病人参考信号、回程信号或病人回程电极信号。在一些实施例中，所述电缆25在一电气屏蔽层中携带两根导线，而所述屏蔽层充当接地导线。所述电缆25中的电气接口23可使所述电缆能够与所述生理监测器17的一连接器20中的电气接口21电气连接。在另一实施例中，所述传感器13与所述生理监测器17无线地通信。

在一些配置中，图1A及图1B中显示的***可进一步包括一存储器、数据库或可由所述处理器19访问的其他数据存储位置(图中未显示)，以存储参考信息或其他数据。具体地说，这样的参考信息可包括一列表的病人类别(比如年龄类别)，以及相关信号、信号标记或信号特征。例如，信号标记或信号特征可包括各种信号振幅、相位、频率、功率谱或它们的范围或组合，以及其他信号标记或信号特征。在一些方面，这样的参考信息可由所述处理器19用来确定特定病人特性，比如明显或可能的病人年龄或其他病人状况或类别。在一些方面，数据获取可根据所选择及/或所确定的病人特性来调整或修改。

现专门参阅图2，其描绘一根据本发明的示例性***200，所述***200可构建为一独立脑监测装置或可携带装置，或可并入作为一现有脑监测装置的一中心部件。从即将到来的描述中将可了解，在与进行各种医疗程序(比如麻醉剂的施用)相关的手术室或重症加护环境中，以及在手术前或手术后评估的情况中，可找到所述***200的有价值用途。

所述***200包括一病人监测装置202(比如一生理监测装置)，其在图2中描绘为一脑电图(EEG)电极阵列。然而，预期所述病人监测装置202可包括多个不同的传感器。明确地说，所述病人监测装置202亦可包括用于监测皮肤电反应(GSR)(例如测量对外部刺激的唤醒度)的机制，或其他监测***(比如心血管监测器及血压监测器)，包括心电图及血压监测器，以及眼球微震监测器。这个设计的一个实现乃利用一具有额外电极的额部“拉普拉斯脑电图”(Laplacian EEG)电极布置，以测量皮肤电反应(GSR)及/或眼球微震。这个设计的另一实现可能并入一具有独立皮肤电反应(GSR)电极、可在后处理中结合的额部电极阵列，以获得被发现能最佳地检测前述脑电图(EEG)信号特征的任何电极组合。这个设计的另一实现可能利用一具有独立皮肤电反应(GSR)电极、以64至256个传感器采样整个头皮表面的高密度布置，以达致源定位目的。

所述病人监测装置202通过一电缆204连接，与一监测***206通信。此外，所述电缆204及类似连接可由组件之间的无线连接取代。所述监测***206可用于接收来自所述病人监测装置202的原始信号(比如由所述脑电图(EEG)电极阵列获取的信号)，并以各种形式收集、处理甚至显示所述信号，包括时序波形、频谱图等等。在一些操作模式中，所述监测***206可设计成获取来自所述病人监测装置202上的传感器的、以生理数据或其他数据为形式的侦察数据，并使用所述侦察数据来识别其中的信号标记或信号特征。例如，可以使用各种适合方法，在侦察数据及其他获取的数据中识别信号振幅、相位、频率、功率谱及其他信号标记或信号特征。此外，可执行多锥分析来识别及列入一动态范围的跨越几个数量级的信号。这样的信号标记或信号特征可以接着由所述监测***206用来确定各种病人特性，包括明显及/或可能的病人年龄。

在一个实施例中，使用所述监测***206来进行生理数据的获取操作，可按照根据侦察数据确定的病人特性来调整或调节。明确地说，所述监测***206可用于确定一与某些已确定的病人特性一致的尺度，以及根据所确定的尺度及/或由用户提供的任何指示来调整后续的数据获取。例如，可以通过调整一个或多个放大器增益以及其他数据获取参数，对数据获取进行调节。此外，在一些方面，所述监测***206可以进一步用于格式化所获取的生理数据，以按所述尺度显示。以这种方式，可根据所述明显及/或可能的病人年龄来确定一年龄适合尺度，而且使用一选定年龄适合尺度的任何后续数据获取将产生并描绘年龄补偿数据。

如图所示，所述监测***206可进一步连接至以专用分析***208。然而，所述监测***206及所述分析***208可以集成或合并成一共同***。所述分析***208可接收来自所述监测***206的脑电图(EEG)波形，以及(如下文所述)分析所述脑电图(EEG)波形及其中的信号特征。然而，预期所述监测***206及所述分析***208的任何分析或处理功能，可按需要或期望共享或单独分配。

在一些方面，涉及一接受特定医疗程序的病人的已确定特性相关的信息可以提供予所述***200的临床医师或操作员。例如，此前发现年老病人更容易在手术室中进入爆发抑制。明确地说，爆发抑制是深度的脑失活，其中脑电活动电的爆发不时***等电周期，故称为抑制。麻醉剂引起无意识的大脑状态，由alpha波(8～10Hz)及慢速波(0.1～4Hz)信号振荡决定，这些大脑状态可以以少于产生爆发抑制所需要的麻醉剂剂量达致。这可能意味着，需将麻醉剂剂量降低至大大低于当前推荐给老年人的水平。由于当前推荐的剂量一般将老年病人置于爆发抑制状态，通过根据实时的脑电图(EEG)监测来滴定麻醉剂剂量，可以实现适当的全身麻醉及减低麻醉暴露的状态。因此，所述***200可根据所确定的病人特性，提供信息用于选择适当的麻醉剂剂量。以这种方式，(例如)可减少接受全身麻醉的老年病人发生手术后认知功能障碍。

在另一例子中，所述监测***206及/或所述分析***208可提供对特定病人(比如年轻、中年、老年以及毒品成瘾病人)的手术前或手术后评估，以确定先验信息，这些先验信息可用于识别及/或预测特定病人状况(包括麻醉剂敏感性)以及任何手术后并发症可能(比如认知功能障碍)。此外，亦可提供麻醉护理、麻醉后护理或重症加护的特殊疗法。

所述***200亦可包括一药物传输***210。所述药物传输***210可连接至所述分析***208及所述监测***206，使所述***200形成一闭环监测及控制***。这样的根据本发明的闭环监测及控制***能够进行范围广泛的操作，但包括用户界面212，以允许用户配置所述闭环监测及控制***、接收来自所述闭环监测及控制***的反馈、以及在需要时重新配置及/或否决所述闭环监测及控制***。在一些配置中，所述药物传输***210不仅能够控制麻醉剂化合物的施用，以便将病人置于受所述麻醉剂化合物影响的意识减退状态(比如全身麻醉或药物镇静)，但亦可实施及反映多个***及方法，以使病人进入或脱离意识增强或意识减弱的状态。

例如，根据本发明的一个方面，甲哌酯(MPH)可以作为多巴胺与多种去甲肾上腺素再摄取转运体的一种抑制剂，并积极诱导异氟醚全身麻醉后苏醒。甲哌酯(MPH)可用于恢复意识、诱导与唤醒度相一致的脑电图变化以及增加呼吸动力。甲哌酯(MPH)对行为及呼吸的影响可被氟哌利多抑制，这支持甲哌酯(MPH)是通过激活多巴胺能唤醒途径来诱导唤醒的证据。体积描记术及血气实验确认甲哌酯(MPH)增加每分钟通气量，而增加每分钟通气量增加麻醉剂从大脑中排出的速率。此外，可使用一控制***(比如上述控制***)，通过提升唤醒度，将乙哌酯(EPH)或其他药剂用于积极诱导从异氟醚、丙泊酚或其他全身麻醉后苏醒。

因此，可以通过包括一具有两个特定子***的药物传输***210，提供一***(比如以上所述与图2有关的***)，以实行麻醉后主动苏醒。因此，所述药物传输***210可包括一麻醉剂化合物给药***224，所述麻醉剂化合物给药***224被设计成将一种或多种麻醉剂化合物的剂量传输至一麻醉对象，而且亦可包括一苏醒化合物给药***226，所述苏醒化合物给药***226被设计成传输会消除全身麻醉或促进一麻醉对象从麻醉中自然苏醒的一种或多种化合物剂量。

例如，甲哌酯(MPH)与类似物及其衍生物通过提升唤醒度及呼吸动力，诱导一麻醉对象从麻醉诱导的无意识中苏醒。因此，所述苏醒化合物给药***226可用于传输甲哌酯(MPH)、***、***、氨酚烷胺或咖啡因，以便在手术结束时消除全身麻醉诱导的无意识及呼吸抑制。所述甲哌酯(MPH)可以是右旋甲哌酯(D-MPH)、外消旋甲哌酯、或左旋甲哌酯(L-MPH)，或可以是相等或不同比率的合成物，比如大约50％：50％、或大约60％：40％、或大约70％：30％、或80％：20％、90％：10％、95％：5％等等。其他药剂可以作为一剂量比用于治疗注意力缺陷障碍(ADD)或注意力缺陷多动障碍(ADHD)的剂量更大的甲哌酯(MPH)给药，比如甲哌酯(MPH)的剂量可以是介于大约10mg/kg与大约5mg/kg之间，以及大约5mg/kg与10mg/kg之间的任何整数。在某些情况下，所述剂量介于大约7mg/kg与大约0.1mg/kg之间，或介于大约5mg/kg与大约0.5mg/kg之间。其他药剂可包括那些吸入给药者。

现参阅图3，其中显示根据本发明的流程300。从流程块302处开始，可以获取任何数量的生理数据，其中所述生理数据代表生理信号，比如脑电图(EEG)信号，这些生理信号乃通过使用(例如)所述病人监测装置202，从一病人处获取。在一些方面，所述生理数据可包括侦察数据，这些侦察数据用于多种目的，包括确定各种病人特性。然后在流程块304处，使用所获取的生理数据来识别或确定信号标记或信号特征。例如，可以使用各种适合方法，在侦察数据及/或其他获取的数据中识别信号振幅、相位、频率、功率谱及其他信号标记或信号特征。

在一些优选实施例中，所述信号标记或信号特征可用于确定病人特性，包括明显及/或可能的病人年龄。此外，流程块304亦可包括确定一与已确定的病人特性一致的尺度的步骤。在一个方面，使用谱估计方法(比如多锥方法)可天然地列入一动态范围的跨越许多数量级的信号。在另一方面，可执行信号振幅的自动估计，以推断一用于可视化尺度以及获取放大器增益的正确年龄群组及护理者设置。

在下一流程块306处，使用根据所述侦察数据确定的信号标记或信号特征，可调整或调节一涉及后续获取的信号数据的数据获取流程。例如，可以通过调整一个或多个放大器增益以及其他数据获取参数，对数据获取进行调节。在一些方面，调节数据获取亦可包括使用一与所确定的病人特性一致的尺度，以及根据所确定的尺度及/或由用户提供的任何指示来调整后续的数据获取。例如，可以在流程块304处使用一根据所述明显及/或可能的病人年龄来确定的年龄适合尺度，而且使用一选定年龄适合尺度的任何后续数据获取将产生年龄补偿数据。

在流程块308处，以一方式获取的数据可用于确定病人的当前或未来大脑状态。例如，使用年龄补偿数据收集的已分析或已处理的脑电图(EEG)波形可用于评估一当前及/或未来的麻醉或药物镇静深度。此外，确定这样的大脑状态亦可包括由一临床医师或用户提供的任何信息，比如涉及一医疗程序的信息。

接着，在流程块310处，产生一报告，例如以打印报告或(优选)实时显示为形式的报告。所述报告可包括原始或经处理的数据、信号特征信息、当前或未来大脑状态的迹象以及涉及病人特异性特点的信息，包括可能及/或明显的病人年龄。已显示的信号特征信息或已确定的状态可能以波形、频谱、相干图、概率曲线等为形式。在一些方面，所述报告可包括按一尺度显示的格式化生理数据。在其他方面，所述报告可指示一麻醉敏感性、一手术后并发症(比如认知功能障碍)概率，以及用于麻醉护理、麻醉后护理或重症加护等的疗法。

参阅图4，其描绘根据本发明的另一流程400的步骤。明确地说，所述流程400在流程块402处开始，在该处使用(例如)此中所述的病人监测***来获采样本数据或侦察数据。在流程块404处，接着使用调整或参考类别来分析所述样本数据，以识别代表所获取的样本数据的病人类别。明确地说，所述步骤可包括识别所述样本数据中的信号标记或信号特征，以及与所述参考类别相关的信号标记户特征进行比较。例如，可以使用各种适合方法，在所述样本数据中检测信号振幅、相位、频率、功率谱及其他信号标记或信号特征。

在流程块404处执行的分析可指示特定病人特性，包括明显或可能的病人年龄。在一些方面，一指示特定病人特性的已识别或明显的类别可选择地在流程块406处显示。此外，在流程块408处，亦可接收一用户输入。

随后，在流程块410处，对各种通信参数进行确定。这包括考虑已确定或已推断的病人特性或类别，以及可选择地考虑一用户输入。例如，在流程块410处，可根据存在于所获取的数据中的已确定病人特性及/或信号、信号标记或信号特征，确定用于所获取数据的年龄适合尺度。接着，在流程块412处，可使用所确定的通信参数来调节后续的数据获取，以获取年龄适合数据。如此中所述，调节数据获取的步骤可包括使用所述通信参数适当地调整或修改各种放大器增益。在一些方面，所确定的通信参数可直接地应用于所获取的样本数据。例如，可应用一年龄适合尺度于所述样本数据，以创建年龄适合数据。

接着，在流程块414处，以一方式获取或处理的数据可用于确定病人的当前或未来大脑状态。例如，使用年龄补偿数据收集的已分析或已处理的脑电图(EEG)波形可用于评估一当前及/或未来的麻醉或药物镇静深度。此外，确定这样的大脑状态亦可包括由一临床医师或用户提供的任何信息，比如涉及一医疗程序的信息。

接着，在流程块416处，产生任何适当形状或形式的报告。在一些方面，所述报告可显示缩放数据或描述所述数据的数据类别。在其他方面，所述报告可指示一麻醉敏感性、一手术或手术后并发症概率、一明显或可能的病人年龄以及涉及本发明的其他信息。

本发明的***及方法可以以许多方式来提供大脑状态指标，这些大脑状态指标可用于评估神经信号，以导出对病人监测有用的指标或其他信息。一个这样的非限制性例子可包括：

确定及/或分析慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)的振幅的特性。

确定及/或分析慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)中的分析能力的特性。

确定及/或分析振荡的峰值频率，包括确定及/或分析慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)的峰值频率。

确定及/或分析不同传感器之间或多个传感器之间的相干性特征，包括确定及/或分析用于慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)的不同传感器之间或多个传感器之间的相干性特性。

确定及/或分析慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)的相位特性。

确定及/或分析用于慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)的不同传感器之间或不同传感器组之间的相位差的特性。

确定及/或分析多次振荡之间的跨频耦合关系的特性，例如一振荡器的相位与其他振荡器的振幅之间的关系，包括确定及/或分析多次慢速振荡(0.1～1Hz)、delta振荡(1～4Hz)、theta振荡(4～8Hz)、alpha振荡(8～12Hz)、beta振荡(12～25Hz)、gamma振荡(25～70Hz)以及高gamma振荡(70Hz及以上)之间的跨频耦合关系的特性。

在一个例子中，为了监测及跟踪全身麻醉或药物镇静，可呈现一跨越所述theta频段(4～8Hz)、alpha频段(8～12Hz)及beta频段(12～25Hz)的振荡的峰值频率。与此同时，可以呈现所述慢速振荡频段(0.1～1Hz)及delta振荡频段(1～4Hz)。这些由在此呈现的***及方法提供的具体指标，可在一***物(比如丙泊酚或七氟醚)施加时，由(例如)一麻醉师或其他临床医师用来辨别对应于慢速振荡/delta及alpha振荡的全身麻醉状态，这些全身麻醉状态不同于表现beta及gamma振荡特征的药物镇静状态。

在另一例子中，为了监测及跟踪全身麻醉或药物镇静，可呈现一慢速(0.1～1Hz)至delta(1～4Hz)振荡的相位与一alpha(8～12Hz)至beta(12～25Hz)振荡的振幅之间的相位振幅耦合。该关系的相位可用于识别一镇静状态(例如“troughmax”调制)或一深度无意识状态(例如“peakmax”调制)之间的差异。所述耦合关系的强度可在闭环控制之下，用来调节维持深度无意识状态而需传输的麻醉剂用量。

在另一例子中，本发明的***及方法可用于表达与感兴趣的生理数据不同的次级生理信号，比如所述脑电图(EEG)及脑电图(EEG)导出量。例如，可将一低频一阶自回归(AR(1))分量添加到所述模型的状态动力学分量，以列入低频漂移。同样地，亦可包括移动诱导的次级生理波形的状态空间模型。同样地，可将一自回归或移动平均分量并入所述模型的状态动力学分量，以表达肌肉诱导的次级生理信号。以这种方式，本发明的***及方法可用于使次级生理信号与感兴趣的生理信号分离。作为不同状态的指标(比如意识或移动)，这些次级生理信号本身可能有用。

在另一操作模式中，所述***可使用来自其它装置的信息，比如使用来自其他生理监测器或来自所述医疗记录***的信息，以告知记录自脑电图(EEG)或生理传感器的大脑状态或生理状态的解释。这个信息可告知特定振荡器模型以及特定次级生理模型的使用，以提供所期望的监测变量。

在另一操作模式中，所述***可使用旨在表达不同临床情景的一数据库的模型及模型参数化。所述***可使用来自医疗记录的信息，比如病人人口统计信息或病史，以及关于施用的***物的信息以及关于正被执行的外科病例类别的信息，以选择可应用于已知情景的一组候选模型。所述***可使用所述脑电图(EEG)或被监测的特定病人的其他生理数据，以选择能最好地应用于已知情景的一个或多个模型。所述数据库的模型可以由许多不同形式，包括不同参数的先验概率分布，这些先验概率分布表示为不同输入变量的函数或表示为查阅表形式或类似概念。

一般而言，一脑电图(EEG)的振荡可表示为在噪声中观测到的不同分量。例如，在丙泊酚麻醉下的一慢速振荡及一alpha振荡可使用以下方程式建模：

y_t＝x_slow，t+x_alpha，t+v_t (1)

以下描述一个更一般的案例。每一振荡分量接着可以使用一个二阶自回归(AR(2))过程来建模，该过程可以以状态空间形式书写。值得注意的是，由于观测的存在(下文将对此进行描述)，在技术上，每一振荡分量的模型乃为自回归移动平均(ARMA)模型，但我们将根据它们的基本状态动力学来称谓这些模型，以强调模型的振荡结构。现在可以假设所观测的信号

是依次表达一慢速分量及一快速分量的两种潜在状态

及

的线性组合。这些二维(2D)潜在状态中的每一潜在状态可以假设为独立，而且它们在一固定步长上的演化建模为一缩放及噪声旋转，因为j＝s,f：

其中，a_j为缩放参数，

为角ω_j(径向频率)的二维旋转，而

为过程噪声方差。图12A-D中显示这个状态空间振荡分解的一个例子。由于所述慢速分量及快速分量时在所述模型下直接估计，这个方法消除对带通滤波的需要。或许更重要的是可以将所述两种振荡各自的分量看作是一相量信号或解析信号的实项及虚项。因此，不再需要使用在先前方法中用来合成来自所观测的实信号的所述虚信号分量的希尔伯特变换(Hilberttransform)。因此，所述潜在向量的极坐标提供所述慢速瞬时相位

及所述快速振荡振幅

的直接表示方式(图12F及图12G)。注意

并取得一标准状态空间表示方式：

其中，

为一由前述的旋转组成的分块对角矩阵，Q为总过程噪声协方差，R为观测噪声协方差，而

将所述观测与所述振荡的第一(实)坐标联系起来。通过使用一“期望最大化(EM)”程序或使用数字计算程序，可取得估计数字(Φ，Q，R)，以最大化对数似然密度或对数后验密度。为了扩展此模型，以添加更多独立振荡分量或列入一振荡的谐波(如在一非线性***中可能遇到者)，以下描述及推导出一更为普遍的公式。

以下描述方程式(3)及(4)体现的状态空间模型的更为具体的推导。对于在F_s采样的时间序列，我们考虑一时窗

而且我们在此部份假设，y_t为来自导致一慢速分量及一快速分量的两种潜在状态

及

的观测噪声及分量的和。因j＝s，f及t＝2..N，每一分量遵守以下过程方程式：

以及

其中a_j∈(0，1)及

乃为一角度为ω_j＝2πf_j/F_s的旋转矩阵。

如先前所述，每一种振荡的相位

及振幅

是使用所述潜在向量极坐标求得：

以及

每一种振荡在频率fj处带有一峰值时，具有一宽带功率谱密度(PSD)。下文将进一步推导这个宽带功率谱密度(PSD)的参数表达式。

设定M＝[1 0 1 0]、

而Q及Φ为分块对角矩阵(它们的分块依次为Q_i及

我们求得方程式(3)及(4)的标准状态空间。

已知所观测的信号y_t，我们的目标是估计这两种隐藏的振荡x_t以及它们的生成参数(Φ，Q，R)。我们是使用包括一“期望步骤”(E-step)及一“最大化步骤”(M-step)的“期望最大化(EM)算法”，下文将描述其一般推导。所述隐藏的振荡x_t在所述“期望最大化(EM)算法”的“期望步骤”(E-step)中，使用“卡尔曼滤波器”(Kalman filter)及“固定区间平滑器”来估计，而所述生成参数则在所述“最大化步骤”(M-step)的每一迭代中估计。

值得注意的是，所述二阶自回归(AR(2))形式可表达具有不同水平的共振或带宽的振荡。此外，在高于所述振荡频率的频率处，所述二阶自回归(AR(2))频谱有一“1/f”轮廓。一般而言，可能有许多不存在这两种振荡中的一种或两种的另类情景。例如，在清醒休息状态下，额部脑电图(EEG)频道将不预期有alpha振荡或慢速振荡。在此情况下，宽带“1/f”动力学可能存在，而且可使用一阶自回归(AR(1))过程建模。可选择地，一阶自回归(AR(1))及二阶自回归(AR(2))动力学的组合亦可以这个状态空间形式表达。

作为简略的表达方式，我们将称方程式(3)及(4)中描述的模型为一“AR(2+2)”模型，既然这个模型中有两个二阶自回归(AR(2))过程。所述AR(2+2)模型可拟合至脑电图(EEG)数据，以识别一alpha振荡及一慢速振荡。一阶自回归(AR(1))与二阶自回归(AR(2))过程的组合可称为一“AR(1+2)”模型。我们将考虑具有以一阶自回归(AR(1))或二阶自回归(AR(2))状态动力学表达的单一动态分量的模型。如以上所述，在技术上，由于方程式(4)中存在观测噪声，这些模型乃为自回归移动平均(ARMA)模型，但我们将根据它们的基本状态动力学来称谓这些模型，以强调模型的振荡结构。可以使用一用于状态空间模型的“期望最大化(EM)算法”(比如“Shumway及Stoffer”期望最大化算法)，将一个或多个模型(如所述AR(2+2)模型)拟合至脑电图(EEG)数据。初始自回归(AR)参数可由一用户设定。例如，一初始慢速振荡峰值频率可设定至1Hz，一初始alpha振荡峰值频率可设定至10Hz，而这两种振荡的极点可初始地设定为距离所述复平面的半径为0.99处。所述“期望最大化(EM)算法”可使用“卡尔曼滤波器”(Kalman filter)、“固定区间平滑器”以及协方差平滑算法来进一步优化估计参数。所述“期望最大化(EM)算法”的期望步骤及最大化步骤在估计参数(比如自回归AR参数及噪声协方差)输出之前，各可迭代一预定次数，比如迭代200次。所述迭代亦可继续，直到满足某些收敛性判别准则为止(包括但不限于满足被估计的参数的收敛性、或对数似然密度或对数后验密度的值中的收敛性)。在已经估计了自回归(AR)参数及噪声协方差之后，可以计算用于不同模型的“赤池信息准则”(AIC)，以比较每一模型拟合一已知脑电图(EEG)的程度优劣。所述“赤池信息准则”(AIC)定义为AIC＝-2log(L)+2p，其中L为所述模型的似然性，而p为所述模型中的参数数目。

在为准确地确定慢速振荡及alpha振荡而确定所述AR(2+2)模型的可行性而进行的某个测试时，在以丙泊酚诱导年龄介于18至36岁的健康志愿者进入无意识状态期间，记录脑电图(EEG)。所述脑电图(EEG)信号在采样频率为5000Hz时记录，并降采样至200Hz。通过对一清醒闭眼基准状态与一丙泊酚诱导无意识状态进行比较，分析了脑电图(EEG)的单一额部频道。在处于所述丙泊酚诱导无意识状态时，所述额部脑电图(EEG)含有大慢速(0.1～1Hz)及alpha(8～12Hz)振荡。在处于所述基准状态时，不存在这些振荡，虽然可观测到以慢速漂移。我们使用了所述“赤池信息准则”(AIC)，在所述AR(1)、AR(2)、AR(1+2)及AR(2+2)模型之中选择一最佳模型。为了列入每一模型的状态空间公式，我们将所述“自回归”(AR)参数、状态噪声方差以及观测噪音方差作为参数，包括在所述“赤池信息准则”(AIC)中，所有这些参数通过所述“期望最大化(EM)算法”进行估计。其他信息-比如“贝叶斯信息准则”(BIC)亦可用于代替所述“赤池信息准则”(AIC)，以选择一最佳模型。

一信号可通过使用带通滤波器分解为频率依赖分量。因此，我们对所述状态空间模型的性能与带通滤波信号进行了比较。为了表达所述慢速波，我们使用了一截止频率为1Hz的100阶汉明窗(Hamming-windowed)低通有限脉冲响应(FIR)滤波器。为了表达所述alpha波，我们使用了一通频带介于8～12Hz之间的100阶汉明窗(Hamming-windowed)低通有限脉冲响应(FIR)滤波器。通过总结带通慢速及alpha分量，创建了一重构带通信号。

以下“表I”显示所述丙泊酚诱导无意识状态的AR(1)模型、AR(2)模型、AR(1+2)模型及AR(2+2)模型的“赤池信息准则”(AIC)评分。表中显示，在所述丙泊酚诱导无意识状态期间，所述AR(2+2)模型超越其他模型，这可由其在所有模型之中为具有最低“赤池信息准则”(AIC)水平者得以证明。对于所述基准案例而言，所述AR(1)模型的“赤池信息准则”(AIC)评分最低；而对于所述丙泊酚案例而言，所述AR(2+2)模型的“赤池信息准则”(AIC)评分最低；这使我们能够为每一情况选择这些模型。图5中显示所述拟合模型的频谱，以及所述观测数据的非参数多锥频谱估计(3Hz频谱分辨率、10秒的数据内26个锥)。图5A显示所述丙泊酚案例中的AR(2+2)模型的频谱，而图5B显示所述休息状态诱导案例中的AR(1)的频谱。

表I

模型	参数数目	休息	丙泊酚
				AR(1)	3	6683.3	6664.6
AR(2)	4	6685.4	3728.8
				AR(1+2)	6	6690.8	3567.8
AR(2+2)	7	6693.3	3407.4

在所述休息状态的案例中，不存在可识别的alpha振荡，这支持选择AR(1)模型。所述建模观测数据与在休息状态中收集的数据之间的残差的大小为4.9348uV或更小。在这个AR(1)模型案例中，有一隐藏状态几乎与所记录的数据相同。用于这个基准状态的AR(1)模型是：

图5A中显示代表所述丙泊酚诱导无意识状态案例的AR(2+2)模型的频谱，以及所述观测数据的多锥频谱估计。所述AR(2+2)模型的慢速振荡分量及alpha振荡分量显示，彼等密切对应在所述多锥频谱中观测到的峰值。所述AR(2+2)模型具有代表所述慢速振荡及alpha振荡的独立隐藏状态，这使得可能在时域中将所述振荡分离。图6显示在所述AR(2+2)模型之下估计的慢速、alpha及组合振荡时间序列。明确地说，图6A显示所述组合振荡时间序列，图6B显示所述慢速振荡时间序列，而图6C显示所述alpha振荡时间序列。所述残差的大小为0.2198uV或更小。

除了方程式(5)、(6)及(7)中描述的状态空间模型的自回归(AR)分量的振荡形式之外，亦可使用传统自回归(AR)模型公式。例如，描述这个丙泊酚案例中的每一分量的自回归(AR)模型可以描述为：

在这种形式的模型之下，所述振荡的特征(比如它们的峰值频率)可以根据所述模型参数直接计算。每一AR(2)分量的以弧度计算的峰值频率可以下式求得：

其中在这种情况下，a₁为第一时滞参数，而a₂为第二时滞参数。在丙泊酚诱导下的慢速波的峰值频率为1.0118Hz，而所述alpha振荡的峰值频率为11.6928Hz。由于一AR(1)过程在零频率时具有一单一实数取值极点，根据一AR(1)过程的定义，在休息状态期间的慢速分量的峰值频率为0Hz。所述AR模型分量的极点决定所述功率谱的形状及所述峰值频率。所述极点的半径决定所述峰值频率的显著性，而且实际上决定所述振荡过程的阻尼。在这个用于丙泊酚诱导之下的慢速波振荡的AR(2)模型中，所述极点的半径是0.99，而在所述alpha模型中，所述极点半径是0.956，这果然显示所述慢速波比alpha波稍微显著一点。

所述状态空间模型的性能与所述更为传统的带通滤波方法进行了比较。图7显示所述带通慢速分量及带通alpha分量的多锥频谱、所述重构信号、所述残差，以及所述状态空间模型的残差。虽然带通滤波识别与使用所述状态空间模型获得的振荡相似的振荡，但在所述带通方法之下的残差显示显著的低频振荡结构。相比之下，来自所述状态空间模型的残差小了大约60dB，而且不甚结构化。在所述带通模型中，所述残差中的功率及结构是所述任意带通截止之外的泄漏造成，而我们的方法通过使用AR模型来表达所述成分振荡来避免这种泄漏。图8显示用于时域中的带通方法的所述成分振荡、重构信号、残差以及与所述状态空间的比较。图8A显示所述带通方法及所述状态空间方法的总信号；图8B显示所述带通方法及所述状态空间方法的慢速振荡分量；图8C显示所述带通方法及所述状态空间方法的alpha振荡分量；而图8D显示所述带通方法及所述状态空间方法的残差。余所述频域分析类似，我们在带通滤波之下看见显著的时域残差结构，这种时域残差结构在所述状态空间模型之下不存在。因此，所述状态空间模型可以用于识别没有某些缺点(比如基于带通方法的相对较大的残差)的振荡。

我们获得的结果说明，由一阶及二阶***的结合体组成、以一特定形式构建的状态空间模型可如何用于识别神经时间序列中的振荡结构。在神经科学分析中采取的一种典型方法是：计算一频谱带内的功率或一带通滤波信号中的方差。这个方法的一个主要局限是，跨越感兴趣频率范围的宽带噪声不能与振荡区分开来。我们的方法通过明确地对潜在振荡分量进行建模，然后接着在具有不同振荡配置的一组模型(即：AR(1)、AR(2)、AR(1+2)及AR(2+2))中进行选择以确定所述信号的组成，从而解决了这个限制。在这种方法下，归因于一特定振荡分量的功率可以在其他宽带分量之外独立计算，而且所述分量振荡时间序列可以分开。

这个状态空间方法通过指定一用于所述自回归(AR)参数的结构来表达两个漂移(具有实极点的AR(1)或AR(2))中的任何一个漂移或振荡项(具有复极点的AR(2))，而不同于传统自回归(AR)或自回归移动平均(ARMA)建模。虽然一更普遍的自回归(AR)或自回归移动平均(ARMA)模型可定可识别这样的分量，但我们的经验是：这样的模型的频谱表达式可能随着模型阶次增加而高度变化。通过为推定的漂移项或振荡认定一特定形式，我们可以减少表达这些特征而需要的模型参数数目，并增强所述模型评估的鲁棒性。此外，在特定振荡分量方面而言，我们提议的模型的结构比一更普遍的自回归(AR)或自回归移动平均(ARMA)结构更容易解释。

相对于某些方法，所述状态空间方法具有一优势，其原因是振荡分量可在时域中分离。由于所述自回归(AR)分量默示的时间结构或脉冲响应特定于振荡器，我们已描述的线性***方法亦可在检测振荡时提供更高的特异性。此外，所述状态空间方法可用于未处理(即：未经带通滤波)的脑电图(EEG)数据，从而避免带通滤波的一些缺点。

用于所述“模型参数”的“先验分布”。在以模型拟合某些脑电图(EEG)数据时，如果所观测的振荡的功率低或无法区别于较宽的频谱景观，一模型可能无法准确地拟合所观测的振荡中的一些振荡。使用先前描述的方法，用户对所述模型的唯一控制是选择初始化参数。在所述拟合程序期间，意在表达一低功率振荡的分量可转移到另一个振荡分量描述的频率范围，尤其是如果该分量大，而且是宽带。当以传统模型选择程序(比如“赤池信息准则”(AIC))评估时，这些过拟合模型的AIC值可能较低。在某种程度上，这是因为“赤池信息准则”(AIC)倾向于偏袒具有较多分量的状态空间模型。以上描述的状态空间模型可能对初始化敏感，而且基于关于正在研究的神经振荡的先验知识，该状态空间模型可能需要附加限制。

作为一第一示例性案例，使用一AR(2+2)状态空间模型来拟合具有强alpha振荡的脑电图(EEG)数据。所述脑电图(EEG)数据的自回归(AR)分量的频谱显示于图9A，而所述估计振荡的频谱则显示于图9B。所述状态空间模型倾向于很好地拟合所述具有强alpha振荡的脑电图(EEG)数据。

相反地，当不存在alpha振荡时，一模型可能未能成功地拟合。在不存在alpha振荡的情况下，所述alpha分量不正确地拟合至所述慢速/delta频率范围。由于没有振荡功率来将所述模型的所述alpha分量(8～12Hz)锚定到预期的频率范围，在所述模型拟合程序期间，这个分量的频率向较大的慢速波分量下降。用于这个模型的“赤池信息准则”(AIC)比仅具有一慢速分量的模型为低，而且将表明应选择所述“alpha+慢速”模型(AR(2+2))。然而，既然所述(非参数)频谱估计中没有明显的alpha振荡，这与我们的直觉不符。因此，所述alpha分量不再存在于所述alpha频率额定值中，而是迁移到所述慢波频率范围中，其峰值频率为1.0872Hz。图10A中显示用于一慢速及alpha模型的自回归(AR)分量的频谱(所述慢速及alpha模型在中心频率上没有先验分布)，下文将对此进行解释。图10B中显示用于所述慢速及alpha模型的估计振荡的频谱(所述慢速及alpha模型具有所述先验分布)。图10C中显示用于一慢速无alpha模型的自回归(AR)分量的频谱。图10D中显示用于所述慢速无alpha模型的估计振荡的频谱。

为了使这类模型更易于识别，以及使它们对不适当的模型拟合更为鲁棒，先验分布可分别用于所述振荡频率及半径或阻尼因数。

对于一单一分量的中心频率上的先验概率分布而言，所述模型可能需要使用一系列的可能频率，并同时保持限定在一预定频率范围，即：将所述分量锚定到一特定频带。

我们选择了一“冯·米塞斯”(Von Mises)分布作为先验概率分布。所述“冯·米塞斯”(Von Mises)分布为具有一跨越0Hz至采样频率的域的对称概率分布。该分布的形式为：

其中ω代表我们希望以所述先验分布限制的振荡的中心频率，I₀为一0阶变形贝塞尔(Bessel)函数，而κ及μ为描述所述分布的形状的超参数。参数μ设定所述先验分布的平均数或中心频率。浓度参数κ决定所述分布的宽度或带宽，并决定对所述预期频带的锚定的强度。亦可考虑使用具有类似特性的其他概率分布，以将一振荡自回归(AR)分量锚定到以特定频率范围。由于频率分量被限制与一有限频率范围，而且因此不能在模型拟合或参数估计过程期间与其他振荡分量重叠，所述“冯·米塞斯”(Von Mises)先验分布可使模型拟合更为鲁棒。以下的测试案例显示没有所述“冯·米塞斯”(Von Mises)先验分布的模型在哪里分解：在没有所述“冯·米塞斯”(Von Mises)先验分布的情况下，“赤池信息准则”(AIC)为此对象选择所述两个组件模型，而所述慢速带中实际上只有一种振荡。值得注意的是，所选择的模型拟合所述慢速带中的两种振荡：在不存在所述“冯·米塞斯”(Von Mises)先验分布时，较高频率的alpha分量漂移入所述慢速带中，以1.0782Hz为中心。

控制中心频率及浓度的所述“冯·米塞斯”(Von Mises)先验参数，可根据所期望的频带选择。例如，可以根据神经科学文献中被普遍接受的频带，确定这些先验参数。例如，以下“表II”中的慢速及alpha模型在所述alpha频率分量上具有“冯·米塞斯”(Von Mises)先验分布(所述alpha频率分量在概率值为0.367时使用一中心频率10Hz及一带宽～1.2Hz)，以帮助将所述alpha分量锚定到所述alpha频率范围。

通过添加所述“冯·米塞斯”(Von Mises)先验分布，所述分量保持在期望范围内，而由于所述数据中缺乏alpha功率，所述分量获分配的功率很小，而且其带宽宽广。当“赤池信息准则”(AIC)实施时，其值高于仅有一种慢速振荡(无alpha振荡)的模型。因此，在使用所述“冯·米塞斯”(Von Mises)先验分布的模型中，“赤池信息准则”(AIC)选择正确的模型。这表明所述“冯·米塞斯”(Von Mises)先验分布的使用如何提高模型鲁棒性，使得模型选择准则(比如“赤池信息准则”(AIC))能够选择正确的模型。如我们在下文展示那样，这个已增强的鲁棒性也使得能够运用算法以迭代地拟合附加的振荡分量。“表II”显示用于所述不具有先验分布的慢速及alpha模型、所述慢速无alpha模型以及具有先验分布的慢速及alpha模型的“赤池信息准则”(AIC)、慢速频率、慢速半径(阻尼因数)、alpha频率以及alpha半径(阻尼因数)。

表II

由于类似原因，其他参数可能对初始化敏感，而且在这些参数上添加先验分布可提高所述模型拟合及模型选择的鲁棒性，使模型选择的自动化能够更完整。

一个这样的参数是所述半径或阻尼因数。所述半径或阻尼因数控制所述个别分量的稳定性，并因此控制整个***的稳定性。随着所述半径或阻尼因数增加超过1，所述振荡时间序列将变得不稳定，而且所述时间序列的振幅将向无穷大增加。为了使所述***保持稳定，我们可以以一beta分布控制所述阻尼因数，以使所述参数保持在0和1之间。所述beta分布的形状是由两个参数控制，并具有从0到1的域，而且可由以下方程式描述：

其中a为所述半径或阻尼因数，而α及β则为决定所述先验分布的形状的超参数。

如以上所述，一种常用于拟合自回归(AR)过程的模型选择程序是“赤池信息准则”(AIC)。所述模型拟合步骤及选择程序在以下用图解表示。所述最大化步骤及所述期望步骤的迭代包括所述“期望最大化(EM)算法”。有必要使用所述“卡尔曼滤波器”(Kalmanfilter)及“固定区间平滑器”来评估隐藏状态的时间序列(亦称估计振荡时间序列)。

【用于“独立及谐波振荡分解”的“期望最大化(EM)算法”】

现在描述一种能够估计独立及谐波振荡的示例性“期望最大化(EM)算法”。所述“期望最大化(EM)算法”可以在模型拟合期间用于迭代地搜索顾及一给定模型的最佳拟合的模型参数。所述模型包括两个主要步骤：最大化及期望。所述期望及最大化步骤可以执行一预定次数(比如两百次)，以达到最优解决方案。所述“期望最大化(EM)算法”可接受一具有初始化自回归(AR)参数的模型，并为所述模型输出估计自回归(AR)参数。

由于所有噪声项被假定为加性高斯(Gaussian)噪声，一个时窗长度N的完整数据对数似然性为：

我们希望最大化关于Θ＝(φ，Q，R)的对数L，但我们未能获知所述隐藏振荡x_t。我们运用所述“期望最大化(EM)算法”交替及迭代地估计(在期望步骤)及最大化(在最大化步骤)所述对数似然性。在迭代r时，在给定一集合Θ_r的情况下，我们使用所述“卡尔曼滤波器”(Kalman filter)来估计x_t，这使我们能够获知：

接着，我们推断Θ_r+1：

在给定所述振荡及所述模型参数的情况下，所述“卡尔曼滤波器”(Kalmanfilter)可用于估计所述隐藏振荡。首先，预测下一个时间点的状态，然后将所述状态与所述观测进行比较。接着，根据最近观测到的数据的更新估计可以生成。在已知所述完整观测时间序列的情况下，我们可运用后向平滑来优化所述更新，以列入所述完整观测时间序列(即：固定区间平滑)。

由于t，t₁，t₂＝1..N，我们注意到：

及

并运用后向平滑算法来计算那些数量如下：

预测：

“卡尔曼”(Kalman)增益：

更新：

以及一集合的后向递归。由于t＝N..1及t₁＜t₂：

反向增益：

平滑：

协方差：

现在描述所述“期望最大化(EM)算法”的期望步骤。以上描述的“卡尔曼滤波器”(Kalman filter)估计方法可包括在所述“期望最大化(EM)算法”的期望步骤中。这里我们描述用于以谐波分量进行状态空间振荡分解的G_r(Θ)。一单一振荡以一旋转矩阵

一缩放参数a及一过程噪声协方差矩阵Q＝σ²I_2×2定义。接下来，我们考虑d独立振荡，这些独立振荡分别与h₁，...，h_d谐波相关。由于j＝1..d，一基本频率为ω_j的振荡是h＝1...h_j、分别由

a_j，h及

定义的谐波的和。我们注意到振荡分量的总数为

由于

j＝1..d及h＝1..h_j，我们注意到与振荡j的

谐波相关的2乘2对角区块V_j，h。Φ及Q是区块对角矩阵，它们的对角区块是

及Q_j，h：

此外，我们将利用

而且由于

我们注意到：rt(U)：＝U₂₁-U₁₂，而且tr(U)＝U₁₁+U₂₂。取用一固定集合的参数Θ_r＝(Φ，Q，R)_r在迭代r时的对数似然性log L的条件期望，我们求得：

其中：

在所述最大化步骤时，G_r可以独立地对R及(Φ，Q)最大化。我们求得：

由于Q为区块对角，我们可写：

A为对称，而且Φ乃一区块对角矩阵(其元素为2×2旋转矩阵)，我们建立方程式(23)，并求得：

对过程噪声协方差

缩放参数a_j，h及基本频率ω_j微分，得到：

我们将(28)的两个上方的方程式代入下方(即第三个)方程式，而且我们注意到：

及

使用三角恒等式，方程式(28)最后可重写为：

总体上，由于j＝1..d，如果h_j＞1，我们利用方程式(30)来数值求解ω_j，并推导h＝1..h_j的a_j，h及

如果h_j＝1，我们立即得到：

参阅图11，其中显示一个以一状态空间模型来拟合脑电图(EEG)数据的示例性流程1100。所述流程1100可包括使用一先验分布来拟合alpha振荡，以及使用如以上所述的能够检测谐波的“期望最大化(EM)算法”。所述状态空间模型可包括几个振荡分量及这些分量中任何分量的谐波振荡，视所述脑电图(EEG)数据(即：一病人是否正在休息、清醒、无意识、等等)而定。所述丙泊酚示例性案例的状态空间模型可包括一慢速波、一alpha振荡以及所述的alpha振荡的谐波振荡。所述状态空间模型亦可包括附加的非谐波振荡，视所述脑电图(EEG)数据而定。例如，如果一病人已被施用丙泊酚，所述流程1100可拟合一慢速振荡分量及一具有或不具有谐波振荡的alpha振荡。所述流程1100可以实施为至少一个存储器上的指令。所述指令可接着由至少一个处理器执行。

在步骤1104，所述流程1100可接收脑电图(EEG)数据。所述脑电图(EEG)数据可接收自一病人监测装置，比如接收自以上连同图2描述的病人监测装置202。所述脑电图(EEG)数据可作为一波形解释。所述流程可接着进入步骤1108。

在步骤1108，所述流程可选择一起始模型。在一些实施例中，所述流程1100可通过拟合一高阶自回归(AR)模型、以最低AIC识别一目标阶次、然后识别所述目标阶次的最显著频率峰值以确定所述AR(2)分量的起始参数，得而选择一起始模型。所述丙泊酚示例性案例的起始参数可包括每一AR(2)分量(即振荡)的起始极点，所述起始极点限定所述初始慢速振荡峰值频率及所述初始alpha振荡峰值频率。所述起始极点包括可以用来确定一峰值频率的信息。在一些实施例中，所述流程1100可选择一具有最少个分量(比如只具有一个代表一慢速振荡的单一AR(2)分量)的模型。在后续的多个步骤中，所述流程1100可确定是否有更多分量，然后将附加分量添加到所述模型。所述AR(2)分量中的任何分量可能在模型参数上有一个或多个先验分布。所述先验分布可让所述模型能够更准确地拟合到一给定振荡。如以上所述，所述alpha振荡范围中的AR(2)分量可在所述振荡频率上有一“冯·米塞斯”(Von Mises)先验分布，以便更好地拟合某些alpha振荡，比如相对弱的alpha振荡。在概率值为0.367时，所述“冯·米塞斯”(Von Mises)先验分布的中心频率可能为10Hz，带宽可能为大约1.2～2Hz。如以上所述，一beta分布可用于代表用于所述模型的阻尼因数的先验分布。所述beta分布可通过具有一从0至1的域，使所述阻尼因数保持在0与1之间。所述模型可接着进入步骤1112。

在步骤1112，所述流程1100可将所述模型拟合至所述脑电图(EEG)数据。所述流程1100可通过固定所述噪声频率来隔离线路噪声，而线路噪声可能alpha频率范围之上发生。所述流程1100可使用一“期望最大化(EM)算法”来拟合所述模型。所述“期望最大化(EM)算法”可以是以上所述的【用于“独立及谐波振荡分解”的“期望最大化(EM)算法”】部份中描述的“期望最大化(EM)算法”。所述“期望最大化(EM)算法”可拟合任何数目的独立振荡，而其中每一独立振荡与一个或多个谐波有关。以这种方式，可以检测所述独立振荡(即：一慢速波及一alpha振荡)的谐波的存在，即使如果所述谐波附近有其他独立振荡。如以上所述，所述“期望最大化(EM)算法”的一期望步骤可迭代一期望步骤，而且所述“期望最大化(EM)算法”的一最大化步骤各可在估计参数(比如自回归AR参数及噪声协方差)输出之前迭代一预定次数，比如迭代200次。所述迭代亦可继续，直到满足某些收敛性判别准则为止(包括但不限于满足被估计的参数的收敛性、或对数似然密度或对数后验密度的值中的收敛性)。接着，最后的估计参数(比如振荡器频率)可由一从业人员或其他流程用于分析所述振荡器。在所述“期望最大化(EM)算法”已经迭代所述期望步骤及最大化步骤所预定的次数之后，所述流程1100可进入步骤1116。

在步骤1116，所述流程1100可确定是否存在更多振荡。在一些实施例中，所述流程1100可确定所述脑电图(EEG)数据是否有不能仅凭白噪声解释的残差或创新。所述流程1100可将所述预测估计时间序列(即：来自上述“卡尔曼滤波器”(Kalman filter)的

)从所述脑电图(EEG)数据的一时间序列中减去，以获得所述创新的一时间序列。所述流程1100可确定所述创新的频谱中是否存在显著的峰值。白噪声通常具有恒定的频谱。所述流程1100可计算所述创新的时间序列的累积功率级图，并确定所述累积功率级图是否在统计学上显著不同于一不断增加的功率图。如果没有找到所述创新与白噪声之间在统计学上的显著差异，则所述流程1100可确定所述脑电图(EEG)数据中没有更多振荡存在，以及应使用在步骤1116中拟合的模型(亦可称作一拟合模型)来近似化所述脑电图(EEG)数据的振荡。如果存在一种多多种在统计学上的显著差异，则所述流程1100可确定所述脑电图(EEG)数据中存在更多振荡。所述流程1100可确定所述创新的累积功率级是否像白噪声那样不断增加，以确定所述创新是否与白噪声一致。在一些实施例中，所述流程1100可简单地为所述模型计算一“赤池信息准则”(AIC)评分，并将所述“赤池信息准则”(AIC)评分与一先前拟合至所述数据(即：先前在步骤1116执行者)的模型的“赤池信息准则”(AIC)评分进行比较。在后续的多个步骤中，如果存在更多振荡，所述流程1100可添加更多振荡分量至所述模型，以便更好地拟合所述模型。如果模型的“赤池信息准则”(AIC)评分一直以附加振荡分量改善，则所述流程1100可继续添加振荡分量，直到“赤池信息准则”(AIC)评分恶化为止。如果当前模型的“赤池信息准则”(AIC)评分优于所述先前模型，所述流程1100可确定所述脑电图(EEG)数据中存在更多振荡。如果当前模型的“赤池信息准则”(AIC)评分劣于所述先前模型，所述流程1100可确定所述脑电图(EEG)数据中不存在更多振荡，并确定所述先前模型应作为所述拟合模型，用于估计所述脑电图(EEG)数据中的振荡。接着，所述流程可进入步骤1120。

在步骤1120，如果所述流程1100确定在步骤1116处不存在更多振荡(步骤1120处的答案为“否”)，所述流程可进入步骤1128。如果所述流程1100确定在步骤1116处存在更多振荡(步骤1120处的答案为“是”)，所述流程可进入步骤1124。

在步骤1124，所述流程1100可估计一附加AR(2)分量及相关参数。一高阶自回归(AR)模型可拟合至所述创新的所述时间序列，而且对促成所述自回归(AR)分量的最高峰值贡献最大的一对复极点可被选择。接着，可以根据该对复极点，确定起始参数。可选择地，可通过估计所述时间序列的功率及根据所估计的功率估计峰值频率，从所述残差的时间序列估计所述附加AR(2)分量的一中心频率。可选择地，可根据所述创新的功率谱，估计所述中心频率。所述附加AR(2)分量的一阻尼因数可初始化为一来自先前数据集合的平均数，或根据所述累积功率级图的形状识别。所述估计中心频率及阻尼因数可提供与该对复极点所提供的信息相同的信息。所述时间序列可在导出所述中心频率及/或阻尼因数之前，以一低通滤波器拟合。可选择地，所述流程可使用Halleret等人描述的“FOOOF算法”，根据所述残差估计所述中心频率及阻尼参数。接着，所述流程可进入步骤1112。

在步骤1128，所述流程1100可输出所述包括估计自回归(AR)参数及噪声参数的拟合模型。所述流程1100可将所述拟合模型输出到一显示器供一用户查看，及/或输出到一存储器供作存储及/或供另一流程使用。所述拟合模型可以是在步骤1120中确定不存在更多振荡之前、在步骤1116拟合的最后一个模型。如果所述流程1100迭代更多分量至所述模型直到“赤池信息准则”(AIC)评分恶化为止，则所述拟合模型可以是一具有最佳“赤池信息准则”(AIC)评分的模型。所述流程1100可根据所述拟合模型产生一报告，并输出所述报告，以使所述拟合模型更容易被人看到。所述流程1100可接着结束。

以上所述的状态空间模型及相关的拟合方法可用于对脑电图(EEG)数据执行跨频分析。跨频分析包括任何对不同振荡之间的关系(例如一第一振荡的相位或振幅与一较高频率的振荡的相位或振幅之间的关系)的分析。相位振幅耦合(PAC)是一常见例子，其中一第一振荡的相位调制一第二振荡的振幅。相位振幅耦合(PAC)可在一频率范围内执行，或在一振荡的相位与一宽带(非振荡)信号的振幅之间执行。我们也提供一新方法来估计一第一振荡(如以上所述，表达为一解析信号)的实部份及虚部份与一第二振荡的振幅、频率范围或宽带信号之间的关系。

用于相位振幅耦合(PAC)的标准方法使用分箱直方图来量化相位与振幅之间的关系，而这是统计效率低下的主要原因。原始信号y_t首先经带通滤波，以隔离慢速振荡及快速振荡。接着，实施一希尔伯特变换(Hilbert transform)，以估计所述慢速振荡的瞬时相位

以及所述快速振荡的瞬时振幅

在时间t处，根据所述慢速振荡相位的瞬时值

将所述alpha振幅

分配予长度为δψ的等距相位分箱(通常18个)的其中之一。所述直方图是在观测的某个时窗T内(比如在a～2分钟时期内)构建的，这产生以下相位振幅调制图(PAM)：

对于一给定时窗T而言，相位振幅调制图(PAM)(T，.)是一概率分布函数，其评估所述快速振荡振幅如何对所述慢速振荡相位分布。接着，通常以均匀分布的KL散度(Kullback-Leibler divergence)来测量所述调制的强度。这产生“调制指数”(MI)：

最后，在这个标准方法下，使用替代数据(比如随机排列)来评估所观测的“调制指数”(MI)的统计显著性。根据一均匀分布(其间隔取决于问题动力学)，绘制随机时移Δt，并使用所述迁移快速振幅

及所述原始慢速相位

来估计相位振幅耦合。接着，为这个排列时间序列计算“调制指数”(MI)，然后重复所述过程，为所述“调制指数”(MI)构建零分布。如果所述原始“调制指数”(MI)大于所述排列值的95％，则该“调制指数”(MI)被视为显著。总体上，这个方法需要所述基本流程在足够长的时间内保持静止，以期能够足够合理地评估所述调制图，而且使得能够排列数量足够的可比数据片段，以评估显著性。

相反地，我们引用跨频分析的一个参数表达式，该参数表达式能够以数量较少的参数量化所述跨频耦合关系。为了这么做，我们考虑一形式为

的约束线性回归问题，该问题最后可重写为：

K^mod控制所述调制的强度，而φ^mod是优选相位，如图12H-J所示，在该优选相位附近，所述快速振荡

的振幅为最大。图12A显示原始脑电图(EEG)数据。图12B显示所述脑电图(EEG)数据的一个六秒时窗。图13C显示在所述六秒时窗中拟合的一慢速振荡。图12D显示在所述六秒时窗中拟合的一快速振荡。图12E显示所述状态空间模型的一示例性形式。图12F显示一慢速振荡相位。图12G显示一快速振荡振幅。图12H显示一估计K^mod及其分布。图12I显示一估计φ^mod及其分布。图12J显示已回归的alpha波振幅。例如，如果K^mod＝1及φ^mod＝0，在所述慢速振荡的峰值处，所述快速振荡最强。另一方面，如果φ^mod＝π，则在所述慢速振荡的低谷或最低点处，所述快速振荡最强。

最后，如图12F-J所示，非但不依赖替代数据来确定统计显著性(那么做会进一步降低效率)，我们的模型表达式允许我们估计所述调制参数的后验分布p(K^mod，φ^mod|{y_t}_t)，并推断所述相关可信区间(CI)。

我们将我们的方法称为“状态空间跨频分析”(SSCFA)方法。由于生理***随时间变化，我们应用这个方法于多个非重叠时窗。在我们的方法的变化中，我们亦可使用一第二状态空间模型，对所述跨时窗的调制参数施加一时间连续性约束，得到我们称为“双状态空间跨频分析”(dSSCFA)估计。

为了证明我们的方法的性能，如图14中所示，我们首先分析来自一接受丙泊酚以诱导镇静及无意识的志愿者的脑电图(EEG)数据。图14A显示供应予所述病人脑电图(EEG)数据的丙泊酚浓度。图14B显示一对应于所述病人脑电图(EEG)数据的响应概率。图14C显示所述病人脑电图(EEG)数据的调制回归值。图14D显示所述病人脑电图(EEG)数据的参数谱。图14E显示所述病人脑电图(EEG)数据的“相位振幅调制图”(PAM)。图14F显示所述病人脑电图(EEG)数据的调制参数。正如预期那样，随着丙泊酚的浓度增高，所述对象对听觉刺激的响应概率减低。如图14D中所示，在此时间期间，所述参数功率谱密度(见方程式(59))改变，在所述响应概率开始减低时产生beta(12.5～25Hz)振荡，紧随其后，在所述响应概率为零时产生慢速(0.1～1Hz)振荡及alpha(8～12Hz)振荡。对于一窗T，如图14F中所示，我们以“双状态空间跨频分析”(dSSCFA)估计所述调制强度K^mod及

(以及可信区间)，而且我们收集所述“相位振幅调制图”(PAM)中的那些估计：图14E的PAM(T，)。对于一给定窗T而言，PAM(T，)是一具有支持[-π，π]的“概率密度函数”(pdf))。该函数评估所述快速振荡的振幅如何对所述慢速振荡的相位分布。当所述响应概率为零时，我们观察到一强“最大峰值”(K^mod≈0.8，

)图案，其中在所述慢速振荡的峰值处，所述快速振荡振幅为最大。在过渡至无意识及从无意识过渡期间，我们观察到一较弱强度的“最大低谷”图案(K^mod≈0.25，

)，其中在所述慢速振荡的低谷处，所述快速振荡振幅为最大。注意由于在耦合强时，我们的模型更好地预测

所述调制关系的决定系数R²反映所述耦合的强度K^mod。

当在冗长、连续且固定的时窗上平均时，传统方法提供对相位振幅耦合(PAC)的良好定性评估。然而，在许多情况下，如果实验条件或临床情况迅速变化，可能有必要在较短的时窗内进行分析。在以前的工作中，已经以相对长的时窗(δt＝120s)，使用传统方法分析相位振幅耦合(PAC)，此举在该情况下合适，这是由于丙泊酚在超过～14分钟的间隔内以固定速率给药。所述“状态空间跨频分析”(SSCFA)及“双状态空间跨频分析”(dSSCFA)方法的增高统计效率使得可能分析短得多的时窗(δt＝6s)；我们以两个对象来说明，其一为强耦合(如图15中所示)，而令一为弱耦合(如图16中所示)。图15A为一响应概率图。图15B为一丙泊酚浓度图。图15C为一使用六秒时窗的标准方法的调制图。图15D为一与图15C相关的调制指数图。图15E为一使用一百二十秒时窗的标准方法的调制图。图15F为一与图15E相关的调制指数图。图15G为一用于“状态空间跨频分析”(SSCFA)方法的调制图。图15H为一与图15G相关的调制指数图。图15I为一用于“双状态空间跨频分析”(dSSCFA)方法的调制图。图15J为一与图15I相关的调制指数图。图16A为一响应概率图。图16B为一丙泊酚浓度图。图16C为一使用六秒时窗的标准方法的调制图。图16D为一与图16C相关的调制指数图。图16E为一使用一百二十秒时窗的标准方法的调制图。图16F为一与图16E相关的调制指数图。图16G为一用于“状态空间跨频分析”(SSCFA)方法的调制图。图16H为一与图16G相关的调制指数图。图16I为一用于“双状态空间跨频分析”(dSSCFA)方法的调制图。图16J为一与图16I相关的调制指数图。

为了这么做，我们根据所述调制图及所述“调制指数”(MI)评估值，以δt＝120s或δt＝6s，对“状态空间跨频分析”(SSCFA)方法、“双状态空间跨频分析”(dSSCFA)方法以及所使用的标准方法进行比较。后者通过测量(对于任何窗T)PAM(T,)与所述平均分布的差异程度，评估所述调制的强度。KL散度(Kullback-Leibler divergence)典型地用于此目的。因此，所述估计“相位振幅调制图”(PAM)中的任何随机波动将使“调制指数”(MI)增高，引起偏差。我们的模型参数化用于推导“相位振幅调制图”(PAM)、“调制指数”(MI)及相关的可信区间(CI)，但标准的非参数分析一般依赖分箱直方图。结果是，它们通过构建替代数据集及报告概率值(p-values)来估计统计显著性。

两个对象都在他们过渡至无意识及从无意识过渡时，展示先前描述的典型相位振幅互相关性轮廓。然而，由于“状态空间跨频分析”(SSCFA)更有效率地估计相位及振幅，并且生成平滑的“相位振幅调制图”(PAM)估计值(即使是在短窗上)，推导自“状态空间跨频分析”(SSCFA)的“调制指数”(MI)估计值显示的偏差少于所述标准方法。出于相同原因，φ^mod估计值显示的方差少于所述标准方法。所述“状态空间跨频分析”(SSCFA)算法在“相位振幅调制图”(PAM)上提供一时间连续性约束，使其可能跟踪相位振幅耦合(PAC)中的随时变化，并进一步减低“相位振幅调制图”(PAM)估计值的方差。最后，我们的参数策略为K^mod、φ^mod及“调制指数”(MI)提供后验分布，使其可能为每一变量估计可信区间(CI)，而且我们的参数策略可在不需诉诸替代数据方法的情况下评估显著性。

为了说明我们的方法在一代表动物模型中的侵入性记录的不同情境中的性能，在一假定涉及theta(6～10H)及低gamma(25～60Hz)振荡的学习任务期间，我们对小鼠数据进行了分析。我们对2秒时窗运用“双状态空间跨频分析”(dSSCFA)，并确定海马体的CA3区域中的theta-gamma耦合随着小鼠学会辨别任务而增强。在我们使用“双状态空间跨频分析”(dSSCFA)方法的分析中，我们并未预先选择感兴趣的频率，而且未指定带通滤波截止频率。情况却是，所述“期望最大化(EM)算法”能够在theta及gamma范围中的初始起点给定的情况下，根据所述数据，为相位及振幅估计特定的基本振荡频率。因此，我们说明了：我们的方法可有效地运用来分析局部场电位(LFP)数据，而且可在不需指定固定频率或频率范围的情况下识别基本振荡结构。

为了以更***的方式对我们的、作为不同调制特征及信噪比水平的函数的算法，我们分析了多个模拟数据集。这些模拟数据是有意地使用与所述状态空间振荡器模型不同的生成过程或模型(即：所述模拟数据生成过程在我们的方法中使用的“模型类”之外)来创建。这里，我们着重于一慢速分量及一alpha分量，以重现我们的主要实验数据案例。我们这么做的目的并未为我们的算法的精度及准确性提供详尽的表征，这是由于其强烈取决于信噪比、信号形状等。相反地，我们是为了说明，在短及噪声时变数据集的情况下，我们的算法如何及为何超过标准分析。

我们首先就具有在多个时间尺度上变化的调制参数的宽带信号，对“双状态空间跨频分析”(dSSCFA)的分辨率及鲁棒性与常规技术进行比较。而且报告了不同生成参数及相关信号的结果(详见下文，

σ_s及σ_f)。虽然标准技术在具有固定耦合参数的窗上平均时鲁棒，但当所述调制参数在不同的窗中迅速变化时，标准技术变得无效。当使用长窗时，所述调制无法解析。然而，如果我们通过减小窗大小来补偿，所述估计值的方差显著地增大。必须凭经验找到一折中方式。另一方面，我们看到，当运用于6秒窗时，即使是在信噪比低的情况下，“双状态空间跨频分析”(dSSCFA)可跟踪振幅调制中的迅速变化。所述“双状态空间跨频分析”(dSSCFA)算法也所述调制参数的后验分布的估计值，使得能够直截了当地构建可信区间(CI)及执行统计推断。相比之下，所述替代数据方法变得不可行，这是由于越来越少数据片段可供进行无序排列。

在一项研究中，一非线性相位振幅耦合(PAC)公式被描述为一驱动自回归(DAR)过程，其中调制信号为慢速振荡的多项式函数。被称为驱动程序的后者在靠近一预设频率处从所述观测被过滤掉，并用于估计驱动自回归(DAR)系数。所述信号参数谱密度随后被推导为所述慢速振荡的函数。所述调制接着根据相位表达，所述快速振荡的功率在靠近所述相位处优先分布。对所述驱动程序进行一网格搜索，在一快速频率范围内，产生每一慢速中心频率的调制图。接着，与最高似然性及/或最强耦合关系相关的频率被选择为最终耦合估计值。

这个参数表示改进效率，尤其是在短信号窗的情况下，但是由于它依赖一初始过滤步骤，而且它也分享常规技术的局限性。如我们所见，虚假跨频耦合(CFC)可从突变信号或非线性出现。此外，这个初始过滤步骤可能污染来自具有宽带慢速振荡的短数据片段的相位振幅耦合(PAC)估计值。

为了将我们的方法与标准技术及所述驱动自回归(DAR)方法比较，我们以不同的感兴趣频率(f_s及f_f)、谱宽

及信噪比(SNR)生成调制信号(方程式(88)，λ＝3，及φ^mod＝-π/3)。用于那些产生参数的典型信号迹线，可以使用方程式(87)来生成。接着，我们对这些方法恢复所述慢速振荡及所述快速振荡的优劣程度、哪个或所述调制相位进行比较。与这里介绍的其他方法相反，“状态空间跨频分析”(SSCFA)并未计算完整的共调制图以选择感兴趣频率，而是通过拟合所述状态空间振荡器模型来识别感兴趣频率。耦合只是在一第二步骤中估计。虽然我们在前面段落中使用有形的先验知识来初始化所述算法，但我们调整一初始化程序，以提供公平的比较。对每种情况，我们生成了400个六秒窗。在必要时，使用

来提取所述驱动程序。

我们发现，在每一案例中，我们的算法更好地检索快速频率，尤其是在所述慢速振荡为宽带时。在估计调制相位时，我们的算法也超过其他方法：在宽带

或弱[(σ_s，σ_s)＝(0.7，0.3)]慢速振荡的情况下，我们的算法稳定；而且在几乎所有考虑的情况下，我们的算法准确地估计φ^mod，异常值很少，标准差较小。

尽管虚假跨频耦合(CFC)在协调神经***中可能扮演的中心角色，虚假跨频耦合(CFC)分析的标准方法受限于许多需注意事项，而这些需注意事项是一持续担忧根源。当基本信号急剧转变或具有非线性时，可能引起虚假耦合。另一方面，如果用于带通滤波的频带选择不适当，真实的基本耦合可能会丢失。这里我们说明，对所有这些局限性，我们的“状态空间跨频分析”(SSCFA)方法如何鲁棒。我们也显示，我们的方法如何能够使用一线性模型-反直觉地提取一信号的非线性特征。

振荡神经波形可能具有不限于窄带的特征，比如急剧变化或不对称性。在这样的情况下，以标准带通滤波器来截短它们的频谱内容可能会扭曲所述信号的形状，而且可能引进人工制造的分量，这些人工制造的分量可能被不正确地解释为耦合。

所述状态空间振荡器模型提供一替代带通滤波的、可接纳非正弦曲线波形的选择。在本段落中，我们扩展所述模型以明确地表示所述慢速振荡信号的谐波，因而允许所述模型更好地表达具有急剧转变的振荡以及那些可能由非线性***产生的振荡。如以下进一步地描述那样，为了这么做，我们对相同的基本频率f_s，优化h振荡。我们进一步将这个模型结合信息准则(“赤池信息准则”(AIC)或“贝叶斯信息准则”(BIC))，以确定：(i)慢速谐波h的数目；以及(ii)一快速振荡的存在或不存在。我们通过最小化ΔIC＝IC-min(IC)来选择最佳模型。虽然“赤池信息准则”(AIC)及“贝叶斯信息准则”(BIC)表现相似，但这里我们只报告基于“赤池信息准则”(AIC)的相位振幅耦合(PAC)估计。当多个慢速谐波受青睐时，我们使用所述基本振荡的相位来估计相位振幅耦合(PAC)。

我们使用一范德波尔(Van der Pol)振荡器来模拟一非对称急剧变化信号(方程式(89)，∈＝5，ω＝5)，我们将观测噪声

添加到所述振荡器。接着，我们考虑了两种情景，一种是带有调制快速正选曲线波(图27A，

及f_f＝10Hz)，而令一种不带调制快速正选曲线波(图27B)。由于我们的模型能够拟合所述急剧转换，“赤池信息准则”(AIC)及“贝叶斯信息准则”(BIC)都识别正确数目的独立分量。因此，当未检测到清晰的快速振荡时，不计算任何相位振幅耦合(PAC)。另一方面，当不存在任何快速振荡时，标准技术基于所述急剧变化的慢速振荡提取一快速分量，导致检测到虚假耦合。

在虚假跨频耦合(CFC)分析中，产生自信号传导谐波的非线性输入是一类似障碍。如果我们考虑一慢速振荡

非线性传导为

我们可写以下二级近似公式：

如果

带通滤波y_t以接近0.9-3.1Hz的频率提取一峰值为f_f＝2Hz的振荡，将产生：

在这样的情况下，标准的虚假跨频耦合(CFC)分析推断显著的耦合，而振荡分解在没有虚假跨频耦合(CFC)的情况下正确地识别一谐波分解。

我们观察到，所述(扩展的)状态空间振荡器比窄带分量更适合建模生理信号。此外，与所述信号内容的先验知识结合的模型选择范例(例如丙泊酚麻醉慢速alpha振荡或啮齿动物海马体theta-gamma振荡)让我们能够以更有原则的方式研究跨频分析。

如果带宽过窄的带通滤波器应用于一调制信号，所述调制结构可能被摧毁。我们的“状态空间跨频分析”(SSCFA)算法使用一状态空间振荡器分解，所述分解并未明确地建模引起所述调制旁瓣的结构关系。在测试期间，我们发现所述调制成功提取，如在拟合时间序列中及所述频谱中所观测的那样。图17A显示使用一标准方法对振荡分量的分量进行分解，而图17D显示为一给定信号使用一“状态空间跨频分析”(SSCFA)方法对振荡分量的分量进行分解。图17B显示使用所述标准方法确定功率谱密度，而图17E显示使用所述“状态空间跨频分析”(SSCFA)方法确定功率谱密度。图17C显示一通过标准方法进行的相位振幅调制估计，而图17F显示一通过所述“状态空间跨频分析”(SSCFA)方法为所给定的脑电图(EEG)信号进行的相位振幅调制估计。所述“状态空间跨频分析”(SSCFA)方法所使用的振荡分解方法与所述“双状态空间跨频分析”(dSSCFA)方法所使用的振荡分解方法相同。通过使所述快速分量的频率响应的宽度足以包围所述旁瓣，所述模型能够实现这个目的。所述算法通过扩大所述噪声协方差

以及缩小a_f，实现这个目的。可能适用一较高阶次的模型(比如像ARMA(4，2)的模型，这种模型将表达2种振荡的乘积，而且其极点在

之中)，或通过一非线性观测，直接建模耦合。然而，在这两种情况下，我们发现：这样的模型难以拟合至所述数据，而且当运用于嘈杂、非平稳、非正弦曲线的生理信号时，这样的模型迅速变得不受约束。相反地，我们发现我们更简单的模型能够鲁棒地提取所述调制高频分量。

虽然所述“期望最大化(EM)”算法确保收敛，但需最大化的对数似然性并非总是凹函数。一由d振荡组成的信号可以以p∈[|d，2d|]阶的最佳自回归(AR)过程的参数初始化。然而，由于所述电生理信号的非周期性分量，这样的程序可能偏移所述初始化。确实地，所述非周期性分量通常以一1/f^χ幂律函数描述，该函数可被所述自回归(AR)过程回归。在这样的情况下，所述初始化可能无法顾及一实际基本振荡。

为了帮助缓解这个潜在问题，我们使Halleret等人的算法适应所述状态空间振荡框架。我们的初始化算法的目的是在以所述振荡的参数功率谱密度(PSD)(方程式(ZZHH))拟合产生的频谱之前，从所述非周期性分量解开所述振荡分量，其将在下文进一步描述。

所述观测数据信号y_t的参数功率谱密度(PSD)可以使用所述多锥方法估计。接着设定频率分辨率r_f(一般设定为1Hz)，因而产生时间带宽乘积

然后选择锥数目K，使得

所述观测噪声R₀(用于初始化R)可以以下列方程式估计：

接着，可以从所述参数功率谱密度(PSD)移除所述观测噪声R₀。

所述非振荡分量可从所述参数功率谱密度(PSD)回归。接着，以dB为单位的所述非周期性信号的参数功率谱密度(PSD)在频率f处以下列方程式建模：

χ控制所述非周期性信号的斜率、偏移g₀以及截止频率f₀。应用一第一带通拟合，以识别对应于非振荡分量的频率：如图13A中所示，只拟合f₀，而χ及g₀依次设定为χ＝2及g₀＝PSD(f＝0)。如图13B中所示，固定一阈值(一般为残差的0.8分位数)，以识别与所述非周期性信号相关的频率。接着，如图13C中所示，只对那些我们据以推断f₀、χ及g₀的频率应用一第二带通拟合。如图13D中所示，我们从以dB为单位的原始参数功率谱密度(PSD)去除g(f)，以生成一修正参数功率谱密度(PSD)用于所述算法的第二步骤。所拟合的参数可接着用于初始化所述“期望最大化(EM)算法”。

根据所述修正参数功率谱密度(PSD)，我们使用方程式(ZZHH)中给予的理论参数功率谱密度(PSD)来拟合一给定数目d₀(例如d₀＝4)的独立振荡。为了这么做，在将一振荡理论频谱拟合在接近需识别的峰值的宽度邻域2r_f处之前，我们识别具有足够宽度(宽于r_f/2)的参数功率谱密度(PSD)峰值。对振荡j，我们推定(f_j)₀、(a_j)₀及

由于

代表在移除所述非周期性分量之后一给定振荡的偏移，我们调整它以估计

接着，如图13E中的DSS初始化拟合线1300所示，将产生的频谱PSD_j减去，然后重复所述过程，直到对所有的振荡的估计完成为止。最后，我们估计在所述邻域(f_j)₀中的振荡j的功率P_j，然后根据下式，估计其对总功率P₀的贡献：

将振荡排序，然后使用产生的参数，以所述d∈[|1，d₀|]第一振荡初始化所述所述“期望最大化(EM)算法”。

为了改进统计效率，我们引入相位振幅耦合(PAC)的参数表达式。对一给定窗，我们考虑以下(约束)线性回归问题：

其中：β＝[β₀ β₁ β₂]^T，

及

如果我们定义：

φ^mod＝tan^-1(β₂/β₁)以及A₀＝β₀ (42)

我们看到方程式(41)相当于：

设定

确保所述模型一致，即：所述调制包络不会超过所述载波信号的振幅。但是实施这个约束，在计算上颇为昂贵。如果所述数据的信噪比高而使得K^mod不太可能偶然大于1，我们亦可选择解决所述无约束问题

在约束解之下，β的后验分布为一截短多变量t-分布：

所述似然性、共轭先验、后验参数β，V，b，v，以及所述归一化常数在以下进一步说明及推导。我们称此估计值为状态空间跨频分析(SSCFA)，而且我们注意到：

所述标准方法依赖替代数据来确定统计显著性，这进一步降低其效率。相反地，我们估计所述后验分布p(β|{y_t}_t)，根据这个估计，我们求得所述调制参数K^mod及φ^mod的可信区间(CI)。为了估计所述后验分布，我们从以下模型给予的后验分布采样：(i)所述状态空间振荡器模型；及(ii)以上所述的参数相位振幅耦合(PAC)模型。

(i)用于以上所述的“期望最大化(EM)算法”的r^th“期望步骤”的“卡尔曼滤波器”(Kalman filter)提供下列瞬间(t，t′＝1..N)：

因此，我们采用以下返程时采样时间l₁序列：

其使用：

其中P为一4N x4N矩阵，其块条目由右式给出：

(ii)对于每个

我们使用方程式(8)来计算重复采样的慢速振荡的相位φ及快速振荡的振幅

我们接着使用方程式(44)，从

采l₂个样本。结果我们生成l₁×l₂个样本供估计：

我们最后使用一L₂范数来构建相关可信区间(CI)，并继而导出K^mod及φ^mod的可信区间(CI)。

我们将所述时间序列分段为多个长度为N的非重叠窗，并对这些非重叠窗应用先前描述的分析。因此，我们生成{β^SSCFA}_T，这是说明所述调制的、以

为单位的一集合向量，其中T代表一长度为N的时窗，或其T代表一长度为N的时窗的调制。

一第二状态空间模型可用于表达所述调制动力学。这里我们跨时窗将一带有观测噪声的p阶自回归(AR)模型拟合至所述调制向量β^SSCFA。这产生所述“双状态空间跨频分析”(dSSCFA)估计值：

如以下所述，我们通过在数字上求解及优化“尤尔-沃克”(Yule-Walker)方程式行进，而且我们以“贝叶斯信息准则”(BIC)选择所述阶次p。最后，在有必要时，我们可使用所述拟合参数来过滤所述l₁×l₂重复采样参数，以便为{β^SSCFA}_T构建一可信区间(CI)。

为了更好地以所述“状态空间跨频分析”(SSCFA)比较标准技术，在我们用于一窗T的参数模型之下，我们导出所述“相位振幅调制图”(PAM)的一个近似表达式：

我们对在所述麻醉剂丙泊酚给药期间发生意识丧失及恢复期间的人类脑电图(EEG)数据进行了分析。简要地，10名健康志愿者(18-36岁)被注入剂量越来越多、跨6个目标效应位浓度(0、1、2、3、4及5μgL^-1)的丙泊酚。输液由电脑控制，而且每一浓度维持14分钟。为了在行为上监测意识丧失及恢复，监测对象每4秒被给予一音频刺激(一咔嗒声或口头命令)，而且必须通过按压按钮进行响应。响应概率及相关的95％可信区间(CI)以蒙特·卡罗方法(Monte Carlo methods)来估计，以将一状态空间模型拟合至这些数据。最后，脑电图(EEG)数据以一抗混叠滤波器预处理，然后降采样至250Hz。

我们使用与振荡分解相关的参数表达式来计算所述脑电图(EEG)的频谱图。用于相位振幅耦合(PAC)分析的标准技术被应用在6秒时窗及120秒时窗，其alpha分量及慢速分量假定为已知，并以频率为0.1～1Hz及9～12Hz的带通滤波器提取。所述标准相位振幅耦合(PAC)分析方法的显著性以200个随机排列评估。

在下文，我们通过以相同的频谱内容建立一自回归移动平均(ARMA)过程，导出一振荡

的功率谱密度(PSD)的参数表达式。为了方便起见，我们接下来将删除指数j。首先，我们注意到一振荡为渐近二阶平稳。让我们计算其自協方差序列。由于

为一旋转矩阵，

且

因此，从方程式(2)，可得：

以及：

因此，我们可写，因t＝1..N，k＝0..N-t，

结果是，可以通过一个二阶平稳过程来近似化一振荡；而且，由于维纳-辛钦定理(Wiener-Khinchintheorem)，其理论功率谱密度为：

我们现在考虑所述ARMA(2，1)模型：

我们对其施加，t＝1...N，k＝0..N-t：

由此得出：

s_k＝φ₁s_k-1+φ₂s_k-2，k＞2

取：

φ₁＝2a cos(ω)

φ₂＝-a² (56)

满足方程式(55)的第一等式。其余条件可接着重写为：

从上式，我们导出2个负根

最后导致相同的自协方差序列。总体上，我们选择：

应用离散傅里叶变换(discrete Fourier transform)于方程式(54)，得到：

由于

是高斯噪声，

因此，我们的ARMA（2，1)PSD是：

最后，一居中于f₀的振荡的功率谱密度(PSD)为：

我们使用

而且我们假设方程式(41)中定义的模型的似然性为：

其中

而

共轭先验为：

我们选择先验超参数

及

以传达尽可能少的关于所述调制的相位及强度的信息。将方程式(61)的τ_β边缘化，产生一截短多变量t-分布：

确保定义所述多变量t方差。该多变量t方差为

接着，我们考虑独立随机变量

及

使得：

我们注意到

使用

然后，我们定义

以及

使得：

以及

此外，我们注意到，如果

而且由于

(其中

代表跨试验或窗的平均数，而<.>_t是跨一给定窗的时间平均值)，

及

因此，可以合理地使用

而且c＝1。总体上：

及

后验参数由以下方程式给定：

其中：

并且：

我们推断所述后验分布：

归一化参数

通过积分求得，为：

其中：

使用参数v、β及b^-1V的多变量t-分布计算。

最后，我们通过将方程式(68)的τ_β边缘化，推导方程式(44)，然后得到：

请注意：对于大样本，使用下列方程式可能有用：

我们现在为一给定自回归(AR)阶次p估计方程式(49)中定义的参数R_β、Q_β及

设

为以方程式(45)估计的调试向量的自协方差序列。我们得到：

其中δ_i，j为克罗内克函数(Kronecker delta)。方程式(72)可重写为：

对于一候选观测噪声

如果我们可反演方程式(73)，我们立即获得

及

使用“卡尔曼滤波器”(Kalman filter)，我们因此推论所述候选模型的似然性。

因此，我们注意到托普利兹矩阵(Toeplitz matrix)C＝(C_|i-j|)_i，j＝0..p的最小本征值

然后对

中的

在数字上优化所述模型似然性，其中我们获知

为满秩。

根据R_β，我们获得Q_β及

然后再次使用“卡尔曼滤波器”(Kalman filter)来估计

我们现在为一居中于τ、长度为δt＝N/F_s的窗推导一近似化参数调制图。我们将使用：

为了清楚起见，我们将无区别地使用t或k，而且我们提醒：

另外，我们假设：

·因k∈Ω_τ，

其中

及

·而且，为了简单起见，对于ψ∈[-π，π]，因所有h：

平滑，

根据中心极限定理，

及

我们因此可得：

但

因此：

另一方面：

但

及l∝N，所以：

使用如以上详述的相同论点，我们得到：

此外：

其中，γ为一函数，使得

由于：

因

“足够小”，我们可写：

最后：

模拟

我们在产生自不同模型的模拟数据集上测试了我们的“状态空间跨频分析”(SSCFA)方法及“双状态空间跨频分析”(dSSCFA)方法。我们通过结合单位方差高斯噪声、一居中于f_s＝1Hz(除非另外说明)的慢速振荡以及一居中于f_f＝10Hz(除非另外说明)的快速振荡，构建了每一个模拟信号。值得注意的是，我们选择使用一与我们用于分析所述数据的状态空间振荡器模型不同的方法或“模型类”，以产生这些模拟信号。对于标准处理，显著性以200个随机排列评估，f_s及f_f被假定为已知，而分量则以通频带设为0.1～1Hz(慢速分量)及8～12Hz(快速分量)的带通滤波器提取。

模拟所述慢速振荡

神经振荡不是完美正弦曲线，而却具有宽带特性。我们通过卷积(过滤)具有以下脉冲响应函数的独立同分布高斯噪声，模拟一宽带慢速振荡：

c(t)＝c₀(t)cos(ω_st) (86)

其中ω_s＝2πf_s，c₀为一

阶的布莱克曼窗(Blackman window)。所述慢速频率带宽

越小，所述信号越接近一正选曲线。当必要时，我们另外使用一π/2相移滤波器

来建模一解析慢速振荡

根据此解析慢速振荡，我们推断所述相位

产生的序列最后归一化，使其标准差设为σ_s。

为了评估我们的方法及所述标准方法的时间分辨率，我们生成了具有不同时变调制速率的模拟数据集。首先，如以上所述，为了构建所述模拟快速振荡，我们构建了一居中于ω_f＝2πf_f并归一化至σ_f的快速振荡，并通过以下方程式模拟所述快速振荡：

这里，

及

随时间变化，而且可遵守依次如图18A及图18B所示的动力学。模拟数据亦可使用以下替代调制函数产生：

其中

突然或急剧转变的信号可导致人工制造的相位振幅互相关性或相位振幅调制。为了评估我们的状态空间相位振幅耦合(PAC)方法在这样的情况下的鲁棒性，我们使用一范德波尔(Van der Pol)振荡器来产生一带突然变化的信号。这里，所述振荡x受以下微分方程式管辖：

使用一带固定时间步长的欧拉方法(Euler method)解方程式(89)。

综上所述，我们的算法的第一阶段，可在于第二阶段中以一自回归模型拟合非线性之前，基于所述调制提取非线性。这个方式的主要后果是扩大所述分量估计中的方差。参看图12G中所述快速振荡估计中的宽大可信区间(CI)的示例。我们继而从一比实际情况更宽的分布重复采样所述快速振荡。虽然这不影响φ^mod的估计值，但其在重复采样K^mod时产生一保守估计，即：所述可信区间宽于它们原来的可能宽度。尽管如此，我们发现我们的方法还是超过先前的方法。

我们已呈献一将振荡的状态空间模型与相位振幅耦合(PAC)的参数公式集成的新方法。在这个状态空间模型之下，我们将每一振荡表示为一解析信号，以直接地估计所述相位或振幅。我们接着使用一具有约束线性回归的参数模型来表征所述相位振幅耦合(PAC)关系。回归系数有效地表示只是几个参数中的耦合关系，其可以并入一个第二状态空间模型，以跟踪所述相位振幅耦合(PAC)中的时变。我们通过从多个应用分析神经时间序列数据，证明了这个方法的效率，并且使用基于不同生成模型的模拟研究，将这个方法的改进统计效率与标准技术进行了比较。最后，我们证明，我们的方法比与标准跨频耦合分析方法相关的许多局限鲁棒。

基于许多因素，我们的方法有效率。首先，所述状态空间解析信号模型提供对被分析的振荡的相位及振幅的直接了解。这个线性模型亦具有卓越能力通过施加根据所述数据估计的“软”频带限制，提取一非线性特征(所述调制)。所述振荡分解因此非常适合分析不受严格频带限制的生理信号。我们也提议了谐波延伸，谐波延伸可表达非线性振荡，使得可能更好地区分真实相位振幅耦合(PAC)与因自带通滤波人工产物而产生的虚假相位振幅耦合(PAC)。所述耦合关系的参数表达式可容纳不同调制形状，并甚至进一步提高所述模型的效率。

总体上，我们解决了与当前用于某些类别的跨频分析(比如相位振幅耦合(PAC)分析)相关的显著限制中的大部份。所述神经时间序列更有效地得到处理，感兴趣频带得以自动选择和提取，并建模。与标准方法相反，我们不需跨时平均化相位振幅耦合(PAC)相关数量，这使所需要的连续时间序列数据的量减少。此外，所述感兴趣信号的后验分布在我们提议的模型之下得以很好地定义。采样自它们绕过对建立替代数据的需求，而建立替代数据会蒙蔽所述数据中的非平稳结构以及低估假阳率。相反地，由于“状态空间跨频分析”(SSCFA)估计调制参数的后验分布，我们报告可信区间(CI)，并提供关于我们的结果的统计显著性以及所述调制的强度和方向的信息。我们对相位振幅耦合(PAC)的动态估计因此使得可能基于短得多的时窗进行推断-短至6秒(对0.1～1Hz的慢速信号而言)。已经提议其他新模型用于表达相位振幅耦合(PAC)，包括驱动自回归(DAR)模型及广义线性模型(GLM)。正如我们之前看到的那样，“状态空间跨频分析”(SSCFA)表现优于所述驱动自回归(DAR)及标准方法，尤其是在信噪比低时。与我们一样，所述广义线性模型(GLM)参数地表达所述调制关系，但以更一般的形式而言，该模型在统计上可能效率较低，而且其使用靴圈来提供置信区间。驱动自回归(DAR)及广义线性模型(GLM)方法都保持对用于信号提取的传统带通滤波的依赖，并因此容易受到这些滤波器引致的关键问题的影响。我们的方法是第一个使用状态空间模型与所述调制的一参数模型结合的方法。

这样的改进可显著地改进对未来涉及跨频耦合(CFC)的研究的分析，而且可使能够进行需要接近实时跟踪跨频耦合(CFC)的医疗应用。一种这样的医疗应用可能是基于脑电图(EEG)的对麻醉剂诱导无意识的监测。在丙泊酚诱导麻醉期间，频带不仅只是很好地定义，但所述相位振幅耦合(PAC)信号特征强烈区分无响应状态(最大峰值)与过渡状态(最大通过)，其最有可能反映丘脑极化水平的基本变化。因此，相位振幅耦合(PAC)可为麻醉剂诱导脑状态提供一敏感而合生理性逻辑的标记，其提供的信息比频谱特征单独提供的信息多。因此，一项研究对频谱特征不能完美地预测接受全身麻醉的病人丧失意识的多个案例进行了报告。在这些相同的数据中，跨频耦合(CFC)的测量值(最高峰值)可更准确地指示病人不能被唤醒的完全无意识状态。在手术室中，药物可以通过推注给药，药物输液速率可能突然变化，而且病人可能被手术自己唤醒，导致病人脑状态在以秒为尺度的时间内发生相应变化。这些状态的快速转变可能使使用常规方法估计的调制图案模糊。更快及更可靠的调制分析因此可导致全身麻醉的管理重大改进。由于常规方法需要多分钟的数据来作出一估计，因此不切实际；相反地，我们的方可在以秒为尺度的时间内估计跨频耦合(CFC)，此与这样的医疗应用相配。

由于跨频耦合(CFC)分析方法首先被引进神经科学，已经由大量数据表明，在健康与疾病中，跨频耦合(CFC)是大脑协调的基本机制。我们的方法解决现有技术面对的许多挑战性问题，并同时显著地改进统计效率及时间分辨率。此改良性能可为重要的新发现铺平道路(分析方法不足曾限制新发现)，而且可提高相位振幅耦合(PAC)分析的可靠性及效率，以使它们可用于医疗用途。

以下提供在一振荡波与一频率范围之间执行跨频分析的方法。虽然下述方法一般使用滤波技术以隔离一慢速频率与一宽带范围的频率，但据了解，以上呈献的状态空间建模技术可以被使用而代替所述滤波技术的至少一部份，以便更好地拟合一振荡及可能更准确地估计一振荡(比如以慢速波)与一范围的频率(比如一大约为5～50Hz的宽带范围)之间的关系。下文亦将解释上述状态空间模型及相关拟合技术可如何用于估计所述慢速波与所述宽带波之间的关系。

在过去的几年间发生了关于前皮质及后皮质在调节意识中扮演的角色的争议。“后皮质热区”假说主张，与前额皮层不同，后感觉皮质及感觉联合皮质是主要意识调节器。相反地，一些群体则认为前皮质后皮质互动是激发有意识意识的关键。所述争议最近已扩大至包括无意识的研究以及前后皮质区域的功能障碍与意识的丧失有关。

不同状态的无意识(比如麻醉、非快速动眼(NREM)睡眠以及昏迷)有截然不同的电生理信号特征。然而，它们的一个共同特征是慢速波活动，这个特征在脑电图(EEG)中表现为以大约1Hz的频率交替的大偏转。这些波被认为是皮质的高态及低态的大规模指标，在这些状态中，神经元在持续的去极化(高态)及超极化(低态)时期之间循环。在处于所述去极化高态时，神经元可能激发，但在处于所述超极化低态时，神经元沉默。

虽然慢速波活动在无意识期间在空间上广泛分布，最近的研究检查慢波动力学的空间分布与无意识状态如何关联。在非快速动眼(NREM)睡眠期间，后皮质区中的慢速波功率在非做梦睡眠期间比在做梦睡眠期间强。相反地，前皮质区已牵连麻醉剂诱导无意识。前额皮质的药理学刺激可恢复被麻醉动物的行为唤醒，与此同时，慢速振荡功率明显减低。麻醉剂诱导的慢速振荡在不同相位调制前额alpha振荡(“最大低谷”及“最大峰值”动力学)，视麻醉深度而定：当前额alpha振荡功率在所述慢速波的峰值(最大峰值)处为最强时，它似乎反映一深度麻醉状态，在这个状态下，监测对象不能从无意识中被唤醒，及时存在疼痛的刺激。这些数据表明，共享相似慢速波功率的不同行为状态可基于所述慢速波对较高频率节律的调制影响，相互分离。从慢速波功率到慢速波调制的这个属性的转变，是否可以用来调和后皮质及前皮质在调节无意识中扮演的角色？

虽然先前的关于无意识期间的跨频调制分析已着重于慢速波对所述alpha节律的影响，将慢速波解释成皮质的高态及低态的作法将暗示：由于高态及低态影响节律性及非节律性的神经活动，在所有频率处的脑电图(EEG)功率应耦合到所述慢速波。遵循这个想法，我们引入一新分析方法，以跟踪所述慢速波对在所述脑电图(EEG)中测量的一广泛范围的频率的调制影响。我们使用这个方法来分析在不同的丙泊酚诱导无意识状态下、在不同皮质区域中的慢速波调制。首先，我们发现前额最大峰值动力学是一个宽带现象，这表明最大峰值动力学反映基本皮质高态及低态。此外，我们发现，后皮质区及前皮质区藉不同无意识状态时的慢速波展示这种宽带调制：后皮质区藉较低程度的无意识状态时的慢速波展示这种宽带调制；而前皮质区则藉较高程度的无意识状态时的慢速波展示这种宽带调制。这个结果支持“麻醉剂诱导无意识并非单一状态而是多种状态，而每种状态显示不同程度的前皮质及后皮质破裂以及与其同时的有意识处理及行为可唤醒性中的相应差异”的想法。

我们使用如图19A中所示的带通滤波慢电压(0.1～4Hz)与带通滤波高频信号直线的互相关性，量化了跨频耦合。这个互相关性提供跨频耦合的一维概要：如果所述互相关性为正，这意味着当所述慢电压为正时，所述高频振幅较高(最大峰值耦合)；如果所述互相关性为负，这意味着当所述慢电压为负时，所述高频振幅较高(最大低谷耦合)。然后，我们可以改变用于所述高频振幅的频率，以查看不同频率如何与所述慢速波相关。以这个方法，我们可查看最大峰值及最大低谷动力学在过渡至无意识状态期间如何跨电极及频率改变。

图19B显示最大峰值及最大低谷动力学如何在一具代表性的监测对象身上发生：所述监测对象每14分钟被施加逐次增加剂量的丙泊酚，同时以一按钮咔嗒声响应来执行一听觉任务。当丙泊酚剂量够高时，所述监测对象对刺激停止反应(“意识丧失”)(LOC)。在施加5次逐增水平的丙泊酚之后减低剂量，最后所述监测对象又开始对刺激做出反应(“意识恢复”)(ROC)。

所述前皮质电极上的基于互相关性跨频耦合再次确定从先前队所述alpha节律所做工作获得的基本结果：在10Hz时，对所述慢速波的耦合作为最大低谷动力学开始，在较高剂量时移至最大峰值，然后在恢复期间向后过渡通过最大低谷。然而，通过观察10Hz以外的数据，我们可发现所述最大低谷动力学远在较高频率(高至大约20Hz)时的“意识丧失”(LOC)之前开始，这些频率显示：在监测对象保持意识的丙泊酚药物镇静期间，功率较高。这与“在最大低谷动力学期间病人可能清醒”的观测一致。

或许所述跨频耦合中最惊人的结果是，所述最大峰值图案扩展到所有频率(5～50Hz)。由于除了alpha波及慢速波之外，所述信号中没有窄带节律，这个结果表明非节律性宽带活动耦合至所述慢速波。这与“最大峰值动力学反映基本皮质高态及低态，其中在所述慢速波处于低谷期间，皮质活动(不论是否为有节律性)被关闭”的解释一致。

比较所述前及后电极之间的跨频耦合，我们可发现：在意识丧失不久之后，与前额最大低谷发生同时，一宽带最大峰值图案在后皮质波段形成。换句话说，即使在所述前皮质波段过渡到最大峰值动力学之前，后皮质波段已经显示与所述慢速波的宽带耦合。在所有监测对象中，后皮质最大峰值动力学在“意识丧失”(LOC)之后及在前皮质最大峰值动力学发动之前开始发动。这个结果在所述8～16Hz宽频耦合的拓扑图中按水平显示：在意识丧失之后不久，一组合的前皮质电极还是有最大低谷耦合，而后皮质电极上的电极显示最大峰值耦合。随着丙泊酚水平增高，所述前皮质电极开始参与所述最大峰值耦合。这表明，在丙泊酚剂量较高时向前皮质扩展之前，由所述慢速波对皮质活动进行的宽带消声便已在后皮质区中首先开始。

在其余的分析中，我们识别了四个感兴趣水平：(1)基线-在施加丙泊酚之前；(2)药物镇静-在监测对象停止反应之前的水平；(3)无意识低剂量-在监测对象停止反应之后的一个水平；以及(4)无意识高剂量-在未导致爆发抑制(一种超越了无意识状态所需要的条件的类昏迷状态)的最高丙泊酚水平。鉴于每一监测对象的“意识丧失”(LOC)及“意识恢复”(ROC)的时程不同，识别这些水平让我们能够结合跨不同监测对象的信息。

由于不同频率的跨频耦合的图案因电极、丙泊酚水平及监测对象而异，我们要识别识别最好地概括了这个活动的耦合图案。为了这么做，我们使用一非居中“主成分分析”(PCA)，以将所述跨频耦合图案分解为图20A中所示的多个主频众数。所述多个主频众数为跨频率的正交图案，它们连续地捕获通过频率、电极、丙泊酚水平及监测对象代表跨频耦合的矩阵中的总能量的多个较小部份。

在所有频率中，第一主频众数为正，因此所述第一主频众数代表所述跨频耦合中的宽带效应，而且其捕获所述总能量的78％。第二主频众数及第三主频众数依次在所述alpha及beta频率范围中有窄频峰值，并捕获少得多的能量。因此，宽带皮层活动对所述慢速波的夹带到目前为止，是丙泊酚麻醉期间的跨频耦合的最大贡献者。图20B中显示每一众数占总能量的百分比。

此显著的第一主频众数特别令人感兴趣，这是由于它捕获在所述个别监测对象概要中观测到的宽带最大峰值图案。尤其是，我们想使用所述第一主频众数来描述在丙泊酚剂量逐步增高期间的脑电图(EEG)的皮质发生器中的宽带最大峰值现象的空间分布。为了这么做，我们对所述源本地化脑电图(EEG)信号执行所述跨频耦合分析，并将产生的跨皮质位置的跨频耦合图案映射到得自以上所述的传感器空间分析的第一主频众数上。图21显示产生的映射，这些映射已变形为Freesurfer平均表面并跨监测对象平均。在这些表面善，一正映射(红色)代表对所述慢速波的宽带最大峰值的耦合。

与所述个别监测对象传感器空间结果一致，对所述慢速波的宽带最大峰值的耦合在所述“无意识低剂量”状态期间存在于后皮质区域。所述耦合在所述“无意识高剂量”状态期间加强，并扩展至前皮质区域。图22为一宽带频率调制图，此宽带频率调制由低频活动在不同状态期间在大脑的不同位置进行。我们可从图22中发现，在所述“无意识低剂量”状态期间，大脑顶叶、颞叶及枕叶都有宽带最大峰值耦合，但大脑额叶几乎没有对所述慢速波的宽带耦合。相反地，在所述“无意识高剂量”状态期间，大脑额叶参与其他脑叶的宽带最大峰值耦合。

这些结果共同表明，如图20中所示，在丙泊酚麻醉期间，对所述慢速波的宽带最大峰值耦合是跨频耦合动力学的主要贡献者。此外，如图21及图22中所示，后皮质区域在意识丧失后首先展示宽带最大峰值耦合，其次是前皮质区域在丙泊酚剂量高时展示宽带最大峰值耦合。

我们提议，我们在所述脑电图(EEG)信号中看到的、对所述慢速波的宽带耦合是区域皮质高态及低态的宏观指标。种群尖峰活动对所述神经功率频谱具有宽带效应，因此在皮质高态及低态期间，种群尖峰活动对所述慢速波的夹带应导致宽带功率在所述慢速波的一特定相位处为最强。这个想法获得大脑皮层电描记术中的宏观证据、以及显示高态与种群尖峰活动及宽带功率(<50Hz)同时增高相关的局部场电位数据的支持。虽然人们通常假设脑电图(EEG)中的慢速功率增加总是对应于皮质放电中的高态及低态，但我们的数据表明，在夹带所述种群尖峰活动的程度足以产生宽带调制之前，所述慢速波必须达到一临界振幅，该临界振幅可能因位置而异。后皮质区域越过这个阈值时的丙泊酚剂量，比前皮质区域越过这个阈值时的丙泊酚剂量来得低。所述慢速波临界阈值的想法与慢速波活动饱和假说一致，在所述假说中，慢速波功率的饱和伴随丘脑与感觉皮质之间的分离。我们假定所述皮质高态及低态本身为这个感觉隔离提供一个机制。根据这个解释，宽带最大峰值动力学表示，局部皮质活动正在被规模足以在所述脑电图(EEG)中被检测到的高态及低态破坏。

在全身麻醉期间对一病人的无意识状态进行监测的活动具有明显的实际效益，但为此活动识别有原则的方法仍然是一个悬而未决的问题。我们可区分病人可由外部刺激唤醒的(可唤醒)及病人不可被唤醒(不可唤醒)的不同无意识状态。一些新近数据间接地表明，当前皮质alpha节律在所述慢速波处于低谷(最大低谷)而为最强时，病人可被唤醒，但当前皮质alpha节律在所述慢速波处于峰值(最大峰值)而为最强时，病人不可被唤醒。这里，我们显示，当前皮质节律耦合至所述慢速波的峰值(最大峰值)时，前皮质及后皮质上的宽带谱功率亦耦合至所述慢速波。因此，不可唤醒性似乎与在前皮质及后皮质上的慢速波的宽带耦合相关。在病人无意识但可能可唤醒的较低麻醉剂剂量时，我们发现，后皮质参与宽带最大峰值动力学，而前额皮质参与(alpha频带)最大低谷动力学。因此，这些跨频耦合的不同空间图案似乎指示不同的无意识状态。对这些状态的知识可由麻醉师用于确定对临床情况适当的不同可唤醒或不可唤醒无意识状态。

所述“后皮质热区”及前-后皮质假说提出后皮质及前皮质电路在调节意识中扮演的截然不同的角色。麻醉对这个问题引入另一个维度，即无意识病人可能或不可能由外部刺激唤醒，视药物剂量而定。沿着这个维度，前额皮质似乎扮演一中心角色，在这个角色中，前额皮质的刺激可跨接意识。我们的结果似乎调合这些观点，因前皮质及后皮质的破坏与不同的无意识状态不谋而合。首先，“后皮质热区”及前-后皮质假说都预测后皮质区域的破坏会降低意识内容。在我们的结果中，当监测对象无意识但可被唤醒时，后皮质区域首先在丙泊酚剂量低时显示宽带最大峰值。这些后皮质区域显然惊人地与那些在非做梦睡眠期间有较高慢速波功率的监测对象相似。这与“后皮质区域调节意识内容，而后皮质区域在无意识期间受交替高态及低态的破坏”的解释一致。其次，前皮质在丙泊酚剂量较高时(处于与不可唤醒相关的状态时)出现宽带最大峰值耦合。这与“从行为唤醒到有意识状态需要前额皮质，而且前额皮质中的交替高态及低态阻扰这种行为唤醒”的解释一致。

本研究的主要结果是，全脑跨频耦合分析捕获可靠的变化，但却是在丙泊酚麻醉下在过渡到无意识及不可唤醒期间的大脑动力学中的空间改变的变化。我们将所述慢速波峰值的宽带耦合解释为宏观皮质活动中破坏皮质功能的高态及低态的宏观指标。这个功能破坏围绕处于不同无意识状态的前皮质及后皮质区域。我们的结果表明，麻醉之下的无意识并非一单一现象，而却是涉及大脑状态中的几个明显转变，这些转变破坏后皮质对意识的“内容”的处理，这些内容与前皮质唤醒及对这些内容的认识截然不同。

我们执行了一项研究，分析先前已经描述的数据。十名年龄介于18至36岁的健康志愿者参与了这项研究。在所述研究之前，取得了每位监测对象的结构磁共振成像(MRI)扫描图(Siemens Trio 3Tesla，T1-加权磁化准备快速梯度回波影像，1.3mm层厚，1.31mm面内分辨率，TR/TE＝2530/3.3ms，7翻转角)。以FreeSurfer图像分析套件执行了皮质重建。

在至少14分钟的闭眼休息之后，监测对象被施加了丙泊酚麻醉，丙泊酚剂量每14分钟增加至目标1、2、3、4或5μgmL^-1的效应位浓度。在苏醒期间，丙泊酚剂量以三个间隔14分钟的步骤减低，目标效应位浓度为比监测对象停止对所述听觉提示时的剂量低0.5μgmL^-1、1.0μgmL^-1及1.5μgmL^-1。使用一高密度(64频道)BrainVision MRI Plus***(BrainProducts GmbH产品)，以每秒5,000样本的采样率收集了脑电图(EEG)记录，而且电极位置数字化(Polhemus FASTRACK 3D)。

在整个研究过程中，监测对象执行了一项任务，他们以手指按压按钮来响应一听觉咔嗒声或听觉词。当一状态空间模型估计了一低于5％的正确反应概率并继续保持至少5分钟之后，将该状态定义为“意识丧失”(LOC)。同样地，当所述状态空间模型达到了一5％的正确反应概率并继续保持至少5分钟之后，将所述第一次苏醒定义为“意识恢复”(ROC)。

脑电图(EEG)预处理：通过目视检查，识别了带有严重、反复出现的人工产物的频道，并移除了这些频道。我们也根据eBridge软件包，移除了显示正被电气桥接的迹象的电极。

以下描述的宽频耦合分析使用已经在所述慢速频带(0.1～4Hz)中带通滤波的信号的版本，以及一组介于4～50Hz之间、每2Hz频带有一个1Hz过渡频带的较高频带。我们使用了一个零相位有限脉冲响应(FIR)滤波器，将整个过程的数据过滤为这些频带，然后降采样至每秒200个样本。我们选择了较高频带的上限50Hz，以避免线路噪声及捕获具有最大脑电图(EEG)功率的频率。

为进行传感器电平分析，我们在滤波之后根据第一近邻，应用一拉普拉斯(Laplacian)参考。.

由于所述丙泊酚剂量的阶梯式结构，在每一剂量施用后，监测对象身上药剂有恒定的预测效应位浓度的时长为大约12分钟。我们称这些时段为“水平”，并定义了四个感兴趣水平：“基线”、“药物镇静”、“无意识低剂量”及“无意识高剂量”。“基线”被定义为丙泊酚发挥作用之前的时段；“药物镇静”被定义为“意识丧失”(LOC)之前的水平；“无意识低剂量”被定义为“意识丧失”(LOC)之后的第一个全水平；而“无意识高剂量”被定义为丙泊酚的最高水平，或爆发抑制之前的水平(适用于进入爆发抑制的监测对象)。为了对这些水平进行分析，我们为每一监测对象选择了距每一水平30s的10个无伪影期。监测对象中有两位在丙泊酚的第一水平期间停止反应，因此他们没有药物镇静期。结果是，跨监测对象汇总的数字代表10位“基线”、“低剂量”及“高剂量”的监测对象，以及8位“药物镇静”的监测对象。

跨频耦合分析：在丙泊酚之下，高频振幅耦合至所述慢速波的峰值或低谷的可能性比其耦合至上升相位或下降相位的可能性高得多。因此，我们使用正(最大峰值：在所述慢速波处于峰值期间，高频振幅较高)或负(在所述慢速波处于低谷期间，高频振幅较高)的度量标准。尤其是，我们使用所述带通高频信号的瞬时振幅与所述带通慢电压(0.1～4Hz)之间的皮尔森(Pearson)相关系数：

其中V为一代表所述慢电压的时间序列的列向量，而A为一代表所述高频带的瞬时振幅的对应时间序列的向量。所述瞬时振幅以所述带通滤波信号的解析信号的大小计算。在计算所述相关性之前，通过减去平均数使每30s时段的瞬时振幅居中。V并未被明确地居中，这是由于其预期值为零。

为了求得跨空间(电极、源位置)的平均相关性及/或时间(多个30s时段)，需平均的量纲垂直地重叠于所述V中及方程式1中的A列向量中。由于所述相关性代表由每一信号的标准差(在所述分母

及

中)归一化的信号(在所述分子V^TA中)之间的协方差，垂直地重叠所述信号具有在执行所述归一化之前分别估计所述协方差及标准差的效果。

注意本分析并未涉及估计所述慢速波的相位，而所述估计是相位振幅耦合分析中的一个典型步骤。我们为了三个原因而选择不估计相位。首先，丙泊酚中的相位振幅耦合在零(慢速波峰，当所述信号为正时)和π(慢速波谷，当所述信号为负时)之间附近聚集。其次，所述慢速波被认为反映高态及低态，不论其是否很好地被描述为一正弦曲线。在丙泊酚诱导无意识的早期，所述慢速波的波形及频率经常是相当不规则。我们选择为所述慢速波使用一较宽频带0.1～4Hz，以捕获所述波形的更多非正弦特征。在丙泊酚麻醉之下，没有窄带delta频率节律活动(1～4Hz)，所以这个频率范围对所述带通信号的贡献是非节律性的，影响所述波形的形状。通过使用互相关性来描述所述高频振幅与所述非正弦低频波形的关系，我们能更好地捕获高频活动与所述宽带信号中明显的高态及低态波动之间的关系(图19A)。最后，将所述问题减为一维(正耦合及负耦合)简化了对一单一位置及频率的分析，使我们能够将所述分析扩展至其他频率及位置。

图19A显示一前皮质电极的计算，其计算在两个时窗期间所述慢电压(低频活动，LFA)与所述alpha频带(α，9～11Hz)之间的互相关性。对左方的时窗而言，alpha振幅在所述慢电压为负时较高，导致一负相关性(最大低谷)。对右方的时窗而言，alpha振幅在所述慢电压为正时较高，导致一正相关性(最大峰值)。我们可在所述过程中为所有的30s时段及为一范围的振幅频率(介于4～50Hz之间的2Hz频带)执行这个分析：这产生一调制图，该调制图为所选择的电极显示在所述过程整个期间每一频率与所述慢速波之间的跨频耦合。同样地，我们可为丙泊酚水平内的每一电极计算互相关性，并显示所述alpha频带与所述慢速波之间的跨频耦合的空间分布为丙泊酚水平的函数。

我们执行了跨振幅频带、时间、电极及监测对象的宽频耦合分析。我们以两种方式来表征时间。第一个方式为一基于时段的分析，其中在所述时段中，所述耦合每30s时间被分别评估一次，并产生调制图，如图19B中所示的前皮质及后皮质概要：这些调制图表达跨几个电极的平均相关性，这些平均相关性在右边的头皮图中指示。用于表征时间的第二个方式基于丙泊酚水平，其中在每一丙泊酚水平，所述耦合在10个30s时间的无伪影期期间被评估。因此，每一耦合值代表超过300秒的数据，而且我们有每一电极、水平及监测对象的这样的值。这个基于水平的分析用于图19B中的拓扑图，并用于以下所述的“主成分分析”(PCA)。

我们使用MNE-python工具箱(27.martinos.org/mne/)，通过最小范数估计执行了源定位。使用高于65Hz的谱功率，根据所述脑电图(EEG)估计，假设噪声协方差为对角。所述源空间为所述皮质表面的Freesurfer重建的的ico-3抽取，而所述源被约束为与所述皮质表面垂直。使用一3层边界元模型计算所述正向模型，消除距颅骨内表面5mm范围内的的任何源。

为了在源空间执行所述跨频耦合分析，使用以上所述方法，对所述慢速频带中的带通信号以及4～50Hz之间的每一2Hz频带执行源定位。对所述皮质表面的表达(图21)而言，将以上所述的跨频耦合分析应用于所述四个感兴趣水平的每一监测对象的源位置。脑叶中的平均跨频耦合(图22)使用了相同的技术，但结合了一脑叶中的所有源的信号：用于所述脑叶中的所有源的10个时间(适用于每一水平)垂直地重叠于所述V中及方程式1中的A列向量中。

非居中“主成分分析”(PCA)与“映射至源空间”：为了识别跨频耦合的图案(它们解释了大多数所观察的特征)，我们对所述传感器空间耦合图案执行了跨丙泊酚水平及监测对象的非居中“主成分分析”(PCA)。所述四个感兴趣水平(基线、药物镇静、无意识低剂量及无意识高剂量)的基于水平的耦合结果被合并为一集合矩阵A(f，i)，其中f指数化所述振幅频带，而i指数化所述丙泊酚水平、电极及监测对象(即：行包含所有重叠的水平、电极及监测对象)。接着产生了所有监测对象的两个电极在所述四个感兴趣水平期间的跨频耦合结果：所述分析的目的在于将这些图案分解为捕获它们跨频率的主要特征的不同分量。主成分分析涉及一奇异值分解：

A＝USV^T (91)

在所述分解中，U为一举证，其列为主要众数。第一列为第一主要众数，意思是在所述A矩阵中捕获最多能量的跨频图案。第二列为第二主要众数，其在所述第一众数被去除之后捕获所述数据中最多的能量。所述主要众数为正交，而且都有单位长度。为了保留从正值到最大峰值的映射及从负值到最大低谷的映射，在必要时将所述主要众数翻转(即乘以-1)，使最大元素必为正数。所述S矩阵为对角矩阵，而且可用于估计总能量百分比，总能量百分比由第j个(j^th)众数解释：

注意我们选择了使用分解总能量的非居中“主成分分析”(PCA)，而不是使用分解方差的居中“主成分分析”(PCA)。在原点具有特殊意义(而如果通过减去所述平均数来居中所述数据，所述原点将失去特殊意义)的情况下，非居中“主成分分析”(PCA)特别有用。以跨频耦合的图案而言，原点代表所述脑电图(EEG)信号不包含在任何频率处与所述慢速波耦合的情况。相反地，所述平均数代表跨监测对象、丙泊酚水平及电极的、与所述慢速波的平均耦合。我们觉得，所述原点比所述平均数更具意义，因此所述非居中“主成分分析”(PCA)比常用的居中“主成分分析”(PCA)更易于解释。

为了量化所述主要众数如何在源空间中表示，我们每一监测对象的源空间图案映射到所述第一传感器空间主要众数上。首先，来自所述源空间的耦合结果结合到一集合A^source(f，i)矩阵，其中i指数化所述丙泊酚水平、源位置(或脑叶)及监测对象。接着，A^source到U的映射为：

P＝U^TA^source (93)

生成的P矩阵的第一行包含源空间耦合图案(每一丙泊酚水平、源位置及监测对象)到所述第一主要众数的映射。由于所述第一主要众数跨频率相对恒定(见图20)，映射到此众数表示宽带的宽频耦合。此外，由于所有频率的第一主要众数都是正数，正映射反映宽带最大峰值耦合(所有频率耦合到所述慢速波的峰值)，而负映射反映宽带最大低谷耦合(所有频率耦合到所述慢速波的低谷)。

所述(subject×level×lobe)(监测对象x水平x脑叶)耦合图案到所述第一主要众数的映射跨监测对象平均，以得到图22中所示的估计值，这些估计值代表所述第一主要众数对该脑叶中的跨频耦合的平均贡献。所述95％置信区间乃使用一跨监测对象的靴圈求得。

为了总结跨监测对象的整个皮质表面的结果，我们使用MNE-python工具箱，将所述(sub ject×level×sourcelocation)(监测对象x水平x源位置)耦合图案到所述第一主要众数的映射变形为Freesurfer-平均表面，而且我们跨监测对象平均了生成的映射图。所述生成的映射图估计所述选择的众数对每一皮质位置的跨频耦合的贡献(见图21)。

如以上所述，所述状态空间模型可用于检测一慢速分量与一范围的宽带频率之间的跨频关系。为了在一基于状态空间模型的跨频耦合分析中使用所述状态空间模型，可以对以下所示的方程式(3)及(4)进行参数替换，以仅仅获得一慢速分量。这里重复方程式(3)及(4)，供作参考：

其中M＝[1 0]，

Q＝Q_s，

接着可以使用一如以上所述的“期望最大化(EM)算法”拟合所述模型。然后，可以在以上所述的跨频耦合估计方法中，将所述向量时间序列

的两个元素(实部及虚部)作为所述带通滤波慢电压(即；所述慢速信号)使用。所述高频带中的振幅还是使用带通滤波及取所述解析信号的大小来计算，即：使用所述希尔伯特变换(Hilbert transform)方法来计算。

接着可使用方程式(90)来量化所述跨频耦合，这里重复方程式(90)，供作参考：

其中V为一列向量，其包含

的第一或第二元素，而A为一代表所述高频带的瞬时振幅的对应时间序列的向量。在计算所述互相关性之前，应通过减去平均数将所述瞬时振幅居中。V不需被明确地居中，这是由于其预期值为零。

接着，可以为

的每一元素产生一类似图19B的图。对应于

的第一元素(即所述实部)的分析将捕获所述高频活动与所述慢速波的峰值(正，最大峰值)及低谷(负，最大低谷)的耦合，而所述图案在品质上非常相似于跨时间、频率及电极的原始分析。对应于

的第二元素(即所述虚部)的分析将捕获所述高频活动与所述慢速波的下降相位或上升相位的耦合。在一些的实施例中，在拟合所述慢速振荡之后剩余的残差可以被滤波，并用于代替所述带通原始脑电图(EEG)信号。使用残差可去除来自所述宽带信号的估计慢速波的任何高频元素。

以上所述的基于状态空间模型的跨频耦合分析方法的一个潜在问题是，如果所述信号包含强振荡或所述慢速范围以外的局部频谱功率，所述“期望最大化(EM)算法”可能执行不理想：例如，其可能将所述频率f_s移到与慢速波相关的范围之外(低于大约1Hz)，及/或增加所述方差

以捕获所述高频功率及所述慢速波。

解决上述问题的一种潜在方法是选择所述状态空间模型的振荡分量数目，并使用生成的隐藏状态(其中之一应对应于所述慢速波，而其他隐藏状态将对应于高频振荡)。接着，可计算所述高频分量的瞬时振幅与所述慢速分量的实部/虚部之间的互相关性。然后，可确定所述模型的残差(反映所述信号的“非振荡”分量)，以及执行所述希尔伯特变换(Hilbert transform)方法，以求得所述高频残差(如以上所述，在一集合的频带中)的瞬时振幅。接着，可计算这些频带与所述慢速分量的实部/虚部之间的互相关性。

另一潜在解决方案是，非但不使用所述希尔伯特变换(Hilbert transform)方法来计算所述高频信号的瞬时振幅，我们可使用平铺更高频率的一集合高频分量，从所述状态空间模型求得所述瞬时振幅。换言之，可以以所述慢速波

的分量以及一集合的高频隐藏状态

(其中i指数化频率-例如5～50Hz之间)，拟合所述标准状态空间方程式(3，4)。由于每一高频分量并非一般存在窄带振荡，在这种情况下，所述“期望最大化(EM)算法”可能表现得不可预测。因此，为了使此方法起作用，用户将需为所述高频带的频率f_i(所述频带的峰值位置)及自回归系数a_j(峰的锐度)选择一先验值，而不在所述“期望最大化(EM)算法”中更新它们。

参阅图2、图11以及图23，图23显示一使用拟合状态空间模型来执行跨频分析的示例性流程2300。所述流程2300可实施为至少一个存储器上的指令。所述指令可接着由至少一个处理器执行。

在步骤2304，所述流程2300可接收脑电图(EEG)数据。所述脑电图(EEG)数据可以接收自一病人监测装置，比如接收自以上连同图2描述的病人监测装置202。所述脑电图(EEG)数据可作为一波形解释。所述流程2300可接着进入步骤2308。

在步骤2308，所述流程2300可拟合一状态空间模型至所述脑电图(EEG)数据。步骤2308可包括流程1100的步骤的子集，例如步骤1108-1128。所述流程2300可接着在随后的步骤中使用所述拟合模型来执行跨频分析。然后，所述流程2300可进入步骤2312。

在步骤2312，所述流程2300可对所述拟合模型的一慢速波分量及一alpha波分量执行跨频分析。所述流程2300可计算两个波分量之间的相位振幅耦合的一个或多个参数表达式。所述流程2300可使用方程式(33)进行“状态空间跨频分析”(SSCFA)。可选择地(或此外)，所述流程2300可使用方程式(49)进行“双状态空间跨频分析”(dSSCFA)。“状态空间跨频分析”(SSCFA)及“双状态空间跨频分析”(dSSCFA)都能提供所述alpha波分量与所述慢速波分量之间的相位振幅耦合的准确估计，更明确地说，它们能在不依赖时间分箱的情况下，估计所述慢速波分量的相位如何调制所述alpha波分量的振幅。如以上所述，所述“双状态空间跨频分析”(dSSCFA)对跨所述脑电图(EEG)时窗的调制参数施加一时间连续性约束。换言之，基于以时间约束进行“双状态空间跨频分析”(dSSCFA)计算。由“状态空间跨频分析”(SSCFA)及“双状态空间跨频分析”(dSSCFA)提供的参数表达式的相位振幅耦合(PAC)提供所述慢速波分量与所述快速波分量之间的最大峰值及/或最大低谷耦合的指标。所述流程2300可接着进入步骤2316。

在步骤2316，所述流程2300可执行所述慢速波分量与所述脑电图(EEG)数据的一范围的频率之间的跨频分析。例如，所述慢速波分量与所述5～50Hz频率之间的关系。所述慢速波分量与所述范围的频率之间互相关性可使用方程式(90)来计算。在方程式(90)中，V为一列向量，其包含

的第一或第二元素，第一或第二元素可从所述拟合模型取得。A为一代表一高频带(即5～50Hz)的瞬时振幅的对应时间序列的向量。可通过应用一带通滤波器使所述高频带通入所述脑电图(EEG)数据，确定A。在一些的实施例中，在拟合所述慢速波分量至所述脑电图(EEG)数据之后剩余的残差可以以一用于带通所述高频带的带通滤波器来带通。使用所述残差来确定A可能去除来自所述宽带信号的估计慢速波的任何高频元素，这可能提供一更佳信号。所述流程2300可接着使用方程式(90)估计所述慢速波与所述高频带之间的互相关性。然后，所述流程2300可进入步骤2320。

在步骤2320，所述流程2300可将所述“状态空间跨频分析”(SSCFA)、所述“双状态空间跨频分析”(dSSCFA)及/或所述慢速波与所述高频带模型之间的互相关性输出到一显示器，供一用户查看，及/或输出到一存储器供存储及/或供另一流程使用。其他流程可进一步处理计算求得的跨频分析值(即：所述“状态空间跨频分析”(SSCFA)、所述“双状态空间跨频分析”(dSSCFA)及/或所述慢速波与所述高频带模型之间的互相关性)，并产生一报告，以使这些值能更容易由用户理解，比如将计算自多个电极度数的多个“状态空间跨频分析”(SSCFA)值映射到一基于每一电极的位置的大脑模型。所述报告可以输出到一显示器，供一用户查看，及/或输出到一存储器供存储及/或供另一流程使用。然后，所述流程2300可结束。

以上呈献的多种配置纯粹为示例，它们绝非意在限制本发明的范围。对本领域的普通技术人员而言，本文中描述的配置的变体是显而易见的，这些变体属于本专利申请的预期范围。尤其是，可以选择以上所述的一个或多个配置中的特征来创建包括可能未在上文明确描述的特征的子组合的替代配置。此外，可以选择及结合以上所述的一个或多个配置中的特征来创建包括可能未在上文明确描述的特征组合的替代配置。对本领域中的技术人员而言，在对本申请书进行整体审查后，适合这样的组合及分组合的特征是显而易见的。本文中描述的及权利要求提及的主题预期覆盖及包含所有合适的技术变化。

Claims

1.一种用于估计一病人大脑状态的至少一个振荡分量的方法，所述方法包括以下步骤：

接收脑电图(EEG)数据；

拟合一状态空间模型至所述脑电图(EEG)数据，所述状态空间模型包括至少一个振荡分量；

根据所述状态空间模型，估计所述脑电图(EEG)数据的跨频耦合的至少一个值；

根据跨频耦合的所述至少一个值，产生一报告；以及

将所述报告输出到一显示器及一存储器的至少其中之一。

2.根据权利要求1所述的方法，其中跨频耦合的所述至少一个值包括所述状态空间模型中包括的一第一波分量及一第二波分量之间的相位振幅调制的一估计值。

3.根据权利要求8所述的方法，其中所述第一波分量是一慢速波分量，而所述第二波分量是一alpha波分量。

4.根据权利要求3所述的方法，其中相位振幅调制的所述估计值是根据所述慢速波分量的一相位及所述alpha波分量的一振幅计算。

5.根据权利要求1所述的方法，其中跨频耦合的所述至少一个值包括所述状态空间模型中包括的一振荡分量与所述脑电图(EEG)数据的一频带之间的一互相关值。

6.根据权利要求1所述的方法，其中所述拟合所述状态空间模型的步骤包括拟合所述至少一个振荡分量中包括的每一振荡分量的谐波频率。

7.根据权利要求1所述的方法，其中所述至少一个振荡分量包括一alpha分量，所述alpha分量与一中心频率的先验分布相关。

8.根据权利要求7所述的方法，其中所述先验分布是一“冯·米塞斯”(Von Mises)先验分布。

9.根据权利要求1所述的方法，其中所述状态空间模型的一阻尼因数以一先验分布约束。

10.根据权利要求9所述的方法，其中所述先验分布是一beta分布。

11.一种用于估计一病人大脑状态的至少一个振荡分量的***，所述***包括：

用于接收脑电图(EEG)数据的多个传感器；

用于接收所述脑电图(EEG)数据及执行下列步骤的一处理器，这些步骤包括：

根据跨频耦合的所述至少一个值，产生一报告；以及

用于显示所述报告的一显示器。

12.根据权利要求11所述的***，其中跨频耦合的所述至少一个值包括所述状态空间模型中包括的一第一波分量及一第二波分量之间的相位振幅调制的一估计值。

13.根据权利要求18所述的***，其中所述第一波分量是一慢速波分量，而所述第二波分量是一alpha波分量。

14.根据权利要求13所述的***，其中相位振幅调制的所述估计值是根据所述慢速波分量的一相位及所述alpha波分量的一振幅计算。

15.根据权利要求11所述的***，其中跨频耦合的所述至少一个值包括所述状态空间模型中包括的一振荡分量与所述脑电图(EEG)数据的一频带之间的一互相关值。

16.根据权利要求11所述的***，其中所述拟合所述状态空间模型的步骤包括拟合所述至少一个振荡分量中包括的每一振荡分量的谐波频率。

17.根据权利要求11所述的***，其中所述至少一个振荡分量包括一alpha分量，所述alpha分量与一中心频率的先验分布相关。

18.根据权利要求17所述的***，其中所述先验分布是一“冯·米塞斯”(Von Mises)先验分布。

19.根据权利要求11所述的***，其中所述状态空间模型的一阻尼因数以一先验分布约束。

20.根据权利要求19所述的***，其中所述先验分布是一beta分布。