方法文章

利用人类新皮层神经求解器揭示神经治疗对脑电图影响的神经机制研究方案

389 次观看

DOI:

10.3791/70618

2026年5月19日

本文内容

摘要

本方案展示了如何利用基于物理的神经模拟来解读神经治疗药物的电生理生物标志物,并揭示其对神经回路的影响,从而为神经治疗药物的开发提供一种具有机制基础的方法。

摘要

脑电图(EEG)和电生理学方法可提供毫秒级分辨率的中枢神经系统疾病生物标志物,广泛用于评估治疗相关效应。然而,对产生这些生物标志物的神经机制理解有限,阻碍了基于这些信号的诊断和治疗方法的发展。人类新皮层神经求解器(Human Neocortical Neurosolver, HNN)是一款开源的生物物理建模软件,可将局部化的脑电图生物标志物与其多尺度神经发生源相联系。本实验方案展示了一种基于假设的工作流程,利用 HNN 通过优化模型参数,使模拟的与实测的电流源波形达到匹配,从而检验神经治疗诱导的脑电图生物标志物的神经机制。随后可对相应的多尺度细胞和环路水平活动进行可视化和量化分析,为后续实证研究中模型预测的验证提供目标。本文提供了一个示例,演示如何研究听觉诱发电位反应中早期事件相关电位成分(P1、N1 和 P2)背后的神经机制,并评估神经环路活动在神经治疗干预后的变化。本方案能够设计模拟实验,生成可检验的预测,将脑电图生物标志物与潜在的神经环路机制联系起来。类似的工作流程也可用于研究疾病机制或其他治疗干预手段。

引言

中枢神经系统(CNS)治疗药物的研发面临独特的挑战,其获批率低于其他疾病领域,凸显了对创新性方法学的迫切需求1,尤其是那些能够揭示治疗对脑动态影响的方法。研究治疗手段对神经活动影响的一种成熟方法是脑电图(EEG)2,3。EEG可提供in vivo条件下脑环路水平动态活动的特征性信号,并具有从啮齿类动物模型向人类临床试验转化的强大潜力,因为生成EEG信号的神经环路在不同物种间具有同源性4,5,6,7,8。在药物研发中,EEG可发挥多种作用,包括提供动物与人类研究之间的转化读数、评估药物安全性、指导化合物筛选、阐明剂量-效应关系、在早期临床阶段验证作用机制,以及实现临床试验的分层和队列富集9,10,11,12,13,14。尽管具有上述优势,EEG信号的解读仍是一项重大挑战,尤其是在试图将观察到的变化与潜在神经机制相联系时。

在中枢神经系统药物研发中,一种可靠的脑电图生物标志物是事件相关电位(ERP)。ERP反映与感觉刺激同步的脑部活动,已被广泛用于研究神经发育和神经精神障碍,包括抑郁症15,16、精神分裂症17,18、自闭症谱系障碍19,20以及阿尔茨海默病21。ERP还被用于评估治疗对脑环路的影响及剂量范围22,23,24,25,其中向健康反应的正常化可能提示治疗有效性26。然而,ERP及其他脑电图生物标志物(如脑节律振荡)的一个关键局限在于,它们与疾病状态或药物效应之间的关联在很大程度上仅为相关性。尽管统计分析能够识别生物标志物与结局之间的关系,但无法揭示特定神经环路组分产生这些信号的具体机制。因此,特定细胞类型和环路机制的因果性贡献仍不明确。阐明脑电图信号的细胞和环路起源,可将观察到的信号特征与潜在生理机制相联系,从而显著提升其应用价值27,28。在本文中,“生物标志物”一词指治疗干预后脑电图信号的可测量变化,符合美国食品药品监督管理局—美国国立卫生研究院生物标志物、终点及其他工具(FDA–NIH BEST)框架的定义29,而非意味着已正式认证可用于特定临床用途30

尽管侵入性电生理记录可提供细胞和神经环路水平的详细信息,但这些方法主要局限于动物模型,难以直接应用于人类研究。其他替代方法(如逆向建模技术)可根据脑电图(EEG)信号估计源活动,但往往缺乏对底层神经环路的明确机制性表征。生物物理模拟通过建模神经环路产生可测量EEG信号的物理过程,提供了一种互补的研究框架31,32,33,34图1)。与纯粹的统计学生物标志物分析或缺乏机制基础的逆向方法相比,生物物理建模能够直接检验神经环路动力学与观测到的电生理信号之间关联的假设。

figure-introduction-1
图1.利用生物物理建模开发和验证药理学脑电图(EEG)生物标志物的机制假说。A)基于不同条件下脑信号差异识别EEG生物标志物。例如,听觉事件相关电位(ERP)在治疗后条件(红色)下较治疗前条件(蓝色)减弱。(B)生物物理建模可用于检验解释EEG生物标志物如何产生以及如何随药物干预而变化的机制假说。研究者提出关于药物诱导神经活动变化的假说,并确定相应的模型参数。(C)默认的人类新皮层神经求解器(HNN)模型可作为检验假说的起点,通过手动调整模型参数或应用自动优化与推断算法进行分析。治疗前与治疗后条件之间参数值的差异对应于基于模型的预测结果。请点击此处查看该图的放大版本。

本方案采用人脑新皮层神经求解器(Human Neocortical Neurosolver, HNN),这是一种开源的生物物理建模框架,用于将治疗相关效应的事件相关电位(ERP)生物标志物与其潜在的细胞和环路层面机制联系起来33图2)。HNN基于这样一个原理:排列整齐的锥体神经元树突中的同步胞内电流流动产生脑电图(EEG)信号背后的主电流偶极子6,35,36,37。该模型代表了一个典型的新皮层柱,包含分布于皮层第2/3层和第5层的兴奋性锥体神经元和抑制性中间神经元。默认的HNN网络每层包含100个锥体神经元和33个抑制性神经元,构成一个简化但具有生物学基础的皮层环路表征。锥体神经元采用多区室树突结构建模,以捕捉关键的形态学特征38,而由于抑制性神经元对细胞外电流的贡献有限,因此被表示为单区室模型33。突触相互作用包括兴奋性的α-氨基-3-羟基-5-甲基-4-异噁唑丙酸(AMPA)受体和N-甲基-D-天冬氨酸(NMDA)受体,以及抑制性的γ-氨基丁酸A型和γ-氨基丁酸B型(GABAB)受体,所有神经元均包含由霍奇金-赫胥黎(Hodgkin–Huxley)动力学调控的主动离子电导。

figure-introduction-2
图 2.HNN 模型示意图。 HNN 模型主要组成部分的可视化图示,包括兴奋性神经元与抑制性神经元之间的局部网络连接,以及称为“近端输入”和“远端输入”的外源性输入通路。请点击此处查看该图的放大版本。

HNN 中的神经活动由代表前馈和反馈通路的外源性输入驱动。前馈“近端”输入对应于来自髓板内丘脑并作用于近端树突的输入,而反馈“远端”输入则代表作用于远端树突的皮层-皮层和非髓板内丘脑输入。这些输入被建模为动作电位序列,可诱发突触电流,并在锥体神经元树突上产生细胞内电流流动。由此产生的群体水平电流偶极子以纳安·米(nanoampere-meters)表示,可直接与方向受限的源定位脑电图(EEG)或脑磁图(MEG)数据进行比较。HNN 的默认参数化基于体感皮层研究的实证数据39,40,41,并已成功应用于听觉42,43,44、视觉45和额叶皮层信号46,其模型推导出的预测结果已在后续实验研究中得到验证7,41,47

HNN 仿真可应用于药物研发的多个阶段,包括靶点验证、药物作用机制比较、剂量优化以及后续实验的假设生成14,48,49,50。这使得研究人员能够将机制性建模整合到实际科研工作流程中,支持关于神经治疗药物如何影响神经环路的假设生成与验证。在本方案中,我们重点关注听觉事件相关电位(ERPs)中的早期 P1、N1 和 P2 成分,因为这些特征已被充分表征,可为假设驱动的建模提供约束条件51。尽管本研究重点在于药物诱导的变化,但该方法也可拓展至其他神经治疗干预手段,如脑刺激或行为训练,以及中枢神经系统疾病的研究。

HNN 的使用遵循一种迭代建模框架,其中模型结构和参数最初由现有数据限定,随后通过与实证观察结果的比较进行优化。大规模神经网络模型包含大量参数,但仅有一部分——称为感兴趣参数——会被调整以检验特定假设。这些参数并非随意选择,而是基于先前的实验证据和文献描述的神经治疗机制潜在作用方式来确定。在本方案中,以外源性输入的时间和强度、局部抑制性连接以及树突离子通道电导为例,说明可能受神经治疗手段影响的、具有生物学可解释性的变量。

从默认模型开始,用户首先通过手动调节与自动优化相结合的方法,将参数拟合到治疗前的事件相关电位(ERP)数据。手动调节可调整全局缩放和输入参数,以近似实际观测的波形,从而直观理解参数变化对模型输出的影响。随后采用自动优化方法,如协方差矩阵自适应进化策略(CMA-ES)、贝叶斯优化和基于线性近似的约束优化,进一步精调参数值并提升拟合效果。在建立治疗前模型后,针对假设能够解释治疗后变化的参数进行调整,以拟合治疗后的ERP数据。

