脑网络和拓扑神经科学的分步教程

文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。

导读

大脑是一个非常复杂的系统，它能够优化整合来自不同脑区的信息以便执行其功能。随着最近技术的进步，研究人员可以使用不同规模和多种形式的神经成像技术从大脑中收集大量的数据。随之而来的是对复杂分析工具的需求。网络神经科学领域一直在努力应对这些挑战，其中，图论技术是研究大脑网络的重要分支之一。最近，拓扑数据分析作为一种替代框架获得了更多的关注，它提供了一套超越成对连接的指标，并具有更好的抗噪声稳定性。本教程的目标是提供使用这些框架的计算工具来探索神经影像数据，并促进该领域的可访问性、数据可视化和新手对该领域的理解。这里将首先介绍这两个框架，然后解释如何计算关于静息态功能磁共振成像(refMRI)中已经建立的和更新的指标。使用开源工具(Python)，并提供一个公开可用的Jupyter Notebook，它使用了1000功能连接组项目数据集，其中有一部分专门用于大脑网络中高阶交互的真实可视化。

前言

神经科学仍然是一个年轻的研究领域，大约在70年前作为一门正式学科出现。此后，这一领域迅速发展，目前关于人脑神经生物学的了解，很大程度上得益于高分辨率和不同规模的大脑研究技术的迅速发展。例如，磁共振成像(MRI)，它使我们能够非侵入性地测量大脑结构的区域特征，也可用于评估大脑区域之间的结构和功能的相互作用。该领域的这种扩展导致数据规模和复杂性呈指数级增长。为了分析和解释这种“大规模数据”，研究人员必须建立健全的理论框架。复杂的网络科学被引入神经科学，并越来越多地用于研究大脑复杂的通信和连接。由此产生的领域——网络神经科学——旨在通过对大脑的元素和相互作用进行映射和建模，并以此来观察大脑。

复杂网络科学中用于建模、估计和模拟大脑网络的主要理论框架之一是图论。一个图由一组相互连接的元素组成，也称为顶点和边。例如，网络中的顶点（也称为节点）可以是大脑区域，而边则表示（也称为连接）顶点对之间的功能连接。重建大脑网络的成像方式有很多种。本文的重点是静息态功能MRI(rsfMRI)。顾名思义，rsfMRI间接测量被试处于休息状态（即不执行任何任务）时的大脑活动。这种类型的数据提供了关于自发性大脑功能连接的信息。功能连接通常通过从结构分离的大脑区域测量信号之间的统计依赖性（通常是Pearson相关系数）来操作。

可以计算出几种描述图指标来形容大脑网络的特征；例如顶点的连接度或总数以及两个顶点之间的路径长度。这些指标使得研究人员能够识别大脑网络的非随机特征。一个突破性的发现是：大脑（与大多数其他现实世界网络一样）遵循一个“小世界网络”架构。这是指这样一种现象，即为了最小化连接成本，同时保持最佳效率和抗扰动的稳定性，大脑网络在执行局部处理（即分离）的能力和在全局水平上整合信息流（即整合）的能力之间保持平衡。因此，网络神经科学提供了一套全面的分析工具，不仅可以研究大脑区域的局部特性，还可以研究它们对整个大脑网络功能的意义。利用图论，研究者们已经收集了许多关于健康和疾病脑神经生物学方面的见解。

拓扑数据分析（TDA）可以通过分析一组顶点之间的相互作用（即高阶相互作用），而不是“简单”的成对连接，来提供关于大脑网络特征的另一视角。使用TDA，可以识别网络的“形状”及其不变属性。因此，TDA通常比图论分析提供更高的抗噪声稳定性，这可能是成像数据中非常重要的一点。尽管TDA最近才被用于网络神经科学，但它在rsfMRI研究上已经显示出令人兴奋的结果。但该框架的一个局限性是，应用TDA和解释结果所必需的数学抽象的复杂性可能会使未经数学训练的临床神经科学家无法使用它。此外，从TDA分析中产生的高阶交互结构通常难以真实地、易于理解地进行可视化。尽管存在技术上的限制，但TDA允许我们适当地处理高阶和庞大的组合编码容量。

通过提供有关如何计算研究大脑网络和现实高阶网络图的不同指标的分步教程，来促进网络神经科学及其组成部分图论和TDA的使用。研究者提供了这些指标的理论和实验背景，并在每个部分中提供相应的代码来解释如何计算不同的指标。此外，还列出了一些附加资源（表1和表2），包括为配合本实践教程而创建的Jupyter Notebook，都可在GitHub和Zenodo上公开获取。

表1.计算资源列表

表2.阅读资源列表

