小鼠神经元记录中的皮层慢动力学分析方案

原文：Shymkiv Y, Yuste R. Protocol for analyzing slow cortical dynamics in mouse neuronal recordings. STAR Protocols, 2026-08-27, 7(3):104804。

DOI: 10.1016/j.xpro.2026.104804 ｜ PMCID: PMC13544393

单位：美国哥伦比亚大学（Yuste 实验室）· CC BY 4.0 开放获取

一句话概要

用三套互补的计算方法——逐 bin 刺激解码器、逐试次变异性分析和内在时间尺度（INT）估计——从群体神经元记录中定量测量"神经网络对最近刺激的记忆能维持多久"。方法以清醒小鼠听皮层双光子钙成像数据和 RNN 合成数据为示范，但只要求格式规整的活动时间序列，不限定脑区与物种。完整使用细节参考原始研究：Shymkiv et al., Neuron 113 (2025) 847–857.e8。

亮点（Highlights）

• 用逐 bin 解码器测量皮层网络中的慢速多神经元活动；
• 刻画慢网络动力学对逐试次变异性的影响；
• 提供分析网络与神经元时间尺度的内在时间尺度（INT）指标。

背景与原理

大脑在短时间尺度上保留信息，是新奇检测（novelty detection）、适应（adaptation）和时间归一化等计算的基础；这些时间信息加工过程受损与精神分裂症等疾病相关。感觉皮层中的"慢动力学"指网络通过延长反应来短期存储感觉输入的记忆，可能表现为单个神经元的持续放电，也可能表现为神经元依次激活形成的状态空间轨迹。

本方案用三种视角测量同一个量——网络活动保留近期刺激记忆的时间尺度——并互相印证：

方法	测什么	回答什么问题
对角线（diagonal）逐 bin 解码器	刺激身份保持可解码的时长	刺激信息在网络里存多久
全（cross-temporal）逐 bin 解码器	跨时间 bin 的解码泛化能力	群体编码是"静态维持"还是"沿轨迹演化"
逐试次变异性分析	相同刺激各试次响应的相关性	前序刺激诱导的持续状态对当前响应的影响（随 ISI 增大而消退）
内在时间尺度（INT）	活动偏离基线的持续时间	单神经元与网络各自的记忆时标

比较单神经元 INT 与网络 INT，还能判断慢网络动力学来自个体神经元的持续放电，还是来自瞬时激活的接力序列。

创新点

1. 用逐 bin 跨时间解码器追踪群体刺激信息的演化，识别演化轨迹中的持续状态；
2. 以逐试次相关随 ISI 增大而恢复的过程来估计记忆时标；
3. 把 INT 从单神经元/合并活动扩展到网络水平——以群体活动向量偏离基线的欧氏距离为标量轨迹，再估计其自相关衰减。这样能检出"单神经元时标都短、但群体靠接力保持长时标"的情形。RNN 提供参照区间：记忆无关任务（频率识别）的网络给出下界，记忆依赖任务（oddball）的网络给出上界，皮层网络落在两者之间。

硬件与软件准备

硬件

• 钙成像分析与 RNN control 数据集：≥16 GB 内存、4 核 CPU 即可；
• RNN oddball 数据集整批载入内存，建议 ≥64 GB；
• 原文验证环境：Windows 10，Intel Core i9-14900K，128 GB RAM。

软件（Timing: <15 min）

1. 安装 Miniconda（anaconda.com）；2. 建立专用环境并安装依赖：

conda create -n slowdyn python=3.12.12
conda activate slowdyn
conda install numpy scipy matplotlib scikit-learn jupyter notebook h5py

3. 下载分析管线：github.com/shymkivy/slow_dynamics_protocol，解压到本地目录。管线含全部三种分析的 demo 脚本（.ipynb / .py）。

关键资源表

