NC：皮层微结构的神经生理特征

在整个皮层中观察到微结构的系统空间变化。这些微结构梯度反映在神经活动中，可以通过神经生理时间序列捕获。自发的神经生理动力学是如何在整个皮层组织的，以及它们是如何从异质皮层微结构中产生的，目前尚不清楚。在这里，我们通过估计来自静息状态脑磁图(MEG)信号的6800多个时间序列特征，广泛地描绘了整个人脑的区域神经生理动力学。然后，我们将区域时间序列概况映射到一个全面的多模式，多尺度的皮质微结构图谱，包括微观结构，代谢，神经递质受体，细胞类型和层流分化。我们发现神经生理动力学的主导轴反映了信号的功率谱密度和线性相关结构的特征，强调了电磁动力学的常规特征的重要性，同时识别了传统上较少受到关注的附加信息特征。此外，神经生理动力学的空间变化与多种微结构特征共定位，包括基因表达梯度、皮质髓鞘、神经递质受体和转运体、氧和葡萄糖代谢。总的来说，这项工作为研究神经活动的解剖学基础开辟了新的途径。

信号以电脉冲的形式，通过多种类型的神经元和神经元群不断地产生、传播和整合。大脑的线路引导信号通过嵌套的多功能神经回路网络传播。由此产生的膜电位和放电率的波动最终表现为神经生理活动的模式。

大量文献证明了皮层微结构与动力学之间的联系。许多研究调查了皮层节律的细胞和层状起源。例如，脑电图和脑磁图(EEG/MEG)信号似乎对来自皮质层II-III和V14、15的锥体细胞的偶极子更敏感。此外，神经元电生理的特定时间序列特征取决于神经元类型、形态和局部基因转录，特别是与离子通道调节相关的基因。然而，以前的研究大多集中在单个或小的感兴趣的特征集上，经常将单个微架构特征映射到单个动态特征。从脑电图中发现8 - 12hz的α节律开始，神经科学中传统的时间序列分析通常集中在典型的电生理节律上。最近，人们对内在时间尺度的研究也越来越感兴趣，这些内在时间尺度显示了从单峰皮层的快速波动活动到跨峰皮层中较慢的上下文信息编码的时间处理层次。神经回路微结构的特定特征如何产生持续的神经生理动力学仍然是神经科学中的一个关键问题。

最近的分析进展为通过计算超越功率谱测量的综合特征集(包括信号幅度分布、熵、分形缩放和自相关的测量)来执行神经生理时间序列表型提供了新的机会。随之而来的成像技术和数据共享的进步提供了前所未有的细节和深度测量大脑结构的新方法，包括基因表达、髓鞘形成、神经递质受体、细胞结构、层流分化、细胞类型组成、代谢和进化扩张。

在这里，我们全面表征神经生理活动的动态特征，并通过整合人类皮层的多个多模态图将其与潜在的微结构联系起来。我们不是手动选择少量感兴趣的特征，而是使用数据驱动的方法使用大量的动态和微架构特征。具体来说，我们首先使用来自人类连接组计划(HCP)。然后，我们应用高度比较的时间序列分析(hctsa)来估计每个大脑区域的综合时间序列特征(图1)。同时，我们构建了皮层的微结构图谱，其中包括微观结构、代谢、神经递质受体和转运体、层流分化和细胞类型(图2)。最后，我们使用多元统计分析将这些广泛的微结构和动态图谱相互映射。

图片

图1 基于特征的神经生理时间序列表示。

图片

图2 多模态脑图

1. 结果

