How to apply LFPykit on more than one neuron

[llandsmeer]

The README states:

the extracellular potentials can be calculated from a distance-weighted sum of contributions from transmembrane currents of neurons

suggesting that LFPykit Is able to handle more than one neuron (which is logical given that we are looking at LFP's). Now, both documentation and examples show LFPykit applied to only one neuron, eg. the form PointSourcePotential(cell=CellGeometry(...)).

Does this mean that LFPykit can only handle a single neuron? Or is there the possibility to apply it to large networks of multicompartmental neurons (eg by just adding all neurons into a single CellGeometry)? Please let me know, we'd really like to apply this our network simulated using Arbor

Kind regards,
Lennart Landsmeer

[github-actions]

Hello @llandsmeer , thank you for submitting an issue!

[llandsmeer]

So reading through some more code and Arbor/LFPy issues, the general process seems to be

1. calculate LFP's for each neuron using something like PointSourcePotential(cell=CellGeometry(...))
2. 'reduce the individual contributions' - which just means summing together?

Is that correct?

[espenhgn]

Hello @llandsmeer;
At this moment we haven't yet prepared any example network simulations incorporating extracellular signal predictions via LFPykit in Arbor.
Right now there is only this single-cell example provided with Arbor: https://github.com/arbor-sim/arbor/blob/master/python/example/single_cell_extracellular_potentials.py (which is just a slightly reworked version of this notebook https://github.com/LFPy/LFPykit/blob/master/examples/Example_Arbor_swc.ipynb).

So in this context what you outline above would be the way to go:

• set up the network simulation in Arbor
• record the transmembrane currents per control volume (CV) using arbor.cable_probe_total_current_cell()
• get information about the 3D geometry/segments of each CV using the arbor.place_pwlin API as in the above example. Right now this is a bit cumbersome as each CV may have a complex geometry and the different segments need to be mapped to the corresponding CV. What I incorporated in the example is that I first map segments to the different extracellular measurement sites; then weigh the different segments into another matrix allowing CV transmembrane currents to the different extracellular measurement sites (see https://github.com/arbor-sim/arbor/blob/2bf5404776df258373c2df5f67482edf4daf9c82/python/example/single_cell_extracellular_potentials.py#L206). I had a recent discussion with the main Arbor devs to improve this procedure here: LFPykit and Arbor arbor-sim/arbor#1806
• Depending on how the Arbor network simulation is distributed using MPI; if each neuron exist on different MPI processes one could sum up all contributions post simulation using say MPI.Reduce(data, op=MPI.SUM, root=0) (see mpi4py docs). If more than one neuron exist per MPI process, first sum over local neuron contributions, followed by summing using MPI.Reduce.

In terms of what we have done with networks and LFPykit; LFPy itself supports computing extracellular signals from networks directly: https://github.com/LFPy/LFPy/blob/master/examples/example_network/example_network.py. In LFPy the procedure above is modified somewhat as it made sense with NEURON specifically. We build linear maps between compartments and measurements for all compartments across cells per MPI process pre simulation. Then at each time point in the simulation build a vector of compartment transmembrane currents and multiply this with the mapping matrix. Post simulation sum up across MPI processes and export the summed data. This avoids recording transmembrane currents altogether (which can be very memory demanding).

Please let me know if you have further questions. I think having an actual use case with networks in Arbor would be very handy.

[espenhgn]

Hi again @llandsmeer. For arbor-sim/arbor#1825 I've now incorporated some improvements mainly to the notebook https://github.com/LFPy/LFPykit/blob/master/examples/Example_Arbor_swc.ipynb. In there I define a couple new classes ArborCellGeometry and ArborLineSourcePotential to simplify the code required to compute extracellular potentials.

[llandsmeer]

Hi @espenhgn, thanks so much for the detailed explanation and for updating the arbor example notebook with reusable classes!

Some information about the network:

• ~1000 cells
• ~200.000 CV's, mostly linear but can be branched
• Most simulation runs are in the order of 5-20 seconds
• Runs completely on a single GPU in a single cell group, no MPI
• Development GPU has 25 GB RAM, which might allow us to record for some time (6000 samples as a rough estimate)
• There are a lot of exactly overlapping cables from different neurons- I hope that's not a problem

Thanks to your messages and the example I got a minimal version working over here: https://github.com/llandsmeer/iopublic/blob/main/Local%20field%20potential.ipynb

Some notes:

• Puling out the morphologies from a recipe is a bit cumbersome at the moment. With a proper API to that, we could AFAICS have a generic LFPykit solution for arbor recipes (the only other thing being that you need to modify the probes() function)
• I do get nan values out of the M matrix transformation. That probably means I'm doing something wrong :)
• Calling ArborLineSourcePotential(...).get_transformation_matrix() for each neuron is slow. At 64x64 points, this step takes 3x as much time as the simulation. It's probably embarrassingly parallel but I still have to figure out how to parallelize it without GIL problems
• Summing up M matrix transforms is also a bit slow
• During postprocessing I went Out Of Memory - that's 126GB

The results are ... not that interesting yet :)

(this is a full network crosssection)

[espenhgn]

Hi @llandsmeer;
Great that you could incorporate LFPykit, although the output is not as expected. It's also great to have an idea of network size/compartment count to get an idea of the magnitude of this problem.

I did try and run your notebook locally on my M1 Mac but failed as a CUDA library is needed during the recipe = iopublic.build_recipe(**kwargs) step. Hence I can mostly speculate what went wrong above:

• In terms of performance, I think this can be explained by the fact that LFPykit is written mostly in pure Python and numpy, resulting in the usual performance issues with nested for loops. Here we appear to have an outer for loop over cells, and in class LineSourcePotential there is another loop over each extracellular measurement site in the function that returns the M matrix:

  LFPykit/lfpykit/models.py

  Line 437 in 5f8bde5

  def get_transformation_matrix(self):

  
  The function that is being called avoids further for loops, but has not seen much in terms of profiling:

  LFPykit/lfpykit/lfpcalc.py

  Line 346 in 5f8bde5

  def calc_lfp_linesource(cell, x, y, z, sigma, r_limit):

  
  I assume these operations can be made faster using Cython or numba for instance. I have little experience with GPU programming, but I assume one could distribute calculations relevant for each cell across threads on the GPU (I guess you also have multiple CPU cores) with relative ease? The GIL could cause issues I know. The problem is indeed embarrassingly parallel.
• Whether or not segments overlap in space should not matter, as there are no mutual interactions between them. Their contributions to the extracellular potential should anyway sum linearly.
• In terms of NaN values, please check that no segments in your segtree objects have zero lengths or zero diameters. These segments should not have any transmembrane currents that could possibly contribute to the extracellular potential anyway.
• You could perhaps concatenate both the list of matrices M and corresponding I_m along their appropriate axes and compute the extracellular potential as M_all @ I_m_all, although both matrices will be huge depending on the number of extracellular points. As a test, what happens if you reduce the number of extracellular points? I assume you could compute this product also on the GPU (https://scikit-cuda.readthedocs.io/en/latest/generated/skcuda.linalg.dot.html)

I see Brent already gave some feedback to my PR, so I may update what I did here with the Arbor example(s)

[espenhgn]

Btw., I thought this plot sort of made sense if this indeed the extracellular potential of one cell (at index 30):

From what you showed above, perhaps you sum over the wrong axis? To me it looks like you should sum over axis=1, not axis=0

[llandsmeer]

Ah sorry I didn't explain that I think. Index 30 here is a spatial index. So we have time on the X-axis and 1-dimensional space on the Y axis. V_e already contains the summed up potentials over all cells. I'm pretty sure I'm summing over the right axes (else the output shape is not correct)

To what extent does it look wrong? I think I forgot to say that the inferior olive cells show threshold oscillations - about 2/3 of the cells in the network always oscillate. They also synchronize, and as far as I can see the global synchronization mode is displayed in the LFP's. But this is also the first time I'm doing this kind of analysis, so I might be doing something wrong of course (we're working on experimental validation of the LFP's currently).

[llandsmeer]

Writing this I just realized I might have made a very obvious mistake, I'll try to get back on this asap

[espenhgn]

Writing this I just realized I might have made a very obvious mistake, I'll try to get back on this asap

No worries!

[llandsmeer]

Yes I had an error - I was placing all segments at position (0, 0, 0). That didn't matter for the cable cell equations but obviously changes the LFPs..

Anyway I fixed that - now I get this, which looks much more logical
(x, y) crosssection

(x, t) crosssection

[espenhgn]

Very nice! Did you manage to address the performance/memory issues you encountered in a good way (beyond reducing the number of locations in space)? I did open an issue here #153. Seems to me that a few operations could be sped up using numba with relative ease.

[llandsmeer]

I just ran some small tests. The results are pretty obvious for you probably, but it was nice to dive into the codebase a bit more. In general, the slow calls are ArborCellGeometry() and lsp.get_transformation_matrix()

1. Directly using a concurrent.futures.ThreadPoolExecutor(max_workers=128 slows things down, with CPU load staying around 1, suggesting some GIL problems.

2. Replacing repeated calls to numpy.row_stack with appending to a list and np.array(, dtype=float) speeds up ArborCellGeometry.__init__ by about 2.5 times

class ArborCellGeometry(lfpykit.CellGeometry):
    def __init__(self, p, cables):
        x, y, z, r = [], [], [], []
        CV_ind = np.array([], dtype=int)  # tracks which CV owns segment
        for i, m in enumerate(cables):
            segs = p.segments([m])
            for j, seg in enumerate(segs):
                x.append([seg.prox.x, seg.dist.x])
                y.append([seg.prox.y, seg.dist.y])
                z.append([seg.prox.z, seg.dist.z])
                r.append([seg.prox.radius, seg.dist.radius])
                CV_ind = np.r_[CV_ind, i]
        x = np.array(x, dtype=float)
        y = np.array(y, dtype=float)
        z = np.array(z, dtype=float)
        d = 2*np.array(r, dtype=float)
        super().__init__(x=x, y=y, z=z, d=d)
        self._CV_ind = CV_ind

3. Seems like we're limited to calling ArborCellGeometry() single threaded for now as it's all python objects. So then I looked at using a ThreadPoolExecutor to speed up just lsp.get_transformation_matrix ... which only slowed things down by a factor of about 2x.

4. Within ArborLineSourcePotential.get_transformation_matrix, the call to LineSourcePotential.get_transformation_matrix() takes up 90% of the time, the CV reshuffling the rest. Within this call, the loop over segments, calling lfpcalc.calc_lfp_linesource(cell, x[j], y[j], z[j], sigma, r_limit) is the bottleneck

5. Within calc_lfp_linesource, the only input argument not being a float or ndarray is cell, which is immediately destructed into its constituent ndarrays. It calls the functions _deltaS_calc, _h_calc, _r2_calc, _linesource_calc_case{1,2,3}`, which all take in numpy arrays and output numpy arrays, suggesting they're fine with numba.jit(nopython)

So indeed, adding numba.jit(nopython=True, nogil=True) to some functions in lfpcalc.py seems like a good solution. Then we either need to destruct cell before calling the numba function or make it a jitclass. I'll try the first solution now in fork