为解决参数估计中的不确定性问题,采用基于模拟的推断(simulation-based inference, SBI)来估计能够重现观测数据的参数值分布52,53。SBI 考虑了不同参数组合可能产生相似输出的可能性,并可实现对参数不确定性的量化。可通过重叠指数(overlap index, OVL)54,55评估治疗前与治疗后参数分布之间的差异,从而揭示潜在的作用机制。

该方法的一个关键优势在于,将模型拟合到特定数据模态时,能够生成跨越多个神经活动尺度的预测结果,包括单细胞放电、特定脑层的局部场电位(LFPs)以及电流源密度(CSD)。这些预测为使用互补技术进行实验验证提供了目标。如果预测结果未得到实证数据支持,则可通过引入新的约束条件来更新模型,从而形成一个假设生成、检验与优化的迭代循环(图3)。

figure-introduction-3
图3. 基于HNN开发与测试ERP生物标志物预测的迭代工作流程 工作流程对应于实验方案的步骤。脑电生物标志物的识别及默认HNN模型的初始化以红色表示(步骤1–2)。通过手动调节与优化,将模型参数拟合至治疗前和治疗后的ERP信号(紫色;步骤3–5)。基于模拟推断(SBI)的不确定性量化以绿色表示(步骤6)。随后,检验模型预测结果并与实验数据进行比较,以验证模型或进一步约束模型参数(橙色;步骤7)。 请点击此处查看此图的放大版本。

本方案专为在诱发反应范式期间采集的具有方向约束且经过源定位的脑电图(EEG)或脑磁图(MEG)数据而设计。可采用标准预处理和源定位方法(例如最小范数估计 [MNE]-Python56)生成所需的输入数据。以纳安·米(nanoampere-meters)表示的源水平信号可直接与HNN的输出结果进行比较。对于快速感觉反应,源水平和传感器水平的信号通常高度相似,因此源定位建模所得的见解可用于指导对传感器水平EEG数据的解释57,58

方案

所有涉及人类数据的操作均按照相关机构的指南和规定进行。本研究使用的数据集来自先前已发表的研究43,无需额外的伦理审批。本方案不涉及任何危险材料或操作。

1. 确定治疗诱导的脑电图事件相关电位生物标志物并定义模型假设

  1. 收集或确定一个包含来自目标受试者(例如,在神经治疗背景下,治疗前和治疗后)实验记录的脑电图(EEG)信号的数据集。在呈现感觉刺激期间记录EEG信号,并同时记录感觉刺激的时间戳,以便将数据分割为多个试验。确保EEG数据以与预处理软件兼容的格式存储(例如,.fif、.set 或 .edf)。
    注意: 相关代码仓库(https://github.com/ntolley/hnn_jove)提供了用于生成代表性结果的数据文件。该仓库包含来自Kohl等(2022)的预处理听觉MEG事件相关电位(ERP),用作治疗前ERP(原始数据可从以下网址获取:https://github.com/kohl-carmen/HNN-AEF)。假设的治疗后ERP是通过对治疗前波形应用高斯 taper 窗函数进行缩放生成的。相应的数据文件位于仓库中的 data/pre-treatment.txt 和 data/post-treatment.txt。由于MEG和EEG信号反映了相似的神经源基础,本方案适用于这两种模态。
  2. 确定一组假设可区分治疗相关效应的候选ERP生物标志物特征(例如,ERP峰值时间点和幅值)。
    注意:在本示例方案中,使用峰值幅值作为感兴趣的生物标志物。
  3. 对EEG数据进行预处理,并提取感兴趣的生物标志物特征。
    注意:多个软件包支持EEG数据的预处理和ERP分析,包括 MNE-Python56、EEGLAB59 和 FieldTrip60。建议进行源定位以建模ERP信号,但非必需。示例工作流程见 https://jonescompneurolab.github.io/hnn-core/stable/auto_examples/workflows/plot_simulate_somato.html。已有若干先前研究详细描述了EEG信号的预处理与分析过程;读者可特别参考56,61 以获得更完整的背景知识。
    1. 使用所有通道的传感器水平信号进行源定位,或选择待分析的EEG传感器。使用源定位数据可与模型输出进行直接比较;传感器水平数据不具备单位对应关系。
      注意:下文所述的一一对应单位关系不适用于传感器水平信号。
    2. 利用感觉刺激的时间戳将记录的EEG数据分割为多个试验。
    3. 计算治疗前和治疗后条件下的试验平均ERP波形。
    4. 从试验平均波形中提取候选ERP生物标志物(例如,计算N1峰值幅值)。在提取前应预先定义峰值检测标准(例如,时间窗口和极性)。
  4. 进行统计检验,以确定哪些ERP特征在不同条件下存在显著差异(例如,治疗前与治疗后)。根据研究设计选择适当的统计检验方法,并在必要时进行多重比较校正(例如,重复测量方差分析后接Tukey HSD事后检验用于多重比较)。
    注意:统计检验的代码示例见 https://mne.tools/stable/auto_tutorials/stats-sensor-space/20_erp_stats.html。
  5. 输出特定的、具有统计学显著性的区分性EEG生物标志物特征(例如,N1幅值的差异)。保存输出结果以供后续步骤使用。
  6. 基于文献提出关于药物作用机制及相应模型参数的假设。查阅先前文献和实验数据,识别神经治疗可能改变的生物物理特性,这些特性可用于解释特征差异。
  7. 确定生物物理神经网络模型(HNN)中哪些参数直接表示或间接关联于步骤1.6中识别出的生物学特性。将这些参数定义为感兴趣参数。结合先前文献和HNN文档,将生物学机制映射到模型参数。
  8. 输出一组与假设可产生已识别EEG特征差异的生物物理特性相对应的模型参数。以默认HNN模型(在步骤2中初始化)作为所有参数值的起点,并保存输出结果用于后续步骤。

2. 初始化默认的 HNN 模型:安装建模软件并设置项目文件夹

注意:本研究中使用的软件版本及其最低系统要求详见材料表。Linux、macOS 和 Windows 系统均提供多种安装方式(例如 pip、conda 和源码安装)。

  1. 下载并安装可用版本的 Anaconda Python。创建并激活一个新的 Python 环境,用于安装所需的软件包。
  2. 根据操作系统特定的安装说明,安装生物物理神经建模 HNN-core 软件,说明详见 https://jonescompneurolab.github.io/textbook/content/01_getting_started/installation.html。
    注意:为了高效安装本研究中使用的软件依赖项,相关代码仓库(https://github.com/ntolley/hnn_jove)使用了 pixi(https://pixi.prefix.dev/latest/)。请按照仓库 README 文件中的说明安装 pixi,并设置本地代码仓库副本。
  3. 在终端中输入以下命令,验证已安装的生物物理神经建模软件版本是否为 0.6.0 或更高版本:pip show hnn_core
  4. 确保已激活 Python 环境且安装成功完成。在终端中输入 hnn-gui 并按 Enter 键,启动图形用户界面(GUI)。
  5. 在计算机文件系统中创建一个新的项目文件夹,用于存储本方案生成的所有数据文件。将该文件夹创建在易于访问的目录中(例如,主目录或当前项目工作目录)。

3. 通过手动调节建立预处理模型拟合

  1. 从标准的HNN ERP模拟及其默认参数开始。手动调整缩放因子和外源性驱动参数,以拟合治疗前ERP(例如,治疗前ERP)。
    注意:HNN 图形用户界面(GUI)会自动加载拟合体感事件相关电位(ERP)的模型参数40,这一点已通过大量研究被证实为一个良好的“典型ERP”起点。本教程重点介绍从该起点出发,调整缩放因子和外源性输入参数的方法。
  2. 将步骤1中获得的预处理经验性ERP波形载入HNN图形用户界面(GUI)图 4A–4F)
    1. 单击位于 GUI 窗口左下部分菜单栏上的“Load data”按钮(图 4D).
      注意:ERP 峰命名的术语在文献中差异较大;P1/N1/P2 标签的使用 图4F 仅用于说明目的,可能与其它研究中使用的命名惯例不符。
    2. 在文件浏览器窗口中,选择一个包含待建模ERP波形(即目标波形)的.csv或.txt文件。确保该文件为逗号分隔,且包含两列数据:第一列为时间(ms),第二列为源定位的实测偶极子波形(nAm)。首行将被视为标题行,不应包含数据值。可选择性地包含具有信息性的列标签(例如,“Time (ms)”和“Dipole (nAm)”)。
      注意:实证数据文件的名称为 预处理.txt 在相关代码仓库中。
    3. 检查在图形面板中自动绘制的波形图4F).
  3. 运行典型事件相关电位的默认模拟
    1. 设置参数值 tstop, dt, 试验, 后端,以及 核心 在“模拟参数”面板中(图4B)调整至所需数值。使用 tstop 以控制模拟时长, dt 以控制积分时间步长,并 试验 用于控制使用相同模型参数值运行的重复模拟次数。选择 后端 无论是串行(Joblib)还是并行(MPI),并指定计算机数量 核心.
      注意:不同试验间的变异性来源于下文步骤 3.5 中所述的外源性诱发驱动时序的标准差。
    2. 点击 运行 按钮(图 4D) 开始典型事件相关电位(ERP)的默认模拟。
  4. 绘制一个将模拟ERP与实证ERP进行比较的图表
    1. 点击 可视化 位于GUI窗口左上角的选项卡图 4A).
    2. 单击标有下拉菜单 用于比较的数据 (未显示)并从步骤 3.2 中选择已加载的目标波形。
    3. 点击 清晰轴 重置图表。
    4. 点击 添加图表 生成一个新图表,叠加显示模拟的初始ERP波形(蓝色)和目标波形(橙色),并附有文字说明,标示出两个波形之间自动计算的相关系数(Corr)和均方根误差(RMSE)。图4F).
      注意:HNN-GUI 提供计算两种拟合优度指标的选项:Corr 和 RMSE。这些指标用于手动调参和优化(步骤 4)。
  5. 调整缩放系数
    1. 通过手动调节缩放系数,使模拟的和实测的偶极子波形幅度大致匹配。将默认值设置为 偶极子缩放 参数(图4C模拟”选项卡中的)图 4A) 至 3000。
      注意:缩放因子对应于对产生脑电图信号的神经元数量的估计值的预测。默认值3000表示200个锥体神经元(HNN模型的大小) × 3000 = 生成 y 轴上以 nAm 为单位所示幅度的诱发电位需要 600,000 个神经元 图4F.
  6. 调整外源性驱动的时间
    1. 通过手动调节外源性输入的均值和标准差,以更精确地拟合实证记录的治疗前ERP峰值(即P1/N1/P2)的时间特征。图5A–5D).
      注意:HNN 自带的默认局部连接性和细胞参数经过调校,可重现健康的单细胞及网络水平活动模式。尽管可以调整局部网络参数,但建议初始时保持预调校的 HNN 新皮层模板模型参数不变,仅通过调整外源性驱动来测试是否能够实现可靠的拟合结果。
    2. 识别哪些模拟的ERP峰值在时间上与实证ERP波形不一致图4).
      注意:本示例假设经验性事件相关电位(ERP)中存在三个早期峰,如默认的典型ERP模拟所示。如需增加峰,可模拟额外的外部驱动。
    3. 点击 外部驱动器 位于GUI窗口左上角的选项卡图 4A图 5A).
      注意:三个预定义的外源性驱动参数可见,分别代表前馈近端驱动(evprox1)、反馈远端驱动(evdist1)和重新出现的前馈近端驱动(evprox2),用于生成默认的标准事件相关电位(ERP)模拟(参见 引言 有关HNN模型和外源性驱动结构的详细信息,请参见)。显示放电时间和计数的直方图如图所示。 图4E.
    4. 点击外源驱动的下拉菜单 平均时间 最接近错位峰。
    5. 修改文本框中的数值 平均时间标准差时间 以更好地匹配目标波形中峰的出现时间和宽度图5B–5D)。调节 平均时间 调整峰时间 标准偏差时间 调整峰宽。
      注意:平均时间(Mean time)和标准差时间(Std dev time)用于控制外源性脉冲激活局部网络时的均值和方差,这些脉冲以近端或远端投射模式呈现(参见直方图) 图4E)。这些参数并不能完全决定事件相关电位(ERP)的峰值潜伏期或波宽。确切的潜伏期和波宽取决于外部刺激输入和神经网络内在活动的共同作用。
      1. 设置 平均时间 将 evprox1 外接驱动器设置为 60 毫秒。
      2. 设定 平均时间 将 evdist1 外接驱动器设置为 100 毫秒。
      3. 设置 平均时间 将 evprox2 外部驱动器设置为 150 毫秒。
  7. 调节外源性驱动的强度
    1. 通过手动调节外源性输入的突触权重(突触后电导),以获得与实证记录的事件相关电位峰值(如 P1/N1/P2)幅度更接近的拟合效果图6A图6B).
    2. 识别哪些模拟的ERP峰值在幅度上与实测ERP波形不一致。
    3. 点击 外部驱动器 位于GUI窗口左上角的选项卡图 4A).
    4. 点击外源驱动的下拉菜单以选择其 平均时间 最接近错位峰。
    5. 修改下方文本框中的数值 AMPA权重NMDA 权重 以调节突触电导。增强对L5和L2/3锥体神经元的近端驱动强度通常会产生更正向的峰电位,而增强远端驱动强度则通常会产生更负向的峰电位。
      注意:与外源性驱动时序类似,事件相关电位(ERP)峰值幅度并不完全由驱动强度决定。神经元放电动力学可能产生非直观效应。建议在1个数量级范围内测试参数变化(例如,将L5_pyramidal的AMPA参数从0.014调整至0.14),并进行迭代优化。 图6 显示设置为 10 的数值× 小于默认模拟尺寸。
      1. 设置 AMPA权重 evdist1 驱动 L5_pyramidal = 0.014243 且 L2_pyramidal = 0.0000007。
      2. 设置 NMDA 权重 evdist1 驱动 L5_pyramidal = 0.0080074 且 L2_pyramidal = 0.0004317。
      3. 设置 AMPA权重 evprox2 驱动 L5_pyramidal = 0.0684013 且 L2_pyramidal = 0.143884。
        注意:生成代表性结果所用的完整参数集可在相关代码仓库中获取(https://github.com/ntolley/hnn_jove;参见 data/opt_baseline_config_correlation_best.json)。建议用户加载该配置文件以及提供的数据文件(data/pre-treatment.txt 和 data/post-treatment.txt),并参考 notebooks/ 目录中的示例工作流程,以复现报告中的模拟结果。
  8. 保存修改后的模拟设置。
    1. 完成步骤 3.5–3.7 的修改后,单击 模拟 标签图 4A)并在其中输入“pre-treatment_handtuned” 名称 文本框图4B).
  9. 运行修改后的模拟
    1. 点击 运行 单击按钮以模拟修改后的参数集。
    2. 检查图层面板中生成的图表(图4F图7A–7D)。使用相应的图表选项卡(例如,“图 1”和“图 2”)访问之前的图表。
  10. 反复手动调节
    1. 继续进行迭代式手动调谐以提高相关系数。
    2. 重复步骤 3.4,以目标波形重新绘制模拟结果并重新计算相关系数。
  11. 保存最终模型输出。
    注意:在保存模拟输出后,可暂停本实验方案。通过将已保存的配置文件加载到软件中以继续实验。
    1. 点击 保存网络 单击按钮,将最佳拟合参数集保存为名为“pre-treatment_handtuned.json”的 .json 文件。
    2. 点击 保存模拟 按钮以保存一个名为“pre-treatment_handtuned.txt”的 .txt 文件,其中包含模拟的偶极子波形图 4D).
    3. 将两个文件移至步骤 2.5 中创建的项目文件夹中。确保文件名与下拉菜单中的模拟名称一致。
      注意:文件将保存到用于运行图形用户界面的网络浏览器的默认下载目录中。请手动移动文件,或临时更改浏览器的下载目录。

