simple internal calcium concentration

[gabrielleg]

Hi,

I've been trying to get my code for internal calcium concentration to work. It should be really simple, but for some reason, it's not working. Even though there is plenty of inward calcium current, cai won't budge from the C0 value (which is kind of like cai0, but not; it's the "background [ca]"). C0=0.5e-3 mM. Cai is initialized to the same value as C0. There is no error when compiling or when running hoc code that uses this mechanism.

I have checked the units and conversion factors. I did a little test in MATLAB to see if the equation was the problem. In MATLAB cai changes, but I used a sin wave to approximate ica.

Something I noticed is that cai = C0, not cai0, at initialization if cai0 and C0 are set to different things, but this makes sense if initialization brings cai to steady state and not to the initial value. So what's wrong here?

Code: Select all

NEURON
  {
  SUFFIX cainT
  USEION ca READ ica WRITE cai
  RANGE tau, F, L, C0
  }

 UNITS
  {
  (mA) = (milliamp)
  (mV) = (millivolt)
  (mol) = (1)
  (molar) = (mol/liter)
  (mM) = (millimolar)
  (S) = (mho)
  (uS) = (microS)
  (mJ) = (millijoule)
  (um) = (micrometer)
  pi = (pi) (1)
  }

PARAMETER  
  {
  F = 0 (mM/mA)
  C0 = 0 (mM)
  }

ASSIGNED
  {
  ica (mA/cm2)
  tau (ms)
  L (um)
  diam (um)
  }

STATE  
  {
  cai (mM)
  }

INITIAL
  {
  cai = cai0
  }

  BREAKPOINT
  {
  SOLVE state_change METHOD cnexp
  }

DERIVATIVE state_change
  {
  cai' = (-(F*ica*pi*diam*L*(1e-8))+C0-cai)/tau
   : F is not Faraday
  }

gabrielle

[ted]

what's wrong here?

Several things, some trivial, some potentially very serious, but all fixable. Let's start with
the trivial.

cai won't budge from the C0 value

because F is 0 (unless you assign it a value during model setup).

tau is also 0, and so is L. This will cause lots of problems because cai will become
"nan" ("not a number"). Are you using hoc statements to assign nonzero values to
these?

The following are more important.

1. There is confusion about the initialization of cai. This is a potentially more serious
problem than tau or L being 0, because cai will have a "real" numerical value . . . but will it
be what the user expects? To remove any ambiguiity about the initialization of cai, add
RANGE cai0
to the NEURON block, and an explicit assignment of a value to cai0 in the PARAMETER
block, e.g.

Code: Select all

cai0 = 0.5e-3 (mM) : if that's the value you want

Then you have control over both the "background" and the iniital value of ca.

2. For several reasons, this is the most problematical aspect of the code:
cai' = (-(F*ica*pi*diam*L*(1e-8))+C0-cai)/tau
This looks like an attempted port to NEURON of something that was previously implemented
for some other simulator or programming language. I say that because it contains what
looks like an explicit calculation of the volume in which ca is distributed. It happens that
the volume is then used incorrectly--as a multiplicative factor, so that increasing
compartment volume _increases_ the effect of ica on cai'--but even if done correctly,
this is a bad idea because it defeats NEURON's strategy for flexible management of
longitudinal spatial discretization. Ignoring for the moment the important issue of units, ica
should be multiplied by a factor that expresses
(surface area)/(volume)
which for a cylinder is simply 4*pi/diam
modlunit will tell you what units your F factor must have.

After making the requisite changes, it would be good to tidy up by eliminating L from the
NEURON and ASSIGNED blocks.

[gabrielleg]

Ok. So the trivial things aren't the problem here. F, tau, L, and C0 are all assigned in hoc code.

