Izhikevichch Model: feedback events

[mongibellili]

The implementation of the model by DifferentialEquations.jl can be found here.

QuantizedSystemSolver

using QuantizedSystemSolver
function izh!(du, u, p, t)
    a, b, c, d, I = p
    du[1] = 0.04 * u[1]^2 + 5.0 * u[1] + 140.0 - u[2] + I
    du[2] = a * (b * u[1] - u[2])
    if u[1] >= 30.0
        u[1] = c
        u[2] += d
    end
    if t >= 50.0
        p[5] += 10.0
    end
end
p = [0.02, 0.2, -50.0, 2.0, 0.0]
u0 = [-60.0, p[2] * -60.0]
tspan = (0.0, 300.0)
prob = ODEProblem(izh!, u0, tspan, p)
sol = solve(prob,liqss2(),reltol=1e-3,abstol=1e-3 );
p1=plot(sol, idxs = 1)
p2=plot(sol, idxs = 2)
@btime solve($prob,liqss2(),reltol=1e-3,abstol=1e-3 )

Jacobian at (u₁, u₂) = (-65, -13)
[[-0.2 -1. ]
[ 0.004 -0.02 ]]

Eigenvalues
λ ≈ { -0.1740, -0.0460 }

The total simulation steps: 1665
The number of events: 27
598.200 μs (79 allocations: 62.29 KiB)

removing the time event:

sol.stats =
The total simulation steps: 1780
The number of events: 36
The number of state steps per Var: [1458, 286]
618.700 μs (64 allocations: 69.01 KiB)

removing all the events gives:

sol.stats =
The total simulation steps: 397
The number of state steps per Var: [335, 62]
120.800 μs (62 allocations: 60.21 KiB)

[mongibellili]

classic solvers (julia)

ContinuousCallback leads to
#= ┌ Warning: dt(5.684341886080802e-14) <= dtmin(5.684341886080802e-14) at t=3.3714630621970634, ...
so DiscreteCallback is used, which produces inaccurate results near events as shown here:

algName = :ABDF2
sol.stats = SciMLBase.DEStats(6120, 0, 641, 4896, 21, 4704, 4, 1488, 744, 442, 0.0)
52.091 ms (217733 allocations: 17.07 MiB)
algName = :QNDF2
sol.stats = SciMLBase.DEStats(9599, 0, 1981, 8246, 180, 6875, 127, 2292, 1146, 840, 0.0)
86.566 ms (677242 allocations: 51.70 MiB)
algName = :Trapezoid
sol.stats = SciMLBase.DEStats(5818, 0, 606, 4759, 58, 4273, 26, 1422, 711, 259, 0.0)
46.629 ms (211526 allocations: 16.67 MiB)
algName = :TRBDF2
sol.stats = SciMLBase.DEStats(5650, 0, 569, 5997, 115, 4505, 49, 754, 377, 115, 0.0)
11.305 ms (37846 allocations: 3.19 MiB)
algName = :Rosenbrock23
sol.stats = SciMLBase.DEStats(2188, 0, 615, 1845, 464, 0, 0, 928, 464, 151, 0.0)
27.626 ms (184214 allocations: 13.56 MiB)
algName = :Rodas5
sol.stats = SciMLBase.DEStats(3804, 0, 406, 3248, 277, 0, 0, 554, 277, 129, 0.0)
16.433 ms (96493 allocations: 8.14 MiB)

adding these benchmarks at 1e-4,1e-5 gives:

@btime solve($prob,Heun(),reltol=1e-4,abstol=1e-5 );#106.868 ms (1051583 allocations: 85.56 MiB)
@btime solve($prob,BS3(),reltol=1e-4,abstol=1e-5 );#18.501 ms (240360 allocations: 20.12 MiB)
@btime solve($prob,Rosenbrock23(),reltol=1e-4,abstol=1e-5 ); #96.980 ms (594070 allocations: 45.88 MiB)
@btime solve($prob,Rodas5(),reltol=1e-4,abstol=1e-5 );#20.499 ms (111503 allocations: 9.18 MiB)
@btime solve($prob,Trapezoid(),reltol=1e-4,abstol=1e-5 )#71.321 ms (270055 allocations: 20.67 MiB)