figure-protocol-1
图4.经典HNN模拟ERP波形与实证治疗前ERP的比较。A)通过图形用户界面(GUI)选项卡可访问的参数类别。(B)控制模拟时长和试验次数的模拟参数。(C)控制波形显示的可视化参数。(D)用于加载数据、运行模拟和保存输出的模拟控制面板。(E)尖峰直方图,显示经典ERP模拟中外部驱动输入的分布情况。(F)经典ERP模拟的偶极子波形(蓝色)与Kohl等人43提供的实证听觉ERP(橙色)叠加。初始模拟与数据拟合不佳,峰值时间和幅度未对齐(相关系数 < 0.95)。用于生成实证ERP的实验范式详见Kohl等人的研究43:将音调(1 kHz,持续时间50 ms,10 ms淡入/淡出)交替呈现于左耳和右耳,刺激间间隔为0.8–1.2 s,强度为受试者听觉阈值以上60 dB。请点击此处查看该图的放大版本。

figure-protocol-2
图 5.调整外源性驱动时序以对齐 ERP 峰值。A)图形用户界面中的“外部驱动”选项卡,用于配置模型的诱发输入。(B–D)调整各个外源性驱动的平均时间参数,以使模拟的 ERP 峰值与实证数据对齐。具体而言,(B)近端驱动 evprox1 对齐至约 60 ms,(C)远端驱动 evdist1 对齐至约 100 ms,以及(D)近端驱动 evprox2 对齐至约 150 ms。调整平均时间参数(已突出显示)可改变模拟峰值的出现时间,从而提高与实证波形的一致性。这些调整有助于改善对齐效果,并提高与目标 ERP 的相关性(见图 7B)。请点击此处查看该图的放大版本。

figure-protocol-3
图6。通过调节外源性驱动强度来调整事件相关电位(ERP)峰值幅度。AB)在图形用户界面(GUI)的“外部驱动”选项卡中调节α-氨基-3-羟基-5-甲基-4-异噁唑丙酸(AMPA)和N-甲基-D-天冬氨酸(NMDA)受体的突触权重。(A)调节远端驱动(evdist1)的突触权重,包括作用于第2/3层(L2/3)和第5层(L5)锥体神经元的AMPA和NMDA电导。(B)调节近端驱动(evprox2)的突触权重,主要影响锥体神经元中的AMPA电导。在此示例中,突触权重相对于默认值降低了10倍,导致ERP峰值幅度减小,并与实测波形具有更好的一致性(见图7C)。请点击此处查看该图的放大版本。

