神经发育过程中大脑临界状态的图神经网络分析方法

摘要

最近的研究表明，静息态功能动态可通过临界点附近的晶格模型(如二维Ising模型)进行建模。Ising温度是决定模型相变的控制参数，能够为大脑宏观动力学提供重要见解，现已被用于分析不同脑状态及神经发育过程。神经发育阶段以微观神经环路重组为特征，这些变化会影响宏观大脑动力学及其功能连接，并可通过功能磁共振成像(fMRI)进行观察。因此，本研究提出了一种新的方法，利用功能连接和图神经网络(在Ising模型网络上训练)，从fMRI数据中估算大脑的Ising温度以探究神经发育。研究结果显示，典型发育儿童(r=-0.48，p＜0.0001)和注意缺陷多动障碍儿童(r=-0.49，p＜0.0001)的年龄与温度之间均存在显著负相关，这表明随着年龄的增长，大脑逐渐远离临界态并向更有序的状态演化。

引言

大脑是由大量神经元构成的复杂网络，受非线性生化相互作用的调控并动态响应环境信息。这种复杂结构产生了从微观到宏观的时空模式，以及认知、情感、运动和感觉等行为。物理学家试图通过建立方程系统来模拟这种自组织动力学，而统计物理与复杂系统理论为此提供了重要工具。最突出的成果之一是：从fMRI BOLD信号中提取的静息态脑网络，与二维Ising模型(统计物理中研究相变的计算模型)提取的网络在统计学上没有显著差异，这表明宏观大脑动力学在功能上接近二维Ising模型的二阶相变临界点。该模型通过集体微观相互作用(由控制参数Ising温度调控)来解释自发磁化等涌现现象。这些涌现模式出现在控制参数的特定临界值区间，此时系统会在多种状态之间发生转换，由此提出了“大脑临界性”概念。

尽管二维Ising模型(基于哈密顿量调控的局部直接相互作用)与真实大脑动力学存在显著差异，但其仍能从统计物理视角为全脑动力学提供见解。需要注意的是，许多研究致力于改进Ising模型(如整合结构数据)，采用成对最大熵模型(PMEM)和最大似然估计等技术从数据中推断交互矩阵J和外磁场h等参数。但当感兴趣区域(ROI)数量较大时，使用PMEM就会存在一定的局限性。因此，本研究聚焦于基础的二维Ising模型，目的是创建一个基于机器学习的轻量级方法，可在333个ROI(大规模)的连接性水平上推断温度，从而实现大数据的可扩展性。深度学习模型能够整合多种模态数据，并且可以在不同的动力学系统上进行泛化。因此，本文提出了一种方法，该方法可以扩展至深度学习模型，能够从多种生物学信息驱动和改进的Ising模型中学习。

本研究采用几何深度学习方法提取网络连接结构中与Ising温度相关的复杂特征。这类专用于非欧几里得空间(如网络)的算法以图神经网络(GNN)为代表，通过模拟二维Ising模型网络训练，可预测产生特定动力学(及连接图)的控制参数(即Ising温度)。该模型经训练能准确预测Ising温度后，被应用于从静息态功能连接MRI(rs-fcMRI)图中估计神经发育期(8-22岁)的温度变化。

这一阶段以全脑(尤其是在大脑皮层)微环路的复杂重构为特征，这些改变会影响宏观脑动力学及其功能关系，反映了网络水平的脑动力学机制变化，代表了更灵活的动力学所支持的多样化神经模式，并对功能网络的拓扑属性产生了影响。由于二维Ising模型假设交互作用恒定且只考虑了最近邻耦合，基于功能连接网络的临界态转变差异只能通过Ising温度进行评估。

在此背景下，本研究利用Ising温度来探讨神经发育过程中临界状态如何通过控制参数的变化与功能网络的形成和演化相关。具体而言，随着年龄的增长，大脑逐渐脱离临界态，从无序阶段过渡到更有序的阶段(即神经活动的随机性降低、可预测性增强)。尽管受到拓扑表征缺失的限制，但本研究方法开创了一种轻量化的温度估算方式，能够处理大规模数据集(例如，适用于包含333个ROI的数据集)，并且可以扩展至从各种动态模型中提取的训练集，用以推断控制参数。

方法

数据描述