类别	资源	来源/标识
数据	钙成像与 RNN 数据	Dryad: doi.org/10.5061/dryad.xsj3tx9q6
软件	Python	v3.12.12（RRID:SCR_008394）
软件	Miniconda	Anaconda（RRID:SCR_018317）
软件	NumPy / SciPy / Scikit-learn	v2.3.5 / v1.16.3 / v1.7.1
软件	Jupyter Notebook / Matplotlib / h5py	v7.5.0 / v3.10.7 / v3.15.1
代码	慢动力学分析管线	GitHub: shymkivy/slow_dynamics_protocol

数据准备（Timing: <30 min）

本方案面向同时记录的多个神经元的时间序列数据（放电率、突触电流等代理信号均可）。示范数据两套：

1. 钙成像：清醒小鼠听皮层（A1/A2/AAF/UF）在 oddball 范式下被动听取声音的双光子记录。oddball 刺激天然依赖对近期刺激的短期记忆——新异刺激放大响应（新奇检测），重复刺激抑制响应（适应）；
2. RNN：三类循环网络——oddball 训练（新奇检测）、频率识别训练（对照任务）、未训练网络；训练后用同分布的新输入测试，测试期单元活动按钙成像数据同等处理（25–50 个单元，时间步长 0.05 s）。

方法与数据的对应关系

分析	数据源	输入范式	分析对象	所需文件
逐 bin 解码器	钙成像	对照（10 个频率）	频率身份	Oddball（"mismatch"），M1–M10
逐 bin 解码器	RNN	对照（50 个频率）	频率身份	*_cont_data.npy + *_params.npy
逐 bin 解码器	RNN	Oddball（50 频率两两组合）	频率身份	*_ob_data.npy + *_params.npy
逐试次变异性	钙成像	变 ISI 对照 run	同一刺激的逐试次相关	变 ISI（"echo"）：M226、M4264、M4265、M4266、M4371、M4372
INT	钙成像	全连续记录（对照和/或 oddball）	神经元与网络时标	Oddball（"mismatch"），M1–M10
INT	RNN	全连续记录（对照或 oddball）	神经元与网络时标	*_cont_data.npy（或 *_ob_data.npy）+ *_params.npy

样本量建议（重要）：这些都是群体水平分析，应尽量最大化同时记录的神经元数与试次数。所需样本量随任务表征的活动维度增长：简单任务（如 10 个感觉刺激）至少需要约 10 个神经元；实操底线是 ≥10 个神经元、每个刺激至少重复 10 次（解码与变异性分析）。INT 分析可用自发活动的连续记录替代；若用刺激试次计算，事件之间最好留出衰减间隔。

分步实验方案

第 1 部分：数据加载与预处理（Timing: <5 min）

目标：把原始数据转成后续所有分析通用的放电率矩阵、试次身份、刺激起始时间。

步骤 1｜启动环境

conda activate slowdyn
jupyter notebook

步骤 2｜导入库

pipeline_dir = 'your analysis pipeline directory goes here'  # edit this

# ---- importing functions ----
import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform

# importing slow dynamics pipeline
sys.path.append(pipeline_dir + '/functions')
import sd_utils as sd
import sd_decoder as dec

步骤 3｜固定随机种子（保证可复现）

seed = 0  # integer seed for reproducible results, None for random

步骤 4｜加载并预处理钙成像数据

data_dir = 'your CaIm data directory goes here'  # edit this

# ---- loading CaIm datasets ----
# firing rates, trial types, and stimuli times
data_ob = sd.load_caim_data_mat(
  data_dir,                      # data directory
  ext_list=['mat'],              # data extension
  tags=['ammn', '_processed_data'],  # data tags
  num_files=20,                  # limit number of files to load
  deconvolution='oasis',         # deconvolution oasis or smoothdfdt
  smooth_std_duration=0.1)       # in seconds

