# From 'Local glutamate-glutamine cycling underlies presynaptic ATP homeostasis'
# implementation of the plain model of Table 3 (Equations 1 and 2)
# by Reinoud Maex, University of Hertfordshire, UK,
# to be run with XPP (Bard Ermentrout).
#
# N.B. This code reproduces the black traces of Figure 3.
# 
# the stimulus (rectangular modulation of the workload w)
rect_pulse(t,from,to,amp) = amp * (heav (t-from) - heav (t-to))
# for left column of Figure 3 use
w = 1 + rect_pulse(t,10,60,1.0)
# for right column use
# w = 1 + rect_pulse (t,10,60,-0.8)
#
# the differential equations (Eqs. 1 and 2)
Glu' = w * k11/k12 * orthoP^2 * Gln - k21/k22 * Glu * ATP
ATP' = k31/k32 * Pyr * Glu * orthoP * ADP / ATP - 0.5 * k21/k22 * Glu * ATP - w * ATP - w * k11/k12 * orthoP^2 * Gln 
#
# the static variables Pi and ADP
orthoP = 0.2 * (9 - ATP)
ADP = orthoP
# 
# auxiliary variables for plotting: 
# workload, rate of ATP production, rate of ATPconsumption for Glu packaging, rate of Glu formation 
aux auxw = w
aux ATPprod = k31/k32 * Pyr * Glu * orthoP * ADP / ATP
aux ATPves = 0.5 * k21/k22 * Glu * ATP
aux Gluprod = w * k11/k12 * orthoP^2 * Gln
#
# the constants (see Table 1)
param Pyr=0.04,Gln=0.4
# 
# the free parameters (see Table 1 or Eqs. 8)
par k11=1,k12=18.3
par k21=1,k22=30.5
par k31=800,k32=30.5
# 
# the initial conditions
init Glu=1,ATP=1.5
#
@ total=120,dt=.1,xhi=10,maxstor=5000000
@ xp=t,yp=ATP,xlo=0,xhi=120,ylo=0,yhi=8
done