# BMB_04.ode
#
# This XPPAUT file contains the program for pancreatic beta-cells, published
# by Bertram and Sherman, Bull. Math. Biol., 66:1313-1344, 2004.

# Variables:
#    v -- voltage
#    n -- activation variable for a delayed rectifier
#    c -- free cytosolic calcium concentration
#    cer -- concentration of free calcium in the endoplasmic reticulum
#    a -- fraction of activated KATP channels

v(0)=-60
n(0)=0
c(0)=0.1
a(0)=0.46
cer(0)=95

# Persistence of slow bursting in tolbutamide and Tg
# gkatpbar=300, gkca=100, r=0.04, kserca=0

# Conductances in pS; currents in fA; Ca concentrations in uM; time in ms

# conductances
par gkatpbar=500
# Ikca (800 = fast , 700 = medium, 100 = slow)
par gkca=700
# ADP parameters (no glucose, r=0, kserca=0.05
                  low glucose, r=0.04, kserca=0.1
#                 med. glucose, r=0.14, kserca=0.4)
par r=0.14, kserca=0.4
# Ca parameters  (sigmav=cyt volume/ER volume)
par sigmav=5, kc=0.2, ip3=0
# ACh current
par gach=0, vach=0

# Ik
number vn=-16, vk=-75, taun=20, gk=3000, sn=5
# Ica
number vca=25, gca=1200
number vm=-20, sm=12
# Ikca
number kd=0.3
# Miscellaneous
number iapp=0, lambda=1.25
number cm=5300
# Calcium Handling
number f=0.01, fer=0.01
number alpha=4.50e-6
# Calcium Handling: ER
par perl=0.0005
number dact=0.35, dip3=0.5, dinh=0.4
number sa=0.1, taua=300000

# Activation variables
ninf = 1/(1+exp((vn-v)/sn))
minf = 1/(1+exp((vm-v)/sm))
ainf = 1/(1+exp((r-c)/sa))
omega = 1/(1+(kd/c)^5)

# Ionic currents
ica = gca*minf*(v-vca)
ikca = gkca*omega*(v-vk)
ikatp = gkatpbar*a*(v-vk)
ik = gk*n*(v-vk)
iach = gach*(v-vach)

% ER functions
ainf_ip3 = 1/(1 + dact/c)
hinfer = 1/(1 + c/dinh)
jerp = kserca*c
binf = ip3/(ip3 + dip3)
o = ainf_ip3^3*binf^3*hinfer^3

% Ca fluxes
jmemtot = -(alpha*ica + kc*c)
jerleak = perl*(cer - c)
jerip3 = o*(cer - c)
jertot = jerleak + jerip3 - jerp

# Differential Equations

v' = (-ica - ik - ikatp  - ikca - iach + iapp)/cm
n' =  lambda*(ninf - n)/taun
c' = f*(jertot + jmemtot) 
a' = (ainf-a)/taua 
cer' = -fer*sigmav*jertot

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=20.0, total=300000, 
@ maxstor=20000,bounds=10000000, xp=tsec, yp=v
@ xlo=0, xhi=300, ylo=-70, yhi=-10

aux tsec=t/1000
aux tmin=t/60000

done