霍奇金-赫胥黎(HH)神经元模型在微秒级双相脉冲刺激下膜电位响应异常咨询

先给你梳理下，你现在遇到的情况既有代码参数设置的硬错误，也有HH模型本身时间尺度匹配的问题，咱们一步步拆解：

一、代码里的关键参数错误（这是导致响应“不合理”的主要原因） #

你的代码里有几个不符合HH模型生理意义的参数设置，直接影响了膜电位的基础稳定性：

• 漏电流反转电位设置完全错误：你把EKleak设为10.6mV，但标准HH模型中，漏电流的反转电位应该接近神经元的静息电位（约-54.4mV）。这个正电位会导致漏电流持续驱动膜电位偏离正常范围，让整个模型的基础状态就不对。
• 初始膜电位与门控变量初始化bug：
  1. 你把起始膜电位startingVoltage设为0mV，但HH模型的正常静息电位在-60~-70mV之间，初始电位设为0mV会让门控变量一开始就处于去极化的异常状态；
  2. 在__init__函数里，你连续调用了两次self.n.setInfiniteState()，漏掉了self.h.setInfiniteState()，导致钠通道失活门h的初始状态是0（类默认值），但静息状态下h门应该处于高开放状态（约0.6），这会让钠通道的初始电导完全异常。
• 绘图时的电位偏移容易误导：你绘图时用了voltages - 70，把实际0mV的膜电位显示为-70mV，这会混淆你对实际膜电位水平的判断。

二、HH模型与微秒脉冲的时间尺度不匹配 #

HH模型的门控变量（m、n、h）的响应时间是毫秒级的：比如静息电位下，m门的时间常数τ大概在0.1_{1ms，`h`门的τ在1}10ms。而你用的是1微秒的脉冲，这个时间尺度远小于门控变量的响应时间——也就是说，在1微秒的脉冲持续时间内，门控变量几乎来不及发生任何变化，膜电位的变化仅仅是细胞膜电容的瞬时充放电（公式：ΔV = I*Δt/Cm）。

这种情况下，你看到的膜电位波动其实是电容对脉冲的直接响应，不会触发动作电位——因为动作电位需要足够长的去极化时间来激活钠通道（让m门打开、h门慢慢关闭），微秒脉冲的时间太短，根本达不到这个条件。

三、修正后的代码建议 #

我把核心错误的参数和初始化逻辑修正了，你可以运行试试：

"""
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, -54.4
    gNa, gK, gKleak = 120, 36, 0.3
    m, n, h = Gate(), Gate(), Gate()
    Cm = 1

    def __init__(self, startingVoltage=-65):  # 修正初始膜电位为静息电位附近
        self.Vm = startingVoltage
        self.UpdateGateTimeConstants(startingVoltage)
        self.m.setInfiniteState()
        self.n.setInfiniteState()
        self.h.setInfiniteState()  # 修复漏掉的h门初始化

    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)
    # 时间单位转换：每个点对应0.01微秒，100个点=1微秒
    times = np.arange(pointCount) * 0.00001
    stim = np.zeros(pointCount)

    start = 0
    num_repeats = 400
    amplitude = 100
    pulse_len = 100  # 对应1微秒
    pause_len = 100  # 对应1微秒
    total_cycle_len = 2 * (pulse_len + pause_len)

    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

    f, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8, 5),
                                 gridspec_kw={'height_ratios': [3, 1]})
    
    # 去掉不必要的偏移，直接绘制实际膜电位
    ax1.plot(times, voltages, 'b')
    ax1.set_ylabel("Membrane Potential (mV)")
    ax1.set_title("Hodgkin-Huxley Spiking Neuron Model with 1us Biphasic Pulses")
    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.show()

四、修正后的预期结果 #

1. 运行修正后的代码，你会看到膜电位稳定在-65mV左右，每次双相脉冲会引起微小的电位波动（先升后降，幅度约0.1mV），这是完全合理的——因为门控变量来不及响应，只有电容电流起作用；
2. 如果你想看到HH模型的“经典动作电位”，需要把脉冲长度改成毫秒级（比如把pulse_len设为10000，对应1ms脉冲），或者增大电流幅度并延长脉冲时间，让门控变量有足够时间响应。

备注：内容来源于stack exchange，提问作者marina05