如何调整solve_ivp输出实现带阈值脉冲的整合发放神经元模型

你之前修改微分方程或输入电流的思路会触发解爆炸，原因是脉冲发放属于硬重置的不连续过程，无法纳入连续微分方程的求解范围，solve_ivp自带的事件检测机制就是专门处理这类场景的，具体实现方案如下：

实现原理 #

通过事件函数捕获膜电位从下往上穿越阈值的事件，自动终止当前段的连续求解，手动插入尖峰电位、重置膜电位后，再重启下一段的求解，完全保留solve_ivp的高精度、自适应步长特性。

完整可运行代码 #

import numpy as np
from scipy.integrate import solve_ivp

# 原有方波输入电流函数
def impulse_func(t, V):
    if 100 <= t <= 200:
        return 0.25
    return 0

# 原有被动膜电位微分方程
def passive_membrane(t, V, El, R, tau, I_app):
    dVdt = (-(V - El) + R*I_app(t, V))/tau     
    return dVdt

# 阈值穿越事件函数：返回值为0时触发事件
def threshold_cross_event(t, V, El, R, tau, I_app, V_thresh=-37):
    return V - V_thresh
# 配置事件属性：触发时终止求解，仅检测膜电位从下往上穿阈值的情况
threshold_cross_event.terminate = True
threshold_cross_event.direction = 1

# 可自由配置的模型参数
total_t_start = 0
total_t_end = 400
V_init = -60.
El = -60.
R = 100.
tau = 10.
V_thresh = -37. # 发放阈值可修改
spike_voltage = 60. # 尖峰电位可修改
spike_duration = 1. # 尖峰时长可修改
reset_V = -65. # 重置电位可修改

# 存储所有分段求解结果
all_t = []
all_V = []

current_t_start = total_t_start
current_V0 = V_init

while current_t_start < total_t_end:
    # 求解当前时间区间的膜电位变化
    sol = solve_ivp(
        fun=passive_membrane,
        t_span=[current_t_start, total_t_end],
        y0=np.array([current_V0]),
        args=(El, R, tau, impulse_func),
        events=threshold_cross_event,
        dense_output=True,
        rtol=1e-8,
        atol=1e-8
    )
    # 保存当前段的求解结果
    all_t.append(sol.t)
    all_V.append(sol.y[0])
    
    # 判断是否触发了阈值发放事件
    if len(sol.t_events[0]) > 0:
        event_t = sol.t_events[0][0]
        # 插入指定时长的尖峰电位
        spike_t = np.linspace(event_t, event_t + spike_duration, 2)
        spike_V = np.full_like(spike_t, spike_voltage)
        all_t.append(spike_t)
        all_V.append(spike_V)
        # 更新下一段求解的起始参数
        current_t_start = event_t + spike_duration
        current_V0 = reset_V
    else:
        # 全时间段求解完成，退出循环
        break

# 拼接所有结果得到完整的时间-电位序列
t_full = np.concatenate(all_t)
V_full = np.concatenate(all_V)

# 如果需要和你原代码一样得到0.1步长的固定采样结果，做插值即可
t_eval = np.arange(0, 400, 0.1)
V_eval = np.interp(t_eval, t_full, V_full)

特性说明 #

• 所有参数（阈值、尖峰幅度、尖峰时长、重置电位等）都支持自由修改，满足可变参数的需求
• 不会出现解爆炸问题，求解精度和效率远高于手动实现的RK2求解器
• 后续如果需要添加 refractory period（不应期）等逻辑，只需要修改重置后的起始时间和电位即可

内容的提问来源于stack exchange，提问作者NoviceConfusion