• 去卷积算法二选一：oasis（CaImAn）或 smoothdfdt（平滑整流一阶导数）；
• 放电率可用高斯核平滑（smooth_std_duration，默认标准差 0.1 s）；
• 返回值 data_ob 是字典列表，每条记录含 firing_rates（神经元 × 帧矩阵）、stim_times（刺激起始帧）、trial_types（试次身份）。外部数据需改写加载步骤，同样组织成字典列表。

步骤 5｜提取帧率（外部数据可手动指定，单位 Hz）

frame_rate = 1000/np.mean(sd.get_values(data_ob, 'volume_period'))  # in Hz

步骤 6｜平行加载 RNN 数据（对照输入测试）

data_dir_rnn = 'your RNN data directory goes here'  # edit this
fname_rnn = 'RNN_test_data_2024_5_24_9h_42m2'

data_rnn = sd.load_rnn_test(
  data_dir_rnn,
  fname_rnn + '_cont_data.npy',
  fname_rnn + '_params.npy',
  max_net_load = 5,
  flatten_runs = True,
  max_trial_types = 10,
  max_trials = 400,
  limit_network_types=['ob trained', 'freq trained'],
  seed=seed)

frame_rate_rnn = 1000/np.mean(sd.get_values(data_rnn, 'volume_period'))
training_type = np.array(sd.get_values(data_rnn, 'training'))

RNN 数据结构与钙成像一致，可直接复用同一套分析函数；frame_rate_rnn 由 0.05 s 时间步换算。

步骤 7｜画示例放电率栅格图（检查数据质量，sd.plot_fig_raster）。

步骤 8｜组装三维刺激触发响应矩阵（神经元 × 帧 × 试次）

# compute trial window frames (assuming all datasets have similar frame rate)
trial_frames, plot_t = sd.get_frames(
  trial_win=[-1,3],     # seconds relative to stim onset (pre, post)
  frame_rate=frame_rate)

# extract stimulus-triggered average matrix (neurons, frames, trials)
stim_trig_resp_all = []
for n_fl in range(len(data_ob)):
  stim_trig_resp = sd.get_stim_trig_resp(
    data_ob[n_fl]['firing_rates'],
    data_ob[n_fl]['stim_times'],
    trial_frames=trial_frames)
  stim_trig_resp_all.append(stim_trig_resp)

• 钙成像试次窗口取 (−1, +3) s（1 s 基线 + 3 s 响应）；
• RNN 慢动力学更持久，窗口取 (−1, +5) s（get_frames + get_stim_trig_resp 同上，得到 stim_trig_resp_all_rnn）。

第 2 部分：逐 bin 解码器（Timing: <30 min）

用 demo 脚本 slow_dynamics_decoder_comb_data_rnn.ipynb 运行。

原理：对试期内每个时间 bin 分别训练并测试一个解码器。对角线（diag）结构回答"刺激信息维持多久"；全（full）结构让某一 bin 训练的解码器在所有 bin 上测试，回答"群体编码是静态还是沿轨迹演化"。

步骤 9｜定义对角线解码器。每 bin 用线性 SVM 做 one-vs-all 分类，5 折交叉验证；同时跑标签打乱（shuffle）对照：

def f_diag_decoder(stim_trig_resp, trial_types, seed=None):
  X_all = [stim_trig_resp,   # data
           stim_trig_resp]   # data repeat for label-shuffled control
  Y_all = [trial_types,
           dec.shuffle_trials(trial_types, seed=seed)]
  dec_data = dec.run_binwise_dec(
    X_all,
    Y_all,
    train_test_method='diag',   # options: full, diag
    pca_var_frac=1,             # reduce data with pca
    num_cv=5,                   # cross validation
    normalize=False,
    add_noise_sigma=1e-5,       # for stability
    get_train_coeffs=True,
    seed=seed)
  return dec_data

步骤 10｜遍历全部钙成像数据集训练（trial_types_use = trial_types<10 选择前 10 类试次）。

步骤 11｜同样训练所有 RNN 数据集（得到 dec_data_all_rnn）。

