Tutorial 1: VERTEX basics


In this tutorial, we will introduce the core methods for setting up and running a VERTEX simulation, and recording LFPs from the network. VERTEX is designed to facilitate the simulation of extracellular potentials generated by activity in spiking neural networks; in particular, spatially-organised networks containing thousands or hundreds of thousands of neurons. VERTEX’s interface and model specification options were designed with this particular task in mind. It is therefore less flexible than other neural simulators (e.g. NEURON, NEST, Brian, GENESIS, Moose), but the limited scope has allowed us to simplify the user interface so that a simulation can be specified simply by setting some parameters and run using a few function calls. The parameters are divided into five categories: neuron group properties, connectivity & synapses, tissue properties, recording settings, and simulation settings. The parameters associated with each category are specified in Matlab structures or structure arrays.

In this tutorial, we will walk through the different parameter settings for creating a single group of neurons that synapses with itself.

Before we start though, we’ll quickly go through VERTEX’s built-in help features.

If you need a list of all the functions available to you in VERTEX, you can run the vertexFunctions() function, which prints all available functions with a brief description of what they do:

vertexFunctions
  initNetwork:  Initialise simulation environment, setup the network and calculate constants for extracellular potential simulation.
  loadResults:  loads the results of a simulation run.
  neuronDynamics:  runs a simulation of a single neuron group.
  resetRandomSeed:  Resets the random number generator seed.
  runSimulation:  Run the simulation given the model generated by initNetwork().
  getGroupConnectivity:  Get the number of connections between neuron groups.
  getSparseConnectivity:  Gets the sparse connectivity matrix for the network.
  groupRates:  Calculate the average firing rate of the neurons in each group.
  neuronRates:  calculates the mean firing rate of each neuron.
  plotSomaPositions:  Plots the positions of the neurons' soma compartments.
  plotSpikeRaster:  Creates a spike raster plot.
  viewMorphology:  Plots a neuron's morphology.

To get more in-depth help on a function, type help and then the function name just like a built-in Matlab function. For example, to see details of for the runSimulation() function, we enter:

help runSimulation
 RUNSIMULATION Run the simulation given the model generated by initNetwork().
    RUNSIMULATION(PARAMS, CONNECTIONS, ELECTRODES) runs the simulation
    given the model generated by initNetwork(). PARAMS, CONNECTIONS and
    ELECTRODES are the PARAMS, CONNECTIONS and ELECTRODES outputs from the
    initNetwork() function. RUNSIMULATION automatically saves the simulation
    results in the directory specified by the user in the recording
    settings structure given to initNetwork().

If you need to check the physical units used in VERTEX, run the vertexUnits() function, which prints out a list of physical quantities with their associated units and a representative example:

vertexUnits
  Units used by VERTEX: 
 
       Electric potential:            milliVolts              (e.g. leak potential E_leak)
              Conductance:           nanoSiemens              (e.g. synaptic conductance weight)
                  Current:             picoAmps               (e.g. synaptic current weight)
      Specific resistance:     Ohm * square centimetres       (e.g. specific membrane resistance R_M)
     Specific capacitance:  micfroFarads / square centimetre  (e.g. specific membrane capacitance C)
  Longitudinal resistance:         Ohm * centimetres          (e.g. intracellular axial resistance R_A)
             Conductivity:          Siemens / metre           (e.g. extracellular conductivity sigma)
                     Time:            milliseconds            (e.g. synaptic time constant tau)
                   Length:            micrometres             (e.g. model width X)
 

If you’re really lost and forget the names of these functions, you can remind yourself by simply running vertexSimulator:

vertexSimulator
 
    Virtual Electrode Recording Tool for EXtracellular potentials
        For a list of functions, enter vertexFunctions
        To display the physical units used by VERTEX, enter vertexUnits
 

Now that we’ve covered the built-in help, we can move on to setting up and running a simulation.

Contents

Tissue Parameters

