#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file morris_lecar_synapse.m
%Morris-Lecar model of synapse onto excitable barnacle muscle fiber
%adapted from Morris and Lecar (1981) Biophysical Journal 35 pp. 193-231
%Figure 8.9


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



#assign parameter values

#mossi-lecar model
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

#synapse activity
gsyn=1 # millisiemens/ cm^2 
Esyn= 150 # millisiemens/ cm^2 
alpha=1 # /ms
beta=0.12 #/ms
Tmax=1
Vsynstar=27 # mV

Iapplied=0;



#assign initial condition.  
#state is Vpre wpre ssyn Vpost wpost
#rest:
S00=[-60.8554, 0.0149, 0, -60.8554, 0.0149]
#super-threshold
S10=[-10, 0.0149, 0, -60.8554, 0.0149]


 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(5) #generate a list to store derivatives
    
    Vpre=S[0]
    wpre=S[1]
    ssyn=S[2]
    Vpost=S[3]
    wpost=S[4]
    
    #local functions:
    mpre_inf = (0.5*(1+np.tanh((Vpre-v1)/v2)));
    wpre_inf = (0.5*(1+np.tanh((Vpre-v3)/v4)));

    mpost_inf = (0.5*(1+np.tanh((Vpost-v1)/v2)));
    wpost_inf = (0.5*(1+np.tanh((Vpost-v3)/v4)));

    if Vpre>=Vsynstar:
        T=Tmax
    else:
        T=0

    dS[0]= (1/C)*(gbarCa*mpre_inf*(ECa-Vpre) + gbarK*wpre*(EK-Vpre) + gleak*(Eleak-Vpre));
    dS[1]= phi*(wpre_inf-wpre)/(tau);
    dS[2]= alpha*T*(1-ssyn) - beta*ssyn;
    dS[3]= (1/C)*(gbarCa*mpost_inf*(ECa-Vpost) + gbarK*wpost*(EK-Vpost) + gleak*(Eleak-Vpost)+gsyn*ssyn*(Esyn-Vpost));
    dS[4]= phi*(wpost_inf-wpost)/(tau);
    
    
    return dS

#run simulation

#set time_grid for simulation
times1=np.linspace(-10, 100, 1000) #generate time-grid list
times2=np.linspace(100, 300, 1000)

S1=odeint(dSdt_original, S00, times1) #run simulation
S2=odeint(dSdt_original, S10, times2)

times=np.concatenate((times1, times2))
S=np.concatenate((S1, S2))

#plot simulation
plt.figure() #generate Figure 8.9A
plt.plot(times, S[:,0], label="$V_{pre}$", linewidth=2)
plt.plot(times, S[:,3], label="$V_{post}$", linewidth=2)
plt.xlabel("Time (msec)")
plt.ylabel("voltage (mV)")
plt.legend()

#plot simulation
plt.figure() #generate Figure 8.9B
plt.plot(times, S[:,1], label="$I_{syn}$", linewidth=2)
plt.xlabel("Time (msec)")
plt.ylabel("Current $(pA/cm^2)$")
plt.legend()