4.1 GUI Tutorial of Event Related Potentials (ERPs) Simulation

Getting Started

In order to understand the workflow and initial parameter sets provided with this tutorial, we must first briefly describe prior studies that led to the creation of the provided data and evoked response parameter set that you will work with. This tutorial is based on results from our 2007 study where we recorded and simulated tactile evoked responses source localized to the primary somatosensory cortex (SI) (Jones et al. 2007).

In our 2007 study, we investigated the early evoked activity (0-175 ms) elicited by a brief tap to the D3 digit and source localized to an an equivalent current dipole in the contralateral hand area of the primary somatosensory cortex (SI) (Jones et al. 2007). The strength of the tap was set at either suprathreshold (100% detection probability) or perceptual threshold (50% detection) levels (see Figure 1, left panel below). Note, to be precise, this data represents source localized event related field (ERF) activity because it was collected using MEG. We use the terminology ERP for simplicity, since the primary current dipoles generating evoked fields and potentials are the same.

We found that we could reproduce evoked responses that accurately reflected the recorded waveform in our neocortical model from a layer specific sequence of exogenous excitatory synaptic drive to the local SI circuit, see Figure 1right panel below. This drive 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). This sequence of drive generated spiking activity and intracellular dendritic current flow in the pyramidal neuron dendrites to reproduce the current dipole signal. This sequence of drive 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. Intracranial recordings in non-human primates motivated and supported this assumption (Cauller and Kulics 1991).

In our model, the exogenous driving inputs were simulated as predefined trains of action potentials (pre-synaptic spikes) that activated excitatory synapses in the local cortical circuit in proximal and distal projection patterns (i.e. feedforward, and feedback, respectively, as shown schematically in Figure 1 right, and in the HNN GUI Model Schematics). 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 (all other model parameters were tuned and fixed based on the morphology, physiology and connectivity within layered neocortical circuits) (Jones et al. 2007). 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 (purple circle, Figure 1, right panel). The dipole units were in nAm, with a one-to-one comparison between data and model output due to the biophysical detail in our model.

Adapted from Jones et al. 2007 (Jones et al. 2007). Comparison of SI evoked response in experiment and neural model simulation. Left: MEG data showing tactile evoked response (ERP) source localized to the hand area of SI. Red: suprathreshold stimulation; Blue: Threshold stimulation (avg. n=100 trials). Right: Neural model simulation depicting proximal/distal inputs needed to replicate the ERP waveform (avg. n=25 trials)

In summary, to simulate the SI evoked response, a sequence of exogenous excitatory synaptic drive was simulated (by creating presynaptic spikes that activate layer specific synapses in the neocortical network) consisting of proximal drive at ~25 ms, followed by distal drive at ~60 ms, followed by a second proximal drive at ~122 ms. Given this background information, we can now walk you through the steps of simulating a similar ERP, using a subset of the data shown in Figure 1.

Tutorial Table of Contents

1. Load/view data

2. Load/view parameters to define network structure and to activate the network

3. Running the simulation and visualizing net current dipole

4. A closer look inside the simulations: contribution of layers and cell types

5. Comparing model output and recorded data

1. Load/view data

An example ERP dataset is provided in the hnn-data GitHub repository. We use the data as an orienting example for where to begin to simulate an ERP.

This dataset represents early evoked activity (0-175 ms) from an equivalent current dipole source localized to the hand area of the primary somatosensory cortex (SI), elicited by a brief perceptual threshold level tap to the contralateral D3 digit (read Getting Started above for details). The example dataset provided was collected at 600Hz and contains only averaged data from 100 trials in which the tap was detected. (Note, when loading your own data, if it was not collected at 600Hz, you must first downsample to 600Hz to view it in the HNN GUI).

To load and view this data, navigate to the main GUI window and on the bottom left corner click: Load data

If you have cloned the hnn-data repository, navigate to hnn-data folder on your desktop and select MEG_detection_data/yes_trial_S1_ERP_all_avg.txt. HNN will then load the data and display the waveform in the dipole window as shown below.

Alternatively, if you have not cloned the hnn-data repository, you can download the file directly by clicking here.

Figure 2