We start with the structure array holding the overall information about the piece of brain tissue we are modelling, which we call TissueParams. The order in which parameters are specified does not matter to VERTEX, so we proceed in what seems to us like a reasonable order, starting with the size of the model.

TissueParams.X = 2500;
TissueParams.Y = 400;
TissueParams.Z = 200;
TissueParams.neuronDensity = 25000;

X, Y and Z specify the dimensions of the model, in micrometres. The model space in this case is cubic, as in the brain slice model described in (Tomsett et al. 2014). Alternatively, if we don’t specify X or Y, but instead specify R, then VERTEX will create a cylindrical model with radius R micrometres and depth Z micrometres. neuronDensity gives the overall density of neurons in the model, in neurons per cubic mm. The number of neurons in the model is then calculated automatically (in this case there will be 2.5 * 0.4 * 0.2 * 25 000 = 5 000 neurons in the model).

Next we specify some more parameters concerning the spatial layout of the model:

TissueParams.numLayers = 1;
TissueParams.layerBoundaryArr = [200, 0];
TissueParams.numStrips = 10;

numLayers is the total number of layers we want to create in the model. Layers can be used to organise both the spatial organisation of the model, and the connectivity of and dynamics of neuron groups in different layers; for example, the model described in (Tomsett et al. 2014) uses five layers to represent neocortical layers 1, 2/3, 4, 5 and 6. In this tutorial, we just want to have one layer. We assume the layer boundaries are defined as x-y planes, so we set their boundaries by specifying their z-depths. layerBoundaryArr contains a list of these depths, going from the top of the highest layer to the bottom of the lowest one. The length of layerBoundaryArr should therefore be equal to numLayers + 1.

Finally, we introduce the concept of ‘strips’. When setting up the neuron positions in a cuboid model, VERTEX will position each neuron’s soma randomly within the boundaries of its containing layer. This often does not matter, but sometimes it is useful to know more precisely where the neurons are positioned – when plotting spike trains to reveal spatial variation in firing patterns, for example. A cuboid model can be divided into two or more spatial ‘strips’, which are ordered from left to right on the x-axis. Neurons are then positioned by strip so that the lowest neuron IDs are in the the furthest left strip, and the highest neuron IDs in the the furthest right strip. The width of each strip is the width of the model, X, divided by the number of strips. Changing numStrips does not affect the model’s dynamics as such; it affects the precision of the placement of neurons along the x-axis.

TissueParams.tissueConductivity = 0.3;
TissueParams.maxZOverlap = [-1 , -1];

tissueConductivity sets the extracellular medium’s conductivity, in Siemens per metre. This parameter is only required when using simulated extracellular electrodes. maxZOverlap sets the maximum distance outside the model’s upper and lower z-boundaries that a neuron’s dendrites can extend. As neurons are positioned according to their soma compartment, if the soma is placed close to the model’s boundaries, then its dendritic compartments will extend outside these boundaries. In the present model, we set the model z-depth to be 200 microns, so pyramidal neuron apical dendrites are bound to extend above the model space (we will be using layer 2/3 pyramidal cell models with an apical dendrite 330 micrometres long). In the current case this is desirable, as we set the 200 micron depth to be the depth of the soma-layer rather than meaning it to be the depth of the whole layer. Therefore we set maxZOverlap to -1, meaning that we do not want to set a maximum z-overlap length. In the cortical slice model described in (Tomsett et al. 2014), though, the layer boundaries are meant to represent the entire depth of the layer and the full size of the model. In that case, we set maxZOverlap to [0, 100]. This prevented the uppermost pyramidal neuron dendrite compartments from extending above the top of layer 1, but allowed dendritic compartments at the bottom of layer 6 to extend below the model by up to 100 microns into the “white matter”.

Neuron Group Parameters

We have finished specifying the TissueParams structure, so now we specify the neuron group information in a structure array we call NeuronParams. We store the parameters for our only neuron group in NeuronParams(1), in case we want to add further groups at a later date. This neuron group is henceforth referred to as neuron group 1.

