Zero sum of currents

[somogyvari]

Hi!

We tried to calculate extracellular potential around a multicompartmental model of L5b neuron from Hay et al. 2011. We tried to collect all the membrane currents, but was unsure that we succeed, because the momentary sum of the currents was not zero. Than we turned to the LFPy, as the Hay model is a part of its tutorial, but the membrane currents still not sum to zero. We consulted with the LFPy guys, and they figured out, that if the synapses are located at the initial (0) or end (1) points of any segments than the zero sum stands, but if the alpha synapses are not located on either of the the ends of the segment, than the sum of all currents does not give zero, however they should.
Could you help with this problem?

Best regards,
Zoltan Somogyvari

[ted]

First let me thank you for using NEURON in your work, and also for asking an excellent question.

We tried to collect all the membrane currents . . . [but] the momentary sum of the currents was not zero

This statement raises several interesting issues about conservation of charge and current.

To start with, NEURON handles conservation of charge and current correctly, regardless of what some careless reader might infer from statements like this

We consulted with the LFPy guys, and they figured out, that if the synapses are located at the initial (0) or end (1) points of any segments than the zero sum stands, but if the alpha synapses are not located on either of the the ends of the segment, than the sum of all currents does not give zero, however they should.

I'm sure you will agree that if axial current is nonzero anywhere in a cell, then somewhere in that cell the net membrane current will be nonzero, at least for a while. And it follows that this is also true if a cell's membrane potential is nonuniform; actually it's true from the very instant at which membrane potential starts to deviate from being uniform over the surface of a cell.

Does "collect all the membrane currents" mean that you integrated all membrane currents over the entire surface of your model cell? Because that would be necessary to check for conservation of charge (and conservation of charge is what we're discussing).

Also, you had to deal with the fact that currents generated by point processes (synaptic currents, currents injected by electrodes) are in nA (nanoamps), but membrane capacitive current and currents generated by "density mechanisms" (like the fast sodium and delayed rectifier potassium currents, H current, voltage- or ligand-gated Ca currents) are in mA/cm2. And you had to deal with the fact that NEURON's sign convention for point processes is that synaptic currents that depolarize a cell are negative, but electrode currents (currents generated by an IClamp or SEClamp) that depolarize a cell are positive.

And you had to decide what to do about the currents generated by point processes attached to "zero area nodes" (the 0 and 1 ends of sections). If no other section is attached to the same location, just add the point process's current to the sum of currents of the segment at the 0 or 1 end (and remember to scale for units: multiply nA by 100 and divide by the segment's area in square microns, and the result is the numerical value of the current in mA/cm2). If you're dealing with a current source that is attached to the junction between two or more sections, and you are particularly interested in the extracellular potential in the vicinity of this junction, you might want to apportion the current among them according to some rule that you'll have to be ready to explain to reviewers. I'd be tempted to take a "current divider" approach based on the input impedances of the joined sections.

Here is a possible problem you may not have considered: when NEURON advances a simulation from t to t+dt using its fixed time step methods, it starts by making sure that the currents and other assigned variables (conductances etc.) are consistent with the state variables (including v) at time t, then calculates what the state variables should be at t + dt -- good so far -- but then TO SAVE TIME it doesn't update the assigned variables (conductances, currents etc.) to be consistent with the new state variables. Consequently the assigned variables are always "a bit behind" and if you try to calculate the net current at the new t, you'll be using the old conductances with the new membrane potentials, and the currents won't add up to 0. You COULD eliminate this problem by executing fcurrent() (h.fcurrent() if you're using Python); read about h.fcurrent() here https://nrn.readthedocs.io/en/latest/py ... l#fcurrent). Or you could eliminate the problem by using any of NEURON's adaptive integrators, because the states and assigned variables are consistent after each step of adaptive integration.

I have created a toy model that illustrates many of these points (a single compartment model cell with two AlphaSynapses--one attached to the 0 end and one attached to the 0.5 location of the soma--plus some informative graphs and a simple user interface for interactive usage) and demonstrates that NEURON handles conservation of charge correctly. I'll make it available online with a bit of documentation later today.

[ted]

Download https://www.neuron.yale.edu/ftp/ted/neu ... charge.zip, unzip it, read the contents of init.hoc, and then execute
nrngui init.hoc
and follow the instructions in init.hoc's header. This program confirms conservation of charge regardless of whether a point current source (synapse, IClamp, SEClamp or other point process) is located at the internal node or one of the zero area nodes of a single compartment model.

[somogyvari]

Dear Ted!

Thank you for your detailed answer, I really appreciate your help!
Regarding your question, yes, we tried to integrate to net membrane current on the whole surface of the neuron, so it is really about the charge conservation. I should definitely check all the points you raised (did we handle well the units, what type of intergrator we used etc.)
Now I realized, that I was wrong: the charge conservation worked for us for synapse positions between 0 and 1 but not for the two ends (I wrote the inverse).
Here you can find our simulation results and the discussion around:
https://github.com/LFPy/LFPy/issues/448
It seems that they were able to reproduce our problem. It may be related to the actual model...

I have a further question: I did not find a "universal" resistive (transmembrane) current variable in the NEURON, which would be the sum of all resistive currents. Or similarly, a net membrane current variable which would be the sum of all currents including capacitive currents as well. In case of a complex model (like the Hay model) it is hard to collect all types of currents just by reading their names (eg. Cai is a current or the internal calcium concentration). There are plenty of variables, so it is not easy to get through all the mechanisms... However, the NEURON knows, which variables stand on the right side of the equation, so it could tell us in one variable. Are there such universal current variables available?

Best regards,
Zoltán

[ted]

More good questions!

Thanks for the pointer to the simulation that was giving you problems. I'll take a look at it.

it is hard to collect all types of currents just by reading their names

True. But there are two ways to get

a net membrane current variable that is the sum of all membrane currents (ionic and capacitive)

One way is to insert extracellular, because that creates a range variable called i_membrane, so each segment has its own i_membrane variable. The only disadvantage is that this adds layers to the cable, and each layer has its own capacitance, which means more ODEs to be solved, and that consumes more memory and run time. But if your model already has extracellular, the i_membrane range variable (units mA/cm2) already exists, so why not use it? Example:
soma.i_membrane(0.5) // hoc
or
soma(0.5).i_membrane # Python
reports the membrane current density (in mA/cm2) for the segment that contains the midpoint of the soma. To convert current density in mA/cm2 to total membrane current in nA, multiply by segment area in um2 and divide by 100. Example:
soma i = area(0.5) * i_membrane(0.5) / 100 // hoc
or
i = soma(0.5).area() * soma(0.5).i_membrane / 100 # Python
and the units of i will be nA.

Note that i_membrane does not include contributions from point processes that generate an ELECTRODE_CURRENT, such as IClamp, SEClamp, or VClamp. I don't recall if i_membrane_ does that too; it would be easy to find out by running a test on a simple model.

The other way is to invoke the CVode class's use_fast_imem()--just execute
cvode.use_fast_imem(1) // hoc
or
h.cvode.use_fast_imem(True) # Python
Then
soma.i_membrane_(0.5) // hoc
or
soma(0.5).i_membrane_ # Python
reports the total membrane current (in nA!!) for the segment that contains the midpoint of the soma.

[somogyvari]

Dear Ted,

Thank you for the idea of inserting extracellular to get the total membrane current. Implicitly, we did this by using LFPy to get the membrane current. However, the results were strange in that sense, that changing only the location of the synapses changed the sum of currents, as it is described in the above github discussion.

Best regards,
Zoltán