figure-protocol-4
图7.手动调节与优化以拟合模型参数。所有模拟均显示5次试验,其中平均ERP为深蓝色,单次试验为浅蓝色。(A)标准ERP模拟(蓝色)与治疗前ERP(橙色)叠加对比。(B)调整外源性驱动时序可改善峰值对齐。(C)降低突触权重可减小峰值幅度。(D)自动优化可实现与实测波形的高度拟合(相关系数 = 1.0),并包含诱发驱动时序的更大变异性。请点击此处查看该图的放大版本。

4. 通过参数优化建立预处理模型拟合

注意: 图形用户界面(GUI)目前不支持对优化过程中的随机种子进行控制。如需进行可重复的优化运行,请使用 Python API。相关代码仓库中包含一个示例实现(参见 code/baseline_optimization.py),可通过向优化函数传递 seed 参数来设置固定的随机种子(例如,optim.fit(..., seed=123))。

注意:本示例展示如何优化目标参数,以使用CMA-ES方法估计出能够与波形高度拟合的单个数值(请注意,这与SBI不同;两者均为拟合模型参数的方法,但SBI的主要输出是参数分布)。关于如何估计可解释波形的参数分布的方法示例,见结果部分。对于治疗前的事件相关电位(ERPs),首先在假设默认HNN新皮层模型中的细胞和局部网络连接参数固定的前提下,优化外源性驱动参数。第7步中描述的HNN所提供的多尺度预测,为此假设的验证提供了目标依据。随着新信息的获取以进一步约束模型预测,HNN框架允许对任意参数集进行估计。

  1. 打开优化设置
    1. 点击图形用户界面(GUI)左上角的“优化”选项卡(图8A)。
    2. 配置优化运行的设置,包括迭代次数、求解器和目标函数。
      注:默认优化设置(目标函数 = “dipole_corr”;求解器 = “cma”)适用于ERP波形。该目标函数最大化模拟波形与实测波形之间的相关系数。若需优化多个参数,可增加最大迭代次数。相关系数为无量纲度量,因此当使用“dipole_corr”时,应在优化后调整缩放因子(步骤4.7.1)。或者,可使用“dipole_rmse”以最小化均方根误差(RMSE),此时缩放因子保持固定。
    3. 点击“最大迭代次数”文本框,并输入100。
  2. 选择用于优化的参数
    1. 点击将要优化其参数的外源性驱动项的下拉菜单(图8A图8B,红色圆圈)。
    2. 通过点击“是否用于优化?”下方的复选框,选择要优化的驱动参数(图8B)。
  3. 定义参数约束条件
    1. 在“约束条件(%)”下的“最小值”和“最大值”文本框中输入数值,以指定优化器搜索的参数值范围(图8B)。
      注:对于已有较高相关系数(Corr > 0.9)的模拟,默认的20%范围是合适的。例如,对均值时间为65.53 ms的参数应用20%的范围,将产生52.42–78.64 ms的边界。对于初始拟合较差的情况,可增大最小值和最大值的百分比;但所需模拟次数可能会显著增加。
  4. 运行优化
    1. 点击“运行优化”按钮(图8A),执行优化程序。
  5. 保存优化结果
    1. 点击“保存优化历史”按钮(图8A)。
    2. 将保存的文件移至步骤2.5中创建的项目文件夹。
      注:优化结果可被存储并重复使用。本方案可在该阶段暂停,并通过加载已保存的优化历史记录继续执行。
  6. 评估优化质量
    1. 评估优化运行的质量。
      注:当使用相关系数作为拟合优度指标时,建议采用 Corr > 0.95 作为停止标准,因为这通常表示模拟波形能够再现目标ERP的主要峰和谷。目前尚不支持提前停止功能,但正在开发中。如果未达到停止标准,但每10次迭代损失仍在持续下降,则应增加迭代次数。
  7. 根据优化结果确定下一步操作
    1. 如果对治疗前ERP实现了良好拟合(即 Corr > 0.95),则重新调整缩放因子以最小化RMSE,并进入步骤5。
      注:如步骤4.1所述,当使用“dipole_corr”作为目标函数时,应在优化后重新调整缩放因子。在此示例中,缩放因子从默认的3000×(图7A7C)降低至1000×(图7D)。
    2. 如果优化未能实现对治疗前ERP的良好拟合,则返回步骤4.2,通过增加最大迭代次数、改进手动调优的起始点,或选择其他可调参数进行故障排查。
      注:有关故障排查步骤的详细说明,请参见“讨论”部分中的“将参数拟合到数据特征时的故障排查”一节。

figure-protocol-5
图8.优化外源性驱动参数以改善与治疗前ERP的拟合度。A)用于配置优化参数的图形用户界面中的优化选项卡。(B)优化所选参数及其约束范围。(C)优化结果示例,显示与Kohl等人43的实证ERP数据的拟合度得到改善。(D)优化损失曲线,显示约经过80次迭代后收敛。请点击此处查看该图的放大版本。

5. 建立治疗后模型拟合

  1. 从优化的治疗前ERP模拟开始。手动调整并优化相关参数,以拟合治疗后的ERP。
  2. 加载治疗后实证ERP波形
    1. 将第1步中得到的治疗后实证ERP波形加载到图形用户界面(GUI)中(操作方法与步骤3.2相同;图9A)。
  3. 加载优化后的治疗前参数
    1. 将第1–4步中获得的优化治疗前ERP参数作为起点进行加载(图9A)。
  4. 执行手动调整与优化
    1. 对第1.7步中确定的关注参数,执行手动调节和参数优化(操作流程与步骤3.2–3.11及步骤4相同)。
    2. 持续进行调节与优化,直至模拟ERP与治疗后ERP之间达到高相关性(相关系数 Corr > 0.95)。
      注:为说明目的,在图9B中,对手动调节应用于一个信号靶向参数(降低局部网络GABAB最大电导),从而更贴近治疗后数据。此处未进行系统优化,以评估该参数变化对数据的解释程度。“代表性结果”部分将描述如何使用SBI方法估计多个假设为治疗后关注参数的参数分布。建议在严谨研究中采用SBI(详见步骤6),因为它可估计能够解释ERP波形的参数分布,从而实现不同参数拟合之间的稳健比较。
  5. 保存模型配置并比较参数
    1. 保存模型配置,并比较关注参数在治疗前与治疗后条件下的优化值(数据未显示)。
    2. 重复步骤3.11,导出模型参数的.json文件,并将该文件移至第2.5步创建的项目文件夹中。
    3. 通过点击“加载外部驱动”(图5A),选择治疗前或治疗后网络配置文件,查看外源性驱动参数。
    4. 通过点击“加载局部网络连接”(图9C),选择治疗前或治疗后网络配置文件,查看局部网络参数。
    5. 识别治疗前与治疗后网络配置之间参数值的变化,并将这些变化解释为基于模型预测的治疗后生物标志物机制。

figure-protocol-6
图9。评估γ-氨基丁酸B型(GABAB)突触强度作为治疗后脑电图生物标志物的机制。A)治疗前优化模拟结果(蓝色)与治疗后事件相关电位(红色)叠加,显示峰值幅度降低。(B)GABAB突触强度的降低导致N1波幅减小,提示其可能为潜在机制。(C)图形用户界面面板,显示局部GABAB突触强度被修改的位置。请点击此处查看该图的放大版本。

6. 使用 SBI 进行不确定性量化,并通过 HNN-Python 应用程序编程接口评估可分性

