#
# WBrealVmD_K.ode
# WB interneuron with K modulation
# 
# Author: Ernest Ho 
# Date: January 31, 2009

#
# Defining the equations
#


f(_V,_N,_H,_J,_L,_Ko,_Cai)=ILEAKCl(_V)+ILEAKNa(_V)+ILEAKK(_V,_Ko)+INA(_V,_H)+IKDR(_V,_N,_Ko)+IAHP(_V,_Ko,_Cai)+ICA(_V)+INAP(_V,_J)+IKM(_V,_L,_Ko)

ILEAKK(_V,_Ko)=gLK*(_V-EK(_Ko))
ILEAKCl(_V)=gLCl*(_V-ECl)
ILEAKNa(_V)=gLNa*(_V-ENa)

INA(_V,_H)=gNa*(minf(_V)**3)*(_H)*(_V-ENa)
INAP(_V,_J)=gNaP*(_J)*(_V-ENa)


minf(_V)=alpham(_V)/(alpham(_V)+betam(_V))
alpham(_V)=-0.1*(_V+35)/(exp(-0.1*(_V+35))-1)
betam(_V)=4*exp(-(_V+60)/18)
		

IKDR(_V,_N,_Ko)=gK*(_N**4)*(_V-EK(_Ko))
IKM(_V,_L,_Ko)=gKm*(_L)*(_V-EK(_Ko))


Ipump(_Ko)=Ipumpmax/(1+K_oeq/_Ko)**2
r_forward(_Ko)=r_f/(1+exp((_Ko-K_oth)/(rise)))
glia_uptak(_Ko,_B)=r_back*(Bmax-_B)/1.1-r_forward(_Ko)*(_Ko)*(_B)


bound_buff(_Ko,_B)=r_back*(Bmax-_B)-r_forward(_Ko)*(_Ko)*(_B)
Cai_func(_V,_Cai)=-0.002*gCa*(_V-ECa)/(1+exp(-(_V+25)/2.5))-_Cai/80

IAHP(_V,_Ko,_Cai)=gAHP*(V-EK(_Ko))*(_Cai)/(1+_Cai)
ICA(_V)=gCa*(_V-ECa)/(1+exp(-(_V+25)/2.5))


EK(_Ko)=26.64*ln(_Ko/133)




func_n(_V,_N)=5*(ninf(_V)-_N)/taun(_V)
ninf(_V)=alphan(_V)/(alphan(_V)+betan(_V))
taun(_V)=1/(alphan(_V)+betan(_V))
alphan(_V)=-0.01*(_V+34)/(exp(-0.1*(_V+34))-1)
betan(_V)=tnzero*exp(-(_V+44)/80)

func_h(_V,_H)=5*(alphah(_V)*(1-_H)-betah(_V)*_H)
alphah(_V)=0.07*exp(-(_V+58)/20)
betah(_V)=1/(exp(-0.1*(_V+28))+1)
gs_inf(_V)=gs_max*tanh((abs(_V-gs_vth)+(_V-gs_vth))/(2.0*gs_vslope))


func_j(_V,_J)=5.0*(jinf(_V)-_J)
jinf(_V)=1/(1+exp(-(_V+42)/5.0))


func_l(_V,_L)=(linf(_V)-_L)/taul(_V)
linf(_V)=alphal(_v)/(alphal(_V)+betal(_V))
taul(_V)=0.34/(alphal(_V)+betal(_V))
alphal(_V)=0.001*(_V+30)/(1-exp(-(_V+30)/9.0))
betal(_V)=-0.001*(_V+30)/(1-exp((_V+30)/9.0))

#
# Poisson spike train
possrate(_rate) = _rate*[DT]

init V=-65
init n=0.1
init h=0.05
init j=0.05
init l=0.05
init Ko=[KOINIT]
init B=[BINIT]
init Cai=0
init s=0
init se1=0
init se2=0
init si1=0
init si2=0
init gs=0
# weiner noise

