人脑功能的几何约束

大脑的解剖结构必然会限制它的功能，但确切的原因仍然不清楚。神经科学的经典和主导范式是，神经元动力学是由复杂的轴突纤维阵列连接的离散的、功能特化的细胞群之间的相互作用驱动的。然而，来自神经场理论（一个用于模拟大规模大脑活动的确立已久的数学框架）的预测表明，与复杂的区域间连接相比，大脑的几何形状可能代表了对动力学更基本的约束。在本研究中，我们通过分析自发和不同任务诱发条件下获得的人体磁共振成像数据，证实了这些理论预测。具体而言，我们表明，皮质和皮质下的活动可以简单地理解为由大脑几何（即其形状）的基本共振模式的兴奋引起，而不是像传统假设的那样，由复杂的区域间连接模式引起。然后，我们使用这些几何模式表明，超过10,000个脑图谱上的任务诱发激活并不像人们普遍认为的那样局限于焦点区域，而是激发了波长超过60㎜的全脑模式。最后，我们证实了预测结果，即几何和功能之间的密切联系可以由波样活动的主导作用解释，表明波动力学可以再现自发性和诱发记录的许多典型时空特性。我们的研究结果挑战了主流观点，并确认了之前未被充分认识的几何形状在功能塑造中的作用，正如一个统一的、有物理原则的全脑动力学模型所预测的那样。

概述

许多自然系统的动力学从根本上受到其底层结构的限制。例如，鼓的形状影响它的声学特性，河床的形态决定水下水流的形状，蛋白质的几何形状决定与之相互作用的分子。神经系统也不例外，其复杂的轴突连接网络支持了解剖分布的神经元群体丰富而复杂的时空动力学。一些研究已经显示了脑连接和活动的各种特性之间的相关性，但神经动力学的时空模式如何受到相对稳定的神经解剖支架的精确限制仍不清楚。
在物理和工程的不同领域中，对系统动力学的结构约束可以通过系统的本征模来理解，本征模是与系统的自然共振模态相对应的基本空间模式。在线性状态下，例如正常（即非癫痫样）条件下的大脑活动，本征模（下文也称为模式）提供了一种特别强大和严格的形式体系，将大脑解剖与塑造活动的物理过程联系起来。通过这个透镜，神经元动力学的时空模式从大脑结构本征模的激发中显现出来，就像拨动的小提琴弦的谐波由它自身的共振模式的振动产生一样。
关键的是，就像小提琴弦的共振频率由它的长度、密度和张力决定一样，大脑的本征模也由它的结构物理、几何和解剖学特性决定。这些特定的结构特性对动力学有重要的贡献吗?在这里，我们测试了两种有影响力且相互对抗的理论，它们对大脑结构的哪些关键元素塑造功能做出了不同的预测。
一个经典的观点代表了神经科学的主导范式，它起源于 Ramon y Cajal 的神经元学说、Brodmann的细胞架构学说以及一个多世纪以来对特定大脑区域功能定位的研究。根据这一观点，神经动力学的时空模式源于离散的、功能特化的细胞群之间的相互作用，这些细胞群由拓扑结构复杂的短程和长程轴突连接阵列连接。在人类中，这些连接可以通过弥散磁共振成像(dMRI)在宏观尺度上进行估计，从而产生一个结构连接矩阵，称为连接组。这种方法已被广泛用于理解大脑的组织和动力学，最近的工作提出，从这样一个离散的连接组衍生出的本征模（这里称为连接组本征模）可以用来重建由功能磁共振成像(fMRI)映射的人类皮层的典型功能网络的空间模式。
这种基于离散连接组的观点的一个局限性是，它依赖于大脑解剖的抽象表示，而不能直接解释其物理属性和空间嵌入（即几何和拓扑）。这些特征被明确地纳入了一大类神经场理论(NFT)中，这些理论描述了超过0.5㎜空间尺度上的平均场神经动力学。一种受生理限制的NFT形式将皮质活动视为通过物理连续的神经组织片层传播的行波的叠加，从而统一了一系列不同的经验现象。在这一理论中，不同皮质位置之间的神经相互作用由一个均匀的空间核来近似，该核大致随距离呈指数下降，这与实验证据一致，即许多物种的神经系统的组织普遍受指数距离规则(EDR)的连接控制。
考虑到波样动力学和EDR样连接，NFT的一个关键预测是，大脑的固有几何形状，不仅在物理上塑造了其所涌现的动态，还对它施加了边界条件。这一观点的一个非平凡的推论是，如果我们优先考虑大脑解剖中的空间和物理限制，那么我们只需要考虑大脑的形状，而不是其拓扑上复杂的轴突相互连接的全部阵列，就可以理解它的空间活动模式。更正式地说，该理论预测，来自脑几何的本征模（以下简称几何本征模）比连接组代表了对动力学更基本的解剖学约束。这一观点与传统观点形成鲜明对比，传统观点认为复杂的区域间解剖连接模式塑造了大脑活动。
在这里，我们测试了这些关于大脑的相互对抗的观点，目的是确定人类大脑动力学的主要结构约束。与NFT的理论预测相一致，我们发现，与从脑连接测量值获得的本征模（连接组本征模）相比，从人类新皮质的自发和任务诱发记录获得的各种fMRI实验数据可以更简洁地用皮质几何特征（几何本征模）来解释。我们进一步证实，刺激诱发的活动是由具有长空间波长的几何本征模的激发所主导，这挑战了将这种活动局限于聚焦的、空间孤立的簇的经典观点。为了直接将这些结构约束与驱动脑动力学的物理过程联系起来，我们使用生成模型来展示在皮质几何上展开的波动力学如何可以解释功能性脑组织的不同特征。最后，我们展示了由本征模捕获的几何和功能之间的密切关系可以推广到非新皮质结构，表明这种联系是大脑组织的普遍属性。

几何模式限制皮质活动

我们首先研究几何本征模在多大程度上可以解释人类新皮质活动的不同方面。为了推导本征模，我们使用了新皮质表面总体平均模板的网格表示（图1a和方法中皮质几何本征模的推导）。然后，我们从该曲面网格构造Laplace–Beltrami算子(LBO)，该算子捕获局部顶点到顶点的空间关系和曲率，并求解特征值问题，

$${\nabla }^2\psi =\Delta \psi =-\lambda \psi \text{\hspace{1em}(1)}$$
其中∇为梯度算子，Δ为LBO， ψ = { ψ 1 ( r ) , ψ 2 ( r ) , … } ψ=\{ψ_1(\mathbf{r}),ψ_2(\mathbf{r}),…\} ψ={ψ1​(r),ψ2​(r),…}为几何本征模族及其对应的本征值族 λ = { λ 1 , λ 2 , … } λ=\{λ_1,λ_2,…\} λ={λ1​,λ2​,…}。本征值按照各模态空间模式的空间频率或波长（图1a）顺序排列，其中${\psi }_1$是波长最长的模态。得到的本征模是正交的，形成了一个完整的基集，将在皮层上展开的时空动力学分解为不同波长模态的加权总和（图1b和方法中的脑活动的模态分解）。除非另有说明，我们在整个研究中使用N=200种模态。

图1：用几何本征模重建新皮质活动

a，通过求解特征值问题，从皮质表面网格导出几何本征模，Δψ=−λψ（方程(1)）。模态${\psi }_1,{\psi }_2,{\psi }_3,\dots ,{\psi }_N$按从低到高的空间频率（长到短的空间波长）排序。负值、零值和正值分别被标记为蓝色、白色和红色。b，脑活动数据的模态分解。这个例子展示了一个空间映射y(r,t)，在给定时间t，可以分解为$a_j$加权的模态${\psi }_j$的总和。c，左，我们使用一系列刺激对比的激活空间图来重建任务诱发数据。右，我们通过分解每个时间帧的空间图来重建自发活动，并生成一个区域到区域的FC矩阵。d, 7个关键HCP任务对比图的重建精度以及静息态FC与模式数量的关系。图中显示了皮质表面的重建，显示了与前10、100和200模态（分别对应于大约120、40和30㎜的空间波长）相关的空间尺度。e，使用7个关键HCP任务对比的10、100和200种模态获得的实验组平均任务激活图和重建(recon.)。黑色箭头表示局部激活模式，当使用短波长模态时，可以更准确地重建。f，实验组平均静息态FC矩阵和使用10、100和200模态的重建。

