# BJ_09beta.ode
#
# This XPPAUT file contains the program for pancreatic beta-cells, published
# by Fendler, Zhang, Satin, and Bertram, Biophys. J., 97:722-729, 2009.
# It was used to make figures 2, 5, 6, and 7.

# Variables:
#    v -- voltage
#    n -- activation variable for a delayed rectifier
#    c -- free cytosolic calcium concentration
#    cer -- concentration of free calcium in the endoplasmic reticulum
#    I -- insulin concentration
#    g6p -- glucose 6-phosphate concentration
#    fbp -- fructose bis-phosphate concentration
#    adpm -- mitochondrial ADP concentration
#    nadhm -- mitochondrial NADH concentration
#    cam -- mitochondrial free calcium concentration
#    psim -- mitochondrial inner membrane potential
#    adp -- cytosolic ADP concentration

v(0)=-63
n(0)=0
c(0)=0.18
cer(0)=381
I(0)=0.57
g6p(0)=327
fbp(0)=0.65

# Mito variables (MGP) -  in mM
adpm(0)=12.0
nadhm(0)=0.14

# Mito Ca - in uM
cam(0)=0.25

# mV (cyto-mito)
psim(0)=160

# Cytosolic ADP concentration - in uM
adp(0)=1509

# ----------------------------------------------------------------

# Fraction of PFK that is M-type.
par Mfrac=1

# Plasma membrane capacitance
num cm=5300

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

# Ikatp
par gkatpbar=12600

# ICa
num vca=25, gca=1000
num vm=-20,sm=12
minf = 1/(1+exp((vm-v)/sm))
ica = gca*minf*(v-vca)

# Ik
num vk=-75, gk=2700
num vn=-16,sn=5,taun=20
ninf = 1/(1+exp((vn-v)/sn))
Ik = gk*n*(v-vk)

# IKATP
# Dissociation constants (from Magnus and Keizer, 1998) 
num kdd=17, ktd=26, ktt=1
par Dz=0

% KATP channel open probability (from Magnus and Keizer, 1998)
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 = (1-Dz)*gkatpbar*katpo*(v-vk)+Dz*gkatpbar*(v-vk)

# ----------------------------------------------------------------
###### Cytosolic Calcium Handling ######
# sigmav=cyt volume/ER volume 
num alpha=4.50e-6, fcyt=0.01, fer=0.01
par kpmca=0.1
num per=0.0002, kserca=0.4, lambdaer=1, epser=1
num sigmav=31
num cbas=0.05

# Ca fluxes (modified 02/13/08)
#new parameters added (modified 02/13/08)
#the following 3 parameters are measured in microM
#IP3=0.3 constant from phantom bursting
par dact=0.35 ,dinact=0.4 ,dIP3=0.5
Oinf(c,IP3)= (c/(dact+c))^3*(IP3/(dIP3+IP3))^3*(dinact/(dinact+c))^3
Jip3   = Oinf(c,IP3)*(cer - c)
Jmem = -(alpha*Ica + kpmca*(c-cbas))
Jserca = kserca*c
Jleak = per*(cer - c)
Jer = epser*(Jleak - Jserca)/lambdaer + Jip3

#input IP3 Pulses
#The following are the parameters and functions for the IP3 pulse
par periodM=5
par durM=0.5
par IP3Amp=0.1

period=periodM*(60000)
dur=durM*(60000)

#Multiple pulses:
IP3 = IP3Amp*(heav(mod(t,period))-heav(mod(t,period)-dur))
#Single pulse:
#IP3 = IP3Amp*heav(t-period)*heav(period+dur-t)
aux ip3=ip3

# -----------------------------------------------------
# Glycolytic components -- muscle type

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

# Parameters
#  Jgk--glucokinase rate (uM/s)
#  Jgk_ms--glucokinase rate (uM/ms)
#  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
#  k1C--Kd for AMP binding, C-type
#  k2C--Kd for FBP binding, C-type
#  k3C--Kd for F6P binding, C-type
#  k4C--Kd for ATP binding, C-type
#  f41--ATP inhibition by weakening binding of AMP
#  f42--ATP inhibition by weakening binding of FBP
#  f43--ATP inhibition by weakening binding of F6P
#  f23--FBP activation by tightening binding of F6P
#  f13--AMP activation by tightening binding of F6P
#  lambda--ratio of PFK activity in inactivated to non-inactivated form

#  famp,etc--Kd amplification factors for heterotropic binding
#  Jgpdh--glyceraldehyde phosphate dehydrogenase rate (uM/ms)
#  kappa--converts from s to ms.
#  Jgpdh_bas--low basal level of GPDH.

