# This computer program was used in the publication "Chronic Stimulation
# Induces Adaptive Potassium Channel Activity That Restores Calcium
# Oscillations in Pancreatic Islets in Vitro", Am. J. Physiol.,
# 318:E554-E563, 2020. Authors are Nathan C. Law, Isabella Marinelli, 
# Richard Bertram, Kathryn L. Corbin, Cara Schildmeyer, and
# Craig S. Nunemaker.
#
# Variables:
#    v -- voltage (mV)
#    n -- activation of delayed rectifier
#    c -- free cytosolic calcium concentration (muM)
#    adp -- cytosolic ADP concentration (muM)
#    cer -- concentration of free calcium in the endoplasmic reticulum (muM)
#    cam -- mitocondrial calcium concentration (muM)
#    f6p -- fructose 6-phosphate concentration (muM)
#    fbp -- fructose 1,6-bisphosphate concentration (muM)
#    adpm -- mitocondrial ADP concentration (muM)
#    nadhm -- mitocondrial NADH concentration (muM)
#    psim -- mitocondrial membrane potential (mV)

# Units:
# time: ms
# volume: l
# conductance: pS
# V: mV
# I: fA
# J: muM/ms
# alpha: mumol/fA/ms

# To simulate chronic tolbutamide application increase gkatp from 19,700 
# to 21,500 pS.

# To simulate acute application of tolbutamide, decrease ktt from 1 to 0.91 uM.
 
# To simulate 30 mM KCl, increase vk from -75 to -72 mV.

# Initial Conditions

v(0)=-62
n(0)=0
c(0)=0.05
cer(0)=360
cam(0)=3
adp(0)=950
f6p(0)=17
fbp(0)=110
adpm(0)=1400
nadhm(0)=450
Psim(0)=170


##################################
## Electrical and Calcium Model ##
##################################

# Ca2+ current, ica
num gca = 1000, vca = 25

num nim = -20, sm = 12
minf = 1/(1+exp((nim-v)/sm))

ica=gca*minf*(v-vca)


# Delayed-rectifying K+ current, ik
num gk=2700, 
par vk=-75
num taun=20

num nin=-16, sn=5
ninf=1/(1+exp((nin-v)/sn))

ik=gk*n*(v-vk)

#  Ca2+-activated K+ current, ikca
num kd=0.5, gkca=150

qinf=c^2/(kd^2+c^2)
ikca=-gkca*qinf*(vk-v) 

# K-ATP channel current, ikatp
num kdd=17, ktd=26, atot=3000
par ktt=1, gkatpbar=19600,

rad=sqrt(-4*adp^2+(atot-adp)^2)
atp=(atot+rad-adp)/2
mgadp=0.165*adp
adp3m=0.135*adp
atp4m=0.05*atp
topo=0.08+0.89*mgadp^2/kdd^2+0.16*mgadp/kdd 
bottomo=(1+mgadp/kdd)^2*(1+atp4m/ktt +adp3m/ktd) 
katpo=topo/bottomo 

ikatp=gkatpbar*katpo*(v-vk)


# Ca2+ fluxes across the plasma membrane
num alpha=5.18E-18, vcyt=1.15E-12, kpmca=0.2

Jmem=-(alpha/vcyt*ica + kpmca*c)

# Ca2+ fluxes into ER
num kserca=0.4, pleak=2.0E-4

Jer=kserca*c - pleak*(cer-c) 


# Uniporter [uM/ms]
num p21=0.013, p22=1.6

Juni=(p21*Psim-p22)*c^2

# Na/Ca exchanger [uM/ms]
num p23=0.0015 p24=0.016

JNaCa=p23*(cam-c)*exp(p24*Psim)

# Ca2+ fluxes into the mitochondria

Jm=Juni-Jnaca


# Differential equations for electrical and calcium model
num Cm=5300, fca=0.01, sigmaer=31, sigmam=290

v'=-(ica + ik + ikca + ikatp)/Cm
n'=-(n-ninf)/taun
c'= fca*(Jmem -  Jm/sigmam - Jer)
cer'=fca*sigmaer*Jer
cam'=fca*Jm


##################################
##        Metabolic Model       ##
##################################


# Glucokinase reaction rate
par Ge=11e3
num Vgk=0.0037 Kgk=19e3
Jgk=Vgk/(1+(Kgk/Ge)^2)


# PFK reaction rate, Jpfk
amp=adp^2/atp

#  Parameters
#  k1--Kd for AMP binding
#  k2--Kd for FBP binding
#  k3--Kd for F6P binding
#  k4--Kd for ATP binding
# alpha=1 -- AMP bound
# beta=1 -- FBP bound
# gamma=1 -- F6P bound
# delta=1 -- ATP bound

num k1=30, k2=1, k3=50000, k4=1000
num famp=0.02, fatp=20, ffbp=0.2, fbt=20, fmt=20
num kpfk =0.06,  vpfk=0.01
# (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

Jpfk= vpfk*(topb + kpfk*topa16)/bottom16

# PDH reaction rate, Jgpdh
num kCaPDH=1.5

sinfty= cam/(cam+kCaPDH)

Jgpdh=sinfty*sqrt(fbp)


# Adenine nucleotide translocator, Jant
num amtot=15000
num p19=0.6, p20=2, 
num FRT=.037

atpm=amtot-adpm
Jant= p19/(1+p20*adpm/atpm)*exp(FRT/2* Psim)

# Hydrolysis, Jhyd 
num khyd=1.8640e-06, khydbas=6.4800e-07

Jhyd = (khyd*c+khydbas)*atp

# The differential equations for the glycolytic model

adp'= Jhyd - Jant/sigmam
f6p'=0.3*(Jgk-Jpfk)
fbp'=(Jpfk-Jpdh/2/sigmam)



##################################
##      Mitochondria Model      ##
##################################

# Phosphorylation. In uM/ms.
num p13=10000,p14=190,p15=8.5,p16=4

b2=(p16*p13)/(p13+ATPm)
JF1F0=b2/(1.0+exp((p14-Psim)/p15))
# Respiration (uM/ms)
num p4=5.28 p5=250,p6=165,p7=5

JO=p4*nadhm/(p5+nadhm)/(1+exp((Psim-p6)/p7))


# Pyruvate dehydrogenase (PDH) (uM/ms)
num nadmtot=10000,t2=1.3, vpdh=0.4

nadm=nadmtot-nadhm

Jpdh= vpdh*Jgpdh/(t2+nadhm/NADm)

% Other dehydrogenase (DH) (uM/ms)
num vdh=1.1  kCaDH=0.08	t21=1.3 
Jdh= vdh*cam/(cam+kCaDH)/(t21+nadhm/NADm)

# H+ leakage through mitochondrial inner membrane (uM/ms)
num p17=0.0014,p18=0.02

JHleak=p17*Psim-p18

# Proton pumping due to respiration  (uM/ms)
num p8=7.4,p9=100,p10=165,p11=5

JHres=p8*nadhm/(p9+nadhm)/(1+exp((Psim-p10)/p11))

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


# The vector field (Mitochondria Model)
num Cmito=180

adpm'= Jant-JF1F0
nadhm' = Jpdh+ Jdh -JO
Psim' = (JHres-JHatp-Jant-JHleak-JNaCa-2*Juni)/Cmito


############################################
#         XPP: numerical details	   #
############################################


@ meth=cvode, toler=1.0e-10, atoler=1.0e-10, dt=10,
@ total=900000,  maxstor=20000,bounds=10000000
@ xp=tmin, yp=v,xlo=0, xhi=15, ylo=-70, yhi=-10

aux tmin=t/60000
aux ratio=atp/adp

done