步骤 12｜画对角线解码结果（sd.plot_fig_diag_decoder）。

预期：刺激只呈现 0–0.5 s（ISI 0.5 s），但解码显著高于打乱对照的时间远超刺激期——刺激相关活动在刺激结束后仍持续；oddball 训练的 RNN 轨迹比对照训练的更持久，说明 oddball 计算需要时间上的信息存储。

步骤 13｜定义全（cross-temporal）解码器（train_test_method='full'，其余参数同上，log=True）。

步骤 14｜对示例钙成像数据集运行全解码器。

步骤 15｜画全解码热图（sd.plot_fig_full_decoder_caim，clim=[0, 0.7]）。解码性能随训练/测试 bin 变化，即群体活动沿轨迹演化。

步骤 16｜对示例 RNN 数据集运行全解码器（分别取 training_type == 'ob trained' 和 'freq trained' 的数据集）。

步骤 17｜在解码空间画图（sd.plot_fig_full_decoder_RNN）：无论 oddball 还是对照训练的 RNN，刺激响应期群体状态都沿轨迹移动，而非停留在单一持续状态。

第 3 部分：逐试次变异性分析（Timing: <5 min）

用 demo 脚本 slow_dynamics_isi_correlation.ipynb 运行。

原理：标准 oddball 及其对照用短 ISI（0.5 s）。若网络保留前序刺激的持续状态，当前响应会被"污染"，逐试次变异性升高；ISI 拉长后该效应应减弱。本节分析四种 ISI（0.5 / 1 / 2 / 4 s）的扩展对照数据集。

步骤 18｜加载变 ISI 对照数据

data_dir = 'your variable ISI data directory goes here'  # edit this

data_echo = sd.load_caim_data_mat(
  data_dir,
  ext_list=['mat'],
  tags=['cont', '_processed_data'],
  num_files=None,
  deconvolution='smoothdfdt',
  smooth_std_duration=0.15)      # in seconds

frame_rate = 1000/np.mean(sd.get_values(data_echo, 'volume_period'))  # Hz
isi_list = np.array(sd.get_values(data_echo, 'isi'))

注意本节平滑参数为 0.15 s（前两节为 0.1 s）。

步骤 19｜计算刺激响应统计、筛出响应细胞

trial_frames_tuning, plot_t_tuning = sd.get_frames(
  trial_win=[-1, 2],   # seconds relative to stim onset (pre, post)
  frame_rate=frame_rate)

resp_cells_all = []
for n_fl in range(len(data_echo)):
  stim_trig_resp = sd.get_stim_trig_resp(
    data_echo[n_fl]['firing_rates'],
    data_echo[n_fl]['stim_times'],
    trial_frames=trial_frames_tuning)
  resp_cells = sd.compute_tuning(
    stim_trig_resp,
    data_echo[n_fl]['trial_types'],
    np.arange(1,11),
    plot_t_tuning,
    num_samp=2000,
    z_thresh = 3,
    sig_resp_win = [0, 1.2],
    seed=seed)
  resp_cells_all.append(resp_cells)

响应细胞判定用自助法峰值检测：对每个细胞、每种刺激，取刺激前后 −1 至 +2 s 窗口的试次平均响应，定位峰值并在峰值附近约 0.25 s 窗口内取均值作为峰幅值；再随机抽 2000 组与该条件等试次数（不区分刺激身份）的试次集，构造峰幅值零分布。实测峰值超过 z = 3 的零分布阈值、且落在刺激后 0–1.2 s 响应窗内，即判为对该刺激有响应。

步骤 20｜计算同刺激试次间的两两相关

trial_frames, plot_t = sd.get_frames(
  trial_win = [-0.05, .95],   # seconds relative to stim
  frame_rate = frame_rate)

