how to overcome division by zero in the Hodgkin Huxley neuron model

I've implemented the Hodgkin Huxley model, which worked fine in my visualization, I got the spiking etc. Then I've build a network with a few thousand neurons and experienced random NAN (not a number) flooding the network; I could track it down to how sodium and potassium are updated. e.g.

VoltageDelta = MembranVoltage - ReferenceVoltage

AlphaNa = 0.01 * (VoltageDelta - 10.0) / (1.0 - exp( -(VoltageDelta-10.0) / 10.0))

If the voltage delta is exactly 10, the exp(0) returns 1, which ends up being a division by zero. I've compared my implementations to other, and despite some reordering of computations, all seem to have this flaw. I haven't found anything about how to deal with that. I could add a check and assume it's 0/0 -> 0, yet I'm not sure whether close to 0 values would lead to numerical inaccuracies and I wonder how this is usually handled. ( I've looked into some github implementations, and haven't seen special handling).

• Can you include the links to the github implementations and the other implementations that you found?
– honi
Feb 18, 2021 at 14:48

You have

AlphaNa = 0.1 * f(0.1*(VoltageDelta-10.0))


where $$f(x)=\frac{x}{1-\exp(-x)}$$ The limit for $$x\to 0$$ can actually be interpreted as difference quotient, you might also appeal to l'Hopital, to find that $$f(x)\to 1$$. So if you want to play it safe, you can employ a 3-tier or 4-tier implementation of this function

def f(x):
if abs(x)<1e-12: return 1
if abs(x)<1e-4: return -x/expm1(-x)
if x<-20: ex = exp(x); return -x*ex/(1-ex)
return x/(1-exp(-x))


The formula as it originally is can still produce non-representable intermediate values if $$-x$$ is moderately large so that $$\exp(-x)$$ produces a floating-point overflow. This should compensate automatically in the division, but to be sure I added the case x < -20.

• To expand on this answer, consider the Taylor series expansion of this function at $x=0$. wolframalpha.com/input/… Feb 19, 2021 at 16:33