利用这种分解，我们评估了几何本征模在捕获任务诱发和自发脑活动（图1c）方面的准确性，这些活动是从人类连接组项目（Human Connectome Project, HCP; 方法中的HCP数据）中的255个健康个体中测量出来的。对于任务诱发的活动，我们从代表不同诱发激活模式的7个不同任务绘制了47个基于任务的对比图。然后，我们使用递增的模态（最多200种）来重建每个个体的激活图（图1d）。对于自发的、无任务的（所谓的静息态）活动，我们重建每个时间帧的活动空间图，然后生成一个区域到区域的功能耦合(Functional Coupling, FC)矩阵，描述每个半球180个离散脑区之间活动的相关性。为了直接比较任务诱发和自发记录，我们对任务诱发数据应用了相同的脑区划分（方法中的皮质区划）。最后，我们通过计算任务诱发的激活图与自发FC矩阵实验和重建之间的相关性来量化重建精度（图1d-f）。
我们观察到，在所有任务对比和静息状态下，重建精度随着模态数量的增加而增加，仅使用N=10个模态就已经实现了r≥0.38（图1d）。不同的任务也会调动不同的大尺度模态，这表明特定的刺激会激发特定的模态（图1e）。超过10种模态后重建精度的提高变得缓慢，大约在N=100种模式时达到r≥0.80，超过此点重建精度仅呈渐近递增趋势。由于前100个模态的波长都超过40㎜，并且包含较短波长的模态只会细化局部模式的重建（图1e中的箭头），我们的发现表明，数据主要由具有较长空间波长的空间模式组成（见下一节的更详细分析）。
这些结果在所有47个HCP任务对比和不同分辨率的区划中都是一致的，但由于较粗略的区划的低通空间滤波效应，较高分辨率的区划数据需要更多的模态才能达到较高的重建精度。我们的结果也不受使用群体平均皮质表面模板（而不是个体特定表面）来推导几何本征模的影响。综上所述，这些发现表明，皮质几何本征模形成了一个紧凑的表征，捕捉了任务诱发和自发皮质活动的不同方面。此外，它们还表明，这种活动是由长波长、大尺度的本征模主导的。
接下来，我们检验一个假设，即几何本征模比基于图的连接组近似得到的本征模提供了更简洁和更基础的动力学描述。为此，我们比较了几何本征模与三种备选连接组衍生的本征模基集的重建精度（示意图见图2a）。第一个基集是从dMRI示踪图映射出的实验连接组，具有顶点级分辨率并经过阈值化，获得0.10%的连接密度，正如之前所做的（在方法中，连接组本征模的推导）。第二个基集是一个由指数距离依赖的连接概率控制的均匀随机连接过程合成的连接组来模拟简单的EDR样连接（方法中，EDR本征模的推导）。由于实验连接组和EDR连接组的连接密度不同，我们还评估了使用1.55%阈值的实验连接组的第三个基集，以匹配EDR连接组的密度。上面描述的连接组、EDR和密度匹配的连接组本征模是由它们各自连接矩阵的图拉普拉斯（LBO的离散形式）导出的（图2b）。

图2：以基于连接组的本征模为基准的几何本征模

a，用于推导皮质几何本征模、连接组和EDR连接组的解剖特性示意图。几何本征模依赖于局部表面网格信息，如相邻表面网格顶点（点）之间的链接（蓝色）和曲率。连接组本征模依赖于从dMRI实验重建的网格顶点（蓝色）之间的局部连接和短、长连接（品红）。EDR本征模依赖于由随机布线过程产生的连接（红色），在这个过程中，顶点之间连接的概率随它们的距离呈指数衰减。b，示例连接组和EDR本征模。负值、零值和正值分别被标记为蓝色、白色和红色。尽管有一些相似之处，但这些模态的空间模式与使用皮质几何导出的模式不同（与图1a相比）。c，静息态FC矩阵重构精度，通过几何、EDR和两种不同的连接组本征模实现：一种是使用之前方法定义的连接组，另一种使用与EDR连接组相同的连接密度来进行公平比较。d，通过几何本征模和其他基集实现的所有47张HCP任务对比图的重建精度差异，如每个面板上方的文本所示。每一行代表不同的任务对比，在这里按大类分组；红色表示几何本征模的优越性能。需要注意的是，在少于10个本征模的重建中，连接组比几何表现更好；但相比100个本征模的重建（不同任务的平均r=0.71），精度通常较低（平均r=0.42）。

综上所述，几何本征模解释了皮质表面的内在曲率和表面网格中的局部顶点到顶点关系；连接组本征模不考虑曲率，而是捕捉网格顶点之间的局部空间关系，以及dMRI测量的短连接和长连接；EDR本征模解释了均匀的、随机的、依赖距离的连接规则的影响，而没有完全捕捉皮质的几何信息（图2a）。通过比较这些不同的基集，我们可以从结构连接中分离出皮质的几何结构对脑动力学的贡献。
对这些不同基集的重建精度的直接比较表明，在自发（图2c）和任务诱发（图2d）数据中，几何本征模一致地显示出最高的重建精度。EDR本征模的表现几乎与几何本征模一样好，而连接组本征模的精度最低。无论使用的区划、用于生成连接组本征模的特定连接密度，以及我们是否使用离散的区划而不是以顶点分辨率生成连接组，这一发现都是正确的。此外，我们还发现几何本征模比功能数据的主成分本身（通过主成分分析(PCA)计算）具有更强的样本外泛化性，并且比Fourier空间基集（方法中的与统计基集的比较）具有更好的表现。
综上所述，这些结果证明了几何本征模作为脑功能基集的简洁性、稳健性和普适性。它们也支持NFT的预测，即大脑活动最好以直接来自于皮质形状的本征模来表示，因此强调了几何在约束动力学方面的基本作用。

长波长支配皮质活动

利用几何本征模对自发数据和任务诱发数据进行的重建表明，大脑活动的空间组织由空间波长约40㎜或更长的模态主导（图1d-f）。这一结果与神经影像学数据分析的经典方法相反，在经典方法中，刺激诱发的激活是通过阈值化统计图来确定活动增强的局部性孤立区域。这种经典的方法建立在这样的假设上，即局部焦点代表被刺激激活的离散脑区，而其他区域的阈下活动可以忽略不计。任务激活数据惊人的长波长内容（图1d-e）表明，经典流程只关注了冰山一角，而掩盖了由任务诱发的潜在的空间扩展的和结构化的活动模式。这些观察结果与NFT的理论预测和之前对任务诱发脑电图(EEG)信号的分析一致。
在这里，我们利用如图1b所示的模态分解来表征任务诱发激活的完整空间模式——整个冰山。为此，我们分析了HCP（方法中任务诱发激活图的模态功率谱）中使用组平均未阈值化激活图的几何模态分解得到的空间功率谱。作为一项独立的重复试验，我们还分析了来自NeuroVault存储库中1,178个独立实验的10,000个未阈值化激活图，从而全面了解了在人脑中绘制的刺激诱发激活模式的多样性。
尽管我们使用了广泛的刺激、范式和数据处理方法来获取这些激活图，但我们观察到，激活图中的很大一部分功率集中在前50种模态中，对应于大于60㎜的空间波长（图3a；对于每个关键的HCP任务对比图，也分别发现了类似的结果）。使用合成的数据，我们证实这些发现不能用典型的fMRI处理流程产生的空间平滑来解释，这种处理流程可以过滤掉短波活动的空间模式。我们进一步观察到，与去除短波模态相比，从小到大依次去除长波模态对重建精度的影响要大得多（图3b和方法中长波和短波模态的贡献）。例如，在7个关键的HCP任务对比中，去除前25%的长波模态（模态1-50）会使重建精度下降约40-60%，而去除前25%的短波模态（模式151-200）只会使重建精度下降约2-4%（图3b，插图）。这些结果表明，在fMRI可达到的时间和空间尺度上，诱发的皮质活动包括大规模的、几乎是全脑的空间模式，这挑战了传统观点，即此类活动应以离散、孤立和解剖学上局部的激活簇来描述。

