如何为经典Hodgkin-Huxley膜电位ODE添加高斯白噪声？

首先直接回答你的问题：直接将高斯白噪声项加入膜电位的微分方程dx[0]中是一种可行的近似实现方式，不过需要注意白噪声在连续时间系统里的特性，以及数值求解时的细节。

为什么这个近似可行？ #

高斯白噪声在连续时间下对应维纳过程（布朗运动）的导数，而标准的odeint是用来求解确定性常微分方程（ODE）的。当你的时间步长足够小（像你代码里用了10000个点覆盖200ms，步长0.02ms，已经比较小），直接把离散的高斯噪声样本加到dx[0]里，能近似模拟膜电位受到的随机扰动。

修改后的代码示例 #

你只需要生成和时间数组T长度一致的噪声，然后在derivatives函数里把噪声项加到dx[0]中即可。注意要根据需求调整噪声的标准差（std），避免噪声过大掩盖HH模型本身的电活动：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

# runtime in milliseconds
t_0 = 0.0
t_1 = 200.0
print("The interval of runtime is:")
print(t_0, t_1)

# Potassium (alpha_n, beta_n) and Sodium (alpha_m, beta_m, alpha_h, beta_h) ion-channel rate functions
def alpha_n(V):
    return (0.01*(10.0-V))/(np.exp((10.0-V)/10.0)-1.0)
def alpha_m(V):
    return (0.1*(25.0-V))/(np.exp((25.0-V)/10.0)-1)
def beta_n(V):
    return 0.125*np.exp(-V/80.0)
def beta_m(V):
    return 4.0*np.exp(-V/18.0)
def alpha_h(V):
    return 0.07*np.exp(-V/20.0)
def beta_h(V):
    return (1.0)/((np.exp((30.0-V)/10.0))+1)

# mostly used parameters
V_NA = 50.0    # Sodium potential (mV)
V_K = -77.0    # Potassium Potential
V_L = -54.4    # Leak Potential
g_k = 36.0     # Potassium channel conductance
g_NA = 120.0   # Sodium channel conductance
g_L = 0.3      # Leak channel conductance
C = 1.0        # Membrane capacitance

# equally distributed time values
T = np.linspace(t_0,t_1,10000)
# 生成高斯白噪声：均值0，标准差可调整，这里用0.5作为示例
noise = np.random.normal(loc=0.0, scale=0.5, size=len(T))

# stimulus function
def stim(t):
    if 0.0 < t < 1.0:
        return 150.0
    elif 35.0 < t < 36.0 :
        return 5000.0
    return 0.0

# steady-state values
def n_infty(V = 0.0):
    return alpha_n(V)/(alpha_n(V)+beta_n(V))
def m_infty(V = 0.0):
    return alpha_m(V)/(alpha_m(V)+beta_m(V))
def h_infty(V = 0.0):
    return alpha_h(V)/(alpha_h(V)+beta_h(V))
def tau_m(V = 0.0):
    return 1.0/(alpha_m(V)+beta_m(V))
def tau_n(V = 0.0):
    return 1.0/(alpha_n(V)+beta_n(V))
def tau_h(V = 0.0):
    return 1.0/(alpha_h(V)+beta_h(V))

# 修改derivatives函数，加入噪声项
def derivatives(x, z, noise):
    # 找到当前时间z对应的噪声索引
    idx = np.argmin(np.abs(T - z))
    dx = np.zeros((4,))
    V = x[0]
    m = x[1]
    n = x[2]
    h = x[3]
    # 把噪声项加到膜电位的微分方程中
    dx[0] = (stim(z)-g_NA*np.power(m,3.0)*h*(V-V_NA)-g_k*np.power(n,4.0)*(V-V_K)-g_L*(V-V_L) + noise[idx])/C
    dx[1] = (m_infty(V)-m)/tau_m(V)
    dx[2] = (n_infty(V)-n)/tau_n(V)
    dx[3] = (h_infty(V)-h)/tau_h(V)
    return dx

# 初始条件
X = np.array([0.0,m_infty(),n_infty(),h_infty()])
# 调用odeint时传递noise参数
V_x = odeint(derivatives, X, T, args=(noise,))

# 绘图：只画膜电位（第一列）
plt.plot(T, V_x[:,0])
plt.xlabel('Time in ms')
plt.ylabel('Membrane Potential in mV')
plt.title('Hodgkin-Huxley Model with Gaussian White Noise')
plt.show()

更严谨的处理方式：随机微分方程（SDE） #

如果你需要更准确地模拟随机扰动，应该把HH模型作为随机微分方程（SDE）来求解，因为白噪声本质上是连续时间的随机过程，odeint不是专门的SDE求解器。你可以使用sdeint库（需要先安装：pip install sdeint），它专门用于求解随机微分方程。

示例代码框架大概是这样：

import sdeint

# 定义确定性部分（原HH方程的导数）
def f(x, t):
    V, m, n, h = x
    dV = (stim(t)-g_NA*m**3*h*(V-V_NA)-g_k*n**4*(V-V_K)-g_L*(V-V_L))/C
    dm = (m_infty(V)-m)/tau_m(V)
    dn = (n_infty(V)-n)/tau_n(V)
    dh = (h_infty(V)-h)/tau_h(V)
    return np.array([dV, dm, dn, dh])

# 定义噪声部分：这里假设只有膜电位受到噪声，噪声强度为sigma
sigma = 0.5
def G(x, t):
    return np.array([sigma, 0, 0, 0])  # 只有V的方程有噪声项

# 求解SDE
V_x_sde = sdeint.itoint(f, G, X, T)

注意事项 #

• 噪声的标准差（scale参数）要合理设置：如果太大，会导致膜电位的随机波动掩盖原本的动作电位；太小的话，扰动效果不明显。
• 用odeint加离散噪声时，时间步长越小，近似效果越好，因为离散噪声样本更接近连续白噪声的特性。
• 如果你需要重复实验，记得设置随机种子（np.random.seed(42)），这样每次生成的噪声是一致的，方便对比结果。

内容的提问来源于stack exchange，提问作者user10674336