# traub sodium and potassium kinetics and it
# cell 1
v1'=-(gna*h1*m1^3*(v1-ena)+gk*n1^4*(v1-ek)+gl*(v1-el) -i1+gsyn2*s2*(v1-vsyn2))/c
m1'=am(v1)*(1-m1)-bm(v1)*m1
h1'=ah(v1)*(1-h1)-bh(v1)*h1
n1'=an(v1)*(1-n1)-bn(v1)*n1
s1'=alpha1*tmax*(1-s1)/(1+exp(-(v1-vt)/vs))-beta1*s1
# cell 2
v2'=-(gna*h2*m2^3*(v2-ena)+gk*n2^4*(v2-ek)+gl*(v2-el) -i2+gsyn1*s1*(v2-vsyn1))/c
m2'=am(v2)*(1-m2)-bm(v2)*m2
h2'=ah(v2)*(1-h2)-bh(v2)*h2
n2'=an(v2)*(1-n2)-bn(v2)*n2
s2'=alpha2*tmax*(1-s2)/(1+exp(-(v2-vt)/vs))-beta2*s2
am(v)=.32*(54+v)/(1-exp(-(v+54)/4))
bm(v)=.28*(v+27)/(exp((v+27)/5)-1)
ah(v)=.128*exp(-(50+v)/18)
bh(v)=4/(1+exp(-(v+27)/5))
an(v)=.032*(v+52)/(1-exp(-(v+52)/5))
bn(v)=.5*exp(-(57+v)/40)
par ek=-100,ena=50,el=-67
par gl=.1,gk=80,gna=100
par c=1,i1=0,i2=0,
par gsyn1=0,gsyn2=0,vsyn1=0,vsyn2=0
par vt=2,vs=5,tmax=3.2
par alpha1=1,beta1=.2,alpha2=1,beta2=.2
init v1=-67,v2=-67,h1=1,h2=1
@ dt=.25,meth=qualrk,total=100,xhi=100,ylo=-85,yhi=25
done