python - Hodgkin Huxley model - incorrect membrane potential? - Stack Overflow

admin管理员组

文章数量:1405323

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

"""
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[i], deltaTms=0.00001)
        voltages[i] = 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:

本文标签： pythonHodgkin Huxley modelincorrect membrane potentialStack Overflow