# BJ_09gly.ode
#
# This XPPAUT file contains the program for glycolytic oscillations in beta-cells,
# published by Fendler, Zhang, Satin, and Bertram, Biophys. J., 97:722-729, 2009.
# It was used to make figures 3 and 4.

# Variables:
# g6p -- glucose 6-phosphate
# fbp -- fructose 1,6-bisphosphate

g6p(0)=200
fbp(0)=0

# 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
#  speed--modifies period of glycolytic oscillations.

#moved to expression:par atp=2000
atot=3000
num speed=5
par Rgk=0.2
num phigk=0.3
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

# Auto parameters
par autoatp=1,atpknot=1300

# Expressions
#Sinusoidal ATP function:
atp = atp0*sin(2*pi*t/period)*(1-osc)+2000
#Periodic Pulsatile ATP function
#atp = 2000-atp0*(heav(mod(t,period))-heav(mod(t,period)-dur))

#new params added
param osc=0
param atp0=50
param periodM=5
param durM=0.5
period=periodM*(60)
dur=durM*(60)

f6p = 0.3*g6p
Rgpdh = 0.2*sqrt(fbp)

# nucleotide concentrations, based on eqs. 14 and 15 from Smolen, ref. 1265.
rad = sqrt(1-4*(1/atp)*(atp-atot))
adp=0.5*(rad-1)*atp
amp = adp^2/atp
ratio = atp/adp

# 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
% rate=(pfkbas*cat*topa16 + cat*topb)/bottom16  
% pfk=rate*(atp/(atp+2))    atp saturates
pfk=(pfkbas*cat*topa16 + cat*topb)/bottom16

# ODEs
g6p' = speed*(Rgk - pfk)
fbp' = speed*(pfk - 0.5*Rgpdh)

aux tmin=t/60
aux ttrans=(t-6000)/60
aux F6P=f6p
aux RGPDH=Rgpdh
aux atp=atp

@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=0.5, total=36000, maxstor=70000
@ bounds=10000000, xp=tmin, yp=fbp, trans=6000
@ xlo=0, xhi=100, ylo=0, yhi=40

done