求解器实际上确实会来回移动,我认为也是因为可变的时间步长。但我认为困难来自于结果f(t, X)
不仅是一个函数t
and X
但是之前对此函数的调用,这不是一个好主意。
您的代码通过替换来工作:
inputspike(t)
g_syn = synapse(t, spiketrain[-1])
by:
last_spike_date = np.max( a[a<t] )
g_syn = synapse(t, last_spike_date)
并通过为“解决时间”设置“旧事件”a = np.insert(a, 0, -1e4)
。这是需要始终有一个last_spike_date
已定义(请参阅下面代码中的注释)。
这是您的代码的修改版本:
我修改了如何找到最后一个尖峰的时间(这次使用 Numpy 函数搜索排序 https://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html以便函数可以向量化)。我还修改了数组的方式a
被建造。这不是我的领域,所以也许我误解了意图。
I used solve_ivp https://scipy.github.io/devdocs/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp代替ode
但仍然有一个BDF求解器 https://scipy.github.io/devdocs/generated/scipy.integrate.BDF.html#scipy.integrate.BDF(但是,它的实现与ode
这是 Fortran 语言)。
import numpy as np # rather than scipy
import matplotlib.pylab as plt
from scipy.integrate import solve_ivp
def synapse(t, t0):
tau_1 = 5.3
tau_2 = 0.05
tau_rise = (tau_1 * tau_2) / (tau_1 - tau_2)
B = ((tau_2 / tau_1)**(tau_rise / tau_1) - (tau_2 / tau_1)**(tau_rise / tau_2)) ** -1
return B*(np.exp(-(t - t0) / tau_1) - np.exp(-(t - t0) / tau_2))
def alpha_m(v, vt):
return -0.32*(v - vt -13)/(np.exp(-1*(v-vt-13)/4)-1)
def beta_m(v, vt):
return 0.28 * (v - vt - 40) / (np.exp((v- vt - 40) / 5) - 1)
def alpha_h(v, vt):
return 0.128 * np.exp(-1 * (v - vt - 17) / 18)
def beta_h(v, vt):
return 4 / (np.exp(-1 * (v - vt - 40) / 5) + 1)
def alpha_n(v, vt):
return -0.032*(v - vt - 15)/(np.exp(-1*(v-vt-15)/5) - 1)
def beta_n(v, vt):
return 0.5* np.exp(-1*(v-vt-10)/40)
def f(t, X):
V = X[0]
m = X[1]
h = X[2]
n = X[3]
# Find the largest value in `a` before t:
last_spike_date = a[ a.searchsorted(t, side='right') - 1 ]
# Another simpler way to write this is:
# last_spike_date = np.max( a[a<t] )
# but didn't work with an array for t
g_syn = synapse(t, last_spike_date)
syn = 0.5 * g_syn * (V - 0)
dVdt = - 50*m**3*h*(V-60) - 10*n**4*(V+100) - syn - 0.1*(V + 70)
dmdt = alpha_m(V, -45)*(1-m) - beta_m(V, -45)*m
dhdt = alpha_h(V, -45)*(1-h) - beta_h(V, -45)*h
dndt = alpha_n(V, -45)*(1-n) - beta_n(V, -45)*n
return [dVdt, dmdt, dhdt, dndt]
# Define the spike events:
nbr_spike = 20
beta = 100
first_spike_date = 500
np.random.seed(0)
a = np.cumsum( np.random.exponential(beta, size=nbr_spike) ) + first_spike_date
a = np.insert(a, 0, -1e4) # set a very old spike at t=-1e4
# it is a hack in order to set a t0 for t<first_spike_date (model settle time)
# so that `synapse(t, t0)` can be called regardless of t
# synapse(t, -1e4) = 0 for t>0
# Solve:
t_start = 0.0
t_end = 2000
X_start = [-70, 0, 1,0]
sol = solve_ivp(f, [t_start, t_end], X_start, method='BDF', max_step=1, vectorized=True)
print(sol.message)
# Graph
V, m, h, n = sol.y
plt.plot(sol.t, V);
plt.xlabel('time'); plt.ylabel('V');
这使:
注意:有一个事件参数solve_ivp
这可能有用。