###########################################################################
# This model has a delayed rectifier, a leak and
# a sodium current based on the famous Traub and Miles model but with changes put in
# by Bard Ermentrout in linearization of F-I curves by adaptation, Neural computation, 
# 10:1721-9, 1998
# 
# Input Currents are in nA, maximal conductances in mS/cm2, membrane potential 
# in mV time in ms, and capacitance in uF/cm2
############################################################################


#adjust defaults for xpp internal variables
@ MAXSTOR=40000
@ YP=v
@ TOTAL=1000
@ DT=0.05
@ BOUND=100000
@ XHI=1000
@ YLO=-100
@ YHI=50
#@ METH=gear
@ BACKGROUND=white

# parameters


par gbarNa=80


par gbarK=80

par gleak=0.2

par C=1., i=1


v' = (i-(gleak*(v-Vl)+gK*(v-vk)+gNa*(v-vna))+eps*(u1-v)/dx)/C+p0*pulse(t-taup)

aux longcur=eps*(u1-v)/dx
#****************************************************
gNa=gbarNa*(m**3)*h

alpham=if (v+54) then (0.32*(V+54)/(1-exp(-(V+54)/4))) else (0.32)
betam=if (V+27) then (0.28*(V+27)/(exp((V+27)/5)-1)) else (0.28)

alphah=0.128*exp(-(V+50)/18)
betah=4/(1+exp(-(V+27)/5))

dm/dt=alpham*(1-m)-betam*m
dh/dt=alphah*(1-h)-betah*h

#**********************************************************

gK=gbarK*(n**4)

alphan=if (V+52) then (0.032*(V+52)/(1-exp(-(V+52)/5))) else (0.032)
betan=0.5*exp(-(V+57)/40)

dn/dt=alphan*(1-n)-betan*n



# auxiliary variables
# sod is sodium current
aux sod=gbarNa*(m**3)*h*(v-vna)
# pot is potassium current
aux pot=gbarK*(n**4)*(v-vk)

#initial conditions
init v=-70
init h=1-.000001,m=.000001,n=.000001

par Vna=50,Vk=-100,Vl=-67,Vh=-65,Vrh=-40,Vsh=6,gld=0.1,gh=0.05,tauh=400



hinf(V)=1/(1+exp((V-Vh)/Vsh))
Ild(V)=V-Vl
Ih(V,y)=gh*y*(V-Vrh)/gld

h[1..50]'=(hinf(u[j])-h[j])/tauh
# NOT TO CONFUSE WITH h GATE IN SOMA!!

# cable equation has HCN current

p p[0..50]=0


u1'=((lambda/dx)^2*(u2-2*u1+v)-Ild(u1)-Ih(u1,h1))/tau +p1*pulse(t-taup)
u[2..50]'= ((lambda/dx)^2*(u[j+1]-2*u[j]+u[j-1])-Ild(u[j])-Ih(u[j],h[j]))/tau +p[j]*pulse(t-taup)
u51=(c1+b1*u50/dx)/(a1+b1/dx)

par lambda=1,tau=10,dx=.1,c1=0,a1=0,b1=1,c0=0,a0=0,b0=1,eps=.025

pulse(t)=heav(t)*heav(sigma-t)
par sigma=.2
par t0=241.3
aux prc=t0-t
taup'=0

i u[1..50]=-60
i h[1..50]=.1

par taur=1,taud=3,thresh=-30, gsyn=.1 Esyn=0
x1'=(-x1+.5*(1+tanh((v-thresh)/3.0)))/taur
y1'=(-y1+x1)/taud
init x1=0,y1=0

d