NeuronParams(1).modelProportion = 1;
NeuronParams(1).somaLayer = 1;
NeuronParams(1).neuronModel = 'passive';

The first line defines the proportion of the total model size made up by neuron group 1. In this case we are only modelling a single neuron group, so the whole model is made up of neurons in this group. The second line specifies which layer of the model this neuron group is in (this must be layer 1, as we only have 1 layer in the model). The third line specifies what model to use to calculate the dynamics of the neurons in group 1. We have set this to 'passive' – so the neurons will only exhibit passive membrane dynamics. These kinds of groups can be useful for simulating the LFP given synaptic input from other neuron groups (see later tutorials), but will not generate any spikes on their own. Instead, we will specify group 1 to consist of neurons randomly firing with Poisson statistics, and a firing rate of 5 Hz:

NeuronParams(1).neuronModel = 'poisson';
NeuronParams(1).firingRate = 5;

Each neuron model in VERTEX has some unique parameters that must be specified when using them. The 'poisson' model’s only extra parameter is called firingRate, and sets each neuron in the group’s firing rate to be the specified number of spikes per second. Other models’ unique parameters can be found in the VERTEX reference.

Next we specify the number of compartments per neuron we want.

NeuronParams(1).numCompartments = 1;

If we specify 1 compartment, then VERTEX will treat this group as consisting of point neurons. No further parameters would be required in the NeuronParams structure array in this case. This is useful for specifying groups of “external” populations that provide external input to the simulation, but do not contribute to the LFP, for example. However, as we want to simulate the LFP in this group, we need to use compartmental neurons (Pettersen et al. 2012). We will use the layer 2/3 pyramidal neuron models described in (Tomsett et al. 2014), which contain 8 compartments.

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';

compartmentParentArr gives the compartmental structure. We assume the neurons’ compartmental structure will always be a tree, so each compartment will have a maximum of one parent compartment – and the soma has no parent. compartmentParentArr lists the parent compartment ID at the index of each child compartment; compartment ID 1 is the soma so has no parent compartment and is set to 0. compartmentDiameterArr gives the diameter of each compartment in microns, again with the indices of the array representing the compartment IDs. The compartment[XYZ]PositionMat values specify the start and end coordinates of each compartment, presuming the soma compartment is positioned at (0,0,0). During network construction, these values are used to calculate the position of each compartment of each neuron after random positioning and rotation. They are also used to calculate the length of each compartment (though a compartmentLengthArr list containing the length of each compartment as defined by the user can be specified to override the automatic calculations). Row indices represent compartment IDs, column 1 holds compartment start coordinates and column 2 holds column end coordinates.

Next we specify the neurons’ passive properties:

NeuronParams(1).C = 1.0 * 2.96;
NeuronParams(1).R_M = 20000 / 2.96;
NeuronParams(1).R_A = 150;
NeuronParams(1).E_leak = -70;

C is the specific membrane capacitance in microFarads per square centimetre, R_M is the specific membrane resistance in Ohms per square centimetre, R_A is the axial resistance in Ohm-centimetres and E_leak is the membrane leak conductance reversal potential in mV.

The next code listing is entirely optional, but useful when specifying the model connectivity parameters. We specify labels for the different compartments in the model, so that we don’t have to remember explicitly the numbers of the compartments that comprise different sections of the neuron. The name of each parameter is chosen to be descriptive, but could be chosen arbitrarily.

NeuronParams(1).basalID = [6, 7, 8];
NeuronParams(1).apicalID = [2 3 4 5];

As we are only specifying one neuron group, we are now finished with the NeuronParams structure array (further groups would have parameters specified in NeuronParams(2), NeuronParams(3) etc.).

Connectivity Parameters

VERTEX treats spikes as discrete events that trigger some postsynaptic event at targeted neurons. Connectivity is specified probabilistically in terms of connections between groups of neurons. We specify the relevant parameters in a structure array called ConnectionParams.