通过将线性约束最小方差(LCMV)波束形成应用于来自人类连接组计划(HCP)，使用Brainstorm软件。高度比较时间序列分析(hctsa)，然后从区域MEG记录中进行大量时间序列特征提取。该过程提供了一种基于特征的时间序列表示，其中给定的时间序列由时间序列特征向量表示。这种时间序列表型分析是一种数据驱动的方法，它使用信号的时间结构的跨学科度量来量化神经活动的动态曲目，并产生每个大脑区域动态特性的全面“指纹”。将时间序列表型应用于区域MEG时间序列，我们从Schaefer-100图谱中估计了100个皮质区域的6880个时间序列特征。hctsa库包含大量跨学科的特征集，这些特征具有潜在的相关值，这些值跨越了各种概念时间序列特征。时间序列特征列表包括但不限于自相关函数、功率谱、幅度分布和熵估计得出的统计数据(图1)。

为了估计一组全面的多模态微结构特征，我们使用了最近开发的neuromaps工具箱、BigBrainWarp工具箱、艾伦人脑图谱(AHBA)和abagen工具箱来转换和编译一组45个特征。包括微观结构、代谢、皮质扩张、受体和转运体、层厚度和细胞类型特异性基因表达的测量(图2)。请注意，微观结构图包括基因表达和神经递质谱的主要梯度，因为它们各自代表了某些分子特性的代理测量。具体来说，基因表达的主要成分(基因表达PC1)提供了细胞类型在皮层分布的潜在代理，神经递质受体和转运体的主要成分(神经递质PC1)提供了18种神经递质受体和转运体的蛋白质密度的综合测量。我们还纳入了个体神经递质受体和转运蛋白图谱以及细胞类型特异性基因表达图谱，分别评估它们的作用。

在随后的分析中，我们首先通过量化静息状态MEG时间序列特性变化的主要模式来评估神经生理动力学的地形组织。然后，我们描述了神经生理动力学的特征，涉及到整个皮层的微观结构属性。最后，我们进行敏感性分析，以调查混杂因素对结果的潜在影响，如信噪比和包裹分辨率(详见敏感性分析)。

1.1 神经生理动力学的地形分布

hctsa时间序列表型过程在每个脑区生成6880个时间序列特征。由于hctsa包含用于量化任何给定时间序列属性的多个算法变体，因此确定的时间序列特征可能捕获相关的动态行为，并包括相关属性组。因此，我们首先试图利用主成分分析(PCA)来确定神经生理动力学的主要宏观模式或梯度。将PCA应用于组平均区域x特征矩阵，我们发现单个主导分量捕获了区域时间序列特征中48.7%方差的证据(图3a)。神经生理动力学(PC1)的主导成分或“梯度”主要横跨后顶叶皮层和感觉运动皮层的一端，以及颞叶前部、眶额和腹内侧皮层的另一端(图3a)。关注内在功能网络，我们发现主导神经生理动力学的地形组织沿着感觉-神经轴变化，从背侧注意、躯体运动和视觉网络到边缘和默认模式网络(图3a)。

图片

图3 神经生理动力学的地形分布

接下来，我们利用每个原始特征图和PC1图(即PCA加载)之间的单变量相关性，研究了第一个分量上的top-loading时间序列特征。使用保留空间自相关的零模型对所有相关性进行统计评估("自旋试验")。图3b显示，许多特征与PC1呈正相关或负相关;完整的特征列表及其相关系数和p值可在在线补充数据集S1中获得。对顶部荷载特征的检查表明，大多数是从功率谱结构或密切相关的测量中得出的统计数据。例子包括不同频段的功率，各种模型的参数与功率谱的拟合，以及相关的测量，如自相关函数的形状和波动分析测量。图3b显示了大脑皮层的功率谱变化情况，每条线代表一个大脑区域。区域根据其在假定的单峰-多峰等级中的位置着色;这种变化在视觉上表明单峰区域显示出更突出的α (8-12 Hz)和β (15-29 Hz)功率峰值。总的来说，这些结果表明，电生理时间序列分析的传统重点是功率谱的统计，这与不同hctsa时间序列特征库捕获的MEG动态的主要变化是一致的。