本研究使用Neuro Bureau提供的ADHD-200竞赛数据集。该公开数据集整合自8个不同研究中心，采用1.5T扫描仪采集静息态fMRI数据。数据集包含776个训练样本和197个测试样本，分为四组：健康对照组、ADHD混合型、ADHD多动冲动型和ADHD注意力缺陷型。为简化分析，本研究将所有ADHD亚型合并为一类，专注于ADHD与健康对照组之间的二元分类。除了基本的fMRI数据外，还提供了每个受试者的各种表型信息，包括年龄、性别、利手、智商和ADHD分型等。根据既往研究标准，设定头动阈值(平移＜3mm或旋转＜3°)，为严格控制头动影响，最终将阈值降至1.5mm和1.5°，最终纳入242名典型发育儿童和111名ADHD症状儿童。数据集使用Athena管道进行静息态fMRI预处理，并使用AFNI/FSL进行基于体素的形态学分析。BOLD信号时间序列数据来自Connectome网站(http://www.preprocessed-connectomes-project.org/adhd200/)。预处理步骤包括：剔除前4个时间点(确保磁化平衡)、时间层校正、头动校正、使用蒙特利尔神经研究所(MNI)空间进行4×4×4体素分辨率配准、带通滤波(0.009Hz＜f＜0.08Hz)以及6mm FWHM高斯平滑。各研究中心使用的扫描设备和参数详见ADHD-200-Webpage(2011)。对于有多次扫描的受试者仅取首次扫描数据，每位受试者统一截取140个时间点。

为控制年龄相关混杂因素，本研究评估了人口统计学信息：性别对年龄的影响(男性12.7±3.8岁，女性12.4±3.5岁)以及ADHD分型对年龄的影响(典型发育儿童12.8±3.7岁，ADHD混合型11.6±3.4岁，ADHD多动冲动型14.4±4.9岁，ADHD注意力缺陷型11.9±3.0岁)。此外，头动分析显示，最大旋转角度(度)与年龄之间呈显著负相关(Pearson's r=-0.13，p=0.01)，而最大平移(毫米)与年龄之间无显著相关性(r=0，p＞0.05)(见图1)。

图1.数据集的人口统计学特征。

为探究ADHD症状对头动的影响，本研究将数据分为典型发育儿童和ADHD症状组。仅ADHD组的最大平移(毫米)与年龄呈显著正相关(Pearson's r=0.27，p=0.001；Spearman's r=0.25，p=0.001)(见图2)。

图2.ADHD儿童(左)与典型发育儿童(右)的最大平移(毫米)和最大旋转(度)与年龄的相关性散点图(含Pearson相关系数和p值)。

图像处理

根据Gordon等人(2016)的研究，将功能数据划分为333个皮层感兴趣区(ROI)，该分区基于神经发育期皮层的剧烈变化及皮层网络的固有临界特性。计算各ROI内体素的平均BOLD信号后，通过Pearson相关系数构建邻接矩阵A(连接矩阵)。脑网络图定义为G=(V，E，A)，其中V={1，…，N}表示节点集合，E为边集合，A∈ℝN×N为加权邻接矩阵，矩阵元素wij表示节点i与j之间的连接强度(Pearson相关系数)。

二维Ising模型模拟

二维Ising模型通过二阶相变描述磁化现象。该模型由自旋网格构成，每个自旋为离散变量(取±1，即“上旋”或“下旋”)，其哈密顿量描述无外磁场时的自旋相互作用：

其中Σi,j表示近邻自旋求和，J为耦合常数。通过H=E可获得网格任意状态的能量信息。

鉴于Ising模型的离散随机特性，本研究采用Metropolis-Hastings蒙特卡洛算法进行模拟。该算法基于马尔可夫链的概念构建。所有模拟均在L=250的网格上进行，初始状态为自旋对齐或随机分布，蒙特卡洛时间步长t=L×L以确保每个自旋至少有一次状态翻转机会。当二维Ising模型通过临界温度Tc=2.269J时，会发生二阶相变(设玻尔兹曼常数k=1，耦合常数J=1)。

热平衡后，对250×250自旋网格进行200个时间步长的波动模拟(每个时间步长允许所有自旋改变其状态)。将网格按13×13自旋块进行平均以获得一个类fMRI信号的连续时间序列(图3)。由于250不能被13整除，因此采用取整法，最终得到19×19(361个ROI)的简化网格。计算全相关矩阵(361×361)后，从矩阵四边各移除14行/列，获得中心区域333×333的目标矩阵以消除边界效应(相关代码可在以下网址获取：https://github.com/Rodrigo-Motta/BRAIN_ISING_GNN)。

图3.二维Ising模型网络构建流程。

本研究在二维Ising模型中分别模拟了200和2000个时间点，发现图神经网络(GNN)在200个时间点的条件下对Ising温度的预测性能略优，因此选择200个时间点进行后续分析。需要说明的是，尽管存在性能差异，但2000个时间点的结果依然具有一致性。为验证无外磁场条件下二维Ising模型模拟流程的正确性与热平衡状态，本研究在[1.6，3.3]kT/J的范围内评估了序参量(单位格点平均磁化强度)。为确保模型的稳健性，针对三种不同晶格尺寸进行验证并与解析解对比，最终选择的晶格尺寸在临界区域表现出较小波动性，且其维度与Gordon的333个ROIs分区的维度接近。

图卷积神经网络

图卷积神经网络(GCN)的运算机制包含卷积操作、非线性激活函数、池化及反向传播四个核心环节。其中图卷积作为GCN的核心运算，通过消息传递机制提取图的拓扑结构。这种本征几何运算使GCN能够实现复杂的表征学习，并从网络中预测全局参数。本研究中，该模型通过图卷积层提取并转换节点的连接特征，进而学习估计图的全局属性。

虽然GNN适用于图结构数据处理，但仍需对输入图进行预处理。本研究采用k近邻算法确定图的连边关系(每个节点代表一个ROI，节点特征为其与其他ROIs的连接强度)，该方法有效避免了噪声拟合和梯度消失问题。具体网络架构包含三层图卷积、全局平均池化以及两个全连接层(见图4)，激活函数采用允许负相关值的Leaky-ReLU，并应用dropout技术和L2正则化的平均绝对误差(MAE)损失函数以防止过拟合。

图4.GCN回归器架构示意图。

建模数据集包含1500个均匀分布在临界点附近(区间=[1.8，2.5])的二维Ising模型模拟数据，其中1000个用于训练、250个用于验证、250个用于测试。采用循环优化技术进行100轮训练(学习率按三角波模式在设定极值间动态调整)，并通过验证集损失阈值筛选最优参数进行微调以提高精度。

基于GNN的Ising温度估计

GNN的输入本质上是图结构，这使其能有效分析具有复杂关联关系的数据。鉴于脑网络与Ising模型图具有相同维度，训练后的GNN可无缝处理这两种类型的数据。这种维度一致性确保了模型在Ising模型图上训练后，能有效估计脑网络的温度参数，实现从模拟数据到真实fMRI数据的迁移学习(方法论框架见图5)。经临界区([1.8，2.4]kBT/J)模拟数据训练及超参数优化后，GNN在测试集上展现出卓越的预测性能(MAE=0.04，R2=0.92)。

图5.方法论框架概览。

统计分析

首先采用Pearson相关和Spearman等级相关评估温度随年龄的整体变化趋势。通过偏差校正和加速bootstrap(1000次重采样)进行组间(ADHD与对照组)相关系数差异检验，每次重采样计算相关系数(平均排除36.8%的受试者)，设定95%置信区间。为了检验不同年龄段可能存在的温度差异，将数据划分为7个两年间隔的年龄组(7-9、9-11、11-13、13-15、15-17、17-19、19-21岁)，各组数据占比分别为：9-11岁(20.17%)、11-13岁(19.75%)、13-15岁(16.39%)、7-9岁(15.55%)、17-19岁(15.13%)、15-17岁(9.66%)、19-21岁(3.36%)。

结果

神经发育过程中Ising温度的估计

当GNN在200个时间点的Ising模型模拟数据上完成训练后，本研究分别估算了典型发育儿童和注意力缺陷多动障碍(ADHD)患者静息态功能磁共振(rs-fMRI)的Ising温度。通过计算估计温度T与临界温度Tc的差值(T-Tc)，Pearson和Spearman相关性分析均显示两组受试者的估计温度与年龄呈显著负相关：典型发育儿童(Pearson=-0.48，p＜0.0001；Spearman=-0.47，p＜0.0001)，ADHD患者(Pearson=-0.49，p＜0.0001；Spearman=-0.47，p＜0.0001)。使用2000个时间点训练的GNN同样显示出与此一致的显著负相关趋势。如图6所示，两组受试者的Ising温度均低于临界温度Tc，并随神经发育呈线性下降趋势。

图6.估计Ising温度与临界温度差值(T-Tc)随年龄变化的散点图。

将年龄划分为7个两年间隔的组别(7-9、9-11、11-13、13-15、15-17、17-19、19-21岁)后，计算各组温度与临界温度的绝对距离|T-Tc|。图7显示，无论是典型发育儿童还是ADHD患者，其与临界点的距离均随神经发育过程逐渐增大。虽然离临界点的中位距离持续增加，但并非所有年龄组均存在统计学差异。双样本t检验显示，年龄跨度较大的组别间存在显著差异(组间比较p值见表1)。

图7.不同年龄组的估计Ising温度与临界温度距离(|T-Tc|)的箱线图(对照组与ADHD组)。

表1.各年龄组间Ising温度与临界温度距离分布的配对t检验p值。

组间分析

通过自助法(bootstrap)比较两组相关系数分布的差异。该统计技术通过有放回地重复抽样构建相关系数分布，从而实现更稳健的组间比较。分析结果显示，ADHD症状对Ising温度随年龄下降的趋势无显著影响(Pearson相关p=0.5；Spearman相关p=0.6)，表明无论ADHD状态如何，大脑在成熟过程中远离临界态的趋势是一致的(图8)。

图8.对照组(蓝)与ADHD组(橙)相关系数的自助核密度估计(高斯核；带宽=1)。

混杂因素分析

考虑到头动对fMRI数据的显著影响，本研究通过线性回归评估了头动对Ising温度估计的潜在干扰。以最大平移(mm)、最大旋转(度)和年龄(岁)为自变量，估计Ising温度为因变量建立模型。对于典型发育儿童组(R2=0.19)，仅年龄与温度呈显著负相关(β=-0.43，双侧p＜0.001)；ADHD组(R2=0.32)也呈现相同的模式(β=-0.37，双侧p＜0.001)，结果如表2所示。

表2.头动对Ising温度估计影响的线性回归结果(线性系数β、p值和R2)。

讨论

本研究基于复杂系统、功能连接和临界态的既往研究，特别是受“大脑处于二阶相变临界点附近”证据的启发，采用Ising模型探究大脑状态。如Fraiman等(2009年)提出的，本框架假设二维Ising模型可描述实证静息功能网络，并将此概念拓展至神经发育领域。因此，本研究聚焦于大脑临界态而非功能网络拓扑结构的变化特征。临界模型是理解集体复杂神经动力学涌现的关键，对麻醉学、睡眠医学、发育行为儿科学和精神病学等多个临床领域具有重要价值。在临界相，模型实现了对齐相互作用与热激发(熵增)的平衡，这种平衡催生了多样化复杂模式的涌现。通过估算rs-fMRI的Ising温度，本研究建立了探究系统在不同条件(如神经发育)下可占据和转换的动态状态范围的桥梁。

研究主要发现，典型发育儿童和ADHD症状儿童的年龄与Ising温度均存在统计学显著负相关，这表明随着年龄的增长，大脑逐渐远离临界态并向有序状态转变。这一结果在二维Ising模型两种不同数量级的时间点条件下均保持一致。既往研究支持静息态脑动力学维持在略低于临界点的亚临界状态，但该结论仍存争议，另有证据表明大脑可能略微偏向顺磁态。不同预处理流程及多重混杂因素可能影响温度估计，因此最具价值的研究启示来自整体趋势。本研究表明大脑从无序相向更有序相转变，神经活动随机性降低而可预测性增强。两组相关系数在统计学上无显著差异，这表明以Ising温度衡量的成熟过程在典型发育儿童与ADHD儿童中遵循相似的轨迹(可能受ADHD患者用药的影响)。通过线性回归分析，本研究排除了头动因素对两组结果的潜在干扰。

现有推断Ising温度的主流框架多基于成对最大熵模型(PMEM)，该模型通过数据约束确定最优概率分布。这类技术还能估计耦合矩阵J(反映自旋间相互作用强度的差异)和外磁场h。另有研究采用扩散张量成像获取的脑结构几何特征来优化Ising模型的生物学合理性。尽管这些框架更为复杂稳健，但其拟合效果受ROI数量和时间点长度的限制，且计算大样本多ROI数据时需要极高的算力。本研究的机器学习框架通过简化Ising模型(耦合矩阵J恒定、仅考虑最近邻相互作用、忽略外磁场h、假设基础模型能描述实证数据)来规避上述问题，但也不可避免地存在局限性。该方法的显著优势在于：直接从连接信息估计控制参数，且在训练完成后能以极低计算成本处理多个ROI的连接矩阵数据，具有良好的可扩展性。

此外，本研究还为功能连接分析提供了一种新的方法：利用图结构模拟数据训练神经网络并应用于实证数据。由于脑网络与Ising网络共享图结构特征，该方法实现了理论模型与实证数据的整合。该框架可进一步扩展至Ising模型和其他临界模型(如描述振荡器/神经元相互作用产生集体态的Kuramoto模型)，使GNN能够学习复杂数据表征并从网络拓扑中提取参数信息。但必须验证训练所用模拟数据能否有效描述实证数据以确保推断的有效性。本研究的一个重要局限是缺乏功能连接网络的拓扑特征分析——尽管二维Ising模型通过临界温度有效捕捉了静息态功能网络，但当前方法可能忽略了影响大脑从无序态向有序态转变的关键网络重构特征。因此，未来的模型应致力于将拓扑指标整合到训练集模拟和事后分析中，以提升分析的完备性。

参考文献：Rodrigo M. Cabral-Carvalho, Walter H. L. Pinaya, João R. Sato; A graph neural network approach to investigate brain critical states over neurodevelopment. Network Neuroscience 2025; doi: https://doi.org/10.1162/netn_a_00451

小伙伴们关注茗创科技，将第一时间收到精彩内容推送哦～