You have chosen a difficult task, and you have done quite well with it. Still, there are
many problems to be fixed before the job is finished. Part of the difficulty is with NMODL,
which has many complexities of its own. But another source of difficulty lies in NEURON's
use of cylindrical sections, and "density" units for membrane current (e.g. mA/cm2) and
other "distributed mechanisms" (mechanisms that are distributed in space). Many
authors who develop their models with Matlab, or C, or other programming tools, make
their model cell have a fixed diameter and length, and specify their variables and
parameters in terms of "absolute" units (e.g. mA). Also, they often assume a spherical
geometry. As a consequence, their model equations and parameter values may be valid
only for a cell with that particular size and shape. In order to port such a model to
NEURON, it is necessary to generalize the model specification so that it can be used
with a cylinder of any diameter or length. Also, if ion accumulation mechanisms (pumps,
buffers, diffusion) are involved, it may be necessary to deal with the differences between
the surface/volume ratios of spherical and cylindrical geometries.
All of the problems with your model are in the er.mod file, and modlunit is the diagnostic
tool for discovering and fixing most of them.
The peculiar result you got when running a simulation may have been caused by failure
of the code in the mod file to initialize Ver, or maybe it was caused by using cnexp (this
mechanism's nonlinearities require derivimplicit, not cnexp).
But it turns out that a NEURON implementation of this mechanism doesn't need the
actual value of Ver--it only needs the ratio of ER volume to cytoplasmic volume.
Here is the mod file, with all fixes so that it passes modlunit, and generates a stable
simulation.
Code: Select all
UNITS {
(mol) = (1/liter)
(mM) = (milli/liter)
: add these so modlunit doesn't complain
(mA) = (milliampere)
(um) = (micron)
FARADAY = (faraday) (coulomb)
PI = (pi) (1)
}
NEURON {
SUFFIX er
USEION ca READ cai WRITE ica
NONSPECIFIC_CURRENT ix
}
STATE {
cai (mM)
caer (mM)
}
PARAMETER {
fer = 0.0025 : ca buffering coefficient in ER from John Rinzel 1997
kerm = 0.2e-3 (mM) : from Huang 2004
kerrel = 3e-12 (/s): from Huang, to be adjusted
kerfila = 0.75e-12 (mM/s) : from John Rinzel 1997
kerfilb = 0.2e-3 (mM) : from John Rinzel 1997
kerlek = 6.15e-14 (/s): from Huang, to be adjusted
rhover = 0.15 : ratio of volume to cytoplasmic volume
: total cell vol = vol_cytoplasm + vol_er
: where vol_er = rhover * vol_cytoplasm
: thus vol_er = v_cell * rhover / (1 + rhover)
caer0 = 0.25 (mM)
}
ASSIGNED {
diam (um)
ica (mA/cm2)
ix (mA/cm2)
}
INITIAL {
: custom initialization may be needed
caer = caer0 : initial calcium concentration in ER
}
BREAKPOINT {
SOLVE store METHOD derivimplicit
ica = -2*FARADAY*(errel(cai,caer,kerrel,kerm)-erfil(cai,caer,kerfila,kerfilb)+erlek(cai,caer,kerlek)) * diam * (1e-7) / 4
ix = -ica
}
DERIVATIVE store {
caer' = -(0.001)*(errel(cai,caer,kerrel,kerm)-erfil(cai,caer,kerfila,kerfilb)+erlek(cai,caer,kerlek))/(rhover/fer)
}
COMMENT
1. Original LHS & RHS units are inconsistent
2. Original RHS is in units of millimoles / second, which doesn't make sense for this mechanism.
Should be mM / second.
3. If RHS is in mM / second, caer' should be insensitive to diam, and only take rhover into account.
Specifically, given that errel is in mM/second,
and that er volume is 0.15 * cytoplasm volume, then if er accumulation of Ca makes
cytoplasmic ca fall at rate cai', unbuffered ca in the er should rise at rate -cai'/0.15
i.e. -cai'/rhover
ENDCOMMENT
FUNCTION errel(cai (mM),caer (mM),kerrel (/s),kerm (mM)) (mM/s) { : from Huang
errel = kerrel*pow((cai/(cai+kerm)),4)*(caer-cai)
}
FUNCTION erfil(cai (mM),caer (mM),kerfila (mM/s),kerfilb (mM)) (mM/s) { : from John Rinzel 1997
erfil = kerfila * pow(cai/(1 (mM)),2) / ( pow(cai/(1 (mM)),2) + pow(kerfilb/(1 (mM)),2) )
}
FUNCTION erlek(cai (mM),caer(mM),kerlek (/s)) (mM/s) { : from Huang
erlek = kerlek*(caer-cai)
}
Comments from the top down:
1. Note the definitions of (mA), (um), FARADAY, and PI in the UNITS block.
2. I think the original units of kerrel, kerfila, and kerlek presume a fixed volume. In any
case, they don't make much sense in a model specification in which a single compartment
can have any arbitrary volume.
3. Note the use of a PARAMETER to specify the initial value of caer. This avoids
embedding "magic numbers" in the body of the program.
4. Note the use of diam. Local diameter is accessible from within a mod file.
5. Note the use of METHOD derivimplicit instead of cnexp.
6. The equation for ica must take surface/volume ratio into account. This is why it
involves diam.
7. Note the use of rhover (i.e. ratio of ER volume to cytoplasm volume) rather than Ver.
8. In the original mod file, the left hand and right hand sides (LHS and RHS) of the
FUNCTIONs were in inconsistent units. The functions should return values whose
units are concentration/time. Otherwise it would be necessary to factor in compartment
volume per micron length. You may have to go back to the formulas, lengths, areas, and
volumes used in the original models, in order to figure out whether the FUNCTIONs
need additional scale factors (e.g. division by actual volume used in an original model
that assumed fixed dimensions).
9. Notice the division by (1 (mM)) so that pow(cai/(1 (mM)), 2) etc. do not generate
"inconsistent units" error messages when tested by modlunit.
Although this mechanism now "works" it is entirely possible that the FUNCTIONs will
need adjustment (e.g. division by actual volume used in the original model) in order to
return numerically correct values.