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

[MRG] Implement beta event modulated ERP model #348

Merged
merged 44 commits into from Jun 24, 2021

Conversation

ntolley
Copy link
Contributor

@ntolley ntolley commented Jun 6, 2021

This PR addresses one of the bullet points of #338. The goal is to create a function that instantiates the network from Law et al. 2021 (which was recently accepted!).

The ultimate goal is to load the network via

net = beta_erp_network() # or beta_erp_model() ?

Setting add_drives_from_params=True will simulate the main result of the paper, a prestimulus beta event inhibiting an ERP.

@codecov-commenter
Copy link

codecov-commenter commented Jun 6, 2021

Codecov Report

Merging #348 (708a67d) into master (3cffe1d) will increase coverage by 0.01%.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #348      +/-   ##
==========================================
+ Coverage   89.91%   89.93%   +0.01%     
==========================================
  Files          15       16       +1     
  Lines        2895     2930      +35     
==========================================
+ Hits         2603     2635      +32     
- Misses        292      295       +3     
Impacted Files Coverage Δ
hnn_core/__init__.py 100.00% <100.00%> (ø)
hnn_core/network.py 91.91% <100.00%> (-1.47%) ⬇️
hnn_core/network_models.py 100.00% <100.00%> (ø)
hnn_core/parallel_backends.py 82.85% <0.00%> (-0.86%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3cffe1d...708a67d. Read the comment docs.

@ntolley
Copy link
Contributor Author

ntolley commented Jun 16, 2021

The tyranny of the incredibly awkward all_to_all_connect function is finally over. Now the network and drives are all defined through net.add_connection.

Now to the original purpose of this PR. Really all that's left for this PR is adding ERP's to the example to show how beta events impact their waveform.

hnn_core/network.py Outdated Show resolved Hide resolved
@ntolley ntolley changed the title WIP: Implement beta event modulated ERP model [MRG] Implement beta event modulated ERP model Jun 19, 2021
@ntolley
Copy link
Contributor Author

ntolley commented Jun 19, 2021

@jasmainak once I can check the example on circleCI, I'm going to push the updated tests. Besides that I am pretty happy with this PR and am fairly confident that the results from the paper are being reproduced. We can iterate on the exact wording in some followup PR's, as I imagine we'll be cleaning everything up before the 0.2 release.

@ntolley
Copy link
Contributor Author

ntolley commented Jun 19, 2021

I think a very necessary discussion will be which examples we are ok with using a longer dt. Especially since this PR is adding two more simulations to the pile. At the moment it is set to 0.025 ms or 40kHz which is right in line with state-of-the-art recording setups. While a larger dt may result in less accurate simulations, most (if any) of the examples will see no appreciable change. Especially for subthreshold simulations where we're only interested in 1-100 Hz oscillations (slight inaccuracies in the ODE solver may be more impactful for spiking dynamics).

@jasmainak jasmainak requested a review from cjayb June 19, 2021 22:39
# but modified to reflect that parameters used in Law et al. 2021.
# Specifically, we are considering the case where a tactile stimulus is
# delivered at 150 ms. 25 ms later, the first input to sensory cortex arrives
# as a proximal drive, followed by one distal and a final late proximal drive.
Copy link
Collaborator

Choose a reason for hiding this comment

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

I am reading the paper as I review the code. A couple of points:

  • it would be helpful to list the drives in temporal order
  • also perhaps explain what each drive represents. E.g., first proximal drive represents the thalamic input, then before adding the next drive add a text block explaining that the second drive represents feedback from SII region, and finally the third proximal is due to thalamocortical loop.
  • I know we haven't done this in the past, but we should start to think of having some examples that are "tutorial style" and others which are more "example style". I am guessing that the ones reproducing a paper will be more tutorial style
  • it might also be helpful to mention why the parameters are different. It's because they have been tuned to match the two conditions in the data

Copy link
Contributor Author

Choose a reason for hiding this comment

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

also perhaps explain what each drive represents. E.g., first proximal drive represents the thalamic input, then before adding the next drive add a text block explaining that the second drive represents feedback from SII region, and finally the third proximal is due to thalamocortical loop.

Along with our discussion to flush out the examples, I think this would belong better in the plot_simulate_evoked.py example when we add more thorough dialogue. I can make a brief reference to it here, but explaining the origin of the ERP is not the main purpose of this paper/tutorial.

doc/api.rst Outdated Show resolved Hide resolved
doc/whats_new.rst Outdated Show resolved Hide resolved
return net


def law_model(params=None, add_drives_from_params=False):
Copy link
Collaborator

Choose a reason for hiding this comment

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

is this function covered by any tests?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Also: law_2021_model to conform with convention

Copy link
Contributor Author

@ntolley ntolley Jun 23, 2021

Choose a reason for hiding this comment

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

Added some tests but I'm not sure how extensive they should be. For example, is it necessary to have a separate dpl.txt for each model we create? Or just a few checks like I have added to make sure things aren't breaking.

@ntolley
Copy link
Contributor Author

ntolley commented Jun 24, 2021

@cjayb @jasmainak I'm ready for the next round of reviews on this PR. I have added more details and references to the examples.

If there are glaring typos and errors then let's definitely fix them here. But for more general improvements to the commentary in the example I would prefer iterate on that in a separate PR. This is mainly so I can focus my efforts on the calcium model for the next few weeks.

Copy link
Collaborator

@jasmainak jasmainak left a comment

Choose a reason for hiding this comment

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

Mostly some minor comments. I'll let @cjayb merge if he's happy

constrained_layout=True)
net_beta.cell_response.plot_spikes_hist(ax=axes[0], show=False)
axes[0].set_title('Beta Event Generation')
plot_dipole(dpls_beta, ax=axes[1], layer='agg', tmin=1.0, show=False)
Copy link
Collaborator