@btime solve($prob,QBDF2(),reltol=1e-4,abstol=1e-5)#199.999 ms (1586020 allocations: 120.08 MiB)
@btime solve($prob,ABDF2(),reltol=1e-4,abstol=1e-5)# 92.802 ms (375751 allocations: 31.36 MiB)
@btime solve($prob,QNDF2(),reltol=1e-4,abstol=1e-5)#  387.156 ms (2603240 allocations: 194.68 MiB)
@btime solve($prob,TRBDF2(),reltol=1e-4,abstol=1e-5)#23.306 ms (64671 allocations: 5.84 MiB)
@btime solve($prob,KenCarp5(),reltol=1e-4,abstol=1e-5)#92.724 ms (83729 allocations: 6.81 MiB)

@btime solve($prob,nmliqss2(),reltol=1e-4,abstol=1e-5); #2.376 ms (74 allocations: 290.77 KiB)

Trapezoid sol.stats = SciMLBase.DEStats(12519, 0, 1076, 10374, 24, 10151, 3, 2868, 1434, 657, 0.0)
Notice that the number of iterations to detect discontinuities or to linearSolve is high.

removing all the events:

algName = :TRBDF2
sol.stats = SciMLBase.DEStats(141, 0, 29, 162, 3, 125, 0, 0, 23, 4, 0.0)
50.200 μs (290 allocations: 21.08 KiB)
algName = :BS3
sol.stats = SciMLBase.DEStats(252, 0, 0, 0, 0, 0, 0, 0, 83, 0, 0.0)
25.900 μs (382 allocations: 33.69 KiB)
algName = :Tsit5

removing the preset time event:

TRBDF2= SciMLBase.DEStats(6000, 0, 581, 6377, 113, 4838, 57, 4501, 400, 121, 0.0)
1.724 ms (3187 allocations: 236.22 KiB)
algName = :Rodas5
sol.stats = SciMLBase.DEStats(3914, 0, 417, 3336, 288, 0, 0, 288, 288, 129, 0.0)
710.200 μs (2826 allocations: 223.42 KiB)

using tstops instead of state events:

sol.stats = SciMLBase.DEStats(5660, 0, 628, 6033, 115, 4512, 50, 398, 398, 121, 0.0)
1.565 ms (3618 allocations: 242.98 KiB)

only discrete callbacks which are not accurate

• Tried all kwargs
• Help asked in julia discourse continuousCallBack

[mongibellili]

Matlab:

function Izh
global fEvalCount
tstart =  0.0; tfinal = 300.0;
y0 = [-60.0,  -12.0];  % initial condition 
refine = 10;
options = odeset('Events',@events,...
   'Refine',refine,'RelTol',1.0e-3,'AbsTol',1.0e-3,...
      'Stats', 'off');  
tout = tstart;  % tout and yout store the output time and state vector
yout = y0;%.';      % for each segment of simulations.
teout = [];     % time at event.
yeout = [];     % state at event.
ieout = [];     % event number at event.
totalSteps = 0;
totalFailed = 0;
totalFuncEvals = 0;
totalEvents = 0;
fEvalCount = 0;
while tstart < tfinal
    [t,y,te,ye,ie] = ode23s(@f,[tstart tfinal],y0,options);
    nt = length(t);
    totalSteps = totalSteps + length(t)-1;        % accepted steps
    % failed steps are not accessible programmatically for ode45, only printed
    totalEvents = totalEvents + numel(ie);       % number of events triggered
    % Append intermediate points if any
    if nt > 1
        tout = [tout; t(2:nt)];
        yout = [yout; y(2:nt,:)];  % Now dimensions match
    else
        % If only one point, append it
        tout = [tout; t(end)];
        yout = [yout; y(end,:)];
    end

    teout = [teout; te];
    yeout = [yeout; ye];
    ieout = [ieout; ie];

    if ~isempty(ie)
        y0 = y(end,:);  % copy last state
        y0(1) = -50.0;
        y0(2) = y0(2) + 2.0;
    end

    tstart = t(end);
end
totalFuncEvals = fEvalCount;
fprintf('Total accepted steps: %d\n', totalSteps);
fprintf('Total function evaluations: %d\n', totalFuncEvals);
fprintf('Total events detected: %d\n', totalEvents);

%plot(teout,yeout(:,1),'ro');
plot(tout,yout);
xlabel('time');
ylabel('Vars');
title('Izh');
odeplot([],[],'done');

% -----------------------------------------------------
function dydt = f(~,y)
    fEvalCount = fEvalCount + 1;
   dydt = [0.04 * y(1)^2 + 5.0 * y(1) + 140.0 - y(2) + 10;0.02 * (0.2 * y(1) - y(2))];   
