The concepts of electromagnetic biophysics are succinctly discussed, for example in Hämäläinen et al. (1993) and Hari and Ilmoniemi (1986). Here, we briefly review the basic framework of forward and inverse modeling and how they relate to HNN (Figure 1).

Figure 1

Download asset Open asset

MEG/EEG signals are created by electrical currents in the brain. The signals are recorded by sensors at least a centimeter from the actual current sources. Given this configuration, a macroscopic scale is employed for the distribution of electrical conductivity. The division between the actual non-ohmic equivalent current sources of activity and the passive ohmic currents is then referred to as primary (J^p) and volume currents (J^V), respectively. As depicted in Figure 1, both MEG and EEG are ultimately generated by the primary currents. The primary currents set up a potential distribution (V) that extends through the brain tissue, the cerebrospinal fluid (CSF), the skull, and the scalp, where it is measured as EEG; the passive J^V is proportional to the electric field (negative gradient of V) and electrical conductivity (σ). MEG, in general, is generated by both the J^p and J^v. The total current is the sum of J^p and J^v: J = J^p + J^v, whence J^p = J – J^v = J + σ∇V. In this definition, the conductivity is considered on a macroscopic scale, omitting the cellular level details. The primary current is the ‘battery’ of the circuit and is nonzero at the active sites in the brain. Furthermore, the direction of the primary current is determined by the cellular level geometry of the active cells. By locating the primary currents from MEG/EEG we locate the sites of activity.

The task of computing EEG and MEG given J^p is commonly called forward modeling, and it is governed by Maxwell’s equations. However, in the geometry of the head, the integral effect of the volume currents to the magnetic field can be relatively easily taken into account and, therefore, modeling of MEG is in general more straightforward than the precise calculation of the electric potentials measured in EEG. Specifically, a first order approximation, the spherically symmetric conductor model (Sarvas, 1987), can often be used (Tarkiainen et al., 2003). In this case, all components of the magnetic field can be computed from an analytical formula which is independent of the value of the electrical conductivity as a function of the distance from the center of the sphere (Sarvas, 1987). In a more complex case with realistically shaped conductivity compartments, the skull and the scalp can be replaced by a perfect insulator (Hämäläinen and Sarvas, 1987; Hämäläinen and Sarvas, 1989).

The availability of these forward models opens up the possibility to estimate the locations (r) and time course of the primary current activity, J^p = J^p(r,t) from MEG and EEG sensor data, that is inverse modeling, or source estimation. However, this inverse problem is fundamentally ill-posed, and constraints are needed to render the problem unique. The different source estimation methods, such as current dipole fitting, minimum-norm estimates, sparse source estimation methods, and beamformer approaches, differ in their capability to approximate the extent of the source activity and in their localization accuracy; there are presently several open-source software packages for source estimation, for example Gramfort et al. (2013). All of these methods are capable of inferring both the location and direction of the neural currents and their time courses. Importantly, due to physiological considerations, the appropriate elementary primary current source in all of these methods is estimated as a current dipole, with units current x distance, that is Am (Ampere x meter).

Thanks to the consistent orientation of the apical dendrites of the pyramidal cells in the cortex this primary current is oriented normal to the cortical mantle and its direction corresponds to the intracellular current flow (Okada et al., 1997; Ikeda et al., 2005; Murakami and Okada, 2006). When source estimation is used in combination with geometrical models of the cortex constructed from anatomical MRI, the current direction can be related to the direction of the outer normal of the cortex: one is thus able to tell whether the estimated current is flowing outwards or inwards at a particular cortical site at a particular point in time. As such, the direction of the current flow can be related to orientation of the pyramidal neuron apical dendrites and inferred as currents flow from soma to apical tuft (up the dendrites) or apical tuft to soma (down the dendrites), see Figure 1.

The focus of HNN is to study how J^p is generated by the assembly of neurons in the brain at the microscopic scale. Currently, the process of estimating the primary current sources with inverse methods, or calculating the forward solution from J^p to the measured sensor level signal, is separate from HNN. A future direction is to integrate the top-down source estimation software with our bottom-up HNN model for all-in-one source estimation and circuit interpretation (see Discussion).

HNN’s underlying neural model contains elements that can simulate the primary current dipoles (J^p) creating EEG/MEG signals in a biophysically principled manner (Figure 1). Specifically, HNN simulates the primary current from a canonical model of a layered neocortical column via the net intracellular electrical current flow in the pyramidal neuron dendrites in a direction parallel to the apical dendrites (see red arrow in Figures 1 and 2, and further discussion in Materials and methods) (Hämäläinen et al., 1993; Ikeda et al., 2005; Jones, 2015; Murakami et al., 2003; Murakami and Okada, 2006; Okada et al., 1997). With this construction, the units of measure produced by the model are the same as those estimated from source localization methods, namely, ampere-meters (Am), enabling one-to-one comparison of results. This construction is unique compared to other EEG/MEG modeling software (see Discussion). A necessary step in comparing model results with source-localized signals is an understanding of the direction of the estimated net current in or out of the cortex, which corresponds to current flow down or up the pyramidal neuron dendrites, respectively, as discussed above. Estimation of current flow orientation at any point in time is an option in most inverse solution software that helps guide the neural interpretation, as does prior knowledge of the relay of sensory information in the cortex, see further discussion in the Tutorials part of the Results section.

Figure 2

Download asset Open asset

By keeping model output in close agreement with the data, HNN’s underlying model has led to new and generative predictions on the origin of sensory evoked responses and low-frequency rhythms, and on the changes in these signals across experimental conditions (Jones et al., 2009; Jones et al., 2007; Khan et al., 2015; Lee and Jones, 2013; Sherman et al., 2016; Sliva et al., 2018; Ziegler et al., 2010) described further below. The macro- to micro-scale nature of the HNN software is designed to develop and test hypotheses that can be directly validated with invasive recordings or other imaging modalities (see further discussion in tutorial on alpha and beta rhythms).

HNN is currently constructed to dissect the cell and network contributions to signals from one source-localized region of interest. Specifically, the HNN GUI is designed to simulate sensory evoked responses and low-frequency brain rhythms from a single region, based on the local network dynamics and the layer-specific thalamo-cortical and cortico-cortical inputs that contribute to the local activity. As such, HNN’s underlying neocortical network represents a scalable patch of neocortex containing canonical features of neocortical circuitry (Figure 2). Ongoing expansions will include the ability to import other user-defined cell types and circuit models into HNN, simulate LFP and sensor level signals, as well as the the interactions among multiple neocortical areas (see Discussion). Of note, users can still benefit from our software if they are working with data directly from EEG/MEG sensor rather than source-localized signals. The primary currents are the foundation of the sensor signal and, as such, can have similar activity profiles (e.g., compare source-localized tactile evoked response in Figure 4 and sensor-level response in Figure 5).