VERTEX treats each structure in the structure array as referring to a presynaptic population of neurons – ConnectionParams(1) contains the connectivity parameters for neuron group one to make connections to other postsynaptic groups. Where parameters need to be specified on a per-postsynaptic-group basis, we use a cell array with each cell index referring to the ID of that postsynaptic neuron group.

ConnectionParams(1).numConnectionsToAllFromOne{1} = 2000;
ConnectionParams(1).synapseType{1} = 'i_exp';
ConnectionParams(1).tau{1} = 2;
ConnectionParams(1).weights{1} = 10;
ConnectionParams(1).targetCompartments{1} = ...
  [NeuronParams(1).basalID, NeuronParams(1).apicalID];

In the above code, we first say that each presynaptic neuron in group 1 (index of the structure array) makes 2000 connections to postsynaptic neurons in group 1 (index of the cell array). We specify that the type of these synapses is 'i_exp' – single exponential, current-based synapses (the other currently available synapse type in VERTEX is 'g_exp', which are conductance-based single exponential synapses; in this case, a reversal potential for the synapse must be specified in the relevant E_reversal parameter cell). We specify the weight of these synapses – for current-based synapses in pA, and for conductance-based in nS – and the exponential decay time constant tau in ms. Finally, targetCompartments lists which compartments on the postsynaptic neurons are allowed to be contacted by the presynaptic neurons (this is where the labels we gave to the compartment numbers comes in handy).

Next, we need to specify the axonal arbour properties.

ConnectionParams(1).axonArborSpatialModel = 'gaussian';
ConnectionParams(1).sliceSynapses = true;
ConnectionParams(1).axonArborRadius = 250;
ConnectionParams(1).axonArborLimit = 500;

The axon arbour of each neuron is defined as a 2D probability distribution in the x-y plane centred at the presynaptic neuron, and defines how the probability of making a connection to a postsynaptic neuron varies in space. Specifying the axonArborSpatialModel to be 'gaussian' means that the connection probability decays away from a presynaptic neuron as a 2D gaussian with standard deviation given by the axonArborRadius parameter. If we instead set axonArborSpatialModel to be 'uniform', the axonArborRadius parameter would represent the maximum extent of the arbor and the connection probability would be constant within this radius. These are the two currently implemented connectivity profiles; we intend to add more complex profiles (e.g. non-isotropic gaussians, patchy projections) in the future. axonArborLimit is an optional parameter for gaussian arbors, and specifies a cutoff point beyond which no connections can be made. If this parameter is not set, then neurons have no limit on their maximum connection distance.

Finally, we need to set the axonal conduction speed, and the constant synaptic release delay for synapses made by this group of neurons:

ConnectionParams(1).axonConductionSpeed = 0.3;
ConnectionParams(1).synapseReleaseDelay = 0.5;

axonConductionSpeed specifies the speed with which action potentials propagate down axons in this neuron group, in m/s. This speed is used to calculate axonal propagation delays after a neuron spikes. synapseReleaseDelay sets a constant time, in ms, for neurotransmitter release after an action potential reaches a synaptic terminal. This time is simply added to the distance-dependent axonal propagation delay.

Tip: if you want to specify a network with constant (not distance-dependent) axonal delay times, you could set synapseReleaseDelay to the desired constant delay time and set axonConductionSpeed = Inf.

Recording settings

In order to analyse the output of a model, we need to be able to record the variables over the time course of the simulation. We specify what we want to record in the recording settings structure. VERTEX will save the spike trains of all neurons by default, but can also record individual neuronal membrane potentials and LFPs.

RecordingSettings.saveDir = '~/VERTEX_results_tutorial_1/';
RecordingSettings.LFP = true;
[meaX, meaY, meaZ] = meshgrid(0:500:2500, 200, 600:-200:0);
RecordingSettings.meaXpositions = meaX;
RecordingSettings.meaYpositions = meaY;
RecordingSettings.meaZpositions = meaZ;
RecordingSettings.minDistToElectrodeTip = 20;