鉴于PC1的地形组织与功率谱特征密切相关，我们直接测试了PC1与常规带限功率谱测量以及固有时间尺度之间的联系(补充图1)。图3c显示了PC1与delta (2-4 Hz)、theta (5-7 Hz)、alpha (8-12 Hz)、beta (15-29 Hz)、lo-gamma (30-59 Hz)和hi-gamma (60-90 Hz)功率图以及固有神经时间尺之间的相关性。我们发现PC1与内在时间标度和高gamma显著相关。当我们使用针对功率谱的非周期分量而不是总功率进行调整的带限功率图时，结果是一致的 (补充图2)。PC1与固有时间标度相关的事实与两者在功率谱中捕获广泛变化的概念是一致的。鉴于固有时间尺度反映了功率谱的非周期分量的特征(这些度量在数学上是相关的;详细信息参见方法)，我们还直接评估了PC1与非周期分量的指数和偏移量之间的关系。指数描述了非周期分量的“曲线”或整体“线”或斜率，偏移量描述了整个功率谱的整体垂直位移(上下平移)。PC1与这两种测量结果均显著相关，这表明PC1捕获的时间序列特征也反映了功率谱的非周期分量的特性(补充图3)。

此外，先前的报告表明，宽带gamma活动也部分反映了非周期性神经生理活动和功率谱的宽带变化。这与我们的研究结果一致，即PC1与gamma功率和内在时间尺度有关，主要捕获功率谱和非周期活动特征的广泛变化。

请注意，我们专注于PC1，因为其他组件(PC2及以上)占时间序列特征方差的10%或更少，并且与hctsa时间序列特征没有显着关联。此外，为了验证数据的明显低维性和识别的PC模式不是由样本数量较少(即大脑区域)而不是特征(即时间序列特征)驱动的，我们进行了敏感性分析，我们随机选择了100个时间序列特征(从原始的6,880个特征列表中)并重新运行PCA(重复1 000次)。识别的PC模式及其相应的方差解释量与使用全套时间序列特征的原始分析一致(补充图4)。

1.2    微结构的神经生理特征

区域神经生理时间序列特征如何映射到多模态微结构特征上?为了解决这个问题，我们实施了多元偏最小二乘分析(PLS)，将多个多模态脑图集成到分析中，并寻求识别时间序列特征的线性组合和相互最佳协变的微结构特征的线性组合。图4a显示，分析确定了多个这样的组合，称为潜在变量(使用稀疏典型相关分析(sCCA)获得了类似的结果;补充图5)。利用空间自相关保持置换检验评估每个潜在变量的统计显著性。第一个潜在变量具有统计学意义，捕获了时间序列和微结构特征之间最大的协方差。

图片

图4 微结构的神经生理特征

图4b显示了时间序列特征的空间地形和第一个潜在变量的微结构得分。这些是原始输入特征的加权和，根据潜在变量确定的权重。分析使评分图之间的相关性最大化。因此，我们试图估计时间序列和微结构特征之间的相同映射是否可以在样本外观察到。我们采用了距离相关的交叉验证程序，随机选择“种子”区域，并选择75%物理上最近的区域作为训练集，而其余25%物理上最远的区域作为测试集(详见方法)。对于每个训练-测试分割，我们将PLS模型拟合到训练集，并将测试集投影到从训练集导出的权重(即奇异向量)上。然后将得到的测试集分数进行相关，以估计样本外相关系数。图4b显示，微结构和时间序列特征得分在训练集(平均rs = 0.75)和测试集(平均rs = 0.5)中是相关的，表明样本外分析的结果是一致的。

接下来，我们检查了相应的时间序列和微结构特征负载，并确定了对第一个潜在变量捕获的空间模式贡献最大的特征(图4c, d)。顶部加载时间序列特征主要与MEG信号的自相关或可预测性测量有关。自相关测度多反映神经生理时间序列的线性相关结构，特别是长滞后自相关。广泛的其他高度加权时间序列特征捕获了信号可预测性的其他方面，包括自相关函数形状的度量(例如，自相关函数越过零的时间滞后)，自相关结构在去除低阶局部趋势后的变化情况(例如，拟合线性模型到滚动5时间步或10 ms窗口的残差)，使用波动分析评估缩放特性(例如，信号方差在时间尺度上的缩放)，以及从小波分解(例如，不同时间尺度上的小波系数)得到的度量。