本文不同于以往的文献，因为本文描述了图论和TDA的核心概念，并提供了一个易于掌握的分步教程。此外，研究者提供了大脑中单纯复形的新3D可视化和TDA指标，这可能有助于这些工具的应用和解释。尽管本教程侧重于rsfMRI，但本文讨论的主要概念和工具可以拓展和延伸至其他成像模式、生物或复杂网络。本文提供的代码只需要具备功能连接矩阵的知识。对于单纯复形的3D可视化，只需要一个给定的大脑图谱节点的坐标。因此，该脚本可以适应不同的数据库，图像模式和大脑图谱。表3是一个简短的词汇表，其中包含了本文中的相关术语。

表3.关键术语词汇表

教程

Python 3工具包是执行以下计算步骤所必需的。附带的Jupyter Notebook可在GitHub（参见表1）或Zenod上找到。

代码示例

邻接矩阵

这里，在rsfMRI中应用图论和TDA的基本单位是邻接或功能连接矩阵（图1G），它表示网络中所有顶点之间的连接。通常，rsfMRI矩阵是对称的并且不指定连接方向（例如，区域A中的活动引发区域B中的活动），从而产生无向网络（图1B和F）。相反，非对称矩阵会产生有向网络。

图1.网络类型。

在计算这些矩阵的任何指标之前，以及在处理连接数据时必须考虑几个关键因素。一个关键决定是是否要保留有关边权重的信息。当边的权重（例如，rsfMRI连接中的相关值）保持不变时，将对网络进行加权（图1D和F）。另一种方法是使用任意阈值或密度，例如，仅保留和二值化20%的最强连接（图1A和B）。目前对于rsfMRI矩阵中的加权问题没有黄金标准，也可能取决于数据集或现有的分析方法。脑网络数据通常使用特定的阈值程序（或标准）进行分析。这些阈值大脑网络通常将大脑的普遍特征显示为一个复杂的系统，例如倾斜度分布、聚类、小世界和短平均路径长度等等。即使在一个阈值正态分布数据的情况下，也可以观察到其中一些被认为是复杂网络特征的属性。然而，一些大脑网络属性对阈值的变化并不稳健。这些结果提高了使用独立于阈值的网络分析方法的意识，例如最小生成树或拓扑数据分析。

关于rsfMRI矩阵的另一个相关讨论是负权重或反相关关系的解释。关于这种负相关在神经生理学中意味着什么的争论仍在继续。有研究表明，它们可以被认为是由全局信号回归或预处理方法引入的伪影，或者仅仅是由大脑区域之间同步信号的巨大相位差引起的。尽管如此，一些研究人员认为，反相关关系可能在长期同步下具有生物学意义，并且在疾病状态下，这些负相关的变化可能反映了网络重组。负权重可以绝对值化，以保留它们可能携带的潜在生物信息。如果决定丢弃它们，那么一些生理信息也可能会丢失。

本教程将使用一个无向、绝对（正）加权矩阵。在GitHub上的Jupyter notebook教程中，研究者提供了一个示例矩阵，它是存储库中所有可用连接矩阵的平均值。为了遵循下面的步骤，假设rsfMRI已经进行了预处理并转换为矩阵(因为预处理步骤不是本教程的目标)。

代码示例

在处理fMRI脑网络数据时，生成的一些图（例如，用于矩阵可视化的热图和边缘权重的分布图）有助于数据探索、理解和标记潜在的伪影。当绘制为log10的概率密度时，研究者期望权重分布具有类似高斯的形态。

代码示例

图论

这里将介绍网络神经科学中最常用的图指标（图2）。首先，需要使用NetworkX包创建一个图形对象并移除自环（即连接矩阵的对角线）。

图2.图论指标。

代码示例

度

顶点度量化了无向二元网络中的顶点连接总数。像rsfMRI矩阵这样的无向加权网络中，顶点度数类似于顶点强度（即，顶点的所有边的总和）并等效于其度的中心性。该指标是网络分析中最基本的指标之一，是对单个顶点连接密集程度的总结。

代码示例

通过从函数中删除参数权重，可以计算所有边为0或1的二值化网络的度（在使用稀疏/非全连接矩阵时很有用）。此更改将通过计算与顶点相邻的边数来得到顶点度，还可以删除指定的顶点来估计所有顶点的度/强度。度/强度分布允许我们通过显示网络是否包含一些高度连接的顶点(即枢纽)来确定一般网络组织的范围。

路径长度

①最短路径是网络中两个顶点之间的边数最少（或总权重最小）的路径。在加权图中，最短路径是通过两个顶点之间的边权重的最小总和计算得到的。它被视为了解网络中信息传播效率的一种测量方法。有几种算法可以计算路径长度，但Dijkstra算法是最经典和最著名的算法之一。但是该算法仅适用于具有非负权重的图。需要记住的是，在相关矩阵的情况下，例如rsfMRI数据，必须通过计算原始权重的倒数将权重转换为“距离”；相关值越大，表示距离越短。这种转换对于以下所有基于路径的指标都是必不可少的。

