Disable instability checks in OrdinaryDiffEq

Specific Domains Modelling & Simulations
1

Hi,

I have a small network of 3 neurons and I’m using OrdinaryDiffEq to solve the ode system describing their interconnected dynamics. Since this problem is stiff, I’m using an adaptive solver that can handle it and this works just fine.

Now, say I don’t care too much about good accuracy because all I’m interested in is the behaviour which can be roughly eyeballed from the activity of the neurons. I want to use a fixed timestep solver, hoping that this might (?) speed up the solution.
However using Euler, even with a small timestep like with dt = 0.0001, gives an instability warning.

So my question is: How to disable the instability checking option if I want to grind through this regardless?

Here’s the working code

using NeuronBuilder, OrdinaryDiffEq

const Area = 0.0628 
const Cm = 10.0 
conv = Cm / Area
synaptic_conv = 1e-3 / Area^2 

defaults = Dict(Voltage => BasicVoltageDynamics(),
    Calcium => Prinz.CalciumDynamics(200.0, 0.05, 14.96 * 0.0628),
    Reversal{Calcium} => Prinz.CaReversalDynamics())

somatic_parameters = Dict(
    Reversal{Sodium} => 50.0,
    Reversal{Potassium} => -80.0,
    Reversal{Leak} => -50.0,
    Reversal{Proton} => -20.0,
    Voltage => -60.0,
    Calcium => 0.05,
    Reversal{Calcium} => 0.0)

function AB_Neuron(num_connections)
    AB2_ch = [
        Prinz.Na(100.0 * conv),
        Prinz.CaS(6.0 * conv),
        Prinz.CaT(2.5 * conv),
        Prinz.H(0.01 * conv),
        Prinz.Ka(50.0 * conv),
        Prinz.KCa(5.0 * conv),
        Prinz.Kdr(100.0 * conv),
        Prinz.leak(0.0 * conv)
    ]
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, AB2_ch, :AB)
    return b(incoming_connections=num_connections)
end

function PY_Neuron(num_connections)
    PY1_ch = [Prinz.Na(100.0 * conv),
        Prinz.CaS(2.0 * conv),
        Prinz.CaT(2.4 * conv),
        Prinz.H(0.05 * conv),
        Prinz.Ka(50.0 * conv),
        Prinz.KCa(0.0 * conv),
        Prinz.Kdr(125.0 * conv),
        Prinz.leak(0.01 * conv)
    ]
    somatic_parameters[Voltage] = -55.0
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, PY1_ch, :PY)
    return b(incoming_connections=num_connections)
end

function LP_Neuron(num_connections)
    LP4_ch = [Prinz.Na(100.0 * conv),
        Prinz.CaS(4.0 * conv),
        Prinz.CaT(0.0 * conv),
        Prinz.H(0.05 * conv),
        Prinz.Ka(20.0 * conv),
        Prinz.KCa(0.0 * conv),
        Prinz.Kdr(25.0 * conv),
        Prinz.leak(0.03 * conv)
    ]
    somatic_parameters[Voltage] = -65.0
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, LP4_ch, :LP)
    return b(incoming_connections=num_connections)
end

AB = AB_Neuron(1)
LP = LP_Neuron(3)
PY = PY_Neuron(3)
neurons = [AB, PY, LP]

ABLP_chol = directed_synapse(AB, LP, Chol(30.0 * synaptic_conv))
ABPY_chol = directed_synapse(AB, PY, Chol(3.0 * synaptic_conv))
ABLP_glut = directed_synapse(AB, LP, Glut(30.0 * synaptic_conv))
ABPY_glut = directed_synapse(AB, PY, Glut(10.0 * synaptic_conv))

LPAB_glut = directed_synapse(LP, AB, Glut(30.0 * synaptic_conv))
LPPY_glut = directed_synapse(LP, PY, Glut(1.0 * synaptic_conv))

PYLP_glut = directed_synapse(PY, LP, Glut(30.0 * synaptic_conv))

connections = [LPAB_glut, ABPY_chol, ABPY_glut, LPPY_glut, ABLP_chol, ABLP_glut, PYLP_glut]

stg = add_all_connections(neurons, connections)
prob = ODEProblem(stg, [], (0.0, 10000.0), []; jac=true)

@time sol = solve(prob, AutoTsit5(Rosenbrock23()))

The output of this block of code with the AutoTsit5(Rosenbrock23()) algorithm is 1.683045 seconds (253.82 k allocations: 78.925 MiB, 1.97% gc time).

But changing the last line to @time sol = solve(prob, Euler(), dt=0.0001) gives

┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/OHiiA/src/integrator_interface.jl:351

Any help is greatly appreciated!

Cheers.

2

unstable_check = (args...) -> false.

https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous

I’m not sure about that.

3

Ah great, I’d found the same issue in another discourse but the link there was broken.
Thanks!

Yeah a fixed step solver with very small step isn’t speeding up the solution, just wanted to check if the instability checks allowed me to try with a simple euler.

Will need to give more thought to the other solving options :slight_smile:

1 Like