#
% Beta-cell model with I_crac. 
% for bi-phasic, start with cer = 0.5 and gkatp=5000; then step gkatp to 150.
% for muscarinic bursting, set ip3=0.6 (be patient - the pattern develops
% slowly)
% for chaotic muscarinic bursting, set ip3=0.8
% Ref: Bertram et al., Biophys. J. 68:2323-2332 1995
%
% jm = 1 - j
%
init v=-61,n=0.0005,jm=0.12,c=0.11,cer=9

params gi=810.0,gs=510.0,gk=3900.0,gcrac=75.0,gkca=1200.0,gkatp=150.0
params gleak=0.0,vca=100.0,vk=-70.0
params ip3=0.0,gamma=.000003606,cm=6157.521,lambda=1.85,f=0.01,kca=0.07
params lambdaer=250.0,perl=0.02,sigmaer=5.0,verp=0.24,kerp=0.1,nerp=2.0
params caerbar=3,vcrac=0.0,kdkca=0.55,hil=5
params dact=0.1,dip3=0.2,dinh=0.4

% current functions
sinf =  1.0/(1.0+exp((-16.0-v)/10.0))
minf =  1.0/(1.0+exp((-20.0-v)/7.5))
ninf =  1.0/(1.0+exp((-15.0-v)/6.0))
# jminf =  1.0 - 1.0/(1.0+exp((v+53.0)/2.0))
jminf =  1.0/(1.0+exp(-(v+53.0)/2.0))
cinf =  1.0/(1.0+exp(cer-caerbar))
taun =  9.09/(1.0+exp((v+15.0)/6.0))
tauj =  50000.0/(exp((v+53.0)/4.0)+exp((-53.0-v)/4.0))+1500.0

% currents
iin = gi*minf*(v-vca)
is = gs*sinf*(1.0-jm)*(v-vca)
icrac = gcrac*cinf*(v-vcrac)
ik = gk*n*(v-vk)
ikatp = gkatp*(v-vk)
ica = iin+is
ileak = gleak*(v-vcrac)
ikca = gkca/(1 + (kdkca/c)^hil)*(v-vk)

% ER functions

ainf = 1/(1 + dact/c)
binf = ip3/(ip3 + dip3)
hinfer = 1/(1 + c/dinh)
# Note error in paper: exponent 3 is missing
o = (ainf*binf*hinfer)^3

% Ca fluxes
jerp = verp/(1 + (kerp/c)^nerp)
jmemtot = -f*(gamma*ica + kca*c)
jerleak = perl*(cer-c)
jerip3 = o*(cer-c)
jertot = jerleak + jerip3 - jerp


# ODEs variables

v' =   -(ica + icrac + ik + ikatp + ileak + ikca)/cm
n' =   lambda*(ninf-n)/taun
jm' =   (jminf-jm)/tauj
c' =   (jertot/lambdaer) + jmemtot
cer' = -jertot/(lambdaer*sigmaer)

@ meth=gear, dt=10.0, total=20000, maxstor=20000
@ xp=t, yp=v
@ xlo=0, xhi=20000, ylo=-70, yhi=-10

done