import matplotlib.pyplot as plt import numpy as np from scipy.signal import find_peaks

def HodHux_ELEC380(I, dt, plot_activation = 0): #Hodgkin Huxley model of a neuron #ELEC 480 Fall 2014 #JTRobinson # #Converted from MATLAB to Python by CATuppen Fall 2021 # #Input: #I = a vector of current values [uA] #dt = time step between I measurements [ms] #plot_activation (optional) = 1 to plot the activation state variables vs V (default 0) # #Output: #V = membrane voltage in mV #n = sodium activation #m = potassium activation #h = 1 - potassium inactivation #Im = membrane current in uA # #State variables: V, n, m, h # #DV = 1/C (I - Ik - Ina - Il) # # Il = gl(V-El) # Ik = gk*n^4(V-Ek) # Ina = gna*m^3*h(V-Ena) # #Dn = (ninf(V) - n)/taun(V) #Dm = (minf(V) - m)/taum(V) #Dh = (hinf(V) - h)/tauh(V) # can move these to outside of function for efficiency if desired global V global m global n global h global gkt global gnat global Ik global Ina global Il #starting membrane potential; Vstart = -65.625 #[mV] #Constants: #reversal potentials for various ions Ek = -77 #[mV] Ena = 50 #[mV] El = -54.4 #[mV] H&H used approx -55 mV #Membrane capacitance: C = 1 #[uF] #maximum conductances [mS] gna = 120 gk = 36 gl = 0.3 #Constants for GHK eq for channel activation #from J. Bossu, J Physiol. 496.2 (1996) #V50n = 57; #kn = 13.5; #from Izhikevich pg. 46 (2007) #K+ Delayed Rectifier 1: Ik = g*n^4*(V-Ek) #activation V50n = -50 kn = 15 Vmaxn = -79 sig_n = 50 Campn = 4.7 Cbasen = 1.1 #from Izhikevich pg. 46 (2007) #Na+ Fast Transient 1: Ina = g*m^3*h*(V-Ena) #activation: V50m = -40 km = 8 #adjusted to match Fig. 2.3 on Pg. 39. probably a typo in the table Vmaxm = -38 sig_m = 30 Campm = 0.46 Cbasem = 0.04 #inactivation V50h = -60 kh = -7 Vmaxh = -67 sig_h = 20 Camph = 7.4 Cbaseh = 1.2 #initialize the vector sizes n = np.zeros(len(I)) m = np.zeros(len(I)) h = np.zeros(len(I)) V = np.zeros(len(I)) Im = np.zeros(len(I)) #set initial conditions: n[0] = 0.3 m[0] = 0.006 h[0] = 0.626 V[0] = Vstart #[mV] #iterate through the input sequence and calculate the response for ii in range(len(I)-1): #display progress (can remove this block if desired) #if (ii % 1000) == 0: #print(ii) #print(V[ii],n[ii],m[ii],h[ii]) #update voltage state variable #DV = 1/C (I - Ik - Ina - Il) #units: #I = uA #V = mV, g = mS; g*V = uA #C = uF #uA/uC*ms = mV #full model V[ii+1] = V[ii] + dt/C*(I[ii] - gl*(V[ii]-El) - gk*n[ii]**4*(V[ii]-Ek) - gna*m[ii]**3*h[ii]*(V[ii]-Ena)) #steady state activation values at the given voltage ninf = 1/( 1 + np.exp( (V50n-V[ii])/kn ) ) minf = 1/( 1 + np.exp( (V50m-V[ii])/km ) ) hinf = 1/( 1 + np.exp( (V50h-V[ii])/kh ) ) #time constants #Izhikevich pg. 45 taun = Cbasen + Campn*np.exp(-(Vmaxn-V[ii])**2/sig_n**2) taum = Cbasem + Campm*np.exp(-(Vmaxm-V[ii])**2/sig_m**2) tauh = Cbaseh + Camph*np.exp(-(Vmaxh-V[ii])**2/sig_h**2) #update activation state variables n[ii+1] = n[ii] + dt*(ninf - n[ii])/taun m[ii+1] = m[ii] + dt*(minf - m[ii])/taum h[ii+1] = h[ii] + dt*(hinf - h[ii])/tauh #define conductances and currents for this v[ii] gkt_ii = gk*n[ii]**4 gnat_ii = gna*m[ii]**3*h[ii] Ik_ii = gkt_ii*(V[ii]-Ek) Ina_ii = gnat_ii * (V[ii] - Ena) Il_ii = gl*(V[ii] - El) #this is wrong lol. Ic_ii = (V[ii+1] - V[ii]) * C / dt #add em up to get total membrane current Im[ii] = Il_ii + Ic_ii + Ina_ii + Ik_ii gkt = gk*n**4 gnat = gna*m**3*h Ik = gkt*(V-Ek) Ina = gnat*(V-Ena) Il = gl*(V-El) #Diagnostics if plot_activation: Vtest = np.arange(-100,11) #steady state activation vs Voltage ninfV = 1/( 1 + np.exp((V50n-Vtest)/kn)) minfV = 1/( 1 + np.exp((V50m-Vtest)/km)) hinfV = 1/( 1 + np.exp((V50h-Vtest)/kh)) line1, = plt.plot(Vtest, ninfV, color = 'b', label = r'$n_{infty}$') line2, = plt.plot(Vtest, minfV, color = 'r', label = r'$m_{infty}$') line3, = plt.plot(Vtest, hinfV, color = 'k', label = r'$h_{infty}$') plt.legend(handles=[line1, line2, line3]) plt.suptitle('steady state activations') plt.show() #time constants vs Voltage taunV = Cbasen + Campn*np.exp(-(Vmaxn-Vtest)**2/sig_n**2) taumV = Cbasem + Campm*np.exp(-(Vmaxm-Vtest)**2/sig_m**2) tauhV = Cbaseh + Camph*np.exp(-(Vmaxh-Vtest)**2/sig_h**2) line1, = plt.plot(Vtest, taunV, color = 'b', label = r'$tau_{n}$') line2, = plt.plot(Vtest, taumV, color = 'r', label = r'$tau_{m}$') line3, = plt.plot(Vtest, tauhV, color = 'k', label = r'$tau_{h}$') plt.suptitle('activation time constants') plt.legend(handles=[line1, line2, line3]) plt.show() return V, Im