Here, we give an overview of the main features that are important to understand in order to begin exploring the origin of macroscale evoked responses and brain rhythms, and we provide details on how these features are implemented in HNN’s template model. Further details can be found in the Materials and methods section, in our prior publications (e.g., Jones et al., 2009), and on our website https://hnn.brown.edu.

Given that the primary electrical current that generates EEG/MEG signals comes from synchronous activity in pyramidal neuron (PN) dendrites across a large population, there are several key features of neocortical circuitry that are essential to consider when simulating these currents. While there are known differences in microscale circuitry across cortical areas and species, many features of neocortical circuits are remarkably similar. We assume these conserved features are minimally sufficient to account for the generation of evoked responses and brain rhythms measured with EEG/MEG, and we have harnessed this generalization into HNN’s foundational model, with success in simulating many of these signals using the same template model (see Introduction). These canonical features include:

(I) A 3-layered structure with pyramidal neurons in the supragranular and infragranular layers whose dendrites span across the layers and are synaptically coupled to inhibitory interneurons in a 3-to-1 ratio of pyramidal to inhibitory cells (Figure 3A). Of note, cells in the granular layer are not explicitly included in the template circuit. This initial design choice was based on the fact that macroscale current dipoles are dominated by PN activity in supragranular and infragranular layers. Thalamic input to granular layers is presumed to propagate directly to basal and oblique dendrites of PN in the supragranular and infragranular layers. In the model, the thalamic input synapses directly onto these dendrites.

Figure 3

Download asset Open asset

The number of cells in the network is adjustable in the Local Network Parameters window via the Cells tab, while maintaining at 3-to-1 pyramidal to inhibitory interneuron ratio in each layer. The connectivity pattern is fixed, but the synaptic weights between cell types can be adjusted in the Local Network menu and the Synaptic Gains menu. Macroscale EEG/MEG signals are generated by the synchronous activity in large populations of PN neurons. Evoked responses are typically on the order of 10 – 100nAm, and are estimated to be generated by the synchronous spiking activity of the order of tens of thousands of pyramidal neurons. Low-frequency oscillations are larger in magnitude and are on the order of 100–1000 nAm, and are estimated to be generated by the subthreshold activity of on the order of a million pyramidal neurons (Jones et al., 2009; Jones et al., 2007; Murakami and Okada, 2006). While HNN is constructed with the ability to adjust local network size, the magnitude of these signals can also be conveniently matched by applying a scaling factor to the model output, providing an estimate of the number of neurons that contributed to the signal.

(II) Exogenous driving input through two known layer-specific pathways. One type of input represents excitatory synaptic drive that comes from the lemniscal thalamus and contacts the cortex in the granular layers, which then propagates to the proximal PN dendrites in the supragranular and infragranular layers and somata of the inhibitory neurons; this input is referred to as proximal drive (Figure 3B). The other input represents excitatory synaptic drive from higher-order cortex or non-specific thalamic nuclei that synapses directly into the supragranular layers and contacts the distal PN dendrites and somata of the inhibitory neuron; this input is referred to as distal drive (Figure 3C). The networks that provide proximal and distal input to the local circuit (e.g., thalamus and higher order cortex) are not explicitly modelled, but rather these inputs are represented by simulated trains of action potentials that activate excitatory post-synaptic receptors in the local network. The temporal profile of these action potentials is adjustable depending on the simulation experiment and can be represented as single spikes, bursts of input, or rhythmic bursts of input. There are several ways to change the pattern of action potential drive through different buttons built into the HNN GUI: Evoked Inputs, Rhythmic Proximal Inputs, and Rhythmic Distal Inputs. The dialog boxes that open with these buttons allow creation and adjustment of patterns of evoked response drive or rhythmic drive to the network (see tutorials described in Results section for further details).

(III) Exogenous drive to the network can also be generated as excitatory synaptic drives following a Poisson process to the somata of chosen cell classes or as tonic input simulated as a somatic current clamp with a fixed current injection. The timing and duration of these drives is adjustable.

Further details of the biophysics and morphology of the cells and of the architecture of the local synaptic connectivity profiles in the template network can be found in the Materials and methods section. As the use of our software grows, we anticipate other cells and network configurations will be made available as template models to work with via open source sharing (see Discussion).

HNN’s template model is a large-scale model simulated with thousands of differential equations and parameters, making the parameter optimization process challenging. The process for tuning this canonical model and constraining the space of parameters to investigate the origin of ERPs and low-frequency oscillations was as follows. First, the individual cell morphologies and physiologies were constrained so individual cells produced realistic spiking patterns to somatic injected current (detailed in Materials and methods). Second, the local connectivity within and among cortical layers was constructed based on a large body of literature from animal studies (detailed in Materials and methods). All of these equations and parameters were then fixed, and the only parameters that were originally tuned to simulate ERPs and oscillations were the timing and the strength of the exogenous drive to the local network. This drive represented our ‘simulation experiment’ and was based on our hypotheses on the origin of these signals motivated by literature and on matching model output to features of the data (see tutorials described in Results). The HNN GUI was constructed assuming ERPs and low-frequency oscillations depend on layer-specific exogenous drives to the network. The simulation experiment workflow and tutorials described below are in large part based on ‘activating the network’ by defining the characteristics of this layer-specific drive. Default parameter sets are provided as a starting point from which the underlying parameters can be interactively manipulated using the GUI, and additional exogenous driving inputs can be created or removed.

Automated parameter optimization is also available in HNN and is specifically designed to accurately reproduce features of an ERP waveform based on the temporal spacing and strength of the exogenous driving inputs assumed to generate the ERP. Before taking advantage of HNN’s automated parameter optimization, we strongly encourage users to begin by understanding our ERP tutorial and by hand-tuning parameters using one of our default parameter sets to get an initial representation of the recorded data. The identification of an appropriate number of driving inputs and their approximate timings and strengths serves as a starting point for the optimization procedure (described in the ERP Model Optimization section below). Hand tuning of parameters and visualizing the resultant changes in the GUI will enable users to understand how specific parameter changes impact features of the current dipole waveform.