为了说明高贡献时间序列特征的空间分布，图4c显示了三个顶部负载特征，它们反映了第一个PLS潜在变量的空间变化。例如，图4c(左)描述了自相关函数的组平均第一个过零点的分布。单峰皮质的自相关函数(用粉色圆圈标记)与跨峰皮质(用紫色圆圈标记)的自相关函数以较低的滞后过零，表明跨峰皮质的自相关衰减更快，相关长度比单峰皮质长。另一个例子是MEG信号在较长时间滞后时的线性自相关。图4c(左)显示了延迟48 ms(24个时间步)的自相关性，显示单峰皮质的自相关性较低，而跨峰皮质的自相关性较高。请注意，顶部加载特征列表包括其他时间滞后的线性自相关，48 ms滞后的自相关只是作为一个说明性的例子(补充图6描述了hctsa中包含的所有时间滞后的线性自相关加载的全部范围)。最后，我们检验了用去趋势波动分析估计的标度指数α，作为对去趋势信号在时间尺度上波动的对数-对数图的线性拟合的斜率。图4c(右)描述了这种跨皮层的缩放指数，其表现出与前两个示例相似的空间模式，表明单峰皮层(粉色圆圈)的自相关性低于跨峰皮层(紫色圆圈)。波动分析的其他变化在顶载特征列表中也占有重要地位，包括线性拟合的拟合优度，多种标度制度的拟合，以及波动大小的不同类型的去趋势和数学公式。

图4d显示了相应的微结构负载。对第一潜变量捕获的空间格局贡献最大的微结构特征是基因表达的主要成分(基因表达PC1;一个潜在的代表细胞类型分布)，T1w/T2w比率(代表皮质内髓鞘)，神经递质受体和转运体的主要成分(神经递质PC1)，以及氧和葡萄糖代谢(强正负荷)。我们还发现特定的神经递质受体和转运体，特别是代谢性血清素能和多巴胺能受体，以及参与髓鞘形成的少突胶质细胞前体(opc)的细胞类型特异性基因表达的高贡献(强负负荷)。当我们使用单变量分析将区域时间序列特征与顶部加载微结构图，特别是基因表达的主成分和T1w/T2w比率联系起来时，得到了一致的结果，这些数据之前被广泛研究为典型的微结构梯度 (补充图7)。这种分析提供了一个全面的图表或“查找表”，说明微结构和时间序列特征图是如何相互关联的。这些结果表明，多种微结构属性的皮质变化表现为神经生理活动的时间序列特性的梯度，特别是反映信号的远程自相关结构的特性。

2. 讨论

在本研究中，我们利用时间序列表型分析，从静息状态MEG信号中全面绘制神经生理活动的动态指纹图谱。然后，我们将得到的动态图谱映射到多模态微结构图谱，以识别皮层微结构的神经生理特征。研究表明，神经生理时间序列特性的皮层变化主要反映功率谱密度，并与信号的内在时间尺度和自相关结构密切相关。此外，神经生理动力学的空间组织遵循神经递质受体和转运体、基因表达和T1w/T2w比值等微结构梯度，反映皮层代谢需求。

之前有大量研究调查了神经振荡及其与大脑神经通讯模式的关系。先前的报告也表明，神经振荡影响行为和认知，并涉及多种神经系统疾病和障碍。神经振荡表现为神经生理信号在频域的功率幅值变化。神经生理信号的功率谱特征，如标准频段的平均功率幅值，已被用于研究大规模大脑活动的潜在机制，并更好地了解大脑功能的个体差异。与功率谱密度相关的其他时间序列特性也被用于研究神经动力学，包括信号的固有时间尺度和自亲和或自相似性的度量(例如，自相关和波动分析)。

