# mean_Cer.ode
 
# Published in Bulletin of Mathematical Biology, 70:1251-1271, 2008.
# This file contains the calcium subsystem only, and voltage is clamped
# to produce a train of depolarizations. An analytical expression for the
# quasi-equilibrium cytosolic calcium (Ceq) is used. There is an ode for the 
# ER calcium concentration (Cer), as well as an expression for the mean ER 
# calcium concentration (mean_Cer).
# 

v(0)=-71
cer(0)=81

# conductance in pS
# currents in fA
# Ca concentrations in uM
# time in ms
# capacitance in fF

# Parameters

# Run control parameters:
# set tg=1 for thapsigargin, otherwise 0:
# set pulses=1 for pulses, pulses=0 for no pulses

par tg=0
par changeG=0, changeCa=0
par pulses=1

# Other pulse protocol parameters:
number toff=6000000, tend=6000000
num vhold=-71, tfirst=5000
num tpulse=5000, period=24000, tchange=300000
num vfirst=31, vpulse=-20

#Ica
par gca=1000, vca=25
num vm=-20, sm=12

# Ca fluxes (in uM ms^(_1))
par pleak=0.0002, kpmca=0.1, fcyt=0.01, fer=0.01
# vcyt_er = v_cyt/v_er; only ratio is needed
par vcyt_er=30
par kserca3=0.1, kserca2b=0.01
# Miscellaneous
number alpha=4.5e-06

# Functions
minf = 1/(1+exp((vm-v)/sm))
gcatot = gca + gca*heav(t-tchange)*changeCa
ica = gcatot*minf*(v-vca)
Jin =-alpha*ica
Jpmca = kpmca*ceq
Jserca = (1 - tg)*(kserca2b + kserca3*ceq)
Jleak = pleak*(cer-ceq)

# Apply the pulse protocol
ts = t-tpulse
thyp = 16000-4000*heav(t-tchange)*changeG
von=vhold+vfirst*(heav(t-tfirst))
voff=vfirst*(heav(t-toff))
vprotocol=von+vpulse*(heav(mod(ts,period))-heav(mod(ts,period)-thyp))*(heav(t-tpulse)-heav(t-tend))-voff

# Equilibrium c  
ceq = (Jin-kserca2b+pleak*cer)/(kpmca+kserca3+pleak)
ceqM = (Jin-kserca2b+pleak*mean_cer)/(kpmca+kserca3+pleak)

# Plateau and nadir
par Jdep=0.0461
par Jrep=0.0131
ceqM_pl=(Jdep-kserca2b+pleak*mean_cer)/(kpmca+kserca3+pleak)
ceqM_n=(Jrep-kserca2b+pleak*mean_cer)/(kpmca+kserca3+pleak)
ceqM_amp=(Jdep-Jrep)/(kpmca+kserca3+pleak)

# ER Ca expressions
omega = (kserca3+pleak)/(kserca3+kpmca+pleak)
par mean_Jin=0.0241
par cer_init=81
b = -fer*vcyt_er*(1-omega)*pleak
ab = -(omega*mean_Jin+(1-omega)*kserca2b)/((1-omega)*pleak)
const = ab+cer_init
mean_cer = const*exp(b*t) - ab

# Equations
v' = 10*(vprotocol - v)*pulses+ (vhold-v)*(1-pulses)
cer' = fer*vcyt_er*(omega*Jin+(1-omega)*kserca2b-(1-omega)*pleak*cer)

# itot in pA
aux tsec=t/1000
aux tmin=t/60000
aux Ceq=ceq
aux Jin=Jin
aux mean_cer=mean_cer
aux CeqM=ceqM
aux CeqM_pl=ceqM_pl
aux CeqM_n=ceqM_n
aux CeqM_amp=ceqM_amp

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=20.0, total=300000, maxstor=200000
@ bounds=10000000, xp=tmin, yp=ceqM, bell=off
@ xlo=0, xhi=5, ylo=0.0, yhi=0.5
done