注意:SBI 需要安装一个独立的 Python 软件包62。有关如何使用 SBI 软件包在 HNN 中运行参数推断的代码示例,请参考相关代码仓库(https://github.com/ntolley/hnn_jove)。该代码的组织结构遵循后续实验方案中的各个步骤。关于将 SBI 应用于 HNN 模型的完整讨论见文献55

  1. 安装 SBI 软件包
    1. 在激活 Python 环境的终端中运行以下命令来安装 SBI 软件包:pip install sbi。
  2. 定义先验参数范围
    1. 确定围绕目标治疗前和治疗后 ERP 参数子集的参数范围,以创建有界的先验分布,用于不确定性量化。
  3. 生成训练数据集
    1. 定义一个参数更新函数(方法与参数优化相同)。
    2. 固定从先验分布中生成样本的随机种子,以确保可重复性。如果使用 NumPy 生成随机样本,请在 Python 脚本中创建一个随机数生成器实例(例如 rng = np.random.default_rng(123)),并使用该生成器进行采样。
      注意:相关代码仓库(https://github.com/ntolley/hnn_jove)在 code/generate_simulations.py 中提供了使用 NumPy 随机生成器的示例。
    3. 从先验分布中采样参数。
      注意: 本研究使用了 10,000 个样本以生成代表性结果。
    4. 使用采样的参数值生成一组模拟 ERP 数据集。
  4. 选择摘要统计量
    1. 选择一种能够表征 EEG 波形的摘要统计量。
      注意:摘要统计量是指任何能够捕捉 EEG 波形关键特征的量化指标。常见的选择包括峰值时间和幅值。本文中使用主成分分析(PCA)提取摘要统计量(即前四个主成分的载荷)。详见55中的完整讨论。
    2. 训练 SBI 网络
      注意: 本教程使用 SBI 软件包中为神经后验估计器对象提供的默认训练参数(例如 density_estimator="maf",training_batch_size=200,learning_rate=0.0005)。训练参数的详细说明见 SBI 文档(https://sbi.readthedocs.io/en/stable/api_reference/_autosummary/sbi.inference.NPE_B.html)。
    3. 在导入 torch 后,在 Python 脚本中加入 torch.manual_seed(0),以设置全局 PyTorch 随机种子,确保训练过程可重复。
    4. 训练 SBI 网络,将参数组合映射到模拟的 ERP 波形。
      注意:训练好的 SBI 网络是一个 Python 对象,它以 EEG 数据的摘要统计量作为输入,并输出参数的分布(后验分布)。如果训练成功,在 HNN 模型中使用该分布生成的参数进行模拟,将产生与实证数据相似的 EEG 波形(后验预测检验 [PPC])。
    5. 生成后验样本并评估拟合效果
    6. 将实验获得的 EEG 波形作为条件输入提供给训练好的网络。
    7. 从基于实验 EEG 波形条件化的后验分布中抽取参数样本。
    8. 对从后验分布中抽取的参数样本进行模拟。
    9. 计算模拟波形与作为输入提供的实验 EEG 波形之间的相似性。
      注意:此过程称为后验预测检验(PPC)。训练良好的网络应产生与实证波形高度匹配的模拟结果(高相关性或低 RMSE)。如果 PPC 未能产生满意的模拟结果,可能存在两种情况:(1)假设的机制无法解释生物标志物,需要提出新的假设并更新先验分布;或(2)SBI 网络未成功训练。此时,应增加训练预算或修改摘要统计量。
    10. 如果从采样参数分布模拟出的 ERP 能够较好拟合治疗前和治疗后的 ERP(PPC 相关系数 > 0.95),则进入第 6.8 步;否则进入第 6.7 步。
  5. 排查 SBI 网络训练问题
    注意:PPC 失败表明 SBI 网络的训练参数需要调整。详见“讨论”部分中的“拟合参数至数据特征时的故障排除”一节的详细说明。
    1. 增加训练数据集的规模。
    2. 修改摘要特征。
    3. 选择不同的 SBI 架构进行训练。
  6. 可视化后验分布并评估可分性
    1. 将第 6.6.2 步得到的参数样本数组传入 pairplot 函数,并为对应于各 ERP 条件的分布分配不同的颜色。
      注意:相关代码仓库展示了绘图功能,可用于复现 图 10
    2. 检查生成的 pairplot 对角线面板中是否存在非重叠的分布。通过计算 OVL 来评估可分性(图 10A)。分布高度分离的参数(OVL < 0.1)对应于预测的神经治疗作用机制,这些机制在治疗后相对于治疗前发生了变化。
      注意:OVL 是一种量化分布可分性的指标,取值范围为 (0,1),其中 OVL = 0.0 表示无重叠,OVL = 1.0 表示完全重叠54,55。相关代码仓库中提供了计算 OVL 的代码。

figure-protocol-7
图10用于参数不确定性量化及神经治疗机制识别的SBI方法 (A) 使用SBI估计的参数分布的成对图可视化。对角线面板(i–iv) 显示各个参数的单变量分布,包括(i丘脑皮层同步性,(ii树突 Km 电导,(iiiγ-氨基丁酸(GABA)B 电导率,以及(iv) 皮质-皮质反馈强度。单位为 (i) 以默认(治疗前)参数值的乘法缩放因子形式表示。单位为 (ii-iv) 以默认(治疗前)参数值的对数尺度上的乘法缩放因子表示。 治疗前(蓝色)和治疗后(红色)条件的分布显示出不同程度的可分性,其中丘脑皮层同步性重叠程度最低(重叠值,OVL = 0.07),表明其治疗相关效应最强。非对角线面板显示参数间的双变量关系。B) 治疗前ERP的后验预测检验(PPC);模拟波形(黑色)与实测数据(蓝色)高度吻合。 (C治疗后ERP的PPC;模拟波形(黑色)与实证数据(红色)高度吻合。 请点击此处以查看此图的放大版本。

7. 进行检查、验证及进一步的模型约束

注意:此步骤提供了如何在图形用户界面(GUI)中可视化模拟活动各要素的示例。这些多尺度细节可作为后续实验中验证和指导模型预测结果的目标7,47。本方案未提供关于选择哪些预测结果最适合用于验证实验的指导,也未说明验证实验应如何进行(即步骤 7.3)。

  1. 加载模型参数并运行模拟
    1. 加载针对治疗前和治疗后条件优化的模型参数,并运行模拟。
      注意:可加载并查看第4–5步中优化得到的参数。相关GitHub代码库中包含了如何通过Python接口导出第6步中SBI生成的网络参数的示例。
  2. 检查多尺度预测结果
    1. 分析来自模拟输出的多尺度预测结果。
    2. 绘制细胞水平的放电活动图
    3. 点击可视化标签页(图4A)。
    4. 点击标记为布局模板的下拉菜单,选择偶极子分层-放电(Dipole Layers-Spikes)。
    5. 在数据集下拉菜单中,选择要绘制的模拟结果。
    6. 点击生成图表以可视化对偶极子波形有贡献的放电活动。
      注意:某些微环路特征(例如,LFP和CSD)仅可通过HNN-Python应用程序编程接口(API)获取。有关这些功能的基于代码的教程可在HNN示例页面获取(https://jonescompneurolab.github.io/hnn-core/stable/index.html)。
  3. 使用实证数据验证模型预测
    1. 识别现有数据集和/或收集新的实证数据(例如,侵入性电生理学、层析MEG/EEG以及磁共振波谱)以检验多尺度模型预测。
    2. 将多尺度模型预测结果与实证数据集进行比较。
    3. 如果多尺度预测结果与实证数据集一致,则认为该模型针对所选微环路特征已通过验证。
    4. 如果多尺度预测结果与实证数据集不一致,则使用新的实证数据约束默认的HNN网络进行更新,并返回第3步。

结果

本部分介绍了一种利用HNN建模软件研究作用机制未知的神经治疗手段的场景。目标是利用治疗前和治疗后的脑电图(EEG)信号,预测该神经治疗手段对神经环路的改变方式。所展示的结果仅为示例,旨在说明如何应用HNN建模来探究神经治疗的作用机制。

建立脑电图事件相关电位生物标志物的机制假说(步骤 1)

在此示例中,采用一种假设的感觉事件相关电位(ERP)范式来研究神经治疗如何改变信号(步骤 1)。图 1A 展示了治疗前的听觉 ERP(蓝色)与假设的治疗后 ERP(红色)(另见 图 9)。治疗前的听觉 ERP 数据来自 Kohl 等人实验记录的源定位结果43,而假设的治疗后 ERP 是通过对治疗前波形应用高斯 taper 窗函数进行缩放生成的。如图所示,与治疗前 ERP 相比,该假设的神经治疗导致 P1、N1 和 P2 成分的幅值显著降低。

请注意,预处理ERP数据来自Kohl等人43的研究,该研究中使用的HNN模拟模型相较于默认的HNN模型,为锥体神经元引入了更为真实的钙通道动力学。因此,Kohl等人43中的模拟结果与本文所示结果略有差异。Kohl等人2020年的模型(以及其他更新的HNN模型)可通过Python API获取(https://jonescompneurolab.github.io/hnn-core/stable/generated/hnn_core.calcium_model.html#hnn_core.calcium_model)。目前,通过图形用户界面(GUI)访问此类扩展模型的功能正在开发中。

接下来,确定代表治疗相关效应的模型参数(即感兴趣参数),这些参数被假定用于解释神经治疗如何降低P1、N1和P2的幅度(步骤1.6–1.8)。候选神经机制(及相应模型参数)的主要类别包括外源性突触输入的时间、局部神经元离子通道电导、局部突触连接性以及外源性突触连接性(图1B)。在此示例中,使用HNN评估每个类别中的候选机制,以分析这些参数的变化如何影响模拟的事件相关电位(ERP)。

关注的参数

  1. 第一(丘脑皮层)近端驱动的标准差(即丘脑皮层同步性),代表初始前馈感觉输入同步性的变异性。
  2. 第5层(L5)锥体神经元中的毒蕈碱型钾(Km)通道电导,调控神经元兴奋性,电导增加时兴奋性降低。
  3. 局部GABAB受体强度,对应于中间神经元向局部网络中所有细胞传递的慢抑制性突触。
  4. 反馈(皮层-皮层)远端驱动的电导强度,代表约100毫秒感觉诱发反馈输入在上颗粒层AMPA和NMDA突触上的强度。

建立治疗前ERP模型拟合度(步骤3–4)

通过执行步骤3–4来模拟治疗前ERP(最终的治疗前模拟结果如图8C所示)。成功的模拟结果表现为模拟波形与实测波形高度吻合,可通过较高的相关系数和较低的均方根误差(RMSE)进行量化评估。

建立治疗后ERP模型拟合(步骤5)

以预处理ERP模型为起点,通过手动调整和参数优化,判断所关注的参数是否能够重现实测的治疗后ERP。若拟合成功,则表明假设的参数足以解释ERP波形中与治疗相关的改变。

使用SBI进行不确定性量化(步骤6)

由于生物物理模型中固有的参数简并性,使用SBI(步骤6)对治疗前到治疗后的参数变化进行预测时,不确定性量化至关重要。SBI的一个关键前提是实现对治疗前和治疗后ERP的精确拟合(步骤3–5)。如果未能实现精确拟合,SBI生成的后验样本可能无法重现实测波形,从而导致预测结果不可靠。

如果在步骤3至5中无法实现成功的拟合,则在应用SBI之前,需重新调整所选参数及其先验范围。

在此示例中,SBI 仅应用于四个感兴趣的治疗后参数,而所有其他参数均保持固定。尽管将 SBI 应用于更大的参数集可提高鲁棒性,但会显著增加计算成本(见讨论)。

SBI 用于估计能够生成与目标波形高度匹配的模拟事件相关电位(ERP)的完整参数分布。简而言之,SBI 是一种贝叶斯推断方法,通过训练神经网络将模型输出映射到模型参数的分布52,53,55。训练后的网络随后被应用于实测波形,以推断与数据一致的参数分布。该过程需要对参数范围作出先验假设。

在此示例中,对四个感兴趣的参数定义了均匀先验分布:丘脑皮层同步性、锥体神经元树突 Km 电导、局部 GABAB 电导以及皮层间反馈强度。先验边界定义为默认值的标量倍数:丘脑皮层同步性为 0–5×,其余参数为 10−1–101×。

图10A 显示了使用配对图(pairplot)可视化展示的治疗前和治疗后ERP的参数分布结果。对角线面板显示单变量分布,非对角线面板显示双变量关系。机制性预测对应于不同条件下参数分布明显分离的情况。

单变量分布的检验显示,丘脑皮层同步性在治疗前后的可分性最大(重叠率OVL最低,为0.07),且在治疗后增加(图10A(iii),红色)。这表明HNN框架预测丘脑皮层同步性的调节可能是其潜在的作用机制。

后验预测验证

使用后验预测检验(PPC)验证推断出的参数分布。从后验分布中生成独立的参数样本,并模拟相应的事件相关电位(ERP)。当模拟的波形与实测ERP高度匹配时,表明PPC成功。

如图所示 图10B图10C,治疗前(图10B,蓝色)和治疗后(图10C,红色)波形与基于后验样本生成的模拟结果(黑色)高度吻合,相关系数分别为0.99和0.96(在10个独立样本上取平均)。这些结果表明,所推断的参数分布能够准确重构波形。

补充图1中提供了一个不成功的后验预测检验(PPC)示例。该示例的结构与图10相同,并使用了相同的经训练的SBI网络;然而,此处采用了一种在训练集中未充分表示的治疗后波形(例如,在N1潜伏期具有正向偏转的ERP波形)。在补充图1C中可观察到失败的PPC,其相关系数较低(例如,Corr < 0.95)。值得注意的是,补充图1A中的后验分布显示出参数分布高度分离。若未进行PPC,这些结果可能被误解释为治疗前与治疗后条件之间存在有意义的差异。此示例强调了在解释后验分布的同时进行PPC的重要性,因为PPC失败的结果不可靠,不应进一步分析。

模型检验与验证(步骤 7)

利用HNN模型,可以直接检查和可视化每个ERP模拟(步骤7.2.2)背后的细胞和环路水平活动,例如放电活动。图11A图11B 展示了从治疗前和治疗后参数分布中采样的模拟ERP,以及相应的细胞特异性放电活动(图11C图11D)。

figure-results-1
图11。脑电生物标志物生成背后的细胞水平放电活动。A)治疗前事件相关电位(蓝色),以及一次后验预测模拟结果(黑色)。(B)治疗后事件相关电位(红色),以及相应的后验预测模拟结果(黑色)。(C)治疗前事件相关电位背后的模拟放电活动。(D)治疗后事件相关电位背后的模拟放电活动。请点击此处查看该图的放大版本。

波形以未平滑的形式可视化,以突出神经元放电时间对电流偶极子的贡献。在实验性脑电图(EEG)信号中,大量神经元群体产生空间平均的信号,因而显得更加平滑。由于HNN模拟的神经元群体较小(200个锥体神经元),需使用平滑处理来近似更大规模的神经活动(>100,000个神经元)。

不同条件下一个显著差异是,治疗后L5锥体神经元的放电活动减少(图11C图11D,红点所示)。请注意,图11 仅显示了后验分布中的一个样本;应分析多个样本以生成可靠的预测结果。这些结果表明,该假设性神经治疗药物改变了多尺度神经环路活动,导致P1–N1–P2波幅降低。

此类预测可通过侵入性电生理学方法(例如,高密度层析探针记录)或其他成像方式(步骤 7.3)直接进行验证。新获取的数据可用于进一步约束模型的预测结果。尽管本方案侧重于拟合宏观尺度的脑电图(EEG)数据以推断微环路活动,但该框架亦可反向应用,即通过拟合微环路数据(例如,神经元放电、局部场电位/电流源密度)来推断宏观尺度的脑电图信号。

补充图1. SBI工作流程中后验预测检验失败的示例。 图表的排列方式与图10相同。治疗前数据(蓝色)与图10一致。假设的治疗后数据生成方式与之前相同(波形乘以高斯 taper 窗口),但经过变换后产生了一个在HNN模拟训练集中未被良好表征的正向峰。(A) 使用SBI估计的参数分布的成对图可视化。对角线子图(i–iv)显示各个参数的一维分布,包括(i)丘脑皮层同步性、(ii)树突Km电导、(iii)GABAB电导和(iv)皮层间反馈强度。治疗前(蓝色)与治疗后(红色)条件下的分布显示出所有参数的高度可分性(OVL < 0.1)。非对角线子图显示参数间的双变量关系。(B) 治疗前ERP的后验预测检验(PPC);模拟波形(黑色)与实测数据(蓝色)高度吻合。(C) 治疗后ERP的后验预测检验(PPC);模拟波形(黑色)与实测数据(红色)差异显著,相关系数Corr < 0.95,表明PPC失败。请点击此处下载该文件。

讨论

对脑电图生物标志物进行计算神经建模,可能有助于深入理解中枢神经系统药物如何重构神经回路,并对治疗效应背后的生物学过程提供预测。本文所展示的工作流程表明,如何将一种常用的脑电图生物标志物——听觉事件相关电位(ERPs),结合使用HNN进行的生物物理建模,作为探究药物影响神经活动机制的窗口。通过将宏观尺度的脑电图测量结果与潜在的细胞和回路水平过程相联系,该方案为机制性解释提供了结构化且基于假设的分析框架。重要的是,该方法不仅限于事件相关电位,还可扩展用于研究其他局部脑电图信号,包括低频神经振荡40,63和瞬态谱事件7,47,64,从而拓宽其在电生理生物标志物和实验范式中的适用范围。

与其他用于脑电图(EEG)神经建模的框架相比,HNN 在模型复杂性与计算效率之间实现了良好平衡,特别有利于进行迭代式假设检验。例如,《The Virtual Brain》能够模拟生成时空脑电图信号的大尺度脑网络34,65。然而,为了实现全脑建模,神经活动通常采用简化的数学表达形式,这省略了诸如锥体神经元形态等详细的细胞特征,限制了模型参数与药物作用细胞机制之间的直接关联能力。相反,具有高度形态学和生理学细节的大规模模型虽能以较高的生物学真实性模拟脑电图信号66,67,68,69,但其计算成本显著,通常需要数小时的计算时间才能模拟几秒钟的神经活动。这种巨大的计算负担可能限制其可及性,并减缓假设生成与验证所需的迭代过程。HNN 处于中间位置(图 2),能够在保持足够生物学细节的同时模拟局部新皮层回路,从而生成细胞层面和回路层面的预测,且仍具备较高的计算效率(即模拟耗时在数秒量级),因此非常适合整合到实验工作流程中。

尽管具有上述优势,在应用脑电图(EEG)和生物物理神经建模来研究脑部疾病及药物作用机制时,仍需考虑若干局限性。产生EEG信号的生物物理细胞与环路特性并不能涵盖药物干预所影响的全部生物学过程。例如,全身性或免疫学反应可能不会直接影响EEG信号,因此在模型输出中可能无法体现。此外,机制性假说通常源于动物研究,而这些研究结果未必能完全外推至人类大脑功能,尤其是在神经精神疾病领域,其临床结局依赖于行为和认知评估70,71。另一个重要挑战在于区分药物的急性与慢性药理效应。虽然急性药物-受体相互作用已有相对明确的特征描述,但长期药物暴露所诱导的适应性变化仍了解甚少,当前的建模框架可能无法充分捕捉这些变化。此外,HNN模型仅代表一个局部的典型新皮层网络,而神经治疗手段及中枢神经系统疾病通常在多个脑区产生广泛分布的作用。尽管可通过改变外源性输入的时间和强度来近似模拟其他脑区的影响,但对这些上游或下游环路的直接实证表征往往有限,这限制了模型结果的解释能力。

参数简并性是所有生物物理神经模型中的一个基本挑战,因为多种参数配置可能产生相似的模型输出。在本方案中,使用基于模拟的推断(SBI)来解决这一问题,通过估计能够生成与实证数据一致的事件相关电位(ERP)波形的参数分布(图10)。该方法能够量化模型参数的不确定性,相较于单一的点估计,为机制性解释提供了更稳健的框架。然而,为了保证计算的可行性,SBI仅应用于与假设药物机制相对应的一组有限参数,而对未估计参数的假设可能会影响最终的网络动力学。通过使用诸如序列化神经后验估计等方法,可将推断扩展到更大的参数空间,该方法通过迭代优化参数估计,实现对高维参数分布的探索52(>10个维度)。除了概率推断外,引入独立的实验约束条件还可进一步降低参数不确定性,并提高模型预测的特异性。由于脑电图(EEG)信号主要反映皮层各层之间的协调活动,因此诸如侵入性层析电生理学等互补技术——包括神经元放电、局部场电位(LFP)和电流源密度(CSD)的测量——可为约束模型解和优化机制性假设提供重要信息。

成功应用本方案取决于对若干关键步骤的仔细执行。在确定ERP生物标志物并安装建模框架后(步骤1–2),主要要求是在工作流程的每个阶段均取得成功结果(图3)。在步骤3–5中,这包括选择和优化假设参数,这些参数可通过手动调节或优化,以实现模拟与实证的治疗前及治疗后ERP之间的良好拟合。如果无法获得满意的拟合效果,应探索并迭代测试其他参数。尽管鉴于HNN先前已成功复现ERP特征,反复失败的可能性较低,但若持续无法拟合,可能表明需要修改默认网络模型或引入更多生物物理细节。步骤6需要仔细配置SBI,包括合理选择参数范围、摘要统计量和训练参数,以确保对参数分布的准确估计。成功完成步骤6后,该方案将同时提供基于模型的预测结果及其相关的不确定性估计。步骤7对于验证这些预测至关重要,但具体的验证策略取决于可用的实验手段。潜在的验证方法包括:使用层析电生理记录评估特定皮层层和细胞的放电活动以及LFP/CSD信号7、进行分层解析的MEG/EEG测量、采用磁共振波谱或正电子发射断层扫描评估神经递质系统,以及利用扩散张量成像评估结构性连接(如丘脑-皮层通路)。

故障排除与自定义调整对于将本方案适配于不同数据集和实验情境至关重要,尤其是在第3–6步中对模型参数进行拟合以匹配实证数据的过程中。在参数优化(第4–5步)过程中,可能出现无法收敛至高相关性(相关系数 Corr > 0.95)的情况,此时可采取多种调整措施。这些措施包括修改优化器的超参数(例如,增加CMA-ES求解器的种群规模以提高鲁棒性,但会增加计算成本)、优化缩放与平滑参数(例如,测试5至60 ms范围内的平滑值),以及扩展外源性驱动参数的取值范围或引入额外的驱动信号,以更准确地捕捉波形特征。在某些情况下,优化后的模拟可能达到较高的整体相关性,但未能复现振幅较低的事件相关电位(ERP)成分(如P1成分);这一问题可通过设置更严格的损失阈值,或对特定时间窗赋予更高权重,以在优化过程中突出这些特征加以解决。对于基于模拟的推断(SBI,第6步),若后验预测检查(PPCs)失败,则表明模拟波形未能充分复现实证数据(补充图1)。此时应通过扩大参数范围或引入新的参数来修订先验参数分布,并可能需要增加训练数据集的规模。进一步改进还可通过调整摘要统计量或选择其他SBI架构实现。最后,若第7步的验证失败,则默认的HNN网络可能需要修改,以加入更多或不同的回路元件。HNN的模块化设计支持此类扩展:用户可通过图形用户界面(GUI)调整突触连接和细胞特性,也可通过Python接口实现更复杂的结构修改。例如,先前的研究已对默认模型进行改进,以纳入前额叶皮层中更精细的中间神经元连接结构46,从而产生新的可检验预测。HNN的开放框架有助于扩展模型的共享与复用,推动其在不同实验情境下的持续优化与验证。

披露

N.T. 和 S.R.J. 是一项关于本研究中所述神经回路模型参数推断方法的待审专利申请的共同发明人。其余作者声明无利益冲突。

致谢

本方案中用于生成所示结果的所有代码均可在以下网址获取:https://github.com/ntolley/hnn_jove。本研究得到了布朗大学生物医学创新与转化奖以及美国国立卫生研究院(NIH; https://www.nih.gov;资助编号 U24NS129945 和 P50MH109429)以及美国国家科学基金会(NSF;https://www.nsf.gov;资助编号 2424101)。资助方在研究设计、数据收集与分析、发表决定或稿件撰写过程中均未发挥任何作用。本研究使用的计算资源由布朗大学计算与可视化中心(CCV)通过NIH S10仪器资助项目S10OD036341(脑科学高性能计算集群)提供支持。

材料

本文使用的材料清单
姓名公司目录编号评论
Anaconda PythonAnaconda, Inc.N.A.Python 发行版;Python 版本 ≥3.9 且 <3.14
计算机工作站N.A.N.A.操作系统:Windows ≥10、Linux 或 macOS。最低推荐硬件配置:≥16 GB 内存,≥8 个 CPU 核心
EEGLABEEGLAB 开发团队N.A.基于 MATLAB 的可选工具箱,用于 EEG 预处理和事件相关电位(ERP)分析
FieldTrip拉德堡德大学唐德斯脑、认知与行为研究所N.A.基于 MATLAB 的可选工具箱,用于 EEG/MEG 分析
人类新皮层神经求解器(HNN-core)HNN 开发团队N.A.生物物理神经建模软件;本研究使用版本 ≥0.6.0
MATLABMathWorksN.A.运行 EEGLAB 和 FieldTrip(如使用)所必需
MNE-PythonMNE 开发团队N.A.用于 EEG 预处理和源定位
NumPyNumPy 开发团队N.A.用于数值计算和随机数生成
Pixi(包/环境管理器)Prefix.devN.A.用于管理相关代码仓库中的依赖项
PyTorchPyTorch 开发团队N.A.用于训练基于模拟的推断(SBI)神经网络和设置随机种子
SBI(基于模拟的推断)软件包SBI 开发团队N.A.用于参数推断和不确定性量化的 Python 软件包
Windows Linux 子系统(WSL2)Microsoft CorporationN.A.仅在基于 Windows 的安装中需要

参考文献

  1. Gribkoff VK, Kaczmarek LK. The need for new approaches in CNS drug discovery: Why drugs have failed, and what can be done to improve outcomes. Neuropharmacology. 2017;120:11-19.
  2. Loo SK, Lenartowicz A, Makeig S. Research review: Use of EEG biomarkers in child psychiatry research—current state and future directions. J Child Psychol Psychiatry. 2016;57(1):4-17.
  3. McLoughlin G, Makeig S, Tsuang MT. In search of biomarkers in psychiatry: EEG-based measures of brain function. Am J Med Genet B Neuropsychiatr Genet. 2014;165(2):111-121.
  4. Douglas RJ, Martin KAC. Neuronal circuits of the neocortex. Annu Rev Neurosci. 2004;27:419-451.
  5. Harris KD, Shepherd GMG. The neocortical circuit: themes and variations. Nat Neurosci. 2015;18(2):170-181.
  6. Murakami S, Okada Y. Contributions of principal neocortical neurons to magnetoencephalography and electroencephalography signals. J Physiol. 2006;575(3):925-936.
  7. Sherman MA, et al. Neural mechanisms of transient neocortical beta rhythms: Converging evidence from humans, computational modeling, monkeys, and mice. Proc Natl Acad Sci U S A. 2016;113(33):E4885-E4894.
  8. Shin H, et al. The rate of transient beta frequency events predicts behavior across tasks and species. eLife. 2017;6:e29086.
  9. De Pieri M, et al. Pharmaco-EEG of antipsychotic treatment response: A systematic review. Schizophrenia. 2023;9(1):85.
  10. Hyun J, Baik M, Kang U. Effects of psychotropic drugs on quantitative EEG among patients with schizophrenia-spectrum disorders. Clin Psychopharmacol Neurosci. 2011;9(2):78-85.
  11. Jobert M, et al. Guidelines for the recording and evaluation of pharmaco-EEG data in man: The International Pharmaco-EEG Society (IPEG). Neuropsychobiology. 2012;66(4):201-220.
  12. Mandema JW, Danhof M. Electroencephalogram effect measures and relationships between pharmacokinetics and pharmacodynamics of centrally acting drugs. Clin Pharmacokinet. 1992;23:191-215.
  13. Leiser SC, Dunlop J, Bowlby MR, Devilbiss DM. Aligning strategies for using EEG as a surrogate biomarker: A review of preclinical and clinical research. Biochem Pharmacol. 2011;81(12):1408-1421.
  14. Wilson FJ, Danjou P. Early decision-making in drug development: The potential role of pharmaco-EEG and pharmaco-sleep. Neuropsychobiology. 2016;72(3-4):188-194.
  15. Klumpp H, Shankman SA. Using event-related potentials and startle to evaluate time course in anxiety and depression. Biol Psychiatry Cogn Neurosci Neuroimaging. 2018;3(1):10-18.
  16. Proudfit GH, et al. Depression and event-related potentials: Emotional disengagement and reward insensitivity. Curr Opin Psychol. 2015;4:110-113.
  17. Luck SJ, et al. A roadmap for the development and validation of event-related potential biomarkers in schizophrenia research. Biol Psychiatry. 2011;70(1):28-34.
  18. Salisbury DF, Collins KC, McCarley RW. Reductions in the N1 and P2 auditory event-related potentials in first-hospitalized and chronic schizophrenia. Schizophr Bull. 2010;36(5):991-1000.
  19. Kang E, et al. Atypicality of the N170 event-related potential in autism spectrum disorder: A meta-analysis. Biol Psychiatry Cogn Neurosci Neuroimaging. 2018;3(8):657-666.
  20. Modi ME, Sahin M. Translational use of event-related potentials to assess circuit integrity in ASD. Nat Rev Neurol. 2017;13(3):160-170.
  21. Horvath A, et al. EEG and ERP biomarkers of Alzheimer’s disease: A critical review. Front Biosci (Landmark Ed). 2018;23:183-220.
  22. Malver LP, et al. Electroencephalography and analgesics. Br J Clin Pharmacol. 2014;77(1):72-95.
  23. Preskorn SH, et al. Normalizing effects of EVP-6124 on event-related potentials and cognition: A randomized trial in schizophrenia. J Psychiatr Pract. 2014;20(1):12-24.
  24. Schwertner A, et al. Effects of subanesthetic ketamine on visual and auditory event-related potentials in humans: A systematic review. Front Behav Neurosci. 2018;12:70.
  25. Visser S, et al. Dose-dependent EEG effects of zolpidem provide evidence for GABAA receptor subtype selectivity in vivo. J Pharmacol Exp Ther. 2003;304(3):1251-1257.
  26. Okoroafor F, et al. Neurophysiologic biomarkers of invasive neuromodulation therapy for epilepsy. Neuromodulation: Technol Neural Interface. 2026;29(3):360-375.
  27. Maki-Marttunen T, et al. Biophysical psychiatry—how computational neuroscience can help understand the complex mechanisms of mental disorders. Front Psychiatry. 2019;10:534.
  28. Murray JD, Demirtas M, Anticevic A. Biophysical modeling of large-scale brain dynamics and applications for computational psychiatry. Biol Psychiatry Cogn Neurosci Neuroimaging. 2018;3(9):777-787.
  29. Cagney DN, et al. The FDA NIH biomarkers, endpoints, and other tools (BEST) resource in neuro-oncology. Neuro Oncol. 2018;20(9):1162-1172.
  30. Cecchi M, et al. Validation of a suite of ERP and QEEG biomarkers in schizophrenia. Schizophr Res. 2023;254:178-189.
  31. Dura-Bernal S, et al. NetPyNE, a tool for data-driven multiscale modeling of brain circuits. eLife. 2019;8:e44494.
  32. Linden H, et al. LFPy: a tool for biophysical simulation of extracellular potentials. Front Neuroinform. 2014;7:41.
  33. Neymotin SA, et al. Human Neocortical Neurosolver (HNN), a new software tool for interpreting MEG/EEG data. eLife. 2020;9:e51214.
  34. Sanz Leon P, et al. The Virtual Brain: a simulator of primate brain network dynamics. Front Neuroinform. 2013;7:10.
  35. Hamalainen M, et al. Magnetoencephalography—theory, instrumentation, and applications. Rev Mod Phys. 1993;65(2):413-497.
  36. Ikeda H, Wang Y, Okada YC. Origins of the somatic N20 and high-frequency oscillations evoked by trigeminal stimulation in the piglets. Clin Neurophysiol. 2005;116(4):827-841.
  37. Okada YC, Wu J, Kyuhou S. Genesis of MEG signals in CNS structure. Electroencephalogr Clin Neurophysiol. 1997;103(4):474-485.
  38. Bush PC, Sejnowski TJ. Reduced compartmental models of pyramidal cells. J Neurosci Methods. 1993;46(2):159-166.
  39. Jones SR, et al. Neural correlates of tactile detection. J Neurosci. 2007;27(40):10751-10764.
  40. Jones SR, et al. Quantitative analysis of MEG mu rhythm. J Neurophysiol. 2009;102(6):3554-3572.
  41. Law RG, et al. Thalamocortical mechanisms regulating beta events. Cereb Cortex. 2022;32(4):668-688.
  42. Fernandez Pujol C, Blundon EG, Dykstra AR. Laminar specificity of the auditory perceptual awareness negativity: A biophysical modeling study. PLoS Comput Biol. 2023;19(6):e1011003.
  43. Kohl C, Parviainen T, Jones SR. Neural mechanisms underlying auditory evoked responses. Brain Topogr. 2022;35(1):19-35.
  44. Lankinen K, Ahveninen J, Jas M, Raij T, Ahlfors SP. Neuronal modeling of cross-sensory visual evoked magnetoencephalography responses in the auditory cortex. J Neurosci. 2024;44(17):e1119232024.
  45. Kaplan L, et al. Modeling cortical dynamics using HNN [poster presentation]. Presented at: Society for Neuroscience Annual Meeting; San Diego, CA, USA; 2025.
  46. Diesburg DA, Wessel JR, Jones SR. Biophysical modeling of ERP generation. J Neurosci. 2024;44(20).
  47. Bonaiuto JJ, et al. Laminar dynamics of beta bursts. Neuroimage. 2021;242:118479.
  48. Ferrante M, Blackwell KT, Migliore M, Ascoli GA. Computational models of neuronal biophysics. Curr Med Chem. 2008;15(24):2456-2471.
  49. Geerts H, et al. Quantitative systems pharmacology for neuroscience drug discovery. CPT Pharmacometrics Syst Pharmacol. 2020;9(1):5-20.
  50. Geerts H, et al. Computational neuroscience and systems pharmacology. J Pharmacokinet Pharmacodyn. 2024;51(5):563-573.
  51. Kappenman ES, Luck SJ. ERP components: brainwave recordings. Oxford Handbook ERP Components. 2012;1:3-30.
  52. Goncalves PJ, et al. Training neural density estimators. eLife. 2020;9:e56261.
  53. Papamakarios G, et al. Normalizing flows for probabilistic modeling. J Mach Learn Res. 2021;22(57):1-64.
  54. Pastore M, Calcagni A. Measuring distribution similarities. Front Psychol. 2019;10.
  55. Tolley N, et al. Estimating parameters in neural models with SBI. PLoS Comput Biol. 2024;20(2):e1011108.
  56. Gramfort A, et al. MEG and EEG analysis with MNE-Python. Front Neuroinform. 2013;7:267.
  57. Sliva DD, et al. Transcranial stimulation and EEG perception. Front Psychol. 2018;9:2117.
  58. Thorpe RV, et al. Distinct neocortical mechanisms underlie human SI responses to median nerve and laser evoked peripheral activation. bioRxiv. 2021; Available at: https://doi.org/10.1101/2021.10.11.463545.
  59. Delorme A, Makeig S. EEGLAB toolbox. J Neurosci Methods. 2004;134(1):9-21.
  60. Oostenveld R, et al. FieldTrip software. Comput Intell Neurosci. 2011;2011:156869.
  61. Puce A, Hämäläinen MS. A review of issues related to data acquisition and analysis in EEG/MEG studies. Brain Sci. 2017;7(6):58.
  62. Tejero-Cantero A, et al. sbi: toolkit for simulation-based inference. J Open Source Softw. 2020;5(52):2505.
  63. Lee S, Jones SR. Distinguishing mechanisms of gamma frequency oscillations in human current source signals using a computational model of a laminar neocortical network. Front Hum Neurosci. 2013;7:869.
  64. Szul MJ, et al. Beta burst waveform motifs. Prog Neurobiol. 2023;228:102490.
  65. Hashemi M, et al. Bayesian virtual epileptic patient. Neuroimage. 2020;217:116839.
  66. Billeh YN, et al. Systematic integration of structural and functional data into multi-scale models of mouse primary visual cortex. Neuron. 2020;106(3):388-403.
  67. Borges FS, Moreira JV, Takarabe LM, Lytton WW, Dura-Bernal S. Large-scale biophysically detailed model of somatosensory thalamocortical circuits in NetPyNE. Front Neuroinform. 2022;16:884245.
  68. Markram H, et al. Reconstruction of neocortical microcircuitry. Cell. 2015;163(2):456-492.
  69. Hagen E, Næss S, Ness TV, Einevoll GT. Multimodal modeling of neural network activity: computing LFP, ECoG, EEG, and MEG signals with LFPy 2.0. Front Neuroinform. 2018;12:92.
  70. Geerts H. Of mice and men in CNS drug discovery. CNS Drugs. 2009;23(11):915-926.
  71. Nestler EJ, Hyman SE. Animal models of neuropsychiatric disorders. Nat Neurosci. 2010;13(10):1161-1169.

重印与许可

标签

脑电图(EEG)生物物理建模EEG生物标志物神经环路活动事件相关电位听觉诱发电位电流源波形