opensourcebrain/HodgkinHuxleyTutorial
0
1import scipy as sp2import numpy as np3import pylab as plt4from scipy.integrate import odeint5import sys6 7class HodgkinHuxley():8 """Full Hodgkin-Huxley Model implemented in Python"""9 10 """ __init__ uses optional arguments """11 """ when no argument is passed default values are used """12 13 def __init__(self, C_m=1, gmax_Na=120, gmax_K=36, gmax_L=0.3, E_Na=50,14 E_K=-77, E_L=-54.387, t_n=450, delta_t=0.01,15 I_inj_amplitude=0, I_inj_duration=0, I_inj_delay=0,16 vc_delay=10, vc_duration=30, vc_condVoltage=-65,17 vc_testVoltage=10, vc_returnVoltage=-65, runMode='iclamp',18 injected_current_plot=True, gating_plot=True, cond_scaling_plot=False,19 cond_dens_plot=True, driving_force_plot=False,20 current_plot=True, memb_pot_plot=True):21 22 self.C_m = C_m23 """ membrane capacitance, in uF/cm^2 """24 25 self.gmax_Na = gmax_Na26 """ Sodium (Na) maximum conductances, in mS/cm^2 """27 28 self.gmax_K = gmax_K29 """ Postassium (K) maximum conductances, in mS/cm^2 """30 31 self.gmax_L = gmax_L32 """ Leak maximum conductances, in mS/cm^2 """33 34 self.E_Na = E_Na35 """ Sodium (Na) Nernst reversal potentials, in mV """36 37 self.E_K = E_K38 """ Postassium (K) Nernst reversal potentials, in mV """39 40 self.E_L = E_L41 """ Leak Nernst reversal potentials, in mV """42 43 self.t = np.arange(0, t_n, delta_t)44 """ The time to integrate over """45 46 """ Advanced input - injection current (single rectangular pulse only) """47 48 self.I_inj_amplitude = I_inj_amplitude49 """ maximum value or amplitude of injection pulse """50 51 self.I_inj_duration = I_inj_duration52 """ duration or width of injection pulse """53 54 self.I_inj_delay = I_inj_delay55 """ start time of injection pulse """56 57 #vclamp parameters58 self.run_mode = runMode59 """default is current clamp"""60 61 self.delay = vc_delay62 """Delay before switching from conditioningVoltage to testingVoltage, in ms"""63 64 self.duration = vc_duration65 """Duration to hold at testingVoltage, in ms"""66 67 self.conditioningVoltage = vc_condVoltage68 """Target voltage before time delay, in mV"""69 70 self.testingVoltage = vc_testVoltage71 """Target voltage between times delay and delay + duration, in mV"""72 73 self.returnVoltage = vc_returnVoltage74 """Target voltage after time duration, in mV"""75 76 self.simpleSeriesResistance = 1e777 """Current will be calculated by the difference in voltage between the target and parent, divided by this value, in mOhm"""78 79 # plotting conditionals80 self.injected_current_plot = injected_current_plot81 self.gating_plot = gating_plot82 self.cond_scaling_plot = cond_scaling_plot83 self.cond_dens_plot = cond_dens_plot84 self.driving_force_plot = driving_force_plot85 self.current_plot = current_plot86 self.memb_pot_plot = memb_pot_plot87 88 self.num_plots = (int(self.injected_current_plot) +89 int(self.gating_plot)+ int(self.cond_scaling_plot) +90 int(self.cond_dens_plot) + int(self.driving_force_plot) +91 int(self.current_plot) + int(self.memb_pot_plot))92 93 self.plot_count = 094 95 def alpha_m(self, V):96 """Channel gating kinetics. Functions of membrane voltage"""97 return 0.1*(V+40.0)/(1.0 - np.exp(-(V+40.0) / 10.0))98 99 def beta_m(self, V):100 """Channel gating kinetics. Functions of membrane voltage"""101 return 4.0*np.exp(-(V+65.0) / 18.0)102 103 def alpha_h(self, V):104 """Channel gating kinetics. Functions of membrane voltage"""105 return 0.07*np.exp(-(V+65.0) / 20.0)106 107 def beta_h(self, V):108 """Channel gating kinetics. Functions of membrane voltage"""109 return 1.0/(1.0 + np.exp(-(V+35.0) / 10.0))110 111 def alpha_n(self, V):112 """Channel gating kinetics. Functions of membrane voltage"""113 return 0.01*(V+55.0)/(1.0 - np.exp(-(V+55.0) / 10.0))114 115 def beta_n(self, V):116 """Channel gating kinetics. Functions of membrane voltage"""117 return 0.125*np.exp(-(V+65) / 80.0)118 119 def g_Na(self, m, h):120 """121 Conductance density (in mS/cm^2)122 Sodium (Na = element name)123 124 | :param m:125 | :param h:126 | :return:127 """128 return self.gmax_Na * m**3 * h129 130 def I_Na(self, V, m, h):131 """132 Membrane current (in uA/cm^2)133 Sodium (Na = element name)134 135 | :param V:136 | :param m:137 | :param h:138 | :return:139 """140 return self.g_Na(m, h) * (V - self.E_Na)141 142 143 def g_K(self, n):144 """145 Conductance density (in mS/cm^2)146 Potassium (K = element name)147 148 | :param n:149 | :return:150 """151 return self.gmax_K * n**4152 153 def I_K(self, V, n):154 """155 Membrane current (in uA/cm^2)156 Potassium (K = element name)157 158 | :param V:159 | :param n:160 | :return:161 """162 return self.g_K(n) * (V - self.E_K)163 164 # Leak165 def I_L(self, V):166 """167 Membrane current (in uA/cm^2)168 Leak169 170 | :param V:171 | :param h:172 | :return:173 """174 return self.gmax_L * (V - self.E_L)175 176 def I_inj(self, t):177 """178 External Current179 180 | :param t: time181 | :return: step up to 10 uA/cm^2 at t>100182 | step down to 0 uA/cm^2 at t>200183 | step up to 35 uA/cm^2 at t>300184 | step down to 0 uA/cm^2 at t>400185 """186 187 """ running standalone python script """188 if __name__ == '__main__':189 return 10*(t>100) - 10*(t>200) + 35*(t>300) - 35*(t>400)190 191 #""" running jupyterLab notebook """192 else:193 return self.I_inj_amplitude*(t>self.I_inj_delay) - self.I_inj_amplitude*(t>self.I_inj_delay+self.I_inj_duration)194 195 def I_inj_vclamp(self,t,v):196 """197 External Current (vclamp)198 199 | :param t: time200 | :return: injector current for voltage clamp201 |202 """203 if t > (self.delay + self.duration):204 current_A = (self.returnVoltage - v) / self.simpleSeriesResistance205 elif t >= self.delay:206 current_A = (self.testingVoltage - v) / self.simpleSeriesResistance207 elif t < self.delay:208 current_A = (self.conditioningVoltage - v) / self.simpleSeriesResistance209 else:210 print('Problem in injection current calculation for voltage clamp...')211 return 0212 213 #convert current to current density (uA/cm^2)214 current_uA = current_A*10**6 #convert ampere to micro ampere215 surface_area = 1000*10**-8 #surface area of 1000 um^2 converted to cm^2216 current_density = current_uA/surface_area217 218 return current_density219 220 @staticmethod221 def dALLdt(X, t, self):222 """223 Integrate224 225 | :param X:226 | :param t:227 | :return: calculate membrane potential & activation variables228 """229 V, m, h, n = X230 if self.is_vclamp():231 dVdt = (self.I_inj_vclamp(t,V) - self.I_Na(V, m, h) - self.I_K(V, n) - self.I_L(V)) / self.C_m232 else:233 dVdt = (self.I_inj(t) - self.I_Na(V, m, h) - self.I_K(V, n) - self.I_L(V)) / self.C_m234 235 dmdt = self.alpha_m(V)*(1.0-m) - self.beta_m(V)*m236 dhdt = self.alpha_h(V)*(1.0-h) - self.beta_h(V)*h237 dndt = self.alpha_n(V)*(1.0-n) - self.beta_n(V)*n238 return dVdt, dmdt, dhdt, dndt239 240 def is_vclamp(self):241 return self.run_mode=='vclamp' or self.run_mode=='Voltage Clamp'242 243 def simulate(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]):244 """245 Main simulate method for the Hodgkin Huxley neuron model246 """247 248 # init_values are the steady state values for v,m,h,n at zero current injection249 X = odeint(self.dALLdt, init_values, self.t, args=(self,))250 V = X[:,0]251 m = X[:,1]252 h = X[:,2]253 n = X[:,3]254 ina = self.I_Na(V, m, h)255 ik = self.I_K(V, n)256 il = self.I_L(V)257 gna = self.g_Na(m, h)258 gk = self.g_K(n)259 260 # Save some of the data to file261 with open('hh_py_v.dat','w') as f:262 for ti in range(len(self.t)):263 f.write('%s\t%s\n'%(self.t[ti],V[ti]))264 265 if not '-nogui' in sys.argv:266 #increase figure and font size for display in jupyter notebook267 268 269 if __name__ != '__main__':270 plt.rcParams['figure.figsize'] = [7, 7]271 #plt.rcParams['font.size'] = 15272 #plt.rcParams['legend.fontsize'] = 12273 plt.rcParams['legend.loc'] = "upper right"274 #275 else:276 plt.rcParams['figure.figsize'] = [10, 7]277 278 plt.close()279 280 fig=plt.figure(figsize=(7, self.num_plots * 2))281 fig.canvas.header_visible = False282 # plt.xlim([np.min(self.t),np.max(self.t)]) #for all subplots283 284 if self.injected_current_plot:285 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)286 plt.title('Simulation of Hodgkin Huxley model neuron')287 if self.is_vclamp():288 i_inj_values = [self.I_inj_vclamp(t,v) for t,v in zip(self.t,V)]289 else:290 i_inj_values = [self.I_inj(t) for t in self.t]291 292 if self.is_vclamp(): plt.ylim(-2000,3000)293 294 plt.plot(self.t, i_inj_values, 'k')295 plt.ylabel('$I_{inj}$ ($\\mu{A}/cm^2$)')296 297 self.plot_count += 1298 299 300 if self.gating_plot:301 try:302 plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)303 except NameError:304 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)305 plt.title('Simulation of Hodgkin Huxley model neuron')306 plt.plot(self.t, m, 'r', label='$m$')307 plt.plot(self.t, h, 'g', label='$h$')308 plt.plot(self.t, n, 'b', label='$n$')309 plt.ylabel('Gating variable')310 plt.legend()311 self.plot_count += 1312 313 if self.cond_scaling_plot:314 try:315 plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)316 except NameError:317 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)318 plt.title('Simulation of Hodgkin Huxley model neuron')319 scale_na = m*m*m*h320 scale_k = n*n*n*n321 plt.plot(self.t, scale_na, 'c', label='$m^{3}h$')322 plt.plot(self.t, scale_k, 'y', label='$n^{4}$')323 plt.ylabel('Cond scaling')324 plt.legend()325 self.plot_count += 1326 327 if self.cond_dens_plot:328 try:329 plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)330 except NameError:331 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)332 plt.title('Simulation of Hodgkin Huxley model neuron')333 plt.plot(self.t, gna, 'c', label='$g_{Na}$')334 plt.plot(self.t, gk, 'y', label='$g_{K}$')335 plt.ylabel('Cond dens ($mS/cm^2$)')336 plt.legend()337 self.plot_count += 1338 339 340 if self.driving_force_plot:341 try:342 ax_here = plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)343 except NameError:344 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)345 plt.title('Simulation of Hodgkin Huxley model neuron')346 ax_here = ax1347 348 dna = V - self.E_Na349 dk = V - self.E_K350 zero = [0 for v in V]351 352 #plt.plot(self.t, dna, 'c', label='$V - E_{Na}$')353 ax_here.fill_between(self.t, dna, color='c', alpha=0.5)354 ax_here.fill_between(self.t, dk, color='y', alpha=0.5)355 356 plt.plot(self.t, dna, 'c', label='$V_{m} - E_{Na}$', linewidth=0.8)357 plt.plot(self.t, dk, 'y', label='$V_{m} - E_{K}$', linewidth=0.8)358 plt.plot(self.t, zero, 'k', linestyle='dashed', linewidth=0.5)359 plt.ylabel('Driving force (mV)')360 plt.legend()361 #if not self.is_vclamp(): plt.ylim(-85,60)362 #plt.ylim(-1, 40)363 self.plot_count += 1364 365 if self.current_plot:366 try:367 plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)368 except NameError:369 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)370 plt.title('Simulation of Hodgkin Huxley model neuron')371 plt.plot(self.t, ina, 'c', label='$I_{Na}$')372 plt.plot(self.t, ik, 'y', label='$I_{K}$')373 plt.plot(self.t, il, 'm', label='$I_{L}$')374 plt.ylabel('Curr dens ($\\mu{A}/cm^2$)')375 plt.legend()376 self.plot_count += 1377 378 if self.memb_pot_plot:379 try:380 plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)381 except NameError:382 ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)383 plt.title('Simulation of Hodgkin Huxley model neuron')384 plt.plot(self.t, V, 'k')385 plt.ylabel('$V_{m}$ (mV)')386 plt.xlabel('Time (ms)')387 if not self.is_vclamp(): plt.ylim(-85,60)388 #plt.ylim(-1, 40)389 self.plot_count += 1390 391 plt.tight_layout()392 plt.show()393 394if __name__ == '__main__':395 396 if '-vclamp' in sys.argv:397 runner = HodgkinHuxley(runMode='vclamp', t_n=50, delta_t=0.0005)398 else: #default mode399 runner = HodgkinHuxley(runMode='iclamp', t_n=450, delta_t=0.01)400 401 runner.simulate()402