Transmembrane currents do not sum to zero

[joseph-tharayil]

Context

We are injecting current into neurons with SEClamp, by playing a random variable into rs. By Kirchoff's law the sum of the transmembrane currents over all segments should be equal to the SEClamp current at each timestep

Overview of the issue

The sum of the transmembrane currents and the SEClamp current is not equal to zero when a random vector is played into the SEClamp conductance with interpolation (i.e., random_vector.play(seclamp._ref_rs,time_vec,1) ). However, the currents do sum to zero when the values are played in without interpolation. This problem seems to be unique to the SEClamp.

Expected result/behavior

The sum of the membrane currents over the neuron and the SEClamp current should be zero

NEURON setup

• Version: 9.0a
• Installed as a module on a CentOS Linux cluster

Minimal working example - MWE

testSEClamp.zip

Logs

See attached notebook

[nrnhines]

Yes. The total extracellular current should be equal to the stimulus current. That seems to be the case, with reasonable numerical accuracy, when I run your model with the variable time step method activated. The above plot is for totalCurrent where totalCurrent = membraneCurrents - svec. Note: I had technical difficulties with plotting using your jupyter notebook and worked around them by translating the ipynb file to a stand-along python file and modifying a copy of that to allow multiple runs with the InterViews GUI.

jupyter nbconvert --to=python testSEClamp_morph.ipynb

However, running with fixed step method, gives results that don't satisfy Kirchoff's law and exhibit several qualitative features that have to be diagnosed. The upper graph shows a very large current imbalance with a magnitude of -950 and the zoomed in lower graph shows way too large current imbalances over the region of 1ms piecewise linear random rs changes when the SEClamp reference potential is 0mV. Note: to get the same results with multiple runs, one must additionally initialize with h.SEClamp[0].rs=1

Turns out, given the existence of current imbalances, the outlier is easy to understand. At 162ms the random value of rs is the minimum of the s vector 0.000481829 and that is a small enough resistance to make the SEClamp almost ideal and therefore cause an almost unlimited current to force soma(.5).v to be close to SEClamp.vc. The next smallest value in the time range of 100-200 is at 135ms with value 0.01588 which would generate 100 fold less current for a given v-vc than the outlier.

My initial hypothesis as to why the current imbalance occurs was incorrect and I need to look more deeply. That is, it turns out that

vcvec = h.Vector()
rsvec = h.Vector()
vcvec.record(stim._ref_vc)
rsvec.record(stim._ref_rs)

# run

zz = vcvec.c().sub(vec).div(rsvec)

generates a zz which is identical to svec in the domain 100-200. I will get back to you when I know more.

[joseph-tharayil]

Dear Dr. Hines,

Thank you for your insight. I just wanted to follow up to see if you had been able to figure out why the error occurs. For our purposes, using a variable-timestep method won't work, since we need to use coreneuron.

Best, Joe

[nrnhines]

Unfortunately the combination of vacation and work on other issues has so far prevented developing a detailed explanation and general recipe for summation of currents to exhibit that charge conservation is satisfied to the expected order of accuracy. I hope to have something by the end of the coming weekend.

[nrnhines]

I carried out all my diagnostic tests with a single compartment model with hh channel inserted. #2733 fixes a bug in the calculation of i_membrane_ when an electode current depends on time varying conductance. With that fix, i_membrane_ and i_membrane * area * 0.01 are now the same in all cases to within accumulated round_off error for a few arithmetic operations (approximately the 1e-14 level).

With respect to the issue of SEClamp.i + total i_membrane_ not summing to 0, the problem is that in the process of fixed step integration of the current balance equations from t to t + dt, current balance is solved at t + dt/2 using the matrix equation G(t+dt/2, v(t)) * (delta_v) = I(t+dt2, v(t)) where G is somewhat analogous to conductance, I is somewhat analogous to current, and delta_v is the solution to that equation. i_membrane and v end up being calculated at t + dt. Vector recording of v is correct to order dt when secondorder=0 and (usually) order dt^2 when secondorder=1. However the value of SEClamp.i that makes it into the Vector.record element at time t, though first order correct, is based on the value of v(t-dt) and SEClamp.rs(t - dt/2) . This is in contrast to the effect of SEClamp on voltage which is rs(t+dt/2)*(vc(t+dt) - v(t + dt)). It is this latter when summed into total i_membrane that sums to 0.

