CoCalc Shared FilesHodgkinsHuxleyManipulation.sagews
Author: Aarushi Solanki
Views : 13
from scipy.integrate import solve_ivp

#Constants
C_m = 0.01
E_Na = 55.17
E_K = -72.14
E_l = -49.42
g_Na = 1.2
g_K = 0.36
g_l = 0.003

#Rate Functions for Potassium Ion Channel

alpha_n(v) = (0.01 * (v+50)) / (1 - e^(-(v+50)/10))
beta_n(v) = 0.125 * e^(-(v+60) / 80 )

#Rate Functions for Sodium Ion Channel

alpha_m(v) = (0.1 * (v+35)) / (1 - e^( -(v+35) / 10))
beta_m(v) = 4.0 * e^(-0.0556 * (v+60))
alpha_h(v) = 0.07 * e^(-0.05 * (v+60))
beta_h(v) = 1 / (1 + e^(-0.1 * (v+30)))

#Initial values for n, m, and h
nnaught = alpha_n(-60) / (alpha_n(-60) + beta_n(-60))
mnaught = alpha_m(-60) / (alpha_m(-60) + beta_m(-60))
hnaught = alpha_h(-60) / (alpha_h(-60) + beta_h(-60))

@interact
def neuron(vnaught=(-60,60), Inaught=(0,1,0.01)):
def diffeqs(t, state):
n,m,h,v = state
nprime = (alpha_n(v) * (1 - n)) - (beta_n(v) * n)
mprime = (alpha_m(v) * (1 - m)) - (beta_m(v) * m)
hprime = (alpha_h(v) * (1 - h)) - (beta_h(v) * h)
vprime = (1/C_m) * (Inaught - (g_Na * m^3 * h * (v - E_Na)) - (g_K * n^4 * (v - E_K)) - (g_l * (v - E_l)))
return (nprime, mprime, hprime, vprime)
t_span = (0,100)
solution = solve_ivp(diffeqs, t_span, (nnaught, mnaught, hnaught, vnaught), method="LSODA")
Vsol = solution.y[3]
p = list_plot(zip(solution.t,Vsol), axes_labels=["t", "Voltage (mV)"], plotjoined=True, ymax=100)
show(p)



#when I = 0, action potential fires when V>=-53.4
#When I = 1, v=60 actually most, periodic

# btwn 0.05 and 0.06 theres the hopf bifurcation
#spiral inward to spiral outward --> hopf bifurcation