Looks like a case of egregious complexity to me, but who am I to judge?
For the moment let us assume that the model specification downloadable from
http://www.cellml.org is complete and correct (chances are slim to none for this being
the case with the specification in the article itself--there are just too many parameters and
equations with fancy notation to escape the inevitable typo).
Regardless of what simulation environment you use, the way to implement this model is
to start small and add complexities incrementally, testing at every step. It doesn't matter
how smart you are, how facile you may be with programming, or how hard you work--
frustration and debugging nightmares are guaranteed if you try to do it all in one pass.
Start by ignoring all the ion accumulation stuff. Do the basic spike currents first:
1. the fast ina
2. pick whichever k current seems most likely to terminate the spike
3. the L current
For the moment ignore complexities such as eq. A41 (ENa depending on {K] as well as
[Na]). You can add that later. At this point, just get the voltage-dependence of the gating
right. Do voltage clamp experiments on each voltage-gated channel and compare the
results with the figures in the paper. Fix any discrepancy you discover.
Add the other ionic currents one at a time. Give each one its own mod file. Make your
life easy by using economical nomenclature--e.g. there are a bunch of K currents, each
of which will have a descriptive suffix. So for IKs,
Code: Select all
NEURON {
SUFFIX ks
USEION k READ ek WRITE ik
RANGE gbar, g, i
. . .
BREAKPOINT {
. . .
g = gbar*etc.
i = g*(v - ek)
ik = i
}
lets you use gbar_ks, g_ks, i_ks instead of the slightly redundant gkbar_ks, gk_ks,ik_ks.
The authors assume specific membrane capacitance of 1 uf/cm2, so their conductance
densities in mS/uF are numerically equal to mS/cm2. You can use the same units for
gbar and g, but be sure to use modlunit to check for units consistency (if your currents
are in mA/cm2, it will tell you that you need to write i = (0.001)*g*(v-ek) ). Comparing
your current plots to theirs will require scaling. Figure 16 uses pA/pF, which is the same
as uA/cm2, i.e. numerically 10^3 times bigger than NEURON's default units for current
density (mA/cm2).
If there are Ca-dependent ion channels, you will need to add them before implementing
a Ca accumulation mechanism. Try to figure out what the "mean" cai and cao are
supposed to be, and assign the values to cao_ca_ion and cai_ca_ion in hoc. Just
remember that eventually you will have to eliminate these assignments (well, at least
the assignment of cai_ca_ion).
Once you have the ionic currents working, it will be time to add ion accumulation.
Electrically, this is a single compartment model, so you'll be using a section with
nseg = 1 and must choose an L and diam to give the surface area and "myoplasmic
volume" specified in table 2.
Add Na and K accumulation first. This can be done with a single mod file that
READs ina and ik, and WRITEs nai and ki (I assume that there is no extracellular
ion accumulation). Then add an Na-K pump. This will be in another mod file that READs
nai and ki (maybe also nao and ko) and WRITEs ina and ik (the pump's fluxes).
Now it makes sense to incorporate the mixed ionic dependence of the driving force for
ina. To that mod file add
Code: Select all
FUNCTION enax() (/mv) {
enax = (R*(celsius+273.15)/FARADAY) * log((0.9*nao + 0.1*ko)/(0.9*nak + 0.1*ki))
}
(I may have omitted a conversion factor here--check with modlunit!),
and in the BREAKPOINT block change
to
Code: Select all
: i = g*(v - ena)
i = g*(v - enax())
Don't forget that this mechanism must now READ nai, nao, ki, and ko, and that you must
be sure to propertly initialize all of these concentrations.
Finally begin to work on Ca accumulation. The first thing to do here is to take care of
ordinary buffering and the Ca pump. You'll have to read the paper to figure out what
all those arrows and boxes mean, and determine what is the basic buffering equation.
Then use one of the ca*mod files that NEURON comes with as a starting point, and
change it accordingly (under MSWin they're in c/nrn58/examples/nrniv/nmodl/ , and on
my Linux box they're in /usr/local/src/nrn-5.9/share/examples/nrniv/nmodl/ ). Be sure
to test the pump and buffering to verify that it works the same as the one described in
the paper. This mechanism will WRITE cai.
Then add the NaCa exchange in its own mod file. This will READ cai and WRITE ina and
ica (the exchanger's fluxes).
Finally tackle the sarcoplasmic reticulum. This can be done in part by elaborating the
basic Ca accumulation mechanism. However, it appears from Fig 1 that icaL enters a
confined region where it accumulates as "Ca_ss" and is subject to exchange with free
Ca in the cytosol and the JSR. Is this interpretation correct?