# general parameters
p C=20
p gL=2.8,gsynE=.10,gsynI=.60
p EsynE=-10,EsynI=-75
p EK=-85,ENa=50
p d=1
Eleak1=-54.5*d 
Eleak2=-59*d
Eleak3=-63.5*d
Eleak4=-75*d
p v12n=-30,v12hp=-59,v12mp=-40,v12hk=-92,v12mk=-22
p kn=-4,khp=10,kmp=-6,khk=10,kmk=-6
p tauhp=5000,gNaP=5

# parameters of output function
p lowerlimE=-50,upperlimE=0
p lowerlimI=-50,upperlimI=0
p a=0

# Activity functions
foute(v)=heav(v-lowerlimE)*(v-lowerlimE)/(upperlimE-lowerlimE)+heav(v-upperlimE)*(1-(v-lowerlimE)/(upperlimE-lowerlimE))

# Initial Conditions
init v1=-60,v2=-60,v3=-60,v4=-60
init h1=0,h2=0,h3=0,h4=0,count=0

# functions
ninf(v)=1./(1+exp((v-v12n)/kn))
mpinf(v)=1./(1+exp((v-v12mp)/kmp))
hpinf(v)=1./(1+exp((v-v12hp)/khp))
tauinf(v)=tauhp/cosh((v-v12hp)/(khp*2))

# currents
Inap(v,h)=gNaP*mpinf(v)*h*(v-ENa)
Il1(v)=gL*(v-Eleak1)
Il2(v)=gL*(v-Eleak2)
Il3(v)=gL*(v-Eleak3)
Il4(v)=gL*(v-Eleak4)


v1'=(-Inap(v1,h1)-Il1(v1)-gsynE*(a*foute(v2)+a*foute(v3)+0*a*foute(v4))*(v1-EsynE))/C
v2'=(-Inap(v2,h2)-Il2(v2)-gsynE*(a*foute(v1)+a*foute(v3)+0*a*foute(v4))*(v2-EsynE))/C
v3'=(-Inap(v3,h3)-Il3(v3)-gsynE*(a*foute(v1)+a*foute(v2)+0*a*foute(v4))*(v3-EsynE))/C
v4'=(-Inap(v4,h4)-Il4(v4)-gsynE*(a*foute(v1)+a*foute(v2)+a*foute(v3))*(v4-EsynE))/C


h1'=(hpinf(v1)-h1)/tauinf(v1)
h2'=(hpinf(v2)-h2)/tauinf(v2)
h3'=(hpinf(v3)-h3)/tauinf(v3)
h4'=(hpinf(v4)-h4)/tauinf(v4)
count'=1

aux sum=v1+v2+v3+v4


@ dt=0.1,total=30000,meth=cvode,xp=t,yp=v1,xlow=0,xhi=1000,ylo=-80,yhi=20.,bound=500000,maxstor=10000001

done