COMMENT

Author: Ching-Lung Hsu (adapted from Hines and Carnevale, 2000)

Calcium ion handling mechanisms with endogenous buffer, radial & longitudinal diffusion and pump
A mobile calcium indicator OGB5n is included

ENDCOMMENT

NEURON {
  SUFFIX cdp5
  USEION ca READ cao, cai, ica WRITE cai
  RANGE ica_pmp, BTC_ca, TotalBuffer, DCa
  GLOBAL vrat, TotalPump
  GLOBAL TotalIndicator
    : vrat must be GLOBAL--see INITIAL block
    : however TotalBuffer and TotalPump may be RANGE
}

DEFINE Nannuli 4

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

PARAMETER {
  DCa   = 0.6 (um2/ms)
  
  k1buf = 570 (/mM-ms) : Yamada et al. 1989 : forward rate        : Kd = 1/KD = k1/k2
  k2buf = 5.7 (/ms) : backward rate
  TotalBuffer = 0.39  (mM)
  
  Dogb = 0.3 (um2/ms) : smaller than DCa (Schmidt et al., 2003; Cornelisse et al., 2007)
                      : An unexpected syntax error -
                      : Don't use "DIndicator" as the parameter name for indicator's diffusion constant.
                      : Because after cdp.mod is compiled into cdp.c, the C code will automatically use DIndicator
                      : as an array whenever the state variable Indicator is involved in any REACTION
                      : (including the radial diffusion and the calcium-indicator binding reaction);
                      : otherwise, the following error message will show up in the shell:
                      : "cdp.c: In fuction _ode_specl:
                      :  error: subscripted value is neither array nor pointer"
  k1ind = 570 (/mM-ms)      : Schmidt et al. (2003)
  k2ind = 20 (/ms)        : Schmidt et al. (2003)
  TotalIndicator = 1 (mM) : take 50% of the pipette concentration as the effective concentration at equilibrium 

  k1    = 1       (/mM-ms)
  k2    = 0.005   (/ms)
  k3    = 1       (/ms)
  k4    = 0.005   (/mM-ms)
  : to eliminate pump, set TotalPump to 0 in hoc
  TotalPump = 1e-11 (mol/cm2)
}

ASSIGNED {
  diam      (um)
  ica       (mA/cm2)
  ica_pmp   (mA/cm2)
  ica_pmp_last   (mA/cm2)
  parea     (um)     : pump area per unit length
  cai       (mM)
  cao       (mM)
  vrat[Nannuli]  (1) : dimensionless
                     : numeric value of vrat[i] equals the volume 
                     : of annulus i of a 1um diameter cylinder
                     : multiply by diam^2 to get volume per um length
  Kd        (/mM)
  B0        (mM)
  Kd_ind    (/mM)    : Ching-Lung
  B0_ind    (mM)     :
  BTC_ca    (mM)
}

CONSTANT { volo = 1e10 (um2) }

STATE {
  : ca[0] is equivalent to cai
  : ca[] are very small, so specify absolute tolerance
  : let it be ~1.5 - 2 orders of magnitude smaller than baseline level
  ca[Nannuli]       (mM) <1e-7>
  CaBuffer[Nannuli] (mM) <1e-5>
  Buffer[Nannuli]   (mM) <1e-5>

  CaIndicator[Nannuli] (mM) <1e-5> :
  Indicator[Nannuli]   (mM) <1e-5> :
  
  pump              (mol/cm2) <1e-15>
  pumpca            (mol/cm2) <1e-15>
}

BREAKPOINT {
  SOLVE state METHOD sparse
  ica_pmp_last = ica_pmp
  ica = ica_pmp
}

LOCAL factors_done