Importantly, the biophysical constraints on the origin of the current dipoles signal (discussed above) will dictate the output of the model and necessarily limit the space of parameter adjustments that can accurately account for the recorded data. The same principle underlies the fact that a limited space of signals are typically studied at the macroscale (ERPs and low-frequency oscillations). A parameter sensitivity analysis on perturbations around the default ERP parameter sets confirmed that a subset of the parameters have the strongest influence on features of the ERP waveform (see Supplementary Materials). Insights from GUI-interactive hand tuning and sensitivity analyses can help narrow the number of parameters to include in the subsequent optimization procedure and greatly decrease the number of simulations required for optimization.

We have applied HNN to study the neural origin of tactile evoked responses localized with inverse methods to primary somatosensory cortex from MEG data (Jones et al., 2007). In this study, the tactile evoked response was elicited from a brief perceptual threshold level tap - stimulus strength maintained at 50% detection - to the contralateral middle finger tip during a tactile detection experiment (experimental details in Jones et al., 2009 and Jones et al., 2007). The average tactile evoked response during detected trials is shown in Figure 4. The data from this study is distributed with HNN installation.

Following the workflow described above, the process for reproducing these results in HNN is as follows.

Load the evoked response data distributed with HNN, ‘yes_trial_SI_ERP_all_avg.txt’. The data shown in Figure 4B will be displayed. Adjust parameters defining the automatically loaded default local network, if desired.

‘Activate’ the local network. In prior publications, we showed that this tactile evoked response could be reproduced in HNN by ‘activating’ the network with a sequence of layer-specific proximal and distal spike train drive to the local network, which is distributed with HNN in the file ‘ERPYes100Trials.param’.

The sequence described below was motivated by intracranial recordings in non-human primates, which guided the initial hypothesis testing in the model. Additionally, we established with inverse methods that at the prominent ~70 ms negative peak (Figure 4D), the orientation of the current was into the cortex (e.g., down the pyramidal neuron dendrites), consistent with prior intracranial recordings (see Jones et al., 2007). As such, in this example, negative current dipole values correspond to current flow down the dendrites, and positive values up the dendrites. In sensory cortex, the earliest evoked response peak corresponds to excitatory synaptic input from the lemniscal thalamus that leads to current flow out of the cortex (e.g., up the dendrites). This earliest evoked response in somatosensory cortex occurs at ~25 ms. The corresponding current dipole positive peak is small for the threshold tactile response in Figure 4D, but clearly visible in Figure 11 for a suprathreshold (100% detection) level tactile response.

The drive sequence that accurately reproduced the tactile evoked response consisted of ‘feedforward’/proximal input at ~25 ms post stimulus, followed by ‘feedback’/distal input at ~60 ms, followed by a subsequent ‘feedforward’/proximal input at ~125 ms (Gaussian distribution of input times on each simulated trial, Figure 4C). This ‘activation’ of the network generated spiking activity and a pattern of intracellular dendritic current flow in the pyramidal neuron dendrites in the local network to reproduce the current dipole waveform, many features of which fell naturally out of the local network dynamics (details in Jones et al., 2007). This sequence can be interpreted as initial ‘feedforward’ input from the lemniscal thalamus followed by ‘feedback’ input from higher-order cortex or non-lemniscal thalamus, followed by a re-emergent leminsical thalamic drive. A similar sequence of information flow likely applies to most sensory evoked signals. The inputs are distinguished with red and green arrows (corresponding to proximal and distal input, respectively) in the main GUI window. The number, timing, and strength (post-synaptic conductance) of the driving spikes were manually adjusted in the model until a close representation of the data was found (see section on parameter tuning above). To account for some variability across trials, the exact time of the driving spikes for each input was chosen from a Gaussian distribution with a mean and standard deviation (see Evoked Inputs dialog box, Figure 4C, and green and red histograms on the top of the GUI in Figure 4D). The gray curves in Figure 4D show 25 trials of the simulation (decreased from 100 trials in the Set Parameters, Run dialog box) and the black curve is the average across simulations. The top of the GUI windows displays histograms of the temporal profile of the spiking activity providing the sequence of proximal (red) and distal (green) synaptic input to the local network across the 25 trials. Note, a scaling factor was applied to net dipole output to match to the magnitude of the recorded ERP data and used to predict the number of neurons contributing to the recorded ERP. This scaling factor is chosen from Set Parameters, Run dialog box, and is shown as 3000 on the y-axis of the main GUI window in Figure 4D. Note that the scaling factor is used to predict the number of pyramidal neurons contributing to the observed signal. In this case, since there are 100 pyramidal neurons in each of layers 2/3 and 5, that amounts to 600,000 neurons (200 neurons x 3000 scaling factor) contributing to the evoked response, consistent with the experimental literature (described in Jones et al., 2009 and Jones et al., 2007).

Based on the assumption that sensory evoked responses will be generated by a layer-specific sequence of drive to the local network similar to that described above, HNN’s GUI was designed for users to begin simulating evoked responses by starting with the aforementioned default sequence of drive that is defined when starting HNN and by loading in the parameter set from the ‘ERPYes100Trials.param’ file, as described above. The Evoked Inputs dialog box (Figure 4C) shows the parameters of the proximal and distal drive (number, timing, and strength) used to produce the evoked response in Figure 4D. Here, there were two proximal drives and one distal drive to the network. These parameters were found by first hand tuning the inputs to get a close representation of the data and then running the parameter optimization procedure described below.

The evoked response shown in Figure 4 is reproduced by clicking the ‘Run Simulation’ button at the top of the GUI, and the RMSE of the goodness of fit to the data is automatically calculated and displayed. Additional network features can also be visualized through pull down menus (Step 5).

Evoked response parameters can now be adjusted, and additional inputs can be created or removed to account for the user-defined ‘simulation experiment’ and hypothesis testing goals (Step 6). With each parameter change, a new parameter file will be saved by renaming the simulation under ‘Simulation Name’ in the ‘Set Parameters’ dialog box (see Figure 4C). From here, other cell or network parameters can be adjusted to compare across conditions (Step 7).

We have applied HNN to study the neural origin of spontaneous rhythms localized to the primary somatosensory cortex from MEG data; it is often referred to as the mu-rhythm, and it contains a complex of (7–14 Hz) alpha and (15–29 Hz) beta frequency components (Jones et al., 2009). A 1 s time frequency spectrogram of the spontaneous unaveraged SI rhythm from this study is shown in Figure 6a. This data is distributed on the HNN website (‘SI_ongoing.txt’), and contains 1000 1 s epochs of spontaneous data (100 trials each from 10 subjects). The data is plotted in HNN through the ‘View → View Spectrograms’ menu item, followed by ‘Load Data’ and then selecting the ‘SI_ongoing.txt’ file. Note that it may take a few minutes to calculate the wavelet transforms for all 1000 1 s trials included. Next, select an individual trial (e.g. trial 32) from the drop-down menu. The dipole waveform from a single 1 s epoch will then be shown in the top.