One can see this with (stim is an instance of SEClamp)

tvec = h.Vector().record(h._ref_t)
vvec = h.Vector().record(s(0.5)._ref_v)
vstimvec = h.Vector().record(stim._ref_vc)
istimvec = h.Vector().record(stim._ref_i)
rsstimvec = h.Vector().record(stim._ref_rs)

def fstim():
    rs = rsstimvec.c().add(rsstimvec.c().rotate(1)).mul(0.5)
    rs[0] = rsstimvec[0]
    vc = vstimvec.c()
    s = vc.sub(vvec).div(rs)
    return s

# after a run
plt("imemvec.c().sub(fstim())")

That yields a balance on the order of 1e-8 for the discontinuity in vc at 1ms and 1e-10 while rs is changing:

in contrast to

plt("imemvec.c().sub(istimvec)")

In summary, all the states and currents are first order correct with secondorder=0 and the error is proportional to dt. However, internally, the current balance equations are satisfied to high accuracy. It is just difficult to see that kind of accuracy without taking into account not only the electrode current variable but the change in electrode current with respect to voltage in order to reconstruct exactly how electrode current affects the total current balance.

My next step is to apply this overall strategy to your above model. (My test model cannot use secondorder=2 since my SEClamp is too stiff and the closeness to numerical instability causes ringing.)

[nrnhines]

test1.py.txt is a modified version of testSEClamp_morph.ipynb that was in your above zip file. I use it with python -i test1.py.txt from the terminal and it runs the model and plots (InterViews GUI) totalCurrent (ideally 0.0 everywhere). In the file, I manually change h.secondorder to be 0, 1, or 2; cv.active(0 or 1), and coreneuron.enable = True or False. Can't have both coreneuron enabled and cvode active. I'm using cv.use_fast_imem(1) and extracellular commented out. (So at the moment, for correct results, one needs to use the hines/fast-imem-fix branch of [email protected]:neuronsimulator/nrn) . Note that extracellular does not work with CoreNEURON ; and with cvode active, it can't work with the CVode ODE solver but switches to the differential algebraic solver. I configured with

cmake .. -G Ninja -DCMAKE_INSTALL_PREFIX=install -DPYTHON_EXECUTABLE=`which python` -DNRN_ENABLE_CORENEURON=ON

For the calculation of totalCurrent I use

    totalCurrent = h.Vector(membraneCurrents - fstim())

where

def fstim():
    if h.cvode.active():
        return svec
    # SEClamp conductance is evaluated (Vector.play) at t + .5*dt
    # But recording is at time t + dt.
    # So use rsstimvec to evaluate clamp resistance at midpoint
    # of each step. That is the resistance for the entire step.
    rs = rsvec.c().add(rsvec.c().rotate(1)).mul(0.5)
    rs[0] = rsvec[0]
    vc = vcvec.c()
    v = vec.c()
    if h.secondorder:
        v = vec.c().add(vec.c().rotate(1)).mul(0.5)
        v[0] = vec[0]

    s = vc.sub(v).div(rs)
    # but nonzero only when SEClamp is on (stim.i != 0)
    for i, x in enumerate(svec):
        if (x == 0.0):
            s[i] = 0.0
    return s

Results for h.secondorder = 0 are

Results for h.secondorder=2 are

Results for h.secondorder=2 and coreneuron.enable = True are

And results for cv.active(1) are

These results differ in part due to different randomness from

randomvals = np.random.rand(300) # Random values for SEClamp rs

since

$ python -c 'import numpy; print(numpy.random.rand(5))'
[0.92735562 0.79503546 0.71969147 0.43453152 0.09363631]
$ python -c 'import numpy; print(numpy.random.rand(5))'
[0.71812691 0.35137599 0.35065615 0.16233652 0.9573502 ]

[jorblancoa]

Hi @nrnhines, I've been collaborating with @joseph-tharayil on integrating online LFP calculations into CoreNEURON. During this process, we encountered an issue where the currents did not sum to zero.

From my understanding, due to the inability to utilize cvode in CoreNEURON, the workaround involves substituting SEClamp with the fstim() function. My query is: Is it feasible to incorporate fstim() into NEURON, maybe within the svclmp.mod file located at https://github.com/neuronsimulator/nrn/blob/c7d200c99e959e2675a38fab7e414b1cc3b227ee/src/nrnoc/svclmp.mod? This would allow for the retrieval of updated fstim() values via the SEClamp mechanism pointer here: https://github.com/neuronsimulator/nrn/blob/c7d200c99e959e2675a38fab7e414b1cc3b227ee/src/coreneuron/io/reports/report_handler.cpp#L255.