# Glycolytic parameters
num phigk=0.3, KGPDH=0.0005
par Jgk=0.5
par k1=30, k2=1, k3=50000, k4=220
par k1C=30, k2C=10, k3C=50000, k4C=100
num f13=0.02, f43=20, f23=0.2, f42=20, f41=20, lambda=0.06
par VmaxPFK=5
num kappa=0.001
num Jgpdh_bas=0.0005

# Glycolytic expressions
f6p = phigk*g6p
Jgpdh = KGPDH*sqrt(fbp) 
Jgk_ms=kappa*Jgk
% aux JGPDH=Jgpdh 

# 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/(f43*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*f42)
topa6=topa5
bottom6=bottom5+weight6

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

# (0,1,1,1)
weight8=(fbp*f6p^2*atp^2)/(k2*k3*k4*f23*f42*f43)
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*f41)
topa10=topa9
bottom10=bottom9+weight10

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

# (1,0,1,1)
weight12=(amp*f6p^2*atp^2)/(k1*k3*k4*f13*f41*f43)
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*f42*f41)
topa14=topa13
bottom14=bottom13+weight14

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

# (1,1,1,1)
weight16=(amp*fbp*f6p^2*atp^2)/(k1*k2*k3*k4*f23*f13*f42*f41*f43)
topa16=topa15+weight16
bottom16=bottom15+weight16

# Phosphofructokinase rate 
pfkM=(lambda*VmaxPFK*topa16 + VmaxPFK*topb)/bottom16


# -----------------------------------------------------
# Glycolytic components -- brain type (C-type)

# Parameters
# alpha=1 -- AMP bound
# beta=1 -- FBP bound
# gamma=1 -- F6P bound
# delta=1 -- ATP bound

# (alpha,beta,gamma,delta)
# (0,0,0,0)
weight1C=1
topa1C=0
bottom1C=1

# (0,0,0,1)
weight2C=atp^2/k4C
topa2C=topa1C
bottom2C=bottom1C+weight2C

# (0,0,1,0)
weight3C=f6p^2/k3C
topa3C=topa2C+weight3C
bottom3C=bottom2C+weight3C

# (0,0,1,1)
weight4C=(f6p*atp)^2/(f43*k3C*k4C)
topa4C=topa3C+weight4C
bottom4C=bottom3C+weight4C

# (0,1,0,0)
weight5C=fbp/k2C
topa5C=topa4C
bottom5C=bottom4C+weight5C

# (0,1,0,1)
weight6C=(fbp*atp^2)/(k2C*k4C*f42)
topa6C=topa5C
bottom6C=bottom5C+weight6C

# (0,1,1,0)
weight7C=(fbp*f6p^2)/(k2C*k3C*f23)
topa7C=topa6C+weight7C
bottom7C=bottom6C+weight7C

# (0,1,1,1)
weight8C=(fbp*f6p^2*atp^2)/(k2C*k3C*k4C*f23*f42*f43)
topa8C=topa7C+weight8C
bottom8C=bottom7C+weight8C

# (1,0,0,0)
weight9C=amp/k1C
topa9C=topa8C
bottom9C=bottom8C+weight9C

# (1,0,0,1)
weight10C=(amp*atp^2)/(k1C*k4C*f41)
topa10C=topa9C
bottom10C=bottom9C+weight10C

# (1,0,1,0)
weight11C=(amp*f6p^2)/(k1C*k3C*f13)
topa11C=topa10C+weight11C
bottom11C=bottom10C+weight11C

# (1,0,1,1)
weight12C=(amp*f6p^2*atp^2)/(k1C*k3C*k4C*f13*f41*f43)
topa12C=topa11C+weight12C
bottom12C=bottom11C+weight12C

# (1,1,0,0)
weight13C=(amp*fbp)/(k1C*k2C)
topa13C=topa12C
bottom13C=bottom12C+weight13C

# (1,1,0,1)
weight14C=(amp*fbp*atp^2)/(k1C*k2C*k4C*f42*f41)
topa14C=topa13C
bottom14C=bottom13C+weight14C

# (1,1,1,0) -- the most active state of the enzyme
weight15C=(amp*fbp*f6p^2)/(k1C*k2C*k3C*f23*f13)
topa15C=topa14C
topbC=weight15C
bottom15C=bottom14C+weight15C

# (1,1,1,1)
weight16C=(amp*fbp*f6p^2*atp^2)/(k1C*k2C*k3C*k4C*f23*f13*f42*f41*f43)
topa16C=topa15C+weight16C
bottom16C=bottom15C+weight16C

# C-type Phosphofructokinase rate 
pfkC=(lambda*VmaxPFK*topa16C + VmaxPFK*topbC)/bottom16C

aux pfkM=pfkM
aux pfkC=pfkC

pfk=Mfrac*pfkM+(1-Mfrac)*pfkC
pfk_ms=kappa*pfk

