VERTEX
The Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX) is a Matlab tool for simulating extracellular potential recordings in spiking neural network (SNN) models. It uses a forward modelling approach to calculate extracellular potentials in a model given the position of the neurons relative to the virtual electrodes. VERTEX makes use of established theory on extracellular potential generation, combined with modern simulation methods and developments in simplified neuron modelling to simulate local field potentials (LFPs) from large neuronal network models encompassing more than 100,000 neurons [1]. The effect of applied electric fields can also be incorporated into the ongoing neuronal dynamics, allowing VERTEX to simulate the effect of electric field stimulation [2]. VERTEX is therefore best suited for simulations that seek to model a particular experimental output using realistic tissue geometry and neuron density. Example simulations may be:
- Modelling the LFP generated by pharmacologically induced neuronal oscillations in an in vitro brain slice preparation.
- Modelling the field potential generated by electrical stimulation in cortical grey matter.
- Modelling the effect of non-invasive brain stimulation techniques such as transcranial direct current stimulation.
Contents |
Overview
VERTEX is specifically set up to simulate layered structures - each neuron type specified by the user has a layer in which its soma will be placed. Synaptic connections are specified between neuron groups, each connection has a synapse type, which may indicate whether the connection should have a form of plasticity. Currently short term plasticity models and a spike timing dependent plasticity model are available. VERTEX saves the results of a simulation to file, these can then be loaded into MATLAB and visualised using the functions provided or analysed using home made functions written by users themselves.
What comprises a VERTEX simulation
To run a VERTEX simulation one must first specify parameters to build the network over which the simulation will run. These parameters can be split into the following groups: neuron parameters, connection parameters, and tissue parameters.
Tissue Parameters
The tissue parameters group includes the dimensions of the tissue we are simulating, as well as the layer boundaries (specifying the number and size of each layer in the structure), neuron density, and conductivity. The tissue parameters also includes any stimulation applied to the tissue, so here we would include the electric field and stimulation timings that we wish to apply. An example of the specification of tissue parameters is shown below. This specifies a single layered slice of tissue that is 2500 x 400 x 200 microns in volume, with a neuron density of 20000 and conductivity of 0.3 Sm^-1.
TissueParams.X = 2500; TissueParams.Y = 400; TissueParams.Z = 200; TissueParams.neuronDensity = 20000; TissueParams.numLayers = 1; TissueParams.layerBoundaryArr = [200, 0]; TissueParams.numStrips = 10; TissueParams.tissueConductivity = 0.3; TissueParams.maxZOverlap = [-1 , -1];
Neuron Parameters
The neuron parameters group includes the details required to simulate the intrinsic dynamics of each cell in the network. Each neuron group has its own set of neuron parameters. These include the proportion of the entire network that this group comprises, the layer that the group is in, the details of the dendrite compartments (how many there are and their dimensions), and the passive properties of the neuron (the membrane resistance, the axial resistance, the capacitance, and the leak potential). The model used to simulate the spiking mechanism of the group is also included here, the default neuron model used in VERTEX is the Adaptive exponential integrate-and-fire model (AdEx), and so this would also require the parameters that this model uses. The following illustrates how one would construct a neuron model in VERTEX, this will create an AdEx neuron group with 8 compartments as the first neuron group in the simulation. Additional neuron groups can be added to the NeuronParams structure as subsequent entries into the array.
NeuronParams(1).neuronModel = 'adex';
NeuronParams(1).V_t = -50;
NeuronParams(1).delta_t = 2;
NeuronParams(1).a = 2.6;
NeuronParams(1).tau_w = 65;
NeuronParams(1).b = 220;
NeuronParams(1).v_reset = -60;
NeuronParams(1).v_cutoff = -45;
NeuronParams(1).numCompartments = 8;
NeuronParams(1).compartmentParentArr = [0, 1, 2, 2, 4, 1, 6, 6];
NeuronParams(1).compartmentLengthArr = [13 48 124 145 137 40 143 143];
NeuronParams(1).compartmentDiameterArr = ...
[29.8, 3.75, 1.91, 2.81, 2.69, 2.62, 1.69, 1.69];
NeuronParams(1).compartmentXPositionMat = ...
[ 0, 0;
0, 0;
0, 124;
0, 0;
0, 0;
0, 0;
0, -139;
0, 139];
NeuronParams(1).compartmentYPositionMat = ...
[ 0, 0;
0, 0;
0, 0;
0, 0;
0, 0;
0, 0;
0, 0;
0, 0];
NeuronParams(1).compartmentZPositionMat = ...
[ -13, 0;
0, 48;
48, 48;
48, 193;
193, 330;
-13, -53;
-53, -139;
-53, -139];
NeuronParams(1).axisAligned = 'z';
NeuronParams(1).C = 1.0*2.96;
NeuronParams(1).R_M = 20000/2.96;
NeuronParams(1).R_A = 150;
NeuronParams(1).E_leak = -70;
NeuronParams(1).somaID = 1;
NeuronParams(1).basalID = [6, 7, 8];
NeuronParams(1).apicalID = [2 3 4 5];
Connection parameters
Using the tissue and neuron parameters VERTEX can create the neuron models required and position each cell randomly within its layer. The connection parameters include all the information required to connect our network up, and to simulate the synaptic dynamics. This includes the number of connections between neuron groups, the neuron compartments onto which the synapses form, and the spatial model used to generate the connections (the default in VERTEX is a Gaussian model, requiring a radius and limit). This allows VERTEX to wire up the network, it uses MATLAB’s built in random number generation and so one can ensure that the network is reproducible by specifying the random seed to use before generating the network. The synaptic model used between each pair of neuron groups is the final connection parameter. This includes the base synapse model (current based, conductance based, exponential, alpha), as well as any plasticity (short term plasticity, spike timing dependent plasticity) that may be applied. Synaptic parameters, such as the weight or decay constants, can be specified as a single value for all synapses in a particular connection or as a distribution, where each synaptic parameter is selected from the the distribution. The following code creates connections from neuron group 1 to neuron group 1 and 2, with a current based exponential synapse, and specifies the particular compartment groups which the synapses will target. Common parameters for all connections from neuron group 1 are also specified, these include the axon arbor spatial model and its parameters (the radius and limit), as well as the conduction speed and release delay.
ConnectionParams(1).numConnectionsToAllFromOne{1} = 1700;
ConnectionParams(1).synapseType{1} = 'i_exp';
ConnectionParams(1).targetCompartments{1} = {'basalID', ...
'apicalID'};
ConnectionParams(1).weights{1} = 1;
ConnectionParams(1).tau{1} = 2;
ConnectionParams(1).numConnectionsToAllFromOne{2} = 300;
ConnectionParams(1).synapseType{2} = 'i_exp';
ConnectionParams(1).targetCompartments{2} = {'dendritesID'};
ConnectionParams(1).weights{2} = 28;
ConnectionParams(1).tau{2} = 1;
ConnectionParams(1).axonArborSpatialModel = 'gaussian'; ConnectionParams(1).sliceSynapses = true; ConnectionParams(1).axonArborRadius = 250; ConnectionParams(1).axonArborLimit = 500; ConnectionParams(1).axonConductionSpeed = 0.3; ConnectionParams(1).synapseReleaseDelay = 0.5;
The simulation process
Having built the network one can then look to simulating its dynamics. This will require the specification of some simulation settings (e.g. the duration of the simulation, the time step), and recording settings (e.g. which variables to record and from which neurons to record them, sampling rate, etc). The simulation loop (the code executed at each time step) involves updating any externally applied electric field, updating all neuron and synapse variables, recording the value of all variables we wish to record, and processing all spikes (this will also involve updates to synaptic variables when plasticity is used). Shown below is an example of how to set up the recording settings structure to record the LFP using a grid of electrodes, as well as recording the membrane potential of every tenth of the first 2000 neurons. The sample rate is set at 5 kHz and the maximum recording time is set at 500 ms indicating after 500 ms of recording a new file will be created. The simulation settings are set up to simulate for 500 ms with a time step of 0.03125 ms, and the simulation will be run in parallel.
RecordingSettings.saveDir = '~/VERTEX_results/'; RecordingSettings.LFP = true; [meaX, meaY, meaZ] = meshgrid(0:1000:2000, 200, 600:-300:0); RecordingSettings.meaXpositions = meaX; RecordingSettings.meaYpositions = meaY; RecordingSettings.meaZpositions = meaZ; RecordingSettings.minDistToElectrodeTip = 20; RecordingSettings.v_m = 1:10:2000; RecordingSettings.maxRecTime = 500; RecordingSettings.sampleRate = 5000;
SimulationSettings.simulationTime = 500; SimulationSettings.timeStep = 0.03125; SimulationSettings.parallelSim = true;
Parallelisation
VERTEX makes use of MATLAB's single program multiple data (SPMD) parallel programming construct to take advantage of multiple cores. This requires the MATLAB parallel computing toolbox and involves assigning each neuron to a worker on which its neural and synaptic dynamics will be computed. Spikes are transmitted between workers during each simulation step.
Neuron and synapse models
Neuron and synapse types are described using inheritance to avoid the duplication of functionality. The abstract NeuronModel class describes the functionality provided by all multi-compartment neurons. It contains the membrane potential, external potential, and axial current (the currents that flow between compartments as a result of the difference between their membrane potentials) properties, as well as the functionality required to integrate these. The integration of the stimulating field is included as an additional step during the calculation of the axial currents and is performed at each time step when the stimulation is turned on. It is part of the core functionality of the abstract Neuron class. Classes with specific mechanisms then inherit from this, e.g. the NeuronModel_passive class provides a simple wrapper on top to allow a neuron with no active channels. The NeuronModel_adex adds the adaptive exponential integrate and fire mechanism to the soma, allowing the cell to generate action potentials. Here each instance of a class would represent a group of neurons in the same layer and of the same type. This allows us to ultilise MATLAB’s vectorised operations when updating variables so that for example: the membrane potential variable (v_m) holds the membrane potentials of all neurons in this group as a matrix. This also allows us to utilise the object oriented design advantages without the overhead that would come from storing each neuron or synapse as its own object. The integration of the axial current involves a loop over all possible neighbouring compartments with an operation vectorised for each compartment. The class hierarchy relevant to conductance based exponential synapses (SynapseModel_g_exp) is shown in Figure 2. Here, multiple inheritance is used to allow many combinations of synapse types to be efficiently described. Synapse models have a base synapse type (defining how the synapse operates without plasticity, e.g. g_exp will be a conductance based exponential synapse), it can then also have short term plasticity (ab for the Abbott model or mt for the Markram and Tsodyks), spike timing dependent plasticity (stdp), or both. The plasticity models are defined as separate classes from which the synapse model can inherit from.
| Synapse Class | Description |
|---|---|
| SynapseModel_i_exp | Current based, exponential synapse |
| SynapseModel_g_exp | Conductance based, exponential synapse |
| SynapseModel_g_exp_mt | Conductance based, exponential synapse with Markram and Tsodyks[4] short term plasticity |
| SynapseModel_g_exp_ab | Conductance based, exponential synapse with Abbott model[3] of short term plasticity |
| SynapseModel_g_exp_stdp | Conductance based, exponential synapse with STDP |
| SynapseModel_g_exp_mt_stdp | Conductance based, exponential synapse with STDP and Markram and Tsodyks[4] short term plasticity |
| SynapseModel_g_exp_ab_stdp | Conductance based, exponential synapse with STDP and Abbott model[3] of short term plasticity |
Electric field stimulation
Incorporating the effect of an electric field generated by a stimulating electrode is relatively straight forward in VERTEX. The first step is to model the electric field, the MATLAB Partial Differential Equation (PDE) toolbox can be used to model most simple electrode setups, and the output of this model attached to the tissue parameter structure, forming part of the input to the VERTEX simulation. The only further information needed is the times at which to apply the stimulation, which are specified as a list of times at which to turn the stimulation on and off. The following code will set the stimulation field and the times at which to apply it. To generate an appropriate field (either as a stationary or time dependent field) follow the guide on the MATLAB PDE Toolbox website. Time dependent fields that are shorter than the duration over which they are applied will be looped, this allows oscillating fields to be generated only for a single period, saving memory.
TissueParams.StimulationField = stimfield; % stimfield is a either a StationaryResults object or TimeDependentResults object from the PDE toolbox. TissueParams.StimulationOn = [600 700];% pulse interval of 100 ms TissueParams.StimulationOff = [600.5 700.5];% pulse width of 0.5 ms
Analysing simulation output
Having run a VERTEX simulation one can then load the simulation results and perform analysis, either using custom code or using functions provided with VERTEX. VERTEX provides functions to plot a spike raster of the activity, or to plot the LFP recorded at a particular electrode. It also provides functions to visualise the weights of synapses in the network, either on a group to group basis, or for each individual connection. The following code will load the results from the directory they were saved in, and plot a spike raster of the activity.
Results = loadResults(RecordingSettings.saveDir); rasterParams.groupBoundaryLines = [0.7, 0.7, 0.7]; rasterParams.title = 'Spike Raster'; rasterParams.xlabel = 'Time (ms)'; rasterParams.ylabel = 'Neuron ID'; rasterFigure = plotSpikeRaster(Results, rasterParams);
The following code will plot the stimulation field stored within the TissueParams struct, and then plot within this the soma positions of the cells in the network.
pdeplot3D(TissueParams.StimulationModel, 'ColorMap', TissueParams.StimulationField.NodalSolution, 'FaceAlpha', 0.8) hold on; pars.colors = rasterParams.colors; plotSomaPositions(Results.params.TissueParams,pars);
Performance
The VERTEX simulator can run simulations at the scale of a cortical column on a modern desktop computer, and at the scale of an in vitro brain slice preparation on a high performance computing node. Small scale simulations Parallelisation significantly improves the performance of large VERTEX simulations. Figure 4 shows the simulation times for a small and large network with increasing numbers of processing cores available, increasing the number of electrodes from which to record an LFP also increases the simulation time.
Comparison with similar software tools
Several simulation packages offer similar capabilities to VERTEX. The Brian simulator, a python based spiking neural network simulator, has good support for synaptic plasticity, both short term plasticity models and spike timing dependent plasticity models. As it focuses on simulating networks of point neurons, it has limited support for including electric field stimulation and LFP calculations using the line source model.
The Neuron simulation environment [5] supports electrical stimulation through its extracellular mechanism, and LFP simulation through the LFPy tool [6]. It also supports a wide variety of synapse types including short term plasticity and spike timing dependent plasticity. However, the Neuron simulator is not considered easy to use as it is designed for simulating single neurons or small networks in detail and it can become cumbersome for those wishing to simulate larger multi-layered networks. VERTEX is specifically designed for modelling multi-layered structures such as the neocortex, and it provides built in support for modelling LFP generation, synaptic plasticity, and electric field stimulation.
The NEST (NEural Simulation Tool) [7] is a mature and easy to use spiking neural network simulator, that allows users to simulate the activity of large scale laminar structures with various forms of synaptic plasticity and a range of neuron and synapse models. It can also be combined with the LFPy tool [8] to simulate the local field potential generated by a large network.
References
1. Tomsett, R. J., Ainsworth, M., Thiele, A., Sanayei, M., Chen, X., Gieselmann, M. A., … Kaiser, M. (2015). Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX): comparing multi-electrode recordings from simulated and biological mammalian cortical tissue. Brain Struct Funct, 220(4), 2333–2353. https://doi.org/10.1007/s00429-014-0793-x
2. Thornton C, Hutchings F and Kaiser M. The Virtual Electrode Recording Tool for EXtracellular Potentials (VERTEX) Version 2.0: Modelling in vitro electrical stimulation of brain tissue [version 1; peer review: 2 approved]. Wellcome Open Res 2019, 4:20 (https://doi.org/10.12688/wellcomeopenres.15058.1)
3. Abbott, L. F., Varela, J. A., Sen, K., & Nelson, S. B. (1997). Synaptic depression and cortical gain control. Science, 275(5297), 221–224.
4. Tsodyks, M., Pawelzik, K., & Markram, H. (1998). Neural networks with dynamic synapses. Neural Computation, 10(4), 821–835. https://doi.org/10.1162/089976698300017502
5. Hines ML, Carnevale NT: The NEURON simulation environment. Neural Comput. 1997; 9(6): 1179–209.
6. Lindén Henrik, Hagen Espen, Leski Szymon, Norheim Eivind, Pettersen Klas, Einevoll Gaute (2014) LFPy: a tool for biophysical simulation of extracellular potentials generated by detailed model neurons. Frontiers in Neuroinformatics , Volume 7
7. Diesmann M and Gewaltig M-O. (2002) NEST: An Environment for Neural Systems Simulations. Forschung und wisschenschaftliches Rechnen, Beiträge zum Heinz-Billing-Preis 2001. Ges. für Wiss. Datenverarbeitung.
8. Hagen, E., Stavrinou, M. L., Linden, H., Tetzlaff, T., van Albada, S., Dahmen, D., … Einevoll, G. T. (2013). Hybrid scheme for modeling LFPs from spiking cortical network models. BMC Neuroscience, 14(Suppl 1), P119. doi:10.1186/1471-2202-14-S1-P119
Acknowledgements
Development of this software was supported by BBSRC (BB/F016980/1), EPSRC (EP/E002331/1, EP/G03950X/1, EP/K026992/1, NS/A000026/1), and the Wellcome Trust (102037).
Internal links
Ted Carnevale (2007) Neuron simulation environment. Scholarpedia, 2(6):1378.
Marc-Oliver Gewaltig and Markus Diesmann (2007) NEST (NEural Simulation Tool). Scholarpedia, 2(4):1430.
Dan F. M. Goodman and Romain Brette (2013) Brian simulator. Scholarpedia, 8(1):10883.
Alain Destexhe and Claude Bedard (2013) Local field potential. Scholarpedia, 8(8):10713.
Rob Schreiber (2007) MATLAB. Scholarpedia, 2(7):2929.
Frances K. Skinner (2006) Conductance-based models. Scholarpedia, 1(11):1408.
External links
VERTEX website [[1]]
VERTEX introduction video [[2]]
Dynamic Connectome website [[3]]
See also
Genesis, NEURON Simulation Environment, XPPAUT, NEST (NEural Simulation Tool)
