# Vclamp.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. The calcium ode is solved (c), as
# well as an analytical expression for the quasi-equilibrium calcium (Ceq).

# Variables:
#
#  v -- clamped voltage
#  c -- cytosolic calcium
#  cer -- ER calcium
#  cer2 -- ER calcium using the quasi-equilibrium cytosolic calcium

v(0)=-71
c(0)=0.06
cer(0)=81
cer2(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
par vhold=-71, tfirst=5000
par tpulse=5000, period=24000, tchange=300000
par vfirst=31, vpulse=-20

#Ica
par gca=1000, vca=25, 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*c
Jserca = (1 - tg)*(kserca2b + kserca3*c)
Jrelease = pleak*(cer-c)

# 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)
omega = (kserca3+pleak)/(kserca3+kpmca+pleak)

# Equations
v' = 10.0*(vprotocol - v)*pulses+ (vhold-v)*(1-pulses)
c' = fcyt*(Jin-Jpmca-Jserca+Jrelease)
cer' = fer*vcyt_er*(Jserca-Jrelease)
cer2' = 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

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