②平均路径长度（或特征路径长度）是网络中所有可能的顶点对的平均最短路径长度。它是衡量网络中信息传输效率和集成度的全局性指标，因著名的Watts-Strogatz模型而广为人知。

代码示例

聚类系数

聚类系数评估顶点的任意两个相邻之间的直接连接（或者在加权情况下更强的连接）趋势，也可以称为小圈子。该指标也用于计算小世界系数（特征路径长度与相对于随机网络的聚类系数之间的比率）。

代码示例

中心性

①特征向量（基于度）中心性衡量了顶点在网络中的重要性，而且还考虑其相邻节点的影响。因此，它同时考虑了顶点连接的数量和质量。对于加权网络，有一定的条件。根据Perron-Frobenius定理，邻接矩阵的最大特征值，必须是唯一且是正值。

②接近度（基于最短路径）中心性衡量一个顶点与网络其余部分的紧密度或“直接”连接程度。如果这个顶点距离网络中的所有其他元素最近，它就有可能快速有效地传播信息。形式上，顶点i的接近中心性是其到所有N-1个其他顶点的平均最短路径长度的倒数。在加权网络中，可以根据Dijkstra算法，通过考虑最短路径的总权重来估计接近中心性。

③中间（基于最短路径）中心性是指网络中通过特定顶点的所有顶点对最短路径的比例。它用于理解顶点对网络中整体信息流的影响。计算顶点i的中间中心性，必须计算两个顶点之间最短路径的比例。对于加权图，边必须大于零，并且考虑权重的总和。

代码示例

最小生成树

最小生成树是网络的主干，即确保路径存在于所有顶点之间而不形成循环所需的最小边集。生成树的算法主要有几种，其中Kruskal的算法在NetworkX中实现。简而言之，该算法对顶点之间的距离进行排序，首先添加距离最小的顶点，然后在添加的每条边上检查是否形成了循环。该算法不会保留导致形成循环的边缘。

代码示例

模块化

模块化说明了网络在不同模块（或社区）中的可分割性。模块的识别由社区检测算法执行。这里将使用Louvain算法。它以两步迭代的方式工作，首先通过优化局部模块来寻找社区，然后连接属于同一模块的顶点。

代码示例

拓扑数据分析

接下来，将在rsfMRI邻接矩阵上使用TDA。TDA可以通过解决图论中使用的成对连接之外的网络高阶结构来识别不同的网络特征。TDA通常使用拓扑和几何方法来研究数据的形状。与其他方法相比，TDA的一个核心特征是能够提供稳健的结果，即使是数据有噪声的情况下。此时，必须区分噪声和系统误差之间的差异，以便正确应用TDA。噪声是由完全随机影响感兴趣变量的测量因素引起的。然而，系统误差并不完全是由偶然因素决定的。它们是由系统地影响感兴趣变量测量的因素引起的（例如，由于观察或测量过程的不准确性）。在rsfMRI（和其他与大脑相关的测量）中，两种类型的噪声都可能存在。

在网络神经科学中使用TDA的好处之一是可以发现网络的全局属性，无论我们以何种方式表示网络，都可以找到这些属性。这些性质就是所谓的拓扑不变量。接下来，将介绍一些基本的TDA概念：过滤、单纯复形、欧拉特征、相变、Betti数、曲率和持久同源性，如图3所示。

图3.拓扑数据分析。

邻接矩阵和过滤

如前面关于图论的部分所述，对于基于rsfMRI的邻接矩阵执行阈值处理的必要性或标准没有达成共识。但是，TDA通过调查网络中所有可能阈值的功能连接性来克服这个问题。这种调查网络属性以寻找所有可能的阈值而不是选择固定阈值的过程称为过滤（图3B）。注意，过滤的概念不仅用于高阶交互，而且还应用于成对关系的图论分析。

代码示例

单纯复形

在TDA中，将网络视为一种称为单纯复形的多维结构。这样的网络不仅由一组顶点（0-单纯形）和边（1-单纯形）组成，而且还由三角形（2-单纯形）、四面体（3-单纯形）和更高的k维结构组成（图3A)。简而言之，k-单纯形是一个k-维的对象，是由网络的k + 1个顶点组成的子集（图4）。

图4.单纯复形。

我们可以通过多种方式将网络编码为单纯复形。然而，这里将专注于仅从大脑网络的‘团’来构建单纯复合形，即所谓的大脑网络的团复合形。在网络中，k-团是具有k个all-to-all连接节点的网络子集。0-团对应于空集，1-团对应于节点，2-团对应于连接，3-团对应于三角形，等等。在团复合形中，每个k + 1团都与一个k-单纯形相关联。本实践教程Jupyter Notebook中提供了可视化团复合形的代码。为了创建3D图，这里使用了Plotly (Inc. 2015) 中提供的网格算法，以及Fan等人(2016)和Bakker等人(2015)提供的整个大脑的网格表面。图5展示了1000功能连接组数据集中3-团的3D可视化示例图。当增加过滤密度d时，我们会获得更多的连接，并且会出现更多的3-团。在图5中，只显示出了3个团；但该方法同样适用于四面体等高维团。

