# Morris-Lecar beta-cell (based on Utah model with kca instead of s)
# parameters chosen to illustrate Chay-Keizer + ER + constant gkatp
# crac added. Developed by Artie.

init v=-65.0, n=0.00016, c=0.12, d=0.5, cer=90.0

# equations
dv/dt=(i+gl*(vl-v)+gk*n*(vk-v)- Ica(v) - Ikca(v,c) + gkatp*(vk-v) -Icrac(v))/cm
dn/dt=(ninf(v)-n)/tau(v)
dc/dt=autoc*(f*(-alpha*ica(v) - kc*c) + jertot/lambdaer) + (1-autoc)*(cknot-c)
dd/dt = ka*((tot-d)-d*exp(r*(1.0-(c/r1))))
dcer/dt = autocer*(-jertot/(lambdaer*sigmaer)) + (1-autocer)*(cerknot-cer)

par autoc=1,cknot=0.1
par autocer=1,cerknot=100
par cstar=1
# par lambdaer=250, perl=0.0005, sigmaer=0.04, kserca=0.4, epser=1
par lambdaer=100, sigmaer=0.12
par perl=0.0008, kserca=0.6
# par lambdaer=200, sigmaer=0.032, kserca=0.6
# par perl0=0.0008, perl1=0.08, tstep=40000
# perl= perl0 + heav(t-tstep)*(perl1 - perl0)
# to remove ER set epser to 0:
par  epser=1

# Icrac
par cerbar=40, gcracbar=20
number vcrac=-30, sloper=3

gcrac = gcracbar/(1+exp((cer-cerbar)/sloper))

icrac(v) = gcrac*(v-vcrac)

jerp = kserca*c
jerleak = perl*(cer - c)
jertot = epser*(jerleak - jerp)

minf(v)=.5*(1+tanh((v-v1)/v2))
ninf(v)=.5*(1+tanh((v-v3)/v4))
tau(v)=1/(phi*cosh((v-v3)/(2*v4)))
#
param gkcastar=100.0, auto=0, gkatpstar=150, freezeatp=1
number vk=-75, vl=-60, vca=25, gk=2700, gl=0, gca=1000, cm=5300
number i=0, kd=5 
number v1=-20.0, v2=24, v3=-16, v4=11.2, phi=0.035
param gkca=2000.0
# Ikca
# ikca(v,c) = auto*gkcastar*(v-vk) + (1-auto)*gkca/(1+(kd/c))*(v-vk)
# linearized:
ikca(v,c) = gkca*((1-auto)*c + auto*cstar)/kd*(v-vk)
# Calcium Handling
par alpha=4.50e-6, kc=0.15
par f=0.01
# Ikatp
par gkatpbar=4000.0
# ADP parameters (from Smolen-Keizer)
# p kone=0.45,ktwo=0.012,ka=5.0e-5,r=0.9,r1=0.35,tot=1.0
number kone=0.45,ktwo=0.012,ka=5.0e-5, r1=0.35,tot=1.0
p r=1

ica(v) = gca*minf(v)*(v-vca)
gkatp = freezeatp*gkatpstar + (1-freezeatp)*gkatpbar*(1.0+d/kone)/(1.0+d/kone+(tot-d)/ktwo)

# arrays for animator; v, c mapped into (0,1),(0,1)

par vmin=-70, vmax=-15, cmin=0.0, cmax=0.3, ttotal=60000
#theta = vmin + t*(vmax - vmin)

zee(v) = (i+gl*(vl-v)+gk*ninf(v)*(vk-v)- Ica(v) + gkatp*(vk-v) -Icrac(v))/(gkca*(v - vk)/kd)

cnull(v) = (-f*alpha*ica(v) + epser*Perl*cer/lambdaer)/(f*kc + epser*(Perl + kserca)/lambdaer)

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

vee[1..51] = (vvec[j] - vmin)/(vmax - vmin)
null[1..51] = (cnull(vvec[j]) - cmin)/(cmax - cmin)
 

vscale = (v - vmin)/(vmax - vmin)
cscale = (c - cmin)/(cmax - cmin)
tscale = t/ttotal

# track some currents
aux Condkca=gkca/(1+(kd/c))
aux Ica=gca*minf(V)*(V-Vca)
aux Ik=gk*n*(V-Vk)
aux Curkca=Ikca(v,c)
aux gkatpact = gkatp
aux gcracact = gcrac
#aux zee = (i+gl*(vl-theta)+gk*ninf(theta)*(vk-theta)- Ica(theta) + gkatp*(vk-theta) -Icrac(theta))/(gkca*(theta - vk)/kd)
#aux null = (-f*alpha*ica(theta) + Perl*cer/lambdaer)/(f*kc + (Perl + kserca)/lambdaer)
#aux theta=theta

@ meth=cvode, atol=1.0e-6, tol=1.0e-6, dt=20.0, total=60000, maxstor=10000
@ xp=t, yp=v, bound=100000000
@ xlo=0, xhi=60000, ylo=-80, yhi=-10

done