Note, the software can be used without loading data. If you wish to play with simulations without data, proceed to Step 2 first.

2. Load/view parameters to define network structure & to “activate” the network

An initial parameter set that will simulate the evoked drives that generate an evoked response in close agreement with the SI data described in Step 1 is distributed in the hnn-data repository. Click on the External drives tab at the top of the GUI and then click the Load external drives button. Navigate to the hnn-data/network-configurations folder on your computer and select ERPYes100Trials.json, or download the parameter file and load it into the GUI.

Figure 3

The template cortical column networks structure for this simulation is described in the “HNN Template Model” page on the hnn.brown.edu website. Several of the network parameter can be adjusted via the HNN GUI (e.g., local excitatory and inhibitory connection strengths) under the Network connectivity tab, but we will leave them fixed for this tutorial and only adjust the inputs that “activate” the network.

The values of the parameters that you loaded to “activate” the network in a manner that will generate an evoked response can now be viewed under the External drives tab. As described in the “Getting Started” section, the evoked response can be simulated with a sequence of exogenous driving inputs consisting of a proximal input at ~26 ms (evprox1), followed by a distal input at ~64 ms (evdist1), followed by a subsequent proximal input at ~137 ms (evprox).

To see the detailed parameter values defining each of these drives click on the dropdown button next to the name of each drive. Note: additional evoked proximal or distal inputs can be added to your simulation for your hypothesis testing goals by using the Add drive button and specifying the drive as “Evoked” and the location as either “proximal” or “distal”. Other types of drives can also be defined including poisson, rhythmic, and tonic, as detailed in other tutorials.

Figure 4

Each evoked input consists of a Gaussian distributed train of presynaptic action potentials that will target all of the cells in the post-synaptic network, with several adjustable parameters, including the mean arrival and standard deviation of the time each spike activates the network (in milliseconds), the number of the driving spikes on each trial of the simulation, and a random number seed that enables reproducibility of simulation results across trials. You can also adjust the postsynaptic conductance of the drive onto the postsynaptic cell. For example, under the “AMPA weights” section in the “evprox1” dropdown menu, the L5_pyramidal field represents the the post synaptic AMPA conductance of the proximal input onto the layer 5 pyramidal neuron at each location targeted by the proximal drive. For further details on the connectivity structure of the network, see the “HNN Template Model” page.

3. Running the simulation and visualizing net current dipole

Now that we have an initial parameter set, we can run a simulation for a set number or trials. Let’s start by defining 3 trials, by clicking on the Simulation tab and defining Trials=3. On each simulated trial, the timings of the evoked inputs (i.e., spikes) are chosen from a Gaussian distribution with mean and stdev (standard deviation) as defined in the “Evoked drives” tap. Histograms of each of the evoked inputs will be displayed at the top of the Figure tab after the simulations run.

Before running the simulation, we’ll first change the simulation name (i.e., the name under which the simulated data will be saved) to a new descriptive-name for the simulation here. Under the Simulation tab, change the name from “default” to “ERPyes-3trials”. (Note that the default simulation is in fact 1 run of the “ERPyes” simulation.) There are several other adjustable simulate parameters in the Simulation tab. These parameters control the duration (stop), integration time step (dt), number of trials (Trials), and the choice of the simulation backend of either MPI (parallel across neurons) or Joblib (parallel across trials), assuming both backends are installed.

Hit the Run button to run the simulation. A simulation log is shown under the Run button that will tell you the status of your simulaion.

Figure 5

Once complete, a new Figure 2 window will appear showing the output of the simulation as in the figure below. The thin blue traces are net current dipole signals from individual trials while the thick blue trace is the average ERP, with histograms of the proximal and distal driving spikes shown above.

Next, let’s create a new figure to view the simulation on top of the loaded data. Click on the Visualization tab and select the [blank] 2row x 1col (1:3) option from the Layout template dropdown menu, then select Make figure.

You will see that our new figure has two subplots defined by ax0 and ax1, both of which are blank. Select the ax1 subplot and set the Simution Data to ERPyes-3trials. Set Data to compare to yes_trial_S1_ERP_all_avg, then click Add plot.

Next, select ax0, set the Type to input histogram, and click Add plot.

