[Python] Hodgkin Huxley model - incorrect membrane potential?

I am using Hodgkin Huxley model of neuron and stimulating it with microsecond pulses. No matter what I try, I can't seem to get "logical" response of neuron aka membrane potential. Is this result logical or the code needs changes? - in this case I'm stimulating neuron with biphasic pulses of length 1 us and pauses are also 1 us long.

"""
Python implementation of the Hodgkin-Huxley spiking neuron model
https://github.com/swharden/pyHH
"""
import matplotlib.pyplot as plt
import numpy as np


class HHModel:
"""The HHModel tracks conductances of 3 channels to calculate Vm"""

class Gate:
"""The Gate object manages a channel's kinetics and open state"""
alpha, beta, state = 0, 0, 0

def update(self, deltaTms):
alphaState = self.alpha * (1-self.state)
betaState = self.beta * self.state
self.state += deltaTms * (alphaState - betaState)

def setInfiniteState(self):
self.state = self.alpha / (self.alpha + self.beta)

ENa, EK, EKleak = 115, -12, 10.6
gNa, gK, gKleak = 120, 36, 0.3
m, n, h = Gate(), Gate(), Gate()
Cm = 1

def __init__(self, startingVoltage=0):
self.Vm = startingVoltage
self.UpdateGateTimeConstants(startingVoltage)
self.m.setInfiniteState()
self.n.setInfiniteState()
self.n.setInfiniteState()

def UpdateGateTimeConstants(self, Vm):
"""Update time constants of all gates based on the given Vm"""
self.n.alpha = .01 * ((10-Vm) / (np.exp((10-Vm)/10)-1))
self.n.beta = .125*np.exp(-Vm/80)
self.m.alpha = .1*((25-Vm) / (np.exp((25-Vm)/10)-1))
self.m.beta = 4*np.exp(-Vm/18)
self.h.alpha = .07*np.exp(-Vm/20)
self.h.beta = 1/(np.exp((30-Vm)/10)+1)

def UpdateCellVoltage(self, stimulusCurrent, deltaTms):
"""calculate channel currents using the latest gate time constants"""
INa = np.power(self.m.state, 3) * self.gNa * \
self.h.state*(self.Vm-self.ENa)
IK = np.power(self.n.state, 4) * self.gK * (self.Vm-self.EK)
IKleak = self.gKleak * (self.Vm-self.EKleak)
Isum = stimulusCurrent - INa - IK - IKleak
self.Vm += deltaTms * Isum / self.Cm

def UpdateGateStates(self, deltaTms):
"""calculate new channel open states using latest Vm"""
self.n.update(deltaTms)
self.m.update(deltaTms)
self.h.update(deltaTms)

def Iterate(self, stimulusCurrent=0, deltaTms=0.05):
self.UpdateGateTimeConstants(self.Vm)
self.UpdateCellVoltage(stimulusCurrent, deltaTms)
self.UpdateGateStates(deltaTms)


if __name__ == "__main__":
hh = HHModel()

pointCount = 170000
voltages = np.empty(pointCount)
times = np.arange(pointCount) * 0.00001
stim = np.zeros(pointCount)

start = 0
num_repeats = 400
amplitude = 100
pulse_len = 100
pause_len = 100
total_cycle_len = 2 * (pulse_len + pause_len) # 21000

for i in range(num_repeats):
idx = start + i * total_cycle_len
stim[idx : idx + pulse_len] = amplitude # + pulse
stim[idx + pulse_len : idx + pulse_len + pause_len] = 0 # pause
stim[idx + pulse_len + pause_len : idx + 2*pulse_len + pause_len] = -amplitude # - pulse
stim[idx + 2*pulse_len + pause_len : idx + total_cycle_len] = 0 # pause


print(f"Initial m.state: {hh.m.state}, n.state: {hh.n.state}, h.state: {hh.h.state}")

for i in range(len(times)):
hh.Iterate(stimulusCurrent=stim, deltaTms=0.00001)
voltages = hh.Vm
# note: you could also plot hh's n, m, and k (channel open states)

f, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8, 5),
gridspec_kw={'height_ratios': [3, 1]})


print(len(times))
print(times[0:3])
print(len(stim))
print(len(voltages))

ax1.plot(times, voltages - 70, 'b')
ax1.set_ylabel("Membrane Potential (mV)")
ax1.set_title("Hodgkin-Huxley Spiking Neuron Model")
ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
ax1.spines['bottom'].set_visible(False)
ax1.tick_params(bottom=False)

ax2.plot(times, stim, 'r')
ax2.set_ylabel("Stimulus (µA/cm²)")
ax2.set_xlabel("Simulation Time (milliseconds)")
ax2.spines['right'].set_visible(False)
ax2.spines['top'].set_visible(False)

plt.margins(0, 0.1)
plt.tight_layout()
# plt.savefig("dev/concept4.png")
plt.show()


The response I get is: and zoomed:

Continue reading...