Membrane currents do not sum to zero

[somogyvari]

Hi!

We simulated the Hay model, with synaptic stimulation and recorded imem, Ipas, icap and vmem from LFPy. According to my understanding, imem contains the sum of all the currents which are necessary for the LFP calculation. According to the theory, imem should sum up to zero along the whole cell at every time instance. But it does not. If I assume that imem contains only the transmembrane currents, then icap should be added to get the total membrane currents, but it is still very far from zero sum. Adding ipas does not help either. I can imagine, that the synaptic currents are still missing, but the time series of the remainder current sum does not resembles to the time course of the synaptic currents.
Could you help me with this issue? Do I misunderstand something? Which variable contains the total membrane current? Why non of them sums up to zero? Could i record and save the synaptic currents as well?

Thank you for your answers, in advance!

Zoltán Somogyvári

[github-actions[bot]]

Hello @somogyvari , thank you for submitting an issue!

[torbjone]

Hi, and thank you for your interest in using LFPy! Your original understanding is correct: cell.imem should already contain all currents needed to calculate extracellular potentials, including capacitive and synaptic currents, and unless there is also a current injection (?), then it should indeed sum to zero at all time steps.

I'm not sure what the issue is, but when I run the code below (similar to https://github.com/LFPy/LFPy/blob/master/examples/LFPy-example-04.ipynb , assuming mod files are compiled and morphology is available):

=== `import numpy as np import LFPy import neuron

cellParameters = { 'morphology' : 'L5bPCmodelsEH/morphologies/cell1.asc', 'templatefile' : ['L5bPCmodelsEH/models/L5PCbiophys3.hoc', 'L5bPCmodelsEH/models/L5PCtemplate.hoc'], 'templatename' : 'L5PCtemplate', 'templateargs' : 'L5bPCmodelsEH/morphologies/cell1.asc', 'passive' : False, 'nsegs_method' : None, 'dt' : 2**-6, 'tstart' : -159, 'tstop' : 10, 'v_init' : -60, 'celsius': 34, 'pt3d' : True, }

neuron.load_mechanisms('L5bPCmodelsEH/mod/')

cell = LFPy.TemplateCell(**cellParameters)

for sec in cell.allseclist: for seg in sec: seg.e_pas = -59.5

cell.simulate( rec_imem=True) print(f'Maximum net membrane currents at any time step: {np.max(np.abs(np.sum(cell.imem, axis=0)))}')`

=== Then I get: "Maximum net membrane currents at any time step: 4.554689958524705e-14" which seems sufficiently close to the expected numerical accuracy to conclude current conservation is enforced. If you can give some more details on this, I'd be happy to help!

Best, Torbjørn

[somogyvari]

Dear Torbjørn,

Thank you for your prompt answer. We should check what we did differently.

Best regards,

Zoltán

[somogyvari]

Dear Torbjørn,

The main difference between your example and our case is that you excited your cell by setting the membrane potential, while we applied synapses (actually 10 alpha synapses) to excite the cell. As a result, the sum of imem (blue line) is different from zero during this initial period. Does imem contain the synaptic currents? (I think it should. Synaptic currents come from the extracellular space, thus they contribute to the EC field, unlike the electrode currents.) Sorry, if I misunderstood something...

Best regards, Zoltán

[espenhgn]

Hi @somogyvari; I might add: You will have to record also all active currents (i_Na, i_K, etc.) present in the Hay model neuron in order to fetch all membrane currents. The LFPy.Cell.simulate method takes another argument rec_variables=['cai', ...] for this purpose. It's not a well-documented feature, so let us know if you encounter problems implementing it. 'cai' would correspond to the recordable variable (RANGE) of a NMODL mechanism (.mod file) as defined by NEURON.

[torbjone]

Hmmm, this is strange, and there might be some bug relating to AlphaSynapses.

When I try to use LFPy.Synapse with AlphaSynapse, like this, synapse_params = { "idx": 0, "syntype": "AlphaSynapse", } syn = LFPy.Synapse(cell, **synapse_params) syn.set_spike_times(np.array([10, 15, 20, 25]))'

