python - 高斯白噪声到霍奇金赫胥黎方程
问题描述
将高斯白噪声添加到经典的 Hodgkin Huxley ODE 中以获得膜电位
嗨,大家好,
我发现这个代码是通过霍奇金赫胥黎方程评估的膜电位。
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
#runtime in milliseconds
#print("Type in the runtime-interval (in milliseconds)")
#t_0 = input()
#t_1 = input()
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
#Sodium potential (mV)
V_NA = 50.0
#Potassium Potential
V_K = -77.0
#Leak Potential
V_L = -54.4
#Potassium channel conductance
g_k = 36.0
#Sodium channel conductance
g_NA = 120.0
#Leak channel conductance
g_L = 0.3
#Membrane capacitance
C = 1.0
#equally ('distributed') time values
T = np.linspace(t_0,t_1,10000)
#please input the stimulus (first of all with a fixed stimulus)
#one spike below
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):
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))/1.0
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()])
V_x = odeint(derivatives,X,T)
print(V_x)
plt.plot(T,V_x)
plt.xlabel('Time in ms')
plt.ylabel('Membrane Potential in mV')
我想添加一个高斯白噪声过程 - 但我不太确定如何,因为我不是一个伟大的 python 用户 -
也许
noise = numpy.random.normal(mean=0.0, std=1.0, size=T)
和
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))/1.0
+ noise
我可以简单地将它添加到微分方程吗?
亲切的问候和非常感谢提前
解决方案
首先,您需要将噪声参数添加到函数中,
def derivatives(x, z,noisyInput):
你已经用正态分布创造了你的噪音。
noisyInput = np.ones(length)*7 + np.random.normal(0,3,length)
您只需在除以容量的 dvdt 项中添加噪声项即可。
dx[0] = (noisyInput + 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))/1.0
你应该把噪音传递给函数。
V_x = odeint(derivatives,X,T,args=(noisyInput,))
你也可以查看我写的代码: https ://github.com/nosratullah/HodgkinHuxely
推荐阅读
- r - R工作室中的“调试位置是近似的,因为源不可用”
- c - 用int初始化ac作为EOF,然后将它用于getchar()?
- c++ - c ++ const_cast从类型'const BIGNUM *'到类型'BIGNUM**'
- ios - ios图像选择器委托未在nativescript中触发
- c - Integer division which results in less than 1
- vba - 将 office 主题设置为 powerpoint VBA
- html - 下拉悬停菜单未显示在 IE/Edge 中
- image - 从 ResourceStream 加载 PNG 图像的函数不返回任何内容
- if-statement - 使用 if 状态级联多个条件,以创建 2 个分类列
- javascript - C# ChromiumWebBrowser 获取图像(验证码)