多变量时序响应函数工具箱：一个用于关联神经信号与连续刺激的MATLAB的工具箱

了解大脑如何在自然环境中处理感觉信号是二十一世纪神经科学的关键目标之一。虽然脑成像和侵入性电生理学将在这一努力中发挥关键作用，但具有高时间分辨率的非侵入式宏观技术，如EEG和MEG也将发挥重要作用。但在确定如何最好地分析这种对复杂的、时变的和多变量的自然感觉刺激的神经反应方面仍存在挑战。应用系统识别技术将神经元放电活动与复杂的感觉刺激联系起来已经有很长的历史，这种技术现在越来越多地应用于EEG和MEG数据。

一个特别的例子涉及到拟合一个滤波器（通常被称为时序响应函数），它描述了感觉刺激的某些特征和神经反应之间的映射。在这里，我们首先简要回顾了这些系统识别方法的历史，并描述了一种用于推导时序响应函数的具体技术，即正则化线性回归。然后，我们介绍了一个新的开源工具箱来进行这种分析。我们描述了如何使用它来推导（多变量）时序响应函数，描述刺激和响应之间在两个方向上的映射。我们还解释了正则化分析的重要性，以及如何为特定的数据集优化这种正则化。然后，我们具体概述了工具箱是如何实现这些分析的，并提供了几个工具箱可以产生的结果类型的例子。最后，我们考虑了该工具箱的一些局限性以及未来发展和应用的机会。本文发表在Frontiers in Human Neuroscience杂志。

介绍：

传统的对人类感觉处理的电生理学研究集中在短暂、独立的刺激这一相当特殊的情况，因为需要对离散的感觉事件进行锁时，以估计事件相关电位（ERPs）。其目的是通过用一个瞬时的、类脉冲的刺激对系统进行卷积，并对几百次时间锁定的反应试验进行平均，以估计感觉系统的脉冲响应函数。这种方法已经被广泛用于研究人脑如何处理各种生态事件，甚至是那些以连续、动态的方式发生的事件（如人类语言）。然而，此类ERP研究中使用的语音刺激通常由单个音素或音节组成，因此并不完全反映自然、有联系、持续且丰富的具有词汇复杂性的语音。最近的研究已经开始使用更自然且更长的语音刺激，将分析重点放在测量同一语段多次重复的神经反应相位上。虽然这种方法揭示了对语音处理神经生理学的新颖而重要的见解，但它并不有助于描述系统的响应功能，而且无论如何，它是一种间接的衡量大脑如何随着时间推移对刺激进行响应的方法。

研究连续刺激神经诱导的直接方法是对一个函数进行数学建模，该函数描述了刺激的特定属性映射到神经反应的方式，这种技术被称为系统识别（System identification，SI）。虽然有几类模型可以实现这一目的，但最直接的一类是线性时不变（Lineartime-invariant，LTI）系统。尽管人脑既不是线性的也不是非时变的，但这些假设在某些情况下是合理的，并允许系统由其脉冲响应表征。 一种被称为 “反向相关（reverse correlation）”的SI方法已经成为神经生理学中表征LTI系统的常用技术，这种方法在视觉和听觉动物电生理学中早已确立。这种技术近似于被研究的感觉系统的脉冲响应，只是它不需要使用离散的刺激。虽然这有点类似于计算ERP，但必须考虑一些重要的区别：

（1）通过反向相关得到的反应函数只反映了系统对实验者定义的特定刺激参数的反应，而不是整个事件；

（2）不像锁时平均，反向相关假设系统的输入-输出关系是线性的；

（3）反向相关比ERP（容易受到时域坏点的影响）收敛于时域上更精确的系统脉冲响应估计。

最简单形式的反向相关可以通过LTI系统的输入和输出之间的直接交叉相关来实现。虽然这种方法已被用于研究语音如何在人脑活动中编码，但它更适合于被高斯白噪声等随机过程调节的刺激。因此，这种方法在动物模型中的大多数实例传统上都使用白噪声刺激。这项工作甚至促进了研究人员研究这种随机信号是如何在人脑中编码的。也就是说，人类的大脑已经进化到可以处理几乎不符合白噪随机过程的生态相关刺激。例如，在人类神经科学研究的背景下，鉴于神经元对更复杂的刺激有不同的反应，要正确理解大脑如何处理自然语音，需要在实验室中使用自然语音作为刺激。因此，由于SI（系统识别）方法的发展，例如normalizedreversecorrelation（NRC）、ridgeregression和boosting，使用动物模型的研究人员已将重点转向使用更自然的刺激研究大脑。每一个技术使用不同的先验，但都收敛于相同的理论解决方案，关键是对非白噪刺激给出一个无偏的脉冲响应估计。这激励了研究人员对各种动物模型中听觉皮层神经元的“频谱接受区”进行定性。因此，对人类语音处理过程感兴趣的研究人员已经开始建立响应函数模型，描述自然语音的属性（如包络或频谱图）与群体（动物和人类）响应之间的线性映射。也有类似的努力来模拟与更自然的视觉刺激属性（如运动）有关的人类神经反应的响应函数，这也是受到以前单位电生理学工作的启发。

