-
Notifications
You must be signed in to change notification settings - Fork 55
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
Conversation
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. |
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) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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 :)
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). |
For validating the "biophysical intuition" a plot like this from may be useful: |
@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? |
I can see about implementing that. Just to make sure I understand, do you mean print statements similar to the |
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. |
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! |
That sounds good! Some of the code here which takes a long time to run can be moved to demo notebooks/gists elsewhere. |
Agree about storing it in |
examples/plot_simulate_lfp.py
Outdated
|
||
h.tstop = 500.0 | ||
lfp_rec = {'lsa': ([], []), 'psa': ([], [])} | ||
for method in ['lsa', 'psa']: |
There was a problem hiding this comment.
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
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() |
There was a problem hiding this comment.
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)
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
How about renaming pos_dict to be |
Love the Also realized an interesting mismatch, the |
def _get_all_sections(sec_type='Pyr'): | ||
ls = h.allsec() | ||
ls = [s for s in ls if sec_type in s.name()] | ||
return ls |
There was a problem hiding this comment.
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?
use_point : bool | ||
Whether to do a point source approximation | ||
for extracellular currents. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Docstring needs updating
# 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) |
There was a problem hiding this comment.
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.
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 |
There was a problem hiding this comment.
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.
# elec.pc.allreduce(elec.lfp_v, 1) | ||
# lfp_py.append(elec_list[pos].lfp_t.to_python()) |
There was a problem hiding this comment.
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.
@ntolley do you have any local changes here? Can I work on this PR if you're busy with other PRs? |
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. |
@jasmainak are you working on this? Looks great! I'd like to propose the following for moving this ahead (and creating a test)
|
I don't think I'll get to this in the coming week. If you or @ntolley want to take over, please feel free! |
Hmm, it's not clear to me how we can proceed to calculate LFP without implementing physical dimensions on |
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 |
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? |
As far as I know x,y position of cells in original model were in 'arbitrary' units primarily used to scale synaptic weights as a function of distance. So original x,y positions should not be considered to be in microns.
…________________________________
From: Nick Tolley ***@***.***>
Sent: Tuesday, April 27, 2021 8:21 PM
To: jonescompneurolab/hnn-core ***@***.***>
Cc: Neymotin, Samuel (NKI) ***@***.***>; Mention ***@***.***>
Subject: Re: [jonescompneurolab/hnn-core] WIP: Add ability to record LFP's (#150)
ATTENTION: This email came from an external source. Do not open attachments or click on links from unknown senders or unexpected emails.
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
@stephanie-r-jones<https://protect2.fireeye.com/v1/url?k=668f8cbf-3914b478-668d758a-000babd9f75c-b2a0885b68cfea2d&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fstephanie-r-jones> @samnemo<https://protect2.fireeye.com/v1/url?k=889af721-d701cfe6-88980e14-000babd9f75c-bf540c1064fc5ea9&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fsamnemo> do you know if there are real units for the x/y position of cells?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<https://protect2.fireeye.com/v1/url?k=e30b4769-bc907fae-e309be5c-000babd9f75c-3131f95b24be9a23&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fjonescompneurolab%2Fhnn-core%2Fpull%2F150%23issuecomment-828047661>, or unsubscribe<https://protect2.fireeye.com/v1/url?k=7e68af16-21f397d1-7e6a5623-000babd9f75c-39f543e44eb1dda9&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FACS5H2ESZ6F4SH2RSOS63S3TK5IJZANCNFSM4QLI44CA>.
________________________________
IMPORTANT NOTICE: This e-mail is meant only for the use of the intended recipient. It may contain confidential information which is legally privileged or otherwise protected by law. If you received this e-mail in error or from someone who was not authorized to send it to you, you are strictly prohibited from reviewing, using, disseminating, distributing or copying the e-mail. PLEASE NOTIFY US IMMEDIATELY OF THE ERROR BY RETURN E-MAIL AND DELETE THIS MESSAGE FROM YOUR SYSTEM. Thank you for your cooperation.
|
Data from that paper used laminar LFP/CSD but afaik not the model. Also, NEURON internally has coordinates for cells but that doesn't mean those positions influence the circuit dynamics.
…________________________________
From: Neymotin, Samuel (NKI) ***@***.***>
Sent: Tuesday, April 27, 2021 9:08 PM
To: jonescompneurolab/hnn-core ***@***.***>; jonescompneurolab/hnn-core ***@***.***>
Cc: Mention ***@***.***>
Subject: Re: [jonescompneurolab/hnn-core] WIP: Add ability to record LFP's (#150)
As far as I know x,y position of cells in original model were in 'arbitrary' units primarily used to scale synaptic weights as a function of distance. So original x,y positions should not be considered to be in microns.
________________________________
From: Nick Tolley ***@***.***>
Sent: Tuesday, April 27, 2021 8:21 PM
To: jonescompneurolab/hnn-core ***@***.***>
Cc: Neymotin, Samuel (NKI) ***@***.***>; Mention ***@***.***>
Subject: Re: [jonescompneurolab/hnn-core] WIP: Add ability to record LFP's (#150)
ATTENTION: This email came from an external source. Do not open attachments or click on links from unknown senders or unexpected emails.
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
@stephanie-r-jones<https://protect2.fireeye.com/v1/url?k=668f8cbf-3914b478-668d758a-000babd9f75c-b2a0885b68cfea2d&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fstephanie-r-jones> @samnemo<https://protect2.fireeye.com/v1/url?k=889af721-d701cfe6-88980e14-000babd9f75c-bf540c1064fc5ea9&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fsamnemo> do you know if there are real units for the x/y position of cells?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<https://protect2.fireeye.com/v1/url?k=e30b4769-bc907fae-e309be5c-000babd9f75c-3131f95b24be9a23&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fjonescompneurolab%2Fhnn-core%2Fpull%2F150%23issuecomment-828047661>, or unsubscribe<https://protect2.fireeye.com/v1/url?k=7e68af16-21f397d1-7e6a5623-000babd9f75c-39f543e44eb1dda9&q=1&e=a8dd2442-b933-4323-949f-ae76ee01f0a2&u=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FACS5H2ESZ6F4SH2RSOS63S3TK5IJZANCNFSM4QLI44CA>.
________________________________
IMPORTANT NOTICE: This e-mail is meant only for the use of the intended recipient. It may contain confidential information which is legally privileged or otherwise protected by law. If you received this e-mail in error or from someone who was not authorized to send it to you, you are strictly prohibited from reviewing, using, disseminating, distributing or copying the e-mail. PLEASE NOTIFY US IMMEDIATELY OF THE ERROR BY RETURN E-MAIL AND DELETE THIS MESSAGE FROM YOUR SYSTEM. Thank you for your cooperation.
|
---------- | ||
electrode_pos : tuple | list of tuple | ||
Coordinates specifying the position for LFP electrodes in | ||
the form of (x, y, z). |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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]) |
There was a problem hiding this comment.
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?
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" ... ? |
1. original x,y positions were used for wiring the neurons in the circuit via synapses (distances were used to set snypatic delays and weights and nothing else in the ORIGINAL model afaik)
2. LFP was never released to public through HNN GUI - needs some work as you're all noticing - part of which involves making sure units are consistent
…________________________________
From: Mainak Jas ***@***.***>
Sent: Wednesday, April 28, 2021 12:39 PM
To: jonescompneurolab/hnn-core ***@***.***>
Cc: Neymotin, Samuel (NKI) ***@***.***>; Mention ***@***.***>
Subject: Re: [jonescompneurolab/hnn-core] WIP: Add ability to record LFP's (#150)
ATTENTION: This email came from an external source. Do not open attachments or click on links from unknown senders or unexpected emails.
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" ... ?
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<https://protect2.fireeye.com/v1/url?k=93df200f-cc441935-93ddd93a-ac1f6b44fec6-5d60d0d64ad4a052&q=1&e=cdd2d0bf-6765-4209-b7a9-850c60417d6f&u=https%3A%2F%2Fgithub.com%2Fjonescompneurolab%2Fhnn-core%2Fpull%2F150%23issuecomment-828604002>, or unsubscribe<https://protect2.fireeye.com/v1/url?k=282ef002-77b5c938-282c0937-ac1f6b44fec6-5e3d108464530bda&q=1&e=cdd2d0bf-6765-4209-b7a9-850c60417d6f&u=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FACS5H2GWMGKZXCQG57DTCGTTLA25JANCNFSM4QLI44CA>.
________________________________
IMPORTANT NOTICE: This e-mail is meant only for the use of the intended recipient. It may contain confidential information which is legally privileged or otherwise protected by law. If you received this e-mail in error or from someone who was not authorized to send it to you, you are strictly prohibited from reviewing, using, disseminating, distributing or copying the e-mail. PLEASE NOTIFY US IMMEDIATELY OF THE ERROR BY RETURN E-MAIL AND DELETE THIS MESSAGE FROM YOUR SYSTEM. Thank you for your cooperation.
|
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 (Murakami & Okada, 2015): doi:10.1016/j.neuroimage.2015.02.003 |
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. |
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. |
closed by #329 |
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.