inject_time = 60 #inject for 60 ms I_inject = 10 #10 uA total_time = 100 #total time of 100 ms def create_I_list(I_inject, dt, inject_time, total_time): ''' Make current injection of I_inject (mA) for inject_time in intervals of dt. Total_time in ms. ''' I = [I_inject] * int(inject_time / dt) # inject at intervals of dt for inject_time I.extend([0]*int((total_time-inject_time) / dt)) #add on inject current of 0 for the rest of the time after injections return I I_1 = create_I_list(I_inject, 0.1, inject_time, total_time) #dt = 0.1 I_01 = create_I_list(I_inject, 0.01, inject_time, total_time) #dt = 0.01 I_001 = create_I_list(I_inject, 0.001, inject_time, total_time) #dt = 0.001 print(len(I_01))

Vm, Im= HodHux_ELEC380(I_1, 0.1, 0) plt.plot(np.arange(0,100, 0.1), Vm) print(len(Vm)) plt.xlabel('time (ms)') plt.ylabel('Membrane Voltage (mV)') plt.title('10 uA current injection with dt = 0.1') plt.savefig('Q1.2dt0.1.png')

Vm, Im = HodHux_ELEC380(I_01, 0.01, 0) plt.plot(np.arange(0,100, 0.01), Vm) plt.xlabel('time (ms)') plt.ylabel('Membrane Voltage (mV)') plt.title('10 uA current injection with dt = 0.01') plt.savefig('Q1.2dt0.01.png')

Vm, Im = HodHux_ELEC380(I_001, 0.001, 0) plt.plot(np.arange(0,100, 0.001), Vm) plt.xlabel('time (ms)') plt.ylabel('Membrane Voltage (mV)') plt.title('10 uA current injection with dt = 0.001') plt.savefig('Q1.2dt0.001.png')

def when_action_potential(I_inject, dt, total_time): ''' check when the first action potential occurs I_inject: constant input current in uA dt: incremental time counter (s) total_time: total time recorded (s) ''' time = 0 mV = -100 while mV < 0: I_list = create_I_list(I_inject, dt, time, total_time) V_list, Im= HodHux_ELEC380(I_list, dt, plot_activation = 0) #needs to be if any V in this list is above zero, not just the whole list. mV = max(V_list) print("mV:" + str(mV) + "," + "time:" + str(time)) time += dt return time when_action_potential(2, 0.1, 10)

