BlueRecording: A pipeline for the efficient calculation of extracellular recordings in large-scale neural circuit models

Abstract

As the size and complexity of network simulations accessible to computational neuroscience grows, new avenues open for research into extracellularly recorded electric signals. Biophysically detailed simulations permit the identification of the biological origins of the different components of recorded signals, the evaluation of signal sensitivity to different anatomical, physiological, and geometric factors, and selection of recording parameters to maximize the signal information content. Simultaneously, virtual extracellular signals produced by these networks may become important metrics for neuro-simulation validation. To enable efficient calculation of extracellular signals from large neural network simulations, we have developed BlueRecording, a pipeline consisting of standalone Python code, along with extensions to the Neurodamus simulation control application, the CoreNEURON computation engine, and the SONATA data format, to permit online calculation of such signals. In particular, we implement a general form of the reciprocity theorem, which is capable of handling non-dipolar current sources, such as may be found in long axons and recordings close to the current source, as well as complex tissue anatomy, dielectric heterogeneity, and electrode geometries. To our knowledge, this is the first application of this generalized (i.e., non-dipolar) reciprocity-based approach to simulate EEG recordings. We use these tools to calculate extracellular signals from an in silico model of the rat somatosensory cortex and hippocampus and to study signal contribution differences between regions and cell types.

Fig 1. Workflow for the BlueRecording Pipeline.

Input and output files are shown in red boxes, while processes are shown in black boxes. Dashed lines around the FEM simulation step indicate that it is only used for reciprocity-based calculations

Fig 2. Extracellular signals recorded in a large homogeneous medium.

A: Signals recorded with recording electrode (not visible in panel C) distant from the neuron. B: Signals recorded with recording electrode (green dot in panel D) near the neuron. C: Difference in per-compartment weight between generalized reciprocity and dipole-based signal calculations, for electrode far from the neuron (adjusted for constant offset, and normalized to the range of compartment weights in general reciprocity approach). D: The same, for electrode close to the neuron.

Fig 3. BlueRecording simulates resting-state EEG in the rat SSCx model.

A. Rat head model. i: 3D view of the surface-mesh of the rat head model. ii: Comparison between the FEM model brain (green) and the BBP somatosensory cortex model. Locations of EEG electrodes over the forelimb region (dark red) and hindlimb region (blue) and ECoG electrode over the forelimb region (bright red) marked (not to scale). iii. Somatosensory cortex model, with approximate electrode locations as in ii. B: Top down view of mean firing rate over the somatosensory cortex at i. 3200 ms ii. 3250 ms iii. 3300 ms iv. 3350 ms v 3400 ms vi. 3500 ms. C: i: Raster plot of firing over the entire SSCx. Arrows indicate snapshots in panel B. Red arrow indicates start of window highlighted in panels C.iii and C.iv, and Fig 5. ii: EEG and ECoG recorded over the forelimb region. Green and red lines as in C.i. We note that due to baseline current noise injection, the signal is nonzero even in the absence of spiking activity. As the time course of recovery from hyperpolarization at the single-cell is longer than that of the action potential, we observe that the extracellular signal peak is broader than the firing burst. iii. Contributions of different regions of SSCx to the EEG recorded in ii. iv: Contributions of different regions of SSCx to the ECoG recorded in ii.

Fig 4. Differences in compartment weights explain differences between EEG and ECoG.

A–B: Weights for EEG and ECoG recordings, respectively, calculated using the reciprocity approach, for a sample of L5 pyramidal cells in the forelimb (compartments represented as circles), hindlimb (compartments represented as triangles, enclosed in pink box) and upper lip (compartments represented as squares, enclosed in red box) regions. Electrodes are represented as red circles. Note the varying color scale ranges. C–E: Histogram of compartment weights for EEG recording, for L5 pyramidal cells in the forelimb, hindlimb, and upper lip regions, respectively. Dashed lines indicate mean values. F–H: Same as C–E, but for ECoG recording. Note that the x-axis range is different in each figure column.

Fig 5. Simplified calculation methods produce inaccurate signals.

EEG (A), ECoG (B), and LFP (C) signals, recorded over, or within, the somatosensory cortex, calculated with the point-source and line-source approximation, with the generalized reciprocity theorem (ground truth), and with the simplified dipole-based approaches. D-F: Contribution of each region to EEG, ECoG, and LFP signals, respectively. Solid lines indicate general reciprocity approach, dashed lines indicate dipole approach. G-I: Weights for EEG, ECoG, and LFP recordings, respectively, calculated using the reciprocity approach, for a sample of L5 pyramidal cells in the forelimb (compartments represented as circles), hindlimb (compartments represented as triangles, and enclosed in pink box) and upper lip (compartments represented as squares, and enclosed in red box) regions. Electrodes are represented as red circles. Note the varying color scale ranges. J-L: As in G-I, but calculated using the dipole approach

Fig 6. BlueRecording simulates whisker-flick EEG.

A: Selected cells from the 7-column subvolume (blue) and activated thalamic projections (orange). B: Firing rate (first column) and EEG (second column) for the original and the disconnected circuit (red and blue traces, respectively) and the difference in the EEG between the two circuits (third column), for both the full circuit (first row) and each of the layers (subsequent rows). C: Correlation matrices between excitatory (first row) and inhibitory (second row) firing rates in each layer, and the differences in EEG contributions from each layer. In each correlation matrix, firing rates are represented along the rows, and EEG differences along the columns. Correlations are calculated for the full window (first column) and for windows 12-45 ms after the stimulus, and 40-200 ms after the stimulus (second, and third columns, respectively). Start and end times of the windows are marked by red arrows in panel B.

Fig 7. BlueRecording simulates hippocampal LFP.

A: Visualization of the hippocampus and cylindrical subvolumes. B: Visualization of recording electrode placement. C: Raster plot of activity in the full hippocampus (i) and the cylindrical subvolume (ii). D: LFP recorded from a representative electrode in the hippocampus. E: Power-spectral density calculated for the signals in panel D. F: Current source density (CSD) calculated in the full hippocampus simulation. G: Current source density calculated in the cylindrical circuit. CSD maps are calculated using the standard CSD method with Vankin correction [26], as implemented by Rimehaug et al. [27]