应用数据驱动的时间序列特征提取分析，我们发现神经生理时间序列特征的地形组织遵循感觉-神经轴，将躯体运动皮层、枕部皮层和顶叶皮层与颞前部皮层、眶额皮层和腹内侧皮层分离开来。这种神经生理活动的动态指纹主要是通过hctsa时间序列特征捕获的MEG信号的线性相关结构来表征的。线性相关结构既表现在功率谱性质上，也表现在自相关函数上。时间序列特征的这种主要空间变化也类似于内禀时间标度的空间分布，内禀时间标度是另一种与功率谱密度特征相关的度量。总之，虽然研究结果强调了代表性不足的时间序列特征，但它们强调了传统方法在表征神经生理活动方面的重要性，以及线性相关结构在脑电信号动力学中的关键作用。

早期的报告发现，区域神经动力学，包括功率谱和固有时间尺度的测量，反映了潜在的环路特性和皮层微结构。神经动力学和皮层微结构之间的关系通常用一个或几个微观结构特征来研究。数据收集和集成的最新进展以及越来越多的数据共享举措为使用广泛的多模态数据集全面研究皮质回路特性和微结构提供了独特的机会。在这里，我们使用这些数据集并编制多个微结构图，包括微观结构、代谢、皮层扩张、受体和转运体、层厚度和细胞类型特异性基因表达的测量，以绘制神经生理动力学和皮层微结构之间的多元关联。

我们的发现建立在先前的报告上，表明神经生理动力学遵循潜在的细胞结构和微观结构梯度。特别地，我们的研究结果证实了脑磁图内在动力学与基因表达和皮质内髓磷脂以及神经递质受体和转运体的异质性分布有关。此外，我们将正在进行的神经生理活动的动态特征与多种代谢属性联系起来;例如，我们发现氧气和葡萄糖代谢更大的区域往往表现出较低的时间自相关性，因此更多的可变瞬间内在活动。这与先前报道的感觉皮层氧气和葡萄糖消耗的高代谢率是一致的。我们还发现少突胶质细胞前体(opc)与细胞类型特异性基因表达的显著关联，可能反映了这些细胞在发育过程中通过产生髓鞘化少突细胞对髓鞘生成的贡献，以及在成人神经回路中髓鞘化轴突的髓鞘调节和代谢支持。最后，我们发现神经活动的主要动态特征与颗粒状皮层第IV层共同变化，这与第IV层接收突出的皮层下(包括丘脑)前馈投影的观点一致。总的来说，我们的发现建立在新兴文献的基础上，这些文献是关于异质微结构特性与宏观尺度网络嵌入(例如，皮质-皮质连接和皮质下投射)如何共同塑造区域神经动力学的。

总之，使用数据驱动的方法，目前的研究结果表明，皮层微结构的神经生理特征在皮层上分层组织，反映了潜在的环路特性。这些发现强调了研究神经生理活动特征的传统方法的重要性，同时也确定了与皮层微结构共变的不太常用的时间序列特征。总的来说，这项工作为研究神经生理活动的解剖学基础开辟了新的途径。

3. 方法简述

3.1 功率谱分析

使用Welch的方法从每个个体的源级时间序列中估计功率谱密度(PSD)，使用长度为4秒的重叠窗，重叠率为50%。然后为每个顶点(即源)计算每个频带的平均功率，作为给定频带频率范围内的平均功率。源级功率数据使用Schaefer-100 atlas在六个典型电生理波段(即δ (δ: 2-4 Hz)， θ (θ: 5-7 Hz)， α (α: 8-12 Hz)， β (β: 15-29 Hz)，低gamma(低γ: 30-59 Hz)和高gamma(高γ: 60-90 Hz))被分割成100个区域。我们将fsLR4k表面上的组平均顶点级功率图贡献给了公开可用的neuromaps工具箱。

3.2 固有时间尺度