上述研究大多是正向对刺激-响应映射函数进行建模（前向建模）。然而，这种映射也可以反向建模（后向建模），为研究刺激特征如何在神经反应措施中编码提供了一种补充方式。与正向模型不同，后向模型参数不容易在神经生理学上解释，但可用于从神经反应中重建或解码刺激特征，这种方法被称为“刺激重建”。

这种方法比正向建模方法有几个优点，特别是在使用多通道系统（如EEG）从集群反应中记录时。首先，由于重构投射到刺激域，它不需要预先选择神经反应通道。事实上，将所有反应通道纳入后向模型是有利的，因为重建方法对不相关的通道给予较低的权重，同时允许模型使用可能被特征选择方法排除的通道来捕捉额外变量。其次，后向模型可以对彼此高相关的反应通道之间的重要信号差异提供更高的敏感性（如脑电图经常出现的情况）。它可以做到这一点，因为分析同时映射所有响应通道的数据（即以多变量的方式），因此它可以识别数据中任何的信道间关联性。相反，在进行正向建模时，每个分析都是单变量的，因此会忽略其他EEG通道的数据。第三，在神经反应中没有明确编码的刺激特征可以从被编码的相关输入特征中推断出来。这可以防止模型将资源分配给冗余的刺激信息的编码。刺激重建法以前被用来研究各种动物模型的视觉和听觉系统。最近，它被用于研究人脑中的语音处理，使用颅内和非侵入式电生理学。

虽然某些研究小组现在经常使用SI（系统识别）来研究人脑中的感觉处理，但这种方法因为与它使用相关的挑战，还没有被整个神经科学界广泛采用。本文的目的是介绍一个最近开发的SI工具箱，它提供了一个直接和灵活的正则化线性（ridge）回归的实现。我们首先总结了这一技术的数学基础，然后提供了一些如何使用该工具箱的具体例子，最后讨论了它的一些应用和重要考虑因素。

正则化线性回归

前向模型：时序响应函数估计

前向模型有时被称为生成或编码模型，因为它们描述了系统如何产生或编码信息。在这里，它们将被称为时序响应函数（Temporal response functions，TRFs）。有许多方法可以从数学上描述一个系统的输入与它的输出之间的关系。一种常用的方法是假设系统的输出通过一个简单的线性卷积与输入相关。在输出由N个记录通道监测的感觉系统的背景下，让我们假设在时间t=[1…T]和通道n采样的瞬时神经反应r(t, n)由刺激属性s（t）与未知通道特定的时序响应函数w(τ, n)的卷积组成（见公式1）。

其中ε(t, n)是模型无法解释的每个通道的响应残差。从本质上讲，TRF可以被认为是一个滤波器，它描述了正在进行的刺激与正在进行的神经反应之间的线性转换。时序响应函数w(τ, n), 描述了相对于刺激特征s(t)的瞬时发生而言，在特定的时间滞后范围τ中的转换。

以语音为例，s(t)可以是对每一时刻的语音包络的测量，r(t, n)可以是通道n处相应的EEG反应。计算w(τ, n)所需的时间滞后范围通常使用捕捉ERP的皮质反应成分范围，例如-100-400ms。

在-100毫秒的TRF的结果值，将反映前100毫秒的语音包络和神经反应之间的关系（显然这应该有一个零振幅），而100毫秒处的TRF将反映语音包络的振幅的单位变化将如何影响100毫秒后的EEG。时序响应函数w(τ, n)，是通过最小化实际神经反应r(t, n)和卷积预测之间的平均标准误（Mean-squarederror，MSE）来估计的(见公式2)。在实践中，我们通常使用反向相关来解决，可以用公式3的矩阵运算来实现。其中，S是刺激属性的时滞序列，定义如公式4。

τmin和τmax分别代表最小和最大的时滞（样本中）。在S中，每个时滞是按列排列的，非零时滞被填充为零以确保因果关系。计算TRF的窗口被定义为 τwindow = τmax- τmin，因此S的大小为T × τwindow。为了在回归模型中包括常数项（y-截距），在S的左边连接了一列1。在方程（3）中，变量r是一个包含所有神经反应数据的矩阵，按列排列通道（即一个T×N矩阵）。由此产生的时序响应函数，是一个 τwindow×N的矩阵，每一列代表从s到每个通道的神经反应的单变量映射。

这里重要的一点是，这种分析明确地考虑到了刺激的自协方差结构。在非白噪刺激中，如自然语音、声学信号的强度受时间调节，这意味着它在非零时滞自相关。

语音包络和相应的神经反应的简单交叉关联将导致脉冲响应函数的时间污染。这里的解决方案（公式3）是将刺激的自协方差结构从模型中划分出来，这样就可以消除不同时间点之间的相关性。这样做的TRF（时序响应函数）方法比简单的交叉相关方法更不容易产生时间上的污染。这在下一节的一个工作例子中得到了证明。

正则化