corr_vals_all = []
stim_trig_resp_all = []
for n_fl in range(len(data_echo)):
  stim_trig_resp = sd.get_stim_trig_resp(
    data_echo[n_fl]['firing_rates'],
    data_echo[n_fl]['stim_times'],
    trial_frames=trial_frames)
  stim_trig_resp_all.append(stim_trig_resp)
  corr_vals = sd.compute_correlation(
    stim_trig_resp,
    data_echo[n_fl]['trial_types'],
    np.arange(1,11),               # trials types to analyze (1-10)
    resp_cells = resp_cells_all[n_fl],
    metric='correlation',          # cosine, correlation
    min_resp_cells=5,              # minimum responsive cells
    seed=seed)
  corr_vals_all.append(corr_vals)
corr_vals = np.vstack(corr_vals_all)

• 相关窗口取 (−0.05, +0.95) s（匹配 0.5 s 刺激 + 0.5 s ISI），窗口不得与相邻试次明显重叠，避免引入邻近刺激的活动；
• 每个频率计算试次-试次相关矩阵并取均值，仅纳入响应细胞数 ≥5（min_resp_cells）的组。

步骤 21｜画相关值随 ISI 的变化（sd.plot_fig_isi_corr_trials），各较长 ISI 与 0.5 s 比较用双侧 Mann–Whitney U 检验（sd.stat_compare）。

预期：ISI 越短，逐试次相关性越低（变异性越高）；ISI 增大，相关性恢复——这就是网络记忆时标的功能学估计。

步骤 22｜画不同 ISI 的逐试次余弦相似度矩阵（compute_correlation_mat，metric='cosine'，subtract_mean=True；示例选 M4372 鼠频率 5）。颜色越亮 = 相似度越高 = 变异性越低。

第 4 部分：内在时间尺度（INT）分析（Timing: <5 min）

用 demo 脚本 slow_dynamics_net_int.ipynb 运行。

步骤 23｜加载数据

• 钙成像：沿用第 1 部分步骤 2–4 的代码；INT 分析不按试次组织数据，因为不需要刺激时刻信息，连续记录即可；
• RNN：重新加载，改用以下参数：

data_rnn = sd.load_rnn_test(
  data_dir_rnn,
  fname_rnn + '_ob_data.npy',
  fname_rnn + '_params.npy',
  max_net_load = 5,
  flatten_runs = False,          # 3D matrix: runs × neurons × time
  cut_zero_trials = True,
  num_initial_trials_skip = 20,  # discard first 20 trials
  limit_network_types=['ob trained', 'freq trained', 'untrained'],
  seed=seed)

跳过前 20 个试次，避免输入序列开头大规模网络状态跃迁的影响。

步骤 24｜计算网络与单神经元 INT

# ---- Compute intrinsic timescales for CaIm data ----
tau_ob_net_all = []
tau_ob_cell_all = []
for n_dset in range(len(data_ob)):
  tau_net, tau_cell = sd.get_network_intrinsic_timescales(
    data_ob[n_dset]['firing_rates'],
    frame_rate)
  tau_ob_net_all.append(tau_net)
  tau_ob_cell_all.append(tau_cell)

# ---- Compute intrinsic timescales for RNN data ----
tau_rnn_net_all = []
tau_rnn_cell_all = []
for n_rnn in range(len(data_rnn)):
  tau_net, tau_cell = sd.get_network_intrinsic_timescales(
    data_rnn[n_rnn]['firing_rates'],
    frame_rate_rnn)
  tau_rnn_net_all.append(tau_net)
  tau_rnn_cell_all.append(tau_cell)

方法要点：

• 单神经元 INT：对每条放电率迹求自相关，估计衰减时间常数。相比指数拟合，本方案采用更稳健的半高法——从零滞后向外找自相关首次降到峰值一半的时间点，对密集激活信号更稳定；
• 网络 INT：先构造网络活动标量迹，再套用同样的单神经元流程。基线向量 b = 各神经元在全记录的均值活动 b = ⟨r(t)⟩_t；网络迹为 d(t) = ‖r(t) − b‖₂（瞬时群体向量与基线的欧氏距离），INT 取 d(t) 自相关的衰减常数；
• 基线取全记录全局均值；本分析用去卷积钙信号，基线漂移很小，其他数据可能需先去趋势。各神经元放电率先归一化，避免个别神经元主导 d(t)。因为 INT 基于自相关，恒定的基线偏移不影响结果；
• 网络 INT 刻意不区分"单一持续状态"与"一串不同状态"——是哪种由逐 bin 解码器回答。

