Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX): comparing multi-electrode recordings from simulated and biological mammalian cortical tissue

Local field potentials (LFPs) sampled with extracellular electrodes are frequently used as a measure of population neuronal activity. However, relating such measurements to underlying neuronal behaviour and connectivity is non-trivial. To help study this link, we developed the Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX). We first identified a reduced neuron model that retained the spatial and frequency filtering characteristics of extracellular potentials from neocortical neurons. We then developed VERTEX as an easy-to-use Matlab tool for simulating LFPs from large populations (>100,000 neurons). A VERTEX-based simulation successfully reproduced features of the LFPs from an in vitro multi-electrode array recording of macaque neocortical tissue. Our model, with virtual electrodes placed anywhere in 3D, allows direct comparisons with the in vitro recording setup. We envisage that VERTEX will stimulate experimentalists, clinicians, and computational neuroscientists to use models to understand the mechanisms underlying measured brain dynamics in health and disease. Electronic supplementary material The online version of this article (doi:10.1007/s00429-014-0793-x) contains supplementary material, which is available to authorized users.


Introduction
Many measurement techniques have been used to study neuronal dynamics, including optical imaging methods (voltage-sensitive dye imaging, calcium imaging, intrinsic signal optical imaging), intracellular electrode recordings of individual neurons, and extracellular recordings using single or multiple electrodes (Brette and Destexhe 2012).
While each modality provides some information about the system's dynamics, it is not always clear how this information is related to the underlying neuronal activity. Intracellular recordings are easiest to interpret because of the strong theoretical foundations of cellular neurophysiology that have arisen over many decades (Johnston and Wu 1995), but the theory linking measurements made by many other methods to neuronal activity are lacking. This deficit in theory, combined with the increasing use of different recording techniques to sample from ever larger M. O. Cunningham and M. Kaiser are joint senior authors.
Electronic supplementary material The online version of this article (doi:10.1007/s00429-014-0793-x) contains supplementary material, which is available to authorized users. neuron populations, has stimulated the idea of ''modelling what you can measure''  to help fill these theoretical gaps.
We aim to contribute to this effort by modelling the measurements made by multi-electrode arrays (MEAs). MEAs record extracellularly, and allow the simultaneous measurement of local population activity across many network locations, providing information about the spatiotemporal properties of network dynamics (Le Van Quyen and Bragin 2007;Buzsáki 2004;Rubino et al. 2006). Such arrays can be used both in vitro (Simon et al. 2014) and in vivo, including in humans, where applications include recording from epilepsy patients for precise localisation and investigation of epileptic foci (Schevon et al. 2009(Schevon et al. , 2012, and for use in brain machine interfaces (Maynard et al. 1997;Andersen et al. 2004). These diverse applications make understanding the link between MEA recordings and the underlying neuronal dynamics particularly important.
To study this link, we have created the Virtual Electrode Recording Tool for EXtracellular potentials (VERTEX). VERTEX is implemented in Matlab (Mathworks Inc., Natick, MA, USA), and makes use of established theory of 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. As most such models implement singlecompartment neurons and may not include spatial information (e.g. Izhikevich 2006;Lumer et al. 1997;Potjans and Diesmann 2012), the LFP can only be estimated by some proxy that will not necessarily preserve the spatial and frequency-scaling features of real LFPs . VERTEX helps to address this issue by simplifying the specification of spatially organised cortical network models, and implementing simplified compartmental models that are computationally inexpensive to simulate, but also preserve the spatial and frequency-scaling properties of LFPs elucidated by previous modelling studies Łęski et al. 2013;Lindén et al. 2010Lindén et al. , 2011.
To illustrate how VERTEX can be used in conjunction with MEA experiments, we implemented a model of a neocortical slice exhibiting persistent gamma oscillations under bath application of the glutamate receptor agonist kainic acid in vitro. The model is designed to reproduce the spiking activity of individual neurons during a persistent gamma (30-80 Hz) frequency oscillation, with the neuronal membrane currents driven by the resulting synaptic activity generating the extracellular potential (Nunez and Srinivasan 2006). The persistent gamma frequency oscillation model has several advantages for our investigation. First, the theory of how neocortical persistent gamma arises in vitro, and how individual neurons participate in the network oscillation, has been comprehensively documented (Ainsworth et al. 2011;Whittington et al. 1995Whittington et al. , 2011Fisahn et al. 1998;Buhl et al. 1998;Draguhn et al. 1998;Roopun et al. 2008;Cunningham et al. 2003;2004a;Traub et al. 2005a, b;Pafundo et al. 2013; Bartos et al. 2007). Second, the slice preparation ensures that all synapses are local, so MEA recordings are influenced only by the local circuit dynamics and not by input from other areas. The slice edges provide natural spatial boundaries for what needs to be included in the simulation. Third, synaptic currents rather than intrinsic active membrane currents drive neuronal firing in persistent gamma, so the previously developed theory of LFP generation in passive neurons (Lindén et al. 2010(Lindén et al. , 2011Pettersen and Einevoll 2008) can be used without modification.
Using VERTEX, we have created the first model of neocortical networks that not only reproduces experimentally observed spike patterns, but also produces a biophysically meaningful LFP signal. To illustrate VERTEX's potential for use in conjunction with experimental data, we directly compared the LFPs generated by the model with those recorded by an MEA in macaque temporal neocortex in vitro, allowing us to identify future research directions to address discrepancies between the theoretically predicted and experimentally observed LFPs.

Overview
We developed the VERTEX simulation tool for simulating LFPs produced by large ([100,000) populations of neurons. We first investigated a suitable neuron model for generating LFPs from such populations while remaining computationally tractable. To illustrate VERTEX's capabilities, we used it to position populations of these neuron models into a neocortical slice arrangement, with neuron positions constrained by cortical layer and slice boundaries, and connected them according to current knowledge about the local anatomy of neocortical circuits (Binzegger et al. 2004). We simulated a persistent gamma frequency oscillation in the network, using a simplified model of spike generation in each neuron to generate the network dynamics (Brette and Gerstner 2005). Finally, we compared the simulated LFPs to experimental MEA recordings from macaque temporal neocortex.

LFP generation
The extracellular potential at a point in brain tissue is given by the sum of all neuronal membrane currents, weighted by their distance from the point (Nunez and Srinivasan 2006) assuming constant tissue conductivity (Logothetis et al. 2007;Nicholson and Freeman 1975). Recent theoretical studies have shown that the spatial and frequency-scaling properties of the LFP are affected by the particular spatial arrangement of neurons' dendrites (Lindén et al. 2010(Lindén et al. , 2011. We therefore looked for a reduced compartmental model that would generate extracellular potentials capturing the spatial and frequency-scaling properties elucidated by Lindén et al.'s investigations using detailed cell reconstructions, while remaining computationally tractable to simulate in large numbers. The reduced compartmental model should create a similar spread of currents across its compartments to an equivalent morphologically reconstructed neuron given the same input. A compartment's membrane current depends on the neuron's axial resistance as well as on its membrane resistance and capacitance. We therefore chose a reduced model that conserved these quantities, while containing a minimal number of compartments. The compartmental reduction method of Bush and Sejnowski (1993) fulfils these requirements, producing compartments with a length equal to the mean length of the compartments they are representing in the full model. This creates a reduced model of the same length as the original reconstruction, but with a smaller membrane area, smaller lateral spread of the dendrites, and fewer than ten compartments (Online Resource, Fig. ESM2).

Validating the reduced LFP generation model
We tested the effects of this reduction on the generated LFP by reproducing the experiments detailed by Lindén et al. (2011). Ten thousand model neurons with passive membrane dynamics and the same morphology were positioned randomly within a 1-mm radius cylinder, with uniform spatial distribution and constant soma depth. One thousand synapses (excitatory, current-based, single-exponential type with time constant 2 ms and fixed amplitude 50 pA) were placed randomly on the compartments of each neuron, with uniform density with respect to membrane area. Each synapse received an independent Poisson spike input train with a rate of 5 Hz. LFPs were calculated at the centre of the population, at five depths. The magnitude of an LFP signal was defined as its standard deviation. The LFP range was calculated by varying the population radius from 0 to 1 mm and measuring the radius at which the LFP magnitude reached 95 % of its value at the maximum 1-mm radius (Lindén et al. 2011). We repeated this procedure for the three neuron types used by Lindén et al.: layer 2/3 (L2/3) pyramidal, layer 4 (L4) spiny stellate, and layer 5 (L5) pyramidal. We compared LFPs generated by the morphological reconstructions of these neuron types described by Mainen and Sejnowski (1996)-hereafter referred to as Mainen cells-with the LFPs from reduced versions of these models created using Bush and Sejnowski's method (1993)-hereafter referred to as Bush cells.
The results of these experiments are shown in Fig. 1. For each neuron type, the LFP range and magnitude in each layer for the population of Bush cells are close to those for the population of Mainen cells. The LFP range is smallest in the soma layer (\250 lm) with the range increasing in the layers above and below the soma, while the LFP magnitude is largest in the soma layer and decreases in the layers above and below the soma. The differences between the results for the L4 spiny stellate models are small, so we concentrate on the pyramidal neuron population results. For the L2/3 pyramidal neurons, the LFP spatial range in the soma layer is very similar between the Bush and Mainen populations, but above and below this layer the discrepancy increases, with the largest difference of 200 lm in L1. The range differences in all other layers are B110 lm. For the L5 pyramidal neurons, the LFP spatial range difference is again smallest in the soma layer, and \100 lm in layers 4 and 1. The largest difference is 320 lm in L2/3.
To see how important these discrepancies were within the context of the general biological variability of neuronal morphology, we repeated the simulations with neuron populations containing pyramidal cells reconstructed from several different real neurons. These were downloaded from the NeuroMorpho.Org database (Ascoli et al. 2007)further details on the models we used are provided in the Online Resource (Table ESM10). We used ten further groups of L2/3 cat pyramidal neurons, and one further group of L5 cat pyramidal neurons (this was the only other cat L5 pyramidal neuron currently available in the database; we did not use L5 pyramidal cells from other species as the size differences in neurons between species could have provided misleading results). The results of these simulations are plotted in Fig. 1b as light-red dashed lines for the extra L2/3 pyramidal populations, and light blue circles for the extra L5 pyramidal population. The extra simulation results show that the LFP range and magnitude in the Bush neuron populations generally fall within a biologically reasonable range; while the reduced models are not ideal substitutes for the morphological reconstructions, the errors incurred by the reduction method are similar to those introduced by neglecting morphological diversity in reconstructed neuron model populations. The general profile of the LFP across the layers, at least, is preserved adequately.
We also checked the power spectra of the simulated LFPs to make sure the Bush model populations reproduced similar frequency-scaling properties to the Mainen cell populations. Figure 1c shows that, in each layer, the 95 % confidence intervals for each model type overlap over the range of frequencies from 2 to 450 Hz (the overlap continues down to 1 Hz; this is not shown to improve the plot resolution at higher frequencies).
The results in Fig. 1 were generated using uncorrelated synaptic inputs over the entire dendritic tree of each neuron in each population, with neurons all positioned at the same height in their respective layers. This simplified setup was used so that a comparison could be made with the results previously reported by Lindén et al. (2011), but we also wanted to check whether the reduced models would still be suitable approximations to use for a more realistic situation, with neurons placed at varying depths within their layer, receiving correlated inputs. As our particular interest was simulating network gamma oscillations, in which pyramidal neurons receive highly correlated inhibitory synaptic input to their perisomatic regions, we repeated the previously described experiments measuring the LFP magnitude and range, but positioned each neuron's 1,000 synapses onto its soma compartment [we only repeated the simulations for the pyramidal neuron morphologies, as the LFP spatial profile for the spiny stellate cells was shown not to change significantly with correlated input (Lindén et al. 2011)]. In the previous experiments with no correlations between synaptic inputs, each synapse was assigned an independent Poisson spike train, for a total of 10,000 9 1,000 = 10 7 independent spike trains at 10 7 synapse locations. To introduce input correlations, we followed the same method as Lindén et al. (2011). Each synapse in the model was now assigned a spike train drawn without replacement from a finite pool of pregenerated spike trains. By reducing the number of Poisson spike trains in the pool so that some synapses shared a common input pattern, we could control the level of input synchrony to the neurons. The resulting input correlation is given by the total number of synapses per neuron divided by the number of independent spike trains (Łęski et al. 2013). To simulate highly correlated input, we used 2,000 independent spike trains, resulting in an input correlation of 1,000/2,000 = 0.5 (i.e. any two neurons share on average 1,000 9 0.5 = 500 common input spike trains).
For these simulations, we also introduced random variability in the soma depth of the neurons. We distributed L2/ 3 pyramidal neuron somas between -334 and -534 lm, and L5 pyramidal neuron somas between -970 and -1,170 lm from the cortical surface. These ranges ensured that the neuron somas remained within the correct layer boundaries, and that their apical dendrites were not positioned above the cortical surface. Figure 2 shows the spatial profiles of the LFP for the different populations. In these simulations, we measured the LFP at 50 intervals, to see how well the Bush models preserved the LFP at this level of detail. We used 11 electrode points in L1 and L2/3 for the L2/3 populations, and 26 electrode points spanning all layers for the L5 populations.
Both the range and magnitude profiles show that the LFP from the Bush population matched the LFP from the Mainen population well, again within the bounds of the LFP profile of the extra comparison populations. The minimum range and magnitude in the L2/3 populations are just above the minimum soma depth, and a few 100 lm above the minimum soma depth in the L5 population. This depth is where the synaptic currents at the soma are approximately balanced by the opposite return currents in the dendrites; below and above this minimum point, the somatic or the apical dendritic currents dominate the LFP signal, respectively. These simulations also show substantial overlap of the 95 % confidence intervals for the power spectra at each electrode (Fig. 2b). The biggest discrepancy between the LFP power spectra for each model occurs around the level of the LFP range minimum. The LFP power up to 100 Hz is reliably reproduced at every measurement point, and up to 450 Hz at all but one point with the L2/3 populations. This point corresponds to the point at which the LFP range and magnitude are lowest. The reduced accuracy at higher frequencies in the L5 models should be taken into account if frequencies above 100 Hz are analysed in models containing L5 pyramidal cells.
Our results suggested that we could use the reduced neuron models in VERTEX simulations with some confidence that the resulting simulated LFPs would be close to LFPs simulated from equivalent morphologically reconstructed neurons, in magnitude, spatial extent, and frequency content.

The VERTEX simulation tool
To simulate large networks, we wrote custom Matlab software to setup neuron populations, position them, connect them together, and simulate their dynamics and the resultant LFPs. We designed this simulation tool to be easily adaptable to create models of any layered brain tissue containing populations of spiking neurons (Fig. 3). Model parameters are specified by the user in Matlab structures, defining: 1. Neuron group properties (for each group: the neurons' compartmental structures, dimensions and positions, electrotonic properties, spiking model parameters) 2. Connectivity (for each presynaptic group: number of efferent synapses per layer per postsynaptic group, allowed postsynaptic compartments to connect to contact, axonal conduction speeds, neurotransmitter release times, synapse dynamics) 3. Tissue properties (dimensions, layer boundaries, neuron density, tissue conductivity) 4. Recording settings (IDs of neurons to record intracellularly, extracellular electrode positions, sampling rate) 5. Simulation settings (simulation length, time-step, number of parallel processes) A model is initialised by positioning the specified number of neurons from each group within the slice and layer boundaries, pre-calculating distances from the neuron compartments to the virtual electrodes, generating each neuron's connections based on its position, axonal arborisation extent in each layer, and expected number of efferent connections, and initialising the synapses (see ''Experimental Procedures''). At this point, the initialised model can, optionally, be saved to disk as MAT files. Functionality to export to NeuroML (Gleeson et al. 2010) is currently under development.
When the simulation is run, recordings (intracellular, LFPs, spike times) are automatically saved to disk at userspecified time intervals. The simulation run can be performed in serial or parallel (requires Matlab Parallel Computing Toolbox). After the simulation is finished, these files are loaded and recombined for analysis. Our design allows the model to be used with minimal programming knowledge, though as Matlab is a high-level, interpreted language, more experienced programmers can make modifications relatively easily.

Simulation speed and memory usage
While Matlab code may run more slowly than equivalent code in compiled programming languages, performance can be dramatically improved through code vectorisation, which minimises the impact of code interpretation overheads (Brette and Goodman 2011). The Matlab Parallel Computing Toolbox allows further performance improvements by providing a simple way to parallelise computations on multicore computers or over networks. These factors, as well as its ease of use, popularity in the neuroscience community, the ability to perform simulations and analysis in the same environment, and the welldeveloped interface for integrating C or Fortran functions for future performance enhancements influenced our decision to write VERTEX in Matlab. To give the user an idea of the performance improvement over using the other current extracellular potential simulation tool LFPy (Lindén et al. 2014)-a Python package for simulating extracellular potentials with NEURON (Hines and Carnevale 1997;Hines et al. 2009)-we performed equivalent simulations using layer 5 Bush pyramidal neurons in LFPy and in VERTEX (no synapses, one random fluctuating current per neuron, 0.03125 ms step size and 32 000 Hz sample rate). LFPy took *278 min to simulate the LFP from 10,000 neurons at 50 electrode points, while VER-TEX running in serial mode took *18 min to simulate the LFP from 10,000 neurons at 50 electrode points (both running on an Intel Xeon E5640 2.66 GHz workstation). While this performance improvement is important for our purposes, it should be noted that LFPy is designed to simulate extracellular potentials from single cells rather than large populations. Indeed, as the code interpretation overhead begins to dominate VERTEX's calculation times in small simulations, running the same model but with only one neuron in the population took *227 s in VER-TEX but \2 s in LFPy. VERTEX is also not suited to running models containing neurons with very many compartments, because the Runge-Kutta integration method becomes unstable as the number of compartments increases (though we aim to address this limitation by implementing implicit integration methods in future releases). LFPy, therefore, remains the superior tool for modelling extracellular potentials around single neurons, while VERTEX's strength lies in simulating LFPs in large-scale networks.
To show how performance improves in parallel mode, we compared the run times for two network models, one large (123,517 neurons with on average 1,835 synapses per neuron) and one small (9,881 neurons with on average 256 synapses per neuron), using VERTEX on a single multicore computer (Fig. 4). Each model contained two populations: layer 5 pyramidal (P5) neurons and layer 5 basket (B5) interneurons. Spike rates in each small model (large model) simulation were *6 Hz (*7 Hz) and *24 Hz (*31 Hz) for the P5 and B5 neurons, respectively. The large model shows linear speed-up with increasing number Fig. 3 Overview of the VERTEX simulation software. a Simulation workflow. The user provides parameters as Matlab structures to setup the neuron populations, position them in layers, connect them together, and simulate their dynamics and the resultant LFPs. Functionality to export to NeuroML is currently under development.
b Example program structure. The main simulation program only requires calls to the initNetwork() function and the runSimulation() function, with the information required to setup the simulation specified in separate script files Brain Struct Funct (2015) 220:2333-2353 2339 of cores for model initialisation and close-to-linear speedup in simulation time. The speed-up for the small model is sub-linear: as the interpretation overhead for a vectorised operation on a small matrix is the same as on a large matrix, this overhead starts to dominate the calculation times below a certain number of neurons (Brette and Goodman 2011). Therefore, splitting already small neuron state matrices between more processes does not significantly improve performance. This limit is not reached in larger models. Figure 4 also shows how increasing the number of virtual electrodes affects simulation speed. Model initialisation times are affected proportionally more than model run times using more electrodes, in both large and small models, though in the small model the proportional impact from adding electrodes to initialisation time was greater than in the large model. This is because the large model not only has more neurons, but also more synapses per neuron. The increase in time spent connecting the neurons is proportional to the number of synapses, while the increase in time spent calculating constants for the LFP measurements is proportional to the number of compartments (roughly proportional to the number of neurons).
The size of the simulated network is limited by the amount of RAM available. As an example, we tested scaled configurations of our neocortical slice model (described below) using single-core and multi-core computers: an iMac with 4 GB RAM supported a serial simulation with *25,000 neurons, a 16 GB Linux machine supported a simulation of *100,000 neurons in both serial and parallel modes, and our Linux server with 120 GB RAM supported a simulation of *700,000 neurons. In addition to increasing the memory on a single machine, VERTEX could be run across a network of computers using the Matlab Distributed Computing Server. On a network of 16 of our 4 GB RAM iMacs, for example, the simulation size could scale to *400,000 neurons. In summary, existing processing environments of experimental and computational laboratories can be sufficient for running detailed simulations of brain tissue activity.

Spike import
Network dynamics can be simulated directly by providing the model neurons with a spiking mechanism-we used the adaptive exponential (AdEx) mechanism (Brette and Gerstner 2005), which we include in VERTEX. Alternatively, previously generated spike times (output from another simulator, for example) can be imported into the simulation. The neurons whose spike times are imported are then specified with purely passive membrane dynamics. We used the spike import feature to run the control experiment to confirm that the AdEx spiking mechanism has a negligible impact on the simulated LFP (Online Resource, Fig. ESM1). Running models using imported spike times is similar to the approach used by Lindén et al. (2011) to link spiking output from a cortical model implemented in the NEST simulator (Gewaltig and Diesmann 2007) to their LFP generating model implemented in LFPy. However, we consider imported spikes to have been emitted by neurons from within the population we are modelling; imported spikes are delivered to target neurons according to the generated connectivity matrix rather than pre-assigned to postsynaptic targets. By contrast, Lindén et al. (2011) considered the spikes from NEST-simulated neurons as external inputs to the neurons in the LFPy simulation, so they were delivered to synapses without a connectivity model within the LFPy-simulated population. The practical effect of this is that our software is better suited to modelling the LFP resulting from intrinsic network dynamics, when connectivity is known or when different spatial connectivity models are to be tested. Input from external populations can be simulated by specifying a population of single-compartment neurons and setting this population's output using the spike import functionality. As single-compartment neurons do not contribute to the extracellular potential (Pettersen et al. 2012), VERTEX ignores them in its LFP calculations. This population can, therefore, be considered as providing ''external'' input from a distant population.

Neocortical slice model
To demonstrate the capabilities of VERTEX for simulating LFPs in large neuron populations, we created a neocortical slice model to use in conjunction with MEA experiments in vitro (Fig. 5). The model comprises fifteen neuron groups, defined in Table 1. It is designed to contain a similar number of neurons to the comparison experimental slice. This was calculated to be 175,421 neurons, based on the slice dimensions and neuron density. The slice has clear spatial boundaries: neurons cannot be positioned outside of the slice edges, and axons cannot 'wrap around' these boundaries. We therefore required a connectivity model that would produce a suitable number of synapses given the large number of neurons, and that took into account each neuron's position in relation to the slice boundaries. We used anatomical data from Binzegger et al. (2004) to specify the numbers of connections between neuron groups, and a 2D Gaussian spatial profile to model the decay in connection probability with increasing distance from a presynaptic neuron (Hellwig 2000). The standard deviation parameter of the Gaussian profile was set using axonal arborisation radius measurements reported by Blasdel et al. (1985); Fitzpatrick et al. (1985), as adapted by Izhikevich and Edelman (2008). These were different for each neuron group in each layer (see Online Resource, Table ESM4). Finally, we modelled the effect of slice cutting on connectivity by reducing the number of connections a presynaptic neuron could make by the proportion of the integral of its Gaussian connectivity profile that fell outside the slice boundaries (Eq. 3). The connectivity generation code in VERTEX implements this connectivity model automatically, though the user can also specify a uniform spatial connection probability and/or ignore slice cutting effects. VERTEX also allows users to specify specific target compartments on postsynaptic neurons that presynaptic neurons are allowed to connect to. We used this feature to incorporate known details about the dendritic regions targeted by different presynaptic neuron types-basket interneurons only make connections with pyramidal cell somas and their two adjacent compartments, for example. We used a similar pattern of connectivity to that described by Traub et al. (2005b); details are provided in the Online Resource (Supplementary Methods: Connectivity and Table ESM7). Incorporating this detail into the model is important, as the locations of synaptic inputs onto the neurons will affect the locations and sizes of the currents that contribute to the simulated LFP. Figure 6 shows the number of connections between neuron groups compared with the original numbers specified by Binzegger et al. (2004). The proportional reduction in synapses is not the same for each connection type because of the varying axonal arborisation radii. These reductions are important to consider when assessing the effect of connectivity changes on dynamics, but they illustrate that the general profile of connections between neuron groups is not substantially altered-connections from P2/3 to P2/3 and P5 neurons remain the most numerous, for example. Modelling thinner slices, or different axon arborisation profiles, could lead to the over-or under-representation of particular connections in the model.

Modelling persistent gamma oscillations
To make a comparison with experimental data, we generated a persistent gamma oscillation in the model by applying random currents to all neurons , and adding an AdEx spiking mechanism to the somatic compartments (see Online Resource). In slice experiments with nanomolar kainate concentrations, this activity regime is driven by L2/3, where neurons receive noisy excitatory drive from the excited axonal plexus of L2/3 pyramidal neurons (Ainsworth et al. 2011;Cunningham et al. 2003Cunningham et al. , 2004b. We simulate this by providing a relatively large noisy current to P2/3 neurons, similar to Ainsworth et al. (2011); . We set synaptic strengths (based on Traub et al. 2005b) and noise currents to match the spiking activity and observed membrane potential fluctuation sizes reported in previous studies in vitro. Model parameters are given in tables ESM1-ESM9.
As described in previous experiments (Ainsworth et al. 2011;Cunningham et al. 2003Cunningham et al. , 2004bTraub et al. 2005a, b), P2/3 neurons spike infrequently, while B2/3 neurons spike on most oscillation periods. Excitatory neurons in L4 do not take part in the oscillation (though still spike infrequently), while L4 interneurons are weakly entrained to the oscillation. In addition to the L2/3 gamma, the comparison slice exhibited increased gamma power in part of the infra-granular layers (see Fig. 9a,electrodes 6,7,16,17,26,27), presumably caused by L5 as described by Ainsworth et al. (2011). We therefore used a relatively high coupling strength of P5 to B5 and NB5 neurons and a larger noisy drive current to L5 neurons to enable the L2/3 gamma to generate gamma in L5. The L5 gamma oscillation also weakly entrained L6 neurons to the oscillation.
The resulting spiking behaviour is shown in Fig. 7, which shows a spike raster for 5 % of the neurons in the model, along with example somatic membrane potential traces for each neuron group. The spike raster reveals that neurons near the slice x-boundaries (neurons nearest the cyan boundary markers in Fig. 7) are less strongly entrained to the oscillation than neurons in the centre of the slice, because they receive fewer inhibitory inputs than more central neurons (neurons closer to the edge of the slice have more connections removed by slice cutting than those towards the middle of the slice, because they lose proportionally more of their axonal arborisation). To demonstrate how the oscillation is generated by the interaction of the excitatory and inhibitory populations, we simulated activity in the model under four different conditions: firstly the original case described above (connection weights in Table ESM5), secondly with P2/3 to B2/3 synapses reduced to 1 % of their original weight, thirdly with B2/3 to P2/3 synapses reduced to 1 % of their original weight, and fourthly with the original synapse weights but increased input current to the B2/3 population (1.5 times the mean and standard deviation used in the original simulation values given in Table ESM9). Simulation results using these different configurations are plotted in Fig. 8, which shows that both P2/3 to B2/3 synapses and B2/3 to P2/3 synapses are necessary for the generation of a population gamma oscillation in the model. Without these connections-or with their strengths severely reduced-no oscillation emerges. This oscillation mechanism is the same as the ''weak'' pyramidal-interneuron network gamma (PING) model described by . Firing in a subset of P2/3 cells, which are densely connected with B2/3 neurons with strong synapses, causes a population spike from the B2/3 cells. This suppresses the network until the P2/3 neurons that receive the most input from the stochastic drive reach threshold. This subset of P2/3 neurons then fires, causing another B2/3 cell population spike, and so the oscillation continues . Figure 8m-p shows that the oscillation is also suppressed in our model when the driving current to B2/3 cells is increased, allowing them to suppress P2/3 cell firing. This is in line with the gamma suppression mechanism described by . Figure 8 also demonstrates that the gamma oscillation in layer 5 is dependent on a gamma oscillation occurring in layer 2/3: layer 5 gamma is suppressed in each of the cases where layer 2/3 gamma is suppressed. Firing rates for each population in each case are given in Table ESM11.
Having verified that the model produced the expected spiking output and that the gamma oscillation was being generated by the correct mechanism, we looked at the simulated LFPs and compared them with those recorded in vitro. Figure 9 shows a comparison over the whole electrode array between the model and the experimental recordings. Figure 9a shows the shape of the experimental neocortical slice with, as predicted by previous research, strong gamma power in the supra-granular layers. The gamma power at each electrode is highly variable, resulting in a patchy power map. This is not captured by the model, whose structure is homogeneous along the x-axis. However, the phase inversion between layer 1 and layer 2, illustrated in Fig. 9b, c, emerges in the model (Fig. 9e, f) from the positioning of current sinks and sources on the P2/ 3 neurons during the gamma oscillation. This is in agreement with the source-sink interaction mechanism of phase inversion demonstrated experimentally in kainate-induced gamma oscillations in entorhinal cortex in vitro (Cunningham et al. 2003). The cross-correlations between electrodes shown in Fig. 9c, f also reveal how the strong gamma oscillation in L2/3 dominates the across the electrodes more than in the experimental recordings. This is, again, a result of the relatively homogeneous activity along the x-axis in the model, meaning that the LFP signal created by the gamma oscillation is not degraded by influences from the non-oscillating areas in the slice as occurs in vitro. Our model, though not capturing all the details of the experimentally measured network dynamics, provides a starting point for further investigations into cortical dynamics on this spatial scale, allowing for better integration of theory and experiment. Boundaries between neuron groups marked in cyan. Note strong persistent gamma oscillation in L2/3, with weaker oscillation in L5. b Example soma membrane potential plots for the various neuron types. Most neurons fire sparsely, while B2/3 and B5 neurons fire on most oscillation periods. Note occasional spike doublet firing in the B2/3 neuron. Spikes are cut-off at V t ? 5 mV in the simulation; we extend them up to 10 mV here for illustrative purposes. c Close-up of P2/3 neuron soma membrane potential (cut-off-45 mV). Scale bar 5 mV Fig. 8 Illustration of the gamma oscillation mechanism in the model. a Spike raster of 250 ms from a simulation of a model with the same parameters as that shown in Fig. 6. For clarity, spikes from only 5 % of the neurons are shown. A gamma oscillation is apparent in layers 2/3 and 5. b Zoomed spike raster showing only neurons in layer 2/3. Spikes from only 1 % of the neurons are shown. c LFP recording from the virtual electrode with the highest gamma power in the LFP. d Power spectrum of the LFP from this electrode, calculated for 1.5 s simulation time, showing a clear gamma peak. e-h same as a-d, but with synaptic weights from P2/3 cells to B2/3 cells reduced to 1 % of their original value. e-f show B2/3 cell firing is greatly reduced, as they are not receiving excitation from the P2/3 cells. No gamma oscillation emerges. i-l same as a-d, but with synaptic weights from B2/3 cells to P2/3 cells reduced to 1 % of their original value. B2/3 cells fire rapidly and randomly: they are driven by the P2/3 cells but they cannot synchronise them as their synapses are too weak. No gamma oscillation emerges. m-p same as a-d, but with the mean and standard deviation of the stochastic input current to the B2/3 cells increased by 50 %. P2/3 cell firing is suppressed by the increased B2/ 3 cell firing, so no gamma oscillation occurs

Discussion
We have developed the VERTEX tool for simulating LFPs generated by large neuronal populations. VERTEX is easily customisable, and makes use of recent developments in simulation techniques and insights from our experiments with simplified neuron models to reduce simulation times for LFPs generated by large networks. To illustrate how VERTEX can be used in conjunction with experimental MEA data, we simulated kainateinduced persistent gamma oscillations in a large-scale neocortical slice model. The model reproduces the spiking activity underlying persistent gamma, and generates the theoretically predicted LFP from this activity. We compared this simulated LFP with Utah array recordings of persistent gamma from macaque temporal neocortical slices. The model predicted the oscillation phase inversion between L2/3 and L1, but not the spatial variation in gamma power within layers, suggesting directions for further research into the cause of the spatial discrepancies between theoretically predicted and experimentally measured LFPs. c Cross correlation of signals from electrodes 41-44 with signal from electrode 42, illustrating phase inversion in the signal from electrode 41. This electrode was identified as being in layer 1 by post hoc histology (not shown). Gamma map and cross-correlations estimated from 18 s of data. d-f as a-c, but for the neocortical slice model (gamma map and cross-correlations estimated from 1.5 s of simulation data) Speed of the VERTEX simulator Parallel computing and code vectorisation allow VER-TEX to simulate network activity and LFPs in reasonable time on hardware that is available to most scientists. We showed typical simulation times and how performance scales with increasing numbers of parallel processes in Fig. 4. However, performance could be improved further by rewriting some of the Matlab code in C or Fortran, which could be incorporated into Matlab via its MEX interface. In particular, the spike queuing and delivery code would benefit from this approach when simulating networks with high spike rates, as it is only vectorised over individual spikes. High spike rates can result in longer simulation times as the spike queue interpretation overhead increases. This is, therefore, a priority for future VERTEX development. However, the pure Matlab versions of VERTEX will continue to be maintained, as some users may not have access to a suitable C or Fortran compiler.

LFP simulation: spatial properties and resolution
We found that the compartmental reduction method described by Bush and Sejnowski (1993) created neuron models that, in a population, reproduced the spatial properties of the LFPs generated by the equivalent full morphological reconstructions to a reasonable degree of accuracy. Where there were large discrepancies, they were close to or fell within the range of the spatial values measured in several further populations of different morphologically reconstructed neurons. The suitability of this reduced model allows VERTEX to simulate LFPs from large networks in reasonable time.
The largest compartment in the reduced models was 400 lm long, which is the inter-electrode distance in a Utah array. New, very high-density MEAs with several 1,000 electrodes can record with such high spatial resolution as to enable the visualisation of individual dendritic tree and synapse activity in detail (Frey et al. 2009), or to record the spiking activity of thousands of neurons (Berdondini et al. 2009), making our reduced neuron models unsuitable for use in conjunction with these experiments. The array described by Frey et al. (2009) is designed to record only from a subset of 126 electrodes concurrently, allowing very high-resolution recordings from small areas, but making it unsuitable for recording the wider population activity that our model is designed to capture. The 4,096 electrode array presented by Berdondini et al. (2009) can record simultaneously from all electrodes, allowing the detailed visualisation of signal propagation through a network. However, this array is designed for capturing the spike times of thousands of individual neurons rather than investigating the properties of extracellular signals. Given the spatial smearing of LFP signals, it would not be appropriate to use this type of array to investigate LFPs across active neural circuits. Additionally, very high-density arrays are new technologies with usage and data analysis techniques still under development. Lower density MEAs will remain useful for studying neuronal population activity for the foreseeable future, especially given the Utah array's approval for use in humans. As higher density arrays become more common, we anticipate that advances in computing speed [through, for example, use of generalpurpose graphical processing unit (GPGPU) computing (Brette and Goodman 2012), already a feature of the Matlab Parallel Computing Toolbox] will permit the simulation of large populations of higher resolution neuron models if desired.

Slice model properties
To demonstrate our simulation approach, we constructed a model of a neocortical slice. We combined the connection probabilities given by Binzegger et al. (2004) with axonal arborisation radii measured in macaque visual cortex Fitzpatrick et al. 1985), and use a Gaussian kernel as suggested by the data from Hellwig (2000) as the decay in connection probability away from the soma. This approach allowed us to calculate the number of connections removed by slice cutting for each neuron, and reduce its number of connections accordingly when initialising the model.
Our anatomical model results in spatially uniform neuron densities and connectivity statistics, with small decreases in connection numbers nearer the slice boundaries. However, the recordings from the experimental slice illustrate substantial inhomogeneities in gamma power between electrodes, even within layers, that are not seen in the model. These could be caused by spatial variations in synapse densities and strengths, neuron group densities, neurons' dynamical properties, gap junction densities and strengths, or axonal plexus properties. While our software does not currently allow specification of gap junctions or axon properties, the other potential inhomogeneities can be investigated further in conjunction with experiments in vitro: VERTEX includes functions to modify parameters in spatially localised regions, allowing spatially inhomogeneous tissue to be modelled. As the results of these modifications can be compared directly with extracellular recordings, theoretical predictions can be tested even when spiking data are lacking. For example, spatial variations in synapse densities may be caused by the ''patchy'' projections made by excitatory neurons (Binzegger et al. 2007;Voges et al. 2010b;Bauer et al. 2012;Douglas and Martin 2004).
Future research could incorporate the patchy projection model of Voges et al. (2010a, b) into our slice model to investigate how patchy connectivity affects network activity and resultant LFP across the slice.
We model the cortical layers as being flat, with boundaries at constant depths below the cortical surface. Neocortex is a folded structure, though, which is apparent even at the small scale of the slice-note the curved shaded regions in Fig. 9a showing the cortical surface and white matter boundaries, as well as the curved profile of gamma power across the MEA. Curves add further complications to the other inhomogeneities discussed above, in terms of neuronal densities, layer thicknesses and axonal arborisation variations. Additionally, the alignment of pyramidal apical dendrites is perpendicular to the cortical surface, so the alignment of the current dipoles arising from synaptic currents on pyramidal dendrites (Lindén et al. 2010;Nunez and Srinivasan 2006) varies across space, with implications for the measured LFP. VERTEX functions for specifying curved layer boundaries are currently under development so that future experiments can investigate the effects of curved surfaces on the measured LFP.

Further considerations for LFP simulation
In its current state, VERTEX is designed for investigating LFPs in medium to large-scale spiking neural networks, as these are most often used for modelling the activity of large neural populations. We have, therefore, only implemented simplified neuron models that do not include realistic active conductances that produce, for example, backpropagating dendritic spikes or sub-threshold membrane oscillations, which would also contribute to the LFP. As gamma oscillations are driven by synaptic interactions between populations, we consider this to be a reasonable simplification for our neocortical slice model. When investigating other dynamical regimes-such as subthreshold oscillations in the absence of spiking (Hutcheon and Yarom 2000)-this simplification may not be appropriate. However, VERTEX will still be useful for investigating many research questions even with these simplifications. For example, most previous spiking neural network models use highly simplified neuron models, for which there is no general, reliable method for estimating the LFP . VERTEX allows researchers to implement similar networks using neuron models that produce a spatially realistic LFP, so that they can directly compare the LFPs produced by the spiking activity in their models to experimental data. Such comparisons may reveal both agreements and discrepancies between model and experiment, which might not have been apparent from comparisons of spiking alone. This was the case for our slice model: we could not directly compare spiking across space as it was massively under-sampled in vitro, but the simulated LFPs based on our prior knowledge of neuronal firing during gamma oscillations revealed that we can account for the observed phase inversion between L2/3 and L1, but cannot account for the spatial variation in gamma power with our current model. Future research to address this discrepancy is discussed above.
Several exciting experimental results have recently shown that neuronally generated electric fields impact on the membrane potentials of nearby neurons without requiring any synaptic contact. Such ''ephaptic'' coupling of neurons was investigated in models (Holt and Koch 1999) and, more recently, confirmed in experiments showing that such interactions could modulate oscillatory network activity (Fröhlich and McCormick 2010), entrain action potentials (Anastassiou et al. 2011) and potentially contribute to the spread of epileptiform activity (Zhang et al. 2014). We have purposefully ignored the contribution of ephaptic interactions in our model for the sake of simplicity, and have not incorporated the simulation of ephaptic coupling into the VERTEX simulator. While the results reported by Fröhlich and McCormick (2010) suggest that endogenous electric fields should be taken into account in models of oscillatory activity, they concentrated on neocortical slow oscillations, which are greater in amplitude than the gamma oscillations we modelled. However, the role of ephaptic interactions on network activity under different conditions must be investigated further. As VERTEX can simulate the LFP at arbitrary locations in a network, it would be possible to incorporate an ephaptic coupling mechanism that depended on the LFP. However, doing this rigorously would entail measuring the LFP near every compartment in the model, which is not feasible. Developing suitable approximation methods for incorporating realistic ephaptic coupling is, therefore, an important direction for future research. Similar methods could also be used for simulating artificially applied electric fields/currents, such as from extracellular stimulating electrodes.
Finally, the VERTEX simulator assumes a purely resistive, constant and homogeneous extracellular conductivity, with no frequency dependence (Pettersen et al. 2012). The extracellular medium's frequency filtering effects are not currently known for certain ): some results have demonstrated an intrinsic lowpass filtering effect (Gabriel et al. 1996;Dehghani et al. 2010) potentially created by ionic diffusion (Bédard and Destexhe 2009), though direct measurements in macaque cortex in vivo found minimal frequency filtering from intrinsic tissue properties (Logothetis et al. 2007). If a frequency-dependent effect of the extracellular medium is confirmed by future studies, Eqs. (1), (2) (see ''Methods'') can be modified to take this into account (Pettersen et al. 2012).

Conclusion
We have described the VERTEX simulation tool for simulating LFPs in large neuronal populations. VERTEX includes functionality for generating spatially constrained networks of several neuron populations, whose parameters are easily specified in Matlab structures. ''Virtual electrodes'' can be positioned at arbitrary locations in the model to simulate the LFP generated by the network. Parallel computing and code vectorisation, as well as the use of reduced compartmental neuron models, allow VERTEX to simulate network activity and LFPs in reasonable time. Finally, we simulated LFPs from a neocortical slice model and compared them with LFPs recorded from macaque neocortex in vitro, illustrating new avenues for research into spatial variations in the LFP signal. We hope that the VERTEX and our neocortical slice model will prove useful to other researchers investigating the relationship between neuronal circuit dynamics and experimental or clinical brain tissue recordings.

Software and simulation methods
Spatial LFP characteristics of each individual compartmental neuron model were tested using LFPy as described in the Results section. LFPy simulations used a 0.125 ms time-step and NEURON's standard implicit Euler numerical integration method. Each simulation was run for 1,250 ms simulation time, and the first 250 ms were discarded to remove simulation start-up effects.
VERTEX is implemented in Matlab. It uses the Matlab Parallel Computing Toolbox for parallelisation, though it can also be run serially. Equations are integrated numerically using a second-order Runge-Kutta method (Press et al. 2007); we used a 0.03125 ms time-step unless otherwise specified. VERTEX incorporates the methods outlined by Morrison et al. (2005) for parallel simulation, and the algorithms and data structures described by Brette and Goodman (2011) for code vectorisation.
In both LFPy and VERTEX, extracellular potentials are calculated by summing the membrane currents of each compartment, weighted by distance from the electrode tips. The line-source method (Holt 1998), as used previously by Lindén et al. (2010Lindén et al. ( , 2011 and Pettersen and Einevoll (2008), is used to calculate the contribution from all dendritic compartments to the LFP, U dend , at a measurement point r and time t: (Holt 1998), where I mem,k is the membrane current from compartment k, r ex is the extracellular conductivity, Ds k is the length of compartment k, q k is the perpendicular distance from compartment k, h k is the longitudinal distance from the end of compartment k, and l k = Ds k ? h k is the longitudinal distance from the start of the compartment. As in Lindén et al. (2010Lindén et al. ( , 2011 and Pettersen and Einevoll (2008), somatic compartments are modelled as point current sources in VERTEX: Nunez and Srinivasan 2006), where r s is the distance between point r and the centre of soma s. The total extracellular potential measurement at point r is then In all simulations, we used a value of r ex = 0.3 S/m (Hämäläinen et al. 1993). LFPy simulations to test the model reduction method (Bush and Sejnowski 1993) were run on an Intel Core-i7 based PC running Ubuntu Linux 11.10 using a pre-release version of LFPy, NEURON 7.1 and Python 2.7.2, and an Intel Xeon E5640 workstation running Linux Mint 16 using LFPy 1.0 with NEURON 7.3 and Python 2.7.5. The LFPy vs. VERTEX performance comparison was run on the same Intel Xeon E5640 workstation, using Matlab 2013a. All other VERTEX simulations were run on a 48-core HP ProLiant server running CentOS Linux 5.8 with Matlab R2012b. Parallel simulations were run on 12 cores unless otherwise specified. The code, as well as documentation and tutorials, will be made available at http://www. dynamic-connectome.org/ upon publication.

Neoortical slice model
The neocortical slice model contained fifteen neuron populations, defined by location, connectivity, morphology, dynamics, and type of neurotransmitter effect (excitatory or inhibitory). We used the naming convention from Binzegger et al. (2004) as adapted by Izhikevich and Edelman (2008), defining the groups listed in Table 1 (full model parameters are given in Tables  ESM1-ESM9). Individual neurons are represented by compartmental models with 7, 8 or 9 compartments, derived from the neuron models given by Mainen and Sejnowski (1996) using the compartmental reduction method of Bush and Sejnowski (1993). Compartmental structure and neuron parameters are given in Fig. ESM2; Tables ESM1 and ESM2. Our connectivity data are from cat visual cortex (Binzegger et al. 2004), so we took parameters for the neuronal density and layer boundaries from the same source. We scaled the layer boundaries to increase the total cortical depth to 2.6 mm, which was approximately the cortical depth in the comparison experimental slice (established by post hoc histology, not shown).
VERTEX is designed for specifying models in 3D space, giving all neuronal compartments 3D start and end coordinates. For the neocortical slice, we defined the z-axis to be the cortical depth from white matter through the layers to the cortical surface, with the border between layer 6 and the white matter set to z = 0 mm. The x-and y-axes ran parallel to the cortical surface, with the y-axis pointing along the thickness of the slice, and the x-axis along the slice width. The boundaries between cortical layers were then defined as x-y planes with constant depth z l . Layer 1 was aneuronal, and layers 2 and 3 were combined. The total model size was then specified by the cortical depth z max , the thickness of the slice y max , the width of the slice x max , and the neuronal density D, with the total number of neurons calculated as N = x max 9 y max 9 z max 9 D. The model slice had dimensions x max = 4.4 mm, y max = 0.4 mm and z max = 2.6 mm, and D = 38,335 neurons/mm 3 , resulting in a model size of 175,421 neurons. We then positioned neurons by placing their somas at random x, y and z values constrained by x max , y max and the z l boundaries of the containing layer, and rotating them by random angles. Pyramidal cells had their apical dendrites aligned parallel to the z-axis.
All neurons could form connections within their group and with neurons from all other neuron groups, according to the values given in Table ESM3. Connections were also constrained by the axonal arborisation radii of the presynaptic neurons, taken from Blasdel et al. (1985) and Fitzpatrick et al. (1985) as adapted by Izhikevich and Edelman (2008). Arborisations were considered in 2D: on the x-y plane on a per-layer basis. We assumed an isotropic Gaussian spatial distribution of connections centred on the presynaptic neuron (Hellwig 2000), setting the arborisation radius equal to two standard deviations of the Gaussian kernel, so that *91 % of connections were contained within the specified arborisation radius. Arborisation radii are given in Table ESM4. When deciding on the targets of a pre-synaptic neuron i in layer l, we calculated the expected number of connections made by i in l remaining inside the slice by multiplying the number of connections specified in Table ESM3 by the ratio f li , defined as the integral of the kernel within the slice boundaries: where a xi is the distance from neuron i to the left edge of the slice, b xi the distance to the right edge, a yi the distance to the front edge, b yi the distance to the back edge, r li is half the arborisation radius of i in layer l (see Fig. 5), and erf the Gaussian error function (this solution is valid provided that a xi and a yi are negative, and b xi and b yi are positive). Pyramidal neuron dendrites span several layers above their soma layer, so we used the connectivity statistics provided per layer for pyramidal neurons by Binzegger et al. (2004). Axonal transmission delays were calculated as the Euclidean distance between the presynaptic and postsynaptic neurons' somas divided by the axonal transmission speed of 0.3 m/s (Hirsch and Gilbert 1991), plus a constant synaptic delay of 0.5 ms to account for the time taken for neurotransmitter release and binding (Katz and Miledi 1965). All these modelling decisions are handled by the initialisation functions in VERTEX, which also allow models to be initialised with uniform spatial connectivity profiles, no synapse reduction, arbitrary delay times, and in a cylindrical shape rather than the cuboid of our slice model. Synapse weights are specified in Table ESM5. We included AMPA and GABA A type conductance-based synapses (the minimal set of synapse types required for generating gamma). When a neuron fired a spike, the synaptic conductance at the contacted target compartments increased by the relevant synaptic weight after the relevant axonal delay time, then decayed exponentially. VERTEX currently includes current-based and conductance-based models of single-exponential and alpha synapses. We stimulated our model to mimic the bath application of kainate, which excites the pyramidal axonal plexus, providing the neurons with excitatory drive. We simulated this by applying independent random input currents to each neuron, modelled as Ornstein-Uhlenbeck processes (similar to Arsiero et al. 2007). Input current parameters are given in Table ESM6. VERTEX can provide random inputs to neurons as either currents or membrane conductance fluctuations.

In vitro experimental methods
All experiments were carried out in accordance with the European Communities Council Directive 1986 (86/609/ EEC), the US National Institutes of Health Guidelines for the Care and Use of Animals for Experimental Procedures, and the UK Animals Scientific Procedures Act.

Surgical preparation
The monkey (Macaca mulatta, male, 8 years old) used in this study was subject to experiments in vivo involving extracellular recording of neural activity and local drug application (iontophoresis). All tissue samples used in this study were taken from intact brain areas that were not the subject of studies performed before tissue extraction. Extraction was performed under general anaesthesia, which was maintained over the course of 4 days. For the anaesthesia, the animal was initially sedated with a 0.1 ml/kg ketamine intra-muscular injection (100 mg/ml). Thereafter, bolus injections of propofol were administered intravenously to allow for tracheotomy and placement of catheters for measuring intra-arterial and central venous blood pressure. During surgery, anaesthesia was maintained by gaseous anaesthetic (2.5-3.9 % sevoflurane) combined with continuous intravenous application of an opioid analgesic (Alfentanil, 120 lg/kg/h), a glucocorticoid (Methylprednisolone, 5.4 mg/kg/h) and saline (50 ml/h). The animal's rectal temperature, heart rate, blood oxygenation and expired CO 2 were monitored continuously during anaesthesia.

Slice preparation
Macaque neocortical samples were routinely obtained from the inferior temporal gyrus. This was confirmed by post hoc anatomical examination of the fixed (paraformaldehyde) whole brain. Following resection, cortical samples were immediately placed in ice-cold sucrose artificial cerebrospinal fluid (ACSF) containing: 252 mM sucrose, 3 mM KCl, 1.25 mM NaH 2 PO 4 , 2 mM MgSO 4 , 2 mM CaCl 2 , 24 mM NaHCO 3 , and 10 mM glucose. Neocortical slices containing all layers were cut at 450 lm (Microm HM 650 V), incubated at room temperature for 20-30 min, then transferred to a standard interface recording chamber at 34-36°C perfused with oxygenated ACSF containing: 126 mM NaCl, 3 mM KCl, 1.25 mM NaH 2 PO 4 , 1 mM MgSO 4 , 1.2 mM CaCl 2 , 24 mM NaHCO 3 , and 10 mM glucose. Persistent gamma frequency oscillations were induced by the application of kainate (400-800 nM) to the circulating ACSF and were deemed stable if there was no change to frequency or power after 1 h. In general, we did not observe spontaneous network activity in the slices before the bath addition of kainate. LFP recordings were taken using multichannel 10 9 10 silicon electrodes with an inter-electrode distance of 400 lm (Utah array, Blackrock Microsystems, Salt Lake City, UT, USA). Time series were digitally sampled at 10 kHz.

Data processing and analysis
Data processing and analysis was performed in Matlab R2012b. We used the same processing chain for both simulated and experimental recordings, except that common average re-referencing, line noise removal and renormalisation were only applied to the experimental recordings. For LFP analysis, recordings were first re-referenced to the common average, then resampled at 1 kHz. We removed line noise and harmonics by band-pass filtering each recording at [49][50][51][99][100][101][149][150][151][199][200][201] Hz and 249-251 Hz (symmetrical Butterworth filter, 8th order) and subtracting the resulting signal from the original signal. The recordings were then band-pass filtered between 2 and 300 Hz (symmetrical FIR filter, Kaiser window, 2,000th order). We restricted our analysis to an 18 s segment of the recording that was identified as artefact-free in all channels by visual inspection of the filtered traces. After filtering, these segments were normalised to zero mean, unit standard deviation to facilitate signal comparison across the MEA.
Power spectra were calculated using the Thomson multitaper method with a time-bandwidth product of 10 (19 tapers) for experimental recordings and 3 (5 tapers) for the shorter simulated recordings, with estimated 95 % confidence intervals calculated using a Chi-squared approach. Total gamma power at each electrode was calculated by taking the integral of the power spectrum between 20 and 40 Hz. Gamma power between electrodes was estimated by bicubic interpolation between electrode locations.