区域内禀时间尺度是使用带有FOOOF(拟合振荡& 1/f)工具箱的频谱参数化来估计的。具体来说，源级功率谱密度用于提取每个顶点和每个个体的神经时间尺度，使用参考文献30中描述的程序。FOOOF算法通过拟合功率谱密度，将功率谱分解为周期(振荡)分量和非周期(1/f-like)分量。该算法识别振荡峰值(周期分量)、控制非周期分量弯曲的“knee参数”k和非周期“指数”χ。然后使用knee参数k来计算“knee频率”，如

图片

，这是功率谱密度中knee或弯曲发生的频率。最后，本征时间标度τ估计为：

图片

3.3 基于hctsa的时间序列特征提取

我们使用高度比较性的时间序列分析工具箱hctsa，对每个参与者的每个大脑区域的预处理时间序列进行大量特征提取。hctsa包使用广泛的时间序列分析方法提取了超过7000个本地时间序列特征。提取的特征包括但不限于数据分布的度量、时间依赖性和相关性、熵和可变性、时间序列模型拟合的参数以及给定时间序列的非线性特性。对每个参与者的分割后的MEG时间序列进行hctsa特征提取分析。考虑到在整个时间序列上应用hctsa的计算成本很高，我们在放弃前30秒后使用80秒的数据进行特征提取。先前的报告表明，相对较短的大约30到120秒的静息状态数据片段足以估计大脑内在活动的强大特性。然而，为了确保我们能够从80秒的数据中稳健地估计时间序列特征，我们使用catch-22工具箱在每个参与者的不同长度的时间序列的后续片段上计算了hctsa特征的子集。具体来说，我们以5秒的增量从5到125秒的短段数据中提取时间序列特征。为了确定估计稳健和稳定特征所需的时间序列长度，我们计算了两个后续片段(例如，从10秒和5秒数据估计的特征)之间的Pearson相关系数。估计特征之间的相关系数在30 s左右的时间序列段开始趋于稳定，与之前的报道一致 (补充图12)。从80秒的时间序列片段中提取特征后，删除产生错误的操作输出，并对每个参与者分别使用离群鲁棒s型变换在节点上对剩余的特征(6880个特征)进行归一化。将归一化的个体水平特征生成群平均区域x特征矩阵。我们还将hctsa分析应用于分割的empty-room记录(80秒)，以使用与静息状态数据相同的程序从噪声数据中估计时间序列特征，确定每个区域每个参与者的6148个特征。时间序列特征在每个参与者的大脑区域被标准化。获得了组平均empty-room特征集，并用于进一步分析。

3.4 neuromaps的微观结构特征

我们使用neuromaps工具箱(https://github.com/netneurolab/ neuromaps)在地图的原生空间中获得微结构和神经递质受体和转运体地图。有关所有地图及其数据源的详细信息可在43中找到。简单地说，在任何表面空间中原始可用的所有数据都使用线性插值方法转换到fsLR32k表面空间进行重新采样，并使用fsLR32k空间中的Schaefer-100地图集划分为100个皮质区域。所有体积数据保留在其原始的MNI152体积空间中，并使用MNI152空间中的体积Schaefer地图集划分为100个皮质区域。微结构图包括T1w/T2w作为皮质髓磷脂、皮质厚度、基因表达主成分、神经递质受体和转运体主成分、突触密度(使用[11C]UCB-J PET示踪剂结合突触囊泡糖蛋白2A (SV2A))、代谢(即脑血流量和体积、氧和葡萄糖代谢、糖酵解指数)、进化和发育扩张、来自费城神经发育队列(PNC)和美国国立卫生研究院(NIH)的异速测量量表。神经递质图谱包括9种不同神经递质系统中的18种不同的神经递质受体和转运体，即5-羟色胺(5-HT1a, 5-HT1b, 5-HT2a, 5-HT4, 5-HT6, 5-HTT)，组胺(H3)，多巴胺(D1, D2, DAT)，去甲肾上腺素(NET)，乙酰胆碱(α4β2, M1, VAChT)，大麻素(CB1)，阿片(MOR)，谷氨酸(mGluR5)和GABA (GABAa/bz)。

参考文献：Neurophysiological signatures of cortical micro-architecture.