Using the timing of the event in the ODE

Question

(This is a follow up question related to Scipy ODE time steps going backward)

I have a system of equations that I am trying to solve with scipy's solve_ivp. Here's a minimal working code:

import numpy as np
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 event(t,X):
    return X[0] + 20
event.terminal = False
event.direction = +1

def f(t, X):
    V = X[0]
    m = X[1]
    h = X[2]
    n = X[3]

    last_inputspike = inputspike[inputspike.searchsorted(t, side='right') - 1 ]
    last_t_event = -100 #Not sure what to put here

    g_syn_in = synapse(t, last_inputspike)
    g_syn_spike = synapse(t, last_t_event)
    syn = 0.5 * g_syn_in * (V - 0) + 0.2 * g_syn_spike * (V + 70)

    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)
inputspike = np.cumsum( np.random.exponential(beta, size=nbr_spike) ) + first_spike_date
inputspike = np.insert(inputspike, 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, events=event)
print(sol.message)

I want to detect when there is a spike (defined as V > 20), and have the timing of the spike affect the syn in the ODE via the changing g_syn_spike, in a similar way that the random input affects it.

In essence, I was wondering if it is possible and how I could go about accessing the last value of the sol.t_events at the given iteration of the solver?

It looks like your f already takes t and V as arguments and you assign g_syn_spike a value with those in scope. Can't you put an "if event(t, X)..." inside f? — GeorgeG, Sep 10 '18 at 17:28

score 1 · Answer 1 · answered Oct 30 '18 at 16:41

I have been looking for a way to simulate discrete events in continuous systems of differential equations as well. Modeling such discontinuities is not trivial, and two (recent) packages out there that can help you coop with this are:

assimulo - https://pypi.org/project/Assimulo/

simupy - https://pypi.org/project/simpy/

(and not simpy, this is only for discrete systems)

I hope this helps, in case you found a different solution already I'd like to hear what that is as well

Using the timing of the event in the ODE

1 Answers1