静息态fMRI中的非线性功能网络连接

在这项工作中，我们关注功能网络中的显式非线性关系。我们介绍了一种使用归一化互信息(NMI)计算不同大脑区域之间非线性关系的技术。我们使用模拟数据演示了我们提出的方法，然后将其应用到Damaraju等人先前研究过的数据集。静息状态fMRI数据包括151名精神分裂症患者和163名年龄和性别匹配的健康对照组。我们首先使用组独立成分分析(ICA)对这些数据进行分解，得到47个功能相关的内在连通性网络。我们的分析显示，大脑功能网络之间存在模块化的非线性关系，在感觉和视觉皮层尤其明显。有趣的是，模块化看起来既有意义又与线性方法所揭示的不同。分组分析发现，精神分裂症患者与健康对照组在显式非线性功能网络连接(FNC)方面存在显著差异，特别是在视觉皮层，在大多数情况下，对照组表现出更多的非线性(即，去掉线性关系的时间过程之间更高的归一化互信息)。某些域，包括皮层下和听觉，显示出相对较少的非线性FNC(即较低的归一化互信息)，而视觉域和其他域之间的联系显示出实质性的非线性和模块化特性的证据。总之，这些结果表明，量化功能连接的非线性依赖性可能通过揭示通常被忽略的相关变化，为研究大脑功能提供一个补充和潜在的重要工具。除此之外，我们提出了一种方法，在增强的方法中捕捉线性和非线性效应。与标准线性方法相比，这种方法增加了对群体差异的敏感性，代价是无法分离线性和非线性效应。

1.简介
在过去的几十年里，功能性磁共振成像(fMRI)已经成为了解大脑功能最广泛的方法之一。在这一领域，功能性网络连接(FNC)已被广泛用于分析不同大脑区域之间的关系。大多数功能连通性研究集中在时间过程之间的线性关系。对功能磁共振数据网络模型的研究表明，线性相关是一种成功的统计工具，可以识别功能磁共振时间过程之间的关系，这也很容易计算和解释该领域的正相关和负相关(例如，默认模式网络往往与其他网络反相关)。然而，一些研究表明，一些大脑活动表现出非线性动态行为。其他研究讨论了fMRI数据中血流动力学响应的非线性效应，重要的是，它也可以随时间(和地点)和不同个体的变化而变化。考虑到这几个非线性效应的例子，很可能，甚至是预期，不同的大脑区域可能在某种程度上是非线性相关的，这是传统线性分析所忽略的。在当前的研究中，我们感兴趣的是评估在功能连接环境中大脑区域之间明显的非线性关系(即去除线性关系后)存在的程度，并确定显著的依赖性。据我们所知，在功能连接中明确研究非线性关系的工作很少。尽管线性相关系数被广泛应用于该领域，但它只测量线性关系而忽略了非线性贡献。其他可以评估非线性的高阶统计工具对功能磁共振时间过程关系不太敏感。因此，我们提出了一个新的统计工具来衡量显式非线性依赖。我们关注的是互信息(MI)的规范化版本，这是一种信息理论方法，具有测量线性和非线性依赖关系的能力。早期的工作评估MI以捕捉更一般的关系。最近，人们探索了包括MI在内的功能性连接的替代指标。然而，据我们所知，我们是第一个评估大脑网络之间显式非线性关系的小组，评估其独特的方面相对于线性关系。总之，本文的贡献如下。我们开发了一种方法，探索功能大脑网络之间的非线性依赖关系。我们的方法通过使用回归方案去除线性关系后，利用剩余依赖性之间的互信息计算非线性依赖性。为了评估非线性关系是否具有潜在的意义，我们首先关注所得到的FNC矩阵是否表现出与功能整合一致的模块关系。其次，我们评估了从精神分裂症患者和健康对照组收集的静息fMRI数据集中，非线性FNC是否显示出有意义的组间差异。最后，我们提出了一种统计方法，它提供了一种选择，在保留线性解释的同时，也考虑了额外的非线性依赖性。

2.材料和方法

2.1 被试和预处理