图3：任务诱发活动激发长波模态

a，47个HCP任务对比图（左）和来自NeuroVault数据库的10,000个对比图（右）的归一化平均功率谱。插图的皮质表面重建显示了与空间波长分别约为60、40和30㎜的前50、100和200模态相关的空间尺度。b, 7个关键HCP任务对比图的重建精度与重建过程中去除的模态（从200个模态中）的百分比有关。实线和虚线分别对应长波和短波模态的依次去除。插图显示组平均实验激活图（数据）及其在去除25%的模态后的重建。负值、零值和正值分别被标记为蓝色、白色和红色。

波动力学桥接几何和功能

通过求解LBO的特征值问题，即Helmholtz方程（式(1)），得到皮质的几何本征模。在物理连续系统中，Helmholtz方程的解对应于更一般的波动方程解的空间投影，这样得到的本征模天然地表征了系统动力学的振动模式或驻波。这一等效性表明，正如NFT所预测的那样，几何本征模在重建不同的脑活动模式方面的优越效能，源于波动力学在这些模式的塑造中发挥的基础性作用。这一预测已经通过脑电图记录的模型得到了证实，但直到最近才在fMRI信号中观察到整个大脑的波，因此目前还缺乏理论解释。在这里，我们使用NFT和几何本征模来表明，波动力学可以在fMRI可及的尺度上为不同的实验和生理现象提供统一的解释。
我们使用无再生的各向同性阻尼NFT波动方程来模拟神经活动（图4a和方法中的NFT波模型）。在这个模型下，活动通过白质连接在新皮质的点之间传播，其强度随着距离的增加而近似指数衰减。为了模拟静息态神经活动，我们使用白噪声输入来模拟非结构化的随机波动（方法中的静息态动力学建模）。我们比较了这种简单的波模型和基于生物物理学的神经质量模型（兴奋-抑制平衡，balanced excitation–inhibition (BEI) 模型）的性能，后者已被广泛用于理解静息态fMRI信号（图4a和方法中的神经质量模型）。神经质量模型与脑功能的经典的、以连接组为中心的观点紧密一致，代表了由离散解剖区域的神经元群相互作用（根据实验测量的连接组耦合）导致的动力学。

图4：波动力学塑造自发和刺激诱发的活动模式

a，对于波模型，位置r和时间t的活动φ(r,t)由具有阻尼率${\gamma }_s$、空间长度尺度$r_s$和输入Q的波动方程控制。对于神经质量模型，区域i的活动$S_i(t)$由函数f描述，该函数依赖于其他区域的活动S、局部群体参数${\theta }_i$和由全局耦合参数G缩放的连接组C。模型动态被用来计算一个模拟的FC矩阵（方法）。$M_{fixed}$和$M_{free}$分别对应每个模型的固定参数和自由参数的数量。b，数据和模型仿真的比较，基于各种度量，从左到右依次为：FC矩阵（用于视觉）、边FC、节点FC和FCD。对于边FC和节点FC，红线表示线性拟合，Pearson相关系数r；对于FCD，利用 Kolmogorov-Smirnov (KS) 统计量比较数据和模型动力学的区域间同步相似性的概率密度函数(pdf)。c，刺激 V1 1㎳ 后，t=1-2㎳活动的波传播。箭头表示传播方向。d，视觉皮层体系中17个区域的活动特征。插图是根据激活谱着色的皮层表面区域的空间位置。e，d中各区域活动谱达峰时间排序号和T1w:T2w值排序号的关系。红线代表排序变量的线性拟合，Spearman相关系数r和单侧$P_{spin}$，从10,000个排列中估计。

首先，我们比较了两种模型捕捉通常研究的不同FC（包括多种自发的、无任务的FC：静态两两FC（边FC），静态节点级平均FC（节点FC）和时间解析的FC动态属性（FCD）（在方法中，静息态动力学建模））特性的效能。在所有基于FC的基准测量中，波模型在实验数据重建方面显示出与神经质量模型相当或更优的性能（图4b）。波模型也比质量模型更准确地捕捉了实验中静息态活动的时滞特性（在方法中，静息态动力学时滞特性的测量）。波模型的简洁性赋予其显著强大的性能：波模型只需要皮质的几何形状（即表面网格）作为输入，并包括一个固定参数和一个自由参数($r_s$)，来拟合数据；而神经质量模型需要一个来自dMRI的区域间解剖连接矩阵，包括15个固定参数和4个自由参数。这表明，波动力学为fMRI捕捉到的宏观自发性皮质动力学提供了更准确和简洁的机制描述。
接下来，我们考虑波模型中刺激诱发的皮质活动。我们分析了初级视皮层(V1)对感觉刺激的皮质反应，因为它会引发一个已明确定义的区域皮质层级反应（在方法中，刺激诱发的动力学建模）。一个1㎳的脉冲输入到V1会产生一个沿背侧和腹侧视觉处理流快速分裂的活动传播波（图4c（箭头）），这与分层视觉处理的主流理解一致。值得注意的是，这一结果表明，对诱发活动的行波的几何约束足以分离背侧和腹侧处理流，而传统上认为这主要是由层特异性连接的复杂模式驱动。此外，在整个视觉系统中，诱发反应的时间分布遵循明确的时间尺度层次，与低阶视觉区域相比，高阶关联区域显示出延迟和延长的峰值反应（图4d）。因此，这些发现表明，这种先前在实验和模型研究中确定的层次顺序，自然地出现在通过皮质介质传播的激发波中。重要的是，区域反应的这种分层时间顺序于另一项独立的解剖学测量（基于髓质结构的无创性评估（T1加权(T1w)和T2加权(T2w)比率）的皮质处理层次）强相关。这种相关性，在视觉处理层次中特别强（r=−0.72，单侧自旋测试P值($P_{spin}$)=0.003；图4e），但考虑到所有皮质区域时也存在（$r=-0.44,P_{spin}=0.037$）。总之，我们的建模结果展现了，在大脑皮层的几何结构上展开的简单的波动力学，如何，为捕捉大脑时空活动的复杂属性，提供了统一的生成机制。

几何形状限制皮质下的活动

到目前为止，我们的分析集中在新皮质的几何和动力学的强耦合上。接下来，我们在非新皮质区域研究这种耦合，重点是丘脑、纹状体和海马，因为这些结构的几何形状很容易通过MRI数据获得，并且它们的功能组织已经被广泛研究。
我们首先将本征模分析推广到三维(3D)体积（在方法中，非新皮质结构的几何本征模的评估），产生通过每个结构的三维空间维度在空间上扩展的几何本征模。接下来，为了全面捕捉这些非新皮质区域的宏观功能组织，我们将广泛使用的流形学习程序应用于体素级FC数据，以获得每个结构中的关键功能梯度（在方法中，描绘非新皮质结构的功能组织）。这些功能梯度描述了FC相似性所支配的空间组织的主轴，从而代表了功能组织变异的主要模态，并根据它们所解释的FC相似性的方差百分比进行排序。
丘脑、纹状体和海马的前三个功能梯度的空间分布（分别占FC相似性方差的24%、50%和47%）显示出与前三个几何本征模的近乎完美匹配（图5a-c；空间相关性r≥0.93）。这种紧密的对应关系可以推广到每个结构的前20个模态和前20个梯度（前20个梯度分别占FC相似性总方差的49%、70%和68%），除纹状体和海马的第20个梯度和第20个模态外，所有的绝对空间相关性|r|>0.5（图5d-f）。这种强烈的关系是惊人的，因为功能梯度是通过一个复杂的处理管道（应用于fMRI衍生的FC测量）产生的；而本征模是简单地从每个结构的几何形状得到的，独立于功能数据。这些发现表明，非新皮质结构的功能组织直接来源于它们的几何本征模。

图5：几何塑造非新皮质的功能

