Neuronal Dynamics：Hodgkin-Huxley model 实验结果展示

for-nothing

已于 2023-01-12 21:43:05 修改

阅读量301

点赞数 1

分类专栏： Neuronal Dynamics 文章标签： python 开发语言

于 2023-01-09 21:20:11 首次发布

本文链接：https://blog.csdn.net/for_none/article/details/128620325

版权

Neuronal Dynamics 专栏收录该内容

13 篇文章 4 订阅

订阅专栏

文章介绍了使用Python实现的Hodgkin-Huxley模型，该模型用于模拟神经元的电生理行为。代码展示了如何通过调整输入信号来观察神经元膜电位u、激活变量m、h、n随时间的变化。不同输入条件下，模型产生的输出包括静息状态、动作电位和不同类型的神经响应，符合电生理学的预期结果。

摘要由CSDN通过智能技术生成

Neuronal Dynamics：Hodgkin-Huxley model 实验结果展示

一、演示效果和代码
二、模型及参数
三、输出解释

一、演示效果和代码

如下代码运行完之后会在本地产生一个显示输出内容的gif文件（4s一次切换）
效果如下,gif包含一页参数图和四页不同的输入-输出图，后文会详细介绍：
请添加图片描述

import numpy as np
import time
from matplotlib import pyplot as plt
from matplotlib.pyplot import *
import matplotlib.animation
rc('font',size=16)
rc('font',family='SimHei')
plt.rcParams['axes.unicode_minus'] = False
class Hodgkin_Huxley():
    def __init__(self):
        self.urest=-65
        self.u=-65
        # self.Ena=55
        # self.Ek=-77
        # self.El=-65
        # self.gna=40
        # self.gk=35
        # self.gl=0.3
        self.Ena = 50
        self.Ek = -77
        self.El = -54.4
        self.gna = 120
        self.gk = 36
        self.gl = 0.3

        self.input=0
        self.C=1
        self.m=self.amu()/(self.amu()+self.bmu())
        self.h=self.ahu()/(self.ahu()+self.bhu())
        self.n=self.anu()/(self.anu()+self.bnu())
    def anu(self):
        vm = -self.urest+self.u
        return 0.01*(10-vm)/(np.exp(1-vm/10)-1)
    def amu(self):
        vm = -self.urest+self.u
        return 0.1*(25-vm)/(np.exp(2.5-vm/10)-1)
    def ahu(self):
        vm = -self.urest+self.u
        return 0.07*np.exp(-vm/20)
    def bnu(self):
        vm = -self.urest+self.u
        return 0.125*np.exp(-vm/80)
    def bmu(self):
        vm = -self.urest+self.u
        return 4*np.exp(-vm/18)
    def bhu(self):
        vm = -self.urest+self.u
        return 1/(1+np.exp(3-vm/10))
    # def anu(self):
    #     return 0.02*(self.u-25)/(1-np.exp(-(self.u-25)/9))
    # def amu(self):
    #     return 0.182*(self.u+35)/(1-np.exp(-(self.u+35)/9))
    # def ahu(self):
    #     return 1/(1+np.exp(-(self.u+62)/6))
    # def bnu(self):
    #     return 0.002*(self.u-25)/(1-np.exp(-(self.u-25)/9))
    # def bmu(self):
    #     return 0.124*(self.u+35)/(1-np.exp(-(self.u+35)/9))
    # def bhu(self):
    #     return 4*np.exp((self.u+90)/12)/(1+np.exp(-(self.u+62)/6))
    def update(self,input,dt):
        self.input=input
        Ina = self.gna*pow(self.m,3)*self.h*(self.u-self.Ena)
        Ik = self.gk*pow(self.n,4)*(self.u-self.Ek)
        Il = self.gl*(self.u-self.El)
        du=(self.input-Ina-Ik-Il)/self.C
        dm=self.amu()*(1-self.m)-self.bmu()*self.m
        dh=self.ahu()*(1-self.h)-self.bhu()*self.h
        dn=self.anu()*(1-self.n)-self.bnu()*self.n
        # print('参数：',self.u,self.m,self.h,self.n)
        # print('参数变化1',du,dm,dh,dn)
        # print('电流参数2',Ina,Ik,Il)
        self.u +=du*dt
        self.m +=dm*dt
        self.h +=dh*dt
        self.n +=dn*dt