INITIAL {
   if (factors_done == 0) {  : flag becomes 1 in the first segment
      factors_done = 1       :   all subsequent segments will have
      factors()              :   vrat = 0 unless vrat is GLOBAL
   }

  Kd = k1buf/k2buf
  B0 = TotalBuffer/(1 + Kd*cai)
  
  Kd_ind = k1ind/k2ind                      
  B0_ind = TotalIndicator/(1 + Kd_ind*cai)

:outermost shell does not have immobile buffer

    ca[0] = cai
    Buffer[0] = 0
    CaBuffer[0] = 0
    Indicator[0] = B0_ind                       
    CaIndicator[0] = TotalIndicator - B0_ind

  FROM i=1 TO Nannuli-1 {
    ca[i] = cai
    Buffer[i] = B0
    CaBuffer[i] = TotalBuffer - B0
    Indicator[i] = B0_ind                       
    CaIndicator[i] = TotalIndicator - B0_ind
  }
BTC_ca=CaIndicator[0]
  parea = PI*diam

: Manually computed initalization of pump
: assumes that cai has been equal to cai0_ca_ion for a long time
  pump = TotalPump/(1 + (cai*k1/k2))   : a better initialization
  pumpca = TotalPump - pump            
: If possible, instead of using formulas to calculate pump and pumpca,
: let NEURON figure them out--just uncomment the following four statements
  ica=0
  ica_pmp = 0
  ica_pmp_last = 0
  SOLVE state STEADYSTATE sparse
: This requires that pump and pumpca be constrained by the CONSERVE
: statement in the STATE block.

}
COMMENT
As suggested by Ted Carnevale,  
the initialization style below may work best

  pump = TotalPump/(1 + (cai*k1/k2))
  pumpca = TotalPump - pump
  SOLVE state STEADYSTATE sparse
ENDCOMMENT


LOCAL frat[Nannuli]  : scales the rate constants for model geometry
PROCEDURE factors() {
  LOCAL r, dr2
  r = 1/2                : starts at edge (half diam)
  dr2 = r/(Nannuli-1)/2  : full thickness of outermost annulus,
                         : half thickness of all other annuli
  vrat[0] = 0
  frat[0] = 2*r
  FROM i=0 TO Nannuli-2 {
    vrat[i] = vrat[i] + PI*(r-dr2/2)*2*dr2  : interior half
    r = r - dr2
    frat[i+1] = 2*PI*r/(2*dr2)  : outer radius of annulus
                                : div by distance between centers
    r = r - dr2
    vrat[i+1] = PI*(r+dr2/2)*2*dr2  : outer half of annulus
  }
}

LOCAL dsq, dsqvol  : can't define local variable in KINETIC block
                   :   or use in COMPARTMENT statement

KINETIC state {
  : COMPARTMENT i, diam*diam*vrat[i] {ca CaBuffer Buffer}
  COMPARTMENT i, diam*diam*vrat[i] {ca CaBuffer Buffer CaIndicator Indicator}
  COMPARTMENT (1e10)*parea {pump pumpca}
  COMPARTMENT volo {cao}
  
  LONGITUDINAL_DIFFUSION i, DCa*diam*diam*vrat[i] {ca}        : Longitudinal diffusion of free Ca
                                                              : Longitudinal diffusion of the Indicator,
                                                              : assuming CaIndicator and Indicator have the same mobility
  LONGITUDINAL_DIFFUSION i, Dogb*diam*diam*vrat[i] {CaIndicator Indicator}

  :pump                                                       : a calcium sink
  ~ ca[0] + pump <-> pumpca  (k1*parea*(1e10), k2*parea*(1e10)) : Nannuli = 0 is the outermost annulus
  ~ pumpca <-> pump + cao    (k3*parea*(1e10), k4*parea*(1e10))
  CONSERVE pump + pumpca = TotalPump * parea * (1e10)
  ica_pmp = 2*FARADAY*(f_flux - b_flux)/parea

  : all currents except pump
  : ica is Ca efflux
  ~ ca[0] << (-(ica - ica_pmp_last)*PI*diam/(2*FARADAY))       : calcium source
  
  FROM i=0 TO Nannuli-2 {                                      : Radial diffusion of free Ca
    ~ ca[i] <-> ca[i+1]  (DCa*frat[i+1], DCa*frat[i+1])        
    ~ Indicator[i] <-> Indicator[i+1]  (Dogb*frat[i+1], Dogb*frat[i+1])      : Radial diffusion of the Indicator,
    ~ CaIndicator[i] <-> CaIndicator[i+1]  (Dogb*frat[i+1], Dogb*frat[i+1])  : assuming CaIndicator and Indicator have the same mobility
  }
  
  dsq = diam*diam                                              : calcium buffering
  FROM i=0 TO Nannuli-1 {
    dsqvol = dsq*vrat[i]
    ~ ca[i] + Buffer[i] <-> CaBuffer[i]  (k1buf*dsqvol, k2buf*dsqvol)
    ~ ca[i] + Indicator[i] <-> CaIndicator[i]  (k1ind*dsqvol, k2ind*dsqvol)
  }
  
  cai = ca[0]
  BTC_ca=CaIndicator[0]
    FROM i=1 TO Nannuli-1 {
    dsqvol = dsq*vrat[i]
BTC_ca=BTC_ca + CaIndicator[i]
  }
}