在计算TRF时，一个重要的考虑因素是正则化，即引入额外的信息来解决不适定估计并防止过拟合。不适定估计与转置自协方差矩阵STS有关。当以有限精度求解时，矩阵转置易数值不稳定。换句话说，如果前者存在不适定估计的情况，STS的微小变化（如离散化引起的取整误差）可能导致w的巨大变化。换句话说，w的估计可以有非常高的方差。当刺激物代表一个随机过程时，这通常并不适用，因为STS将是满秩的。然而，像语音这样的非白噪刺激的自相关特性意味着它更有可能是奇异的（即行列式为0）。通常情况下，对不良条件矩阵的数字处理包括通过增加一个偏置项或 "平滑解 "来减少估计的方差。具体来说，由于总体估计误差是由偏置项（即估计的期望值与真实值之间的差异）和方差项组成的，因此可以故意增加偏置项，以减少估计值的（高）方差，从而降低总体估计误差。

加入这个平滑项也解决了过拟合的问题。反向相关分析无法得知它被要求建模的数据生物学性质。因此，如果没有正则化，产生的TRF将是用于拟合的特定数据集的最佳拟合标准（如最小二乘法误差）。

而且，鉴于这些数据是 "嘈杂的"，TRF可以显示出生物学上难以置信的特性，如非常高的频率波动。使用这个TRF来预测新数据可能会导致次优的性能，因为高频波动不一定与新数据中的 "噪声"很好地对应。换句话说，TRF已经被 "过拟合"到训练中使用的特定数据集。正则化的作用是防止在低方差维度对这种高频、特定数据集的噪声进行过度拟合。例如，它可以通过规定相邻的TRF值之间的巨大差异来做到这一点，从而迫使TRF变得更加平滑。这使得TRF对用于拟合的数据不过于具体，可以帮助它更好地泛化到新的数据。

在实践中，通过在求逆前对STS的对角线进行加权，可以同时解决不适定问题和过拟合问题，这种方法被称为Tikhonov正则化或岭回归。

其中是I特征矩阵，λ是平滑常数或 "岭参数"。岭参数可以使用交叉验证来调整，以使r(t, n)和之间的相关性最大化。TRF优化将在下一节中详细描述。虽然这种形式的岭回归通过添加惩罚因子限制TRF值与零的距离来对所产生的模型进行平滑性约束，但另一种选择是对W的每两个相邻项之间的差异进行四分法惩罚。

Tikhonov正则化（公式5）通过对TRF估计的平滑化来减少过拟合，这种方式对感兴趣的信号振幅不敏感。然而，二次方法（公式6）在保留信号振幅的同时减少了非样本误差。因此，与Tikhonov正则化相比，这种方法通常会导致对系统响应的更好估计。

多变量分析

上一节重点介绍了将单一的、单变量的输入刺激特征（例如，语音刺激的包络）分别与多个记录通道中的每一个联系起来的具体情况。然而，自然界中大多数复杂的刺激都不是作为简单的单变量特征来处理的。例如，当听觉语音进入耳朵时，信号被耳蜗转化为由多个频段组成的频谱表征，沿听觉通路投射。听觉系统将每一个频段映射到皮层水平上测得的神经表征。这个过程可以用TRF的多变量形式（即mTRF）来模拟。

实际上，有可能定义一个mTRF，将一个多变量刺激特征线性地映射到每个记录通道上。使用上述例子，让s(t, f )代表f=1...F频段的语音信号的频谱图。为了得出mTRF，刺激滞后矩阵S（公式4）被简单地扩展，每一列被替换为F列，每列代表一个不同的频段（一个T × Fτwindow矩阵）。由此产生的多变量时序响应函数w(f , τ , n)，将是一个Fτwindow × N的矩阵，但可以很容易地分解，从而使每个自变量被表示为一个单独的维度（即一个F × τwindow × N矩阵）。这里，常数项是通过连接S左边的F列而被纳入的。

多变量TRF分析的一个重要考虑是使用哪种正则化方法。方程（6）中的二次正则化项被设计为强制执行平滑度约束并保持沿时间维度的信噪比。对于λ高值，这种方法会导致跨频率的污染；因此它不会产生每个频段的TRF的准确表示。在这种情况下，通常最合适的做法是使用特征矩阵进行正则化（公式5），以避免在mTRF的非时间维度上强制执行平滑度约束—尽管在某些情况下，这实际上可能是我们所希望的。

后向模型：刺激重建

前面几节描述了如何对刺激和神经反应之间的线性映射进行正向建模。虽然这种方法可以扩展到适应多变量的刺激特征，但它是次优的，因为它把每个神经反应通道当作一个独立的单变量特征。另一方面，后向模型通过利用多变量背景下的所有可用神经数据，得出反向刺激-反应映射。后向模型有时被称为识别或解码模型，因为它们试图通过从神经反应中解码刺激特征来反向生成数据。这里，它们将被简单地称为解码器。

解码器的建模方式与TRF大致相同。假设解码器g( τ, n)表示从神经反应r(t, n)回到刺激物s(t)的线性映射。这可以用离散时间表示为公式7。其中，s(t)是重构的刺激物属性。在这里，解码器在指定的时滞τ的范围内整合神经反应。理想情况下，这些时滞将捕获神经数据的窗口，以优化刺激物的重建。通常，用于重建的最有信息量的滞后期与用于捕捉前向TRF的主要成分的滞后期相称，但在反向上除外，因为解码器有效地在时间上向后映射。为了扭转先前TRF例子中使用的滞后期（ τmin= -100ms，τmax = 400ms），τmin和τmax的值被调换，但它们的符号保持不变，即 τmin= -400ms，τma = 100ms。解码器g(τ, n)是通过最小化s（t）和预测s(t)的MSE来估计的（见公式8）。