在本研究中，我们使用了fBIRN数据集，该数据集之前已进行过分析。最终策划的数据集包括163名健康参与者(平均年龄36.9岁，117名男性;46名女性)和151名年龄和性别匹配的精神分裂症患者(平均年龄37.8岁;男性114人，女性37人)。闭眼静息状态fMRI数据收集于美国各地的七个地点。扫描前根据相关机构的内部审查委员会取得所有受试者的知情同意。其中一个部位的成像数据由通用电气Discovery MR750扫描仪捕获，其余六个部位则由西门子Tim Trio系统收集。使用AFNI、SPM和GIFT等工具对数据进行预处理。使用SPM中的INRIAlign工具箱进行刚体运动校正，以校正头部运动。为了去除异常值，执行AFNI3s 3dDespike算法。然后将fMRI数据重采样为3 mm³的各向同性体素。然后采用AFNI3s BlurToFWHM算法对数据进行半最大平滑至6 mm，并对每个体素时间过程进行方差归一化处理。运动较大的受试者被排除在分析之外，以减轻运动对策展过程的影响。

2.2 后处理

使用组级空间ICA的GIFT (http://trendscenter.org/software/gift)实现来估计100个功能网络作为ICA成分。首先使用个体特定的数据简化步骤，利用主成分分析将162个时间点数据简化为100个最大变异性的方向。其次，信息最大化方法用于从组PCA约简矩阵中估计100个最大独立分量。为了保证估计的稳定性，对ICA算法进行多次重复，并选取最中心的运行作为代表。最后，将聚合空间映射估计为组成簇的模式。受试者特定空间图(SMs)和时间过程(TC)是通过时空回归重建方法获得的。为了标记各成分，我们获得了每个特定空间图的峰值激活簇区域。ICA处理后，为获得峰值激活簇的区域，对所有受试者的每个SM进行一个样本t检验图，然后进行阈值;并计算了相应TC的平均功率谱。如果其峰值激活簇位于灰质范围内，且与已知的血管、易感性、心室和头部运动对应的边缘区域重叠较少，则可以识别作为内在连接网络(ICNs)的成分集。这导致100个独立成分中有47个ICN。运行20次ICASSO后，除一个ICNs外，所有ICNs的集群稳定性/质量指数都非常高。考虑预估的受试者运动参数，对TC进行去趋势和正交后，采用AFNI3s 3dDespike算法检测尖峰，并用三阶样条拟合值代替尖峰。更多细节请参见Allen等人(2012)和Damaraju等人(2014)。经过处理，fBIRN数据集得到159个时间点47个成分的矩阵，314名受试者，其中163名对照组受试者和151名SZ受试者。

2.3 互信息方法

虽然线性相关是描述相关性最广泛使用的度量，但它可能完全忽略非线性相关性。说明这种不足的一个例子是Anscombe的四重奏，它显示了四个不同的非随机数据点的图具有相同的相关系数，尽管它们的依赖结构截然不同。为了测量一对TC之间的显式非线性关系，本研究采用的方法是去除线性相关并计算残差相关性。

时间过程x,y的皮尔逊积矩相关系数ρ为：

SxSy分别是样本标准误差，Cov(x,y)是样本协方差。

相关系数主要衡量两个分布之间的线性相关性。然而，非线性的相关性并没有被捕捉到。最近的统计方法已被提出，以衡量相关性而不低估非线性依赖性。其中一种方法是互信息(MI)，它测量线性和非线性依赖关系。这个公式决定了MI的值为：

H(x),H(y)是边缘熵，H（x,y）是联合熵。然而，MI单元不是标准化的，这使得它很难在个体和数据集之间进行比较。Kvalseth(2017)讨论了NMI的一些规范化因素，包括多种选项，如：min(H(x),H(y)),H(x)+H(y),max((H(x),H(y))等。在本研究中，我们使用后者作为(Horibe, 1985)证明它是一个(归一化的)相似度度量。归一化MI (NMI)公式为：

在这项工作中，我们的目标是只计算依赖关系的非线性成分。为此，我们在去除线性相关性后测量数据的互信息依赖性。

2.4 模拟实验

我们将所提出的方法应用于模拟数据，以说明其应用。在这个实验中，我们从一个1000x1大小的x向量开始，它的分量来自[0 1]上的随机均匀分布。接下来，我们形成三个向量y1,y2,y3，每一个向量都与x有特殊关系。三种不同类型的关系如下：情形一，x与y1与向量有纯线性关系。情形二，x与y2我们定义为二次关系。情形三，x与y3向量具有线性和非线性的组合依赖。我们还添加了零均值高斯噪声(图1)。

图1 两个向量之间线性和非线性关系的三个仿真案例。向量的分量x随机地来源于[0 1]上的均匀分布。
我们使用皮尔逊相关和标准化互信息方法测量x-y1,x-y2,x-y3的关系。皮尔逊相关性取值从-1到1。简单地说，1表示完全线性负相关，1表示完全线性正相关。我们在这项工作中使用的归一化互信息的范围是[0,1]。NMI = 0表示无依赖性，随着两个分布的依赖性增加，NMI值上升到最大值1。在计算相关性和NMI之前，我们实现了前面解释的过程来去除线性相关性。接下来，我们计算每一对的Pearson相关性和归一化互信息，如表1所示。

表1

2.5 量化fMRI中的非线性连接

每个科目有47个ICA时间过程，长度为159。对于每一对x,y，我们计算传统的FNC(即所有对和之间的线性相关)。接下来，计算314个受试者的平均FNC矩阵(图2a)。

图2
然后我们拟合一个线性模型来估计x和y之间的线性相关。在此之后，我们删除线性效应，通过更新为来研究剩余的依赖性。接下来，我们通过NMI计算功能磁共振数据中功能网络成分之间的剩余依赖关系。这将为每个个体生成一个47 * 47的矩阵，其中条目中的值显示了由NMI方法计算出来的x,y的非线性依赖关系。在此之后，我们计算了总体受试者的平均值(图2b)。然后，为了评估NMI在整个大脑中是否显示出显著的变化，我们进行了t检验，将每个单元的平均值与最小值的平均值进行比较，以确定给定单元的平均值显著大于最小平均单元的单元。此外，我们使用随机矩阵分析方法来检验得到的显式非线性依赖矩阵的模块化。根据图2中线性FNC和显式非线性FNC的结果，我们选择了两个极端单元对(成分23和成分38)，它们的线性相关性低，显式非线性相关性高;对(成分2和成分3)，它们的线性相关性高，显式非线性相关性低。对于每一对，TC及其频谱如图3所示。

图3 两对ICN分量的线性和显非线性关系。左边的图与配对(23,38)有关，右边的图与配对(2,3)有关。

我们还比较了精神分裂症患者和对照组之间的非线性依赖关系。每一组的线性效应被取消，并计算所有受试者的平均NMI。然后，我们实施了t检验以确定显著的组间差异。在错误发现率(FDR)校正中，所有p值均采用Benjamini Hochberg校正法进行校正，阈值为p<0.05。

2.6 提高的方法
虽然我们强调非线性中包含的独特信息，但未来的研究可能希望同时利用线性和非线性信息。皮尔逊相关由于计算简单，能更灵敏地捕捉功能磁共振成像时间过程间的主要线性依赖关系，被广泛应用于功能连通性研究。此外，线性相关的符号揭示了有关大脑功能的有意义的信息。我们提出了一种方法，在保留线性解释的同时，也考虑了额外的非线性依赖性。这种增强的方法是皮尔逊相关和修改互信息的组合，用于量化2.3节中所述的非线性依赖关系。我们将增强的方法定义为：

通过这种技术，我们使用非线性信息来增强皮尔逊相关方向的线性效应，反映出更强的依赖性。因此，线性和非线性关系都被考虑了，线性效应的方向(在非线性情况下没有很好地定义)被保留了。我们提出的非线性和增强方法允许多种用途。我们可以只关注非线性效应，就像我们在本文中展示的那样，它本身可能很有趣。第二，我们可以集中精力捕捉线性和非线性的影响。这似乎增加了对标准线性分析之外的群体差异的敏感性，尽管它没有分别考虑线性和非线性效应。注意NMI捕捉线性和非线性效应;不过，NMI对线性效应的敏感性有所降低。此外，NMI不能捕获方向关系，因为NMI总是正的。改进的方法保留了线性关系的敏感性和方向性，这也为非线性关系提供了一些敏感性。我们评估了精神分裂症患者和健康对照组分的线性相关(Pearson相关)、非线性依赖(2.3节介绍的修改的互信息)以及线性和非线性依赖(增强)。然后分别对每种方法进行t检验，并采用Benjamini Hochberg校正方法对p值进行校正，阈值p<0.05。

2.7 联合分布
为了进一步可视化所识别的非线性关系，我们选择了在HC-SZ的非线性相关性的群体差异t检验中p值最显著的5个成分对。然后我们构建了每对时间过程中关节分布的差异，比较患者和对照组。

3.结果

3.1 模拟实验

三类依赖关系:线性、非线性以及线性和非线性的组合被检验。表1测量并报告了每个病例去除线性相关性前后的Pearson相关性和互信息。皮尔逊相关的范围为-1 ~ +1，独立分布的归一化互信息为0，完美依赖性为1。在情况I中，两个分布只有线性相关，在去除线性效应之前，皮尔逊相关接近于1。去除后，皮尔逊互信息和归一化互信息都接近于零。在情况II中，两个分布具有二次关系，皮尔逊相关显示出低但非零相关。相比之下，归一化互信息计算表明两个分布之间有相当的相关性。在去除线性效应后，Person相关性有效为零，而在去除线性效应前后的互信息几乎相同。在情形III中，两个分布之间既有线性关系又有非线性关系，在去除线性效应之前，皮尔逊相关是显著的。消去线性相关后，归一化互信息消失，而去线性相关后，归一化互信息仅略有下降。这简要地说明了皮尔逊相关并不捕获纯粹的非线性依赖，而互信息同时考虑线性和非线性依赖。类似地，如果我们去除线性效应，相关性将变为零，而互信息将捕获两个分布之间的真正残余非线性依赖关系。

3.2 fMRI数据

我们测量了静息功能磁共振数据估计的47个ICN时间TC之间的线性相关性，结果见图2。图2a是314名受试者的平均FNC。

对于相同的成分，使用我们提出的NMI方法计算所有个体的非线性依赖的平均值(图2b)。对随机性的评估证实了显式非线性FNC表明相对于随机矩阵具有高度显著的模块化。然后，对一个样本进行t检验，以确定平均值与最小平均值有显著差异的对。FDR阈值化后，大多数对被检测到显著大于最小平均值。这向我们展示了NMI差异在FNC矩阵中的显著不同程度，并展示了视觉(VIS)和躯体运动(SM)成分与其他成分之间的模块化非线性依赖关系。有趣的是，尽管皮质下(SC)和听觉(AUD)之间有很高的线性相关，但观察到的非线性依赖率很低。对一个受试者，绘制出配对(23,38)和配对(2,3)的时间轨迹(图3a)，并用快速傅里叶变换(FFT)的振幅计算出它们的频谱(图3b)。在所有受试者中，配对(23,38)显示出相对较高的显式非线性相关性(0.0236)和较低的线性相关性(0.0127)。(2,3)对总体被试的平均线性相关性较高(0.2194)，但显式非线性相关性较低(0.0208)。

接下来，使用我们的NMI方法评估健康对照组和精神分裂症患者之间非线性依赖性的变化，并使用双样本t检验进行比较。在图4a中，下面的三角形显示阈值之前乘以FDR之前的初始p值。阈值化后，上三角形与下三角形相同。用不同颜色显示的条目表明组间在非线性依赖性上有显著差异的对。这一结果表明，精神分裂症(SZ)患者相对于健康对照组(HC)，视觉(VIS)成分对听觉(AUD)、视觉(VIS)、躯体运动(SM)、认知控制(CC)和默认模式(DM)等其他成分的非线性依赖性有显著差异。面板b，连接图，是t检验结果的另一种表示。如果HC和SZ之间的非线性依赖性差异显著，则空间ICN连接一条线(黄色表示HC>SZ和蓝色为SZ>HC)。

图4 (a)上三角:去除线性相关后NMI中的群差(HC-SZ)。p值由FDR调整，阈值(p<0.05)。数值被绘制为-log10§ x sign(t-statistics)。下三角:与上三角相同，除了p值没有阈值化，而且这些值乘以FDR之前的初始p值。(b)连接图显示HC-SZ中具有显著不同非线性关系的成分是连接的。除一种情况外，所有情况的结果都显示，对照组的非线性明显更多，大部分与视觉领域有关。
3.3 增强方法

健康对照和精神分裂症患者中各成分之间的相关性有三种评估方法:(1)皮尔逊相关性，其中只强调线性相关性;(2)修改互信息，如2.3节所述，只量化非线性相关性;(3)增强方法，如2.6节所述，通过获取非线性相关性来增强线性相关性。接下来，在每种方法中，应用t检验来比较两组之间的差异。FDR调整的p值为阈值(p<0.05)。Pearson相关法具有显著差异的对数为530对，修正互信息法为17对，提升法为537对。有5对线性相关无显著差异，但非线性依赖显著。注意，线性分析和显式非线性分析是相互补充的。显式非线性NMI方法给出了17个显著对，而不管它们的线性依赖是否显著。也就是说，线性相关方法完全忽略了非线性变化中的这些信息和差异，重要的是，非线性群体差异的全脑模式是高度模块化和结构化的。增强方法提供了一个选项来捕获依赖性的线性和非线性成分。结果还表明，增强方法识别了线性或全NMI方法都无法捕获的关系。

3.4 联合分布

我们还对时间TC对之间的关系的可视化感兴趣，这显示了患者和对照组之间的非线性差异，图5显示了健康对照组和精神分裂症患者之间关节分布的差异，这五对显示了最大的组差异。数值从左到右，从下到上递增。面板a显示了第26和第6分量联合分布的差异。第26个成分属于认知控制(CC)。第6个成分属于听觉(AUD)域。面板b是第44和第13分量的联合分布差。第44个成分为DM默认模式，第13个成分为VIS域。面板c表示第7和第16分量的联合分布差异。第7部分是听觉(AUD)，第16部分是视觉(VIS)。图d显示了第14和第7分量联合分布的差异，图e显示了第14和第17分量联合分布。第14和第17分量都属于视觉(VIS)域。面板b、d和e有一些相似之处，包括两个成分之间的负相关关系。

图5 5个对(26,6)，(44,13)，(7,16)，(14,7)，(14,17)的HC和SZ关节分布的差异在对的非线性依赖中表现出最大的群差异。每个分布中的值从左到右，从下到上递增。第26部分属于认知控制(CC)。第6和第7个成分是听觉(AUD)。44为默认模式(DM)。所有第13、14、16和17个分量都在视觉(VIS)域中。我们观察到关节分布中一些有趣的差异，患者通常表现出不同的一个网络较高的活性和配对网络较低的活性。

在图5a中，我们可以观察到，不管第26个分量的值是多少，控制组相对于SZ在第6个听觉分量的低水平上花费的时间更多。从图5b中可以看出，在默认模式网络(#44)和可视成分(#13)中，健康的成分都比SZ中更加活跃。在c组中，我们注意到与SZ组相比，健康对照组在听觉(#7)和视觉(#16)两个成分上的活动都更少。从图5d,e可以解释为，与患者相比，健康对照组的视觉(#14)部分相对于更靠后的视觉(#17)和听觉(#7)部分显示出更高的激活水平。

对于这5对，非线性FNC比SZ表现出超连通性，即两个分布的非线性部分在HC中具有更强的依赖性。然而，联合分布的群体差异图可以说明在HC中表现出高度依赖性的两对在SZ中如何不同地分布。面板b、d和e显示健康的成分大多处于较高的激活级别(红色拖到右上角，蓝色拖到左下方)。相反，面板c显示健康成分往往显示较低的激活水平(红色拖到左下方，蓝色拖到右上方)。

4. 讨论
   在这项初步工作中，我们强调了研究功能连通性中的非线性的好处。以往的功能连通性研究都是基于相关系数，只评估线性相关，从而忽略了非线性的贡献。我们建立了一种评估不同大脑区域之间显式非线性依赖的方法，首先消除线性依赖。我们首先在模拟数据上证明了我们的方法如预期的那样工作(图1)。随后评估了314名受试者的47个时间过程之间的非线性依赖关系(图2)。我们应用了类似的方法来估计精神分裂症患者大脑平均不同区域对功能连通性的非线性贡献有多不同(图4)。研究了HC-SZ中非线性依赖群差异最大的5对的联合分布(图5)。非线性依赖的可能原因有很多，包括(1)非线性血流动力学效应。对神经元活动、氧代谢和血流动力学反应之间关系的研究表明，在任务数据中，神经元活动和血流动力学反应大小之间的联系表现出线性和非线性效应。血流动力学响应随时间的变化可诱发非线性。血液动力学在不同的研究对象和群体之间也会有所不同。其他结果表明，神经活动的电生理测量与血流动力学反应之间存在很强的非线性关系。(2)受试者内部和组间血流量、血氧和血容量的差异。实验表明，获得性血管空间占用(VASO)、动脉自旋标记(ASL)灌注和BOLD信号对刺激持续时间的响应呈非线性。(3)主体运动。即使是轻微的头部姿势变化也可能导致大脑中相当复杂的空间场变化。虽然我们不能完全排除运动，但我们仔细整理了数据，以关注低运动的受试者，此外，两组之间没有显著的运动差异。

线性和非线性模块化中的不同模块化模式表明了非线性和线性关系的互补性。在未来的研究中利用这些差异可能很重要。我们的研究结果表明网络之间存在有趣的差异。例如，如图5b所示，在视觉(VIS)、躯体运动(SM)域和认知控制(CC)和默认模式(DM)域之间观察到显著的非线性依赖关系。听觉(AUD)网络在线性依赖方面表现出很强的差异(a)，但非线性不太明显，而视觉和感觉运动网络在域内表现出很强的非线性依赖(b)。此外，在皮层下(SC)和听觉(AUD)与其他成分之间观察到相对较低的非线性依赖率。

我们还发现患者和对照组之间的非线性关系有显著差异。SZ和HC之间的非线性FNC两两比较如图5a所示。在大多数情况下，对照显示出相对于患者较高的非线性依赖性，大多与视觉域有关。SZ患者和HC患者的视觉(VIS)成分之间以及视觉(VIS)成分与听觉(AUD)、躯体运动(SM)、认知控制(CC)和默认模式(DM)等其他成分之间的非线性关系存在显著差异。我们观察到，大多数患者/对照组的差异涉及视觉和听觉部分。鉴于现有证据表明，一些精神分裂症症状可能与视觉系统有关。话虽如此，视觉幻觉等视觉症状在SZ相当罕见，比听觉和触觉异常更罕见。此外，一些研究表明，先天失明可能可以防止精神分裂症的发展，其特征是不可避免的嘈杂的知觉输入，导致错误的推论。这些发现表明，当一个人从出生就看不见东西时，他们会更多地依赖其他感官。因此，由此产生的世界模型更能抵抗错误的解释。

虽然提出的结果显示关注功能磁共振成像数据中显式非线性依赖的潜在效用，但仍有许多工作要做。未来的工作应该集中于仔细评估非线性关系的可能来源。定量功能磁共振成像研究可用于分离血氧、血容量和血流量的非线性。此外，高场层特异性fMRI研究可用于评估输入和输出层的非线性。在未来的工作中，还可以评估各种生理变量(如呼吸、二氧化碳、心率和运动)的贡献。

此外，我们的结果提供了证据，表明功能磁共振成像时间过程之间有意义和显著的非线性依赖性。我们已经有证据表明FNC中有意义的(模块化的和不同群体的)超线性效应，这主要暗示了在精神分裂症中视觉皮层被破坏。我们提出了两种方法，专注于显式线性效应或增强方法，在一个度量中捕获线性和非线性效应。未来的工作应该研究非线性关系中包含的信息。它可以通过更快的采集进行研究，与多模态成像(如并发脑电图数据)相关联，并复制我们在这项工作中显示的结果。

参考文献：Nonlinear functional network connectivity in resting functional magnetic resonance imaging data