a-c，丘脑(a)，纹状体(b)和海马(c)的前3个几何本征模和基于FC的功能梯度。模态和梯度在三维坐标空间中分别用蓝色，白色和红色表示为负值，0和正值。D，背侧；P，后端；R，右侧。散点图显示了模态和梯度之间的关系，红线表示Pearson相关系数r的线性拟合。d-f，下半，丘脑(d)、纹状体(e)和海马(f)的前20个几何本征模和功能梯度的绝对相关性(|r|)。上半，每个功能梯度的最大|r|值（灰色条，考虑到几何本征模的顺序翻转），以及每个功能梯度解释的方差百分比（蓝线）。Max表示最大。

讨论

许多物理系统的动力学受其几何形状的限制，可以理解为对相对较少数量的结构模态的激发。在这里，我们表明，与基于连接组的其他模型相比，仅来源于大脑几何的结构本征模能更紧凑、准确和简约地表示其宏观活动。这种基于模态的大脑视图进一步表明，fMRI捕捉到的自发和诱发的脑活动都由波长相对较长的大尺度本征模主导，其动力学来自生物物理学激发的波动方程。这些发现挑战了经典的神经科学范式，在这一范式中，离散的、特殊的神经元群体之间的区域间连接的拓扑复杂模式被视为动力学的关键解剖学基础。相反，我们的结果表明，这种将大脑视为一个连续的、空间嵌入式系统的，基于物理的方法，为理解宏观神经元功能的不同方面的结构约束提供了一个统一的框架。
对几何本征模与其他解剖学（连接组和EDR本征模）和统计学（PCA和Fourier）基集进行的广泛比较表明，几何本征模在捕捉宏观新皮质活动方面的优越性能并非由基集展开的一般数学特性所驱动。相反，这一结果表明，几何代表了动力学的基础解剖约束。此外，从一个人工合成网络导出的EDR本征模的强大性能表明，同质的、距离依赖的近指数形式的连接，代表了对活动的另一个重要的解剖学限制。EDR样连通性在数学上内含于方程(1)中的Helmholtz方程，因此这种连通性的作用被几何本征模隐式地捕获。
连接组本征模相对较差的性能表明，为了获得能够准确解释fMRI测量的皮层活动的时空模式的本征模，考虑那些存在于简单EDR之外的拓扑复杂连接的额外收益很小。因此，我们的发现与传统观点相反，传统观点强调复杂的解剖连接模式是协调动力学的主要驱动因素。事实上，最近的研究表明，长程皮质连接是罕见的——因此，它们可能代表了EDR样连接主导的效应的一个相对较小的扰动。尽管如此，这种连接的拓扑中心性、代谢成本和紧密的遗传控制表明，它们提供了比波样动力学更重要的功能和进化优势。dMRI和fMRI数据预处理流程的有限分辨率和敏感性，使得充分揭示这些连接的功能变得复杂。高质量的动物追踪和电生理数据可能有助于这方面的研究。
在新皮质和非新皮质结构中，几何和动力学之间的紧密耦合都是明显的，这表明新皮质以外区域的功能组织也由距离依赖的解剖连接和波动力学主导，正如在最近的实验中发现的那样。这些观察表明，几何本征模提供的描述（关于非新皮质结构中公认的功能组织梯度）更简洁、包含更多机制信息，相比于目前文献中使用的复杂的流形学习程序。这是因为这些程序是现象学的，提供的是对数据中方差的主要来源的统计描述；而对结构本征模的研究则来自一个生成过程。
几何模态分解为脑激活图的空间特性提供了独特的见解。经典的脑图分析侧重于超过统计阈值的孤立空间位置集群的反应。相比之下，我们的方法与严格建立的物理学和工程学结果一致，在这些结果中，空间连续系统的扰动会引起全系统的响应，就像小提琴弦的音调是由其整个长度的振动而不是由一个受限部分的行为产生的一样。因此，几何本征模的使用表明，在来自基于任务的fMRI研究的10,000多张不同图谱中，任务主要与波长约60㎜及更长的模态的激发相关。这一结果与在实验脑电图和诱发反应电位数据中对长波兴奋的类似观察结果一致，并表明依赖于阈值化的点统计图的经典分析掩盖了任务实际诱发的空间扩展和复杂的活动模式。
我们的建模结果提供了对物理过程的洞察，这是观测到的几何和功能之间的密切联系的基础。特别是，波模型在捕捉自发fMRI动力学的各个方面的相对简单性和优越表现表明，该模型比复杂的神经质量模型提供了更精简的解释，后者将大脑视为通过连接组（边）耦合的离散解剖区域（节点）的图。这一发现与在人类和动物fMRI数据中对波动力学的实验观察一致。未来的工作可以探索，在波模型中引入空间异质性或复杂的结构输入，是否可以进一步提高其解释各种实验现象的准确性。
应用波模型模拟视觉刺激表明，从刺激部位传播的波沿着经典的背侧和腹侧视觉通路分离，并且对扰动的局部反应符合一个已被充分描述的时间尺度层次，从快速反应的单模态区域到较慢的跨模态区域。分层视觉处理的这些典型特性已经被广泛研究了几十年，通常认为是由特定层间区域连接的复杂模式驱动的，但我们的分析表明，通过皮质几何的波足以产生分离的、分层的处理流。因此，虽然我们的研究结果不能排除复杂的区域间连接的作用，但它们确实表明，这种连接对这些宏观动态的出现不是必需的。
几何本征模的优越性能提供了直接的实际好处，因为仅使用感兴趣结构的网格表示就可以估计模态，而利用T1w解剖图像的成熟的自动化处理管道可以很容易地推导出模型。相比之下，连接组本征模需要：一个基于图的宏观区域间连接模型，该模型通过应用于T1w和dMRI图像的复杂数据处理管道生成；图节点的定义，这是一个充满争议的问题；以及应用阈值程序去除推定的虚假连接，我们自己的分析表明，这可能会影响研究结果。事实上，这些选择并不需要获得几何本征模，这意味着它们可以在人类和其他物种的不同实验环境中稳健和灵活地应用，开辟了新的研究途径。例如，可以研究几何本征模如何随着神经发育而变化或在临床疾病中被破坏。事实上，我们所确定的几何和功能之间的密切联系暗示了物种间在时空动力学上的差异可能很大程度上是由大脑形状的差异所驱动的。研究物种内部和物种之间的脑几何变化如何塑造大脑功能，对于理解神经元活动的物理和解剖学约束至关重要。

方法

皮质几何本征模的推导

如果大脑结构在时间上可以近似为常数，那么由此产生的空间和时间动力学可以通过本征模分解分别处理，类似于处理其他物理系统。特别地，其空间方面满足Laplacian特征值问题，也称为Helmholtz方程，定义在式(1)中。
对于大脑皮质，我们认为它是嵌入在3D欧氏空间中的二维(2D)模型，方程(1)中的LBO捕捉了固有几何形状，其中包括皮质表面的曲率，一般定义为：