Figure 6

Download asset Open asset

Figure 6—source data 1

S1 pre-stimulus activity.

https://cdn.elifesciences.org/articles/51214/elife-51214-fig6-data1-v2.txt

Download elife-51214-fig6-data1-v2.txt

The corresponding time-frequency spectrogram is automatically calculated and displayed at the bottom, as seen in Figure 6A. The default colormap indicating spectral power is the ‘jet’ scheme with blue representing low power and red representing high power. This is the same colormap used in the spectrograms of prior publications using the HNN model for studying low-frequency oscillations. However, HNN allows the user to plot the spectrogram using different standardized colormaps that are perceptually uniform, meaning they have a lightness value that increases monotonically (Pauli, 1976). The configuration option for changing the colormap and the colormap options are shown in Figure 6B. The spectrogram resulting from choosing the perceptually uniform ‘viridis’ colormap is shown in Figure 7C. Note that updating the spectrogram colormap requires repeating the ‘View → View Spectrograms’ menu selection.

Figure 7

Download asset Open asset

Notice that this rhythm contains brief bouts of alpha or beta activity that will occur at different times in different trials due to the spontaneous, non-stationary nature of the signals. Such non-time-locked oscillations are often referred to as spontaneous and/or ‘induced’ rhythms. When the non-negative spectrograms are averaged across trials, the intermittent bouts of high power alpha and beta activity accumulate without cancellation, and bands of alpha and beta activity appear continuous in the spectrogram (data not shown, see Jones, 2016; Jones et al., 2007) and will create peaks at alpha and beta in a power spectral density (Figure 8C). Since the alpha and beta components of this rhythm are not time locked across trials, it is difficult to directly compare the waveform of the recorded data with model output. Rather, to assess the goodness of fit of the model, we compared features of the simulated rhythm to the data (see Jones et al., 2009), including peaks in the power spectral density, as described below. Since we can not directly compare the waveform of this rhythm with the model output, rather than first loading the data, we begin this tutorial with Step 3, ‘activating’ the network, using the default local network defined when starting HNN.

Figure 8

Download asset Open asset

‘Activate’ the local network. In prior publications, we have simulated non-time-locked spontaneous low-frequency alpha and beta rhythms through patterns of rhythmic drive (repeated bursts of spikes) through proximal and distal projection pathways. These patterns of drive were again motivated by literature and by tuning the parameters to match features of the model output to the recorded data (see Jones et al., 2009; Sherman et al., 2016).

We begin by describing the process for simulating a pure alpha frequency rhythm only, and we then describe how a novel prediction for the origin of beta events emerged (Sherman et al., 2016). Motivated by a long history of research showing alpha rhythms in neocortex rely on ~10 Hz bursting in the thalamus, we tested the hypothesis that ~10 Hz bursts of drive through proximal and distal projection pathways (representing lemniscal and non-lemniscal thalamic drive) could reproduce an alpha rhythm in the local circuit. The burst statistics (number of spikes and inter-burst interval chosen from a Gaussian distribution), strength of the input (post-synaptic conductance), and delay between the proximal and distal input were manually adjusted until a pure alpha rhythm sharing feature of the data was found. We showed that when ~ 10 Hz bursts of proximal and distal drives are subthreshold and arrive to the local network in anti-phase (~50 ms delay) a pure alpha rhythm emerges (Jones et al., 2009; Ziegler et al., 2010).

The parameters of this drive are distributed with HNN in the file ‘Alpha.param’, loaded through the Set Parameters From File button and viewed in the Set Parameters dialog box under Rhythmic Proximal and Rhythmic Distal inputs (Figure 7A). Note that the start time mean of the ~10 Hz Rhythmic Proximal and Rhythmic Distal Inputs are delayed by 50 ms. The HNN GUI in Figure 7B displays the simulated current dipole output from this drive (middle), the histogram of the proximal and distal driving spike trains (top), and the corresponding time-frequency domain response (bottom). This GUI window is automatically constructed when rhythmic inputs are given to the network, and HNN is designed to easily define rhythmic input to the network via the Set Parameters dialog box. A scaling factor was also applied to this signal (via Set Parameters, Run dialog box) and is shown as 300,000 on the y-axis of the main GUI window example in Figure 7B. The 300,000 scaling factor predicts that 60,000,000 PNs (300,000 × 200 PNs) contribute to the measured signal.

The alpha rhythm shown in Figure 7B is reproduced by clicking the ‘Run Simulation’ button at the top of the GUI, Additional network features, including power-spectral density plots, can also be visualized through the pull down menus (Step 5).

Rhythmic input parameters can be adjusted to account for the user defined ‘simulation experiment’ and hypothesis testing goals.

The goal in our prior study was to reproduce the alpha/beta complex of the SI mu-rhythm. By hand tuning the parameters we were able to match the output of the model to several features of the recorded data, including symmetric amplitude modulation around zero and PSD plots as shown in Figures 7 and 8 (see further feature matching in Jones et al., 2009 and Sherman et al., 2016), we arrived at the hypothesis that brief bouts of beta activity (‘beta events’) non-time locked to alpha events could be generated by decreasing the mean delay between the proximal and distal drive to 0 ms and increasing the strength of the distal drive relative to the proximal drive. This parameter set is also distributed with HNN (‘AlphaandBeta.param’) and viewed in Figure 8A. With this mechanism, beta events emerged on cycles when the two stochastic drives hit the network simultaneously and when the distal drive was strong enough to break the upward flowing current and create a prominent ~50 ms downward deflection (see red box in Figure 8B). The stronger the distal drive the more prominent the beta activity (data not shown, see Sherman et al., 2016). This beta event hypothesis was purely model derived and was based on matching several features of the SI mu rhythm between the model output and data (detailed in Jones et al., 2009 Sherman et al., 2016).

Importantly, due to the non-time locked nature of this spontaneous rhythm, the waveform can not be directly compared by overlaying the waveform of the model and recorded oscillations as in the evoked response example (e.g. Figure 4). However, one can quantify features of the oscillation and compare to recorded data (Jones et al., 2009; Sherman et al., 2016). One such feature is the amplitude of the oscillation waveform, where a scaling factor can be applied to the model to predict how many cells are needed to produce a waveform amplitude on the same order as the recorded data, as described in Step three above. Additionally, the PSD from the model and data can be directly compared. This can be viewed in the HNN GUI though the ‘View PSD’ pull down menu (see Figure 4, Step 5), where this data (‘SI_ongoing.txt’ - provided with HNN) can be automatically compared to the model output in the PSD window (Figure 8C).

