HodgkinsHuxleyManipulation.sagews

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)
    

    
    

Interact: please open in CoCalc

#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