# Reproduces Fig. 10, Zhang et al, Biophys. J. 84(5):2852:2870
# (stochastic simulation of dynamic-clamp conversion of a noisy single cell
# to a medium burster)

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


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


# Dclamp params:
# Gmaxon in nS, K in sec for compatibility with Dclamp
# (converted to pS and msec below):
par Gmaxon=0.01, K=2, vp=-30, sp=7.5, vr=100
par tonD=20000, toffD=800000
 
# set Dclamp conductance (stepped on at ton, off at toff):
Gmax = Gmaxon*(heav(t - tonD) - heav(t - toffD))
zinf = 1/(1.0+exp((vp-v)/sp)) 

# 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=1450
# 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)
iclamp = Gmax*1000*z*(v-vr)

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

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

# Stochastics

par alphas=1, taus=100, nkatp=500, 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 - iclamp)/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
z' = K*0.001*(zinf - z) 

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

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

done