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).
Edit: some implementations you can find: https://github.com/swharden/HHSharp/blob/master/src/HHSharp/HHModel.cs https://github.com/openworm/hodgkin_huxley_tutorial/blob/master/Tutorial/Source/HodgkinHuxley.py