First we set the directory in which we want to save the files that VERTEX generates. Then we set LFP to true, so VERTEX knows we want to record LFPs. We then use Matlab’s meshgrid function to calculate x, y and z coordinates (in microns) of our multi-electrode array. VERTEX expects electrode coordinates in the format output by meshgrid , but it isn’t necessary to use the meshgrid function to specify the electrode positions – they can be placed arbitrarily. Model coordinates in VERTEX use standard x, y and z axis directions, so we specify the z-coordinates of the MEA in descending order so that they are numbered from top to bottom as in the experimental MEA in (Tomsett et al. 2014). This example produces a 2.5 mm by 0.6 mm MEA with 24 electrodes, at a constant y-coordinate of 0.2 mm, with an inter-electrode spacing of 0.5 mm in the x-plane and 0.2 mm in the z-plane. Finally we set minDistToElectrodeTip – this value ensures that no neuronal compartment can be positioned too close to an electrode tip and therefore unrealistically dominate the LFP. We choose 20 microns, which is the default value if this parameter is not set by the user.

We can also choose which neurons we want to record the membrane potential from, if any. We provide this in a list of neuron IDs. While the IDs are only generated during model initialisation, we know roughly where in the model space the neuron of each ID will be because of the method of using strips to constrain the x-coordinates of the neurons. We can sample membrane potentials of one neuron from each strip, for example, to look at differences in membrane potentials across the space of the model. We will specify:

RecordingSettings.v_m = 500:500:4500;

to tell VERTEX to record from neuron IDs 500, 1500, 2500, … 4500, ensuring that we record from one neuron in each strip. The membrane % potential is always recorded at the soma compartment.

Next we set the sample rate and the maximum simulation time between saves:

RecordingSettings.maxRecTime = 100;
RecordingSettings.sampleRate = 1000;

maxRecTime specifies, in milliseconds, the amount of time that VERTEX records for before saving the recordings to disk and starting recording again. This chunks the recordings into files of length maxRecTime * sampleRate / 1000, and ensures that for large MEAs the file sizes and memory usage do not grow too large. We have found that setting this to 100-200 ms provides a good compromise between memory usage and slowing the simulation down due to excessive file save operations, though the optimal value will depend on your model details. In smaller models with few electrodes, it may be suitable to set this to be equal to the total simulation length. sampleRate is simply the rate at which VERTEX samples the LFP, in Hz. High values can severely slow simulation time in large models, or for large numbers of electrodes. 1000 Hz provides a good resolution for most LFP investigations (note: in the next section we will describe how to set the simulation time step. If the simulation time step and specified sample rate are incompatible, VERTEX will automatically sample at the closest (higher) frequency that is compatible with the time step).

General simulation settings

Finally we need to specify some overall simulation settings.

SimulationSettings.simulationTime = 1000;
SimulationSettings.timeStep = 0.03125;
SimulationSettings.parallelSim = false;

We set the total simulation time to 1000 ms, and the global integration time step to 0.03125 ms. We then tell VERTEX to run in serial mode by setting parallelSim to false. If you want to run a parallel simulation, set parallelSim to true, and provide the Matlab pool size parallel profile to use:

SimulationSettings.poolSize = 4;
SimulationSettings.profileName = 'local';

If poolSize and profileName are not set, then VERTEX will ask Matlab to start a parallel job with Matlab’s default settings. If you do not have the Matlab Parallel Computing Toolbox installed, then VERTEX will just run the simulation in serial mode.

Generate the network

Now that we have set all the required parameters for our simulation in the relevant structures, we can generate the network.

[params, connections, electrodes] = ...
  initNetwork(TissueParams, NeuronParams, ConnectionParams, ...
              RecordingSettings, SimulationSettings);