I get this error: NEURON: No NET_RECEIVE in target PointProcess: AlphaSynapse[0] near line 0 {gCa_HVAbar_Ca_HVA(1.0000000000) = 0.0000555000} ^ NetCon(..., ...) Traceback (most recent call last): File "/home/tone/work/playground/test_sum_to_zero.py", line 37, in syn = LFPy.Synapse(cell, **synapse_params) File "/home/tone/work/LFPy/LFPy/pointprocess.py", line 142, in init self.hocidx = int(cell.set_synapse(idx=idx, syntype=syntype, File "/home/tone/work/LFPy/LFPy/cell.py", line 584, in set_synapse nc = neuron.h.NetCon(self._hoc_netstimlist[-1], syn) RuntimeError: hocobj_call error

which I have not seen before, so this needs to be looked more into.

However, when I make an AlphaSynapse like this, I still have current conservation: for sec in cell.allseclist: for seg in sec: syn = neuron.h.AlphaSynapse(seg.x) syn.onset = 10 syn.gmax = 0.2 break cell.simulate( rec_imem=True) How do you implement your AlphaSynapses? If you could supply a minimal example that reproduces the error, that would be helpful!

[78furu]

Dear @torbjone and @espenhgn,

I am a colleague of @somogyvari and we are working on the same data. Our code roughly looks like this:

Importing much more than this, but these are the essentials:

import numpy as np
import LFPy
import neuron
h = LFPy.cell.neuron.h

import matplotlib.pyplot as plt # for plotting

import sys, os
if sys.version < '3':
    from urllib2 import urlopen
else:    
    from urllib.request import urlopen

Compiling the necessary files:

if not os.path.isfile('L5bPCmodelsEH/morphologies/cell1.asc'):
    #get the model files:
    u = urlopen('http://senselab.med.yale.edu/ModelDB/eavBinDown.asp?o=139653&a=23&mime=application/zip',
                context=ssl._create_unverified_context())
    localFile = open('L5bPCmodelsEH.zip', 'wb')
    localFile.write(u.read())
    localFile.close()
    #unzip:
    myzip = zipfile.ZipFile('L5bPCmodelsEH.zip', 'r')
    myzip.extractall('.')
    myzip.close()

#compile mod files every time, because of incompatibility with Mainen96 files:
if "win32" in sys.platform:
    pth = "L5bPCmodelsEH/mod/"
    warn("no autompile of NMODL (.mod) files on Windows.\n" 
         + "Run mknrndll from NEURON bash in the folder L5bPCmodelsEH/mod and rerun example script")
    if not pth in neuron.nrn_dll_loaded:
        neuron.h.nrn_load_dll(pth+"nrnmech.dll")
    neuron.nrn_dll_loaded.append(pth)
else:
    os.system('''
              cd L5bPCmodelsEH/mod/
              nrnivmodl
              ''')
    neuron.load_mechanisms('L5bPCmodelsEH/mod/')

Defining the bases of the cell and the recording sites on a grid:

# define cell parameters used as input to cell-class
cellParameters = {
    'morphology'    : 'L5bPCmodelsEH/morphologies/cell1.asc',
    'templatefile'  : ['L5bPCmodelsEH/models/L5PCbiophys3.hoc',
                       'L5bPCmodelsEH/models/L5PCtemplate.hoc'],
    'templatename'  : 'L5PCtemplate',
    'templateargs'  : 'L5bPCmodelsEH/morphologies/cell1.asc',
    'passive' : False,
    'nsegs_method' : None,
    'dt' : 2**-6,
    'tstart' : 0,
    'tstop' : 50,
    'v_init' : -70,
    'celsius': 34,
    'pt3d' : True,
}

#Generate the grid in xz-plane over which we calculate local field potentials
X, Y, Z = np.mgrid[-4:5:1, 1:2, -4:5:1] * 20
# define parameters for extracellular recording electrode, using optional method
electrodeParameters = {
    'sigma' : 0.3,          # extracellular conductivity
    'x' : X.flatten(),      # x,y,z-coordinates of contacts
    'y' : Y.flatten(),
    'z' : Z.flatten(),
    'method' : 'root_as_point',  #sphere source soma segment
    'N' : np.array([[0, 1, 0]]*X.size), #surface normals
    'r' : 2.5,              # contact site radius
    'n' : 20,               # datapoints for averaging
}

Creating the cell:

# delete old sections from NEURON namespace
LFPy.cell.neuron.h("forall delete_section()")

# Initialize cell instance, using the LFPy.Cell class
cell = LFPy.TemplateCell(**cellParameters)
cell.set_rotation(x=4.729, y=-3.166)

# create extracellular electrode object for LFPs on grid
electrode = LFPy.RecExtElectrode(cell, **electrodeParameters)

Generating the Alphasynapses:

all_sections = list(h.allsec())

chosen_section_ids = [97, 98, 105, 121, 135, 137, 176, 190, 191, 192]

alphasynapses = {}

for secid in chosen_section_ids:
    f = h.AlphaSynapse(all_sections[secid](1))
    f.onset = 10
    f.tau = .1
    f.gmax = 1
    f.e = 40
    alphasynapses[secid] = f

Creating recordings and simulating cell

alphsyn_recs = [h.Vector().record(alphasyn._ref_i) for alphasyn in alphasynapses.values()]

t = h.Vector().record(h._ref_t)                     # Time stamp vector
cell.simulate(rec_icap=True, rec_vmem=True, rec_ipas=True, rec_imem=True, probes=[electrode])

Converting to numpy and plotting:

alphsyn_recs_np= np.array([v.to_python() for v in alphsyn_recs])

fig, axes = plt.subplots(2,1, figsize=(15,10), facecolor = 'w')

for ax in axes.ravel():

    ax.plot(t, np.sum(cell.imem, axis=0), label = 'imem')
    ax.plot(t, np.sum(cell.icap, axis=0), label = 'icap')
    ax.plot(t, np.sum(cell.ipas, axis=0), label = 'ipas')

    for c, v in enumerate(alphsyn_recs_np):
        ax.plot(t, v, ".-", label = f"i_alphasyn #{c}")

    #ax.plot(t, np.sum(alphsyn_recs_np, axis=0), ls = '--', label = f"alphasyn #{c}")

    ax.axvline(10, label = 'Stimuli on drendites', color = 'k')

    ax.set_xlabel('Time [ms]')
    ax.set_ylabel('Current [nA]')

    ax.grid()
    ax.legend(loc = "upper right")

axes[1].set_xlim( (9.9, 10.5))

fig.suptitle('Comparison of imem, icap, ipas and alphasynapses', fontsize=20, fontweight = 'bold')
#fig.savefig('cap_imem_alphasyn_comparison.png')

Which should result in the following figure:

As one can see, the imem after stimulation does not sum up to zero. If one would take the sum of the 10 alphasynapes' currents, that would overcompensate.

Note: as opposed to the previous figure, added by @somogyvari I removed the np.abs() when plotting the currents, hence the difference in the orange line around 13 and 25 ms.

My personal question after surfing the documentation is the following: why has the rec_isyn = True parameter been removed in cell.simulate()?

Thank you for your answers and efforts, Kristof

[espenhgn]

@torbjone; I could reproduce your error. Not sure what this may be - certainly not the definition of the AlphaSynapse mechanism, which has not changed in >16 years (https://github.com/neuronsimulator/nrn/commits/05deaaf12022887762606421e612b18045f5270c/src/nrnoc/syn.mod). It should be made into a separate ticket, unrelated to the issue above by @somogyvari.

As for that issue, I messed around a bit in this notebook with summing all currents. At first glance, the sum of synaptic currents is well balanced by the summed capacitive currents (as expected), but things go a bit wrong during the generated AP.

However, based on responses made by Ted Carnevale in this thread there are summation issues that arise due to the numerical integration in NEURON itself.

To get a better hold on this, I'd suggest investigating a 1-compartment model with only a single synapse, and a passive membrane before including active ion channels in the model.

[espenhgn]

Thanks for providing the example @78furu. Have a look at the notebook I linked to above: With this active cell model, ionic currents have to be included in addition to capacitive, leaky, and synaptic currents, but differences will occur due to the numerical integration scheme in NEURON (to my understanding) - so for the purpose of computing extracellular potentials etc., please use the rec_imem=True option :)

I removed the "rec_isyn" argument at one point as it was redundant with the LFPy.PointProcess argument record_current (via https://github.com/LFPy/LFPy/pull/37)

[torbjone]

Hi,

it seems there really is a problem with current conservation here, and I think the problem is at this line: f = h.AlphaSynapse(all_sections[secid](1))

here, the synapse is at location 1 (the endpoint of the segment), and locally, I also do not get current conservation if i use exactly 0 or 1, but any number in between results in current conservation. So the quick fix I guess is to use 0.5, so that the synapse is at the middle of the segment, since this is more common anyway. Is that OK for your purposes? With regards to why it does not work for 0 or 1, I have to look more into this :-)

[torbjone]

Btw Espen: As I understand it, the question is not how to sum up current components in the correct way, but why cell.imem does not sum to zero over the cell.

[torbjone]

This actually seems to be a NEURON issue, and not an LFPy issue, as I can reproduce it in NEURON without LFPy: In the script below, there is current conservation for any value of 0 < x < 1, but not for x=0 or x=1. I still don't know why though.

import numpy as np
import matplotlib.pyplot as plt
import neuron

h = neuron.h # neuron imports hoc and does a h = hoc.HocObject()

cvode = neuron.h.CVode()

cvode.use_fast_imem(1)

soma = h.Section(name='soma')
soma.nseg = 1
soma.insert("pas")

dend = h.Section(name='dend')
dend.nseg = 1
dend.insert("pas")

dend.connect(soma(1), 0)

all_sec = h.SectionList()

for sec in h.allsec():
    all_sec.append(sec = sec)

x = 1
f = h.AlphaSynapse(soma(x))
f.onset = 10
f.tau = .1
f.gmax = 0.5
f.e = 0

t = neuron.h.Vector()
v = neuron.h.Vector()

t.record(neuron.h._ref_t)
v.record(soma(0.5)._ref_v)

imem_ = []
for sec in h.allsec():
    for seg in sec:
        i = neuron.h.Vector()
        i.record(seg._ref_i_membrane_)
        imem_.append(i)

##################################################################
# Simulation control
##################################################################
neuron.h.dt = 0.05          # simulation time resolution
tstop = 15.        # simulation duration
v_init = -70        # membrane voltage(s) at t = 0

neuron.h.finitialize(v_init)
neuron.h.fcurrent()

while neuron.h.t < tstop:
    neuron.h.fadvance()

# Extract membrane currents
imem = []
comp_idx = 0
for sec in h.allsec():
    for seg in sec:
        imem.append(np.array(imem_[comp_idx]))
        comp_idx += 1
imem = np.array(imem)

##################################################################
# Plot simulated output
##################################################################
fig, axes = plt.subplots(2)

for comp_idx in range(imem.shape[0]):
    axes[0].plot(t, imem[comp_idx], 'gray', lw=1.,
                 label=f'comp idx: {comp_idx}')

axes[0].plot(t, np.sum(imem, axis=0), 'r', lw=1.5,
             label="net sum over membrane")

axes[0].set_ylabel('current (nA)')
axes[0].legend()

axes[1].plot(t, v, 'k', lw=1.5)
axes[1].set_ylabel('voltage (mV)')
axes[1].set_xlabel('time (ms)')

plt.show()

[78furu]

Dear @torbjone and @espenhgn,

Thank you for your comments, setting the location from 1 to 0.5 is sufficient for our calculations.

Best regards, Soma & Kristóf

[somogyvari]

Dear @torbjone and @espenhgn,

Many thanks! This quick fix is perfect for us now! The exact location of the synapses do not matter, so we relocate them accordingly. We really appreciate your quick help!

Best reagards! Zoltán

[espenhgn]

This actually seems to be a NEURON issue, and not an LFPy issue, as I can reproduce it in NEURON without LFPy: In the script below, there is current conservation for any value of 0 < x < 1, but not for x=0 or x=1. I still don't know why though.

Good that the issue is resolved. I think this is because the section locations at $x\in\{0, 1\}$, strictly speaking, do not belong to the membrane, but the node or terminus points connecting section elements that make up the morphology. It may be documented in the NEURON book or in its documentation.