def plot_line_pic(x,y,label='u-t变化图',ylim1=0,ylim2=0):
    plt.plot(x, y, '--r',linewidth=1.6, markeredgecolor='k', markersize=8)
    title(label,fontsize=40)
    if ylim1!=0 or ylim2!=0:
        ylim(ylim1,ylim2)
def caculate_data(in_list,dt_big = 0.1,dt_small = 0.001,all_time = 50):
    model = Hodgkin_Huxley()
    u_list = []
    m_list = []
    h_list = []
    n_list = []
    t_list = []
    for i in range(int(all_time / dt_big)):
        u_list.append(model.u)
        m_list.append(model.m)
        h_list.append(model.h)
        n_list.append(model.n)
        t_list.append(i * dt_big)
        for j in range(int(dt_big / dt_small)):
            model.update(in_list[i], dt_small)
    plt.subplot(222)
    plot_line_pic(t_list, in_list,'input-t变化图',-50,300)
    plt.subplot(221)
    plot_line_pic(t_list, u_list,'u-t变化图',-100,100)
    plt.subplot(234)
    plot_line_pic(t_list, m_list,'m-t变化图',0,1)
    plt.subplot(235)
    plot_line_pic(t_list, h_list,'h-t变化图',0,1)
    plt.subplot(236)
    plot_line_pic(t_list, n_list,'n-t变化图',0,1)
if __name__ == '__main__':
    model = Hodgkin_Huxley()
    u0_list = []
    mu_list = []
    hu_list = []
    nu_list = []
    mtu_list = []
    htu_list = []
    ntu_list = []
    for i in range(-100,50):
        u0_list.append(i)
        model.u=i
        mu_list.append(model.amu()/(model.amu()+model.bmu()))
        hu_list.append(model.ahu()/(model.ahu()+model.bhu()))
        nu_list.append(model.anu()/(model.anu()+model.bnu()))
        mtu_list.append(1/(model.amu()+model.bmu()))
        htu_list.append(1/(model.ahu()+model.bhu()))
        ntu_list.append(1/(model.anu()+model.bnu()))
    global fig
    fig,ax = plt.subplots(figsize=(28,16))
    def update(i):
        if i==0:
            plt.subplot(111)
            plt.cla()
            plt.subplot(231)
            plot_line_pic(u0_list,mu_list,'m-u变化图',0,1)
            plt.subplot(232)
            plot_line_pic(u0_list,hu_list,'h-u变化图',0,1)
            plt.subplot(233)
            plot_line_pic(u0_list,nu_list,'n-u变化图',0,1)
            plt.subplot(234)
            plot_line_pic(u0_list,mtu_list,'tm-u变化图')
            plt.subplot(235)
            plot_line_pic(u0_list,htu_list,'th-u变化图')
            plt.subplot(236)
            plot_line_pic(u0_list,ntu_list,'tn-u变化图')
            plt.suptitle('parameters',fontsize=60)
        elif i==1:
            dt_big = 0.1
            all_time = 50
            in_list=[]
            for i in range(int(all_time / dt_big)):
                if 10<=i*dt_big and i*dt_big<=10.5:
                    in_list.append(250)
                else:in_list.append(0)
            plt.subplot(111)
            plt.cla()
            plt.suptitle('ouput sample1',fontsize=60)
            caculate_data(in_list,dt_big = 0.1, dt_small = 0.001, all_time = 50)
        elif i==2:
            dt_big = 0.1
            all_time = 50
            in_list=[]
            for i in range(int(all_time / dt_big)):
                if 10<=i*dt_big and i*dt_big<=40:
                    in_list.append(50)
                else:in_list.append(0)
            plt.subplot(111)
            plt.cla()
            plt.suptitle('ouput sample2',fontsize=60)
            caculate_data(in_list,dt_big = 0.1, dt_small = 0.001, all_time = 50)
        elif i==3:
            dt_big = 0.1
            all_time = 50
            in_list=[]
            for i in range(int(all_time / dt_big)):
                if 10<=i*dt_big and i*dt_big<=40:
                    in_list.append(250)
                else:in_list.append(0)
            plt.subplot(111)
            plt.cla()
            plt.suptitle('ouput sample3',fontsize=60)
            caculate_data(in_list,dt_big = 0.1, dt_small = 0.001, all_time = 50)
        elif i==4:
            dt_big = 0.1
            all_time = 50
            in_list=[]
            for i in range(int(all_time / dt_big)):
                if 10<=i*dt_big and i*dt_big<=20:
                    in_list.append(-10)
                else:in_list.append(0)
            plt.subplot(111)
            plt.cla()
            plt.suptitle('ouput sample4',fontsize=60)
            caculate_data(in_list,dt_big = 0.1, dt_small = 0.001, all_time = 50)
    ani = matplotlib.animation.FuncAnimation(fig, update, frames=5,interval=4000, repeat=True)
    ani.save("Hodgkin_Huxley_model.gif", writer='imgemagick')