Note you can change what is shown in either of these subplots by selecting clear axis, picking the Type of data from the pulldown menu, and clicking Add plot.

Figure 6

Importantly, note that a scaling factor of 3000.00 was multiplied by the net dipole produced by the model, as seen on the y-axis scale in the figure above. This scaling factor can be adjusted to match the magnitude of the recorded data; the value of 3000 is the default value for the loaded parameter set. In this case, since the template model contains 200 pyramidal neurons (PNs), the simulation predicts that the number of cells that contribute to the signal is 600,000 (200 x 3000) PNs.

Also note that in the ERP simulation shown, the raw dipole signal was smoothed using a Hamming filter using a window size of 30 milliseconds, in order to reduce noise in the ERP signal generated by this reduced network model. The level of smoothing can be adjusted by changing the value of the Dipole Smooth Window (ms). The longer the window, the more smoothing will occur. To turn off smoothing entirely, set the window size to 0. Below, we provide an example of the same simulation with smoothing turned off entirely. Note the higher-frequency content compared to the ERP simulation with smoothing turned on.

Figure 7

4. A closer look inside the simulations: contribution of layers and cell types

One of the main advantages of simulating neocortical activity is that we can dive into the details of the simulation to investigate the contribution of different components in the network; e.g., layers, cell types, etc. HNN currently enables the viewing of the following:

1. Layer-specific dipole activity
2. Spiking activity in each individual neuron population

4.1 Viewing layer specific current dipoles

From the Visualization tab, select “Dipole Layers (3x1)” from the Layout template dropdown menu. From the Dataset dropdown, select the simulation data you would like to visualize. Next, click Make figure to generate a new figure.

Figure 8

The new figure shows the dipole contributions from Layer 2/3 (top), Layer 5 (middle), and the aggregate (bottom). Note the different features in Layer 2/3 vs Layer 5 dipole signals, allowing you to tease apart how the different cortical layers contribute to different net waveform features. In this figure, the light blue traces are from individual trials (n=2), and the dark blue trace is the average across trials. The same dipole scaling factor (3000.0) is applied.

4.2 Viewing network spiking activity

From the Visualization tab, select “Drive-Spikes (2x1)” from the Layout template dropdown menu. From the Dataset dropdown, select the simulation data you would like to visualize. Next, click Make figure to generate a new figure.

Figure 9

This window shows the spiking activity produced in each population in response to the evoked inputs. The top two panels show histograms of distal evoked inputs (green) and proximal evoked inputs (red) provided to the neurons. The large third panel shows a raster plot of the spiking activity generated by the individual neurons, with different populations in different colors as labeled (x-axis: time in ms; y-axis: neuron identifier). The neuron identifiers are arranged vertically by layer, with top representing supragranular layers and the bottom representing the infragranular layers. Individual neuron types are drawn in the different colors shown in the legend. The dotted lines in the bottom panel show a time-series of summed activity per population (these use the same color code as the individual spikes; you can turn these lines off or on by selecting View> Toggle Histograms). The initial view shows the aggregate spiking activity across trials. To see spiking activity generated by a single trial, select the trial number using the combination box at the bottom of the window. This spike viewer window also provides the standard save/navigation functionality through the matplotlib control at the top.

4.3 Viewing ERP Spectrograms

HNN also allows you to generate a time-frequency representation of ERP dipole signals.

Figure 10 spectrogram

5. Comparing model output and recorded data

When loading data, the root-mean-squared-error (RMSE) is calculate between the average simulated ERP waveform and the loaded waveform. Our ERPyes-3trial simulation yields an RMSE of 5.51.

Let’s run a simulation with 10 trials and see how the RMSE changes.

Notice that when we simulated more trials, the RMSE between the data and the simulated average ERP changed to 5.59. Depending on the number of trials you run and what adjustments you make to the parameter values, you may be able to reduce the RMSE.

Figure 11

6. Adjusting parameters

Parameter adjustments will be key to developing and testing hypotheses on the circuit origin of your own ERP data. HNN is designed so that many of parameters in the model can be adjusted from the GUI (see Tour of the GUI). Here, we’ll walk through examples of how to adjust several “External drives” parameters to investigate how they impact the evoked response.