Choose a reason for hiding this comment

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

why not smooth here as well? For ease of comparison with the other graphs and real data?

Copy link
Contributor Author

@ntolley ntolley Jun 24, 2021

Choose a reason for hiding this comment

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

Eh I think smoothing here really obstructs your ability to see the unique features of the beta event (sharp distal trough, long asymmetric tail). Also this better reflects the paper where no smoothing was used to plot the beta events (and a mixture of smoothed/unsmoothed for the ERP's)

syn_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
'L5_pyramidal': 0.1}
net.add_evoked_drive(
'evdist1', mu=70.0 + stimulus_start, sigma=0.0, numspikes=1,
Copy link
Collaborator

Choose a reason for hiding this comment

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

this is just for my knowledge. Why is numspikes=1 ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This is how it was simulated in the current study and the original 2009 paper https://journals.physiology.org/doi/full/10.1152/jn.00535.2009

EVOKED RESPONSE INPUTS.
These were simulated as in Jones et al. (2007), with a sequence of FF input, followed by FB, followed by a reemergent late feedforward (LFF) input (schematically drawn in Fig. 9A). The timing of the input sequence was fixed as in Table 2. Each driving spike train consisted of a single presynaptic spike on each trial and the weights were distributed uniformly across the SI networks as in Table 2.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

As for theoretically why it's a single spike, it's a good question. Presumably there is some stochasticity in arrival time at the cortex, and plot_simulate_evoked.py is likely the more realistic simulation of an evoked input.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I see ... I have to say I am super confused what's the interaction between sigma=0, sync_within_trials=True, numspikes=1 and how it affects the histogram plot. I'll raise a separate issue though

# :ref:`evoked example <sphx_glr_auto_examples_plot_simulate_evoked.py>`,
# but modified to reflect the parameters used in Law et al. 2021.
# Specifically, we are considering the case where a tactile stimulus is
# delivered at 150 ms. 25 ms later, the first input to sensory cortex arrives
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think it would be nice to change the time axis so that 150 ms = 0 when plotting the dipole. ERP researchers normally name these waveforms ... N20, P50 etc. based on the latencies. I was trying to figure out which one you're trying to simulate but I'm not sure ... perhaps @cjayb has more opinions

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I can do this in a hacky matplotlib hacky way for this example, but maybe this is an API issue? Like in #315 setting a tstart that is negative should allow for this sort of plotting.

Copy link
Collaborator

Choose a reason for hiding this comment

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

okay yes, let's avoid hacky solutions. Could you raise an issue and I'll see if it can be addressed as part of #315 ?

Comment on lines +65 to +72
# Test succeful creation of a new cell type
net_new = net.copy()
n_cells_orig = net_new.n_cells
new_cell_type = {'new_type': net_new.cell_types['L2_basket']}
net_new.cell_types.update(new_cell_type)
net_new._add_cell_type('new_type', pos=[(0, 0, 0)])
net_new._update_cells()
assert net_new.n_cells == n_cells_orig + 1
Copy link
Collaborator

Choose a reason for hiding this comment

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

rebasing issue?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah the test below caused a merge conflict

Copy link
Collaborator

Choose a reason for hiding this comment

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

I think this needs to be fixed ... this test has moved elsewhere now

def test_add_cell_type():

target_type = conn['target_type']
if src_type not in net.cell_types.keys():
assert set(conn['gid_pairs'].keys()) == set(
net.external_drives[src_type]['conn'][target_type]['src_gids'])
Copy link
Collaborator

@jasmainak jasmainak Jun 24, 2021

Choose a reason for hiding this comment

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

we should unify the two at some point ...

@cjayb cjayb merged commit 572b3fe into jonescompneurolab:master Jun 24, 2021
@cjayb
Copy link
Collaborator

cjayb commented Jun 24, 2021

Fantastic contribution @ntolley ! I'm very impressed by the tutorial, which almost doubled out CircleCI build time ;)

@cjayb
Copy link
Collaborator

cjayb commented Jun 24, 2021

@jasmainak will you pick up on your final review comments as new issues as you see fit?

@jasmainak
Copy link
Collaborator

yeah, sounds good. Final example did turn out quite nice. Great work @ntolley !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants