# Reproduces approximately Fig. 5, Zhang et al, Biophys. J. 84:2852-2870, 2003.
# Erratum: gkatp=55, not 60 as stated in legend in paper.

v(0)=-54.19088431173661 
n(0)=0.0009122736985496961
c(0)=0.03115585096668747
cer(0)=217.2405145917104
css(0)=0.3782924543159534

# Units: conductances in pS; currents in fA; Ca concentrations in uM; time in ms

# Calcium Handling: cytosol 
par alpha=4.5e-06, kpmca=0.18

# Calcium Handling: ER
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


# Parameters for Currents
# IKv
par vn=-15, sn=5.6, vk=-70, taun=10.8, gkv=2500
# ICa
par vca=30, vm=-13, sm=8, gca=1450
# IKCa
par kd=0.7, q=8, gkca=1200
# ILeak
par gleak=14, vleak=-30
# IKatp
par gkatp=55

# Miscellaneous
par cm=5300

# 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 = gkatp*(v-vk)
IKv = gkv*n*(v-vk)
ILeak = gleak*(v-vleak)

# ER fluxes
Jserca = kserca*c
Jrelease = per*(cer - css)

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

# Equations

v' = (-ICa - IKv - IKatp  - IKCa - ILeak)/cm
n' =  (ninf - n)/taun
c' =   fcyt*(Jx  - kpmca*c - Jserca - alpha*ICa)
cer' = fer*(vcytver*Jserca - Jrelease)
css' = fss*(vervss*Jrelease - vcytvss*Jx)

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=5.0, total=100000, 
@ maxstor=1000000, bounds=10000000, xp=tsec, yp=v
@ xlo=0, xhi=100, ylo=-80, yhi=-10

# Auxiliary quantities to plot:  
aux tsec=t/1000
# average cytosolic Ca2+ (weighted average of css and c):
aux cavg = (vss*css+vcyt*c)/(vss+vcyt)

done