4302.sagews

from scipy.integrate.odepack import odeint
import pylab
import numpy as np


# Roman, everything is in SI, except for  I convert to mV to 
# find alpha/beta rate constants, since this is what neuron uses.

# Surface Area:
sa = 1e-5 # in cm2

# Capacitance (1uF/cm2)
C = 1e-6 * sa

# Reversal Potentials:
erev_lk = -51.e-3
erev_k = -81.e-3
erev_na = 50.e-3

# All conductances in mS/cm2 
glk = 0.25e-3 * sa * 1.2
gna = 10e-3 *  sa 
gkf = 2.5e-3 * sa * 1.0
gks = 2.0e-3 * sa * 1.5


# Alpha-Beta Coefficients:
h_alpha = [1.64, 0.00, 0.49, 67.30, 15.91 ] 
h_beta  = [0.45, 0.00, 0.12, 0.08, -10.08 ]

m_alpha = [ 7.95, 0.00, 0.93, 0.04, -11.66, ]
m_beta =  [ 1.64, 0.00, 0.43, -0.03, 9.03]

kf_alpha = [5.06, 0.07, 5.12, -18.40, -25.42]
kf_beta =  [0.50, 0.00, 0.00, 28.69, 34.62]

ks_alpha = [0.46, 0.01, 4.59, -4.21, -11.97]
ks_beta  = [0.09, -0.00, 1.62, 210656.27, 332762.27]



# Form from Sautois et al 07
def AlphaBetaStdEqn(v,A):
    a,b,c,d,e = A
    return (a+b*v)/(c+np.exp((v+d)/e) )
                    
def getInfTau(V, alpha_coeffs,beta_coeffs):
    v_mv = V * 1000.
    alpha = AlphaBetaStdEqn(v_mv, alpha_coeffs )
    beta = AlphaBetaStdEqn(v_mv, beta_coeffs )
    
    inf = alpha / (alpha + beta)
    tau = 1.0 / (alpha + beta) * 1e-3 # in ms
    
    return inf,tau


def stdDerivative( V, state,  alpha_coeffs, beta_coeffs):
    inf,tau = getInfTau(V, alpha_coeffs=alpha_coeffs, beta_coeffs=beta_coeffs)
    
    return (inf-state)/tau



        
    


def ode_func(y,t):
    V, m, h, kf, ks = y
    
    # Injected Current: 
    # 120pA, with a 'rebound pulse of 
    if 120e-3 < t < 140e-3:
        iInj = 70e-12
    elif t > 100e-3:
        iInj = 120e-12
    else:
        iInj = 0
        
    
    
    
    iLk = glk* (erev_lk - V ) 
    iNa = gna* (erev_na - V ) * m*m*m*h  
    iKf = gkf* (erev_k - V ) * kf*kf*kf*kf
    iKs = gks* (erev_k - V ) * ks*ks
    
    i = iLk+iNa + iKf + iKs
    
    dV = 1/C * ( i + iInj )
    
    #dV = 0
    dm =  stdDerivative(V, m, m_alpha, m_beta )
    dh =  stdDerivative(V, h, h_alpha, h_beta )
    dkf = stdDerivative(V, kf, kf_alpha, kf_beta )
    dks = stdDerivative(V, ks, ks_alpha, ks_beta )
    

    
    d = [dV, dm, dh, dkf, dks]
    return d
    


    
y0 =[ -51e-3, 0, 0, 0 ,0,]    
t = np.arange(0, 0.300, 0.1 * 1e-3)
res =odeint( ode_func, y0,t )

 
v = res[:,0]
m = res[:,1]
h = res[:,2]
kf = res[:,3]
ks = res[:,4]


f = pylab.figure()
ax1 = f.add_subplot(2,1,1)
ax2 = f.add_subplot(2,1,2)

ax1.plot(t, v )
ax2.plot(t, m, label="m" )
ax2.plot(t, h, label="h" )
ax2.plot(t, ks, label="kf" )
ax2.plot(t, kf, label="ks" )
ax2.legend()

f.savefig('.')

pylab.show()