Hello,
I have the following system of equations:
using DifferentialEquations
using Plots
u_rest = -50.
t_ref = 2.
u_reset = -75.
u_thr = -55.
J = 1.
last_spikes = [1., 1.] .* (- t_ref + eps())
function decay(du, u, p, t)
du[1] = (u_rest-u[1]) * 0.5
du[2] = (u_rest-u[2]) * 0.5
end
function spike_condition(u, t, integrator)
prod(u .- u_thr)
end
function refractory_condition(u, t, integrator)
prod([(integrator.t - last_spike) for last_spike in last_spikes])
end
function spike_affect!(integrator)
above_threshold = integrator.u .>= u_thr
integrator.u[.! above_threshold] .= integrator.u[.! above_threshold] .+ J
integrator.u[above_threshold] .= u_reset
last_spikes[above_threshold] .= integrator.t
end
function refractory_affect!(integrator)
in_refractory = [(integrator.t - last_spike <= t_ref) for last_spike in last_spikes]
integrator.u[in_refractory] .= u_reset
print(integrator.u)
end
spike_cb = ContinuousCallback(spike_condition, spike_affect!);
cb_refractory = ContinuousCallback(refractory_condition, refractory_affect!);
vector_cbs = CallbackSet(spike_cb, cb_refractory);
u0 = [-70.,-61.];
tspan = (0.0, 10.0);
prob = ODEProblem(decay, u0, tspan);
sol = solve(prob, Tsit5(), callback = spike_cb);
plot(sol, size = (1500,500),
xtickfontsize = 15,
ytickfontsize =15,
legendfontsize = 15,
title = "When one curves reaches threshold, the other is offset",
titlefontsize = 20,
)
hline!([-55.], linestyle=:dash, label = "threshold")
which creates the above plot in the figure.
This code has a problem for me: the components of u, when they reach the threshold, trigger an event (i.e. reset of component to 0 and small increment on other component). However, since they interact, whenever a component of u is above the threshold (which should not happen), no event is triggered.
Is there a way to use a continuous callback with a condition that checks if the components of u are greater (or smaller) than 0?
I have tried solving this with a DiscreteCallback, but as stated in the documentation the event is not triggered because it needs to be triggered at a specific time. Here is a minimal working example:
using DifferentialEquations
using Plots
u_rest = -50.
u_thr = -55.
u_reset = -75.
t_refractory = 2.
last_spike = -(t_refractory + eps())
#callback for spiking event
function decay(du, u, p, t)
du[1] = (u_rest - u[1])
end
function spike(u, t, integrator)
u_thr < u[1]
end
function spike_affect!(integrator)
integrator.u[1] = u_reset
last_spike = integrator.t
end
spike_callback = DiscreteCallback(spike, spike_affect!)
u0 = [-70.];
tspan = (0.0, 5.0);
prob = ODEProblem(decay, u0, tspan);
sol = solve(prob, Tsit5(), callback = spike_callback);
p = plot(sol, size = (1500,500),
xtickfontsize = 15,
ytickfontsize =15,
legendfontsize = 15,
title = "Dynamical variable is allowed to cross threshold",
titlefontsize = 20,
fmt = :png
)
hline!([-55.], linestyle=:dash, label = "threshold")
but here the problem is that the dynamical variable can cross the threshold (see below plot):
So my question is: how can I get a DiscreteCallback that checks at every simulation step? Or is there an ‘inequality’ ContinuousCallback?