二、模型及参数

模型核心为以下四个式子：
$C_m\frac{du}{dt}=-g_{Na}m^3h(u-V_{Na})-g_Kn^4(u-V_K)-g_L(u-V_L)+I(t)$

$\frac{dn}{dt}=α~n~(u)(1-n)-β(u)n$

$\frac{dm}{dt}=α~m~(u)(1-m)-β(u)m$

$\frac{dh}{dt}=α~h~(u)(1-h)-β(u)h$

x(u),x=m,n,k为 : $\frac{α_x(u)}{α_x(u)+β_x(u)}$
τ_x(u),x=m,n,k为: $\frac{1}{α_x(u)+β_x(u)}$

书上给的参数是
C_m=1μF,U_rest=-65mV

x	E_x[mV]	g_x[mS/cm²]
Na	55	40
K	-77	35
L	-65	0.3

x	α_x(u/mV)[ms^-1]	β_x(u/mV)[ms^-1]
n	0.02(u-25)/[1-e^-(u-25)/9]	-0.002(u-25)/[1-e^-(u-25)/9]
m	0.182(u+35)/[1-e^-(u+35)/9]	-0.124(u+35)/[1-e^-(u+35)/9]
h	1/[1+e^-(u+62)/6]	4e^(u+90)/12/[1+e^-(u+62)/6]

但是我感觉这个参数有问题，实际带入模型算不出来，会出错（尤其是n(u)和m(u)是常数，这显然有问题）
所以我采用了网络上找到的如下参数：
C_m=1μF,U_rest=-65mV,V=u-U_rest

x	E_x[mV]	g_x[mS/cm²]
Na	50	120
K	-77	36
L	-54.4	0.3

x	α_x(u/mV)[ms^-1]	β_x(u/mV)[ms^-1]
n	0.01(10-V)/[e^(10-u)/10-1]	0.125e^-V/80
m	0.1(25-V)/[e^(25-V)/10-1]	4e^-V/18
h	0.07e^-V/20	1/[e^(30-V)/10+1]

三、输出解释

在这里插入图片描述
第一张图是根据上述式子计算出的x(x=m,h,n)和tx的图像

第二张图是模拟输入是一个冲击信号的时候输出的电压u以及各个参数m，n，h的变化
图像中u随时间t先升高后降低到静息电位以下的值然后缓慢恢复，模型结果符合结论
在这里插入图片描述
第三张图是模拟输入是一个较大的矩形脉冲信号的时候输出的电压u以及各个参数m，n，h的变化
图像中u随时间t进行阻尼振荡，模型结果符合结论

第四张图是模拟输入是一个很大的矩形脉冲信号的时候输出的电压u以及各个参数m，n，h的变化
图像中u先产生一个动作电位，由于输入很大，没有产生振动而是维持在一个较高的水平，知道输入消失，才下降至静息电位以下缓慢恢复，模型结果符合结论
在这里插入图片描述
第五张图是模拟输入是一个抑制类型的矩形脉冲信号的时候输出的电压u以及各个参数m，n，h的变化
图像中u先缓慢的下降到一定值，当足够长时间的抑制输入消失的时候，反弹产生了一个动作电位，模型结果符合结论