# ----------------------------------------------------------------
###### Mitochondrial Ca2+ handling #######
# delta transforms  (mito volume) -> (cyto volume)/ms
# gamma transforms uM to mM.

num fmito=0.01
!delta=3.9/53.2
num gamma=0.001

### Uniporter [uM/ms]
par p21=0.04,p22=1.1
Juni=(p21*psim-p22)*c^2
# aux juni=Juni

### Na/Ca exchanger [uM/ms]
par p23=0.01
par p24=0.016
JNaCa=(p23/c)*cam*exp(p24*Psim)
# aux JNC=JNaCa

# [uM/ms] , mito -> cyto
Jmito=JNaCa-Juni

# ----------------------------------------------------------------
###### Mitochondrial fluxes (JTB, 243:575-586, 2006) ######

# Mitochondrial adenine and pyridine nucleotide conservation. All in [mM].

num Amtot=15, NADmtot=10
NADm=NADmtot-NADHm
ATPm=Amtot-ADPm
RATm=ATPm/ADPm

# Pyruvate dehydrogenase (PDH) (uM/ms)
par p1=400,p2=1,p3=0.01
JPDH=(p1/(p2+NADHm/NADm))*(Cam/(p3+Cam))*(Jgpdh+Jgpdh_bas)
# aux jpdh=JPDH

# ----------------------------------------------------------------
###### Mitochondrial membrane components ######
# Inner membrane capacitance
num Cmito=1.8

#H+ leakage through mitochondrial inner membrane (uM/ms)
par p17=0.002,p18=-0.03
JHleak=p17*Psim+p18

# Respiration (uM/ms)
par p4=0.6
par p5=0.1,p6=177,p7=5
MM1=p4*NADHm/(p5+NADHm)
JO=MM1/(1+exp((Psim-p6)/p7))

aux jo=JO

# Proton pumping due to respiration  (uM/ms)
par p8=7,p9=0.1,p10=177,p11=5
MM2=p8*NADHm/(p9+NADHm)
JHres=MM2/(1+exp((Psim-p10)/p11))

# Phosphorylation. In uM/ms.
par p13=10,p14=190,p15=8.5,p16=35
b2=(p16*p13)/(p13+ATPm)
JF1F0=b2/(1.0+exp((p14-Psim)/p15))
# aux jf1f0=JF1F0

# Proton flux due to ATPase (uM/ms)
JHatp=3*JF1F0

# ADP/ATP translocator. In uM/ms.
par p19=0.35, p20=2
FRT=96480/(310.16*8315)
Jant=p19*(RATm/(RATm+p20))/exp(-0.5*FRT*Psim)
aux JANT=Jant
aux RATm=RATm
aux ATPm=ATPm

# ----------------------------------------------------

###### Cytosolic nucleotide concentrations #####
par khyd=0.00005, JhydSS=0.00005
num amp=500,atot=2500
atp = atot-adp
Jhyd=(khyd*c+JhydSS)*ATP 
mgadp = 0.165*adp
adp3m = 0.135*adp
atp4m = 0.05*atp
RATc = atp/adp
aux Jhyd=Jhyd

# ----------------------------------------------------------------

# Insulin secretion component

num taui=1000
num ki=0.15
num deltaI=2

Iinf = (c^deltaI)/(ki^deltaI + c^deltaI)

# ------------------------------------------------------
# Differential equations

v' = -(Ik + Ica + Ikca + Ikatp)/cm
n' = (ninf-n)/taun
c' = fcyt*(Jmem + Jer + delta*Jmito)
cer' = -fer*sigmav*Jer
g6p' = Jgk_ms - pfk_ms
fbp' = pfk_ms - 0.5*Jgpdh
I' = (Iinf-I)/taui

### Mitochondria 
adpm'= gamma*(JANT-JF1F0)
cam' = -fmito*Jmito
nadhm' = gamma*(Jpdh-JO)

# Mitochondrial membrane potential [mV] (cyto-mito)
Psim'=(JHres-JHatp-JANT-JHleak-JNaCa-2*Juni)/Cmito

# Cytosolic ADP concentration
adp' = -delta*JANT + Jhyd

@ meth=83, toler=1.0e-9, atoler=1.0e-9, dt=30, total=6000000,
@ maxstor=300000,bounds=10000000, xp=tmin, yp=c
@ xlo=0, xhi=50, ylo=0, yhi=0.5, bell=off

# aux itot=(ica + ik + ikatp  + ikca + icrac)/1000
aux tsec=t/1000
aux tmin=t/60000
# aux GKATP=gkatpbar*katpo
aux Atp=atp/1000
% aux ratio=atp/adp
# aux PFK=pfk
# aux amp=amp

### average of last qqq ms
# num qqq=60000
# @ delay=60000
# storI'=I/qqq
# aux smoothI=(storI-delay(storI,qqq))

done