Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Add ability to record LFP's #150

Closed
wants to merge 40 commits into from

Conversation

jasmainak
Copy link
Collaborator

@jasmainak jasmainak commented Aug 26, 2020

cc @blakecaldwell

is this the right file?

closes #68

I just copied it for now. I'll work on it to improve and understand it.

@jasmainak jasmainak changed the title Start LFP WIP: Start LFP Aug 26, 2020
@blakecaldwell
Copy link
Member

Yes, this is the right file for LFP. Calculating CSD is here too. However, plotting CSD (time vs. depth) is kind of complicated and will probably look very different. Perhaps this PR can just be for LFP.

@jasmainak
Copy link
Collaborator Author

okay, I improved the implementation of PSA and LSA here. Code looks much cleaner if you use numpy.

@blakecaldwell do we have some ground truth? I don't know if we're getting the right result when I run the code.

hnn_core/lfp.py Outdated
"""Enables fast calculation of transmembrane current (nA) at
each segment."""
h.cvode.use_fast_imem(1)
self.bscallback = h.cvode.extra_scatter_gather(0, self.callback)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had to replace the line here with this because my version of Neuron doesn't have the function. Maybe @samnemo do you know what's the right thing to do here?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jasmainak I'm guessing you saw this post?
https://www.neuron.yale.edu/phpBB/viewtopic.php?p=15713#p15713

cvode.extra_scatter_gather seems to replace the cvode.beforestep_py functionality? Anyway, I can't find it in the neuron documentation anymore.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah, I figured from the API documentation that this must be the function :)

@jasmainak
Copy link
Collaborator Author

Here is the plot of the LFP for the two methods

Screen Shot 2020-09-01 at 12 14 41 AM

they match exactly. Does it look correct to people? If so, I'll go ahead and add tests and example

Feel free to pull the branch and play with it.

@blakecaldwell
Copy link
Member

blakecaldwell commented Sep 8, 2020

they match exactly. Does it look correct to people? If so, I'll go ahead and add tests and example

The correctness of these depends on the coordinates of the electrode and cells around it (i.e. the LFP above and below the soma will be drastically different, opposite I think). Probably the best way to verify correctness is to calculate LFP at various points around a SINGLE neuron. Many Einevoll papers did this and it'd be excellent to do the same with our code. This is large enough to be a separate pull issue (perhaps a hackathon issue).

Regarding point source vs. line source: it makes sense that they would be the same with a compartmental model. Line source is said to be more appropriate for compartmental models. However, I don't know which cases the point source approximation would be different. It'd be a worthwhile learning experience to investigate this (for me or anyone).

@ntolley
Copy link
Contributor

ntolley commented Sep 8, 2020

For validating the "biophysical intuition" a plot like this from may be useful:
source
image

@jasmainak
Copy link
Collaborator Author

@ntolley looks like you had a very productive hackathon. I just added the code as a commit here so as to not lose it. Would love to learn from you what you did. Any way we can add some print statements so I know how far the simulation went?

@ntolley
Copy link
Contributor

ntolley commented Sep 15, 2020

I can see about implementing that. Just to make sure I understand, do you mean print statements similar to the simulate() function in parallel_backends.py?

@jasmainak
Copy link
Collaborator Author

yep ... something similar :) You can push straight to this branch. My patience runs out after a script runs for 30 seconds, specially if I see no output.

Also before we merge, do make sure to remind me to credit everyone involved in writing the code.

@jasmainak
Copy link
Collaborator Author

@ntolley let's try to revive this PR once you are done with #276

@ntolley
Copy link
Contributor

ntolley commented Mar 11, 2021

Read my mind, I was actually thinking about this PR today as well. We can discuss later, but I think the most straightforward implementation will be passing a list of tuples containing xyz coordinates. Something like:

electrode_list = [(x1, y1, z1), (x2, y2, z2), ...]
net.record_lfp(electrode_list, method='lsa')

No idea how complicated it may get with parallelization, but we'll see!

@jasmainak
Copy link
Collaborator Author

That sounds good! Some of the code here which takes a long time to run can be moved to demo notebooks/gists elsewhere.

@ntolley ntolley changed the title WIP: Start LFP WIP: Add ability to record LFP's Mar 27, 2021
@ntolley
Copy link
Contributor

ntolley commented Mar 27, 2021

Very rough initial pass at this. Currently electrode positions are specified as a list of tuples under net.lfp_electrode_pos. After running a simulation the output is stored under net.lfp:
image

While the Network class is getting pretty bloated, this seems like a natural place to store electrode information at the moment. Absolutely no idea if this is being calculated correctly, we'll need to think of some good tests to validate the code.

@jasmainak
Copy link
Collaborator Author

Agree about storing it in net. However, you might want to "reset" it before the next simulation. See discussion with @cjayb and @blakecaldwell about simulation times


