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)
1.8
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????????