# Model of somato-dendritic burster
# From Toporikova and Butera, J. Comp neuroscience, DOI 10.1007/s10827-010-0274-z
# Model with following parameters reproduce Figure 1B of the paper above
# setting gcan=0, convert model to Pruvis and Butera 2007
# units: V = mV; Kca = uM; Cm = uF/cm^2; g = uS/cm^2; phi = 1/s;
#Parameter fi is normalized with respect to intracellular volume, so it's units fi=1/pL; 

# PARAMETERS

# conductance of CaN current
par gcan=1.5

# conductances of persistent sodium current
par gnaps=2

#IP3 concentration
par I=1

#Total Ca
par Ct=1.25

# Somatic Parameters (same as Butera and Pruvis 2007)
num Cms=21
nun vna=50,Vk=-85
num gk=11.2, gna=28
num vm=-34, vn=-29, vmp=-40,vh=-48
num sm=-5, sl=6, sn=-4, smp=-6, sh=5
num taunb=10, tauhb=10000
par gls=2.3
par vleaks=-57
par Iaps=0

# The ER parameters
par fi=0.0001	
num Vi=4
num LL=0.37
num P=31000
num Ki=1.0
num Ka=0.4
num Ve=400
num Ke=0.2
num A=0.005 
num Kd=0.4
num sigma=0.185

# Other parameters 
num Cmd=5, gc=1,k=0.3
num Kcan=0.74, ncan=0.97
  

# SOMATIC FUNCTIONS

minfs=1/(1+exp((vs-vm) /sm))
ninfs=1/(1+exp((vs-vn) /sn))
minfps=1/(1+exp((vs-vmp)/smp))
hinfs =1/(1+exp((vs-vh) /sh))

tauns=taunb/cosh((vs-vn)/(2*sn))
tauhs=tauhb/cosh((vs-vh)/(2*sh))

I_nas=gna*minfs^3*(1-ns)*(vs-vna)
I_ks=gk*ns^4*(Vs-Vk)
I_naps=gnaps*minfps*hs*(vs-vna)

i_ls =gls*(vs-vleaks)
I_ds =gc*(Vd-Vs)/k

#DENDRITIC FUNCTIONS

J_ER_in=(LL + P*( (I*C*l)/( (I+Ki)*(C+Ka) ) )^3 )*(Ce - C)
J_ER_out=Ve*C^2/(Ke^2+C^2)

caninf =1/(1+(Kcan/C)^ncan)
I_can=gcan*caninf*(vd-vna)
I_sd =gc*(Vs-Vd)/(1-k)

# SOMATIC EQUATIONS

Vs'= (-I_ks - I_nas-I_naps-I_ls-Iaps-I_sd)/Cms
ns'= (ninfs-ns)/tauns
hs'= (hinfs-hs)/tauhs

# DENDRITIC EQUATIONS

Vd'= (-I_can-I_ds)/Cmd
C' = fi/Vi*( J_ER_in- J_ER_out)
l' = A*( Kd - (C + Kd)*l )

#Calculate ER Ca

Ce = (Ct - C)/sigma

# INITIAL CONDITIONS

Vs(0)=-50
Vd(0)=-48

ns(0)=0.004
hs(0)=0.33
C(0)=0.03
l(0)=0.93

aux Ce=Ce
aux ican=I_can
aux inaps=I_naps

@ TOTAL=10000,DT=.21,xlo=0,xhi=10000,ylo=-60,yhi=10
@ xplot=t,yplot=vs, dt=0.11
@ dsmin=1e-5,dsmax=.1,parmin=-.5,parmax=.5,autoxmin=-.5,autoxmax=.5
@ autoymax=.4,autoymin=-.5
set vvst {xplot=c,yplot=h,xlo=0,xhi=1.5,ylo=0,yhi=1.5,total=100 \
	dt=0.01,meth=qualrk}
done