# Equivalent to Fig. 2D, Zhang et al, in preparation.
# (stochastic simulation of a single cell - Cell 902 -with small noise which is a slow
# burster)
# For Fig. 2C (medium burster; Cell 854), change gca to 1465.0
# For these cells, coupling is needed only for synchrony.

v(0)=-21.36600412872696
n(0)=0.141675036335063
c(0)=0.05161719728197939
cer(0)=191.0
css(0)=0.3497928101993529
s(0)=0.1850975132832939


wiener w

# Conductances in pS; currents in fA; Ca concentrations in uM; time in ms

# conductances
par gkatpbar=300
par gkca=1200
# leak
par gleak=14, vleak=-30


# Ca parameters  

par kserca=0.2
par per=0.001

# volume and buffer parameters (*10^-6 ul)
par vcyt=5.0, ver=0.2, vss=2.0
par fcyt=0.01, fer=0.005, fss=0.04
par px=0.025

# Ik
par vn=-15, sn=5.6, vk=-70, taun=20, gk=2500
# Ica
par vca=30, vm=-13, sm=8
par gca=1439.1
# Ikca
par kd=0.7, q=8

# Miscellaneous
par cm=5300 lambda=1.85
# Calcium Handling: cytosol 
par alpha=4.5e-06, kpmca=0.2
# Calcium Handling: ER

# Functions

vervss=ver/vss
vcytvss=vcyt/vss
vcytver=vcyt/ver

ninf = 1/(1+exp((vn-v)/sn))
minf = 1/(1+exp((vm-v)/sm))
omega = 1/(1+(kd/css)^q)

ica = gca*minf*(v-vca)
ikca = gkca*omega*(v-vk)
ikatp = gkatpbar*s*(v-vk)
ik = gk*n*(v-vk)
ileak = gleak*(v-vleak)

# ER fluxes
jerp = kserca*c
jrelease = per*(cer - css)

# Subspace fluxes
jxfer = px*(css-c)

# Stochastics

par alphas=1, taus=100, nkatp=50000, noise=1, s0=0.2
betas = alphas*(1/s0 - 1)

varnoise=sqrt((alphas*(1-s)+betas*s)/(taus*nkatp))

# Equations

v' = (-ica - ik - ikatp  - ikca - ileak)/cm
n' =  lambda*(ninf - n)/taun
c' =   fcyt*(jxfer - kpmca*c - jerp - alpha*ica)
cer' = fer*(vcytver*jerp - jrelease)
css' = fss*(vervss*jrelease - vcytvss*jxfer)
s' = (alphas*(1-s)-betas*s)/taus + noise*w*varnoise

@ meth=euler, toler=1.0e-10, atoler=1.0e-10, dt=0.1, total=200000, 
@ maxstor=1000000, bounds=10000000, xp=tsec, yp=v
@ xlo=0, xhi=200, ylo=-60, yhi=-10, nout=20

aux tsec=t/1000.0
aux cavg=(vss*css+vcyt*c)/(vss+vcyt)

done