Solvers do not model uncertainty

[michaeldeistler]

Hi there!

Thanks for creating this toolbox!

I just tried the solvers for the first time and I'm running into some problems. My ODE is a Hodgkin-Huxley model (a model in neuroscience). Some specs (full code below):

• the model has 4 states
• it's non-linear
• it's non-homogenous

When simulating it with the EK0 solver and low tolerance, the solution is the same as if I used a simple Euler() solver. The EK0 solution is the light blue curve in the figure below. However, when I increase the tolerance in the EK0 solver, the solution changes, but all samples are the same. These are the other five traces in the figure (you can see only the dark green one because they are almost, but not bitwise, identical).

• Have you tried running the Hodgkin-Huxley model?
• Should I expect different / better results when using https://github.com/probabilistic-numerics/probnum?
• Have you observed this failure to model uncertainty?

Thanks! Michael

[michaeldeistler]

Here's my ODE:

function hh_pospischil(du, u, p, t)
    E = u[1]
    m = u[2]
    h = u[3]
    n = u[4]

    V_T = -63.0

    g_Na, g_Kd, g_leak = p

    alpha_act = -0.32 * (E - V_T - 13) / (exp(-(E - V_T - 13) / 4) - 1)
    beta_act = 0.28 * (E - V_T - 40) / (exp((E - V_T - 40) / 5) - 1)
    du[2] = (alpha_act * (1.0 - m)) - (beta_act * m)

    alpha_inact = 0.128 * exp(-(E - V_T - 17) / 18)
    beta_inact = 4 / (1 + exp(-(E - V_T - 40) / 5))
    du[3] = (alpha_inact * (1.0 - h)) - (beta_inact * h)

    alpha_kal = -0.032 * (E - V_T - 15) / (exp(-(E - V_T - 15) / 5) - 1)
    beta_kal = 0.5 * exp(-(E - V_T - 10) / 40)
    du[4] = (alpha_kal * (1.0 - n)) - (beta_kal * n)

    area = 20000.0
    I_Na = -(E - 50) * g_Na * area * (m^3) * h
    I_K = -(E - (-90)) * g_Kd * area * n^4
    I_leak = -(E - (-65)) * g_leak * area

    I_ext = 1.0 * area

    du[1] = (I_leak + I_K + I_Na +I_ext) / (1 * area)

end

This is my setup:

p = [70.0, 10.0, 0.05];

u0 = [-70, 0.0, 1.0, 0.1];
tspan = [0., 30];

I generature the accurate solution with:

prob = ODEProblem(hh_pospischil, u0, tspan, p);
sol1 = solve(prob, EK0(), abstol=1e-3, reltol=1e-4);
u_ = [sol1.u[i][1] for i = 1:size(sol1.u)[1]];
plot(sol1.t, u_, label="`sol` with accurate tol.")

and the inaccurate solution with

prob = ODEProblem(hh_pospischil, u0, tspan, p);
sol = solve(prob, EK0(), abstol=1e-1, reltol=1e-0);
samples = ProbNumDiffEq.sample(sol, 5);
for i in 1:5
    plot!(sol.t, samples[:, 1, i], label="Sample with inaccurate tol.")
end

[nathanaelbosch]

Hi Michael, thanks for opening this issue! It's very helpful to see what problems people are working on and to get some feedback on these solvers.

Thanks also for providing all the code. I looked into the example and might have a few pointers on how to get more useful output, but let me first answer your questions.

• Have you observed this failure to model uncertainty?

Unfortunately, yes. Model calibration is not perfect, and especially for high-tolerance regimes the model can become overconfident. For low tolerances however my current experience rather suggests underconfident uncertainties.

• Should I expect different / better results when using https://github.com/probabilistic-numerics/probnum?

I don't think you should expect different results with ProbNum. It does provide some additional functionality, such as more priors (IOUP and Matern) and an unscented Kalman filtering-based solver, but I also think that these should not influence the qualitative behaviour of the calibration. I still encourage you to have a look if you're interested in this topic, the documentation is much more polished and the tutorials might be helpful if you want to learn more about these solvers.

• Have you tried running the Hodgkin-Huxley model?

I have not previously considered this model, but I suspect that it might be one of the more challenging settings for the EK solvers (compared to sampling-based solvers), especially with high tolerances, due to the fact that they are only able to model Gaussian uncertainties.

But, I was able to obtain a more useful output: The EK0 with default diffusionmodel=:dynamic models the same uncertainty for each dimension, but the scales in the Hodgkin-Huxley example vary a lot between the dimensions. With EK0(diffusionmodel=:dynamicMV), the solver models the uncertainties in each dimension independently from each-other, resulting in larger uncertainties for the first component:

Now, the plot clearly shows that the solver is indeed very uncertain about the result. It is still not perfect and the samples look very different from the actual solution, but that might be the point where you would need non-Gaussian multi-modal distributions, e.g. provided by sampling-based methods.

I hope that helped, let me know if you have more questions! And of course, please feel free to open issues about anything that is unclear or that doesn't work.

[michaeldeistler]

Hi Nathanael,

thanks for your reply and your efforts. It's pretty cool that there is uncertainty now! You are raising a good point about the Gaussianity...I'm pretty sure the "true" uncertainty-aware likelihood is non-Gaussian. When adding Brownian motion to the differential equation (which, from my layman's perspective seems very similar to sampling-based approaches), the marginal histogram over samples is highly non-Gaussian.

Regarding dynamicMV: I guess the scenario of having dimensions on different scales is pretty common. Is there a reason to not make dynamicMV the default? Or at least give a warning if you detect that the scales are different?

Anyways, thanks a lot for your efforts, this is a really cool project!

Michael

[nathanaelbosch]

Hi Michael, I'm glad to hear that my reply was helpful.

The dynamicMV diffusion model unfortunately does not work with the EK1, which is typically my preferred solver, therefore I made the dynamic diffusion model the default. But this is a great example on why I might want to re-evaluate this choice. A warning is also a great idea to communicate these model differences better to users.

So, thanks again for trying out our library and for opening this issue! It provided very helpful input for us :) Feel free to reach out if you have any more questions or issues.