Thanks! Jorge

[nrnhines]

workaround involves substituting SEClamp with the fstim() function.

I'm missing something. My belief is that: Electrode currents do not appear in LFP calculations. The fstim value is no more accurate than SEClamp.i, it is just a time adjusted value to correspond to the time at which i_membrane_ is calculated so that the sum is 0. I regret that the currents calculated in a BREAKPOINT block are left in an intermediate state on return from fadvance and it is difficult to adjust them to the fadvance exit value without calling the relatively expensive h.fcurrent() function (which would have to be called prior to filling in the vector.record elements).

[jorblancoa]

Probably is me the one missing something... So right now, your fix here should be enough to fix the summation reports of i_membrane (nt.nrn_fast_imem->nrn_sav_rhs) and SEClamp(get_var_location_from_var_name(mech_id, var_name.data(), ml, j);) not summing 0 in CoreNEURON? https://github.com/neuronsimulator/nrn/blob/master/src/coreneuron/io/reports/report_handler.cpp#L220

[nrnhines]

That is a good question. The fix will cause i_membrane_ to be the same as extracellular i_membrane * area * .01 . If that suffices to fix the summation report, then ...

[jorblancoa]

After trying your fix, I noticed a small change in the summation reports compared with the previous version...however really small and not even close to 0. lets take the first timestep: current - 0.00158077 vs 0.00498575 - old and with cvode -3.73781e-14

I am out of ideas, maybe trying to implement the fstim() method every timestep before writing the report directly in coreneuron? However @nrnhines you mentioned before that this wont change IClamp value at all? cc: @pramodk @joseph-tharayil

[nrnhines]

@jorblancoa I've lost the context. I don't know how to reproduce those numbers. I've generally been running either the testSEClamp.zip or a single compartment model based on

s = h.Section("s")
s.L = 10
s.diam = 10
s.nseg = 1
s.insert("hh")

with the addition of an SECLamp with varying vc and rs. Somehow I need to get onto the same page you are on. Do you have a zip file that generates your numbers?

[jorblancoa]

The numbers are just to show the change between versions. I got them by running a small 1 cell neurodamus simulation that generates a summation report adding nt.nrn_fast_imem->nrn_sav_rhs and the IClamp. I am not sure I have the expertise to create a small example with plain neuron...

[nrnhines]

Ok. Perhaps we will eventually have to diagnose the neurodamus simulation. But before that, it might be useful to run the new test that is a part of the just merged pull request. It is nrn/test/hoctests/tests/test_imem.py and runs a variety of SEClamp and IClamp tests to verify that i_membrane_ - electrode current == 0 By modifying that file so that

diff --git a/test/pytest_coreneuron/test_imem.py b/test/pytest_coreneuron/test_imem.py
index 902dc6178..7cdc7294f 100644
--- a/test/pytest_coreneuron/test_imem.py
+++ b/test/pytest_coreneuron/test_imem.py
@@ -100,7 +100,7 @@ def tst(xtrcell, secondorder, cvactive, corenrn, seclmp):
         return s
 
     def plt(vec, vecstr):
-        return  # comment out if want to see plots
+#        return  # comment out if want to see plots
         g = h.Graph()
         glist.append(g)
         vec.label(vecstr)
@@ -168,4 +168,5 @@ def test_imem():
 
 
 if __name__ == "__main__":
-    test_imem()
+    glist = tst(0, 2, 0, 0, 0)
+

and running with

$ python3 -i test_imem.py

A number of graphs pop up. The tst(0, 2, 0, 0, 0) refers to

No extracellular
secondorder = 2
cvode is not active
coreneuron is not active
using IClamp as the electrode

First, the i_membrane_ - stim_i is very close to 0.0

Second, the IClamp stimulus (stim.i) current is time varying in the form of a triangle wave and elicits an "identical" i_membrane

Third, the voltage response is