h.tstop = 500.0
lfp_rec = {'lsa': ([], []), 'psa': ([], [])}
for method in ['lsa', 'psa']:
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this script could go into @cjayb 's repo: https://github.com/cjayb/hnn-core-examples ? Because it runs pretty slow ... could be a nice place to host some of the longer examples

hnn_core/lfp.py Outdated
Comment on lines 210 to 235
from hnn_core.pyramidal import L5Pyr
from hnn_core.network_builder import load_custom_mechanisms

load_custom_mechanisms()
cell = L5Pyr()

h.load_file("stdgui.hoc")
h.cvode_active(1)

ns = h.NetStim()
ns.number = 10
ns.start = 100
ns.interval = 50.0
ns.noise = 0. # deterministic

nc = h.NetCon(ns, cell.synapses['apicaltuft_ampa'])
nc.weight[0] = 0.001

h.tstop = 2000.0

for method in ['lsa', 'psa']:
elec = LFPElectrode([0, 100.0, 100.0], pc=h.ParallelContext(),
method='psa')
elec.setup()
elec.LFPinit()
h.run()
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could be factored into a unit test. If you read the reference paper, I think you might find some corner cases in which the LSA and PSA method give the exact same result ... (I suspect)

Copy link
Contributor

@ntolley ntolley Apr 5, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this the sort of setup you're thinking should go in test_lfp,py for testing the code directly?

I think it will be good to offload the net.add_electrode() tests there as well. However, for a full simulation test I think I will add LFP recordings to test_parallel_ backends.py and compare between joblib and MPI.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah one test like that which checks that LSA and PSA give roughly equal results. I think they would be more equal the further you go away from the cell.

And then one where you could just look at the vres and see if it's computed correctly. I'm worried that lines like this might be incorrect (see also this figure). You can give different positions and see if the vres is computed correctly -- e.g., if it's above a segment vs if it's outside the segment etc.

@jasmainak
Copy link
Collaborator Author

How about renaming pos_dict to be pos_cells (slight misnomer but I can live with it ...) and then calling the LFP position pos_lfp. That way it's more discoverable. In Ipython, when you type net.pos and press tab, you'll see both options.

@ntolley
Copy link
Contributor

ntolley commented Mar 28, 2021

Love the net.pos_lfp suggestion!

Also realized an interesting mismatch, the times recorded for the dipole is one element longer than the LFP recordings (see L55) I believe it has to do with the way the LFP is calculated. Anyone against just copying the first element's value twice?

Comment on lines +39 to +42
def _get_all_sections(sec_type='Pyr'):
ls = h.allsec()
ls = [s for s in ls if sec_type in s.name()]
return ls
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe an opportunity to replace these with better names?

Comment on lines +52 to +54
use_point : bool
Whether to do a point source approximation
for extracellular currents.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Docstring needs updating

Comment on lines +182 to +187
# verify sum i_membrane_ == stimulus
# s = pc.allreduce(imem_vec.sum(), 1)
# if rank == 0: print pc.t(0), s

# sum up the weighted i_membrane_. Result in vx
# rx.mulv(imem_vec, vx)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The allreduce method is typically more efficient than gathering and then summing. Maybe the commented code won't work anymore (b/c switch to imem_ptrvec), but something to consider.

Comment on lines +1062 to +1068
assert sigma > 0.0
_validate_type(method, str, 'method')
_check_option('method', method, ['psa', 'lsa'])
if isinstance(electrode_pos, tuple):
electrode_pos = [electrode_pos]
for e_pos in electrode_pos:
assert len(e_pos) == 3
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

More assert statements that should be converted to Exceptions.

Comment on lines +573 to +574
# elec.pc.allreduce(elec.lfp_v, 1)
# lfp_py.append(elec_list[pos].lfp_t.to_python())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe you meant to remove these? Also, realize that you may have been testing out allreduce where my comment above about efficiency. Please ignore if that's the case. Main point is just that commented lines could be removed.

@jasmainak
Copy link
Collaborator Author

@ntolley do you have any local changes here? Can I work on this PR if you're busy with other PRs?

@ntolley
Copy link
Contributor

ntolley commented Apr 24, 2021

No local changes! Mostly paused this to focus on the connectivity PR's. I think the last item on the list was adding tests to ensure the LFP calculation was correct? My initial thought would be to instantiate a point neuron and/or cylinder with nrn.

@cjayb
Copy link
Collaborator

cjayb commented Apr 26, 2021

@jasmainak are you working on this? Looks great! I'd like to propose the following for moving this ahead (and creating a test)

  • use new connectivity API to prune all connections (non-connected cells, no spiking, no IPSPs)
  • send (weak) synchronous EPSPs to L5 distal apical dendrites (for test: 3 x 3 grid or even a single cell would suffice)
  • plot LFP as a function of depth near the cell (should go from negative to positive at soma)
  • calculate far-field LFP (~10x distal-to-apical distance), which should reproduce a dipole field (see attached, calculated analytically, could be used in a test)

