# 26/09/06 /Users/connie/new_ode/ican_2cell2_w_par.ode
# edited parameters to correspond to ican_2cell2synburst.set
# wrote an expression for leak current

# ican_2cell2:  make cabase nonzero, include VGCC in v equation - 7/17/06

#parameters
# h=0
p gsyn=2.2,gton=0.5,el1=-73,el2=-72
p gna=28.,ena=50.0,thetam=-34,sigmam=-5 
p gk=11.2,ek=-85,thetan=-29.0,sigman=-4,tn=4
p gnap=2.8,thetamp=-40.0,sigmamp=-6.0,thetah=-48
p sigmah=6.0,th=10000
p esyn=0.0,iapp=0,c=21, gl1=2.8,gl2=2.8
p k=1,taus=5,sigmas=-5.0,thetas=-10
p thetaca=-39,sigmaca=-8,kca=22.5,eps=5e-05

#functions
minf(v)=1./(1.+exp((v-thetam)/sigmam))
minfp(v)=1./(1.+exp((v-thetamp)/sigmamp))
hinf(v)=1./(1.+exp((v-thetah)/sigmah))
ninf(v)=1./(1.+exp((v-thetan)/sigman))
taun(v)=tn/(cosh((v-thetan)/(2*sigman)))
tauh(v)=th/(cosh((v-thetah)/(2*sigmah)))
sinf(v)=1./(1.+exp((v-thetas)/sigmas))
cainf(v)=1./(1.+exp((v-thetaca)/sigmaca))
#f(s)=1/(1.+exp((v-thetaf)/sigmaf))

#currents
ina(v,n)=gna*(minf(v))^3*(1-n)*(v-ena)
ik(v,n)=gk*n^4*(v-ek)
inap(v,h)=gnap*minfp(v)*h*(v-ena)
iton(v)=gton*(v-esyn)
isyn(v,s)=gsyn*s*v
ican(v,ca)=gcan*ca/(ca+kcan)*(v-vcan)
iahp(v,ca)=gahp*(v-ek)*(ca^hill)/(ca^hill+kahp^hill)
ica(v)=gca1*((cainf(v))^2)*(v-vca)
il(gl,v,el)=gl*(v-el)

aux can1=gcan*ca1/(ca1+kcan)*(v1-vcan)
aux ahp1=gahp*(v1-ek)*(ca1^hill)/(ca1^hill+kahp^hill)
aux can2=gcan*ca2/(ca2+kcan)*(v2-vcan)
aux ahp2=gahp*(v2-ek)*(ca2^hill)/(ca2^hill+kahp^hill)
aux vgcc1=gca*gca1*((cainf(v1))^2)*(v1-vca)
aux vgcc2=gca*gca1*((cainf(v2))^2)*(v2-vca)

p vcan=0,kahp=.7,hill=2,kcan=.3,vca=140.
p gcan=3,gahp=7,gca=0.07,ip3=50,gca1=1,cabase=.1

#initial conditions
i v1=-44.9397839434262  
i n1=0.3814575263697881  
i s1=0
i ca1=0.1322046339284345
i v2=-50
i n2=0.05
i s2=0.09710818195169828
i ca2=0.02

v[1..2]'=-(il(gl[j],v[j],el[j])+ina(v[j],n[j])+ik(v[j],n[j])+ica(v[j])+ican(v[j],ca[j])+iahp(v[j],ca[j])+iton(v[j])+isyn(v[j],s[j])-iapp)/21.0
n[1..2]'=(ninf(v[j])-n[j])/taun(v[j])
s1'=((1-s1)*sinf(v2)-k*s1)/taus
s2'=((1-s2)*sinf(v1)-k*s2)/taus
ca[1..2]'=eps*(-gca*ica(v[j])+ip3*s[j]-kca*(ca[j]-cabase))
#h'=const*(hinf(v)-h)/tauh(v)

#rhs(h,v)=(hinf(v)-h)/tauh(v)
#p const=1

@ dt=.1,total=10000,meth=qualrk,xp=t,yp=v,xlo=0,xhi=10000,ylo=-55,yhi=5,bound=5000,maxstor=1000001

done