# glycophan.ode
#
# This file contains the program for a beta-cell model coupled to glycolysis.
# Glucose modulates [Ca2+]i oscillations in pancreatic islets via ionic and glycolytic mechanisms.
#  This file gives figure 11, regime change


# Variables:
#    v -- voltage
#    n -- activation of delayed rectifier
#    c -- free cytosolic calcium concentration
#    adp -- cytosolic ADP concentration
#    cer -- concentration of free calcium in the endoplasmic reticulum
#    g6p -- glucose 6-phosphate concentration
#    fbp -- fructose 1,6-bisphosphate concentration

# These work but have convergence errors:

# v(0)=-61.5
# n(0)=0.0001
# c(0)=0.1
# adp(0)=790
# cer(0)=120
# g6p(0)=210
# fbp(0)=0.01

# These work but also have convergence errors:

#v(0)=-61.45112250316146
#n(0)=0.0001127171676807812
#c(0)=0.05971173349625042
#adp(0)=790.3985348487028
#cer(0)=119.6392837195939
#g6p(0)=207.7355380115209
#fbp(0)=0.009034527652506757

# These work without convergence problems (used last ICs from 2nd fast burst)

v(0)=-62.6618066014772
n(0)=8.848548358296786e-05
c(0)=0.06095304038160815
adp(0)=797.7195667416995
cer(0)=128.0924185806713
g6p(0)=224.2627494899293
fbp(0)=0.01326094549014909


# Parameter vlaues for various behaviors:
#
# Compound bursting: rgk=0.2, gkatp=25000, gkca=600, kg=10
# Glycolytic bursting: rgk=0.2, gkatp=27000, gkca=100, kg=10
# Simple bursting: rgk=0.4, gkatp=25000, gkca=600, kg=10
# Subthreshold oscillations: rgk=0.2, gkatp=30000, gkca=100, kg=10
# Accordion bursting: rgk=0.2, gkatp=23000, gkca=600, kg=10

# ----------------------------------------------------------------
# Membrane and Ca components

# sigmav=cyt volume/ER volume
num pleak=0.0002, sigmav=31
num kserca=0.4
num lambdaer=1, epser=1
num taun=20

# where
minf = 1/(1+exp(-(20+v)/12))
ninf = 1/(1+exp(-(16+v)/5))

# conversion parameter for glycolytic subsystem
num lambda=0.005

num vk=-75, vca=25
num gk=2700, cm=5300
par gca=1000
num kd=0.5

param gkca=25

# Ikca
ikca = gkca/(1+(kd/c)^2)*(v-vk)

# Calcium Handling
num alpha=4.50e-6, kpmca=0.2, fcyt=0.01, fer=0.01

# Ikatp
par gkatpbar=27000

# ICa
ica = gca*minf*(v-vca)

# Ik
Ik = gk*n*(v-vk)

# Ca fluxes
Jmem = -(alpha*Ica + kpmca*c)

Jserca = kserca*c
Jleak = pleak*(cer - c)
Jer = epser*(Jleak - Jserca)/lambdaer

# -----------------------------------------------------
# Glycolytic and Keizer-Magnus components

# Parameters
#  Rgk--glucokinase rate
#  atot--total adenine nucleotide concentration (micromolar)
#  k1--Kd for AMP binding
#  k2--Kd for FBP binding
#  k3--Kd for F6P binding
#  k4--Kd for ATP binding
#  famp,etc--Kd amplification factors for heterotropic binding
#  Rgpdh--glyceraldehyde phosphate dehydrogenase rate

# Glycolytic parameters
par kg=10

# rgk=if(t<960000)then(0.1)else(0.15)
rgk=0.02 + 0.18*heav(t - 600000) + 0.4*heav(t - 28*60000)
r  = 1.0

num atot=3000
num k1=30, k2=1, k3=50000, k4=1000
num famp=0.02, fatp=20, ffbp=0.2, fbt=20, fmt=20, pfkbas=0.06
number cat=2
num katpase=0.0003

# Glycolytic expressions
f6p = 0.3*g6p
Rgpdh = 0.2*sqrt(fbp)

# Iterative calculation of PFK
# alpha=1 -- AMP bound
# beta=1 -- FBP bound
# gamma=1 -- F6P bound
# delta=1 -- ATP bound

# (alpha,beta,gamma,delta)
# (0,0,0,0)
weight1=1
topa1=0
bottom1=1