The model derived predictions on mechanisms underlying alpha and beta where motivated by literature and further refined by tuning the parameters to match the output of the model with various features of the recorded data. While the mechanisms of the alpha rhythm described above were motivated by literature showing cortical alpha rhythms arise in part from alpha frequency drive from the thalamus and supported by animal studies (Hughes and Crunelli, 2005; for example, see Figure 2 in Bollimunta et al., 2011), the beta event hypothesis was novel. The level of circuit detail in the model led to specific predictions on the laminar profile of synaptic activity occurring during beta events that could be directly tested with invasive recordings in animal models. One specific prediction was that the orientation of the current during the prominent ~50 ms deflection defining a beta event (red box, Figure 8B) was down the pyramidal neuron dendrites (e.g. into the cortex). This prediction, along with several others, were subsequently tested and validated with laminar recordings in both mice and monkeys, where it was also confirmed that features of beta events are conserved across species and recording modalities (Sherman et al., 2016; Shin et al., 2017).

Gamma rhythms can encompass a wide band of frequencies from 30 to 150 Hz. Here, we will focus on the generation of so-called ‘low gamma’ rhythms in the 30–80 Hz range. It has been well established through experiments and computational modeling that these rhythms can emerge in local spiking networks through excitatory and inhibitory cell interactions where synaptic time constants set the frequency of the oscillation, while broadband or ‘high gamma’ rhythms reflect spiking activity in the network that creates sharp waveform deflections (Lee and Jones, 2013). The period of the low gamma oscillation is determined by the time constant of decay of GABAA-mediated inhibitory currents (Buzsáki and Wang, 2012; Cardin et al., 2009; Vierling-Claassen et al., 2010), a mechanism that has been referred to as pyramidal-interneuron gamma (PING). In normal regimes, the decay time constant of GABA_A-mediated synapses (~25 ms) bounds oscillations to the low gamma frequency band (~40 Hz). In general, PING rhythms are initiated by ‘excitation’ to the excitatory (PN) cells, and this initial excitation causes PN spiking that, in turn, synaptically activates a spiking population of inhibitory (I) cells. These (I) cells then inhibit the PN cells, preventing further PN activity until the PN cells can overcome the effects of the inhibition ~25 ms later. The pattern is repeated, creating a gamma frequency oscillation (~40 Hz).

We have applied HNN to determine if features in the current dipole signal could distinguish PING-mediated gamma from other possible mechanisms such as exogenous rhythmic drive or spiking activity that creates ‘high gamma’ oscillations (Lee and Jones, 2013). Here, we describe the process for generating gamma rhythms via the canonical PING mechanisms in HNN. First, to demonstrate the basic mechanisms of PING, we simulate a robust large amplitude and nearly continuous gamma rhythm (Steps 3–5, Figure 9). Second, we adjust the simulation parameters in order to directly compare to experimental data, where on single trials induced gamma rhythms are smaller in amplitude and less continuous emerging as transient bursts of activity (Pantazis et al., 2018) (Steps 6 and 7, Figure 10). Both parameter sets for creating the examples below are distributed with the software.

Figure 9

Download asset Open asset

Figure 10

Download asset Open asset

Figure 10—source data 1

Visual cortex gamma activity from 1 trial.

https://cdn.elifesciences.org/articles/51214/elife-51214-fig10-data1-v2.txt

Download 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

Download elife-51214-fig10-data2-v2.txt

To demonstrate the robust PING-mechanisms and its expression at the level of a current dipole, we begin this tutorial with Step 3, ‘activating’ the network using a slightly altered local network configuration as described below.

‘Activate’ the local network by loading in the parameter set defining the local network and initial input parameters ‘gamma_L5weak_L2weak.param’. In this example, the input was noisy excitatory synaptic drive to the pyramidal neurons. Additionally, all synaptic connections within the network are turned off (synaptic weight = 0) except for reciprocal connections between the excitatory (AMPA only) and inhibitory (GABA_A only) cells within the same layer. This is not biologically realistic, but was done for illustration purposes and to prevent pyramidal-to-pyramidal interactions from disrupting the gamma rhythm. To view the local network connections, click the ‘local network’ button in the Set Parameters dialog box. Figure 9B shows the corresponding dialog box where the values of adjustable parameters are displayed. Notice that the L2/3 and L5 cells are not connected to each other, the inhibitory conductance weights within layers are stronger than the excitatory conductances, and there are also strong inhibitory-to-inhibitory (i.e., basket-to-basket) connections. This strong autonomous inhibition will cause synchrony among the basket cells, and hence strong inhibition onto the PNs.

To reproduce the ~40 Hz gamma oscillation described by the PING mechanism above, we drove the pyramidal neuron somas in L2/3 and L5 with noisy excitatory AMPA synaptic input, distributed in time as a Poisson process with a rate of 140 Hz. This noisy input can be viewed in the ‘Set Parameters’ menu by clicking on the ‘Poisson Inputs’ button (see Figure 9A). Setting the stop time of the Poisson drive to −1, under the Timing tab, keeps it active throughout the simulation duration.

The gamma rhythm shown in Figure 9C is reproduced by clicking the ‘Run Simulation’ button at the top of the main HNN GUI. The top panel shows a histogram of Poisson distributed times of input to the pyramidal neurons, the middle panel the net current dipole across the entire network, and the bottom the corresponding time-frequency spectrogram showing strong gamma band activity. Additional network features, including spiking activity in each cell in the population (Figure 9D), somatic voltages (Figure 9E), and PSD plots for each layer and the entire network (Figure 9F) can also be visualized through the ‘View’ pull down menu (Step 5). Notice the PING mechanisms described above in the spiking activity of the cells (Figure 9D), where in each layer the excitatory pyramidal neurons fire before the inhibitory basket cells. The line plots, which show spike counts over time, also demonstrate rhythmicity. The pyramidal neurons are firing periodically but with lower synchrony due to the Poisson drive (orange histogram at the top), which creates randomized spike times across the populations (once the inhibition sufficiently wears off). Notice also that the power in the gamma band is much smaller in Layers 2/3 than in Layer 5 (Figure 9F). This is reflective, in part, of the fact that the length of the L2/3 PNs is smaller than the L5 PNs, and hence the L2/3 cells produce smaller current dipole moments that can be masked by activity in Layer 5. This example was constructed to describe features in the current dipole signal that could distinguish PING mediated gamma rhythms but has not been directly compared to recorded data (see Lee and Jones, 2013 for further discussion).