total_time = 60 dt = 0.01 Istim = 2 time = np.arange(0, total_time, dt) I = create_I_list(Istim, dt, when_action_potential(Istim, dt, total_time), total_time + dt) print(len(I)) Vm, Im = HodHux_ELEC380(I, dt, 0) fig, (ax1, ax3, ax4, ax5) = plt.subplots(4, sharex=True) fig.suptitle('Questions 1.3 to 1.6') #plot vm ax1.plot(time, Vm, color="red") ax1.set_ylabel("Vm",color="red",fontsize=14) #plot im on same graph ax2 = ax1.twinx() ax2.plot(time, Im) ax2.set_ylabel("Im", color="blue", fontsize=14) ax3.plot(time, n) ax3.plot(time, m) ax3.plot(time, h) ax3.legend(["n", "m", "h"]) ax4.plot(time, gkt) ax4.plot(time, gnat) ax4.legend(['gK', 'gNa']) ax5.plot(time, Ina) ax5.plot(time, Ik) ax5.legend(['INa', 'IK']) ax5.set_xlabel("time (ms)") #axs[1].plot(x, -y)

currents = np.arange(0,70,0.2) I_and_freq_list = [] def find_firing_rate(driving_current): I_list = create_I_list(driving_current, 0.01, 120, 120) Vm, Im = HodHux_ELEC380(I_list, 0.01, 0) peaks, _ = find_peaks(Vm, height=0) if len(peaks) == 0: frequency = 0 else: frequency = (len(peaks) / 60) * 1000 #convert to Hz return frequency for driving_current in currents: I_and_freq_list.append((driving_current, find_firing_rate(driving_current))) driving_currents = [i[0] for i in I_and_freq_list] frequencies = [i[1] for i in I_and_freq_list] plt.plot(driving_currents, frequencies) plt.xlabel('current (mA)') plt.ylabel('firing frequency (Hz)')

I_A = create_I_list(15, 0.01, 80, 80) time = np.arange(0, 80, 0.01) V_mem_A, I_mem_A = HodHux_ELEC380(I_A, 0.01, plot_activation = 0) I_B = (V_mem_A / 0.15) / 1000 #units are weird here... V_mem_B, I_mem_B = HodHux_ELEC380(I_B, 0.01, plot_activation = 0) plt.plot(time, V_mem_A) plt.plot(time, V_mem_B) plt.title('2.1 Part 1') plt.ylabel('Membrane Potential (mV)') plt.xlabel('time (ms)') plt.legend(['V_mem_A', 'V_mem_B'])

avmin, avmax = min(V_mem_A), max(V_mem_A) for i, val in enumerate(V_mem_A): V_mem_A[i] = (val-avmin) / (avmax-avmin) bvmin, bvmax = min(V_mem_B), max(V_mem_B) for i, val in enumerate(V_mem_B): V_mem_B[i] = (val-bvmin) / (bvmax-bvmin) plt.plot(time, V_mem_A) plt.plot(time, V_mem_B) plt.title('2.1 Part 2') plt.ylabel('Normalized Membrane Potential (mV)') plt.xlabel('time (ms)') plt.legend(['norm_V_mem_A', 'norm_V_mem_B'])

I_A = create_I_list(15, 0.01, 80, 80) time = np.arange(0, 80, 0.01) V_mem_A, I_mem_A = HodHux_ELEC380(I_A, 0.01, plot_activation = 0) old_I_B = (V_mem_A / 0.15) / 1000 #units are weird here... new_I_B = (V_mem_A / 0.3) / 1000 #units are weird here... old_V_mem_B, old_I_mem_B = HodHux_ELEC380(old_I_B, 0.01, plot_activation = 0) new_V_mem_B, new_I_mem_B = HodHux_ELEC380(new_I_B, 0.01, plot_activation = 0) plt.plot(time, old_V_mem_B) plt.plot(time, new_V_mem_B) plt.title('2.3 Proof of Concept') plt.ylabel('Membrane Potential (mV)') plt.xlabel('time (ms)') plt.legend(['old_V_mem_B', 'new_V_mem_B']) #wtf????????