Screenshot 2021-04-26 at 17 29 56

@jasmainak
Copy link
Collaborator Author

I don't think I'll get to this in the coming week. If you or @ntolley want to take over, please feel free!

@cjayb
Copy link
Collaborator

cjayb commented Apr 26, 2021

Hmm, it's not clear to me how we can proceed to calculate LFP without implementing physical dimensions on x and y, in addition to the current z? Relative distances appear in the equations for "transfer resistance". Has anyone considered this? I believe NetPyNE moves the HNN network cells around to enable LFP calculations. What's the principled way of doing it, what's the "lateral extent" of the HNN network?

@cjayb
Copy link
Collaborator

cjayb commented Apr 27, 2021

Are our in-plane coordinates actually in physical (micron) units?! I just assumed they weren't I guess... There's something weird going on with the psa distance calculation, but I can't quite put my finger on it... I think I need to sit down and follow the construction of the neuron sections.

@ntolley
Copy link
Contributor

ntolley commented Apr 28, 2021

Are our in-plane coordinates actually in physical (micron) units?! I

I see why in plane units may be arbitrary for calculating the current dipole, but HNN modeling has previously published with LFP/CSD analysis: https://www.pnas.org/content/113/33/E4885 (figure 9)

@stephanie-r-jones @samnemo do you know if there are real units for the x/y position of cells?

@samnemo
Copy link
Member

samnemo commented Apr 28, 2021 via email

@samnemo
Copy link
Member

samnemo commented Apr 28, 2021 via email

----------
electrode_pos : tuple | list of tuple
Coordinates specifying the position for LFP electrodes in
the form of (x, y, z).
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think @cjayb is right that we should specify the units here. Should be micrometer?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also (x, y, z) is ambiguous: NEURON objects are (x, z, y) in my mind (apical dendrite oriented along z)

plt.figure()
trial_idx = 0
plt.plot(times, net.lfp[0]['data'][trial_idx])
plt.plot(times, net.lfp[1]['data'][trial_idx])
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add units in label ... millivolt?

@jasmainak
Copy link
Collaborator Author

So original x,y positions should not be considered to be in microns.

humm ... okay then I'm confused too. If the x, y positions are not in microns, doesn't that mean that the LFP calculation will be in the wrong units/scale? Maybe the influence of other cells is negligible because they are "far enough" ... ?

@samnemo
Copy link
Member

samnemo commented Apr 28, 2021 via email

@cjayb
Copy link
Collaborator

cjayb commented Apr 28, 2021

Thanks @samnemo . It's more than "just" ensuring consistency. We have to define the physical distance between cells within the "HNN column". Given that our "cells" actually represent quite large assemblies, we should probably not use biological cell-to-cell distances (on the order of microns), but at least an order of magnitude or two larger. We could base it on (Murakami & Okada, 2015)'s 1-2 nAm/mm2 for maximal neocortical tissue dipole moment. I haven't done the math yet to figure out what that implies for us, assuming the usual 10 nAm dipole moment said to be detectable using MEG.

I really like the inclusion of LFP-like signals into hnn-core, but these should not be confused with those obtainable in more realistic models. I'm struggling to get the psa-method to produce the expected dipolar field patterns, though. I need to backtrack the imem-currents coming from NEURON. It's probably "just" a matter of getting the distance calculation right, since the target equation is super-simple: , where n are the sections. But we will need some physical units for all our coordinates, and to keep track of the xyz-coordinates in various places.

(Murakami & Okada, 2015): doi:10.1016/j.neuroimage.2015.02.003

@jasmainak
Copy link
Collaborator Author

jasmainak commented Apr 29, 2021

I think we should just use physical distances based on known numbers from microscopic measurements. I'm not sure we can simply scale the LFP signal if the cell density changes ...

Note that if you get something less than 1 nAm / mm^2, it could simply mean that the cells are not firing synchronously. So, I'm not sure you can calculate backwards the cell density from the signal amplitude. The Okada constant is a maximum. Rather, given a certain signal magnitude from our simulated signals, you can say the signal must be due to x neurons (50,000 for 10 nAm?) firing synchronously. This means that a patch of the cortex that is at least (x / cell density) mm^2 is activated. Interestingly, you can come to the same conclusions from the "scaling factor" in HNN.

@cjayb
Copy link
Collaborator

cjayb commented May 9, 2021

Quick shoutout here: I think I've found a bug in the current implementation. I need to investigate further, but will return with more soonest. I've already taken this branch and rebased against master, hope to have an analysis (and PR) ready by Tuesday.

@cjayb cjayb mentioned this pull request May 10, 2021
12 tasks
@jasmainak
Copy link
Collaborator Author

closed by #329

@ntolley ntolley closed this Jun 17, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
Development

Successfully merging this pull request may close these issues.

add CSD plot for comparing against animal studies
6 participants