Next, we describe how parameter adjustments to the canonical PING mechanism above can account for induced gamma rhythms recorded with MEG in a prior published study (Pantazis et al., 2018). In this study, subjects were presented with the following visual stimuli: stationary square-wave Cartesian gratings with vertical orientation and black/white maximum contrast with a frequency of 3 cycles/degree. These stimuli create non-time locked induced gamma rhythms in the visual cortex; a reliable phenomena observed in many studies. The data represents time-domain activity from sources localized in the pericalcarine region of the visual cortex from 0 to 700 ms after the presentation of the stimulus (Figure 10A,B), and is located in HNN’s data/gamma_tutorial folder.

To view the data in HNN, click on HNN’s ‘View’ pull down menu, then click on ‘View Spectrograms’. Then from the Spectrogram Viewer, click on the ‘File’ pull down menu, click ‘Load Data’, and select data/gamma_tutorial/100_trials.txt. Figure 10A shows a single trial of the source-localized current dipole signal, with corresponding wavelet-transform spectrogram below. Here, visual gratings were provided from time = 0–700 ms. On single non-averaged trials, the visual stimulus creates non-time-locked induced gamma rhythms at ~50 Hz that are not continuous but transient in time and ‘burst-like’, with each burst lasting ~50–100 ms. HNN’s Spectrogram Viewer automatically averages the wavelet spectrograms across trials when the data is loaded (Figure 10B). To view a single trial as in Figure 10A, use the drop-down menu at the bottom of the Spectrogram Viewer to select a single trial. The averaged spectrogram from 100 trials makes the oscillation appear nearly continuous in the 40–50 Hz gamma band (Figure 10B).

We hypothesized that in order to change the nearly continuous gamma oscillation from HNN’s canonical PING model (Figure 9) into a bursty lower amplitude gamma rhythm on single trials, the noisy excitatory AMPA synaptic input to the pyramidal neurons that creates the PING-rhythm, as described in Step 3, would need to be reduced. We tested this hypothesis by reducing the rate of the Poisson synaptic inputs to L2/3 and L5 pyramidal neurons, here from 140 Hz to 4 Hz and 5 Hz, respectively. To see the parameter values used, load the ‘gamma_L5weak_L2weak_bursty.param’ file, then click on ‘Set Parameters’ and ‘Poisson Inputs’ (parameter values not shown in Figure 10). Note, that in addition to the change and rate of drive, the NMDA synaptic inputs now have a small positive weight (previously zero), which was required to compensate for the reduced AMPA activation of the pyramidal neurons.

Click on ‘Run Simulation’ to produce the results shown in Figure 10C. The top panel of Figure 10C shows the Poisson inputs, which have a noticeably lower rate than the default PING simulation. The current dipole time-course is shown in the middle panel. The corresponding wavelet spectrogram in the bottom panel shows intermittent gamma bursts recurring with high power, similar to the features seen in the experiment. Note, in this data set, the amplitude of the current dipole signal is small,<0.1 nAm. As such, a small dipole scaling factor of 5 was applied to the output of HNN to compare to the amplitude of the experimental data. This predicts that a highly localized population of ~1000 pyramidal neurons contribute to the recorded gamma signal.

To see the effects of averaging across trials in HNN, click on ‘Set Parameters’ and ‘Run’, enter 100 trials, and then click on ‘Run Simulation’. Output as shown in Figure 10D will be produced. The individual current dipole waveforms do not have consistent phase across trials (Figure 10D middle). However, averaging the spectrograms across trials produces a more continuous band of gamma oscillation throughout the simulation duration, similar to the effects observed in the experiment.

We can now use HNN to take a closer look at the underlying neuronal spiking activity contributing to the observed dipole signal, as shown in Figure 10E. The firing rates of L2/3 pyramidal neurons (green) and L5 pyramidal neurons (red) are significantly lower and more sparse than in the continuous PING simulation shown in Figure 9. This lower pyramidal neuron activation leads to fewer, or no basket interneurons firing on any given gamma cycle (Figure 10E white, blue points have fewer participating interneurons, and do not always occur). This lower interneuron activity produces lower-amplitude bouts of feedback inhibition, which are sometimes nearly absent, and thus creates transient bursty gamma activity. Note that the firing rates of L2/3 pyramids are lower than that of L5 pyramidals, consistent with experimental data (Naka et al., 2019; Schiemann et al., 2015) and previous data-driven modeling (Dura-Bernal et al., 2019; Neymotin et al., 2016). Layer specific current dipole responses in Figure 10F confirm that due to their dendritic length L5 pyramidals are once again the major contributors to the aggregate dipole signal. This example is presented as a proof of principle of the method to begin to study the cell and circuit origin of source localized gamma band data with HNN. The circuit level hypotheses have not yet been published elsewhere or investigated in further detail.

Similar to steps 1–4 above, first load the supra-threshold experimental data file ‘S1_SupraT.txt’ via the ‘Load data file’ menu option and the example starting parameters to activate the network provided in the parameter file ‘ERPYes100Trials.param’ via the ‘Load parameter file’ menu option. Note that in this example, the network is also ‘activated’ by a sequence of three exogenous inputs defined in the parameter file. The parameters for these inputs serve as a baseline for model optimization. The supplied parameter file (used above) runs 100 trials by default for each simulation. For model optimization, this can be reduced to three trials. Click on the ‘Set Parameters’ button, then the ‘Run’ button, and replace 100 trials with 3. In the previous Set Parameters dialog box change the simulation name to ‘ERPYes3Trials’ to reflect this change (Figure 4C). By clicking the ‘Run Simulation’ button the evoked response using this initial parameter set as in Figure 11A will be displayed. As described above, in practice with user defined data, users should apply their own hypotheses related to the number, timing and synaptic input strengths of the exogenous inputs that activate the network to obtain an initial representation of the recorded waveform before beginning the parameter estimation process.

Before running the optimization, rename the simulation to ‘ERPYes3Trials_opt’ in the Set Parameters dialog box as described above, so that the parameter results of the optimization will be saved in a new file.

In the Simulation pull down menu, choose the ‘Configure Optimization’ option. This option is only selectable once data and parameter files have been loaded. A new dialog box pre-populated with values from the parameter file will appear, as shown in Figure 11B. All parameters describing the timing and strength of defined exogenous inputs will be available for optimization. Users can generate their own evoked response parameter files with as many exogenous inputs as desired and they will be automatically populated into the ‘Configure Optimization’ dialog box.