图5.单纯形 3D 可视化。

欧拉特征

欧拉特征是拓扑不变量的一个例子：具有不依赖于特定图表示的网络属性。在3D凸多面体中，欧拉特征定义为所考虑的多面体的顶点数减去边数加上面数。对于与球体同构的没有空腔的凸多面体，欧拉特征始终为2。欧拉特征告诉我们关于多面体拓扑结构及其类似表面的一些信息。换句话说，如果我们有一个曲面并且我们对其进行离散表示，则无论采用何种方式，它的欧拉特性将始终相同。接下来，将欧拉特征的定义推广到任意维的单纯复形。

代码示例

拓扑相变

在比较患者和健康人群时使用TDA，可以通过研究每个组的拓扑相变曲线来识别这些人群的特性。该策略已应用于调查对照组和神经胶质瘤脑网络之间的群体差异，以及典型发育中的儿童和患有注意力缺陷/多动障碍儿童之间的差异。为了研究大脑网络中的拓扑相变，首先需要可视化欧拉熵。在特定假设下，当欧拉特征为零时，复杂网络中会发生拓扑相变。在网络理论中，大规模分量跃迁与网络变化有关，从较小的连接团簇到大规模分量的出现。从理论上讲，拓扑相变与单纯复形的大规模分量跃迁的变化有关。

Betti数

另一组拓扑不变量是Betti数(β）。假设单纯复形是高维结构，βk表示单纯复形中k维孔的数量。对于每个k≥0，这些拓扑不变量对应于单纯复形中线性独立的k维孔的数量。图6展示了k维孔的示意图。在单纯复形中，可能有许多这样的k维孔，对它们进行计数得到Betti数β。

图6.Betti数和每个k维孔的示意图。

代码示例

图7说明了从随机网络获得单纯复形的这一性质。随着概率的增加，网络的密度越高，研究者发现了一个显性序列βk，从k=0开始，每当发生拓扑相变时 ，该序列（即k增加一个单位）就会发生变化。

图7.Betti数和欧拉特征近似。

曲率

曲率是一种TDA指标，可以将上述全局网络属性与局部特征联系起来。在处理脑网络数据时，这将使我们能够计算全脑顶点集的拓扑不变量，并了解特定的单个节点或子网络几何特性对脑网络全局属性的贡献。为了将节点曲率与网络的全局属性联系起来，将使用网络的高斯-博内定理，通过该定理可以连接网络的局部曲率及其欧拉特征。

图8.曲率3D图。

持久同调

同调性是一个拓扑分支，它通过研究对象的孔/环来研究其形状。持久同调性使我们能够识别是否有同调类在许多过滤中“持续”存在。重要的是，要计算持久同调性，需要使用距离矩阵，这是下面代码中的第一步。然后，可以计算单纯复形的持久性并将其绘制为条码或持久性图（图3C和E）。这里，将使用到Gudhi包来实现这些步骤。

代码示例

讨论

本教程解释了与网络神经科学两个分支（图论和TDA）相关的主要指标，提供简短的理论背景和代码示例，并附有公开可用的Jupyter Notebook。将这些子领域的现成代码和大脑中单纯复形的可视化相结合，并且通过分步解释，希望能够降低神经科学家熟悉这些新分析方法的门槛，特别是应用于rsfMRI数据中的新方法。本文的主要目标是提供一个循序渐进的计算教程，以在脑成像数据（尤其是rsfMRI）上使用图论和TDA，并对每个指标进行深入解释。将这些分析框架应用于大脑数据的核心思想是，这两种框架都可以定量地结合大脑的两个明显的基本特征：大脑不仅在特定大脑区域的局部水平上工作，而且还包含对其功能至关重要的全局属性。

网络神经科学是理解大脑组织和功能的关键。图论是迄今为止使用最多的框架，但随着网络神经科学领域的发展，TDA等更新的方法开始运用到研究中。为了进一步改进该领域，特别是临床网络神经科学领域，当务之急是要使已开发指标的计算易于访问、易于理解、可视化和更加高效。此外，研究人员必须意识到在执行数据分析时必须做出的关键决定，以及这些决定如何影响研究结果及其再现性。希望本教程能够促进对网络和拓扑神经科学的某些方面的理解，以及某些指标的计算和可视化。

原文：A hands-on tutorial on network and topological neuroscience.

https://doi.org/10.1007/s00429-021-02435-0