# glycophan.ode
#
# This file contains the program for a beta-cell model coupled to glycolysis.
# This was published by Bertram, Satin, Zhang, Smolen, and Sherman in
# Biophysical Journal, 87:3074-3087, November, 2004. 

# State variables:
#    V    -- membrane potential
#    n    -- activation of delayed rectifier
#    Ca   -- free cytosolic calcium concentration
#    ADP  -- cytosolic ADP concentration
#    Caer -- concentration of free calcium in the endoplasmic reticulum
#    G6P  -- glucose 6-phosphate concentration
#    FBP  -- fructose 1,6-bisphosphate concentration

# Initial conditions:
V(0)=-61
n(0)=0.0001
Ca(0)=0.085
ADP(0)=840
Caer(0)=197
G6P(0)=219
FBP(0)=0.015

# Parameter sets for various behaviors:
#
set compound {R_GK=0.2, gKATPbar=25000, gKCabar=600, k_gamma=10}
set slow {R_GK=0.2, gKATPbar=27000, gKCabar=100, k_gamma=10}
set fast {R_GK=0.4, gKATPbar=25000, gKCabar=600, k_gamma=10}
set subthresh {R_GK=0.2, gKATPbar=30000, gKCabar=100, k_gamma=10}
set accordion {R_GK=0.2, gKATPbar=23000, gKCabar=600, k_gamma=10}

# ----------------------------------------------------------------
# Channel properties

num Cm=5300

num VK=-75, gK=2700, taun=20
ninf = 1/(1+exp(-(16+V)/5))

num gCa=1000, VCa=25
minf = 1/(1+exp(-(20+V)/12))

par gKCabar=600
num Kd=0.5

# Ionic currents:

# IK
Ik = gK*n*(V-VK)

# ICa
ICa = gCa*minf*(V-VCa)

# IKCa
gKCa = gKCabar/(1+(Kd/Ca)^2)
IKCa = gKCa*(V-VK)

# IKATP (see below, after calculation of nucleotide concentrations)
par gKATPbar=23000

# Calcium Handling
num alpha=4.50e-6, kPMCA=0.2, fcyt=0.01, fer=0.01
# sigmav=cyt volume/ER volume = V_cyt/V_er
num pleak=0.0002, sigmav=31
num kSERCA=0.4

Jmem = -(alpha*Ica + kPMCA*Ca)
JSERCA = kSERCA*Ca
Jleak = pleak*(Caer - Ca)
Jer = Jleak - JSERCA

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

# Parameters
#  R_GK--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
#  famp corresponds to f13 in paper
#  fmt  corresponds to f41
#  ffbp corresponds to f23
#  fbt  corresponds to f42
#  fatp corresponds to f43
#  R_GPDH--glyceraldehyde phosphate dehydrogenase rate

# Glycolytic parameters
num K1=30, K2=1, K3=50000, K4=1000
num famp=0.02, fmt=20, ffbp=0.2, fbt=20, fatp=20

# Glycolytic expressions
F6P = 0.3*G6P

# nucleotide concentrations used for R_PFK
num Atot=3000
rad = sqrt((ADP-Atot)^2-4*ADP^2)
ATP = 0.5*(Atot-ADP+rad)
AMP = ADP^2/ATP

# Iterative calculation of R_PFK (cf. Smolen95, Eq. 12)
# 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
# lambda, Vmax as in Smolen95, Eq. 3
num lambda=0.06, Vmax=2
R_PFK=Vmax*(lambda*topa16 + topb)/bottom16

# GPDH flux:
R_GPDH = 0.2*sqrt(FBP)


# KATP channel
num Kdd=17, Ktd=26, Ktt=1

% KATP channel open probability  (cf. Magnus and Keizer, 1998)
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)
oinf = topo/bottomo
gKATP = gKATPbar*oinf
IKATP = gKATP*(V-VK)

# glycolytic input to mitochondrial ADP equation:
par k_gamma=10
num v_gamma=2.2
gamma = v_gamma*(R_GPDH/(k_gamma+R_GPDH))

par r=1
num taua=300000, r1=0.35

# conversion parameter for glycolytic subsystem
# kappa erroneously called lambda in paper; renamed kappa for
# consistency with Smolen, JTB, 1995 and Pedersen et al, 2005.
num kappa=0.005
par R_GK=0.2


# Differential equations

V' = -(IK + ICa + IKCa + IKATP)/Cm
n' = (ninf-n)/taun
Ca' = fcyt*(Jmem + Jer)
ADP' = (ATP-ADP*exp((r + gamma)*(1-Ca/r1)))/taua
Caer' = -fer*sigmav*Jer
G6P' = kappa*(R_GK - R_PFK)
FBP' = kappa*(R_PFK - 0.5*R_GPDH)

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=20.0, total=600000,
@ maxstor=100000,bounds=10000000, xp=tmin, yp=V
@ xlo=0, xhi=10, ylo=-70, yhi=-10

# Auxiliary quantities for plotting:
aux tsec=t/1000
aux tmin=t/60000
aux Atp=ATP/1000
aux ratio=ATP/ADP
aux mito_inpt=r+gamma
aux gkca=gKCa
aux gkatp=gKATP
aux f6p=F6P

done