$$\Delta :=\frac{1}{W}\sum _{i,j}\frac{\partial }{\partial x_i}(g_{ij}W\frac{\partial }{\partial x_j})\text{\hspace{1em}(2)}$$
其中$x_i,x_j$是局部坐标，$g_{ij}$是内积度量张量的倒数$g_{ij}:=⟨\frac{\partial }{\partial x_i},\frac{\partial }{\partial x_j}⟩,W:=\sqrt{det(G)}$，det表示行列式，$G:=(g_{ij})$。
我们使用安装在MASSIVE高性能计算设备中的 LaPy python 库来推导人类皮层的几何本征模。具体而言，我们使用了人类皮质厚度中层表面的三角形表面网格表示，每个半球包含32,492个顶点，这些顶点来自FreeSurfer的平均人群模板的下采样、左右对称版本。该模板与我们所有分析中使用的数据样本无关，因此避免了任何关于循环论证的问题。
注意，连续的LBO作用于曲面的Riemannian流形，而不是直接作用于网格顶点。LaPy采用曲面网格的三次有限元方法，在插值平滑流形上实现方程(1)的数值解。这区别于离散的图Laplacian，后者不编码点之间的空间关系。我们所有的分析都集中在单侧半球的本征模，但我们的方法可以很容易地推广到整个大脑，因为双侧半球的本征模可以表示为来自每个半球的本征模的对称或反对称组合：对称组合对应于横跨矢状面中部的镜像对称，非对称组合对应于半球具有相同空间结构但正负号相反的情况。
方程(1)的本征值解按各本征模空间图案的空间频率或波长顺序排列，即$0\le {\lambda }_1\le {\lambda }_2\le \dots$。注意，第一个特征值${\lambda }_1$近似等于零（波长远大于大脑大小），而对应的特征模${\psi }_1$是一个没有波节线的常数函数（函数的零集）。在我们的研究中，我们使用了前200个模态（包括常数模态，${\psi }_1$）进行分析，因为当使用越来越多的模态时，重建精度的提高越来越小（图1d）。
每个本征模由具有特定空间波长的空间图案组成。我们使用理想化的球形情况来近似本征模波长，因为它在拓扑学上与人类皮质相似。通过在球体上求解方程(1)，简并解的存在使得某些本征模具有相同的本征值和空间波长——这类似于量子物理中的球谐波。事实上，由于本征模在皮质折叠消失的极限下会趋近于球面谐波，因此前者可以组合成一个本征群，具有空间波长：

$$wavelength=\frac{2\pi R_s}{\sqrt{l(l+1)}}\text{\hspace{1em}(3)}$$
式中$R_s$为球的半径（对于本研究使用的群体平均模板，$R_s\approx 67.0㎜$），l为本征群数（原子物理中的角动量量子数）。

脑活动的模态分解

我们使用几何本征模分解时空fMRI数据，在空间位置r和时间t测量，为每个个体，作为模态的加权总和：

$$y(\mathbf{r},t)=\sum _{j=1}^N{a}_j(t){\psi }_j(\mathbf{r})\text{\hspace{1em}(4)}$$
其中，$a_j$为解释数据时j阶模态的幅值，${\psi }_j$为第j阶模态，N为所用模态数；我们使用N=200进行分析。对于时空数据，即无任务fMRI记录的自发动力学数据，将每个时间帧的数据代入公式(4)，得到每个模态${\psi }_j$的时间依赖性振幅$a_j(t)$。对于纯粹的空间数据，即任务诱发的激活图，其振幅与时间无关，如$a_j(t)\to a_j$。在这两种情况下，振幅都可以通过对皮质表面的积分得到：

$$a_j(t)=\int y(\mathbf{r},t){\psi }_j(\mathbf{r})d\mathbf{r}\text{\hspace{1em}(5)}$$
可由式(4)利用本征模的正交性导出。如果没有足够的测量来计算积分，振幅也可以通过统计一般线性模型估计。
得到振幅后，利用公式(4)对重建数据进行计算。我们通过计算实验数据和重建数据之间的相关性来量化这种重建的准确性。对于时空无任务数据，我们首先按照领域内惯例，对离散的区划内的数据取平均值，将实验数据和重建数据分割，然后通过计算成对区划时间序列的Pearson相关系数，构建区域间FC矩阵。对于任务诱发的数据，我们对激活图应用了相同的分割，以允许直接比较。最后，利用实验和重建FC上三角元素的相关性计算了无任务数据的重建精度。对于任务诱发的数据，重建精度是由实验的和重建的区划图的空间相关性来计算的。然后我们计算了它们的平均重建精度。

HCP数据

使用预处理后的 HCP fMRI 数据。我们没有执行任何额外的预处理步骤，例如全局信号去除，因为第一个本征模（被认为是全局、常数模态）已经显式地捕获数据中的全局偏差，从而允许其他模态捕获功能相关的非全局活动。我们分析了255名无亲缘关系的健康人（年龄22 ~ 35岁，132名女性，123名男性）的数据，这是排除双胞胎或兄弟姐妹后最大的HCP样本，所有参与者均完成了任务诱发和无任务静息态数据。所有参与者均为志愿者并签署知情同意书。开放获取数据由 WU-Minn HCP 联盟在获得当地监管伦理委员会批准的情况下获得，并根据HCP的数据使用条款与作者共享。我们的所有程序均按照这些数据使用条款设定的方案实施。
在HCP数据集中，我们分析了已经过HCP预处理的7个任务域的任务诱发fMRI和无任务静息态fMRI。两组数据均映射到 fsLR-32k CIFTI 空间上，每个半球有32,492个顶点。我们还分析了来自dMRI数据的单个连接组。连接组代表每个半球大小为32,492×32,492的高分辨率加权矩阵。

皮质区划

这项研究的主要结果分析了任务诱发的激活图和划分成离散区域的无任务FC数据。我们使用每个半球180个区域的HCP-MMP1区划（我们称之为Glasser360分区）来展示结果，这反映了基于皮质结构、功能、连接和解剖的清晰区域边界。区划在HCP提供的fsLR-32k空间上进行。为了测试我们结果的稳稳性，我们还使用了Schaefer等人在fsLR-32k空间上提供的不同分辨率（两侧大脑半球的100、200、400、600、800和1,000个区划，我们分别将它们称为Schaefer100、Schaefer200、Schaefer400、Schaefer600、Schaefer800和Schaefer1000区划）。

连接组本征模的推导

根据以往的方法推导出连接组本征模，以便与以往的研究结果进行比较。注意，在以前的研究中，连接组本征模被称为连接组谐波，但我们在这里使用本征模，因为谐波意味着整数频率比，这对于脑源性模态并不一定能保证。
我们获得了dMRI示踪成像测量的高分辨率连通性图，其中皮质表面网格中32,492个顶点的连通性通过追踪从每个点出发直到终止于另一个点的流线来估计。在不需要归一化的情况下，将顶点（视为节点）之间的连接权值估计为互连流线的数量。从对HCP个体进行纤维束成像的数据中，我们组合了大小为32,492×32,492的个体加权连接矩阵，以生成组平均连接组($W_{connectome}$)，其权重代表流线的平均数量。然后，我们生成了一个二进制邻接矩阵$A_{local}$，该矩阵捕捉了皮质表面网格模型（通过连接网格中直接相邻的两个顶点构建）中点之间的局部空间关系的离散表示。这些连接旨在捕捉传统dMRI示踪成像无法解决的局部、极短程连接。
组平均加权连接组($W_{connectome}$)被阈值化，以去除最小的权重，使连接数量大于$A_{local}$的4倍。将阈值化的矩阵再二值化，得到组邻接矩阵$A_{connectome}$。最后，我们生成了一个合并邻接矩阵$A_C=A_{local}||A_{connectome}$（矩阵大小32,492×32,492；||是逻辑或算子），它捕捉了局部顶点到顶点的连接，以及实验测量的复杂的短连接和长连接。注意，由上述阈值过程产生的邻接矩阵$A_C$的连接密度为0.10%。
为了得到连接组的本征模，求解特征值问题：

$$L^{\prime }\psi =-\lambda \psi \text{\hspace{1em}(6)}$$
其中L’是归一化的图Laplacian，是LBO的离散对应。归一化的图Laplacian算子与非归一化的图Laplacian算子L有关，$L^{\prime }=D^{-1/2}L{D}^{-1/2}$，根据之前的工作，L定义为：

$$L=\frac{1}{2}[(D-A_C)+(D-A_C{)}^T]\text{\hspace{1em}(7)}$$
其中D为对角线度矩阵，上标T为矩阵转置。与几何本征模一样，方程(6)的本征值排成序列$0\le {\lambda }_1\le {\lambda }_2\le \dots$。还需要注意的是，使用高分辨率的、顶点级的连接组可以得到的连接组本征模，横跨维数（模态的数量）等于顶点数的空间，允许与几何本征模进行公平的比较。

EDR本征模的推导

