Note
Go to the end to download the full example code. or to run this example in your browser via Binder
02. Record extracellular potentials¶
The main output of HNN simulations is the ‘dipole’ waveform, i.e., the net intracellular current flowing in pyramidal cell apical dendrites. At the large distances between cells and M/EEG sensors, this ‘primary’ current is the main contributor to the measured fields. Close to the cells, the local field potential (LFP) is the result of intracellular current leaking into the extracellular medium through active and passive membrane channels. Under some simplifying assumptions, we may approximate the LFP at virtual electrodes placed in and around the HNN network model.
# Authors: Christopher Bailey <cjb@cfin.au.dk>
# Mainak Jas <mainakjas@gmail.com>
# Nick Tolley <nicholas_tolley@brown.edu>
# sphinx_gallery_thumbnail_number = 3
import matplotlib.pyplot as plt
from hnn_core import jones_2009_model, simulate_dipole
from hnn_core.network_models import add_erp_drives_to_jones_model
The default network model defined in Jones et al. (2009) [1] consists of a square grid of pyramidal cells. The in-plane distance between pyramidal cell somas on the grid can be set by the user, which will have an influence on the extracellular potentials (but not on the calculated net intracellular dipole moment). In this example, we’ll simulate a network of model cells spaced 30 um apart. To drive the network dynamics, we’ll use three evoked ‘ERP’ drives; see the event-related potential (ERP) example for details.
net = jones_2009_model()
add_erp_drives_to_jones_model(net)
net.set_cell_positions(inplane_distance=30.)
Extracellular recordings require specifying the electrode positions. It can be useful to visualize the cells of the network to decide on the placement of each electrode.
![plot record extracellular potentials](../../_images/sphx_glr_plot_record_extracellular_potentials_001.png)
<Figure size 640x480 with 1 Axes>
The default network consists of 2 layers (L2 and L5), within which the cell
somas are arranged in a regular grid, and apical dendrites are aligned along
the z-axis. We can simulate a linear multielectrode array with 100 um
intercontact spacing [2] by specifying a list of (x, y, z) coordinate
triplets. The L5 pyramidal cell somas are at z=0 um, with apical dendrites
extending up to z~2000 um. L2 pyramidal cell somas reside at
z~1300 um, and have apical dendrites extending to z~2300 um. We’ll place the
recording array in the center of the network. By default, a value of
0.3 S/m is used for the constant extracellular conductivity and the
‘point source approximation’ for calculations; see
hnn_core.Network.add_electrode_array()
for details.
depths = list(range(-325, 2150, 100))
electrode_pos = [(135, 135, dep) for dep in depths]
net.add_electrode_array('shank1', electrode_pos)
The electrode arrays are stored under Network.rec_arrays
as a dictionary
of hnn_core.extracellular.ElectrodeArray
objects that are now
attached to the network and will be recorded during the simulation. Note that
calculating the extracellular potentials requires additional computational
resources and will thus slightly slow down the simulation.
Using MPI will
speed up computation considerably.
print(net.rec_arrays)
net.plot_cells()
dpl = simulate_dipole(net, tstop=170)
![plot record extracellular potentials](../../_images/sphx_glr_plot_record_extracellular_potentials_002.png)
{'shank1': <ExtracellularArray | 25 electrodes, conductivity=0.3, method=psa (no data recorded yet)>}
Joblib will run 1 trial(s) in parallel by distributing trials over 1 jobs.
Building the NEURON model
[Done]
Trial 1: 0.03 ms...
Trial 1: 10.0 ms...
Trial 1: 20.0 ms...
Trial 1: 30.0 ms...
Trial 1: 40.0 ms...
Trial 1: 50.0 ms...
Trial 1: 60.0 ms...
Trial 1: 70.0 ms...
Trial 1: 80.0 ms...
Trial 1: 90.0 ms...
Trial 1: 100.0 ms...
Trial 1: 110.0 ms...
Trial 1: 120.0 ms...
Trial 1: 130.0 ms...
Trial 1: 140.0 ms...
Trial 1: 150.0 ms...
Trial 1: 160.0 ms...
For plotting both aggregate dipole moment and LFP traces, we’ll use a 10 ms smoothing window, after which both data can be decimated by a factor of 20 from 40 to 2 kHz sampling rates (note that decimation is applied in two steps). Decimation speeds up plotting significantly.
trial_idx = 0
window_len = 10 # ms
decimate = [5, 4] # from 40k to 8k to 2k
fig, axs = plt.subplots(4, 1, sharex=True, figsize=(6, 8),
gridspec_kw={'height_ratios': [1, 3, 3, 3]})
# Then plot the aggregate dipole time series on its own axis
dpl[trial_idx].smooth(window_len=window_len)
dpl[trial_idx].plot(ax=axs[0], decim=decimate, show=False)
# use the same smoothing window on the LFP traces to allow comparison to dipole
net.rec_arrays['shank1'][trial_idx].smooth(window_len=window_len).plot_lfp(
ax=axs[1], decim=decimate, show=False)
axs[1].grid(True, which='major', axis='x')
axs[1].set_xlabel('')
# Add spike raster to subplot
net.cell_response.plot_spikes_raster(ax=axs[2], show=False)
# Finally, add the CSD to the bottom subplot
net.rec_arrays['shank1'][trial_idx].smooth(window_len=window_len).plot_csd(ax=axs[3], show=False)
plt.tight_layout()
plt.show()
![Aggregate (L2 + L5)](../../_images/sphx_glr_plot_record_extracellular_potentials_003.png)
/home/circleci/project/examples/howto/plot_record_extracellular_potentials.py:103: UserWarning: The figure layout has changed to tight
plt.tight_layout()
References¶
Total running time of the script: (1 minutes 6.091 seconds)