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, k1buf, k2buf, Dogb, TotalIndicator, TotalPump, k1ind, k2ind
GLOBAL vrat
: 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) : diffusion coefficient
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) : indicator radial diffusion coefficient (Laila)
: 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) Kon
k2ind = 20 (/ms) : Schmidt et al. (2003) Koff
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]
}
}