指数距离规则的本征模也通过求解方程(6)得到。然而，在这种情况下，非归一化的图Laplacian算子L，是使用一个人工构造的32,492×32,492邻接矩阵$A_E$来定义的，它遵循一个随机EDR。为了构建$A_E$，我们使用了上一节中的组平均、非阈值化加权连接组($W_{connectome}$)。然后，我们用$e^{-\alpha d}$形式的指数函数（其中α是一个尺度指数参数）拟合权值方差，作为皮质表面顶点之间欧氏距离d的函数。在 MATLAB 2019b 中使用非线性最小二乘法进行拟合，得到最优经验参数值${\alpha }_{empirical}=0.12$，与之前基于连接概率与距离函数的估计一致。然后，我们按照随机EDR布线过程生成一个随机的二进制邻接矩阵，其中两个顶点以概率$e^{-\alpha d}$连接，$\alpha ={\alpha }_{empirical}=0.12$。
EDR模型生成的邻接矩阵$A_E$的连接密度为1.55%，显著高于$A_C$的0.10%密度。因此，我们构建了另一个版本的连接组本征模（称为密度匹配的连接组本征模），它对组平均实验连接组进行阈值化，使最终获得的$A_C$具有与$A_E$匹配的密度，从而允许公平比较。

与统计基集的比较

几何、连接组和EDR本征模基集都来自一个生成模型，该模型解释了大脑功能如何从解剖学中产生。这种方法与文献中常用的统计基集不同，后者可以有效地总结数据，但不能洞察潜在的生成过程。我们评估了几何本征模在两个统计基集上的性能，一个来自功能数据本身的PCA，另一个基于Fourier空间基集。

任务诱发激活图的模态功率谱

为了研究任务诱发激活图的频谱内容，我们使用方程(4)中的模态分解并取振幅a的绝对平方来计算它们的模态功率谱。这类似于从Fourier分析中计算时间功率谱密度。然后我们将j模态的功率相对于所有模态的总功率归一化，这样：

$$P_j=\frac{|a_j{|}^2}{{\sum }_{j=1}^N|a_j{|}^2}\text{\hspace{1em}(8)}$$
我们计算了两组任务诱发数据的模态功率谱。第一组包括来自HCP任务诱发数据的非阈值激活图。在进行频谱分析之前，我们对255个个体的组平均激活图进行分析，然后分析每个任务对比的组平均激活图的功率谱。图3a显示了47个HCP任务对比图的平均功率谱。
第二组包括来自NeuroVault的1178个独立实验的10,000个激活图。我们使用python模块Nilearn从NeuroVault检索无阈值且带有fMRI-BOLD模态标签的激活图。我们使用Nilearn（通过函数Nilearn .surface.vol_to_surf）将容积激活图投射到 fsLR-32k CIFTI 空间，以匹配HCP数据。然后我们分析了每个激活图的功率谱。图3a中的结果显示了10,000个NeuroVault图的平均功率谱。
然后，我们比较了实验图的功率谱与具有不同平滑程度的人工合成随机图的功率谱，以进一步研究长波模态的相关性。特别地，我们在体积空间中生成了10,000个随机图，我们用尺寸在0到50㎜的全宽半最大值范围内的核对这些图进行了平滑，并将其投射到 fsLR-32k CIFTI 空间上。然后，我们计算了实验图和合成图之间功率谱的均方对数误差(MSLE)。

长波和短波模态的贡献

为了了解长波和短波几何本征模在任务诱发激活图重建中的作用，我们在进行任务诱发激活图重建之前先将其去除。具体来说，我们首先用200个几何本征模重建了7个关键的HCP对比图，并计算了重建精度（即实验图和重建图之间的相关性），作为我们的基线。然后，我们从长波模态开始进行递增的、顺序的模态移除（即移除模态1、模态1-2、模态1-3、模态1-4、…、模态1-200），并计算每一个台阶的重建精度。我们从短波模态开始重复相同的程序（即去除模态200、模态199-200、模态198-200、模态197-200、…、模态1-200)。

NFT波模型

如上所述，几何本征模的优越表现的一个暗示是，神经元活动由NFT预测的波动力学主导。为了研究波动力学是否能解释神经活动的复杂时空模式，我们实现了一个简单的NFT波模型，这种模型的动力学，由一个无再生的各向同性阻尼波动方程描述：

$$\frac{\mathrm{d}{S}_i}{\mathrm{d}t}=f(\mathbf{S},{\theta }_i,C,G)\text{\hspace{1em}(9)}$$
其中$\phi (\mathbf{r},t)$是位置r和时间t的神经活动，Q是外部输入，${\gamma }_s$是阻尼率，$r_s$是波传播的空间长度尺度（概念上与EDR本征模推导中随机EDR布线过程的α长度尺度有关）。这种形式告诉我们，脉冲输入将产生一种以${\gamma }_s$的速率耗散并以${\gamma }_s{r}_s$的速率传播的活动。在这里，我们将${\gamma }_s$作为一个固定参数，其值为$116s^{-1}$，而$r_s$作为一个自由参数。我们将模型应用于每个半球有32,492个顶点的皮质厚度中层表面网格，以求解每个顶点的活动。请注意，点与点之间的活动传播由它们的白质连接控制，其强度随距离呈近似指数衰减。将式(9)转换为其等效积分形式时，这种距离依赖性更加明显。

神经质量模型

我们比较了简单的NFT波模型和生物物理大尺度神经质量模型的性能，在该模型中，通过实验解剖连接耦合的神经群体（即神经质量）的相互作用产生了介观动力学。在一个典型的神经质量模型中，每个脑区i都有自己的平均场群体动力学，其时间活动$S_i$由通用方程定义：

$$\frac{d{S}_i}{dt}=f(\mathbf{S},{\theta }_i,C,G)\text{\hspace{1em}(10)}$$
其中f是描述区域活动演化的函数。该函数依赖于其他区域的活动S、局部群体参数${\theta }_i$、区域之间的解剖连接C（顶点分辨率连接组的一个区划版本，$W_{connectome}$，用于推导连接组的本征模）和缩放区域之间连接的全局耦合参数G。
文献中有几种我们可以使用的全脑神经质量模型，从简单的相位振荡器模型（如Kuramoto模型）到更复杂的生物物理群体模型（如Wilson-Cowan模型）。所有这些模型都遵循方程(10)的形式，特别是它们依赖于解剖区域间连接矩阵C。在这里，我们重点关注一种广泛使用的神经质量模型，BEI模型。BEI模型使用平均场方法来近似每个脑区的局部群体动态，通过来自dMRI的解剖连接矩阵进行耦合。每个脑区i由兴奋性(E)和抑制性(I)神经元相互作用的群体组成，受以下非线性随机微分方程控制：

$$\frac{\mathrm{d}{S}_i^{(E)}}{\mathrm{d}t}=-\frac{S_i^{(E)}}{{\tau }_E}+\left(1-S_i^{(E)}\right)\gamma r_i^{(E)}+\sigma {\nu }_i(t)\text{\hspace{1em}(11)}$$

$$\frac{\mathrm{d}{S}_i^{(I)}}{\mathrm{d}t}=-\frac{S_i^{(I)}}{{\tau }_I}+r_i^{(I)}+\sigma {\nu }_i(t)\text{\hspace{1em}(12)}$$

$$r_i^{(E)}=H^{(E)}\left(I_i^{(E)}\right)=\frac{a_E{I}_i^{(E)}-b_E}{1-\exp \left[-d_E\left(a_E{I}_i^{(E)}-b_E\right)\right]}\text{\hspace{1em}(13)}$$

$$r_i^{(I)}=H^{(I)}\left(I_i^{(I)}\right)=\frac{a_I{I}_i^{(I)}-b_I}{1-\exp \left[-d_I\left(a_I{I}_i^{(I)}-b_I\right)\right]}\text{\hspace{1em}(14)}$$

$$I_i^{(E)}\left(t\right)=I^{\mathrm{ext}}+W_E{I}_0+w_{EE}{S}_i^{(E)}\left(t\right)+GJ\sum _j{C}_{ij}{S}_j^{(E)}\left(t\right)-w_{IE}{S}_i^{(I)}\left(t\right)\text{\hspace{1em}(15)}$$

