#bursting with three slow variable

#units: V=mV; t=ms; g=pS; I=fA
#Reference:  Bertram and Sherman. 
#Calcium-based model for pancreatic islets
#Figure 13A...slow bursting generated with the interaction of 3 slow variables
#burst period is about 275 seconds
#
#Ica- calcium current
#Ik- delayed rectifier K+ current
#Ik(Ca)- Ca2+ dependent K+ current
#Ik(ATP)- nucleotide-sensitive K+ current
#c - cytosolic free Ca2+ concentration
#cer - ER Ca2+ concentration
#initial conditions
init v=-64.0, c=0.087, n=6.75e-5, cer=80, a=0.483

set basic {freezeatp=1, epser=0, astar=0.46, gkca=300, fcyt=0.01}
set basicmedium {freezeatp=1, epser=0, astar=0.46, gkca=300, fcyt=0.0005}
set erfast  {freezeatp=1, epser=1, gkca=900, astar=0.455, fcyt=0.01}
set ermedium  {freezeatp=1, epser=1, gkca=700, astar=0.46, fcyt=0.01}
set erslow  {freezeatp=1, epser=1, gkca=300, astar=0.46, fcyt=0.01}
set atpfast {freezeatp=0, epser=1, gkca=1000, fcyt=0.01}  
set atpmedium {freezeatp=0, epser=1, gkca=700, fcyt=0.01}  
set atpslow {freezeatp=0, epser=1, gkca=100, fcyt=0.01}  

#parameters
par gca=1200, gkca=100, gk=3000
par vca=25, vk=-75, cm=5300
par taun=16, alpha=4.5e-6
par fcyt=0.01, kpmca=0.2, kd=0.3
par vn=-16, vm=-20, sn=5, sm=12
par kserca=0.4, dact=0.35, dinact=0.4
par fer=0.01, pleak=0.0005, dip3=0.5, vcytver=5
par ip3=0, gkatp=500, sa=0.1, r=0.14 taua=300000
par epser=1, freezeatp=0, astar=0.46

# ionic currents
ica(v)=gca*minf(v)*(v-vca)
ik(v)=gk*n*(v-vk)
ikca(v)=gkca*w*(v-vk)
ikatp(v)=gkatp*a*(v-vk)

#activation functions
minf(v)=1.0/(1.0+exp((vm-v)/sm))
ninf(v)=1.0/(1.0+exp((vn-v)/sn))
ainf(c)=1.0/(1.0+exp((r-c)/sa))

#fraction of K(Ca) channels activated by cytosolic Ca2+
w=c^5/(c^5+kd^5)

#flux of Ca2+ through the membrane
jmem=-(alpha*Ica(v)+kpmca*c)

#Ca2+ influx into the ER via SERCA 
jserca=kserca*c

#efflux out of the ER has two components
# 1. Ca2+ leak is proportional to gradient between Ca2+ and ER
jleak=pleak*(cer-c)

# 2. Ca2+ efflux through the IP3R
jip3=oinf*(cer-c)

#fraction of open channels
oinf=(c/(dact+c))*(ip3/(dip3+ip3))*(dinact/(dinact+c))

#net Ca2+ efflux from the ER
jer=jleak+jip3-jserca

#differential equations
v'=-(ica(v)+ik(v)+ikca(v)+ikatp(v))/cm
n'=(ninf(v)-n)/taun
c'=fcyt*(jmem+epser*jer)
cer'=-epser*fer*(vcytver)*jer
a'=(1 - freezeatp)*(ainf(c)-a)/taua + freezeatp*(astar-a)

aux tsec=t/1000.0
aux W=w

# Arrays for animator

par vmin=-70, vmax=-15, wmin=-0.3, wmax=0.3, ttotal=250000

zee(v) = ( -ica(v) - gk*ninf(v)*(v - vk) - ikatp(v) )/( gkca*(v - vk) )

cnull(v) = ( -alpha*ica(v) + epser*pleak*cer )/( kpmca + epser*pleak + epser*kserca )

wnull(v) = cnull(v)^5/( cnull(v)^5 + kd^5 )

vvec[1..51] = vmin + (vmax - vmin)*([j] - 1)/50
zvec[1..51] = (zee(vvec[j]) - wmin)/(wmax - wmin)

vee[1..51] = (vvec[j] - vmin)/(vmax - vmin)
null[1..51] = (wnull(vvec[j]) - wmin)/(wmax - wmin)

vscale = (v - vmin)/(vmax - vmin)
wscale = (w - wmin)/(wmax - wmin)
tscale = t/ttotal



@ meth=cvode, dtmax=1, dt=50, total=250000, maxstor=1000000
@ bounds=1000, xp=tsec,  yp=v, toler=1.0e-7, atoler=1.0e-7
@ xlo=0, xhi=250, ylo=-80, yhi=-20

done