1.
cai0 was also set in hoc code (actually it's set in a template code, but this is almost the same thing). I tried setting it explicitly in the mod code like you said and this did something strange. The voltage trace is a flat line whereas it used to have some spike-like behavior before. This is weird because cai still initializes to C0 and stays put even though cai0 is explicitly set in the mod file. But cai initializes to cai0 if I use the run control panel which uses the init() function. This makes sense because my initialization code is meant to bring things to steady state.

Code: Select all

proc initialize() {
	finitialize(v_init)
	t = -1e5
	dt = dt1
	while (t<-dt/2) {
		fadvance()
		}

	fcurrent()
	frecord_init()
}

2.
I suspected that the main equation was the culprit. I am confused about the units to use here though. I'm trying to implement someone else's model where they have calcium concentration governed by this:

tau*d[ca2+]/dt = -F*ica - [ca2+] + C0
F is a factor that they gave in units of uM/nA. I set it in hoc code in terms of mM/mA.

So in my code I had:

Code: Select all

cai' = (-F*ica+C0-cai)/tau 

But the units didn't check out and since ica is in mA/cm2, I multiplied it by the surface area to keep cai in terms of mM:
cai' = (-(F*ica*pi*diam*L*(1e-8))+C0-cai)/tau
Are you saying that cai actually has units of mM/cm3? Or that F needs to be multiplied by volume while ica is multiplied by surfArea/volume?

p.s. thanks for your help.

[ted]

cai0 was also set in hoc code

Doesn't matter. The cai0 used in your mod file is not visible to hoc. Nothing you do to a
variable named cai0 at the hoc level will have any effect on the cai0 in your mod file. You
can repair that disconnection by declaring
RANGE cai0
in the NEURON block, and including a statement of the form
cai0 = somevalue (milli/liter)
in the PARAMETER block, where somevalue is some value > 0 (it is not a good idea
for a default initial ionic concentration to be 0, because the resulting driving force will be
nonsense). Then you will be able to control the mechanism's initial cai by hoc statements
of the form
cai0_cainT = somevalue

my initialization code is meant to bring things to steady state.

proc initialize() is seriously broken. It doesn't restore t to 0, or dt to its original value.
It won't work if adaptive integration is being used. If it is based on something you found
in NEURON's documentation, please let me know so we can fix that.

I'm trying to implement someone else's model where they have calcium concentration governed by this:

tau*d[ca2+]/dt = -F*ica - [ca2+] + C0
F is a factor that they gave in units of uM/nA.

The original model was almost certainly not implemented with NEURON. One often sees
such code in models that assume a fixed geometry. Such implementations almost always
specify membrane currents in absolute units (net current) rather than current density
(e.g. mA/cm2), and channel conductances as net conductance or resistance (siemens
or ohms) rather than conductance density (siemens/cm2). To port such a model to
NEURON, you have to first discover the surface area and volume of the original model,
then choose a cylindrical geometry with the same surface/volume ratio. It is also
necessary to use this information to recast all conductances and currents in density
units.

But the units didn't check out and since ica is in mA/cm2, I multiplied it by the surface area to keep cai in terms of mM:
cai' = (-(F*ica*pi*diam*L*(1e-8))+C0-cai)/tau

Now I see what you're doing. My mistake--I was thinking that there had been a couple
of typos--that the original had been of the form F*ica/(pi*L*diam^2). Wrong way to go.

Time for a fresh start. Can you point me to the definitive statement of the original model?

[gabrielleg]

Thanks for your reply. My initialization is based on something I found in the NEURON book, section 8.4.2. Writing initialization routines in NEURON is still a mystery to me. I haven't been able to get a handle on the details of how they work. In particular, changing things like the size of dt and the t expression in the while loop can sometimes give vastly different results.

As for the original model, they do specify conductances, resistances, and capacitances in absolute units. I had to use the total capacitance and axial resistance to extract the geometry. Here is the original statement for calcium concentration:

tau*d[ca]/dt = -F*Ica - [ca] + C0

This comes from Soto-Trevino, et al, 2005. Here's a link to it.

http://jn.physiology.org/cgi/content/abstract/94/1/590

Thanks for your help!

g

[ted]

Please discard your custom initialization. Put aside all consideration of custom initialization
for the moment. Focus just on making a correct implementation of the anatomical and
biophysical specification of the model.

The first step in creating a NEURON implementation of the calcium accumulation
mechanism described by the authors is to determine the volumes and surface areas of the
compartments that contain this mechanism. Let us focus on the soma of the AB cell; you
will then be able to deal with the soma of the PD cell in a similar manner.

Surface area is easy. Start by assuming that specific membrane capacitance is 1 uf/cm2,
which means that 100 um2 of membrane has total capacitance 1 pf = 1e-3 nf.
AB soma has total capacitance = 9 nf, i.e. 9000 times larger.
Therefore AB soma surface area is 9e5 um2.

Volume is easy too.
In the absence of buffering or transport, 1 A of Ca current into a volume of 1 liter makes
[Ca] change at a rate of (1/2/FARADAY)) Molar/s = (1/2/FARADAY) mM/ms.
So if Ca current is 1 nA and volume is 1 um3, the rate is (1/2/FARADAY)*(1e6) mM/ms
(note that 1 liter = 1e15 um3).
And if the volume in um3 is vol, the rate of change is 1e6/(2*FARADAY*vol) mM/ms.

In the authors' formula tauca Cai' = -F*ICa - Cai + C0
the 2nd and 3rd terms on the RHS represent the effects of buffering and transport.
Ignoring those, tauca Cai' = -F*ICa
where tauca is in ms, F is in uM/nA, and ICa is in nA
A net Ca current of 1 nA makes Cai change at a rate of F/tauca (uM/ms)
= (1e-3) F/tauca (mM/ms)
so
1e6/(2*FARADAY*vol) = (1e-3) F/tauca
and vol = (1e9)*tauca/(2*F*FARADAY)

For a cylinder, volume = PI*L(diam/2)^2 and surface area = PI*L*diam so diam = 4*volume/area

Plugging in numbers we get AB soma area = 900000 um2, volume = 3756428.8 um3,
diam = 16.695239 um, and L = 17159.317 um.
So now we have the physical dimensions of the soma compartment.