$$I_i^{(I)}\left(t\right)=W_I{I}_0+w_{EI}{S}_i^{(E)}\left(t\right)-S_i^{(I)}\left(t\right)\text{\hspace{1em}(16)}$$
其中$S_i^{(E,I)},r_i^{(E,I)}$和$I_i^{\left(E,I\right)}$分别代表E群和I群的突触门控变量、放电率和总输入电流。参数${\tau }_{E,I}$为时间常数，γ为动力学速率常数，$v_i(t)$为随时间变化的随机高斯输入，标准差σ。函数$H^{(E,I)}$是将总输入电流$I_i^{(E,I)}$转化为放电速率$r_i^{(E,I)}$的Sigmoidal神经元反应函数，由增益因子$a_{E,I}$，阈值电流$b_{E,I}$和曲率参数$d_{E,I}$进行参数化。在方程(15)和(16)中，$I_i^{\mathrm{ext}}$是外部输入，$I_0$是本地输入电流，由$W_E$和$W_I$缩放，分别代表兴奋和抑制性群体，$W_{EE}$是兴奋-兴奋强度，$W_{EI}$是兴奋-抑制强度，$W_{IE}$是抑制-兴奋强度，G是全局耦合参数，J是有效的N-甲基-D-天冬氨酸(NMDA)电导，$C_{ij}$是i和j区域之间的结构连接强度，根据dMRI估计。
总体而言，BEI模型有15个固定参数和4个自由参数。为了进行直接比较，我们还使用离散的连接组数据来定义区域间结构连接矩阵C。这些连接组数据来自334名无关的HCP受试者的最小预处理dMRI数据，并通过 FMRIB Software Library (FSL) 的概率示踪工具构建。当使用我们的连接组数据的区划版本时，我们的结果没有改变。我们注意到，NFT方程的数值解（如方程(9)）也在空间上离散了皮质，但在将每个点的时间动力学与连续极限中正确的连接性相结合之前，将其整合到一个非常精细的点阵列中。因此，神经质量模型近似于更一般的NFT方法。我们还重申，局部动力学不影响NFT波模型的空间特征函数。

血液动力学的模型

为了模拟fMRI数据，我们使用成熟的Balloon-Windkessel血流动力学模型，将NFT波和BEI神经质量模型产生的神经活动转换为血氧水平依赖(Blood Oxygen-Level Dependent, BOLD)信号。需要注意的是，虽然这个模型是BOLD信号下更详细的生理血流动力学过程模型的简单近似，但我们在这里使用这个近似是为了与文献中的绝大多数建模研究进行直接比较。每个顶点或脑区i的BOLD-fMRI信号由以下微分方程控制：

$$\frac{\mathrm{d}{z}_i}{\mathrm{d}t}=N_i(t)-\kappa z_i(t)-\gamma [f_i(t)-1]\text{\hspace{1em}(17)}$$

$$\frac{\mathrm{d}{f}_i}{\mathrm{d}t}=z_i(t)\text{\hspace{1em}(18)}$$

$$\frac{\mathrm{d}{v}_i}{\mathrm{d}t}=\frac{1}{\tau }[f_i(t)-v_i^{1/\alpha }(t)]\text{\hspace{1em}(19)}$$

$$\frac{\mathrm{d}{q}_i}{\mathrm{d}t}=\frac{1}{\tau }\left\{\frac{f_i(t)}{\rho }[1-{(1-\rho )}^{1/f_i(t)}]-v_i^{1/\alpha -1}(t)\right\}\text{\hspace{1em}(20)}$$

$$\frac{\mathrm{d}{y}_i}{\mathrm{d}t}=V_0\left\{k_1\left[1-q_i(t)\right]+k_2\left[1-\frac{q_i(t)}{v_i(t)}\right]+k_3\left[1-v_i(t)\right]\right\}\text{\hspace{1em}(21)}$$
其中$z_i,f_i,v_i,q_i,y_i$分别为血管舒张信号、入血量、血容量、脱氧血红蛋白含量和BOLD信号变量。变量$N_i$代表由波和BEI神经质量模型产生的神经活动。对于波模型，N(t)为ϕ(t)；对于BEI神经质量模型，$N_i(t)$为$S_i^{(E)}(t)$。模型参数及其值从以前的工作得到，如下：$\kappa =0.65s^{-1}$是信号衰减率，$\gamma =0.41s^{-1}$是流依赖的消除速率，τ=0.98s是血液动力学的穿行时间，α=0.32是Grubb指数，ρ=0.34是静息氧气提取分数，$V_0=0.02$是静息血液体积分数，$k_1=3.72,k_2=0.53,k_3=0.53$是 3T fMRI 参数。

静息态动力学建模

我们使用NFT波和BEI模型，以及血流动力学模型来评估各种自发FC特性。对于波模型，我们根据先前的研究，用白噪声输入来模拟没有任何结构性刺激的情况，求解方程(9)。我们将得到的结果与方程(17-21)中的血流动力学模型相结合来模拟BOLD-fMRI信号。将模拟BOLD信号降采样，采样间隔为0.72s，采样时间为1200帧，与静息态HCP数据进行匹配。最后，我们使用HCP-MMP1区划的方法对模拟的BOLD信号进行分区，并计算出代表模型FC的相关矩阵。
对于BEI模型，我们可以使用上述方法从方程(11)的每个脑区的数值解开始计算模型FC。然而，研究表明，我们可以将BEI和血流动力学模型的方程线性化，以获得模型FC的解析近似。由于BEI模型中有大量的模型参数，我们在本研究中使用了这种解析近似，因为它允许更全面和计算效率更高的模型拟合。
通过拟合每个模型的自由参数，我们将模型FC与所有分析中使用的相同HCP参与者的实验FC（也使用HCP-MMP1区划）进行了拟合。我们将参与者样本分为训练集和测试集，每组125人，以实现对模型性能的样本外评估和避免过拟合。具体来说，我们在训练集数据上拟合模型参数，并使用拟合的模型参数预测测试集数据。模型拟合和性能评估基于三个广泛使用的FC相关指标：边FC、节点FC和FCD。
利用z变换模型和实验FC的上三角元素（即FC边的强度）的Pearson相关性计算边FC度量值。我们只取上三角元素是因为FC值对于对角线是对称的。相关性越高，表示模型与数据的拟合度越好。
采用模型中各脑区平均FC强度与实验FC的Pearson相关性计算节点FC度量值。定义一个脑区i的平均FC强度为$\frac{1}{n}{\sum }_{j=1}^n{\mathrm{F}\mathrm{C}}_{ij}$，其中n为脑区总数。同样，相关性越高，模型拟合度越好。
FCD指标捕捉了静息态活动的时空统计量，计算方法如下。每个区域i的时间序列在0.04~0.07Hz之间使用二阶Butterworth滤波器进行滤波，选取这个频段是由于它与大脑的功能相关。然后对时间序列进行 Hilbert (H) 变换，计算$y_i\left(t\right)=x_i\left(t\right)+j{H}_i(t)$，其中j为虚数。然后计算瞬时复参量${\theta }_i\left(t\right)={\tan }^{-1}\left[H_i(t)/x_i(t)\right]$。区域i和区域j在t时间点的同步水平Δ(i,j,t)计算如下：

$$\Delta \left(i,j,t\right)=\cos \left[{\theta }_i\left(t\right)-{\theta }_j(t)\right]\text{\hspace{1em}(22)}$$
请注意，我们没有像有些文献中那样将${\theta }_i\left(t\right)$视为相位，因为我们的信号是宽带的，而这种解释只适用于几乎是单频的信号。然后计算了两个时间实例${\tau }_u$和${\tau }_v$的全局同步度的相似度${\phi }_{uv}$，定义为：

$${\varphi }_{uv}=\frac{1}{d_u{d}_v}\sum _{i>j}\Delta (i,j,{\tau }_u)\Delta (i,j,{\tau }_v)\text{\hspace{1em}(23)}$$
其中

