# Continuously tracking callback with boolean condition

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?

I solved this particular problem using rootfind = SciMLBase.RightRootFind when constructing the continuous callback, as stated in the documentation.

However, the question remains: How can I continuously check if some condition if true, without relying on a specific time point?

This would be extremely useful to implement e.g. a refractory period.