与TRF方法类似，解码器是通过公式9的矩阵运算来计算的。其中R是响应矩阵的时滞序列r。为简单起见，我们将为单通道响应系统定义R。

如前所述，这可以扩展到N个通道系统的多变量情况，方法是将R的每一列替换为N列（每列代表一个单独的记录通道）。对于N个通道，R的尺寸将是T×Nτwindow。常数项是通过将R左边的N列1连接起来而纳入的。在语音的背景下，刺激变量s代表一个列矢量（如包络）或一个T×F矩阵。由此产生的解码器g，将分别是一个Nτwindow样本的向量或一个Nτwindow×F矩阵。虽然对解码器权重的解释不像TRF那样直接，但人们可能希望将其分开（例如，N×τwindow×F），以检查每个通道在特定时间滞后的相对权重。通道权重代表了每个通道为重建提供的信息量，也就是说，信息量大的通道得到的权重更大，而提供信息量小或没有的通道得到的权重更接近于零。在公式（9）中，使用了Tikhonov正则化，因为它假定神经反应数据是多变量的。如上所述，在重建方法中消除了神经反应通道之间的相关性的任何偏差。在实践中，这是通过除掉神经反应的自协方差结构来实现的（见方程9）。因此，通道加权变得更加局部，因为通道间的冗余不再被编码在模型中，这使得它比正向TRF方法和交叉相关方法更具优势。

MTRF 工具箱: 实施和功能

本节概述了如何在MATLAB中使用mTRF工具箱实现正则化线性回归（https://sourceforge.net/projects/aespa/). 具体来说，它描述了如何在单变量和多变量数据集上进行训练和测试，以及应如何为特定目的优化所产生的模型。

训练

在mTRF工具箱中，使用一个简单的函数mTRFtrain来实现对给定数据集的刺激-反应映射的建模。该函数计算单变量或多变量岭回归，如上节所述（公式5、6和9）。该模型可以通过两种不同的方式在数据集上进行训练：（1）在每个试次上单独训练并在M个模型上取平均值，或者（2）在试次的串联上训练。这两种方法产生相同的结果，因为数据是用线性假设建模的。在此，我们将考虑前一种方法，因为它具有某些优势。首先，通过为每个试次（共M个）生成单独的模型，某些需要重复 "试次"的降噪算法可以应用于模型系数，即使它们是在不同的刺激物上建模的，例如联合去相关。其次，试次之间的不连续造成的伪迹不是问题。第三，交叉验证的效率要高得多，因为在少量的数据上训练模型并在不同的试验中取平均值要比串联大量的数据并对其进行训练快得多。

对于一个给定的试验，mTRFtrain函数同时训练所有的数据特征（例如，频带、响应通道）（见图1）。唯一的要求是刺激和反应数据具有相同的采样率（以Hz为单位），并且在时间上具有相同的长度。如前一节所述，向量和矩阵的组织方式应使所有的特征都以列向量排列。映射方向被指定为 "1"（前向建模）或"-1"（后向建模）。最小和最大的时滞以毫秒为单位输入，并根据输入的采样率转换为样本。算入额外的时滞（如预刺激时滞）来达到可视化的目的往往有效。

而且，人们还应该注意到在所产生的模型的任何一个极端都存在回归假象。然而，当为预测目的优化模型时，最好只使用刺激相关的时滞。岭回归中使用的滞后矩阵是由lagGen函数生成的。 如果用户指定向后映射，滞后就会自动反转，算法也会从公式（5）变为（9）。如果输入的刺激物是单变量（即一个矢量），算法将自动切换到方程6来使用更好的二次岭惩罚。最后必须指定岭参数λ。为了实现模型系数的可视化，λ可以根据经验选择最低值，这样任何增加都不会导致绘图估计值的明显改变。为了优化模型性能，应该采用更系统的方法，如交叉验证，照下文所述。

优化

刺激-反应映射的优化可以通过交叉验证实现（mTRFcrossval函数）。具体来说，目标是确定优化该映射的岭参数的值。这里，整个数据集被一起输入，M个刺激和M个反应矩阵被安排在两个胞元阵列中。不要求各个试验的时间长度相同（尽管出于优化的原因，这是最理想的情况）。优化交叉验证的另一个重要因素是输入和输出数据的标准化。通过对数据进行Z分数标准化，可以大大减少进行综合参数搜索所需的数值范围，使该过程更效率。衡量验证的岭值可以作为一个单个向量输入。其他参数的输入方式与mTRFtrain中相同。mTRFcrossval中的验证方法是留一法交叉验证，尽管这也可以被描述为M倍交叉验证。首先，对每一个指定的岭值的M个试次进行单独的模型拟合。试次被旋转M次后，每个试次会被排除或用作测试集，其余的M-1个试次被分配为训练集（见图1）。测试的实际模型通过对分配给每个训练集的单试次模型进行平均而得到。如前所述，这种方法比串联M-1个试次并对这些数据进行模型拟合更有效。然后，每个平均模型与相应测试集的数据进行卷积，以预测神经反应（前向建模）或刺激信号（后向建模）。这个过程对每个不同的岭值都要重复。模型的验证是通过比较预测的估计值和相应的原始数据来评估的。使用两个不同的验证指标：皮尔逊相关系数和均方误差。一旦得到验证指标，它们就应该在试次水平上取均值。这种方法是可取的，因为每个模型在理论上都需要相同的岭值来进行正则化，因为它们彼此共享M-2个试次的数据。这确保了模型在新数据上很好地泛化，并且不会对测试集过拟合。然而，如果所有的试次都是相同的长度，这种方法效果最好。最佳岭值被确定为产生平均最高的r值或最低的MSE分数。