步骤 25｜画 INT 分布（sd.plot_fig_tau_networks_comb，do_log=True；小提琴图 + 中位数/四分位距箱线，组间用双侧 Mann–Whitney U 检验）。

预期：听皮层网络时间常数长于神经元时间常数——早期感觉皮层的群体沿轨迹演化（瞬时神经元的接力把群体维持在偏离基线的状态）。网络时间常数还可以与 RNN 对照：频率识别（无记忆需求）给出下界、oddball 训练（优化为保留刺激历史）给出上界，皮层网络恰落在两者之间，符合通用皮层网络的中间记忆需求。

预期结果

三种分析只需格式规整的活动数据，不限定脑区与物种（原文应用于小鼠初级与次级听皮层）。解码与变异性分析聚焦刺激触发响应；INT 适用于任意活动迹。整体目标是找到网络活动延长维持的证据（可能支撑短期记忆和时间计算，如新奇检测）：对角线解码超越打乱对照的时间窗、全解码热图呈现的时间特异性轨迹、逐试次相关随 ISI 的恢复曲线、以及网络 INT 显著长于神经元 INT 的分布。

局限性

• 解码器只能给出慢动力学的间接读出，不能揭示具体的神经环路或机制；解码成功也不意味着该信息真的被网络读取或使用；
• INT 依赖自相关衰减估计，对数据质量敏感：密集且高度相关的活动会产生不稳定的时间常数估计（本方案的半高法即为缓解此问题）。网络 INT 以"距基线的欧氏距离"压缩群体活动，把与基线等距的不同状态视为等价——这是有意的设计，状态的"身份/方向"交给解码器刻画；
• 双光子钙成像受钙信号时间分辨率与帧率（约 60 Hz）限制，指示剂信号是经过时间滤波的间接测量，去卷积还会引入平滑；关键结论最好用更高时间分辨率的记录复现；
• 生物学验证目前限于小鼠听皮层；推广到其他皮层区或物种时需按时间分辨率与活动统计差异调整参数；
• 生物记录与 RNN 的对比需谨慎解读：RNN 的架构、训练目标与输入统计未必完全复刻生物环路，应理解为"共享刺激结构下的功能比较"，而非 RNN 离散时间步与生物秒的一一对应。

Troubleshooting

问题	相关步骤	可能原因与解决办法
解码器运行慢或不收敛	9–14	① 先标准化/归一化数据（每神经元基线归零、最大值归一，或 z-score）；② 数据集过大——用 PCA 降维，保留约 70% 总方差（pca_var_frac=0.7）；③ 试次类型太多导致 one-vs-all 过复杂——可简化为两类问题
逐试次变异性看不出结构或差异	20–21	试次数不足。最大化试次数，直到估计维度或解释方差达到平台；可用 PCA 检查随试次增加解释方差的变化，出现平台说明采样充分
INT 估计出过长/过短、无生理意义的时间常数	24	INT 对放电率迹的类型敏感，误差主要来自自相关衰减常数的估计。理想迹应大部分静息、激活稀疏（因此常用自发活动计算）。对密集激活数据，先画自相关迹检查衰减常数是否可可靠估计，必要时调整算法

资源与代码

• 数据（钙成像 + RNN）：Dryad，doi.org/10.5061/dryad.xsj3tx9q6
• 分析管线：GitHub，github.com/shymkivy/slow_dynamics_protocol
• 联系人：Yuriy Shymkiv（ys2605@columbia.edu）