#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file Hodgkin_Huxley.m
%Hodgkin-Huxley model 
%reviewed in Rinzel (1990) Bulletin of Mathematical Biology 52 pp. 5-23.
%Figure 1.9 and problem 8.6.4



@author: bingalls
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint



#assign parameter values

E_N=55
E_K=-72
g_N_bar = 120.0
g_K_bar = 36.0
g_leak = 0.30
E_leak = -49.0
C_M = 1.0

#set time_grid for simulation
times=np.linspace(0, 100, 1000) #generate time-grid list


#assign initial condition.  
#state = [V m h n]
Sinit=[-59.8977, 0.0536, 0.5925 , 0.3192]



 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(4) #generate a list to store derivatives
    
    V=S[0]
    m=S[1]
    h=S[2]
    n=S[3]
    
    m_inf = (0.10*(V+35)/(1.0 - np.exp(-(V+35)/10.0)))/((0.10*(V+35)/(1.0 - np.exp(-(V+35)/10.0)))+(4.0*np.exp( -(V+60)/18.0)))
    t_m   = 1.0/((0.10*(V+35)/(1.0 - np.exp(-(V+35)/10.0)))+(4.0*np.exp( -(V+60)/18.0)))

   
    h_inf = (0.07*np.exp(-(V+60)/20))/((0.07*np.exp(-(V+60)/20))+(1/(np.exp(-(V+30)/10)+1)))
    t_h   = 1.0/((0.07*np.exp(-(V+60)/20))+(1/(np.exp(-(V+30)/10)+1)))


    I_N = g_N_bar*(V-E_N)*m**3*h;  

    t_n = 1.0/((0.01*(V+50)/(1.0 - np.exp(-(V+50)/10.0)))+(0.125*np.exp( -(V+60)/80)))
    n_inf = (0.01*(V+50)/(1.0 - np.exp(-(V+50)/10.0)))/((0.01*(V+50)/(1.0 - np.exp(-(V+50)/10.0)))+(0.125*np.exp( -(V+60)/80)))


    I_K = g_K_bar*(V-E_K)*n**4

    I_leak = g_leak*(V-E_leak)
   
  
    Istim=0
    if t>20:
        Istim=-6.65/1.65
    if t>21:
        Istim=0
    if t>60:
        Istim=-6.85
    if t>61:
        Istim=0

  
    
    
    dS[0]= (- I_N - I_K - I_leak - Istim)/C_M
    dS[1]= (m_inf - m)/t_m
    dS[2]= (h_inf - h)/t_h
    dS[3]= (n_inf - n)/t_n
    
    
    
    return dS


S=odeint(dSdt_original, Sinit, times, hmax=0.005) #run simulation

#plot simulation
plt.figure() #generate figure 7.17A
plt.plot(times, S[:,0], label="Voltage", linewidth=2)
plt.xlabel("Time (msec)")
plt.ylabel("membrane voltage (mV)")
plt.legend()