end 
% -----------------------------------------------------
function [value,isterminal,direction] = events(t,y)
    value = y(1)-30.0; 
    isterminal =1 ; 
    direction = 1 ; 
end 
end % main function

Results:

Rel,abs=-3,-3
ode45
Total accepted steps: 2310
Total function evaluations: 2068
Total events detected: 27
Ode23
Total accepted steps: 5810
Total function evaluations: 2845
Total events detected: 27
Ode23t
Total accepted steps: 9000
Total function evaluations: 2497
Total events detected: 28
Ode23s
Total accepted steps: 4990
Total function evaluations: 3019
Total events detected: 27

[mongibellili]

Chatgpt talk: (not checked)

Why QSS (Quantized State Systems) Performs Well for the Izhikevich Model

1. QSS vs. Traditional ODE Solvers

Traditional solvers (e.g., Heun(), Rodas5()) work by:

• Approximating the solution using small fixed or adaptive time steps.
• Solving the differential equations at every step, even if the state changes slowly.
• Performing frequent, expensive floating-point computations and memory allocations.

In contrast, QSS solvers like nmliqss2() work differently:

• They discretize the state (not just time) and update variables only when needed.
• Instead of global time-stepping, QSS uses event-driven updates.
• Each state variable is updated independently, leading to efficient handling of sparse dynamics.

2. Why QSS Works Well for the Izhikevich Model

The Izhikevich model exhibits sharp voltage spikes due to neuron firing. After ( t = 50 ), the injected current ( I ) increases, leading to frequent action potentials (spikes).

🔹 Sparse & Asynchronous Updates

• Traditional solvers update all variables at every step.
• QSS updates only variables that significantly change, reducing unnecessary computation.
• During the resting phase (before ( t = 50 )), updates are rare.
• After ( t = 50 ), QSS efficiently detects and processes only the sharp changes.

🔹 Event-Driven Execution

• The model has two types of events:
  1. A single external event at ( t = 50 ) (increase in ( I )).
  2. Multiple threshold-crossing events when ( u_1 \geq 30 ) (neuron spikes).
• Traditional solvers handle these by reducing step size (causing slowdowns).
• QSS treats them as discrete events, updating only the affected variables efficiently.

🔹 Adaptive Quantization

• In smooth regions (resting potential), QSS takes large time steps.
• During rapid fluctuations (spiking), QSS refines updates dynamically.
• This reduces unnecessary calculations compared to fixed-step solvers.

3. Computational Efficiency

The benchmarks confirm this:

• nmliqss2() solves the problem in 2.376 ms, while explicit solvers like BS3() take 18.5 ms (8× slower).
• Traditional solvers allocate megabytes of memory, whereas nmliqss2() uses only 290 KiB.

This efficiency is due to:

1. Sparse, asynchronous updates instead of global time-stepping.
2. Event-driven execution, efficiently handling discontinuities.
3. Adaptive quantization, refining resolution only when needed.

Since the Izhikevich model has frequent sharp changes (spikes), QSS naturally aligns with the problem structure, leading to huge computational savings.

[mongibellili]

24-11-2025:
ft 300, rel,abs= -3,-3

No preset event

(typeof(alg)).name.name = :Tsit5
sol.stats = SciMLBase.DEStats(1716, 0, 0, 0, 0, 0, 0, 2499, 215, 66, 0.0)
385.000 μs (2234 allocations: 204.23 KiB)
(typeof(alg)).name.name = :BS3
sol.stats = SciMLBase.DEStats(2331, 0, 0, 0, 0, 0, 0, 6394, 576, 191, 0.0)
830.400 μs (2519 allocations: 209.66 KiB)
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(6594, 0, 0, 0, 0, 0, 0, 38833, 3234, 48, 0.0)
2.699 ms (13750 allocations: 1.15 MiB)