Note that for an IClamp electrode, I directly subtract the stim_i```` from the i_membrane_and do not need the adjusted current computed bydef fstim()```

[pramodk]

We can investigate the neurodamus+coreneuron issue separately.

@joseph-tharayil: Maybe you can better understand Michael's response + changes and see if that addresses issue in your original reproducer?

[joseph-tharayil]

My understanding of the current situation is as follows:

Both the fix in #2733 and updating the SEClamp current values with fstim() are necessary in order to ensure that currents sum to zero at each time step (to clarify, in our use case, the SEClamp currents are not a true electrode current but represent virtual synaptic input from other regions, so we do need total current to sum to zero at each time step to obtain accurate LFP values)

However, implementing something in NEURON to obtain the SEClamp currents at the same time point as the membrane currrent would be infeasible, because of the following:

I regret that the currents calculated in a BREAKPOINT block are left in an intermediate state on return from fadvance and it is difficult to adjust them to the fadvance exit value without calling the relatively expensive h.fcurrent() function (which would have to be called prior to filling in the vector.record elements).

This leaves us with a few different options for calculating LFP, which we can discuss internally.

[nrnhines]

to clarify, in our use case, the SEClamp currents are not a true electrode current but represent virtual synaptic input from other regions, so we do need total current to sum to zero at each time step to obtain accurate LFP values

I don't think so. My present opinion is that i_membrane_ + SEClamp ==? 0 is just a consistency check. i_membrane_ is correct by itself as the current source for LFP.

Consider the following which I think is counterfactual with regard to "virtual synaptic input from other regions". Ignore if it just seems to muddy the waters. Imagine copying the svclmp.mod file to a virtsyn.mod file and replace

        POINT_PROCESS SEClamp
        ELECTRODE_CURRENT i
...
                i = (vc - v)/rs

with

        POINT_PROCESS VirtSynCur
        NONSPECIFIC_CURRENT i
...
                i = (v - vc)/rs

Several things to note about this. VirtSynCur.i is positive when it is an outward current (positive charge leaving the cell interior across the membrane) In VirtSynCur, v is the membrane potential (not the internal cell potential) VirtSynCur.i is already component of i_membrane_

Contrast this to SEClamp.i is positive when it flows into the interior of the cell. In SEClamp, v is the internal potential of the cell (in the interpreter this is v + vext[0]) SEClamp.i is not included in i_membrane_

For VirtSynCur, please consider carefully the proper sign of vc with respect to how it should affect the LFP value. But note that in a single compartment model, i_membrane_ = 0 no matter what vc is.

[joseph-tharayil]

I think I'm misunderstanding something about this example. The way I see it, with a VirtSynCur, i_membrane_ would indeed be the current source for the LFP, but that is because VirtSynCur.i is already included in i_membrane. Since SEClamp.i is not included in i_membrane_, we need to account for it separately

[nrnhines]

In my opinion, i_membrane_ is the LFP source current whether or not SEClamp.i is present in the model. I.e. i_membrane_ includes the contribution of SEClamp.i to the LFP source current. Are you thinking that there is another source of current for the LFP than i_membrane_.

I am happy to discuss this via zoom. I'm available now and the next few mornings from 7am (1pm Geneva time) onward.

[joseph-tharayil]

Yes, I think it would be a good idea to discuss this over Zoom. I am also now available now and at any point after 1pm Geneva time today or tomorrow. My email is [email protected]; you can send me a Zoom invite for whenever works best for you

[nrnhines]

https://yale.zoom.us/j/99660131601

---

From: joseph-tharayil @.> Sent: Monday, February 26, 2024 8:52 AM To: neuronsimulator/nrn @.> Cc: Hines, Michael @.>; Mention @.> Subject: Re: [neuronsimulator/nrn] Transmembrane currents do not sum to zero (Issue #2686)

Yes, I think it would be a good idea to discuss this over Zoom. I am also now available now and at any point after 1pm Geneva time today or tomorrow. My email is @.@.>; you can send me a Zoom invite for whenever works best for you

— Reply to this email directly, view it on GitHubhttps://github.com/neuronsimulator/nrn/issues/2686#issuecomment-1964197050, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABBRFUMULYHMT7BOELP7RHLYVSHQBAVCNFSM6AAAAABCJGVKW2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNRUGE4TOMBVGA. You are receiving this because you were mentioned.Message ID: @.***>

[joseph-tharayil]

We decided to implement a new mechanism that replicates SEClamp, but as a NONSPECIFIC_CURRENT rather than an ELECTRODE_CURRENT. Since the injected current is then included in i_membrane, that fixes the issue with currents summing to zero