Select which parameters to treat as free variables for optimization; parameters that will be fixed in the optimization process are grayed out. By default, all parameters are selected, but it may be desirable to limit the number of free parameters to only the most influential set based on a parameter sensitivity analysis. Fixing non-influential parameters will decrease the complexity of an optimization step, and increase the likelihood of the optimization algorithm converging on parameter estimates after a relatively low number of simulations. Results of a sensitivity analysis using Uncertainpy (Tennøe et al., 2018) on this example data are provided in Supplementary file 1 and may help guide model optimization for similar data. Sensitivity analysis is not yet included in HNN (see Future Directions).

The number of simulations per optimization step is configurable in the top section of the ‘Configure Optimization’ dialog box (Figure 11B). The default values shown in Figure 11B were based on results from our studies where the fit obtained was significantly improved from a single optimization. This value can be decreased as the number of free parameters is reduced.

The parameter ranges defining the bound constraints given to the optimization algorithm are shown in the ‘Defined range’ column of the dialog box in Figure 11B. The displayed range is calculated as plus or minus a specified number of standard deviations for input start time or plus or minus a percentage of the initial value for all other parameters. The user may customize the range by inputting their own ‘Range specifier’ or using the interactive slider bar to define new minimum or maximum values. If a parameter has an initial value of 0, its range is defined by a user-specified maximum value rather than percentage. The ‘Reset Ranges’ button will update ranges using the ‘Range specifier’ and discard custom values set by the slider. Parameters can be fixed during the optimization process by unchecking the ‘Optimize’ checkbox.

Click the ‘Run Optimization’ button to start the stepwise optimization process. After each input has been optimized in sequence followed by a final optimization step that adjusts all input parameters, the final optimized fit will be shown in gray in the main HNN window along with the lowest obtained RMSE (Figure 11C/D).

To perform a second optimization using the results of the first procedure as a starting point, select the optimized simulation parameter set drop-down menu. This will update the values in the Configure Optimization dialog box and pressing Run Optimization will start a new optimization process. For this example, the RMSE improved from 15.54 (Figure 11D) after the first optimization to 10.02 after a second round (data not shown).