测试

一旦用交叉验证法调整了模型参数，就可以用mTRFpredict函数在新的数据上测试优化后的模型。这可以在交叉验证程序之外的数据上进行（更好），也可以在交叉验证所用的同一测试数据上进行（图1）。如前所述，由于上述交叉验证程序采取的是跨试次的验证指标的平均值，所以模型不偏向用于交叉验证的测试数据。因此，基于这些数据报告模型的性能是合理的，因为在新的数据上测试很可能会产生相同的结果。虽然mTRFpredict函数输出的性能指标与mTRFcrossval相同，它还输出预测信号供进一步评估。当预测像EEG这样的多变量信号时，每个特征（即EEG通道）都会计算出一个性能指标，允许用户根据他们认为最相关的特征来评价模型。

图1. mTRF工具箱实现的前向和后向建模方法示意图。

刺激重建（即后向建模）可用于从记录的神经反应数据中解码特定的刺激特征，以估计这一信息在大脑中的编码的准确性。时序响应函数估计（即正向建模）可以用类似的方式来预测对新刺激的神经反应，也可以详细检验刺激特征是如何在大脑中编码的。

实例

本节介绍的例子使用了一项已发表的研究的数据，该研究测量了被试对自然、连续语音的EEG反应。被试听的是由一位男性演讲者用美式英语朗读的经典有声小说。音频分为28段（每段持续时间为155秒），其中5段的子集被用于本章的例子中。EEG数据是用一个128通道的ActiveTwo系统（BioSemi）记录的，采样率为512Hz。离线时，数据在1到15Hz之间进行数字滤波，降采样到128Hz，并以左右乳突为重参考。记录的128个通道中只有32个在分析中使用，但关键是均匀地分布在整个头部。进一步的细节可以在原始研究中找到。

本节详细介绍了几个例子，展示了如何利用mTRF工具箱以各种不同的方式将神经数据与感觉刺激联系起来。这些例子包括：

1.单变量TRF估计

2.优化和预测

3.多变量TRF分析

4.刺激物重建

5.多模态TRF估计

6.TRF与交叉相关

虽然这些例子都与语音刺激期间收集的EEG数据有关，但如前所述，这些方法都可用于其他类型的感觉刺激。

单变量TRF估计

这里的目的是估计时序响应函数，它将语音包络的单变量表示映射到每个通道记录的脑电信号上。语音信号的宽带包络（图2A）的计算方法如公式11所示：

其中xa(t)是由x(t)原始语音和其希尔伯特变换预测x(t)之和得到的复分析信号。包络被定义为xa(t)的绝对值。然后，在应用零相偏移滤波器后，将其降采样至与EEG数据相同的采样率。TRF是在-150和450毫秒的滞后期计算的，允许在两端有额外的50毫秒用于回归伪影。对五次试次中的每一次都单独计算估计值，然后取平均值。岭参数是根据经验选择的，以保持成分振幅。

首先通过计算32个通道的TRF方差来估计全局场功率（GFP）的大小（图2B）。GFP是对整个头皮在每个时滞的反应强度的参照独立测量。根据GFP测量的时间曲线，在50、80和140毫秒处有三个明显的TRF成分。图2C显示了这些成分中后两个的头皮拓扑图。它们的潜伏期和极性与典型的（乳突参考）听觉诱发反应的N1和P2成分相似。N1-P2复合体的地形图表明，这两个成分在中前部的FCz位置最强。图2D显示了在FCz通道计算的平均TRF和作为比较的在枕部位置Oz测得的TRF。

图2.单变量TRF估计。

（A）一个30秒的宽频带语音包络段.（B）在每个时滞测量的全局场功率.（C）发生在80毫秒和140毫秒的主要TRF成分的头皮地形图。黑色标记表示前中部通道FCz和枕部通道Oz的位置.（D）FCz（蓝色）和Oz（红色）的平均TRFs

优化与预测

这里的目的是使用TRF模型来预测新数据的EEG响应。这一次，我们用一种更系统的方法进行模型参数的调整，即前面描述的交叉验证程序的方法。具体来说，在每个独立的试次中，对一定范围的岭值计算TRF，（λ=20，22，...，220）。对于每个岭值，TRFs在每个四试次组合中取平均值，并用于预测第五个试次的EEG。在这里，建模数据在0至200毫秒之间的时滞，因为这些滞后反映了全局TRF反应中的大部分信息（图2B）。纳入额外的时滞（刺激前或刺激后）并没有改善模型的性能。