(typeof(alg)).name.name = :Rodas5
sol.stats = SciMLBase.DEStats(3924, 0, 418, 3344, 289, 0, 0, 3272, 289, 129, 0.0)
1.184 ms (3104 allocations: 242.88 KiB)
(typeof(alg)).name.name = :ABDF2
sol.stats = SciMLBase.DEStats(7122, 0, 751, 5702, 20, 5515, 1, 9327, 846, 536, 0.0)
2.555 ms (5147 allocations: 417.66 KiB)
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(10831, 0, 2307, 9267, 214, 7696, 148, 14571, 1320, 977, 0.0)
5.370 ms (27758 allocations: 2.06 MiB)
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(6298, 0, 733, 5139, 72, 4512, 31, 8405, 759, 297, 0.0)
2.319 ms (4891 allocations: 399.69 KiB)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(6000, 0, 581, 6377, 113, 4838, 57, 4501, 400, 121, 0.0)
1.984 ms (3189 allocations: 236.59 KiB)
(typeof(alg)).name.name = :Rosenbrock23
sol.stats = SciMLBase.DEStats(2298, 0, 645, 1935, 489, 0, 0, 5450, 489, 156, 0.0)
1.158 ms (4071 allocations: 298.98 KiB)

complete izh:

(typeof(alg)).name.name = :Tsit5
sol.stats = SciMLBase.DEStats(1698, 0, 0, 0, 0, 0, 0, 2645, 210, 68, 0.0)
381.200 μs (2270 allocations: 207.33 KiB)
(typeof(alg)).name.name = :BS3
sol.stats = SciMLBase.DEStats(2250, 0, 0, 0, 0, 0, 0, 6735, 557, 183, 0.0)
860.900 μs (2531 allocations: 210.70 KiB)
(typeof(alg)).name.name = :Rodas5
sol.stats = SciMLBase.DEStats(3826, 0, 408, 3264, 280, 0, 0, 3450, 280, 128, 0.0)
1.181 ms (3127 allocations: 244.78 KiB)
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(9303, 0, 1976, 7953, 175, 6623, 116, 13819, 1148, 849, 0.0)
5.096 ms (27814 allocations: 2.06 MiB)
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(6026, 0, 677, 4898, 55, 4417, 24, 8890, 737, 305, 0.0)
2.310 ms (4870 allocations: 399.14 KiB)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(5482, 0, 569, 5822, 112, 4395, 53, 4659, 381, 102, 0.0)
1.954 ms (3211 allocations: 238.14 KiB)
(typeof(alg)).name.name = :Rosenbrock23
sol.stats = SciMLBase.DEStats(2200, 0, 618, 1854, 467, 0, 0, 5674, 467, 151, 0.0)
1.157 ms (4076 allocations: 299.55 KiB)

liqss2() The total simulation steps: 1677
The number of events: 39
The number of state steps per Var: [1393, 245]
606.900 μs (66 allocations: 87.80 KiB)
qss2() The total simulation steps: 2616
The number of events: 37
The number of state steps per Var: [2275, 304]
870.100 μs (56 allocations: 87.40 KiB)

no events (I=0.0)

(typeof(alg)).name.name = :BS3
sol.stats = SciMLBase.DEStats(252, 0, 0, 0, 0, 0, 0, 0, 83, 0, 0.0)
24.800 μs (380 allocations: 33.62 KiB)
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(129, 0, 26, 101, 4, 71, 0, 0, 21, 10, 0.0)
49.400 μs (506 allocations: 38.34 KiB)
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(143, 0, 20, 109, 2, 99, 0, 0, 29, 0, 0.0)
42.100 μs (295 allocations: 21.77 KiB)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(141, 0, 29, 162, 3, 125, 0, 0, 23, 4, 0.0)
50.700 μs (290 allocations: 21.08 KiB)
(typeof(alg)).name.name = :Rosenbrock23
sol.stats = SciMLBase.DEStats(99, 0, 24, 72, 24, 0, 0, 0, 24, 0, 0.0)
27.400 μs (280 allocations: 20.97 KiB)
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(245, 0, 0, 0, 0, 0, 0, 0, 121, 0, 0.0)
33.800 μs (529 allocations: 44.41 KiB)

liqss2() The total simulation steps: 40
The number of state steps per Var: [25, 15]
14.800 μs (62 allocations: 19.71 KiB)
qss2()
The total simulation steps: 260
The number of state steps per Var: [246, 14]
73.200 μs (52 allocations: 19.30 KiB)

dense=false for DE.jl
complete izh:

