#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file morris_lecar.py
%Morris-Lecar model of excitable barnacle muscle fiber
%adapted from Morris and Lecar (1981) Biophysical Journal 35 pp. 193-231
%Figures 8.6 and 8.7


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



#assign parameter values

C= 20  #microfarad/cm^2 
gbarCa=4.4 # millisiemens/ cm^2 
ECa=120 #millivolts
gbarK=8 # millisiemens/ cm^2 
EK=-84 #%millivolts
gleak=2 # millisiemens/ cm^2 
Eleak=-60 #millivolts
v1=-1.2 #millivolts
v2= 18 #millivolts
v3= 2 #millivolts
v4= 30 #millivolts
phi= 0.04 # per millisecond
tau=0.8

Iapplied=0;


#set time_grid for simulation
times1=np.linspace(0, 100, 1000) #generate time-grid list
times2=np.linspace(100, 250, 1000) 
times3=np.linspace(250, 400, 1000) 
times4=np.linspace(400, 550, 1000) 
times5=np.linspace(550, 700, 1000) 


#assign initial conditions 
#state vector S=[V, w]
#rest:
Sinit1=[-60.8554, 0.0149]
#sub-threshold
Sinit2=[-35, 0.0149]
#sub-threshold
Sinit3=[-15, 0.0149]
#super-threshold
Sinit4=[-13, 0.0149]
#super-threshold
Sinit5=[15, 0.0149]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(2) #generate a list to store derivatives
    
    V=S[0]
    w=S[1]
    
    #local functions:
    m_inf = (0.5*(1+np.tanh((V-v1)/v2)))
    w_inf = (0.5*(1+np.tanh((V-v3)/v4)))
   
    dS[0]=(1/C)*(gbarCa*m_inf*(ECa-V) + gbarK*w*(EK-V) + gleak*(Eleak-V)+Iapplied)
    dS[1]=phi*(w_inf-w)/(tau)
   
    return dS


S1=odeint(dSdt_original, Sinit1, times1) #run simulation
S2=odeint(dSdt_original, Sinit2, times2)
S3=odeint(dSdt_original, Sinit3, times3)
S4=odeint(dSdt_original, Sinit4, times4)
S5=odeint(dSdt_original, Sinit5, times5)

times=np.concatenate((times1, times2, times3, times4, times5))
S=np.concatenate((S1, S2, S3, S4, S5))

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


#next, figure 8.6B

times=np.linspace(0,300,1000)

Sp1=odeint(dSdt_original, [-35, 0.0149], times)
Sp2=odeint(dSdt_original, [-15, 0.0149], times)
Sp3=odeint(dSdt_original, [-13, 0.0149], times)
Sp4=odeint(dSdt_original, [5, 0.0149], times)


plt.figure() #generate figure 8.6B
plt.plot(Sp1[:,0], Sp1[:,1],  linewidth=2)
plt.plot(Sp2[:,0], Sp2[:,1],  linewidth=2)
plt.plot(Sp3[:,0], Sp3[:,1],  linewidth=2)
plt.plot(Sp4[:,0], Sp4[:,1],  linewidth=2)
plt.xlabel("Potassium gating variable (w)")
plt.xlabel("membrane voltage V (mV)")

#generate nullclines

delta = 0.025
# xrange = arange(0, 2, delta)
# yrange = arange(0, 2, delta)
xrange = np.arange(-80, 40, delta)
yrange = np.arange(-0.15, 0.5, delta)
Vn, wn = np.meshgrid(xrange,yrange)


F = (1/C)*(gbarCa*(0.5*(1+np.tanh((Vn-v1)/v2)))*(ECa-Vn) + gbarK*wn*(EK-Vn) + gleak*(Eleak-Vn))
G = phi*((0.5*(1+np.tanh((Vn-v3)/v4)))-wn)/tau

#plot nullclines
plt.contour(Vn, wn, F, [0])
plt.contour(Vn, wn, G, [0])


####figure 8.7A####

#initial condition: subthreshold perturbation
Sinita=[-35, 0.0149];

timesa=np.linspace(0, 200, 1000)

Iapplied=20
Sa1=odeint(dSdt_original, Sinita, times)
Iapplied=150
Sa2=odeint(dSdt_original, Sinita, times)
Iapplied=400
Sa3=odeint(dSdt_original, Sinita, times)

#plot simulation
plt.figure() #generate figure 8.7A
plt.plot(timesa, Sa1[:,0], label="$I_{applied} = 20$", linewidth=2)
plt.plot(timesa, Sa2[:,0], label="$I_{applied} = 150$", linewidth=2)
plt.plot(timesa, Sa3[:,0], label="$I_{applied} = 400$", linewidth=2)
plt.xlabel("Time (msec)")
plt.ylabel("membrane voltage (mV)")
plt.legend()

####figure 8.7B####

Iapplied=150

timesb=np.linspace(0, 300, 1000)

Sp1=odeint(dSdt_original, [-40, 0], timesb)
Sp2=odeint(dSdt_original, [60, 0.4], timesb)
Sp3=odeint(dSdt_original, [20, 0.75], timesb)
Sp4=odeint(dSdt_original, [-70, 0.5], timesb)
Sp5=odeint(dSdt_original, [0, 0.35], timesb)
Sp6=odeint(dSdt_original, [-10, 0.45], timesb)

#plot simulation
plt.figure() #generate figure 8.7A
plt.plot(Sp1[:,0], Sp1[:,1],linewidth=2)
plt.plot(Sp2[:,0], Sp2[:,1],linewidth=2)
plt.plot(Sp3[:,0], Sp3[:,1], linewidth=2)
plt.plot(Sp4[:,0], Sp4[:,1], linewidth=2)
plt.plot(Sp5[:,0], Sp5[:,1], linewidth=2)
plt.plot(Sp6[:,0], Sp6[:,1], linewidth=2)
plt.ylabel("Potassium gating variable (w)")
plt.xlabel("membrane voltage V (mV)")
plt.legend()

#generate nullclines

delta = 0.025
# xrange = arange(0, 2, delta)
# yrange = arange(0, 2, delta)
xrange = np.arange(-70, 60, delta)
yrange = np.arange(0, 0.8, delta)
Vn, wn = np.meshgrid(xrange,yrange)



F = (1/C)*(gbarCa*(0.5*(1+np.tanh((Vn-v1)/v2)))*(ECa-Vn) + gbarK*wn*(EK-Vn) + gleak*(Eleak-Vn)+Iapplied);
G = phi*((0.5*(1+np.tanh((Vn-v3)/v4)))-wn)/tau;

#plot nullclines
plt.contour(Vn, wn, F, [0])
plt.contour(Vn, wn, G, [0])