图3A显示了基于原始和预测的EEG反应之间的相关系数（Pearson's r）的交叉验证结果。重要的是，r值在五次试验中取平均值，以防止模型对测试数据的过拟合。r值在32个通道上平均，以便全面优化模型性能。另外，我们也可以只在指定的最高百分位数内或基于特定位置的通道上取平均值。图3B显示了基于平均平方误差的交叉验证的结果。使用同样的平均过程来确定最佳岭值。岭值的选择是为了使原始EEG和预测EEG数据之间的相关最大化。请注意，使用MSE作为交叉验证的标准也会产生同样的结果。图3C显示了使用优化的TRF模型在每个通道获得的相关系数。Pearson'sr的地形分布与TRF主导成分的分布非常相似（图2C）。事实上，该模型在反应最强的通道上表现最好。图3D显示了在FCz的脑电图反应的2s片段和优化的TRF模型所预测的相应估计。

图3.为EEG预测优化TRFs.

(A) 基于原始和预测的EEG响应之间的相关性的模型交叉验证（Pearson'sr在各通道和试次中平均）。被填充的标记表示最高的r值，即最佳岭值。

(B) 基于平均平方误差（MSE）的交叉验证。最佳岭值由最低的MSE分数来确定。

(C) 优化的TRF模型的测试显示了每个通道的相关系数。黑色标记表示通道FCz的位置。

(D) FCz处的脑电图反应的2秒分段（蓝色）和优化的TRF模型预测的相应估计（红色）。

多变量TRF分析

这里的目的是估计语音的多变量（频谱图）表征的TRF，即mTRF。频谱图表征（图4A）是根据Greenwoods方程，首先将语音刺激过滤成250和8kHz之间的16个对数间隔的频段而得到的。以对数方式滤波是试图模拟听觉过程进行的频率分析。每个频段上的能量都是用希尔伯特变换来计算的（公式11）。

mTRFs在-150-450ms的滞后期之间计算，模型参数根据经验进行调整。图4B显示了250-8000Hz之间所有频段的FCz通道的mTRF响应。图4B的视觉检查表明，主导的N1成分TRF和P2成分TRF几乎在每一个频段都编码了语音信息，直到6kHz，这是语音信号中包含的大部分信息（图4A）。对各频段的mTRF进行平均，将得到一个单变量的TRF测量值，它与使用宽带语音包络计算的TRF非常接近（图2D）。为了用mTRF模型预测EEG响应，采用了与之前相同的方法。尽管交叉验证产生的结果（图4C,D）与单变量TRF方法（图3A,B）相似，但mTRF方法似乎对岭值的变化更敏感。进一步调查发现，这不能归因于在单变量和多变量分析中使用不同的正则化惩罚因子。尽管如此，在整个头皮上，优化的mTRF模型的性能与单变量TRF模型相近（图4E,F）。虽然已经证明mTRF模型在预测EEG反应方面优于单变量TRF模型，但必须考虑到多变量TRF分析对正则化更敏感，并可能涉及相当多的计算。

图4.多变量TRF估计和EEG预测.

(A) 语音频谱的30s分段。

(B) 通道FCz的平均mTRF。

(C) 基于原始和预测的EEG反应之间的相关性的模型交叉验证（跨通道和试次平均的Pearson'sr）。被填充的标记表示最高的r值，即最佳岭值。

(D) 基于平均平方误差（MSE）的交叉验证。最佳岭值由最低的MSE分数来确定。

(E) 优化的mTRF模型的测试显示了每个通道的相关系数。黑色标记表示通道FCz的位置。