V'=-ge1zero*se1*(V-Eexc)-ge2zero*se2*(V-Eexc)-gi1zero*si1*(V-Ein)-gi2zero*si2*(V-Ein)-f(V,n,h,j,l,Ko,Cai)+iext-frac*Ipump(Ko)+ipulse*(heav(t-pstart)-heav(t-pend)) 

# last term is a 3ms pulse

# ge'=-(ge-gezero)/tau_e+sqrt(2*sigmae^2/tau_e)*noise
n'=func_n(V,n)
h'=func_h(V,h)
j'=func_j(V,j)
l'=func_l(V,l)
Ko'=(50/96485.0)*(ILEAKK(V,Ko)+IKDR(V,n,Ko)+IKM(V,l,Ko)+IAHP(V,Ko,Cai)-Ipump(Ko))+glia_uptak(Ko,B)
gs' = (gs_inf(V) - gs)/((gs_max - gs_inf(V))*tau_gs)
s' = -s/tau_s
global 1 V {s=s+s_max}
se1' = -se1/tau_e
si1' = -si1/tau_i
se2' = -se2/tau_e
si2' = -si2/tau_i


global 1 ((ran(1)<possrate(rate_e1))-0.5) {se1=se1+1.0}
global 1 ((ran(1)<possrate(rate_e2))-0.5) {se2=se2+1.0}
global 1 ((ran(1)<possrate(rate_i1))-0.5) {si1=si1+1.0}
global 1 ((ran(1)<possrate(rate_i2))-0.5) {si2=si2+1.0}




Cai'=Cai_func(V,Cai)
B'=bound_buff(Ko,B)
aux EKrev=EK(Ko)
#aux ELeak=EL(Ko)

aux glia_curr=glia_uptak(Ko,B)
aux pota_curr=IKDR(V,N,Ko)
aux pump_curr=Ipump(Ko)
aux KB=Bmax-B
aux KI=ILEAKK(V,Ko)+IKDR(V,n,Ko)-Ipump(Ko)
aux IKCA=IAHP(V,Ko,Cai)
aux m_curr=IKM(V,l,Ko)
aux nap_curr=INAP(V,j)
aux ca_curr=ICA(V)
aux DEL=delay(t,1)

# Initial conditions

par s_max=1.0
par gs_max = 1.0
par tau_s = [TAUS] 
#3.0
par tau_gs = [TAUGS] 
#3.0
par gs_vth = -35.0
par gs_vslope=10.0
par tau_e=3
par tau_i=6
par tnzero = 0.125
par gLK=0.035
par gLCl=0.1
par gLNa=0.025
par gK=[GK]
par gKm=0.00
par gNa=35
par gNaP=0.0
par gAHP=[GAHP] 
# 5 or 6
# gAHP=0.8 (IN)
par gCa=0.03
par ENa=70
par ECl = -74
par ECa = 120
par Eexc = 0
par Ein = -73
par ge1zero=0.00
par ge2zero=0.00
par gi1zero=0.00
par gi2zero=0.00
par rate_e1=0.5
par rate_e2=0.0
par rate_i1=0.0
par rate_i2=0.5
# par sigmae=0.00
par r_back=0.0008
par r_f=0.0008
par rise=[SLOPE] 
#-1.15
par Bmax=500
par K_oeq=[KOEQ] 
#3
par K_oth=[KOTHETA]
#15
par iext=[IEXT]
#0
par Ipumpmax=[IPUMPMAX]
#7
par frac=[IPUMPFRAC]
#0.5
par ipulse=0
par pstart=90000
par pend=90005

@ total=[DURATION],dt=[DT]
@ meth=stiff
#@ meth=euler
@ ylo=-90,yhi=60
@ xlo=0, xhi=1000
@ maxstor=200000000
@ bound=100000000
@ autoeval=0
done