import scipy as sp import numpy as np import pylab as plt from scipy.integrate import odeint import sys class HodgkinHuxley(): """Full Hodgkin-Huxley Model implemented in Python""" """ __init__ uses optional arguments """ """ when no argument is passed default values are used """ def __init__(self, C_m=1, gmax_Na=120, gmax_K=36, gmax_L=0.3, E_Na=50, E_K=-77, E_L=-54.387, t_n=450, delta_t=0.01, I_inj_amplitude=0, I_inj_duration=0, I_inj_delay=0, vc_delay=10, vc_duration=30, vc_condVoltage=-65, vc_testVoltage=10, vc_returnVoltage=-65, runMode='iclamp', injected_current_plot=True, gating_plot=True, cond_scaling_plot=False, cond_dens_plot=True, driving_force_plot=False, current_plot=True, memb_pot_plot=True): self.C_m = C_m """ membrane capacitance, in uF/cm^2 """ self.gmax_Na = gmax_Na """ Sodium (Na) maximum conductances, in mS/cm^2 """ self.gmax_K = gmax_K """ Postassium (K) maximum conductances, in mS/cm^2 """ self.gmax_L = gmax_L """ Leak maximum conductances, in mS/cm^2 """ self.E_Na = E_Na """ Sodium (Na) Nernst reversal potentials, in mV """ self.E_K = E_K """ Postassium (K) Nernst reversal potentials, in mV """ self.E_L = E_L """ Leak Nernst reversal potentials, in mV """ self.t = np.arange(0, t_n, delta_t) """ The time to integrate over """ """ Advanced input - injection current (single rectangular pulse only) """ self.I_inj_amplitude = I_inj_amplitude """ maximum value or amplitude of injection pulse """ self.I_inj_duration = I_inj_duration """ duration or width of injection pulse """ self.I_inj_delay = I_inj_delay """ start time of injection pulse """ #vclamp parameters self.run_mode = runMode """default is current clamp""" self.delay = vc_delay """Delay before switching from conditioningVoltage to testingVoltage, in ms""" self.duration = vc_duration """Duration to hold at testingVoltage, in ms""" self.conditioningVoltage = vc_condVoltage """Target voltage before time delay, in mV""" self.testingVoltage = vc_testVoltage """Target voltage between times delay and delay + duration, in mV""" self.returnVoltage = vc_returnVoltage """Target voltage after time duration, in mV""" self.simpleSeriesResistance = 1e7 """Current will be calculated by the difference in voltage between the target and parent, divided by this value, in mOhm""" # plotting conditionals self.injected_current_plot = injected_current_plot self.gating_plot = gating_plot self.cond_scaling_plot = cond_scaling_plot self.cond_dens_plot = cond_dens_plot self.driving_force_plot = driving_force_plot self.current_plot = current_plot self.memb_pot_plot = memb_pot_plot self.num_plots = (int(self.injected_current_plot) + int(self.gating_plot)+ int(self.cond_scaling_plot) + int(self.cond_dens_plot) + int(self.driving_force_plot) + int(self.current_plot) + int(self.memb_pot_plot)) self.plot_count = 0 def alpha_m(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 0.1*(V+40.0)/(1.0 - np.exp(-(V+40.0) / 10.0)) def beta_m(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 4.0*np.exp(-(V+65.0) / 18.0) def alpha_h(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 0.07*np.exp(-(V+65.0) / 20.0) def beta_h(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 1.0/(1.0 + np.exp(-(V+35.0) / 10.0)) def alpha_n(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 0.01*(V+55.0)/(1.0 - np.exp(-(V+55.0) / 10.0)) def beta_n(self, V): """Channel gating kinetics. Functions of membrane voltage""" return 0.125*np.exp(-(V+65) / 80.0) def g_Na(self, m, h): """ Conductance density (in mS/cm^2) Sodium (Na = element name) | :param m: | :param h: | :return: """ return self.gmax_Na * m**3 * h def I_Na(self, V, m, h): """ Membrane current (in uA/cm^2) Sodium (Na = element name) | :param V: | :param m: | :param h: | :return: """ return self.g_Na(m, h) * (V - self.E_Na) def g_K(self, n): """ Conductance density (in mS/cm^2) Potassium (K = element name) | :param n: | :return: """ return self.gmax_K * n**4 def I_K(self, V, n): """ Membrane current (in uA/cm^2) Potassium (K = element name) | :param V: | :param n: | :return: """ return self.g_K(n) * (V - self.E_K) # Leak def I_L(self, V): """ Membrane current (in uA/cm^2) Leak | :param V: | :param h: | :return: """ return self.gmax_L * (V - self.E_L) def I_inj(self, t): """ External Current | :param t: time | :return: step up to 10 uA/cm^2 at t>100 | step down to 0 uA/cm^2 at t>200 | step up to 35 uA/cm^2 at t>300 | step down to 0 uA/cm^2 at t>400 """ """ running standalone python script """ if __name__ == '__main__': return 10*(t>100) - 10*(t>200) + 35*(t>300) - 35*(t>400) #""" running jupyterLab notebook """ else: return self.I_inj_amplitude*(t>self.I_inj_delay) - self.I_inj_amplitude*(t>self.I_inj_delay+self.I_inj_duration) def I_inj_vclamp(self,t,v): """ External Current (vclamp) | :param t: time | :return: injector current for voltage clamp | """ if t > (self.delay + self.duration): current_A = (self.returnVoltage - v) / self.simpleSeriesResistance elif t >= self.delay: current_A = (self.testingVoltage - v) / self.simpleSeriesResistance elif t < self.delay: current_A = (self.conditioningVoltage - v) / self.simpleSeriesResistance else: print('Problem in injection current calculation for voltage clamp...') return 0 #convert current to current density (uA/cm^2) current_uA = current_A*10**6 #convert ampere to micro ampere surface_area = 1000*10**-8 #surface area of 1000 um^2 converted to cm^2 current_density = current_uA/surface_area return current_density @staticmethod def dALLdt(X, t, self): """ Integrate | :param X: | :param t: | :return: calculate membrane potential & activation variables """ V, m, h, n = X if self.is_vclamp(): dVdt = (self.I_inj_vclamp(t,V) - self.I_Na(V, m, h) - self.I_K(V, n) - self.I_L(V)) / self.C_m else: dVdt = (self.I_inj(t) - self.I_Na(V, m, h) - self.I_K(V, n) - self.I_L(V)) / self.C_m dmdt = self.alpha_m(V)*(1.0-m) - self.beta_m(V)*m dhdt = self.alpha_h(V)*(1.0-h) - self.beta_h(V)*h dndt = self.alpha_n(V)*(1.0-n) - self.beta_n(V)*n return dVdt, dmdt, dhdt, dndt def is_vclamp(self): return self.run_mode=='vclamp' or self.run_mode=='Voltage Clamp' def simulate(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]): """ Main simulate method for the Hodgkin Huxley neuron model """ # init_values are the steady state values for v,m,h,n at zero current injection X = odeint(self.dALLdt, init_values, self.t, args=(self,)) V = X[:,0] m = X[:,1] h = X[:,2] n = X[:,3] ina = self.I_Na(V, m, h) ik = self.I_K(V, n) il = self.I_L(V) gna = self.g_Na(m, h) gk = self.g_K(n) # Save some of the data to file with open('hh_py_v.dat','w') as f: for ti in range(len(self.t)): f.write('%s\t%s\n'%(self.t[ti],V[ti])) if not '-nogui' in sys.argv: #increase figure and font size for display in jupyter notebook if __name__ != '__main__': plt.rcParams['figure.figsize'] = [7, 7] #plt.rcParams['font.size'] = 15 #plt.rcParams['legend.fontsize'] = 12 plt.rcParams['legend.loc'] = "upper right" # else: plt.rcParams['figure.figsize'] = [10, 7] plt.close() fig=plt.figure(figsize=(7, self.num_plots * 2)) fig.canvas.header_visible = False # plt.xlim([np.min(self.t),np.max(self.t)]) #for all subplots if self.injected_current_plot: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') if self.is_vclamp(): i_inj_values = [self.I_inj_vclamp(t,v) for t,v in zip(self.t,V)] else: i_inj_values = [self.I_inj(t) for t in self.t] if self.is_vclamp(): plt.ylim(-2000,3000) plt.plot(self.t, i_inj_values, 'k') plt.ylabel('$I_{inj}$ ($\\mu{A}/cm^2$)') self.plot_count += 1 if self.gating_plot: try: plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') plt.plot(self.t, m, 'r', label='$m$') plt.plot(self.t, h, 'g', label='$h$') plt.plot(self.t, n, 'b', label='$n$') plt.ylabel('Gating variable') plt.legend() self.plot_count += 1 if self.cond_scaling_plot: try: plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') scale_na = m*m*m*h scale_k = n*n*n*n plt.plot(self.t, scale_na, 'c', label='$m^{3}h$') plt.plot(self.t, scale_k, 'y', label='$n^{4}$') plt.ylabel('Cond scaling') plt.legend() self.plot_count += 1 if self.cond_dens_plot: try: plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') plt.plot(self.t, gna, 'c', label='$g_{Na}$') plt.plot(self.t, gk, 'y', label='$g_{K}$') plt.ylabel('Cond dens ($mS/cm^2$)') plt.legend() self.plot_count += 1 if self.driving_force_plot: try: ax_here = plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') ax_here = ax1 dna = V - self.E_Na dk = V - self.E_K zero = [0 for v in V] #plt.plot(self.t, dna, 'c', label='$V - E_{Na}$') ax_here.fill_between(self.t, dna, color='c', alpha=0.5) ax_here.fill_between(self.t, dk, color='y', alpha=0.5) plt.plot(self.t, dna, 'c', label='$V_{m} - E_{Na}$', linewidth=0.8) plt.plot(self.t, dk, 'y', label='$V_{m} - E_{K}$', linewidth=0.8) plt.plot(self.t, zero, 'k', linestyle='dashed', linewidth=0.5) plt.ylabel('Driving force (mV)') plt.legend() #if not self.is_vclamp(): plt.ylim(-85,60) #plt.ylim(-1, 40) self.plot_count += 1 if self.current_plot: try: plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') plt.plot(self.t, ina, 'c', label='$I_{Na}$') plt.plot(self.t, ik, 'y', label='$I_{K}$') plt.plot(self.t, il, 'm', label='$I_{L}$') plt.ylabel('Curr dens ($\\mu{A}/cm^2$)') plt.legend() self.plot_count += 1 if self.memb_pot_plot: try: plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1) except NameError: ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1) plt.title('Simulation of Hodgkin Huxley model neuron') plt.plot(self.t, V, 'k') plt.ylabel('$V_{m}$ (mV)') plt.xlabel('Time (ms)') if not self.is_vclamp(): plt.ylim(-85,60) #plt.ylim(-1, 40) self.plot_count += 1 plt.tight_layout() plt.show() if __name__ == '__main__': if '-vclamp' in sys.argv: runner = HodgkinHuxley(runMode='vclamp', t_n=50, delta_t=0.0005) else: #default mode runner = HodgkinHuxley(runMode='iclamp', t_n=450, delta_t=0.01) runner.simulate()