(F）FCz处的EEG响应的2秒片段（蓝色）和由优化的TRF模型预测的相应估计值（红色）。

刺激重建

刺激重建的目的是产生一个解码器，在后向（即从EEG到刺激）对数据进行建模，并使用它来估计重建单变量刺激输入。这种方法相对于正向建模技术的优点在介绍和后向模型部分有概述。使用与TRF模型描述的相同的交叉验证技术调整模型参数。具体来说，在0-200ms的时滞下，对相同范围的岭值计算解码器，（λ=20，22，...，220）这里的区别是，EEG被视为输入，刺激被视为输出，而且滞后的方向相反，即-200到0。图5A显示了使用原始和重建的语音包络之间的相关系数衡量交叉验证的结果，而图5B表示基于MSE的模型脊参数的验证。为防止对测试数据的过拟合，这两个指标在不同的试次中都是平均的。32个EEG通道都纳入在模型验证程序中以优化性能。后向建模方法比前向建模的优势体现在以相关值为指标的残差更小。这主要是由于解码器可以同时利用整个头部的信息（即在多变量意义上）来确定语音估计，而在正向建模时，预测的EEG估计是基于刺激和该特定通道的EEG反应之间的单个单变量映射。此外，前向建模方法的预测是在EEG域中评估的，低信噪比对预测准确性有负面影响。相比之下，后向建模方法将这些措施的估计转移到由实验者定义的刺激域。因此，在低信噪比的刺激物（如语音包络）的情况下，后向建模的拟合质量可能会更高。

虽然解码器通道权重在神经生理学意义上不容易解释，但它们的权重反映了对重建刺激信号贡献最大的通道。图5C显示了110至130毫秒之间的平均解码器权重（由GFP来表示权重最大的区域）。

与TRF拓扑图（图2C）相比，模型权重的分布更加局部。因为解码器不需要像TRF那样在整个头皮的每个通道上编码信息，它可以有选择地只加权对重建重要的通道，同时通过赋予较低的权重来忽略不相关的和噪音通道。图5D中可以看到一个重建的2s估计样本。除了必须对每个特征（例如频带）分别进行模型性能评估，或通过对不同特征进行平均，多变量刺激的刺激重建以同样的方式进行。以前的研究描述了一个能够对后向模型系数进行神经生理学解释的程序。具体来说，他们提出了一种确定性的方法，将之前拟合的线性后向模型系数转化为线性前向模型系数。这一过程使神经生理学能够解释线性后向模型的参数，否则可能存在误导和错误。mTRFtransform函数专门为使用岭回归技术导出的后向模型实现了这一程序（例如图5E,F）。

图5.刺激重建。

(A) 基于原始和重建的语音包络之间的相关性的模型交叉验证（Pearson's r跨试次的均值）。填充的标记表示最高的r值，即最优岭值。

(B) 基于平均平方误差（MSE）的交叉验证。最优岭值是由最低的MSE分数确定的。

(C) 解码器通道权重在110和130毫秒之间的平均时滞。

(D) 原始语音包络（蓝色）和由优化解码器重建的相应估计（红色）的2秒分段。

(E）使用Haufe等人描述的逆程序将解码器通道权重转换到正向模型空间。黑色标记表示中前部通道FCz和枕部通道Oz的位置。

(F) 通过将解码器权重倒置到FCz（蓝色）和Oz（红色）的正向模型域而得到的时序响应函数。

多模态TRF估算

除了提取对单模态刺激的神经反应外，TRF方法还可用于分离多模态信号（或同一模态中的多个信号，如鸡尾酒会场景）。这可以用在自然视听语音期间记录的EEG来证明。这里介绍的数据发表在一项研究中，该研究调查了视觉语音对听觉语音皮质表征的影响。被试听了15分钟的自然视听语音，这里使用了其中7分钟的子集。听觉刺激的特征和以前一样是宽带包络（图6A），而视觉刺激的特征是计算视频的帧与帧之间的运动（图6B）。对于每一帧，使用 "Adaptive Rood Pattern Search "block匹配算法计算出运动向量矩阵。通过计算每一帧的所有运动矢量长度的总和，得到全局运动流的测量值。然后将其从RGB表示转换为相对亮度，并从30到128赫兹升采样，以匹配脑电数据的采样速率。

图6C显示了声学包络映射到EEG时在通道Fz和Oz的TRF，而图6D显示了视觉运动信号映射到相同EEG数据时在相同通道的TRF。当使用听觉信号时，Fz的TRF振幅比Oz大得多；反之，视觉信号则是Oz的振幅更高。这也可以从地形图中看出，语音包络TRF在额部头皮上有优势反应（图6E），运动TRF在枕部头皮上有优势反应（图6F）。虽尽管在两种情况下分析了相同的脑电图数据，但通过简单地映射每种感觉模态的特征，可以提取来自不同感觉皮层区域的响应。为了测量多感官整合，函数可以用来拟合一个基于单感官模型系数代数和的 "additive model（加性模型）"。在多感觉神经反应数据上对加性模型进行测试，并将其性能与多感觉模型进行比较，获得客观的衡量标准。

图6. 多模态TRF估计。

(A) 宽带语音包络的30s分段。

(B) 相应的帧与帧之间的视觉运动的30s分段。

(C) Fz（蓝色）和Oz（红色）的平均包络TRFs。

(D) 在Fz（蓝色）和Oz（红色）的平均运动TRFs。

(E) 发生在78毫秒的主导包络TRF成分的头皮地形。

(F)发生在117毫秒的主导运动TRF成分的头皮地形。

TRF与交叉相关

如前所述，LTI（线性时不变（Linear time-invariant，LTI））系统的脉冲响应可以很容易地通过输入和输出的简单的相关来近似。虽然这种方法比使用归一化反向相关或岭回归等技术更直接，但它只适用于符合随机过程的输入信号。

为了证明这一点，我们用语音和白噪声作为刺激输入信号，对每一种方法进行了比较。语音数据与前面例子中的数据相同，使用的非语音数据是在另一项研究中发表的。该研究调查了用于估计听觉系统对高斯白噪响应的TRF函数。受试者听了10段120秒的不间断噪声刺激，本例中使用了其中的6段子集。刺激是高斯宽带噪声，能量限制在0-22.05kHz的带宽内，用0-30Hz的功率均匀的高斯噪声信号进行调节。为了说明听觉刺激强度感知的对数性质，这些调制信号的值X随后使用公式12的指数关系被映射到音频刺激的振幅x′。

EEG的记录和处理与前面例子中描述的程序完全相同。实验中使用的语音和噪音刺激分别显示在图7A、B中。每个刺激的自相关显示，语音刺激在多个时滞区间自相关（图7C），而噪声刺激只在零时滞自相关（图7D）。图7F显示了使用TRF方法和交叉相关（XCOR）方法在通道FCz计算的白噪刺激的脉冲响应。视觉检查表明，交叉相关法和TRF法对系统响应函数的估计大致相同。然而，对于语音刺激来说，与TRF方法相比，交叉相关方法会导致脉冲响应估计的时间污染（图7E）。这是因为刺激动态在多个重叠的时滞映射到EEG信号上。这证明了TRF技术在描述感觉系统对缓慢调制的自然刺激(如人类语言)的响应方面的效用。

图7. 时间响应函数（TRF）和交叉相关（XCOR）方法的比较.

(A) 宽带语音包络的30秒分段。

(B) 30秒的振幅调节噪声分段。

(C) 语音包络的自相关。

(D) 噪声信号的自相关。

(E) 使用TRF方法（蓝色）和交叉相关方法（红色）估计的通道FCz处的语音脉冲响应。

(F) 使用TRF方法（蓝色）和交叉相关方法（红色）估计的通道FCz上白噪的脉冲响应。

讨论

在此，我们描述了一个新的基于MATLAB的SI（系统识别）工具箱，用于对神经信号和自然、连续的刺激间的关系进行建模。上述例子展示了这个多功能工具箱如何应用于单变量和多变量的数据集，以及单感官和多感官的数据集。重要的是，它还可以用于向前和向后向的映射，分别进行反应函数估计和刺激重建，提供互补的分析技术。

应用

mTRF工具箱在感觉神经科学方面有许多应用，其中最重要的是研究自然语音在人脑中是如何处理的。前向TRF方法以前曾被用来证明神经系统如何使用颅内和非侵入性记录技术，以精确的时间分辨率提取人类对不间断语音的神经响应。使用这种方法的后续研究产生了几个关键的发现，涉及到大脑如何在鸡尾酒会场景中选择性地关注单个语音流，以及时频和语音信息如何在听觉皮质活动中表征。

该工具箱的其他应用包括使用后向和前向模型来研究视听语音处理和视觉语音处理，即语音阅读。替代SI技术（最终产生相同的解决方案）也被用来研究听觉场景分析、噪音中的语音、语音的显性和隐性皮质表征以及从颅内记录中提取语音的详细谱图重建。

除了研究语音，前向TRF方法还被应用于视觉研究，研究人脑如何处理对比度随时间变化的刺激。这种特殊的方法也被用于临床研究，研究自闭症谱系障碍儿童和精神分裂症成人患者的视觉处理缺陷。最近，它被修改为研究大脑如何处理更自然的视觉刺激，如连贯的运动。除了描述视觉刺激特征和脑电图记录之间的映射之外，研究人员最近使用相同的正则化线性回归方法从表面EMG（肌电信号）重建了手指运动，进一步证明了该技术的多功能性。

注意事项

反向相关法的线性假设对其解释有影响。刺激特征和神经反应振幅之间的线性关系的假设可能导致反应度量反映了皮质细胞子集的前馈活动。因此，这种方法可能对以非线性方式与刺激相关的皮质反应不敏感，包括横向和反馈贡献，这可能对研究高阶认知过程的影响存在干扰。这与区分无数的前馈、侧向和反馈贡献给锁时平均ERP所涉及的挑战形成鲜明对比。事实上，需要解决线性假设，以便准确地描述神经元群体以非线性的方式对复杂刺激做出反应。也就是说，以前的一项研究实现了线性TRF方法的二次扩展，用于对对比度刺激的视觉反应建模，但没有发现相对于线性模型性能的任何明显改善。

然而，随后的研究将相同的二次模型应用于听觉系统，确实证明了对白噪刺激的模型性能的改善。使用支持向量回归机器学习技术，也曾探索过TRF模型向高阶的扩展，但同样，只取得了十分有限的改进。虽然这种非线性回归技术可以使模型的性能略有提高，但在性能和计算时间之间有一个相当大的权衡，这往往使它们不切实际。然而，对于集群数据（如脑电图），非线性模型的性能只比线性模型好一点；但非线性模型对于单个单位数据（如ECoG）的建模似乎更有利，这一事实可能意味着关于EEG记录的一些基本性质。由于容积传导的空间污染效应，每个EEG电极都能检测到来自大型皮质群体（107-109个神经元）的神经活动。因此，在最大的神经群体中普遍存在的激活模式将对头皮记录的信号贡献最大。由于神经元间非线性反应的多样性，这种活动很可能被编码在小的神经元亚种群中，而线性反应则可能被编码在更宏观的水平上。因此，容积传导的影响可能导致许多非线性活动在由此产生的EEG记录中被掩盖。事实上，如果是这样的话，就可以解释为什么线性回归技术在脑电图反应建模方面的表现要比非线性技术好。为了支持这一概念，其他EEG/MEG建模算法，（如SPoC，它将神经振荡的振幅与刺激特征或行为相关参数联系起来）专门使用了线性模型，因为已知这种振荡的叠加是线性的和瞬时的。