6.1 Changing the synchrony of the evoked inputs

In the example evoked response simulation described in Step 3, the time that each exogenous driving spike was provided to each cell in the local network was chosen independently from a Gaussian distribution, for both the proximal and distal inputs; this variability in timing created variability in the response of each cell in the network. Asynchronous exogenous drive is the default configuration in HNN. HNN also provides the capability to provide exogenous driving inputs to all cells in the network synchronously (0ms lag); this will reduce variability in timing of evoked inputs, producing a stronger response, and may or may not provide a better fit to the data.

To change the evoked inputs to contact the cells in network synchronously, first change the simulation name something like “ERPYes2TrialsSync” in Simulation. This will save the simulation data to a file with this new name.

Then, to provide synchronous exogenous drives to all cells in the network, go to External drives and set the Std dev time for all drives to 0.

Figure 12

After the simulation has completed, you’ll see the following output. Although the model replicates some gross features of the experimental data, the fit to the data is now substantially worse (RMSE = 14.15). Notice also that there is significantly lower variability of the input times in the histograms at the top of the figure (compare to evoked response inputs shown in Step 5), predicting (in this case) the evoked responses are more likely to be non-synchronous. Remember, however, that this simulation is only based on two trials.

Figure 13

6.1.1 Exercises for further exploration

Try adjusting the following. How do they affect the simulation? - Standard deviation and AMPA weights in External drives - Data Diople Scaling in Visualization - Dipole Smooth Window in Visualization

6.2 Simulating the Suprathreshold response

For this part of the tutorial, we’ll load a different experimental data set into the GUI. The new experimental data represents the evoked response from supra-threshold level stimuli in the experiment. We will adjust the parameters to better fit our simulation to the supra-detected threshold level data.

First, let’s load the data from the supra-threshold detected trials S1_SupraT.txt and the data from the threshold-level detected trials yes_trial_S1_ERP_all_avg.

Figure 18

Notice that the two waveforms displayed have very different features, including altered timing, amplitude, and sharpness of the peaks in this new data set.

If we also load in the previous simulation of threshold-level detected trials, ERP_yes_3trials, we can see that HNN automatically calculate the model fit to both datasets loaded (i.e, the RMSE between the average model response and the loaded data) based on the threshold-level simulation, showing that ERPYes100Trials.json gives a poor fit to the supra-threshold data (RMSE = 30.31).

Figure 19

What might account for these differences between our simulation and the data for the supra-threshold detected trials? We can formulate some hypotheses and use HNN to test it.

Given the timing and magnitude differences we see in the plot, we can hypothesize that these difference can be accounted by 1. An increase in the strength of the inputs that created the evoked response, and 2. A reduction in the delay in the arrival time of the inputs to the network.

To test these hypotheses, we’ll adjust the parameters as described. For simplicity, we have also created a parameter file that accounts for such changes and accurately reproduces the new data. Load the ERPYesSuparT.json file into External Drives.

You should see the new values in the dropdowns.

Next, let’s run the simulation using the Run button.

Figure 20

Let’s take a look at the output. First, note that the previous threshold-level simulation was not removed to allow comparison to the new supra-threshold simulation. For the new suprathreshold detection parameter set, there is now a better fit to the experimental data (RMSE = 24.39), with all major waveform features in the experiment and simulation in agreement. One interpretation of these results is that on supra-threshold detected trials, both the feedforward and feedback inputs are stronger, and both early and late feedforward inputs arrive earlier. These trials have overall stronger inputs, which then promotes tactile detection.

6.2.2 Exercises for further exploration

• Run 5 trials using the parameters from above. How does that impact the RMSE?
• Try manually adjusting the parameters to further reduce the RMSE.
• Try running the simulation for longer and adding an additional proximal or distal drive; how does this affect the simulation?

7. Have fun exploring your own data!

Follow steps 1-6 above using your data and parameter adjustments based on your own hypotheses.

8. Additional documents

8.1 Parameter files All parameter files for this tutorial can be found by clicking here or by visiting the following URL: https://github.com/jonescompneurolab/hnn/tree/master/param.

8.2 Video walkthrough A video walkthrough of an abridged version of the tutorial can be found by clicking here.

References