import scipy as sp
import pylab as plt
from scipy.integrate import odeint
from scipy import stats
import scipy.linalg as lin

## Hodgkin-Huxley Model For L-Type Voltage-dependent Ca2+ Channel EGL-19
## Jospin M. et.al V.S. Boyle and Cohen

# Constants
g_Ca1 = 220.0 # maximum conducances, in mS/cm^2
g_Ca = 199
#g_Ca = 127
a_Ca = 0.283
E_Ca1 = 49.1 # Nernst reversal potentials, in mV
E_Ca = 50.0
#E_Ca = 59.0
E_half_Ca = 0.9
#E_half_Ca = 6.4
E_half_e = -3.4
E_half_f = 25.2 #5.2 is much better
Ca_half_h = 64.1 * sp.power(10,-9)
k_Ca = 4.9 #mV
#k_Ca = 7.9 #mV
k_e = 6.7
k_f = 5.0 #
k_h = -0.01 #mM
Ca_Con = 2.39 * sp.power(10,-6)


# Channel currents (in mA/cm^2)
def I_Ca_boyle(V): return g_Ca1 * e(V)**2 * f(V) * (1 + (h - 1) * a_Ca) * (V - E_Ca1) / 1000
def I_Ca(V): return g_Ca * (V - E_Ca) / ((1 + sp.exp((E_half_Ca - V) / k_Ca)) * 1000)

# The voltage to integrate over
V = sp.arange(-70, 80, 1)

# Gating Kinetics
def X_inf(V, V_half_x, k_x): return 1 / (1 + sp.exp((V_half_x - V) / k_x))

def e(V): return X_inf(V, E_half_e, k_e)
def f(V): return X_inf(V, E_half_f, k_f)

h = X_inf(Ca_Con, Ca_half_h, k_h)

icab = I_Ca_boyle(V)
ica = I_Ca(V)

plt.figure()

plt.subplot(1,1,1)
plt.plot(V, ica, 'r', label='$I_{Ca}$')
plt.plot(V, icab, 'y', label='$I_{Ca}-Boyle$')
plt.ylabel('I (A/F)')
plt.xlabel('Em (mV)')
plt.legend()

plt.show()