# Six Subunit Model of CaMKII
#
# Graupner, M. and Brunel, N., STDP in a bistable synapse model based on CaMKII and associated signaling pathways, PLoS Comput Biol, 3(11), e221, 2299-2323 (2007). 
#
#   Please note that this file allows to compute the steady-states of the CaMKII 
#   phosphorylation level with respect to calicum. The parameter used here allow 
#   to reproduce the data shown in Fig.3C by the blue line (p.2303) in the above 
#	menioned paper. 
#   The steady-state diagram consits of two separate branches which have to be computed
#   separately. This is the case since the initial point (specified by 'init') has to 
#   be a fixed point on the respective branch. This file allows to compute the lower branch 
#   including the DOWN state. The computation starts at Ca_0 = 0.01 \mu M.
#   
#   Note however that all the dynamic simulations of the model were not done with xppaut. 
#   The dynamics of the CaMKII-system has been implemented in a C++ code. Please contact
#   the authors for further informations.
#
#  this file is set to run:
#  1. start xppaut and load file
#     $ xppaut 
#  2. lauch auto
#     click -> File -> AUTO
#  3. run auto
#     click -> Run -> Steady State 
#  and you will get the fix-points of the system with the bistability
#
# 
# auto parameters
@ NPR=400, NMAX=40000, DSMAX=0.01, DS=.01, PARMIN=0, PARMAX=2
@ AUTOXMIN=0, AUTOXMAX=2, AUTOYMIN=0, AUTOYMAX=210, AUTOVAR=Ta
#
# note that total receptor pop is conserved
# so p0+p1+...+p10 is constant
# this leads to a zero eigenvalue, so we set the total
# receptor population to be p0i=20 and eliminate p0
# this allows AUTO to do its thing without 
# choking
#
#
# inital conditions to start at Ca=0.01
# required to compuate the fix-points including the DOWN state
init B1=0,B2=0,B3=0
init PP1=0.001108,I1P=0.0359
# 
# parameters
param Ca=0.01
param b0i=33.3
param K5=0.1, CaM=0.1
param L1=0.1, L2=0.025, L3=0.32, L4=0.40
param k6=6, k7=6
param PP10=0.2
param k12=6000
param KM=0.4
param k11=500, km11=0.1
param I10=1
param Kdcan=0.053, ncan=3, kcan0=0.1, kcan=18
param Kdpka=0.11, npka=8, kpka0=0.00359, kpka=100


# occupied receptors
rr=sum(0,12)of(shift(B1,i'))
# p0 is whats left from total
B0=b0i-rr

# total activated and inactivated subunit concentrations
tact= B1 + 2*(B2 + B3 + B4) + 3*(B5 + B6 + B7 + B8) + 4*(B9 + B10 + B11) + 5*B12 + 6*B13

# kinetic equations
phossum=B1 + 2*(B2 + B3 + B4) + 3*(B5 + B6 + B7 + B8) + 4*(B9 + B10 + B11) + 5*B12 + 6*B13
#PP1=Ca^3/(KL^3 + Ca^3)
#PP1=base + kpp1*Ca^3/(KL^3 +  Ca^3)*KH^4/(KH^4 + Ca^4)
k10=k12*PP1/(KM + phossum)
#
C=CaM/(1 + L4/Ca + L3*L4/(Ca^2) + L2*L3*L4/(Ca^3) + L1*L2*L3*L4/(Ca^4))
gamma=C/(K5+C) 
vPKA=kpka0 + kpka/(1 + (Kdpka/C)^npka)
vCaN=kcan0 + kcan/(1 + (Kdcan/C)^ncan)

# at last the equations

B1' = 6*k6*gamma^2*B0 - 4*k6*gamma^2*B1 - k7*gamma*B1 - k10*B1 + 2*k10*(B2 + B3 + B4)
#
B2' = k7*gamma*B1 + k6*gamma^2*B1 - 3*k6*gamma^2*B2 - k7*gamma*B2 - 2*k10*B2 + k10*(2*B5 + B6 + B7)
B3' = 2*k6*gamma^2*B1 - 2*k7*gamma*B3 - 2*k6*gamma^2*B3 - 2*k10*B3 + k10*(B5 + B6 + B7 + 3*B8) 
B4' = k6*gamma^2*B1 - 2*k7*gamma*B4 - 2*k6*gamma^2*B4 - 2*k10*B4 + k10*(B6 + B7)
#
B5' = k7*gamma*B2 + k7*gamma*B3 + k6*gamma^2*B2 - k7*gamma*B5 - 2*k6*gamma^2*B5 - 3*k10*B5 + k10*(2*B9 + B10)
B6' = k6*gamma^2*B2 + k6*gamma^2*B3  + 2*k7*gamma*B4 - k6*gamma^2*B6 - 2*k7*gamma*B6 - 3*k10*B6 + k10*(B9 + B10 + 2*B11)
B7' = k6*gamma^2*B2 + k7*gamma*B3 + 2*k6*gamma^2*B4 - k6*gamma^2*B7 - 2*k7*gamma*B7 - 3*k10*B7 + k10*(B9 + B10 + 2*B11)
B8' = k6*gamma^2*B3 - 3*k7*gamma*B8 - 3*k10*B8 + k10*B10
#
B9' = k7*gamma*B5 + k6*gamma^2*B5 + k7*gamma*B6 + k7*gamma*B7 - k6*gamma^2*B9 - k7*gamma*B9 - 4*k10*B9 + 2*k10*B12
B10'= k6*gamma^2*B5 + k6*gamma^2*B6 + k7*gamma*B7 + 3*k7*gamma*B8 - 2*k7*gamma*B10 - 4*k10*B10 + 2*k10*B12
B11'= k7*gamma*B6 +  k6*gamma^2*B7 - 2*k7*gamma*B11 - 4*k10*B11 + k10*B12
#
B12'= k7*gamma*B9 + k6*gamma^2*B9 + 2*k7*gamma*B10 + 2*k7*gamma*B11 - k7*gamma*B12 - 5*k10*B12 + 6*k10*B13
#
B13'= k7*gamma*B12 - 6*k10*B13
#
PP1'= -k11*I1P*PP1 + km11*(PP10 - PP1)
I1P'= -k11*I1P*PP1 + km11*(PP10 - PP1) + vPKA*I10 - vCaN*I1P


# dummy to get steady-state value of total phosphate - this can be plotted now
# in AUTO!
ta'=-ta+tact
aux act=tact
#@ total=2000,dt=5,meth=cvode
#@ total=100,dt=0.001
@ total=1000,dt=0.001
@ bound=100000
@ maxstor=100000
@ njmp=10

done