What about the details of the Ca accumulation mechanism itself? This will do:

Code: Select all

COMMENT
NEURON implementation of calcium accumulation mechanism described in
Soto-Trevino et al.
J Neurophysiol 94: 590 – 604, 2005.
ENDCOMMENT

NEURON {
  SUFFIX cacumst
  USEION ca READ ica WRITE cai
  RANGE tau, C0
}

UNITS {
  (mA) = (milliamp)
  (mol) = (1)
  (molar) = (mol/liter)
  (mM) = (millimolar)
  (uM) = (micromolar)
  (um) = (micrometer)
  FARADAY = (faraday) (coulombs)
}

PARAMETER {
  tau = 303 (ms) : AB soma; 300 ms for PD soma
  C0 = 0.5 (uM)
}

ASSIGNED {
  ica (mA/cm2)
  diam (um)
}

STATE {
  cai (mM)
}

INITIAL {
  cai = (1e-3)*C0
}

BREAKPOINT {
  SOLVE states METHOD cnexp
}

COMMENT
Surface area of a cyl of length len is PI*len*diam
so net ca influx is -(ica/(2*FARADAY))*PI*len*diam
Volume is PI*len*diam^2/4
so in the absence of buffering or transport
rate of change of conc is -(ica/(2*FARADAY))*4/diam
ENDCOMMENT

DERIVATIVE states {
  cai' = -(1e4)*2*ica/(FARADAY*diam) + (1e-3)*(C0 - (1e3)*cai)/tau
}

Note the numerous scale factors (constant multipliers inside paired parentheses)
which are needed to reconcile the authors' units with NEURON's preferred units.
Also note carefully the form of the first expression on the RHS of the statement in the
DERIVATIVE block, and its explanation in the preceding COMMENT block.

The following program, which I call testAB.hoc, tests this mechanism:

Code: Select all

load_file("nrngui.hoc")

create soma
access soma

soma {
  diam = 16.695
  L = 17159
  insert cacumst
  tau_cacumst = 1e9 // eliminate transport and buffer for this test
}

/*
F = 0.418 uM/nA
tauca = 303 ms
so 1 nA Ca influx x 303 ms should make cai increase by 0.418 uM
ditto for 303 nA Ca influx x 1 ms
*/

objref capp
capp = new CaPP(0.5) // delivers specified Ca current
capp.del = 1 // ms
capp.dur = 1 // ms
capp.amp = -303 // nA

finitialize(-65)
c0 = cai
print "Initial cai ", c0, "mM"
run()
print "Final cai ", cai, "mM"
print "Difference is ", 1e3*(cai-c0), "uM"

print "Soma capacitance is ", area(0.5)*cm*1e-5, "nF"

The physical dimensions of the soma compartment are those that were derived above
from the authors' model description, rounded off to 4 place precision. Making
tau_cacumst very large reduces the effect of buffer and transport to a negligible value.
CaPP is a custom point process that delivers a specified calcium current over a
specified time interval. It is patterned after the code in stim.mod, which defines the
IClamp class (see nrn/src/nrnoc/stim.mod).

Code: Select all

NEURON {
	POINT_PROCESS CaPP
	USEION ca WRITE ica
	RANGE del, dur, amp
}

UNITS {
	(nA) = (nanoamp)
}

PARAMETER {
	del (ms)
	dur (ms)
	amp (nA)	: negative is inward
}

ASSIGNED {
	ica (nA)
}

BREAKPOINT {
	at_time(del)
	at_time(del + dur)

	if (t > del && t < del + dur) {
		ica = amp
	}else{
		ica = 0
	}
}

Executing testAB.hoc generates this output:
Initial cai 0.0005 mM
Final cai 0.00091801976 mM
Difference is 0.41801976 uM
Soma capacitance is 8.9997049 nF
Note that total soma capacitance and the change of cai are both as expected.

[figoyouwei]

Hi Ted,

For this topic, I understand that one can initialize ion concentration in the mod file. In calcium case, e.g. parameter it <cai0 = some value> and <cao0 = soma value>.

It works fine, but that means if there are ten different calcium mod files, one has to specify the parameters in all these ten files ?

I noticed that some public models defined a <cai0_ca_ion = some value> or <cao0_ca_ion = some value> in the "hoc" file,
my question is:
I don't understand how this hoc variable is called by all mod files ? Could you please elucidate ? thanks !

[ted]

if there are ten different calcium mod files, one has to specify the parameters in all these ten files ?

Only if you aren't happy with the default values of the parameters.

How to specify initial concentrations of ions is discussed in chapter 8 of The NEURON Book. If you don't have the book, get this early version of chapter 8:
https://www.neuron.yale.edu/ftp/ted/boo ... xedref.pdf

You'll also want to read this in the Programmer's Reference:
http://www.neuron.yale.edu/neuron/stati ... #ion_style

[figoyouwei]

thanks Ted, I got it now, it works :) and I appreciated your sharing of resource.