$$d_x=\sqrt{\sum _{i>j}{[\Delta (i,j,{\tau }_x)]}^2}\text{\hspace{1em}(24)}$$
这里${\phi }_{uv}$是FCD，它是一个对称的时间×时间矩阵，解释为${\tau }_u$和${\tau }_v$之间全局同步的相似性。然后，我们比较了模型的上三角元素和实验FCD评估值的分布，将所有个体或模型实现情况串联起来，以更好地捕捉数据的总体波动。采用KS统计量比较模型与数据的FCD分布，KS统计量越低，模型拟合越好。
波模型有一个自由参数，即波传播的空间长度尺度($r_s$)，经过优化以拟合实验fMRI数据。具体来说，我们构造了一个向量，其中20个$r_s$值均匀地间隔在10到100㎜之间。对于每个$r_s$值，计算上述FC指标。我们取了使 FCD KS 统计量最小化的$r_s$值，因为已经发现这个指标是上面讨论的三个指标中模型-数据比较最严格的基准。从这一步骤得到了$r_s=28.9㎜$的优化值，这比之前在脑电数据分析中获得的值要小，可能是由于血流动力学过程限制了神经活动向邻近区域的传播，或者我们这里关注的是皮质-皮质动力学。
BEI神经质量模型有15个固定参数和4个自由参数($w_{EE},w_{EI},w_{IE},G$)，经过优化以拟合数据。优化后的参数$w_{EE}=9.80,w_{EI}=1.48,w_{IE}=7.13,G=6.87$。

静息态动力学时滞特性的测量

除了用于模型拟合和模型性能评估的基于FC的指标外，我们还研究了波和神经质量模型是否也可以捕获传播活动的时间特性。具体来说，我们分析了静息态BOLD-fMRI时间进程的滞后结构（或滞后线程）。请注意，还有其他几种方法可以表征静息态活动的时空特性，但目前选择的方法对于我们的目的已经足够了。

刺激诱发的动力学建模

我们还评估了波模型捕捉诱发神经反应经典特性的程度。具体来说，我们使用了优化的波模型($r_s=28.9㎜$)来研究刺激左V1时的神经动力学。具体来说，刺激$Q(\mathbf{r},t)$是一个1㎳的脉冲(t =1-2㎳)，幅度为$20s^{-1}$（结果对幅度的变化是不敏感的），限制在皮层表面的顶点，这些顶点位于HCP-MMP1区划定义的V1区域内。我们在100㎳的时间内进行了模拟，分辨率为0.1㎳。最后，我们使用HCP-MMP1区划将活动$\phi (\mathbf{r},t)$划分到180个区域。
我们比较了每个脑区的活动情况（振幅与时间），重点关注活动达到峰值幅度的时间（简称达峰时间）。具体地，我们研究了活动的时间顺序是否遵循人类视觉皮层的层次，从视觉皮层到额叶皮层，包括以下17个脑区：V1, V4, 7m, 7Am, TE1p, 7AL, 24dd, 2, 24dv, 8BM, 10r, 10v, 8BL, 10pp, 10d, 9-46v 和9-46d。这些区域与之前使用轨迹追踪数据和非线性网络模型在猕猴新皮层视觉层次中发现的区域非常相似。
我们还将峰值响应时间的区域估计值与T1w:T2w值进行了比较，T1w:T2w值是一种对皮质内髓鞘含量敏感的非侵入性测量方法，也是皮质层次排序的良好替代指标。从HCP数据集获得fsLR-32k空间的髓鞘图谱，然后使用HCP-MMP1区划进行划分。我们通过Spearman位序相关量化了局部达峰时间值与髓鞘含量之间的关系。相关性的统计显著性通过与使用空间约束自旋测试方法获得的10,000个相关值的零分布进行比较来评估。这种方法计算一张图和另一张图的随机空间旋转版本之间的相关性，从而保持图中区划的空间关系。得到的P值，即$P_{spin}$，是零相关值大于实验相关值的频率。最后，我们重复了上述包括所有脑区（即不限于视觉层次脑区）的统计检验，以确保我们的发现不是由我们对感兴趣区域(ROI)的特定选择驱动的。这是一个更为保守的测试，因为并非所有的大脑区域都对视觉刺激表现出强烈的诱发反应。

非新皮质结构的几何本征模的评估

我们将本征模分析扩展到新皮质之外的区域，特别关注皮质下（丘脑和纹状体）和旧皮质（海马）。与皮质带可以被模拟成二维薄片不同的是，这些结构是立体的三维物体。因此，我们使用四面体网格而不是基于表面的三角形网格来计算几何本征模，以考虑非新皮质结构的完整3D几何形状，如下所述。
我们首先使用Harvard-Oxford皮质下概率图谱在每个半球生成丘脑、纹状体和海马的体积二元蒙版。属于这些结构的概率在25%以上的体素被包含在蒙版中。我们使用FreeSurfer的mri_mc函数（该函数实现了一个移动立方体算法），通过对体积蒙版进行镶嵌来构建一个2D表面，然后使用Gmsh软件将2D表面转换为3D四面体网格。然后我们使用 LaPy python 库在每个非新皮质结构的四面体网格上求解方程(1)，得到本征模及其对应的特征值。最后，我们通过插值将四面体空间中的本征模投影回自然体积空间。因此，产生的本征模在空间上随包含每个非新皮质结构体积的3D体素变化。在本部分研究中，我们丢弃了第一个常模态，并在分析中使用了接下来的20个模态。

描绘非新皮质结构的功能组织

fMRI非新皮质结构的信噪比普遍弱于新皮质。此外，由于fMRI分辨率的限制以及常规fMRI处理方法引入的空间平滑性，精细粒度的任务激活往往难以在这些较小的结构中进行解析。因此，为了有效地映射丘脑、纹状体和海马的功能组织，我们使用了静息态fMRI信号的连接组绘图来获得它们的主要功能模态（通常称为梯度，尽管该模式捕捉了FC谱的点相似性的空间变化）。这种技术和相关的程序在过去的工作中被广泛地用于研究这些结构的功能组织。 我们将连接组绘图应用于HCP个体的容积体素静息态fMRI数据。具体来说，对于每个个体和非新皮质ROI（即丘脑、纹状体和海马），我们构建了一个ROI时间序列数据矩阵，A，大小为T×N，其中T是时间帧的数量，N是ROI中的体素数量。类似地，我们为每个个体构建了大小为T×M的灰质时间序列数据矩阵B，其中M是ROI外的灰质体素的数量。因为M一般较大（1000以上），我们使用奇异值分解来降低B的维数，以构造矩阵$\widetilde{B}$，尺寸T×(T-1)。然后使用A和$\widetilde{B}$的Pearson相关性来计算ROI内的每个体素的连接指纹，得到的大小为N×(T-1)的矩阵C。利用$η^2$系数计算每对ROI体素之间连接指纹的相似性，得到大小为N×N的矩阵S。$η^2$系数代表一个连接谱的方差占另一个连接谱的方差的比例。然后取所有个体S的平均值。利用Laplacian特征映射算法对每个ROI的S矩阵进行非线性流形学习，计算其特征向量和特征值。特征向量代表ROI的功能模式，称为功能梯度，根据它们解释的FC相似性的方差排序。我们只分析了前20个非常数梯度，以便与几何本征模进行直接比较。基于FC的梯度分析的典型应用很少考虑超过前5个梯度。 由于模态和梯度的符号是任意的，我们在分析每个非新皮质结构的几何本征模和功能梯度之间的对应关系时，取它们的绝对空间相关性。此外，同一特征群内的几何本征模序是可以改变的（例如，第一个特征群可以在2-4模态之间发生顺序翻转），因此不能保证几何本征模的序号与函数梯度一一对应。因此，我们检查了模态和梯度之间所有可能的成对相关性，并评估了每个梯度与最大相关的模态的对应关系。丘脑、纹状体和海马的几何模态和功能梯度的最大阶差分别为1、8和5（共20阶）。

报告总结

关于研究设计的更多信息可在与本文链接的自然作品集报告总结中获取。

获取数据

原始和预处理的HCP数据
NeuroVault的数据
复现研究结果的源数据可以在https://github.com/NSBLab/BrainEigenmodes和https://osf.io/xczmp/上公开获取。

获取代码

用于计算本征模、分析结果和再现研究数据的计算机代码

引用