modeling ER calcium store

[along82]

Hi, ted
Thank you for your reply, and I need your suggestion for another issue:
Adding a single intracellular endoplasmic recticulum based on a single compartment neuron model to simulate an calcium store.
My idea now is to calculate a current component and add it to the ica current when calculating cai, and subtrate the volume of the ER from the neuron volume. This seems an ugly way, so what is your perfect suggestion?

[ted]

You have the basic idea. Here are some tips.

1. Specify the properties of this mechanism in its own mod file.

2. This should be a density mechanism, so its NEURON block should say something like
SUFFIX ER

3. It must change cai so it must have a
USEION ca . . . WRITE ica
statement. But it should be electrically neutral (contribute 0 to membrane current), so it
must also have a
NONSPECIFIC_CURRENT ix
statement. The BREAKPOINT block will contain a statement that calculates ica; it must also
contain this statement
ix = -ica

subtrate the volume of the ER from the neuron volume

Every ion accumulation mechanism contains a statement in which the net ion flux is
divided by compartment volume (or a statement that does something equivalent to this).
If the ER volume is a significant fraction k of total volume, the calcium accumulation
mechanism should be revised so that net ca flux is divided by
(compartment volume * (1 - k)).

[along82]

Hi, ted!
Actually I have recently sent u an email but probably you didn't get it. Forget about that email and now I state my problem here again.
Now I've got a strange problem on my codes. With these code, when I tried to draw the Voltage graph, current graph, state graph or whatever, the line just go back to the minus zone of the x-axis which is actually Time.
I've already met this problem in another mod file, but I fixed that one. Problem in that file was likely in the FUNCTION block, where I mistakely used a global variable (like FARADAY, Vcyt) within the FUNCTION block without passing it through the function argument. I'm not sure if I've got the point but the problem vanished in that mod file anyway.
But I still can't get out the mistakes in this attached mod file, seems there are more one mistakes within this simple mod file, some straight forward, some implicit, but all beyond my vision.
The MOD file

Code: Select all

TITLE Endoplasmic Reticulum

UNITS {
	(mol) = (1/liter)
	(mM) = (milli/liter)
}

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     (liter/s): from Huang, to be adjusted
	kerfila = 0.75e-12 (mM*liter/s) : from John Rinzel 1997
	kerfilb = 0.2e-3 (mM) : from John Rinzel 1997
	kerlek = 6.15e-14 (liter/s): from Huang, to be adjusted
}

ASSIGNED {
	FARADAY
	Ver                 : ER volume
	ica (mA/cm2)
	ix (mA/cm2)
}

BREAKPOINT {
	SOLVE store METHOD cnexp
	ica = -2*FARADAY*(errel(cai,caer,kerrel,kerm)-erfil(cai,caer,kerfila,kerfilb)+erlek(cai,caer,kerlek))*1e-3
	ix = -ica
}

INITIAL {
	caer = 0.25 : initial calcium concentration in ER
}

DERIVATIVE store {  
	caer' = -(errel(cai,caer,kerrel,kerm)-erfil(cai,caer,kerfila,kerfilb)+erlek(cai,caer,kerlek))/(Ver/fer)
}

FUNCTION errel(cai (mM),caer (mM),kerrel (liter/s),kerm (mM)) (mM*liter/s) { : from Huang 
	errel = kerrel*pow((cai/(cai+kerm)),4)*(caer-cai)
}

FUNCTION erfil(cai (mM),caer (mM),kerfila (mM*liter/s),kerfilb (mM)) (mM*liter/s) { : from John Rinzel 1997
	erfil = kerfila*(pow(cai,2)/(pow(cai,2)+pow(kerfilb,2)))
}

FUNCTION erlek(cai (mM),caer(mM),kerlek (liter/s)) (mM*liter/s) { : from Huang
	erlek = kerlek*(caer-cai)
}

And the HOC file

Code: Select all

load_file("nrngui.hoc")

create soma
access soma

soma {
	nseg = 1
	diam = 20                           // um
	L = 22.5                            // um
	Ra = 123.0
	Area = PI*diam*L*(1e-8)             // Cell surface (cm2)
	Vcell = PI*diam*diam/4*L*(1e-15)    // Cell volume (liter)
	Vcyt = Vcell/(1+0.15)               // cytoplasmic volume
        Ver = 0.15*Vcyt                     // ER volume

//-- ER calcium current
        insert er
     }
	
objectvar stim
soma stim = new SEClamp(0.5)

stim.dur1 = 100
stim.amp1 = -65
stim.dur2 = 100
stim.amp2 = 40
stim.dur3 = 100
stim.amp3 = -65
stim.rs = 1e-5 // actually 10 ohm

tstop = 800

[ted]

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.