The initNetwork function performs a number of operations. Firstly it sets up the parallel environment if a Matlab parallel pool is not already running. Secondly, it creates the neuron groups, positions them in space, and if necessary distributes them between parallel processes. Thirdly, it calculates the model connectivity based on the connection statistics and neuron positions. Finally, if recording the LFP, it pre-calculates the constants required for the line-source extracellular potential calculation based on the neuron and electrode positions.

Run the simulation

Now that we have generated the network, we can run the model to simulate the dynamics in the network and record the resultant LFPs.

runSimulation(params, connections, electrodes);

This function will run the dynamics in the simulation according to the specified parameters and previously generated network, storing the variables that have been specified to record from in the location set in RecordingSettings. The timestep will be printed every 5 milliseconds so you can keep track of how long the simulation has left to run. Once the simulation has run (this may take some time depending on your computer and the size of the model), you can load the results for analysis.

Results = loadResults(RecordingSettings.saveDir);

The loadResults function loads variables saved from a previous simulation in the input directory specified as its input argument. loadResults automatically reconstructs the data saved by runSimulation for both serial and parallel simulations. It returns a structure which holds the spike, LFP and membrane potential (if requested) recordings. Recordings.spikes contains all neurons’ spike times, listed so that each row contains [neuron ID, spike time in ms]. To plot all the spikes in a spike raster, we can do:

figure(1)
plot(Results.spikes(:, 2), Results.spikes(:, 1), 'k.')
axis([0 500 -50 5050])
set(gcf,'color','w');
set(gca,'YDir','reverse');
set(gca,'FontSize',16)
title('Tutorial 1: spike raster', 'FontSize', 16)
xlabel('Time (ms)', 'FontSize', 16)
ylabel('Neuron ID', 'FontSize', 16)

 

Tutorial 1 spike raster

Recordings.LFP is a matrix in which each row contains the LFP recording from the electrode of that row index (see above). So, to plot the LFP from electrode 3, for example, we can do:

figure(2)
plot(Results.LFP(3, :), 'LineWidth', 2)
set(gcf,'color','w');
set(gca,'FontSize',16)
title('Tutorial 1: LFP at electrode 3', 'FontSize', 16)
xlabel('Time (ms)', 'FontSize', 16)
ylabel('LFP (mV)', 'FontSize', 16)

tutorial_1_02

Finally, Recordings.v_m contains the soma membrane potential of each neuron we specified should be recorded from in RecordingSettings, ordered in ascending ID order. Again, each row is a time series for a neuron. Given that we set RecordingSettings.v_m = 500:500:4500 above, if we want to plot the membrane potential of neuron 1500 we would do:

figure(3)
plot(Results.v_m(3, :), 'LineWidth', 2)
set(gcf,'color','w');
set(gca,'FontSize',16)
title('Tutorial 1: membrane potential for neuron 1500', 'FontSize', 16)
xlabel('Time (ms)', 'FontSize', 16)
ylabel('Membrane potential (mV)', 'FontSize', 16)

tutorial_1_03

If all has gone according to plan, you should have set up a network model containing 5000 layer 2/3 pyramidal neurons, spiking with Poisson spike trains at 5 Hz, simulated the dynamics of this model, recorded LFPs, spike trains and membrane potentials from the simulation, loaded these recordings, and plotted them. In the next tutorial, we will look at a network containing two neuron groups, and give the neurons proper spiking dynamics to generate network activity.

If you have experienced any problems when trying to run this tutorial, or if you have any suggestions for improvements, please contact us using the contact form.

References

Tomsett RJ, Ainsworth M, Thiele A, Sanayei M, Chen X et al. (2014) Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX): comparing multi-electrode recordings from simulated and biological mammalian cortical tissue, Brain Structure and Function. doi:10.1007/s00429-014-0793-x

Pettersen KH, Lindén H, Dale AM, Einevoll GT (2012) Extracellular spikes and current-source density. In: Brette R, Destexhe A (eds) Handbook of Neural Activity Measurement. Cambridge University Press