(typeof(alg)).name.name = :BS3
sol.stats = SciMLBase.DEStats(2250, 0, 0, 0, 0, 0, 0, 6735, 557, 183, 0.0)
480.000 μs (709 allocations: 70.73 KiB)
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(9303, 0, 1976, 7953, 175, 6623, 116, 13819, 1148, 849, 0.0)
5.305 ms (23749 allocations: 1.75 MiB)
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(6026, 0, 677, 4898, 55, 4417, 24, 8890, 737, 305, 0.0)
1.841 ms (2491 allocations: 204.02 KiB)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(5482, 0, 569, 5822, 112, 4395, 53, 4659, 381, 102, 0.0)
1.702 ms (1928 allocations: 137.23 KiB)
(typeof(alg)).name.name = :Rosenbrock23
sol.stats = SciMLBase.DEStats(2200, 0, 618, 1854, 467, 0, 0, 5674, 467, 151, 0.0)
909.900 μs (2518 allocations: 178.83 KiB)
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(6594, 0, 0, 0, 0, 0, 0, 38833, 3234, 48, 0.0)
2.742 ms (3516 allocations: 382.58 KiB)

no events
(typeof(alg)).name.name = :BS3
sol.stats = SciMLBase.DEStats(252, 0, 0, 0, 0, 0, 0, 0, 83, 0, 0.0)
18.600 μs (130 allocations: 14.02 KiB)
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(129, 0, 26, 101, 4, 71, 0, 0, 21, 10, 0.0)
48.300 μs (443 allocations: 33.75 KiB)
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(143, 0, 20, 109, 2, 99, 0, 0, 29, 0, 0.0)
39.300 μs (208 allocations: 15.42 KiB)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(141, 0, 29, 162, 3, 125, 0, 0, 23, 4, 0.0)
48.700 μs (221 allocations: 16.05 KiB)
(typeof(alg)).name.name = :Rosenbrock23
sol.stats = SciMLBase.DEStats(99, 0, 24, 72, 24, 0, 0, 0, 24, 0, 0.0)
25.500 μs (208 allocations: 15.72 KiB)
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(245, 0, 0, 0, 0, 0, 0, 0, 121, 0, 0.0)
24.100 μs (165 allocations: 16.48 KiB)

[mongibellili]

no preset event
abstol, reltol = 1e-4,1e-3
Tsist5
sol.stats = SciMLBase.DEStats(3450, 0, 0, 0, 0, 0, 0, 4926, 439, 131, 0.0)
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(6023, 0, 593, 6428, 110, 4913, 40, 4539, 404, 141, 0.0)
1.766 ms (3220 allocations: 238.27 KiB)
IZH_err = 0.017101896553298235
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(11127, 0, 2378, 9492, 197, 8036, 131, 15537, 1408, 1023, 0.0)
4.545 ms (29059 allocations: 2.15 MiB)
IZH_err = 0.009846227254257517
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(6601, 0, 704, 5437, 64, 4872, 25, 8519, 769, 300, 0.0)
1.982 ms (4865 allocations: 398.72 KiB)
IZH_err = 0.023151439286182587
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(7904, 0, 0, 0, 0, 0, 0, 41530, 3773, 164, 0.0)
2.979 ms (15310 allocations: 1.27 MiB)
IZH_err = 0.0015998401417987708

[mongibellili]

complete izh
tf300, rel,abs=-3,-4
(typeof(alg)).name.name = :TRBDF2
sol.stats = SciMLBase.DEStats(5921, 0, 599, 6298, 112, 4789, 49, 4773, 390, 130, 0.0)
1.745 ms (3232 allocations: 239.20 KiB)
IZH_err = 0.010307278052920876
(typeof(alg)).name.name = :QNDF2
sol.stats = SciMLBase.DEStats(10337, 0, 2228, 8806, 172, 7564, 125, 15995, 1329, 961, 0.0)
4.586 ms (28797 allocations: 2.13 MiB)
IZH_err = 0.005327517546096974
(typeof(alg)).name.name = :Trapezoid
sol.stats = SciMLBase.DEStats(6096, 0, 655, 4981, 55, 4510, 23, 8884, 736, 293, 0.0)
1.993 ms (4913 allocations: 402.02 KiB)
IZH_err = 0.015370492666921958
(typeof(alg)).name.name = :Heun
sol.stats = SciMLBase.DEStats(7560, 0, 0, 0, 0, 0, 0, 43310, 3607, 158, 0.0)
3.012 ms (15330 allocations: 1.27 MiB)
IZH_err = 0.0011634768638637014

liqss2
The total simulation steps: 1274
The number of events: 37
The number of state steps per Var: [1039, 198]
456.500 μs (70 allocations: 65.93 KiB)
IZH_err = 0.03180163947407906

qss2()
The total simulation steps: 1996
The number of events: 41
The number of state steps per Var: [1708, 247]
657.600 μs (60 allocations: 65.52 KiB)
IZH_err = 0.032831876122108514