# (0,0,0,1)
weight2=atp^2/k4
topa2=topa1
bottom2=bottom1+weight2

# (0,0,1,0)
weight3=f6p^2/k3
topa3=topa2+weight3
bottom3=bottom2+weight3

# (0,0,1,1)
weight4=(f6p*atp)^2/(fatp*k3*k4)
topa4=topa3+weight4
bottom4=bottom3+weight4

# (0,1,0,0)
weight5=fbp/k2
topa5=topa4
bottom5=bottom4+weight5

# (0,1,0,1)
weight6=(fbp*atp^2)/(k2*k4*fbt)
topa6=topa5
bottom6=bottom5+weight6

# (0,1,1,0)
weight7=(fbp*f6p^2)/(k2*k3*ffbp)
topa7=topa6+weight7
bottom7=bottom6+weight7

# (0,1,1,1)
weight8=(fbp*f6p^2*atp^2)/(k2*k3*k4*ffbp*fbt*fatp)
topa8=topa7+weight8
bottom8=bottom7+weight8

# (1,0,0,0)
weight9=amp/k1
topa9=topa8
bottom9=bottom8+weight9

# (1,0,0,1)
weight10=(amp*atp^2)/(k1*k4*fmt)
topa10=topa9
bottom10=bottom9+weight10

# (1,0,1,0)
weight11=(amp*f6p^2)/(k1*k3*famp)
topa11=topa10+weight11
bottom11=bottom10+weight11

# (1,0,1,1)
weight12=(amp*f6p^2*atp^2)/(k1*k3*k4*famp*fmt*fatp)
topa12=topa11+weight12
bottom12=bottom11+weight12

# (1,1,0,0)
weight13=(amp*fbp)/(k1*k2)
topa13=topa12
bottom13=bottom12+weight13

# (1,1,0,1)
weight14=(amp*fbp*atp^2)/(k1*k2*k4*fbt*fmt)
topa14=topa13
bottom14=bottom13+weight14

# (1,1,1,0) -- the most active state of the enzyme
weight15=(amp*fbp*f6p^2)/(k1*k2*k3*ffbp*famp)
topa15=topa14
topb=weight15
bottom15=bottom14+weight15

# (1,1,1,1)
weight16=(amp*fbp*f6p^2*atp^2)/(k1*k2*k3*k4*ffbp*famp*fbt*fmt*fatp)
topa16=topa15+weight16
bottom16=bottom15+weight16

# Phosphofructokinase rate
pfk=(pfkbas*cat*topa16 + cat*topb)/bottom16

# KATP channel
num kdd=17, ktd=26, ktt=1

# r=if(t<960000)then(0.7)else(0.9)

par vg=50
num taua=300000, r1=0.35

# nucleotide concentrations
rad = sqrt((adp-atot)^2-4*adp^2)
atp = 0.5*(atot-adp+rad)
amp = adp^2/atp
ratio = atp/adp

% KATP channel open probability 
mgadp = 0.165*adp
adp3m = 0.135*adp
atp4m = 0.05*atp
topo = 0.08*(1+2*mgadp/kdd) + 0.89*(mgadp/kdd)^2
bottomo = (1+mgadp/kdd)^2 * (1+adp3m/ktd+atp4m/ktt)
katpo = topo/bottomo
ikatp = gkatpbar*katpo*(v-vk)

# glycolytic feedback onto adp
y = vg*(Rgpdh/(kg+Rgpdh))
fback = r+y

# Differential equations

v' = -(Ik + Ica + Ikca + Ikatp)/cm
n' = (ninf-n)/taun
c' = fcyt*(Jmem + Jer)
adp' = (atp-adp*exp(fback*(1-c/r1)))/taua
cer' = -fer*sigmav*Jer
g6p' = lambda*(Rgk - pfk)
fbp' = lambda*(pfk - 0.5*Rgpdh)

@ meth=cvode, dt=20.0, total=2700000
@ toler=1.0e-10, atoler=1.0e-10
@ maxstor=1000000,bounds=10000000, xp=tmin, yp=Ca
@ xlo=0, xhi=35, ylo=0, yhi=300

aux tsec=t/1000
aux tmin=t/60000
aux Ca=c*1000
aux gkatpact=gkatpbar*katpo
aux Fback=fback
aux Y=y
aux RGPDH=Rgpdh
aux Atp=atp/1000
aux ratio=atp/adp
aux gkcaact=gkca/(1+(kd/c)^2)

done