Human Neocortical Neurosolver (HNN), a new software tool for interpreting the cellular and network origin of human MEG/EEG data
- Brown University, United States
- Nathan S. Kline Institute for Psychiatric Research, United States
- Yale University, United States
- Massachusetts General Hospital, United States
- Harvard Medical School, United States
- Providence VAMC, United States
Figures
Figure 1
Overview of the biophysical origin of MEG/EEG signals and the relationship between HNN and forward/inverse modeling.
(A) HNN bridges the ‘macroscale’ extracranial EEG/MEG recordings to the underlying cellular- and circuit-level activity by simulating the primary electrical currents (JP) underlying EEG/MEG, which are generated by the postsynaptic, intracellular current flow in the long and spatially-aligned dendrites of a large population of synchronously-activated pyramidal neurons. (B) A zoomed in representation of the relationship between HNN and EEG/MEG forward and inverse modeling. Inverse modeling estimates the location, timecourse and orientation of the primary currents (Jp), and HNN simulates the neural activity creating Jp, at the microscopic scale. Adapted from Jones (2015).
Figure 2
A schematic illustration of a canonical patch of neocortex that is represented by HNN’s underlying neural model.
(Left) 3D visualization of HNN’s model (pyramidal neurons drawn in blue, interneurons drawn in orange. (Right) Commonly measured EEG/MEG signals (ERPs and low frequency rhythms) from a single brain area that can be studied with HNN.
Figure 3
Schematic illustrations of HNN’s underlying neocortical network model.
(A) Local Network Connectivity: GABAergic (GABAA/GABAB; lines) and glutamatergic (AMPA/NMDA; circles) synaptic connectivity between single-compartment inhibitory neurons (orange circles) and multi-compartment layer 2/3 and layer five pyramidal neurons (blue neurons). Excitatory to excitatory connections not shown, see Materials and methods. (B) Exogenous proximal drive representing lemniscal thalamic drive to cortex. User defined trains or bursts of action potentials (see tutorials described in text) are simulated and activate post-synaptic excitatory synapses on the basal and oblique dendrites of layer 2/3 and layer five pyramidal neurons as well as the somata of layer 2/3 and layer five interneurons. These excitatory synaptic inputs drive current flow up the dendrites towards supragranular layers (red arrows). (C) Exogenous distal drive representing cortical-cortical inputs or non-lemniscal thalamic drive that synapses directly into the supragranular layer. User defined trains of action potentials are simulated and activate post-synaptic excitatory synapses on the distal apical dendrites of layer 5 and layer 2/3 pyramidal neurons as well as the somata of layer 2/3 interneurons. These excitatory synaptic inputs push the current flow down towards the infragranular layers (green arrows). (D) The full network contains a scalable number of pyramidal neurons in layer 2/3 and layer 5 in a 3-to-1 ratio with inhibitory interneurons, activated by user defined layer specific proximal and distal drive (see Materials and methods for full details).
Figure 4 with 1 supplement
An example workflow showing how HNN can be used to link the macroscale current dipole signal to the underlying cell and circuit activity.
The example shown is for a perceptual threshold level tactile evoked response (50% detected) from SI (Jones et al., 2007; see ERP Tutorial text for details). (A) Steps 1 and 2: load data and define the local network structure. (B) Step 3: activate the local network, starting with a predefined parameter set; shown here for the parameter set for perceptual threshold-level evoked response (ERPYes100Trials.param) (C) Step 3 and 4: adjust the evoked input parameters according to user defined hypotheses and simulation experiment, and run the simulation. (D) Step 5: visualize model output; the net current dipole will be displayed in the main GUI window and microcircuit details, including layer-specific responses, cell membrane voltages, and spiking profiles (E and F) are shown by choosing them from the View pull down menu. Parameters can be adjusted to hypothesized circuit changes under different experimental conditions (e.g. see Figure 5).
-
Figure 4—source data 1
Average S1 ERP from detected threshold-level stimulus.
- https://cdn.elifesciences.org/articles/51214/elife-51214-fig4-data1-v2.txt
Figure 4—figure supplement 1
Sensitivity analysis results of the perceptual threshold level evoked response example showing the relative contribution of each input’s parameters on variance.
Total Sobol indices at each point have been weighted by the std. deviation scaled from 0 to 1.
Figure 5
Application of HNN to test alternative hypotheses on the circuit level impact of tACS on the somatosensory tactile evoked response (adapted from Sliva et al., 2018).
(A) The early tactile evoked response from above somatosensory cortex before and after 10 min of 10 Hz alternating current stimulation over SI shows that the ~70 ms peak is more prominent in the post-tACS condition. Note that the timing of this peak in the sensor level signal is analogous to the 70 ms peak in the source localized signal in Figure 4B, since the tactile stimulation was the same in both studies and the early signal from SI is similar both at the source and sensor level. (B) HNN was applied to investigate the impact of several possible tACS induced changes in local synaptic efficacy and identify which could account for the observed evoked response data. The parameters in HNN were first adjusted to account for the pre-tACS response using the default HNN parameter set (solid blue line). The synaptic gains between the different cell types was then adjusted through the Set Parameters dialog box to predict that 2x gain in the local inhibitory synaptic weights best accounted for the post-tACS evoked response. (C) Simultaneous viewing of the cell spiking activity further predicted that there is less pyramidal neuron spiking at 70 ms post-tACS, despite the more prominent 70 ms current dipole peak.
Figure 6
Example spontaneous data from a current dipole source in SI showing transient alpha (~7–14 Hz Hz) and beta (~15–29 Hz) components (data as in Jones et al., 2009).
The data file (‘SI_ongoing.txt’) used to generate these outputs is provided with HNN and plotted through the ‘View → View Spectrograms’ menu item, followed by ‘Load Data’, and then selecting the file. (A) The spectrogram viewer with the default ‘jet’ colormap. (B) The configuration option to change the colormap to other perceptually uniform colormaps, where the lightness value increases monotonically. (C) The spectrogram viewer updated with the perceptually uniform ‘viridis’ colormap.
-
Figure 6—source data 1
S1 pre-stimulus activity.
- https://cdn.elifesciences.org/articles/51214/elife-51214-fig6-data1-v2.txt
Figure 7
An example workflow for simulating alpha frequency rhythm (Jones et al., 2009; Ziegler et al., 2010; see Alpha and Beta Rhythms Tutorial text for details).
(A) Here we are using the default HNN network configuration and not directly comparing the waveform to data, so begin with Step 3: activate the local network. Motivated by prior studies (see text), in this example alpha rhythms were simulated by driving the network with ~10 Hz bursts (presumed to be generated by thalamus) to the local network through proximal and distal projection pathways. The parameter set describing these burst is provided in the Alpha.param file and loaded through the Set Parameters From File button. Adjustable burst drive parameter are shown and here were set with a 50 ms delay between the ~10 proximal and distal drive (red boxes). (B) Step 4: running the simulation with the ‘Run Simulation’ button, shows that a continuous alpha rhythm emerged in the current dipole signal (middle dipole time trace; bottom time-frequency representation). Green and red histograms at the top display the defined distal and proximal burst drive patterns, respectively. (C) Step 5: additional network features, including layer specific power spectral density plots as shown can be visualized through the ‘View’ pull down menu, and compared to data (here compared to the spontaneous SI data shown in Figure 6A). Features of the burst drive can be adjusted (panel A) and corresponding changes in the current dipole signals studied (Steps 6 and 7, see Figure 8).
Figure 8
An example workflow for simulating transient alpha and beta frequency rhythm as in the spontaneous SI rhythms shown in Figure 6A (Jones et al., 2009; Ziegler et al., 2010; see Alpha and Beta Rhythms Tutorial text for details).
(A) Here we are using the default HNN network configuration and not directly comparing the waveform to data, so begin with Step 3: activate the local network. In this example, a beta component emerged when the parameters of two ~ 10 Hz bursts to the local network through proximal and distal project pathways, as described in Figure 7, were adjusted so that on average they arrived to the network at the same time (see red boxes). This parameter set is provided in the ‘AlphaAndBeta.params’ file. (B) Step 4: running the simulation with the ‘Run Simulation’ button, shows that intermittent and transient alpha and beta rhythms emerge in the current dipole signal (middle dipole time trace; bottom time-frequency representation). Green and red histograms at the top display the defined distal and proximal burst drive patterns, respectively. Due to the stochastic nature of the bursts, on some cycles of the drive, the distal burst was simultaneous with the proximal burst and strong enough to push current flow down the dendrites to create a beta event (see red box). This model derived prediction reproduced several features of the data, including alpha and beta peaks in the corresponding PSD that were more closely matched to the recorded data (C). Model predictions were subsequently validated with invasive recordings in mice and monkeys (Sherman et al., 2016; see further discussion in text).
Figure 9
An example workflow for simulating canonical pyramidal-interneuron gamma (PING) rhythms (Lee and Jones, 2013; see Gamma Rhythms Tutorial text for details).
(A) Here we are using the default HNN network configuration (with some parameter adjustments as shown in panel B) and we are not directly comparing the waveform to data, so begin with Step 3: activate the local network. Motivated by prior studies on PING mechanisms (see text), in this example PING rhythms were simulated by driving the pyramidal neuron somas with noisy excitatory synaptic input following a Poisson process. The parameters defining this noisy drive are viewed and adjusted through the Set Parameters button as shown, see text and Materials and methods for parameter details. This parameter set for this example is provided in the ‘gamma_L5weak_L2weak.param’ file. In this example, all synaptic connections within the network are turned off (synaptic weight = 0), except for reciprocal connections between the excitatory (AMPA only) and inhibitory (GABAA only) cells within the same layer. The local network connectivity can be viewed and adjusted through the Set Network Connection button or pull down menu, as shown in (B). (C) Step 4: running the simulation with the ‘Run Simulation’ button, shows that a ~ 50 Hz gamma rhythm is produced in the current dipole signal (middle dipole time trace; bottom time-frequency representation). The black histogram at the top displays the noisy excitatory drive to the network. (D–F) Step 5: additional network features, including cell spiking responses, somatic voltages, and layer specific power spectral density plots as shown can be visualized through the ‘View’ pull down menu.
Figure 10
An example simulation describing how parameter adjustments to the canonical PING mechanism can account for source localized visually-induced gamma rhythms recorded with MEG.
(A) An example single trial current dipole waveform generated from a stationary square-wave grating, presented from 0 to 700 ms. MEG data was source localized to the pericalcarine region of the visual cortex (top). The corresponding time-frequency representation (bottom) shows transient bursts of induced gamma activity (bottom). Data is as in Pantazis et al. (2018) and described further in text. (B) Analogous waveforms from 100 trials overlaid (top) and corresponding averaged time-frequency spectrograms (bottom). Averaging in the spectral domain creates the appearance of a more continuous oscillation. (C) Parameters used to simulate the canonical PING rhythm in Figure 9 were adjusted (Steps 6 and 7) to test the hypothesis that a reduced rate of Poisson synaptic input to the pyramidal neurons could account for the bursty nature of the observed gamma data. The parameter set for this example is provided in the ‘gamma_L5weak_L2weak_bursty.param’ file. In this example, the Poisson input rates to the Layer 2/3 and Layer 5 cells were reduced to 4 Hz and 5 Hz, respectively, and the NMDA synaptic inputs were set to a small positive weight (dialog box not shown). Running the simulation will display a histogram of the Poisson drive (top), the current dipole waveform (middle), and the corresponding time-frequency spectrogram (bottom). With this parameter adjustment,~50 Hz bursty gamma activity is reproduced. A scaling factor of 5 was applied to the net current dipole to produce a signal amplitude comparable to the data (<0.1 nAm), suggesting that only ~1000 pyramidal neurons contributed the recorded signal. (D) Poisson input over 100 simulations (top), current dipole waveforms overlaid (middle), and the corresponding spectrogram average (bottom). Similar to the recorded data, averaging in the spectral domain creates the appearance of a more continuous oscillation. (E–F) Microcircuit details including cell spiking responses and layer specific current dipoles shows Layer 5 activity dominates the net current dipole signal, see text for further description.
-
Figure 10—source data 1
Visual cortex gamma activity from 1 trial.
- https://cdn.elifesciences.org/articles/51214/elife-51214-fig10-data1-v2.txt
-
Figure 10—source data 2
Visual cortex gamma activity from 100 trials.
- https://cdn.elifesciences.org/articles/51214/elife-51214-fig10-data2-v2.txt
Figure 11 with 2 supplements
Example of the ERP parameter optimization procedure for a suprathreshold tactile evoked response.
(A) Source localized SI data from a suprathreshold tactile evoked response (100% detection; purple) is shown overlaid with the corresponding HNN evoked response (black) using the threshold level evoked response parameter set detailed in Figure 4, as in initial parameter set. The RMSE between the data and the model is initially high at 30.53. (B) To improve the fit to the data, a procedure for sequentially optimizing the strengths of the proximal and distal drive inputs generating the evoked response can be run. The ‘Configure Optimization’ option is available under the ‘Simulation’ pull down menu. A dialog box allows users to choose and set a range over free optimization parameters, see text for details. (C) The GUI displays an intermediate fit after the first optimization step, specific to the first proximal drive. (D) The final fit is displayed once the optimization is complete. Here, the simulation from the optimized parameter set for the suprathreshold evoked response is shown in gray with an improved RMSE of 15.54 compared to 30.53 for the initial model. See Figure 11—figure supplements 1 and 2 for a further description of the optimization routine and parameter sensitivity analysis.
-
Figure 11—source data 1
Average SI ERP from suprathreshood stimulus.
- https://cdn.elifesciences.org/articles/51214/elife-51214-fig11-data1-v2.txt
Figure 11—figure supplement 1
Weighted scoring of the stepwise optimization procedure for suprathreshold level evoked response example, see text for details.
(A) Probability distribution of input timing for each input. (B) Corresponding cumulative distribution as a contribution index. (C) Unique Contribution Index that is reduced by subsequent input’s CDFs. (D) Extended Contribution Index in which the subtracted CDFs are subject to a decay factor.
Figure 11—figure supplement 2
Sensitivity analysis results of the suprathreshold level evoked response example showing the relative contribution of each input’s parameters on variance.
Total Sobol indices at each point have been weighted by the std. deviation scaled from 0 to 1.
Tables
Key resources table
| Reagent type (species) or resource | Designation | Source or reference | Identifiers | Additional information |
|---|---|---|---|---|
| Software, algorithm | HNN | HNN | RRID:SCR_017437 | https://hnn.brown.edu |
| Software, algorithm | NEURON | NEURON | RRID:SCR_005393 | https://neuron.yale.edu |
Table 1
Length (), diameter () of compartments in each modeled neuron type.
Adend (Bdend) represent apical (basal) dendrite. Note the following connectivity for compartments of PNs. The non-oblique Adends of PNs are connected vertically along the Z axis (cortical layer axis from supra- to infragranular layers) from soma → Adend trunk → Adend1 → Adend2 → Adend tuft. The smaller L2/3 PNs do not have an Adend2 compartment. PN Adend oblique are connected to the soma and perpendicular to the Z axis. Bdend1 connects to the soma along the Z axis, and Bdend2 and Bdend3 branch from Bdend1 at a 45 degree angle from the Z axis. L2/3 and L5 basket interneurons have a single somatic compartment. N/A indicates non-applicable, since that specific compartment not present in the neuron type. Geometry illustrated in Figure 3 above.
| Type | Soma | Adend trunk | Adend1 | Adend2 | Adend tuft | Adend oblique | Bdend1 | Bdend2 | Bdend3 |
|---|---|---|---|---|---|---|---|---|---|
| L2/3 PN | 22.1, 23.4 | 59.5, 4.3 | 306.0, 4.1 | N/A | 238.0, 3.4 | 340.0, 3.9 | 85.0, 4.3 | 255.0, 2.7 | 255.0, 2.7 |
| L2/3 basket | 39.0, 20.0 | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A |
| L5 PN | 39.0, 28.9 | 102.0, 10.2 | 680.0, 7.5 | 680.0, 4.9 | 425.0, 3.4 | 255.0, 5.1 | 85.0, 6.8 | 255.0, 8.5 | 255.0, 8.5 |
| L5 basket | 39.0, 20.0 | N/A | N/A | N/A | N/A | N/A | N/A | N/A | N/A |
Table 2
Space constant (arbitrary units) for synaptic connection strengths and delays between different populations of neurons (rows are pre-synaptic type, columns are post-synaptic type).
| (From ↓, To →) | L2/3 PN | L2/3 Basket | L5 pn | L5 basket |
|---|---|---|---|---|
| L2/3 PN | 3 | 3 | 3 | 3 |
| L2/3 Basket | 50 | 20 | 50 | N/A |
| L5 PN | 3 | N/A | 3 | 3 |
| L5 Basket | N/A | N/A | 70 | 20 |
Additional files
-
Supplementary file 1
Summary of weighted total Sobol sensitivity index values averaged across the entire simulation for each exogenous driving input in two sensory evoked response models: suprathreshold (Figure 11) and 50% detected (Figure 4).
The ranking of each input’s parameters is similar between models.
- https://cdn.elifesciences.org/articles/51214/elife-51214-supp1-v2.xlsx
-
Transparent reporting form
- https://cdn.elifesciences.org/articles/51214/elife-51214-transrepform-v2.docx
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Human Neocortical Neurosolver (HNN), a new software tool for interpreting the cellular and network origin of human MEG/EEG data
eLife 9:e51214.
https://doi.org/10.7554/eLife.51214
