我有一个问题是关于使用ODE求解器在Python中实现Leaky integrate and fire模型。
我有一个使用Euler方法和for循环的实现：

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
tau_m = 10  # membrane time constant (ms)
v_rest = -70  # resting membrane potential (mV)
v_thresh = -55  # spike threshold (mV)
v_reset = -80  # reset potential (mV)
i_e = 3  # input current (nA)
dt = 0.1  # time step (ms)
t_max = 100  # simulation time (ms)

# Initialize variables
v = v_rest  # membrane potential (mV)

# Simulate LIF model
for t in np.arange(0, t_max, dt):
    dv_dt = (-(v - v_rest) + i_e) / tau_m  # membrane potential differential equation
    v += dv_dt * dt  # update membrane potential
    if v >= v_thresh:  # spike condition
        v = v_reset  # reset membrane potential

现在的问题是，我有一个更大的模型使用solve_ivp（）实现，我需要将LIF模型集成到其中。
定义膜电位的微分项（dv_dt）是可行的，我可以这样做：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define LIF model
def lif(t, y, tau_m, v_rest, v_thresh, v_reset, i_e):
    dv_dt = (-(y - v_rest) + i_e) / tau_m  # membrane potential differential equation
    return dv_dt

# Define parameters
tau_m = 10  # membrane time constant (ms)
v_rest = -70  # resting membrane potential (mV)
v_thresh = -55  # spike threshold (mV)
v_reset = -80  # reset potential (mV)
i_e = 3  # input current (nA)
t_span = (0, 100)  # simulation time span (ms)
y0 = v_rest  # initial membrane potential (mV)

# Simulate LIF model
sol = solve_ivp(lambda t, y: lif(t, y, tau_m, v_rest, v_thresh, v_reset, i_e), t_span, [y0])

# Extract results
t = sol.t
v = sol.y[0]

但我不确定如何实现去极化后膜电位的硬变化。

v = v_reset if v >= v_thresh  # spike condition

因为当您调用求解器solve_ivp（）时，函数的唯一输出是微分项dv_dt。我想过要做：

dv_dt = (v_reset - y) if y >= v_thresh else (-(y - v_rest) + i_e) / tau_m

因此，我在下一次迭代中通过dv_dt实现对y的改变，但尽管这降低了膜电位，但它不会将其重置为v_reset。
任何想法将不胜感激！