使用NEMOS进行钙成像数据分析教程

前言

本教程将介绍如何使用NEMOS工具包对钙成像数据进行建模分析。钙成像技术是现代神经科学研究中的重要手段，能够记录大量神经元的活动情况。我们将通过一个实际案例，展示如何建立Gamma-GLM模型来分析小鼠后下托(head-direction)神经元的钙信号。

数据准备

数据来源

我们使用的数据集来自自由活动小鼠的1光子钙成像记录，采样频率为30Hz，使用GCaMP6f钙指示剂。记录区域为后下托(postsubiculum)，该区域包含头部方向细胞(head-direction cells)。

数据加载与初步处理

首先，我们需要加载NWB格式的数据文件：

import nemos as nmo
import pynapple as nap

# 加载数据文件
path = nmo.fetch.fetch_data("A0670-221213.nwb")
data = nap.load_file(path)

数据探索

数据中包含多个字段，我们主要关注RoiResponseSeries字段，它包含了钙信号时间序列：

transients = data['RoiResponseSeries']
print(transients.time_support)  # 显示记录时长

我们可以可视化部分神经元的原始钙信号：

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].plot(transients[:, 4].get(0,250))
ax[0].set_title("神经元4")
ax[1].plot(transients[:, 35].get(0,250))
ax[1].set_title("神经元35")
plt.tight_layout()

数据预处理

时间降采样

为了减少计算复杂度，我们将数据从30Hz降采样到10Hz：

Y = transients.bin_average(0.1, ep)  # 100ms时间窗平均

头部方向调谐曲线

我们可以计算神经元的头部方向调谐曲线：

tcurves = nap.compute_1d_tuning_curves_continuous(transients, data['ry'], 120)

模型构建

基函数选择

我们使用两种基函数的组合：

1. 循环B样条基(Cyclic BSpline)用于编码头部方向
2. 对数间隔的升余弦基(Raised Cosine)用于神经元间的耦合关系

heading_basis = nmo.basis.CyclicBSplineEval(n_basis_funcs=12, label="heading")
coupling_basis = nmo.basis.RaisedCosineLogConv(3, window_size=10, label="coupling")
basis = heading_basis + coupling_basis

Gamma-GLM模型

由于钙信号是非负连续值，我们选择Gamma分布作为观测模型，并使用softplus作为链接函数：

model = nmo.glm.GLM(
    regularizer="Ridge",
    regularizer_strength=0.02,
    observation_model=nmo.observation_models.GammaObservations(
        inverse_link_function=jax.nn.softplus
    )
)

模型训练与评估

数据划分

我们将数据分为训练集和测试集：

train_ep = nap.IntervalSet(start=X.time_support.start, end=X.time_support.get_intervals_center().t)
test_ep = X.time_support.set_diff(train_ep)

Xtrain = X.restrict(train_ep)
Ytrain = Y.restrict(train_ep)
Xtest = X.restrict(test_ep)
Ytest = Y.restrict(test_ep)

模型训练

选择目标神经元进行训练：

neu = 4
selected_neurons = jnp.hstack((jnp.arange(0, neu), jnp.arange(neu+1, Y.shape[1])))
model.fit(Xtrain, Ytrain[:, neu])

与传统线性回归比较

我们与scikit-learn的线性回归模型进行对比：

from sklearn.linear_model import LinearRegression

mdl = LinearRegression()
valid = ~jnp.isnan(Xtrain.d.sum(axis=1))
mdl.fit(Xtrain[valid], Ytrain[valid, neu])

结果分析

预测效果可视化

比较三种预测结果：

yp = model.predict(Xtest)
ylreg = nap.Tsd(t=yp.t, d=mdl.predict(Xtest), time_support=yp.time_support)

plt.figure()
plt.plot(Ytest[:,neu].restrict(ep_to_plot), "r", label="真实值")
plt.plot(yp.restrict(ep_to_plot), "k", label="Gamma-GLM")
plt.plot(ylreg.restrict(ep_to_plot), "g", label="线性回归")
plt.legend()
plt.show()

调谐曲线比较

计算并比较不同模型的调谐曲线：

real_tcurves = nap.compute_1d_tuning_curves_continuous(transients, data['ry'], 120, ep=test_ep)
gamma_tcurves = nap.compute_1d_tuning_curves_continuous(yp, data['ry'], 120, ep=test_ep)
linreg_tcurves = nap.compute_1d_tuning_curves_continuous(ylreg, data['ry'], 120, ep=test_ep)

plt.figure()
plt.plot(real_tcurves[neu], "r", label="真实值")
plt.plot(gamma_tcurves, "k", label="Gamma-GLM")
plt.plot(linreg_tcurves, "g", label="线性回归")
plt.legend()
plt.show()

技术要点总结

1. 数据特性处理：钙信号是非负连续值，Gamma分布比高斯分布更合适
2. 基函数选择：循环基函数适合周期性变量(如头部方向)，升余弦基适合描述神经元间耦合
3. 模型比较：Gamma-GLM能保证预测值非负，更符合生理实际
4. 正则化：使用Ridge正则化防止过拟合

应用前景与局限

虽然Gamma-GLM在钙成像数据分析中显示出潜力，但这一方法仍处于探索阶段。未来研究可以关注：

• 不同链接函数的比较
• 更复杂的耦合关系建模
• 与其他去卷积方法的结合

本教程展示了NEMOS在钙成像数据分析中的应用，为神经科学研究提供了新的建模思路。