动力学平均场模型
动力学平均场模型(Dynamic Mean Field model,简称 DMF 模型)是一类用群体平均变量描述大量相互作用单元随时间演化的模型。在全脑计算神经科学中,它将每个脑区表示为包含兴奋性和抑制性神经元的局部回路,再按照脑区之间的解剖连接耦合这些回路,从而研究局部突触过程如何产生全脑活动及其相关结构。本文主要介绍这一全脑神经动力学语境下的 DMF 模型[1]。
平均场近似把许多神经元的共同作用概括为少量群体变量。“动力学”强调这些变量保留时间演化:模型既能计算稳态,也能模拟刺激后的恢复、噪声引起的波动和脑区之间的相互影响。借助血流动力学观测模型,DMF 的神经活动还可转化为模拟的血氧水平依赖信号(blood-oxygen-level-dependent signal,BOLD 信号),与功能磁共振成像(fMRI)数据比较[1]。
基本信息
| 项目 | 内容 |
|---|---|
| 英文名称 | Dynamic Mean Field model |
| 常用缩写 | DMF |
| 研究尺度 | 神经元群体、脑区及全脑网络 |
| 典型状态变量 | 兴奋性与抑制性平均突触门控变量 |
| 数学形式 | 相互耦合的非线性常微分方程或随机微分方程 |
| 连接信息 | 扩散磁共振成像与纤维追踪估计的结构连接矩阵 |
| 主要输出 | 群体放电率、突触活动、模拟 BOLD 信号及功能连接 |
| 代表性研究 | Wong 与 Wang(2006);Deco 等(2013、2014);Mediano 等(2025) |
历史溯源
Wong 与 Wang 在 2006 年研究知觉决策时,将较复杂的脉冲神经网络约化为少量群体变量,保留慢速 NMDA 受体介导的循环兴奋。该工作说明,突触时间尺度可以为神经回路提供跨越数百毫秒的证据积累机制,也是后来神经 DMF 模型的重要方法来源[2]。
Deco 等在 2013 年将动力学平均场约化用于由真实解剖连接约束的全脑网络,并进一步分析稳态附近的均值与协方差。他们发现,在所研究的参数范围内,结构连接引导的噪声传播和接近失稳时的动力学减慢,可以产生类似静息态功能连接的相关模式[1]。2014 年的研究进一步明确保留局部兴奋与抑制群体,并引入反馈抑制控制(feedback inhibition control,FIC),用于调节不同脑区的局部抑制强度,使网络耦合增加时仍能维持较低的自发兴奋性放电率[3]。
DMF 随后被用于比较清醒、麻醉和脑损伤相关的全脑动力学,以及研究信息冗余、协同和跨脑区传递等现象[4][5]。
模型概述
从局部回路到全脑网络
设全脑被划分为 [math]\displaystyle{ N }[/math] 个区域。每个脑区包含兴奋性群体 E 和抑制性群体 I:兴奋性活动促进自身和抑制性群体的活动,抑制性活动则降低局部回路的净兴奋输入。脑区之间主要通过长程兴奋性连接相互作用。
结构连接矩阵记为 [math]\displaystyle{ \mathbf{C}=(C_{ij}) }[/math],其中 [math]\displaystyle{ C_{ij} }[/math] 表示脑区 [math]\displaystyle{ j }[/math] 对脑区 [math]\displaystyle{ i }[/math] 的连接权重。常见实现将 [math]\displaystyle{ C_{ii} }[/math] 设为零,因为脑区内部连接已由局部回路描述。连接组(connectome)指这种对脑区及其连接的整体表示。全局耦合参数 [math]\displaystyle{ G }[/math] 统一调节脑区间作用强度;其数值必须结合矩阵归一化方式理解[1][3]。
NMDA 和 AMPA 是介导兴奋性突触作用的受体,GABA 受体介导抑制性作用。NMDA 的慢动力学使兴奋作用能够持续积累;AMPA 的快速过程在常见约化中被视为近似即时响应。不同 DMF 版本保留的突触过程和状态变量数目可以不同;下文采用 Deco 等 2014 年的双群体形式[3]。
输入电流与群体放电率
兴奋性与抑制性群体的输入电流可写为:
-
[math]\displaystyle{ \begin{aligned} I_i^E={}&W_E I_0+w_+J_{\mathrm{NMDA}}S_i^E\\ &+GJ_{\mathrm{NMDA}}\sum_{j=1}^{N}C_{ij}S_j^E-J_iS_i^I+u_i(t). \end{aligned} }[/math]
(1)
-
[math]\displaystyle{ \begin{aligned} I_i^I={}&W_I I_0+J_{\mathrm{NMDA}}S_i^E-J_{II}S_i^I\\ &+\lambda GJ_{\mathrm{NMDA}}\sum_{j=1}^{N}C_{ij}S_j^E. \end{aligned} }[/math]
(2)
其中,[math]\displaystyle{ S_i^E }[/math] 和 [math]\displaystyle{ S_i^I }[/math] 为平均突触门控变量,即描述群体突触活动强度的状态变量;[math]\displaystyle{ I_0 }[/math] 为背景输入;[math]\displaystyle{ W_E,W_I }[/math] 为背景输入缩放系数;[math]\displaystyle{ w_+ }[/math] 为局部循环兴奋权重;[math]\displaystyle{ J_{\mathrm{NMDA}} }[/math] 为兴奋耦合系数;[math]\displaystyle{ J_i }[/math] 为抑制群体作用于兴奋群体的反馈强度;[math]\displaystyle{ u_i(t) }[/math] 为外部刺激电流。式(2)中 [math]\displaystyle{ J_{II}=1\,\mathrm{nA} }[/math],对应原文以 nA 表示电流时抑制自反馈项的单位系数。[math]\displaystyle{ \lambda=0 }[/math] 表示长程输入仅进入兴奋群体,[math]\displaystyle{ \lambda=1 }[/math] 则同时考虑长程输入对抑制群体的前馈作用[3]。
电流通过输入—输出函数转化为平均放电率:
-
[math]\displaystyle{ \begin{aligned} r_i^p&=H_p(I_i^p),\qquad p\in\{E,I\},\\ H_p(u)&=\frac{a_pu-b_p}{1-\exp[-d_p(a_pu-b_p)]}. \end{aligned} }[/math]
(3)
参数 [math]\displaystyle{ a_p,b_p,d_p }[/math] 分别控制输入敏感度、阈值位置和阈值附近的平滑程度。[math]\displaystyle{ H_p }[/math] 概括带噪声神经元的电流—放电关系:低输入时放电率较低,超过阈值后逐渐升高。分子和分母同时为零处是可去奇点,其极限为 [math]\displaystyle{ 1/d_p }[/math],数值计算时需使用该极限或稳定的等价表达式[3]。
突触动力学与单位约定
为明确单位,以下以毫秒(ms)为时间单位、赫兹(Hz)为放电率单位。兴奋性和抑制性门控变量满足:
-
[math]\displaystyle{ \frac{dS_i^E}{dt}=-\frac{S_i^E}{\tau_E} +\frac{\gamma}{1000}(1-S_i^E)r_i^E+\sigma_E\xi_i^E(t). }[/math]
(4)
-
[math]\displaystyle{ \frac{dS_i^I}{dt}=-\frac{S_i^I}{\tau_I} +\frac{r_i^I}{1000}+\sigma_I\xi_i^I(t). }[/math]
(5)
这里 [math]\displaystyle{ \tau_E }[/math] 和 [math]\displaystyle{ \tau_I }[/math] 为突触时间常数,[math]\displaystyle{ \gamma }[/math] 为兴奋门控动力学系数;[math]\displaystyle{ \xi_i^p(t) }[/math] 表示零均值、彼此独立的高斯白噪声,[math]\displaystyle{ \sigma_p }[/math] 控制群体波动幅度。白噪声的数值幅度应随时间单位及离散化约定保持一致。式(4)、(5)中的 [math]\displaystyle{ 1/1000 }[/math] 将“每秒”的放电率转换为“每毫秒”。如果全式改用秒,时间常数和噪声约定也必须相应转换。
常见参数为 [math]\displaystyle{ \tau_E=100\,\mathrm{ms} }[/math]、[math]\displaystyle{ \tau_I=10\,\mathrm{ms} }[/math] 和 [math]\displaystyle{ \gamma=0.641 }[/math]。兴奋门控中的 [math]\displaystyle{ 1-S_i^E }[/math] 描述开放过程的饱和效应;抑制变量采用线性化的门控描述。该约化模型直接追踪突触变量,放电率由式(3)得到[3]。
动力学机制
兴奋—抑制反馈与全局耦合
局部循环兴奋增强时,兴奋群体的输入和放电率随之增大;它也带动抑制群体活动,后者反过来限制兴奋活动。这种相互作用与兴奋-抑制平衡密切相关。增加 [math]\displaystyle{ G }[/math] 则增强跨脑区输入,可能使低放电稳态接近失稳,也可能在一定参数范围内形成多个稳定活动状态,即多稳态[1][3]。
FIC 是一种参数校准过程:针对给定的连接矩阵和 [math]\displaystyle{ G }[/math],逐脑区调整 [math]\displaystyle{ J_i }[/math],将兴奋性平均放电率维持在约 3 Hz 的目标附近。它让研究者能够比较低放电背景下不同网络耦合的作用。是否使用 FIC、以及改变 [math]\displaystyle{ G }[/math] 时是否重新校准 [math]\displaystyle{ J_i }[/math],都会影响模型行为,应明确记录[3]。
稳态、线性波动与相关结构
将所有门控变量组成状态向量:
-
[math]\displaystyle{ \begin{aligned} \mathbf{s}&=(S_1^E,\ldots,S_N^E,S_1^I,\ldots,S_N^I)^{\mathsf T},\\ \frac{d\mathbf{s}}{dt}&=\mathbf{f}(\mathbf{s})\qquad\text{(无噪声时)}. \end{aligned} }[/math]
(6)
稳态满足 [math]\displaystyle{ \mathbf{f}(\mathbf{s}_*)=\mathbf{0} }[/math]。在其附近线性化,可得到:
-
[math]\displaystyle{ \begin{aligned} \delta\dot{\mathbf{s}}&=\mathbf{A}\delta\mathbf{s}+\boldsymbol{\eta}(t),\\ \mathbf{A}&=\left.\frac{\partial\mathbf{f}}{\partial\mathbf{s}}\right|_{\mathbf{s}_*}. \end{aligned} }[/math]
(7)
[math]\displaystyle{ \mathbf{A} }[/math] 是雅可比矩阵,[math]\displaystyle{ \boldsymbol{\eta} }[/math] 为线性近似中的噪声向量。如果 [math]\displaystyle{ \mathbf{A} }[/math] 的全部特征值实部均为负,稳态局部稳定;当某个实部趋近零时,相应方向的扰动衰减变慢,可能放大低频波动。若稳定条件成立,门控波动的平稳协方差矩阵 [math]\displaystyle{ \mathbf{P} }[/math] 满足连续李雅普诺夫方程:
-
[math]\displaystyle{ \mathbf{A}\mathbf{P}+\mathbf{P}\mathbf{A}^{\mathsf T}+\mathbf{Q}=\mathbf{0}. }[/math]
(8)
其中 [math]\displaystyle{ \mathbf{Q} }[/math] 为噪声强度矩阵。式(8)说明,群体活动之间的相关结构同时依赖局部响应、脑区间连接和噪声。这是将解剖结构与功能相关联系起来的一种近似分析方法;在远离所选稳态或波动很强时,需要直接研究非线性随机动力学[1][3]。
从神经活动到 BOLD 信号
fMRI 的 BOLD 信号反映神经活动引起的血流、血容量与血氧变化。DMF 常与 Balloon–Windkessel 血流动力学模型结合,将区域神经活动作为观测模型的输入。不同实现可能采用兴奋性放电率或兴奋性突触活动作为驱动量,复现时应遵循所用版本的定义[3][4]。
设模拟 BOLD 时间序列为 [math]\displaystyle{ y_i(t) }[/math],常用的静态功能连接为:
-
[math]\displaystyle{ \mathrm{FC}_{ij}=\mathrm{corr}\bigl(y_i(t),y_j(t)\bigr). }[/math]
(9)
这里“功能连接”指活动的统计相关,“结构连接”指解剖通路的估计。两者之间由神经动力学、网络传播和血流动力学共同联系。还可计算滑动时间窗内的功能连接,以研究功能连接动态(functional connectivity dynamics,FCD),并用经验 FCD 与模拟 FCD 的相似程度校准参数[4]。
数值示例:局部回路对短时输入的响应
图 1 使用式(1)—(5)的确定性单脑区模型,[math]\displaystyle{ G=0 }[/math]、[math]\displaystyle{ \lambda=0 }[/math]、[math]\displaystyle{ J_i=J_{II}=1\,\mathrm{nA} }[/math],其余参数如下[3]:
| 参数 | 数值 | 含义 |
|---|---|---|
| [math]\displaystyle{ I_0 }[/math] | 0.382 nA | 背景电流 |
| [math]\displaystyle{ W_E,W_I }[/math] | 1,0.7 | 背景电流缩放 |
| [math]\displaystyle{ w_+ }[/math] | 1.4 | 局部循环兴奋权重 |
| [math]\displaystyle{ J_{\mathrm{NMDA}} }[/math] | 0.15 nA | 兴奋耦合系数 |
| [math]\displaystyle{ a_E,a_I }[/math] | 310,615 nC−1 | 输入敏感度 |
| [math]\displaystyle{ b_E,b_I }[/math] | 125,177 Hz | 输入—输出函数参数 |
| [math]\displaystyle{ d_E,d_I }[/math] | 0.16,0.087 s | 阈值附近形状参数 |
| [math]\displaystyle{ \tau_E,\tau_I }[/math] | 100,10 ms | 突触时间常数 |
| [math]\displaystyle{ \gamma }[/math] | 0.641 | 兴奋门控动力学系数 |
先求无输入稳态,再用 DOP853 自适应积分模拟 1200 ms,每 0.1 ms 保存一次结果。将最大积分步长从 0.5 ms 减半为 0.25 ms 后,两次计算的放电率最大差异小于 [math]\displaystyle{ 10^{-6}\,\mathrm{Hz} }[/math]。
本示例的初始兴奋性和抑制性放电率约为 3.08 Hz、3.92 Hz。刺激期间,兴奋活动持续积累,放电率最高约 9.62 Hz,抑制性放电率随后升高至约 5.87 Hz。撤去刺激时,直接输入立即减小,突触状态则连续变化,随后逐渐恢复。[math]\displaystyle{ \tau_E=100\,\mathrm{ms} }[/math] 表示兴奋突触的固有时间尺度;整个回路的恢复速度还受到循环兴奋和抑制反馈的共同影响。
研究案例:PNAS 的全脑信息整合实验
论文与实验设计
Mediano 等在 2025 年发表于《美国国家科学院院刊》(PNAS)的论文《Toward a unified taxonomy of information dynamics via Integrated Information Decomposition》中,用全脑 DMF 检验不同信息整合指标的解释。这项工作提出整合信息分解(Integrated Information Decomposition,ΦID),将系统过去对未来的预测信息区分为冗余、独有、协同及它们之间的转换[5]。
其图 6 的实验将 Human Connectome Project 的扩散成像连接组与全脑分区结合;示意图标注为 Lausanne-83。每个脑区采用局部兴奋—抑制平均场动力学,区域神经活动经血流动力学模型转换为 BOLD 信号。研究者扫描全局耦合 [math]\displaystyle{ G }[/math],计算所有脑区对的 [math]\displaystyle{ \Phi^{\mathrm{WMS}} }[/math] 与 [math]\displaystyle{ \Phi^{\mathrm R} }[/math],再分析脑区对之间的平均值。论文采用与 Herzog 等(2020)相同的模型设置,完整模拟和估计细节列于其补充材料[5]。
两种信息指标为何不同
对任意一对脑区,设 [math]\displaystyle{ \mathbf{x}_t=(x_t^1,x_t^2)^{\mathsf T} }[/math] 表示两脑区在时刻 [math]\displaystyle{ t }[/math] 的活动,[math]\displaystyle{ \Delta\gt 0 }[/math] 为所选时间滞后。[math]\displaystyle{ \mathcal{I}(\cdot;\cdot) }[/math] 表示互信息,即两个变量之间共享的统计信息。常用的“整体减去各部分”指标(whole-minus-sum,WMS)为:
-
[math]\displaystyle{ \begin{aligned} \Phi^{\mathrm{WMS}}={}&\mathcal{I}(\mathbf{x}_t;\mathbf{x}_{t+\Delta})\\ &-\sum_{k=1}^{2}\mathcal{I}(x_t^k;x_{t+\Delta}^k). \end{aligned} }[/math]
(10)
第一项衡量两个脑区的联合过去对联合未来的预测信息,第二项减去各脑区自身的预测信息。当两个脑区都携带同一份信息时,这份冗余信息可能在第二项中被重复减去,使结果变成负值。ΦID 将这一现象定位为同时在过去和未来两端冗余的成分,即“双重冗余”[5]。
在论文采用的最小互信息冗余定义(minimum mutual information,MMI)下,修正指标为:
-
[math]\displaystyle{ \Phi^{\mathrm R}=\Phi^{\mathrm{WMS}} +\min_{k,\ell\in\{1,2\}}\mathcal{I}(x_t^k;x_{t+\Delta}^{\ell}). }[/math]
(11)
式(11)补回被重复扣除的双重冗余项。论文将修正后的整合理解为协同与独有信息传递的贡献:协同指联合观察两脑区才能获得的预测信息,独有信息传递指一个脑区独有的过去信息进入另一个脑区的未来。式(11)对应 MMI 定义;采用其他冗余定义时,应使用相应的双重冗余项[5]。
实验结果与含义
图 2 中,平均放电率在 [math]\displaystyle{ G\approx2 }[/math] 附近陡升。同一区域内,[math]\displaystyle{ \Phi^{\mathrm{WMS}} }[/math] 明显下降至负值,而 [math]\displaystyle{ \Phi^{\mathrm R} }[/math] 出现强正峰值,并在所扫描范围内保持非负。作者指出,这一区域与相关先前研究中较好拟合清醒受试者数据的耦合范围相符[5]。
该实验展示了 DMF 的一种用途:在连接结构和局部机制明确的模拟系统中,比较分析指标如何响应同一次动力学变化。[math]\displaystyle{ \Phi^{\mathrm{WMS}}\lt 0 }[/math] 在此反映指标对冗余信息的重复扣除,不能直接解释为大脑具有“负的信息整合”。修正指标则揭示耦合变化附近增强的协同与传递成分[5]。
读图时还需区分三个层次:图中的信息曲线是脑区对指标的平均值;[math]\displaystyle{ G\approx2 }[/math] 是该连接组、参数和归一化方式下的结果;该模拟用于研究信息动力学与指标解释,其结果本身不构成个体意识状态的诊断标准。更换连接矩阵、局部参数、抑制控制或信号估计流程,峰值位置和数值都可能改变。
神经科学意义与典型应用
结构与功能的联系
DMF 可以在固定连接组上改变局部时间常数、兴奋增益、抑制强度和噪声,检验哪些机制能产生经验功能连接。它还能通过改变网络边或对局部脑区施加输入,研究扰动如何传播。结构提供传播路径,局部动力学决定响应与稳定性,观测模型决定这些活动如何体现在 fMRI 信号中[1][3]。
神经调质、麻醉与脑损伤
区域受体分布可以用于调节局部群体增益,使模型纳入脑区之间的生物学差异。Luppi 等(2022)结合 GABA-A 受体分布与 DMF 研究丙泊酚麻醉,并通过替换或扰动结构连接研究脑损伤相关的意识障碍。其模拟将受体作用变化和解剖连接损伤作为不同干预,用以解释它们如何产生相近的宏观活动特征[4]。这类工作为机制假设提供可计算的检验方式,具体参数的生物学解释仍需结合实验和模型验证。
参数拟合与复现要点
全脑研究通常先确定脑区分区与结构矩阵,再校准耦合和局部参数,生成与经验数据时长及采样方式可比的信号。报告结果时应说明以下信息:
- 连接定义:分区方案、矩阵方向、对角线处理、权重归一化,以及是否包含传导时延。
- 局部模型:保留哪些群体和突触变量,是否使用 FIC,改变耦合时是否重新校准抑制。
- 数值计算:时间单位、积分方法、步长、初始条件、丢弃暂态的时长,以及噪声定义和随机种子。
- 观测与比较:BOLD 驱动量、血流动力学参数、滤波与采样流程,以及功能连接或 FCD 的拟合目标。
- 信息分析:用于计算的活动变量、时间滞后、互信息估计方法、冗余定义和脑区对汇总方式。
对于 PNAS 案例,完整复现应使用原文的连接数据和补充材料所列流程;图 1 的单脑区刺激示例用于说明基本方程,两项实验的条件与目的分别如各自图注所述。
优势、局限与其他模型的关系
DMF 的优势在于,它以少量群体状态保留突触时间尺度、兴奋—抑制回路和解剖连接,计算成本通常低于逐神经元、逐突触的全脑脉冲模拟。它也提供稳态与线性波动分析的途径,使网络结构与统计相关之间的关系更易研究[1]。
平均场约化会压缩群体内部的异质性和精细脉冲时序;连接成像及脑区分区也带有测量与尺度选择的影响。拟合良好的宏观活动可能由多组参数产生,因此需要多个观测量、参数敏感性分析或外部验证来约束机制解释。FIC、时延、噪声和血流动力学设置均可能影响结果,应随模型版本一并说明。
| 模型 | 典型变量与尺度 | 主要用途及与 DMF 的联系 |
|---|---|---|
| 霍奇金-赫胥黎模型 | 单神经元膜电位与离子通道门控 | 描述动作电位的膜机制;DMF 将大量神经元的作用约化到群体层面 |
| FitzHugh-Nagumo模型 | 快变量与慢恢复变量 | 研究兴奋性与振荡的简化动力学 |
| 鲁尔科夫映射 | 离散时间的快慢变量 | 以较低计算成本模拟放电、簇发与网络活动 |
| Wilson-Cowan模型 | 兴奋与抑制群体平均活动 | 研究群体回路及振荡;常见神经 DMF 更明确联系突触门控、脉冲网络约化和全脑观测模型 |
| 动力学平均场模型 | 突触门控、派生放电率与脑区网络 | 研究解剖连接约束下的全脑活动、功能相关与信息动力学 |
这些模型可在不同研究尺度上互相补充。模型的选择取决于需要解释的是单神经元放电、局部群体响应,还是跨脑区的宏观活动。
参考文献
- ↑ 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Gustavo Deco; Adrián Ponce-Alvarez; Dante Mantini; Gian Luca Romani; Patric Hagmann; Maurizio Corbetta (2013). "Resting-State Functional Connectivity Emerges from Structurally and Dynamically Shaped Slow Linear Fluctuations". The Journal of Neuroscience. 33 (27): 11239–11252. doi:10.1523/JNEUROSCI.1091-13.2013.
- ↑ Kong-Fatt Wong; Xiao-Jing Wang (2006). "A Recurrent Network Mechanism of Time Integration in Perceptual Decisions" (PDF). The Journal of Neuroscience. 26 (4): 1314–1328. doi:10.1523/JNEUROSCI.3733-05.2006.
- ↑ 3.00 3.01 3.02 3.03 3.04 3.05 3.06 3.07 3.08 3.09 3.10 3.11 Gustavo Deco; Adrián Ponce-Alvarez; Patric Hagmann; Gian Luca Romani; Dante Mantini; Maurizio Corbetta (2014). "How Local Excitation–Inhibition Ratio Impacts the Whole Brain Dynamics". The Journal of Neuroscience. 34 (23): 7886–7898. doi:10.1523/JNEUROSCI.5068-13.2014.
- ↑ 4.0 4.1 4.2 4.3 Andrea I. Luppi; Pedro A. M. Mediano; Fernando E. Rosas; et al. (2022). "Whole-brain modelling identifies distinct but convergent paths to unconsciousness in anaesthesia and disorders of consciousness". Communications Biology. 5: 384. doi:10.1038/s42003-022-03330-y.
- ↑ 5.0 5.1 5.2 5.3 5.4 5.5 5.6 Pedro A. M. Mediano; Fernando E. Rosas; Andrea I. Luppi; Robin L. Carhart-Harris; Daniel Bor; Anil K. Seth; Adam B. Barrett (2025). "Toward a unified taxonomy of information dynamics via Integrated Information Decomposition". Proceedings of the National Academy of Sciences. 122 (39): e2423297122. doi:10.1073/pnas.2423297122.
署名
贡献者:杨明哲。