TITLE Calcium ion accumulation without diffusion and buffering and calcium induced calcium release
 
                 
NEURON {
	SUFFIX cabalcicr
	USEION ca READ cai, ica WRITE cai
	USEION cask READ caski, icask WRITE caski VALENCE 2
	USEION cacicr READ cacicri, icacicr WRITE cacicri VALENCE 2
	
	:USEION adp READ adpi WRITE adpi VALENCE 0
	
	THREADSAFE :but we make it so (this SHOULD not be a problem so long as each neuron or compartment is handled by one thread)
	RANGE  cainit, icapump,icapumpmax,km, TotalBuffer, k1buf, k2buf, SCALE, DCa, basal, cac, cast, tog, dsk, dcicr, constrict, dense,cacicri,caski
	RANGE imetamax, proxy, kadp, castore
	POINTER ju_p, jcicr_p
}

DEFINE NCOMP 3

UNITS {
	(molar) = (1/liter)
	(mM) =  (millimolar)
	(um) =  (micron)
	(mA) =  (milliamp)
	FARADAY = (faraday) (coulomb) 
	PI = (pi) (1)
}

PARAMETER {
        cainit = 0.0001 (mM)
        dt    (ms)
        :celsius = 35  (degC)
        icapumpmax  = 0.00191  (mA/cm2) : 191
        km = 0.000500         (mM)
        DCa = 0.6 (um2/ms) : diffusion constant for calcium, set to 0 to remove radial diffusion
        k1buf = 100 (/mM-ms) : Yamada et al
        k2buf = 0.1 (/ms)
        TotalBuffer = 0.03  (mM)
        basal = 0 (mM) 
        SCALE = 1
        tog = 1
		
		constrict = 0.5
		cas = 2 (mM)		
		
		dsk = 0.1
		dcicr = 0.01
		dense = 10 : nano domains per unit length

        proxy = 1 : upscaling/downscaling of adp creation using ca-atpase activity as proxy.
		imetamax = 0.00191 (mA/cm2) : technically not amps
		kadp = 0.0005 (mM)
         }

ASSIGNED {
  diam  (um)
  ica   (mA/cm2)
  icapump (mA/cm2)
  icask (mA/cm2)
  icacicr (mA/cm2)
  Kd  (/mM)
  B0  (mM)
  cai (mM)
  caski (mM)
  cacicri (mM)
  frat[NCOMP]
  vrat[NCOMP]
  :castore (mM)
  
  imeta (mA/cm2) : its not actually a current, but its easier to treat as one
  :adpi (mM)
  cast (mM)
  :cac (mM)
  
  ju_p
  jcicr_p
}

STATE {
  ca[NCOMP] (mM) <1e-10>
  CaBuffer[NCOMP] (mM) <1e-10>
  Buffer[NCOMP] (mM) <1e-10>
  castore (mM) <1e-10>
  :adp (mM) <1e-8>
}

BREAKPOINT {
	SOLVE states METHOD sparse
}

INITIAL{
  factors()
   
  cai=cainit
  caski = cai
  cacicri = cai
  
  Kd = k1buf
  B0 = TotalBuffer/(1+Kd*cainit)
  
  FROM i=0 TO NCOMP-1 {
	ca[i] = cai
	Buffer[i] = B0
	CaBuffer[i] = TotalBuffer - B0
  }
  castore = cas
  :adpi = 0
}

:LOCAL frat[Nanp2]
PROCEDURE factors() {
  frat[0] = 0: dsk*PI/4.0 : 
  vrat[0] = constrict*diam*diam*PI/4.0 : constrict is the fraction of the compartment that is not occupied by orgenelles
  
  frat[1] = diam*dense*dsk*PI/2.0 : note smaller diam = higher concentrations 
  vrat[1] = diam*dense*dsk*dsk*PI/8.0 : semicircle of diam dsk 
  : dense is the number of microdomains per unit length
  
  frat[2] = diam*dense*dcicr*PI/2.0 : surface area multiplied by number
  vrat[2] = dense*diam*dcicr*dcicr*PI/8.0 : semicircle of diam dcicr
  }
 
LOCAL dsq, dsqvol, casq, juinf, jcinf



KINETIC states {
  COMPARTMENT i, vrat[i] {ca CaBuffer Buffer}
  COMPARTMENT vrat[0] {castore}
  COMPARTMENT vrat[0] {adp}
  icapump = icapumpmax*(1/(1 + km/ca[0])) :calcium only exits from outer layer
   
  FROM i = 0 TO NCOMP-2 {
    ~ca[i] <-> ca[i+1] (DCa*frat[i+1], DCa*frat[i+1]) : shell-shell diffusion
  }
  :~ca[0] <-> ca[1] (DCa*frat[1], DCa*frat[1])
  :~ca[0] <-> ca[2] (DCa*frat[2], DCa*frat[2]) : shell-shell diffusion

  
  
  FROM i=0 TO NCOMP-1 { : may want to consider not buffering the nano-domain with CaT
    dsqvol = vrat[i] : scaling is in 
    ~ca[i] + Buffer[i] <-> CaBuffer[i] (k1buf*dsqvol, k2buf*dsqvol)
    :~ca[i] <-> castore (tog*ju_p*dsqvol,tog*jcicr_p*dsqvol) : release, reuptake from ER stores
  }
  dsqvol = vrat[2]
  ~ca[2] + Buffer[2] <-> CaBuffer[2] (k1buf*dsqvol, k2buf*dsqvol)
  : ica is the non-interacting calcium component
  ~ca[0] << (-SCALE*(ica+icapump)*PI*diam*(1e4)/(2*FARADAY)) : this is correct because the current is going in to the dendrite with d=diam
  ~ca[1] << (-SCALE*(icask)*PI*diam*(1e4)/(2*FARADAY))
  ~ca[2] << (-SCALE*(icacicr)*PI*diam*(1e4)/(2*FARADAY)) : this is correct because the current is going in to the dendrite with d=diam
  if (tog > 0) {
	~castore <-> ca[2] (tog*jcicr_p*vrat[2],tog*ju_p*vrat[2]) : cicr nanodomain
  
  	:~castore <-> ca[1] (0.01*tog*jcicr_p*vrat[1],tog*ju_p*vrat[1])

 	:~castore <-> ca[0] (0,tog*ju_p*vrat[0])
  } else {
    castore = castore
  }
  
  :~adp << ((proxy*icapump-imeta)*PI*diam*(1e4)/(2*FARADAY)) : while the adp current is actually 0, it still 'works' here as if it were valence 2
  
  
  
  
  FROM i=0 TO NCOMP-1 { : protect against integration edge cases that can lead to non-physiological states
	if(ca[i] < 0){ca[i] = 1e-8}
	if(Buffer[i] < 0){Buffer[i] = 1e-10}
	if(CaBuffer[i] < 0){CaBuffer[i] = 1e-10}
  } 
  :if (adp < 0) {adp = 1e-12}
  :adpi = adp
  cai = ca[0]
  caski = ca[1]
  cacicri = ca[2]
  :castore = cas : don't actually change stores
}