The template neocortical model provided with HNN is based on prior publications using the model without a graphical user interface, as described in Jones et al. (2009), and available on ModelDB (https://senselab.med.yale.edu/modeldb). All current parameter files included in the software, and described in the tutorials above, are based on this model, except for the gamma tutorial whose parameter file has local network modifications as described above.

HNN’s underlying neocortical model is simulated using the NEURON simulation environment with the Python interpreter (see Key Resources Table). HNN’s model is simulated across multiple cores in parallel using the message-passing interface (MPI). HNN’s Run Parameters dialog box can be accessed through the GUI and provides access to commonly used simulation parameters, including integration time-step (dt), simulation duration (milliseconds), number of trials, neuronal firing threshold (mV), and number of cores over which to parallelize the model.

The model represents a canonical neocortical circuit. It contains multi-compartment pyramidal neurons (PN) in supragranular and infragranular layers (layers 2/3 and 5, respectively), whose apical dendrites are spatially aligned and span the cortical layers. In both layers, the PNs have two basal, one oblique, and one apical dendrite branch, and the layer 2/3 PNs have shorter apical dendrites than layer 5 PNs.

The PNs are synaptically coupled to each other and to a subset of inhibitory neurons in each layer, and are included in the model in a 3/1 PN-to-interneuron ratio, with a scalable number of PNs. The inhibitory neurons are simulated with single compartments representing fast spiking basket cells, and are shown in yellow in (Figure 2). Note that the granular layer is not explicitly included in the template circuit. This design choice was based on the fact that macroscale current dipoles are dominated by PN activity in supragranular, and infragranular layers, due to their alignment (see Calculation of Primary Electrical Current). Thalamic input to granular layers is presumed to propagate directly to basal and oblique dendrites of PN in layer 2/3 and 5.

The morphology of the PN in each layer (see Figure 3, and Table 1) were adapted from the morphology reduction procedure in Bush and Sejnowski (1993). In Bush and Sejnowski (1993), digitized HRP-filled pyramidal neurons from layers 2 and 5 of cat visual cortex consisting of 400 compartments were reduced to 8 and 9 compartments, respectively. This procedure is important for HNN’s purposes because the axial conductance is the basis of the primary current dipole calculation (see below).

The reduction procedure was based on conserving the axial conductance between compartments in the cells rather than the surface area of the dendritic tree. This reduction retained the general morphology of the neurons and allowed for accurate position of synaptic inputs and ionic conductances on individual cells, and construction of spatially accurate network models. The reduced model contained active conductances in the somatic compartents and retaining faithful electrical responses to injected current. In HNN, a scaling factor of 1.3 was applied to the length and diameters of the dendritic compartments to account for increases in dendritic length and volume in human somatosensory neurons, as predicted by larger cortical thickness (Fischl and Dale, 2000; Geyer et al., 1997) and an increase in the number of dendritic spines and arborization (Elston et al., 2001). Additional active ionic conductances were added to the dendritic compartments, as detailed below, and the neurons were embedded in a cortical column network with single compartment inhibitory neurons. The radial geometry of the inhibitory neuron dendrites does not contribute to the primary current dipole calculation. The morphology of each cell was as follows:

• Layer 2/3

  • PN: 8 compartments including four apical dendrites, three basal dendrites, one soma

  • Inhibitory basket neurons: single compartment (soma)

• Layer 5

• • PN: 9 compartments including five apical dendrites, three basal dendrites, one soma.

  • As shown below, L5 PNs have longer dendrites than L2/3 PNs. L5 PN somas are based in L5 with long apical dendrites reaching into L2/3.

  • Inhibitory Basket neurons: single compartment (soma).

  • L2/3 and L5 basket interneurons are identical but their synaptic parameters and local circuit connectivity differs.

Evoked inputs are trains of synaptic inputs to the local network during a sensory stimulus that creates an event related potential (ERP). Parameter choices for defining these inputs are shown in Figure 4A. The following parameter values are used to define each proximal or distal evoked input:

• Start time mean (ms) - average start time

• Start time stdev (ms) - standard deviation of start time

• Number spikes - number of inputs provided to each synapse

• L2/3 Pyr weight AMPA/NMDA ($\mu S$) - weight of AMPA/NMDA synaptic inputs to layer 2/3 pyramidal neurons

• L2/3 Basket weight AMPA/NMDA ($\mu S$) - weight of AMPA/NMDA synaptic inputs to layer /32 basket cells

• L5 Pyr weight AMPA/NMDA ($\mu S$) - weight of AMPA/NMDA synaptic inputs to layer 5 pyramidal neurons

• L5 Basket weight AMPA/NMDA ($\mu S$) - weight of AMPA/NMDA synaptic inputs to layer 5 basket cells (only used for proximal inputs)

Each evoked input also has a ‘Synchronous Inputs’ option, indicating whether for a specific evoked proximal/distal input each neuron receives the input at the same time, or if instead each neuron receives the evoked input events independently drawn from the same distribution. Increment input (ms) indicates whether to increment the Start time of all evoked inputs on each trial. In the studies described above, the evoked input strengths are suprathreshold generating action potentials in the local network.

Rhythmic Inputs are typically bursts of action potentials that drive the local network rhythmically. Parameter choices for defining these inputs are shown in Figure 7A. Each rhythmic input is defined as a series of ‘population bursts’, consisting of a set number of ‘burst units’ which drive post-synaptic conductances in the local network with a set frequency and mean delay between proximal and distal projections. Rhythmic proximal and distal inputs target different cortical layers, as described above. HNN allows setting proximal and distal rhythmic synaptic input start/stop times and frequencies using the following specification:

• Start time mean (ms) - specifies the average start time for rhythmic inputs

• Start time stdev (ms) - specifies the standard deviation of start times for rhythmic inputs

• Stop time (ms) - specifies when the rhythmic inputs should be turned off

• Burst frequency (Hz) - average frequency of bursts

• Burst stdev (ms) - standard deviation of input events

• Spikes/burst - provides n synaptic events at each selected time

• Number bursts - number of times the full Burst sequence is repeated (each repeat adds variability and more inputs)

In addition, HNN’s Rhythmic Input dialog box allows setting the weights of the rhythmic synaptic inputs (units of conductance) to individual neuron types in layers 2/3 and 5, and adding synaptic delays (ms) before the neurons receive the synaptic inputs. In the studies described above, rhythmic inputs are set to sub-threshold synaptic strengths, and therefore do not lead to neuronal action potentials.

Tonic inputs are modeled as somatic current clamps with a fixed current amplitude (nA). These clamps can be used to adjust the resting membrane potential of a neuron, and bring it closer (with positive amplitude injection) or further from firing threshold (with a negative amplitude injection). Parameter choices for defining these inputs are shown in Figure 9A and include setting the current clamp amplitude, and start/stop time for each modeled neuron type separately.

Noisy Inputs are trains of action potentials that follow a Poisson Process and create excitatory AMPA or NMDA synaptic inputs to the somata of all neurons of a given type. Parameter choices for defining these inputs are shown in Figure 9A and include, setting the average frequency of the Poisson drive, synaptic strength to somatic AMPA or NMDA synapses, and start/stop times of all Poisson inputs.

To reduce the computational demands of performing model optimization of HNN ERP simulations, we used variance based sensitivity analysis to identify parameters that were less significant to the simulated dipole waveform. As discussed above, HNN’s model optimization feature focuses on estimating parameters of the exogenous driving inputs. Of those parameters, ones that did not vary model output significantly, as determined by sensitivity analysis, could be excluded from parameter estimation.

The method of variance based sensitivity analysis through Monte Carlo estimation (Sobol′, 2001) provides Sobol sensitivity indices that can be used to explain the relative contribution of individual parameters on model variance. The total Sobol sensitivity index for each parameter serves as a measure that represents that parameter’s contribution to the variance, and also the contributions resulting from interactions with other parameters being varied (Homma and Saltelli, 1996). So a parameter with a low total Sobol sensitivity index can be characterized as an insignificant contributor to variance and can be fixed at its default value during model optimization.

We used the Python package Uncertainpy (Tennøe et al., 2018) to perform sensitivity analyses of parameters belonging to the exogenous driving inputs in the perceptual threshold-level ('yes_trial_SI_ERP_all_avg.txt') and suprathreshold-level ('ERPYesSupraT.txt') evoked response examples provided with HNN and described in the Results (Figure 4 and Figure 11). These analyses were performed using a modified simulation interface to run the simulations in parallel on a high-performance computing cluster, which is not currently included with HNN distribution. The results from our sensitivity analyses are shown in Figure 4—figure supplement 1, Figure 11—figure supplement 2, and Figure 11—figure supplement 1. Each analysis consisted of varying all parameters (except input timing standard deviation) of a driving input over 55,000 simulations using a quasi-Monte Carlo method that sampled from parameter distributions we specified. The input time distribution was defined as a normal distribution with the same mean and standard deviation in the corresponding HNN parameter file. The values for various synaptic weights were chosen from a uniform distribution ranging from the default value plus or minus 500%. For synaptic weight parameters with a default value of 0, the uniform distribution ranged from 0 to 1.0.

We were interested in comparing the contribution of each parameter to the dipole waveform (in units of nAm) across the entire simulation time. However, since the calculation of the total Sobol sensitivity index is relative to the variance at each point of the time series, when variance changes over time, it is not appropriate to average the total Sobol sensitivity indices across the entire simulation (Alexanderian et al., 2020). Instead we computed a weighted total Sobol index at each point in time (weighted by std. deviation scaled to have a range from 0 to 1). The plots in Figure 4—figure supplement 1 and Figure 11—figure supplement 2 show weighted total Sobol indices for each parameter over the duration of the simulation. Supplementary file 1 ranks the parameters with the greatest contribution to model output using the arithmetic means of weighted total Sobol indices across the entire simulation, for each driving input.

The results from our sensitivity analyses of sensory evoked response examples illustrate that there are several candidate parameters for excluding from model optimization. Not surprisingly input timing is an important parameter to optimize. We also found that NMDA weights typically have a greater contribution to variance than AMPA weights, as do connections to Layer 5 neurons compared to Layer 2/3 neurons.

• Samuel A Neymotin
• Dylan S Daniels
• Blake Caldwell
• Robert A McDougal
• Nicholas T Carnevale
• Mainak Jas
• Christopher I Moore
• Michael L Hines
• Matti Hämäläinen
• Stephanie R Jones

• Samuel A Neymotin
• Dylan S Daniels
• Blake Caldwell
• Robert A McDougal
• Nicholas T Carnevale
• Mainak Jas
• Christopher I Moore
• Michael L Hines
• Matti Hämäläinen
• Stephanie R Jones

• Samuel A Neymotin

• Samuel A Neymotin