Computing the Local Field Potential (LFP) from Integrate-and-Fire Network Models — TEI-tagged view

Alberto MazzonipersName, Henrik LindénpersName, Hermann CuntzpersName, Anders LansnerpersName, Stefano PanzeripersName, Gaute T. EinevollpersName (2015date). “Computing the Local Field Potential (LFP) from Integrate-and-Fire Network Models.” PLOS Computational Biologytitle·j.

Object of study — proxies for the local field potentialterm·dom computed from point-neuron LIFterm·mth network simulationsterm·mth

Research question — how best to use LIFterm·mth network models to model and understand LFPterm·dom recordings

Contribution — a fixed linear combination of LIFterm·mth synaptic currentsterm·dom that accurately predicts the LFPterm·dom time course

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Computing the Local Field Potentialterm·dom (LFPterm·dom) from Integrate-and-Fire Network Models Abstract Leaky integrate-and-fireterm·mth (LIFterm·mth) network models are commonly used to study how the spiking dynamics of neural networks changes with stimuli, tasks or dynamic network states. However, neurophysiological studies in vivo often rather measure the mass activity of neuronal microcircuits with the local field potentialterm·dom (LFPterm·dom). Given that LFPsterm·dom are generated by spatially separated currents across the neuronal membrane, they cannot be computed directly from quantities defined in models of point-like LIFterm·mth neurons. Here, we explore the best approximation for predicting the LFPterm·dom based on standard output from point-neuron LIFterm·mth networks. To search for this best “LFP proxyterm·key”, we compared LFPterm·dom predictions from candidate proxies based on LIFterm·mth network output (e.g, firing ratesterm·key, membrane potentialsterm·key, synaptic currentsterm·dom) with “ground-truth” LFPterm·dom obtained when the LIFterm·mth network synaptic input currents were injected into an analogous three-dimensional (3D) network model of multi-compartmental neurons with realistic morphology, spatial distributions of somata and synapses. We found that a specific fixed linear combination of the LIFterm·mth synaptic currentsterm·dom provided an accurate LFP proxyterm·key, accounting for most of the variance of the LFPterm·dom time course observed in the 3D network for all recording locations. This proxy performed well over a broad set of conditions, including substantial variations of the neuronal morphologies. Our results provide a simple formula for estimating the time course of the LFPterm·dom from LIFterm·mth network simulationsterm·mth in cases where a single pyramidalterm·key population dominates the LFPterm·dom generation, and thereby facilitate quantitative comparison between computational models and experimental LFPterm·dom recordings in vivo .

Introduction Models of recurrently connected networks of leaky integrate-and-fireterm·mth (LIFterm·mth) neurons are well established tools for studying brain function [ 1 , 2 ]. The equations describing the single LIFterm·mth neuron are simple and can be easily adapted to generate complex dynamics [ 3 , 4 ]. Despite their simplicity, LIFterm·mth network models have proved able to describe a wide spectrum of different cortical dynamics and cortical functions, from the emergence of up and down states [ 5 – 7 ], working memory [ 8 – 10 ], attention [ 11 , 12 ], decision making [ 13 ], rhythmogenesis [ 14 ], and sensory information coding [ 11 , 15 , 16 ]. In some cases it is possible to describe the dynamics of LIFterm·mth networks analytically [ 17 , 18 ], thus providing deeper insights into how spiking neuronal networks may implement the basic cerebral computational mechanisms [ 19 ].

Models can only be properly tested against experimental evidence when they can predict empirical measures quantitatively. Local cortical activity is often recorded in vivo or in vitro using the local field potentialterm·dom (LFPterm·dom), a measure obtained by low-pass filtering (below a few hundred hertz) the electrical potential recorded from extracellular electrodes. The LFPterm·dom signal reflects mass neural activity arising within a few hundred micrometers or more from the recording electrode [ 20 – 25 ]. This spatial scale is highly relevant for LIFterm·mth network models, which typically aim to describe the activity of thousands or tens of thousands of cells. The recording of LFPsterm·dom has a prominent role in systems neuroscience, and such recordings have been used extensively to investigate cortical network mechanisms involved in sensory processing [ 26 ], motor planning [ 27 ], and higher cognitive processes [ 28 ].

LFPterm·dom is generated by transmembrane currentsterm·dom in the neurons in the vicinity of the recording electrode [ 23 ] and depends on morphological features of the contributing cells, the positioning of synapses, as well as the correlation level of synaptic inputs [ 20 , 21 , 29 , 30 ]. Under reasonable assumptions about the extracellular milieu the cellular LFPterm·dom contributions can be computed as a weighted sum of the transmembrane currentsterm·dom in multi-compartmentterm·key neuron models [ 31 – 34 ]. This allows for detailed numerical investigations of spatial, as well as spectral features of the LFPterm·dom signals [ 35 ]. In particular, such simulationsterm·mth of large populations of morphologically detailed neurons have provided insight into how the neuronal activity at the population level influences the spatial reach and laminar variation of the LFPterm·dom signal in vivo [ 20 , 21 , 33 , 34 , 36 ] the relative importance of active and passive currents [ 37 ], and the population LFPterm·dom signal measured from cortical slices in microelectrode arrays (MEAs) [ 38 , 39 ].

However, it has been unclear how best to use LIFterm·mth networks to model and provide understanding of LFPterm·dom recordings. This is because extracellular potentials arise in biological tissue due to a spatial separation of inward (sinks) and outward (sources) transmembrane currentsterm·dom of the neurons, and neuron models used to compute an LFPterm·dom signal must thus have a minimum of two spatially separated compartments in order to generate a potential [ 32 ]. In LIFterm·mth models, however, a single compartment is typically used as an approximation of the entire neuron, including the spatially extended dendritic structure, and individual cells within a population are not assigned to a specific spatial position.

One possible way to compute LFPsterm·dom from LIFterm·mth network is to project the spiketerm·dom times generated by the LIFterm·mth network under consideration onto morphologically detailed 3D neuron models and then compute the field that the currents flowing through these 3D networks generate. However, this approach would require the modeler to set up a cumbersome and computationally expensive network model based on multi-compartmentterm·key model neuron. As a much simpler alternative, we here instead search for a general and easy-to-use proxy to predict the time course of the LFPterm·dom based on variables available directly from the LIFterm·mth network simulationsterm·mth. Here we investigate and evaluate different strategies to compute an LFP proxyterm·key directly from the output of standard LIFterm·mth network simulationsterm·mth without the use of multi-compartmentterm·key neuronal morphologies.

Our approach is as follows: we first simulate an LIFterm·mth point-neuron network model and record the output spiking activity, membrane potentialsterm·key, and synaptic currentsterm·dom. Next, we compute a realistic ground-truth estimate of the LFPterm·dom that the same LIFterm·mth network activity would generate. We do this by injecting distributed synaptic currentsterm·dom corresponding to the stored LIFterm·mth synaptic events, onto a population of multi-compartmentterm·key neurons with realistic distributions of dendrites and synapses (we call this population the “3D network”). We then compare this simulated ground-truth LFPterm·dom signal to a number of LFPterm·dom proxies computed directly from measures of activity of the point-neuron LIFterm·mth network. These proxies include those previously proposed in the literature (e.g., the average firing rateterm·key [ 11 , 14 , 40 ], the average membrane potentialterm·key [ 24 , 41 – 44 ], the sum of synaptic currentsterm·dom [ 7 , 45 ], and the sum of absolute values of synaptic currentsterm·dom [ 15 ]), as well as others proposed here. By separating the spiking dynamics generated by the LIFterm·mth network from the LFPterm·dom generated by the 3D network, we are also able to investigate how different assumptions regarding cell morphology, synaptic distributions and recording positions influence the accuracy of the different LFPterm·dom proxies.

We find that a simple linear combination of excitatory (AMPAterm·key) and inhibitory (GABAterm·key) synaptic currentsterm·dom extracted from the point-neuron LIFterm·mth network provides a proxy for the LFPterm·dom that closely matches the temporal features of the signals resulting from the morphologically realistic LFPterm·dom model generated by the 3D network. Even with a small set of fixed parameters this LFP proxyterm·key is able to account for the LFPterm·dom signal with a high degree of precision under most investigated conditions.

Results A morphological model for a cortical LIFterm·mth network Our goal was to understand how to compute a simple yet accurate approximation (denoted as “proxy” in the following) of the LFPterm·dom that would be generated by the time series of synaptic activity of an LIFterm·mth network if its neurons had a realistic spatial structure and arrangement. We therefore first simulated an LIFterm·mth network (known to reproduce several features of cortical dynamics). Next, we injected the synaptic activity it generated into a synthetic three-dimensional network model (3D network) of a layer of a cortical column that employed multi-compartmental neurons with realistic morphology, spatial distributions of somata and synapses, and computed the extracellular potentials generated by this synaptic activity.

We selected an LIFterm·mth network (adapted from [ 14 ] and refined in [ 15 , 16 , 46 , 47 ]) that has been shown to reproduce a number of important features of the dynamics of visual primary cortical neural population recorded in vivo during naturalistic sensory stimulation, including a realistic spectrum of cortical dynamics and of its modulation with the visual stimuli, including low-frequency (1–12 Hz) and gamma (50–100 Hz) oscillations [ 15 , 46 ]. Moreover, when using a simple proxy (which is demonstrated below to perform well) to compute an LFPterm·dom from synaptic currentsterm·dom, this LIFterm·mth network reproduced quantitatively several important properties of recorded extracellular potentials, including LFPterm·dom power spectra and spectral information content [ 15 ], and cross-frequency and spike-field relationships [ 16 , 46 ]. Thus, the LIFterm·mth network seemed to generate a sufficiently realistic dynamics to provide synaptic input for the generation of biologically plausible LFPsterm·dom in the 3D network.

The LIFterm·mth network model ( Fig 1A ) was composed of 4000 excitatory and 1000 inhibitory LIFterm·mth neurons that were randomly connected with a pair-wise connection probability of 0.2 (for further details see Methods ). The LIFterm·mth network received two kinds of external inputs: a “thalamic” synaptic input thought to carry the information about the external stimuli and a stimulus-unrelated input representing slow ongoing fluctuations of activity. Synaptic dynamics and parameters are reported in Tables 1 and 2 , and further details can be found in the Methods. Importantly, as is the case for most LIFterm·mth network models to date, our LIFterm·mth network did not have any spatial structure: the individual neurons were not assigned to a specific spatial position and consequently the connectivity had a random and sparse structure.

LIFterm·mth network and 3D morphological network. (A) Sketch of the leaky integrate-and-fireterm·mth (LIFterm·mth) network. A population of 1000 interneurons with GABAterm·key synapses (blue) and a population of 4000 pyramidalterm·key neurons with AMPAterm·key synapses (red) receive recurrent inputs (random connectivity with 20% probability) and two kinds of external inputs: global ongoing cortical activity (Ornstein-Uhlenbeck process) and a regular thalamic stimulation. (B) Sketch of the morphological 3D network made of two stacked cylinders with 250 μm radius and 250 μm height. A representative interneuron and pyramidalterm·key cell are depicted. Interneuron dendrites remain in the lower cylinder while the pyramidalterm·key neuron dendrites reach out to the upper cylinder. Dendrites in the lower cylinder receive both AMPAterm·key and GABAterm·key synapses while dendrites in the upper cylinder receive only AMPAterm·key synapses. (C) Graphical rendering of a subset of the 3D network composed of 10 interneurons and 40 pyramidalterm·key neurons. (D) Raster plot of the spiking activity of the 10 interneurons (blue, top) and the 40 pyramidalterm·key neurons (red, bottom) with the highest spiking activity in the LIFterm·mth network, for a thalamic stimulation of 1.5 spikesterm·dom/ms. (E) Depth-resolved LFPterm·dom signal as simulated by injecting the spikesterm·dom generated by the whole network during thalamic stimulation of 1.5 spikesterm·dom/ms into the 3D network. Black lines show LFPterm·dom for 100 and -100 μm depth.

Summary of leaky integrate-and-fireterm·mth (LIFterm·mth) network model.

Parameters for the two cell types used in the LIFterm·mth network model.

The LFPterm·dom signal that would result from the time series of spikesterm·dom generated by the LIFterm·mth network provided the postsynaptic neurons had biologically plausible dendritic structures, was computed by injecting the LIFterm·mth synaptic activity into a 3D network of morphologically detailed multi-compartmental model neurons ( Fig 1B , see Methods ). A summary of the properties of the 3D network is reported in Table 3 , while the synaptic parameters are listed in Table 4 (see Methods for further details). In order to set up the 3D network we were required to make additional assumptions regarding the spatial positioning of cells, the shape and size of their dendritic structures, as well as the synaptic distributions. We focused on computing the LFPterm·dom generated by one cortical layer (in terms of soma positions) that comprised both inhibitory and excitatory neurons. In our default setting, we assumed all neurons in the 3D network to be inside two cylinders with 250 μm radius and 250 μm height that were stacked one above the other to resemble the vertical structure of layer 2/3 ( Fig 1B ). Note that this spatial scale is similar to the size of the neuronal pool contributing to the recorded LFPterm·dom, the so-called spatial reach, in the case of uncorrelated synaptic activity driving a neuronal population [ 20 , 21 , 35 ], and resulted in a neuronal density consistent with known estimates of 50000 neurons per mm 3 in the cortex [ 48 ]. While our two model populations most directly resemble a pair of excitatory and inhibitory populations in cortical L2/3, we show in subsection “Dependency of the LFPterm·dom signal on dendritic morphology” that our results also pertain to the LFPterm·dom generated by neuron morphologies found in other cortical layers.

Summary of 3D network model of multi-compartmental neurons.

Parameter values for the 3D network model of multi-compartmental neurons.

Given these geometrical constraints, we created the multi-compartmental cell models in the 3D network in the following way: soma locations for all cells were homogeneously distributed within the lower cylinder ( Fig 1B ). Next, we placed artificial straight axons that were distributed at random cortical depths and random orientations within both cylinders. They served as targets in an algorithmic generation of dendrites, through which pyramidalterm·key cell dendrites were connected to all axons within a specified reach distance while optimizing the following wiring conditions: short conductions times, short total cable length and synaptic democracy (i.e., equal impact of synaptic inputs at the site of dendritic integration [ 49 , 50 ]). This procedure has been shown previously to reproduce pyramidal-cell-like dendrites [ 51 ]. The number of axons and their length were set so that the resulting cell morphologies matched the membrane surface distribution of real cortical layer 2/3 pyramidalterm·key cell reconstructions [ 52 ] within the constraints of the simplified columnar arrangement that was chosen for this study. This procedure also provided good matches for total cable lengths and number of branch points (compare membrane surface distribution in S1 Fig and see Methods for more details). Note that the virtual axons used for the generation of the morphologies were subsequently discarded. Since the membrane area (and consequently the transmembrane currentterm·dom) of the axons is very small compared to the dendrites, we expected them to have a negligible contribution to the present 3D network LFPterm·dom generation.

Stellate cell dendrites were generated in a similar manner, but were only connected to axons in the lower cylinder. This resulted in stellate cell morphologies with realistic bush-like dendrites. Fig 1C illustrates the overall structure of the resulting 3D network.

To further validate the simulationterm·mth results obtained with these morphologies, we also built an alternative 3D network with anatomically reconstructed morphologies (see Methods ) and checked that the results were essentially the same as for algorithmically grown morphologies (see subsection “ Performance of LFPterm·dom proxies in different dynamic network states ”).

Finally, AMPAterm·key synapses were homogeneously distributed over the whole neuronal surface while GABAterm·key synapses were located only in the lower cylinder, closer to the soma ( Fig 1B , see Methods for details). Alternative synaptic distributions are explored in the “Dependency of the LFPterm·dom signal on the distribution of synapses” subsection.

For each neuron in the LIFterm·mth network we randomly assigned a multi-compartmental neuron model with a unique dendritic structure in the 3D network. The connectivity of the LIFterm·mth network determined which postsynaptic spikesterm·dom in the LIFterm·mth network simulationterm·mth should serve as input spikesterm·dom for each multi-compartmentterm·key neuron. We then used these spiketerm·dom times together with the external input (see above) to activate synaptic currentsterm·dom in the 3D network (see Methods ). In this way we assured that the synaptic input in a multi-compartmentterm·key neuron was identical to its LIFterm·mth neuron counterpart. The synaptic dynamics in the 3D network was identical to that in the LIFterm·mth network. In a subsequent step, we took into account all transmembrane currentsterm·dom in the neurons of the 3D network to compute the LFPterm·dom by means of well-established volume conduction theory and the so called line-source method [ 31 , 34 ] (see Methods ).

Spatial distribution of simulated LFPterm·dom signal Fig 1 shows a half-second excerpt of results for an example simulationterm·mth using the spiking activity generated by the LIFterm·mth network ( Fig 1D ) in response to a 1.5 spikesterm·dom/ms stimulus (see Methods for details) to calculate the corresponding LFPterm·dom signal along the vertical axis of the cylinder at different electrode depths from the 3D network ( Fig 1E ). The temporal fluctuations of the LIFterm·mth signal were strongly correlated across depth, albeit with a sign shift around the depth just between the two cylinders (which we from now on will refer to as the inversion point). The sign of the baseline (DC) LFPterm·dom was negative above the inversion point while it was positive below it. This reflects that the LFPterm·dom was dominated by the perisomatic inhibitory synapses generating a net source current close to the soma and sink return currents in the apical branches. The excitatory synapses contributed less due to their homogeneous distribution ( Fig 1B ), giving only a weak current dipole [ 29 ], as will be discussed in more detail in the next sections.

We defined the amplitude of LFPterm·dom fluctuations at each depth as the standard deviation of the signal over time, and further assigned it the same sign as the LFPterm·dom baseline, i.e., negative/positive above/below the inversion point. The magnitude of the LFPterm·dom amplitude was largest around the middle of each cylinder ( Fig 2A ), decreased steeply close to the inversion point and more smoothly beyond the vertical boundaries of the network. The decrease of the amplitude of the LFPterm·dom fluctuations when the electrode was moved away from the center of the 3D network is shown in Fig 2B for all depths. This decrease in LFPterm·dom power was consistent with results of [ 20 ]: inside the 3D network ( X / R < 1, where X is the displacement of the electrode from the center and R is the radius of the cylinder) differences were small, but when the electrode was placed outside the 3D network ( X / R > 1) the decrease was steep. Note that the region around the inversion point where the potential is very small, broadened with the distance from the center.

Simulated local field potentialsterm·dom (LFPsterm·dom) as a function of depth and lateral position of the electrode from the center of the 3D network. (A) Amplitude of LFPterm·dom signal generated by the morphological 3D network at 50 μm spaced depths, when the dynamics were driven by a 10 seconds thalamic input of 1.5 spikesterm·dom/ms. The amplitude was measured as the standard deviation of the LFPterm·dom signal over the entire time course with the same sign as the baseline (see Fig 1E ). Dashed lines indicate network boundaries (-250 μm < d < 250 μm) (B) Amplitude of LFPterm·dom signal at different depths and distances from the center of the 3D network. Distances were measured in units of 3D network radius R (= 250 μm). The dashed lines separate the area inside the network ( X / R ≤1 & -250 μm < d < 250 μm) and outside the network. (C) Power spectral density (PSD) of LFPterm·dom signal in the center of the 3D network for different depths. Dashed lines indicate network boundaries -250 μm < d < 250 μm). (D) Power spectral density (PSD) of LFPterm·dom signal at different distances from the center of the 3D network at a reference depth of 100 μm. Dashed line indicates network boundary ( X / R = 1).

We observed that all power spectra recorded outside this noise-dominated region had similar shapes ( Fig 2C and 2D ), suggesting that LFPterm·dom fluctuations could be roughly approximated by the same time series rescaled by the numerical value of the LFPterm·dom amplitude shown in Fig 2B . The observation that LFPsterm·dom recorded in different spatial positions had similar temporal behavior and differed mainly by a scaling factor, suggested that a single LFP proxyterm·key could work for recordings at different depths and positions in the horizontal plane, provided that it is properly scaled. Such a factorization of spatial and temporal dimensions can be expressed (see [ 30 ]) as L F P p r o x y ( r , d , t ) = f p r o x y ( r , d ) * g p r o x y ( t ) (1) where d is the depth and r the distance from the population center. The term f proxy ( r , d ) then gives the amplitude of the signal as a function of the electrode position (as in Fig 2B ) while the dimensionless g proxy ( t ) has variance equal to one and describes the temporal features of the LFPterm·dom signal.

We first focused on finding the optimal g proxy ( t ) for an LFPterm·dom signal recorded at selected depths along the central vertical axis ( X/R = 0 ) of the 3D network. However, we found (see subsection “ New class of LFPterm·dom proxies ”) that the identified optimal LFP proxyterm·key was applicable also to other depths and radial distances of the populations (given an appropriate overall scaling of the signal amplitudes, cf. Fig 2B ).

The contribution to the LFPterm·dom signal from synaptic inputs onto the interneurons (and their associated return currents) was negligible both in amplitude ( Fig 3A ) and in determining the LFPterm·dom spectrum ( Fig 3B ). This was due to the different morphologies of the two types of neurons: consistently with what was shown previously for stellate cells with symmetrically placed synapses [ 20 ] (i.e., a so-called close-field arrangement [ 53 ]), the contribution from the interneurons to the LFPterm·dom was negligible ( Fig 3 ). Further, the associated power spectrum of this contribution was closer to colored noise and did not display gamma fluctuations. We investigate this in detail in the subsection “Dependency of the LFPterm·dom signal on dendritic morphology” .

Contribution from individual neuron types to simulated LFPterm·dom signal. Decomposition of LFPterm·dom obtained in same conditions as Fig 2 into contributions from currents through the membrane of interneurons and pyramidalterm·key neurons. (A) Depth-resolved amplitude of LFPterm·dom signal generated by all neurons (black), by the pyramidalterm·key neurons (red), and by the interneurons (blue). (B) Corresponding LFPterm·dom power spectra for the three sets depicted in (A) at a depth of 100 μm.

Since we obtained a very similar LFPterm·dom when we only simulated the contribution from synaptic inputs onto the pyramidalterm·key neurons, all the results shown in the following will, unless otherwise stated, consider only the contributions from pyramidalterm·key neurons to LFPterm·dom. Likewise, the LFPterm·dom proxies will be based only on input onto excitatory neurons (as done previously [ 15 ]). However, the inhibitory neurons obviously play a key role (i) in generating the dynamics and (ii) in providing the GABAterm·key currents of synapses onto pyramidalterm·key neurons that contribute strongly to the LFPterm·dom.

Performance of different LFPterm·dom proxies We first tested six LFP proxyterm·key candidates ( Fig 4A ): AMPAterm·key currents, GABAterm·key currents, the average firing rateterm·key FR, the average membrane potentialterm·key V m , and the sum of these synaptic currentsterm·dom ∑ I as well as their absolute values ∑| I |. Note that the "AMPAterm·key currents" and "GABAterm·key currents" proxies are defined as the sum of the post-synaptic currents for each type of synapse over all pyramidalterm·key neurons (see Table 1G ). These currents have depolarizing and hyperpolarizing effects, respectively, on the postsynaptic neurons. We thus here use the convention that assigns a positive sign to AMPAterm·key currents and a negative sign to GABAterm·key currents.

Performance of candidate LFPterm·dom proxies. (A) Illustrations of predictions of LFPterm·dom time courses from candidate LFPterm·dom proxies. From top to bottom: firing rateterm·key (FR), membrane potentialterm·key (V m ), AMPAterm·key currents, GABAterm·key currents (note: these have a negative sign), sum of absolute values of AMPAterm·key and GABAterm·key currents ∑| I |, sum of AMPAterm·key and GABAterm·key currents ∑ I . Results are shown for a thalamic stimulation of 1.5 spikesterm·dom/ms, and the proxies are normalized to have variance equal to one (see text),. (B-C) time course of the LFPterm·dom signal (black) for reference depths 100 μm (B) and -100 μm (C) compared to the best matching proxy, ∑| I | (magenta). The title indicates the fraction of variance explained (85% in both cases). (D-E) Cross-correlation in time between the LFPterm·dom and the ∑| I | proxy for the two depths. Note that the peaks corresponding to the highest cross-correlation magnitudes corresponded to a lag of 1 ms, i.e., the LFPterm·dom was best predicted by the value of ∑| I | one millisecond in the past. (F-G) Fraction of LFPterm·dom signal variance explained by different LFPterm·dom proxies with optimal delay (same color code as (A)). as a function of depth. The sum of absolute values of the synaptic currentsterm·dom ∑| I | was the best proxy, followed by the use of GABAterm·key alone. The firing rateterm·key FR was a poor proxy, and the other three were moderately good proxies.

Because of the opposite signs assigned to the AMPAterm·key and GABAterm·key currents, the sum of the absolute values of the currents ∑| I | is equivalent to the difference between the currents. For several reasons, i.e., synaptic delay and dendritic filtering, we expected the best proxy for the LFPterm·dom time course to possibly involve time-delayed measures of LIFterm·mth network variables. To assess the best values of these delays we first computed the cross-correlation function between the ground-truth LFPterm·dom and the considered LFP proxyterm·key obtained from the LIFterm·mth network, and found the delay at which the absolute value of the correlation was largest (for half of the recording depths the correlation is negative due to LFPterm·dom inversion). The LFP proxyterm·key that we chose was the z-scored (i.e., baseline-subtracted and normalized to have variance equal to one) and time-shifted LIFterm·mth network variable that maximized the fraction of variance explained, R 2 . Finding the best delay and rescaling factor was done separately for each depth, but we found that the differences in the observed best values of the delay across depth, were minor (see S2B Fig ). Fig 4B and 4C shows the comparison between the 3D network LFPterm·dom signal at two different electrode depths and the LFP proxyterm·key given by the sum of absolute values of the synaptic currentsterm·dom ∑| I |, that given our sign convention simply becomes the difference between the currents, i.e., L F P ∑ | I | ( r , d , t ) = f ∑ | I | ( r , d ) * N o r m [ ∑ p y r A M P A ( t − τ ) − ∑ p y r G A B A ( t − τ ) ] (2) where Norm [] indicates the mean-subtracted, normalized version of the time series between square brackets. Fig 4D and 4E shows the cross correlation between the 3D network LFPterm·dom signal and proxy for the two depths.

A comparison of the average fraction of variance explained by all the LFPterm·dom proxies displayed in Fig 4A across different depths ( Fig 4F and 4G ) shows that the best one was the sum of absolute values of synaptic currentsterm·dom ∑| I | (<R 2 > = 0.83) followed by the negative of the GABAterm·key currents (<R 2 > = 0.81) and then the AMPAterm·key currents (<R 2 > = 0.78). The negative of the sum of synaptic currentsterm·dom ∑ I and membrane potentialterm·key V m performed in a similar way (<R 2 > = 0.69), while the firing rateterm·key FR gave a poor fit (<R 2 > = 0.51). The R 2 is slightly larger for depths about 100 μm from the inversion point, probably due to stronger synaptic and return currents.

We found two results to be of particular interest. The first was that a proxy based on GABAterm·key currents alone gave clearly a better match for the simulated LFPterm·dom signal than the AMPAterm·key currents alone. The second was that the ∑| I | gives the best fit which suggests that the magnitude of the AMPAterm·key currents locally sums with the magnitude of the GABAterm·key return currents. Thus the two types of synaptic currentsterm·dom contribute to the LFPterm·dom with the same sign. This feature is partly due to the fact that AMPAterm·key synapses are distributed over the whole surface of pyramidalterm·key neurons, while GABAterm·key synapses are located only in the lower cylinder close to the soma ( Fig 1B ). This will be further investigated in the “Dependency of the LFPterm·dom signal on the distribution of synapses” subsection.

The fits above were computed by averaging the time-varying variables over the set of excitatory neurons in the LIFterm·mth network. However, we also tested the quality of the fit obtained by averaging over all the neurons in the LIFterm·mth network or only over inhibitory neurons. The results for each variable and depth are shown in S2A Fig , together with the associated optimal delays. The relative ranking of the candidate proxies remains unaltered. Further, proxies obtained by averaging the firing rateterm·key, the membrane potentialterm·key, or the synaptic input currents over the excitatory neurons (as above) performed better than proxies obtained by averaging the same variables over the inhibitory neurons set and roughly the same as proxies obtained averaging over all neurons ( S2A Fig ).

New class of LFPterm·dom proxies Since AMPAterm·key and GABAterm·key currents contributed differently to the LFPterm·dom signal we investigated a novel proxy, the weighted sum between AMPAterm·key and GABAterm·key currents (WS), that uses a linear combination of AMPAterm·key and GABAterm·key synaptic currentsterm·dom where we introduce a factor α describing the relative contribution of the two currents and a specific delay for each type of current: L F P W S ( r , d , t ) = f W S ( r , d ) * N o r m [ ∑ p y r A M P A ( t − τ A M P A ) − α ( ∑ p y r G A B A ( t − τ G A B A ) ) ] (3) Note that the two proxies ∑| I | and ∑ I are particular cases of the above equation in which the delays are the same, and α is equal to 1 and -1 respectively.

We first tested the WS proxy with the electrode located in the center of the 3D network for different depths. The optimal value of α was always positive, but varied across depths ( Fig 5A ). The optimal delays were always in the range [ 5 – 7 ] ms for τ AMPAterm·key ms and in the range [-1 1] ms for τ GABAterm·key . This implies that the optimal LFP proxyterm·key was achieved by subtracting the GABAterm·key PSCs (postsynaptic currentsterm·dom) from the AMPAterm·key PSCs occurring around 6 ms in the past. Performance was very high for all depths (up to 93% of variance explained, see Fig 5B ). Since the optimal values of α , ( Fig 5A ) τ AMPAterm·key and τ GABAterm·key ( S2B Fig ) were relatively stable across depths, we defined a new proxy: the reference weighted sum LFP proxyterm·mth (RWSterm·mth). The structure of the RWSterm·mth proxy is the same as the WS proxy but the variables are fixed: α is set to the average accross depths of the optimal values for WS (1.65, see Fig 5A ) and the delays to τ AMPAterm·key = 6 ms and τ GABAterm·key = 0 ( S2C Fig ). This results in L F P R W S ( r , d , t ) = f R W S ( r , d ) * N o r m [ ∑ p y r A M P A ( t − 6 m s ) − 1.65 ( ∑ p y r G A B A ( t ) ) ] (4) New proxy explaining more than 90% of the variance in the LFPterm·dom signal. (A) Value of relative contribution of AMPAterm·key and GABAterm·key currents (α parameter in Eq ( 3 )) optimizing correlation between WS proxy and ground truth LFPterm·dom in the same conditions as Fig 4 . Dashed line indicates average value over depths α = 1.65 used for RWSterm·mth proxy ( Eq 4 ). (B) Same as Fig 4F and 4G including also WS proxy (black) and RWSterm·mth proxy (blue). (C) Fraction of LFPterm·dom variance explained, averaged over all depths, by different proxies for different distances from the center of the 3D network. From best to worst: weighted sum (WS), reference weighted sum (RWSterm·mth), sum of absolute values of the synaptic currentsterm·dom (∑| I |), membrane potentialterm·key (V m ), sum of synaptic currentsterm·dom (∑ I ), firing rateterm·key (FR). Error bars are not displayed since they would not be visible in the figure. (D) Mean and standard deviation across depths of optimal coefficients of α in the WS proxy as a function of electrode lateral position. Dashed line indicates the fixed coefficient of the RWSterm·mth proxy that was not optimized but kept constant for all depths and distances.

We found that the performance of this proxy was almost indistinguishable from the single-depth optimized values across depths ( Fig 5B ) and largely outperformed all other proxies. Moreover, we found the performance of a proxy with α = 1.65 to be very good (>80% of variance explained) for a broad range of other AMPA- and GABA-current delays ( S2C Fig ).

We next tested the performance of the proxies for different distances of the electrode from the center of the 3D network: Fig 5C compares the fraction of variance explained by WS, RWSterm·mth and the other proxies mentioned above for LFPsterm·dom measured at different distances from the center of the 3D network. The depicted results are found from averaging across all depths. The Standard Error of the Mean of R 2 across depths was <1% for all proxies and all lateral displacements and is not displayed in the figure since it would not be visible. Values for explained variance were very stable for different lateral electrode positions: in particular, for all lateral displacements RWSterm·mth performances were similar to WS and outperformed all other proxies ( Fig 5C ). The average optimal value of α across depths was always close to the reference value 1.65 ( Fig 5D ). Given that the RWSterm·mth proxy was much simpler than WS (see below) and able to explain more than 90% of the variance of the LFPterm·dom time course at a wide range of electrode recording positions, we tentatively propose this as the best proxy for the LFPterm·dom signal computable directly from LIFterm·mth network variables.

The proxies given by the combination of two synaptic parameters (WS and RWSterm·mth) have four free parameters (scale as described by the function f in Eq 1 and following, AMPAterm·key and GABAterm·key delays, relative amplitude of AMPAterm·key and GABAterm·key contribution) while the other proxies have only two free parameters (scale, delay). We assessed by means of the Bayesian Information Criterion (BIC, [ 54 ], see Methods for details) whether the benefit in terms of improved performance of the models based on the linear combinations of synaptic currentsterm·dom was worth the increase in model complexity due to the higher number of parameters. We found that, according to this model selection criterion, RWSterm·mth outperforms all previous proxies and WS outperforms RWSterm·mth and all other proxies ( S3 Fig ), demonstrating the power of the RWSterm·mth and WS models. However, the optimal WS parameters are by construction different for each recording position and, as we will see in the following, for various network structures and states. Thus comparison of LFPterm·dom predictions from use of the WS proxy requires detailed knowledge about recording position as well as the characteristics of the underlying network, and will thus have limited practical use. On the other hand, as the parameters of the RWSterm·mth proxy are fixed, it can be used directly for all locations in space. As seen in the following, the RWSterm·mth proxy performs well for a broad set of conditions (input intensity, neuron morphology, synaptic distribution), and this means crucially that the proxy can be used also under weak assumptions about the spatial structure of the underlying network. We thus conclude that RWSterm·mth is the best LFP proxyterm·key based on LIFterm·mth network variables. In the following we will test its robustness for different dynamic network states, spatial architectures and synaptic properties.

Performance of LFPterm·dom proxies in different dynamic network states So far we investigated the LFPterm·dom proxies using LIFterm·mth networks in a state exhibiting weakly synchronized oscillations in the spiking dynamics, stimulating the LIFterm·mth network at a relatively low intensity (1.5 spikesterm·dom/ms). However, LIFterm·mth networks can generate a variety of different dynamic network states when the frequency of external inputs is varied [ 14 , 15 , 18 ]. In order to test LFPterm·dom proxies in different dynamic network states, we stimulated the LIFterm·mth network with a wide range of input intensities, covering both much higher and much lower intensities than the one tested in above.

Fig 6A shows, from left to right, a raster plot of a subset of neurons in the LIFterm·mth network for a low-intensity input (0.5 spikesterm·dom/ms), the default input level (1.5 spikesterm·dom/ms), and a high-intensity input (6 spikesterm·dom/ms). Shown below ( Fig 6B ) is the LFPterm·dom signal generated in the 3D network at the reference depth of 100 μm for these three cases together with their corresponding WS fits. For external stimulation with 0.5 spikesterm·dom/ms, recurrent activity in the LIFterm·mth network was almost absent, with all pyramidalterm·key neurons and most interneurons being silent. The LFPterm·dom amplitude was very small and the signal very noisy. For an input of 1.5 spikesterm·dom/ms, firing was sparse with coexisting slow and high-frequency LFPterm·dom fluctuations, and for 6 spikesterm·dom/ms the dynamics were dominated by high-frequency LFPterm·dom gamma oscillations also visible in the LIFterm·mth network spiking activity.

Effects of dynamic network states of the LIFterm·mth model on the simulated LFPterm·dom signal. (A) Raster plots of 50 interneurons (blue, top) and 200 pyramidalterm·key neurons (red, bottom). Neurons displayed are those with the highest number of spikesterm·dom fired in the considered interval. Each panel corresponds to a different stimulation frequency: from left to right: 0.5 spikesterm·dom/ms, 1.5 spikesterm·dom/ms (the stimulation used in Figs 1 – 4 ), 6 spikesterm·dom/ms. Note that in the selected interval all pyramidalterm·key neurons and most interneurons were silent for 0.5 spikesterm·dom/ms stimulation. (B) LFPterm·dom signal (black line) for 100 μm depth and corresponding best fit with the WS proxy (red) for these three stimulation frequencies. The titles show the fraction of variance explained over the whole 10 second simulationterm·mth period. Note the different vertical scales. (C) Average fraction of LFPterm·dom variance explained over all depths by different proxies for different thalamic input frequencies. Error bars are not displayed since they would not be visible in the figure. Same proxy arrangement as Fig 5C . (D) Mean and standard deviation across depths of optimal coefficients of α in the WS proxy as a function of thalamic input. Dashed line indicates the fixed coefficient of the RWSterm·mth proxy.

With an input frequency of 0.5 spikesterm·dom/ms, none of the candidate proxies was able to account for the LFPterm·dom ( Fig 6C ). This was presumably because in these low-firing conditions, randomly occurring, uncorrelated synaptic inputs onto the neurons close to the electrode dominated the LFPterm·dom signal. Such activity does not give a strong dipolar LFPterm·dom pattern [ 21 ] and is apparently more difficult to capture with the global LIFterm·mth network variables considered in the proxies. For the larger inputs ranging from 1 to 6 spikesterm·dom/ms, however, the WS proxy was able to explain more than 91% of the variance. RWSterm·mth was able to explain 88–91% of the variance between 1 to 3 spikesterm·dom/ms with a small decrease to 87% for an input of 6 spikesterm·dom/ms ( Fig 6C ). For inputs of 1 spikesterm·dom/ms or more the sum of the absolute values of the synaptic currentsterm·dom explained 81–85% of the variance, the membrane potentialterm·key 70–79%, the sum of the currents 70–79%, and the firing rateterm·key only 51–60%. Overall, the ranking of the proxies regarding their R 2 values remained the same for all dynamic network states and the RWSterm·mth provided an excellent proxy in all cases. As shown in Fig 6D the relative weighting between AMPAterm·key and GABAterm·key currents as given by the parameter α for the WS proxy was stable and close to the reference values 1.65 chosen for RWSterm·mth for input stimulus intensities, except for the case of very low-intensity input in which the LFPterm·dom signal is almost absent and the fit is poor.

A key property of the LIFterm·mth network is that it exhibits a prominent gamma-band activity (30–100 Hz) in the overall firing activity when the input intensity is increased as indicated by an increased peak in the power spectral density (PSD) [ 15 ]. We therefore investigated how this is reflected in our simulated LFPterm·dom signal and how well the LFPterm·dom proxies capture these properties of the LFPterm·dom signal. Fig 7A shows the power spectra for three different input frequencies. All proxies except for the membrane potentialterm·key tended to underestimate the low frequency LFPterm·dom fluctuations and to overestimate frequencies in the gamma range. WS and RWSterm·mth proxies both produced a nearly perfect fit of the LFPterm·dom spectrum in the gamma-band range while exhibiting the smallest error in the low frequency components among all proxies. In the 1–3 spikesterm·dom/ms input range the modulation of the LFPterm·dom gamma power was well approximated by all proxies, while for 6 spikesterm·dom/ms input, WS and RWSterm·mth underestimated it ( Fig 7B ). All proxies essentially predicted the correct peak LFPterm·dom gamma frequency ( Fig 7C ) for all input levels above 1 spikesterm·dom/ms.

Spectral analysis of LFPterm·dom signal. (A) Power spectra of the LFPterm·dom signal at recorded at the 100 μm depth position, and predictions from candidate LFPterm·dom proxies for input stimulation frequencies of 1 spikesterm·dom/ms (left), 1.5 spikesterm·dom/ms (center), 3 spikesterm·dom/ms (right). Similar results are found for recordings at all depths. Note that for low inputs the power decreased almost monotonously with frequency, while for sufficiently strong input, LFPterm·dom gamma fluctuations appeared and increased in power and frequency with stimulus intensity. (B) Variation of gamma (30–100 Hz) peak power with stimulus input intensity for the LFPterm·dom signal and all proxies (measured as relative increase compared to the power at 1.5 spikesterm·dom/ms). (C) Gamma peak frequency as a function of input frequency for the LFPterm·dom signal and all proxies.

Dependency of the LFPterm·dom signal on dendritic morphology We hypothesized that the negligible contribution of inhibitory neurons was due to the weak dipole moment created by the symmetrically placed synapses on the dendrites of stellate cells [ 29 ]. To test this hypothesis we investigated in the following the effect of neuron shape on the LFPterm·dom generation by systematically altering the morphology of the interneuron population while keeping its inputs fixed. This manipulation also tested the robustness of the LFPterm·dom proxies to the specific choice of the neuronal morphology. We started with two overlapping cylinders (distance = 0 μm) describing the stellate cell morphology. Then we progressively increased their “pyramidalness”, i.e., the distance between the two dendritic bushes and generated a new interneuron population for each cylinder distance ( Fig 8A ; see Methods for details). The generated morphologies ranged from pure stellate cells (the interneuron used in the reference case), to cells corresponding to layer 2/3 pyramidalterm·key cells where the two cylinders were juxtaposed (the pyramidalterm·key neuron used in the reference case), to cells where the two areas were parted by several hundred micrometers (as in layer 5 pyramidalterm·key neurons). In all cases GABAterm·key synapses were distributed only on dendrites located inside the lower cylinder, while AMPAterm·key synapses were distributed over the entire dendritic tree ( Fig 1B ).

Effect of neuronal morphologies on neural signal. (A) Manipulation of the relative position of the two cylinders in which the dendrite arborizes, 50 μm steps. Each relative position induces a specific ‘pyramidalness’ of the model cortical cells (see Methods for details). Red and blue lines indicate the center of the upper and lower cylinder respectively. When the two cylinders are completely superimposed (first panel of top row), the structure corresponds to a stellate cell. When the two cylinders are on top of each other (last panel of top row, first of bottom row), the morphology corresponds to a layer 2/3 pyramidalterm·key cell. When the boundaries of the two cylinders are separated by 250 μm (last panel of bottom row), the cell morphology resembles a layer 5 pyramidalterm·key cell. (B) LFPterm·dom amplitude as a function of depth and distance between the two cylinders. (C) Average absolute amplitude (standard deviation) over depths of LFPterm·dom fluctuations as a function of distance between cylinders. (D) Amplitude of current dipole moment as a function of distance between cylinders. (E) Average fraction of LFPterm·dom variance over all depths explained by different proxies for different distances between cylinders. Error bars are not displayed since they would not be visible in the figure. Same proxy arrangement as in Fig 5C . (F) Mean and standard deviation across depths of optimal coefficients of α in the WS proxy as a function of distance between cylinders. Dashed line indicates the fixed coefficient of the RWSterm·mth proxy. Note that since for distances below 100 μm the performance of the fit was poor (see panel (E)), the fitted value of the relative weight of AMPAterm·key and GABAterm·key currents in contributing to the LFPterm·dom signal has little significance.

We found that the LFPterm·dom signal from the 1000 interneurons was very weak for cylinder distances less than about 100 μm, corresponding to a 40% overlap between the two cylinders (see Fig 8B and 8C ). The amplitude of the LFPterm·dom signal increased with the cylinder distance together with the current dipole moment ( Fig 8C and 8D ; see Methods ). The “transition distance” of about 100 μm is seen to be associated with the appearance of an inversion point in the LFPterm·dom ( Fig 8C ) and with the establishment of a sizable dipole moment ( Fig 8D ). Above this transition distance the LFPterm·dom became larger with larger cylinder separations, yet saturating somewhat for distances above about 250 μm, corresponding to our reference model of layer 2/3 pyramidalterm·key cell. This demonstrates that the lack of a sizable contribution to the overall LFPterm·dom from our interneurons in the reference model was due to their stellate morphologies.

Below the inter-cylinder transition distance all proxies performed poorly with average fraction of variance explained across depths smaller than 70% (for 100 μm the range was <R 2 > between 0.37 and 0.64), but <R 2 > quickly saturated as soon as the dipole appeared ( Fig 8E ). <R 2 > was smaller for all proxies compared to the reference case (since the noise was larger due to the smaller number of neurons, i.e., 1000 neurons versus 4000 neurons for the reference case), but the ranking of performances for different proxies remained roughly the same: above the transition distance the fraction of variance explained by WS was 83%, RWSterm·mth and the sum of absolute values of currents both explained 80%, the membrane potentialterm·key and sum of synaptic currentsterm·dom 59%, while firing rateterm·key explained only 47% of the variance. Note that for inter-cylinder distances above the transition distance, the stable performance of the proxies were accompanied by stable values of the optimal coefficient α ( Fig 8F ). This result implies that the RWSterm·mth we have found for populations of layer-2/3-like pyramidalterm·key cells, likely also can be applied to pyramidalterm·key cell populations with different morphologies, as long as they produce a dipolar LFPterm·dom.

In order to verify that the assumptions we made to algorithmically construct the neuronal morphologies in the 3D network did not bias the results, we also did simulationsterm·mth using realistic morphologies obtained from anatomical reconstructions (see Methods ). The spatiotemporal dynamics of these LFPterm·dom signals was found to be qualitatively very similar to the one previously shown, and the agreement with proxies was even higher, with RWSterm·mth reaching R 2 = 0.95 ( S4 Fig ). This result indicates that our conclusions are not strongly dependent on the detailed branching patterns within the basal and apical dendritic bushes.

Dependency of the LFPterm·dom signal on the distribution of synapses In the reference case ( Fig 1B ) GABAterm·key synapses were distributed only in the lower cylinder while AMPAterm·key synapses were distributed homogeneously across all dendrites. In order to test how our results depended on this distribution we therefore evaluated all LFPterm·dom proxies for a variety of synaptic distribution patterns. Fig 9A illustrates the three main different synaptic distributions tested: (1) a case where all synapses were distributed homogeneously, (Hom.) (2) the reference case (Ref.), and (3) a case where AMPAterm·key synapses were located only in the upper cylinder (AM Up), leading to a complete separation between AMPAterm·key and GABAterm·key synapses. We further considered two conditions where (4) AMPAterm·key synapses were located only in the lower bush leaving the upper bush empty (AM down) and where (5) AMPAterm·key cortical synapses were located in the upper bush while thalamic AMPAterm·key inputs were distributed homogeneously (AM r Up).

LFPterm·dom signal and synaptic distribution. (A) Example cases for different synaptic distributions. Left: both AMPAterm·key and GABAterm·key synapse distributed over the entire surface of the cell. Center: GABAterm·key synapses distributed only in the lower cylinder, with AMPAterm·key synapses distributed over the entire cell. Right: GABAterm·key synapses distributed only in the lower cylinder and AMPAterm·key synapses only in the upper cylinder. (B) LFPterm·dom time course for different synaptic distributions: the three configurations presented in (A) correspond to black, red and green lines, respectively. Additional configurations were tested where GABAterm·key synapses were located in the lower cylinder, thalamic synapses were in both cylinders and the cortical AMPAterm·key synapses were only in the upper cylinder (blue line, AM r Up), and where all the AMPAterm·key synapses were located in the lower cylinder (cyan line, AM Down). (C) LFPterm·dom amplitude as a function of depth (similar to Fig 2A ) for different synaptic distributions. Blue and green markers were superimposed, illustrating that changing the position of the thalamic synapses does not alter the amplitude of the LFPterm·dom (only but its mean value, cf. panel B). (D) Average LFPterm·dom absolute amplitude over depths versus dipole moment (standard deviation over time) for the different synaptic distributions. (E) Contribution of AMPAterm·key and GABAterm·key currents to LFPterm·dom fluctuation amplitudes for different synaptic distributions. (F) Average fraction of LFPterm·dom variance explained by different proxies for different cylinder distances. Error bars are not displayed since they would not be visible in the figure. Same proxy arrangement as Fig 5C . (G) Mean and standard deviation across depths of optimal coefficients of α in the WS proxy as a function of synaptic distribution. Dashed line indicates the fixed coefficient of the RWSterm·mth proxy. Note that since for homogeneous synaptic distribution the performance of the fit was low (see panel (F)), the fitted value of the relative weight of AMPAterm·key and GABAterm·key currents in contributing to the LFPterm·dom signal has little significance.

Even though the parameters in the LIFterm·mth network and thus the output activity remained precisely the same as before in these different situations, the corresponding LFPterm·dom signal was dramatically altered by the choices of synaptic distributions ( Fig 9B ). The amplitude of the fluctuations was strongly affected, while the spatiotemporal features were only moderately altered. Note, however, that the position of the thalamic synapses only marginally affected the LFPterm·dom fluctuations, and only the mean value of the LFPterm·dom was affected. As a rule of thumb, we found that the more spatially segregated AMPAterm·key and GABAterm·key synapses are, the larger are the LFPterm·dom fluctuations ( Fig 9C ). We further observed that the variation of the LFPterm·dom amplitude on the synaptic distribution directly reflected changes in the magnitude of the current dipole moment ( Fig 9D ).

The individual contributions to the LFPterm·dom from AMPAterm·key and GABAterm·key synapses were strongly dependent on the spatial distributions ( Fig 9E ): when synapses were distributed homogeneously, the contribution of their currents to the LFPterm·dom signal was small as compared to when the synapses were segregated. Moreover, the AMPAterm·key contribution was larger when synapses were confined to the upper than to the lower cylinder. When the synapses were distributed homogeneously, the LFPterm·dom signal was very weak resulting in poor performances for all LFPterm·dom proxies ( Fig 9F ). When the cortical AMPAterm·key synapses were confined to the upper bush, the performance of the WS proxy was not affected, but a small decrease of 0.07 in the <R 2 > value was observed for both RWSterm·mth and the sum of the absolute values of synaptic currentsterm·dom. For the same situation there was a larger decrease of 0.17 in the <R 2 > value to a global value of only 0.51 for both the membrane potentialterm·key and the sum of synaptic currentsterm·dom. However, in the configuration in which AMPAterm·key synapses were confined to the lower bush and the LFPterm·dom amplitude was small, the <R 2 > for the membrane potentialterm·key and the sum of synaptic currentsterm·dom rose to 0.81 and 0.79 respectively, a value comparable to results for the WS and RWSterm·mth proxies (0.80 and 0.78). This suggests that the advantage of using the WS and RWSterm·mth proxies over, e.g., a membrane-potential proxy is particularly large when the AMPAterm·key and GABAterm·key synapses are spatially separated so that a large current dipole moment and a large amplitude LFPterm·dom is generated (Figs 8C and 9D ).

The best coefficients for WS strongly depended on the synaptic distribution ( Fig 9G ): When AMPAterm·key synapses were confined to the upper cylinder forming a strong current dipole moment, the optimal AMPAterm·key coefficients became larger than the GABAterm·key ones. Therefore, although the R 2 value of RWSterm·mth was still 0.82 under these conditions, a better result could be achieved with a proper tuning of the coefficients.

Difference between current-based and conductance-based synapses for the LFPterm·dom signal To keep the consistency with the LIFterm·mth network in which the synapses were current-based (see Methods ), all LFPterm·dom simulationsterm·mth considered until now were done using current-based synapses in the 3D network. However, some neuronal features may be better approximated by conductance-based synaptic models in which the postsynaptic currentsterm·dom (PSCs) depend on the local membrane potentialterm·key and do not have a fixed shape as in the case of current-based synapses. To test this situation, we repeated our simulationsterm·mth by introducing conductance-based synapses in the 3D network. Synaptic time constants were left unchanged, while the peak conductance values were scaled to obtain PSC amplitudes equivalent to current-based synapses for the reference stimulus intensity 1.5 spikesterm·dom/ms [ 47 ].

While the simulated LFPterm·dom amplitude was smaller when using conductance-based instead of current-based synapses (compare the three panels in Fig 10A with the three panels in Fig 6B and note the different y-axis scales), the time course was similar. We found that the explanatory power of the proxies was similar or better in all cases compared to the situation with LFPsterm·dom computed with current-based synapses ( Fig 10A ): the R 2 values for the RWSterm·mth were in the range 0.91–0.93 for inputs between 1 and 3 spikesterm·dom/ms, and 0.88 for 6 spikesterm·dom/ms. We hypothesize that the main reason for the increase in performance was that the LFPterm·dom contributions from different neurons were more correlated when synapses were conductance-based [ 47 ]. Note that in the case with conductance-based synapses, the performance of the membrane potentialterm·key proxy is in the very low 0.5–0.6 range for R 2 for all stimuli above 1 spikesterm·dom/ms. This can be understood given that the membrane potentialterm·key no longer depends linearly on synaptic input currents as in the case with current-based synapses.

Effects of modulation of inputs with conductance-based synaptic model. (A) LFPterm·dom (black line) for 150 μm depth position when the stimulation frequency was 0.5 spikesterm·dom/ms (left), 1.5 spikesterm·dom/ms (middle), 6 spikesterm·dom/ms (right) and corresponding best fit with WS proxy (red). The depth was the one for which WS proxy performance was highest. The titles indicate the fraction of variance explained. Note the different vertical scales. (B) Average fraction of LFPterm·dom variance over all depths explained by different proxies for different thalamic inputs. Error bars are not displayed since they would not be visible in the figure. Same proxy arrangement as in Fig 5C . (C) Mean and standard deviation across depths of optimal coefficients of α in the WS proxy as a function of thalamic input. Dashed line indicates the fixed coefficient of the RWSterm·mth proxy.

The WS proxy coefficients for 1 spikesterm·dom/ms inputs were rather similar to the current-based case, but when the input frequency was increased, the optimized value of the coefficient α, describing the ration of GABAterm·key to AMPAterm·key currents in the WS proxy, increased ( Fig 10C ). This likely reflects that for stronger stimuli the neurons were more depolarized, so that the average membrane potentialterm·key was closer to the AMPAterm·key reversal potential and further away from GABAterm·key reversal potential. Consequently, the GABAterm·key versus AMPAterm·key PSC-amplitude ratio increased. Nevertheless, the RWSterm·mth still performed well for all inputs ( Fig 10B ).

Discussion The main aim of this work was to develop an accurate, robust and an easy-to-use method to synthesize the LFPterm·dom signal from output from a model network of LIFterm·mth neurons. We simulated a biophysically realistic LFPterm·dom signal using a population of morphologically detailed multi-compartmental neuron models and compared this LFPterm·dom signal with several LFP proxyterm·key candidates extracted from the LIFterm·mth network simulationsterm·mth. We found that a linear combination of summed and time-shifted AMPAterm·key and GABAterm·key currents in the LIFterm·mth network explained a large fraction of the variance of the LFPterm·dom of the 3D morphologically accurate network in nearly all conditions considered. Specifically, we identified a specific set of parameters ( Eq 4 ), the so-called reference weighted sum LFP proxyterm·mth (RWSterm·mth), which accurately predicted the LFPterm·dom time course for all considered electrode positions ( Fig 5 ), and across different dynamic network states ( Fig 6 ).

The fraction of LFPterm·dom variance explained by the RWSterm·mth proxy was only moderately affected by changes in neural morphology ( Fig 8 , S4 Fig ), in synaptic distribution ( Fig 9 ), or in synaptic dynamics ( Fig 10 ). This LFP proxyterm·key was found to be very accurate for every condition considered in which the dipole generated by the synaptic currentsterm·dom was sizable and hence the amplitude of the LFPterm·dom substantial. This LFP proxyterm·key only performed poorly in situations where the amplitude of the LFPterm·dom signal itself was very small, i.e., at the inversion point or when the resulting current dipole moments from synaptic activation are small (homogeneous synaptic distributions ( Fig 9 ), low activity ( Fig 6 )). Furthermore, we showed that despite the complexity of our LFPterm·dom simulationterm·mth setup (with 5000 different morphologies with realistic dendritic structures) the temporal evolution of the LFPterm·dom was well captured by the RWSterm·mth proxy based on synaptic currentsterm·dom with adjustment of only three parameters, the relative weight of the contributions from the two synaptic currentsterm·dom and the two synaptic delays. Our results further suggested that the RWSterm·mth proxy can be used for a wide range of LIFterm·mth network models and pyramidal-neuron morphologies to synthesize biophysically plausible LFPterm·dom signals that can be compared with experimental LFPterm·dom recordings. Table 5 describes how to properly use the proxy in a variety of conditions and modeling approaches.

Summary of results for proxies and suggested adaption of results to other situations.

Thanks to its robustness, our proxy can expectedly be applied to models of any brain area in which the LFPterm·dom is likely to be generated by one dominant population (as in the hippocampus with a single population of pyramidalterm·key cells). When there are two or more populations giving a significant contribution to the LFPterm·dom (as is likely in whole cortical column model taking into account several pyramidalterm·key neuron populations, e.g., layer 2/3, layer 5, layer 6), the total LFPterm·dom can be approximated as a suitable linear combination of individual contributions if information on the depth positions of the populations relative to the recording electrode is available. Comparison of the model LFPterm·dom with experimental results might then be used to estimate the relative weights of the LFPterm·dom contributions from the different populations.

Comparison of candidate proxies A major difference between the accurate LFPterm·dom proxies using synaptic currentsterm·dom (sum of currents, WS, RWSterm·mth) compared to the less accurate proxy based on firing ratesterm·key is that a spiketerm·dom is a very local event in time, while the postsynaptic currentterm·dom following after a spiketerm·dom (as well as the contribution to the LFPterm·dom) lasts for many milliseconds. So an instantaneous firing rateterm·key proxy like the ones we are considering based on firing ratesterm·key cannot be expected to perform well (even with a fixed delay). In laminar population analysisterm·mth (LPAterm·mth, [ 30 ]) the LFPterm·dom time course was rather assumed to be given by the measured firing ratesterm·key convolved with a suitable (i.e., delayed exponential) kernel, the rationale being that spikesterm·dom causally drive synaptic currentsterm·dom which in turn set up the LFPterm·dom. The present RWSterm·mth proxy is similarly constructed, effectively corresponding to a suitable weighted sum of exponentially convolved presynaptic spiketerm·dom rates corresponding to excitatory and inhibitory synaptic currentsterm·dom. The postsynaptic soma membrane potentialsterm·key following presynaptic spiking is more low-pass filtered than the synaptic currentsterm·dom (and also the transmembrane return currents in the case of multicompartmental models) [ 29 ], and LFPterm·dom proxies based on this dynamical variable will generally fail to predict the most rapid temporal changes in the LFPterm·dom.

Spectral properties An interesting result is that all the proxies tested here displayed largely the same modulation of the LFPterm·dom gamma power as a function of input intensity, both in terms of relative power modulation and peak frequency ( Fig 7D and 7E ). This is encouraging since we did not specifically aim to find a good prediction of the power spectrum when constructing the LFPterm·dom proxies and estimating their parameters. We note however that no proxy is fully able to account for the low-frequency end of the spectrum ( Fig 7A and 7B ), which is overestimated by the membrane potentialterm·key proxy and underestimated by the other proxies. If one is interested in a highly detailed reproduction of the whole LFPterm·dom spectrum, preliminary results hint to the possibility of designing a WS fit optimized to match the spectrum instead than the spatiotemporal features and to define an LFP proxyterm·key that slightly differs from the RWSterm·mth discussed above. However, the fraction of spectral variance explained by the RWSterm·mth is already 0.91 (average over all stimulus intensities above 0.5 spikesterm·dom/ms, standard morphology and synaptic condition) which likely is sufficient for most purposes.

Model assumptions and limitations In the present work we have focused on how the relationship between LIFterm·mth variables and ground-truth LFPterm·dom change when the 3D model features change, keeping the LIFterm·mth model fixed. While different LIFterm·mth networks would generate different activity and hence different synaptic currentsterm·dom, we expect roughly the same relationship between these synaptic currentsterm·dom and the generated LFPterm·dom. Therefore, for any LIFterm·mth network generating enough correlated activity to result in a sizeable LFPterm·dom, we expect RWSterm·mth to be a good proxy.

Our strategy had the advantage that we could vary the assumptions in the LFP-generating model, e.g., the distribution of synapses or neuronal morphologies, without affecting the spiking dynamics. The disadvantage of this approach is, however, that the 3D network does not match the LIFterm·mth network in every respect; for instance, even though the synaptic input currents were identical in the two models, the resulting soma potentials in the multi-compartmental neurons were not identical to those in the LIFterm·mth neurons (due to passive dendritic filtering). It is, however, unlikely that imposing identical somatic potentials, or identical currents entering the soma, in the two models would result in a more realistic LFPterm·dom since large synaptic currentsterm·dom would be needed to overcome the passive filtering for distant synapses. Instead one could consider changing the synaptic weight distribution in the LIFterm·mth network simulationterm·mth to make the two models match better. Our focus here was to use LIFterm·mth models as commonly used in the literature (typically using homogeneous weight distributions), but it would be an interesting topic for future studies to extract effective point-neuron synaptic weight distributions from the multi-compartmental population and use these in the LIFterm·mth network simulationsterm·mth in order to make the two simulationterm·mth environments even more similar.

We did not test different LIFterm·mth network architectures or sizes, but we expect the RWSterm·mth proxy to be applicable as long as the network displays a sufficient level of correlation. We have found in previous modeling studies [ 20 , 21 ] that correlated synaptic activity is necessary to create a sizable LFPterm·dom signal, and in this case all cells in the dominant LFP-generating population will contribute. Making the network size larger or altering its connectivity would therefore likely not qualitatively change the form of the best LFP proxyterm·key (as long as a sufficient level of spiking correlations is maintained in the network). The LFPterm·dom generated by larger populations, however, should be tested in further studies taking into account the summed effect of several cortical populations, across layers as well as heterogeneous spatial structure in the horizontal direction.

A limitation of the presented simulationterm·mth setup is that it models only AMPAterm·key and GABAterm·key synapse contributions. However, most of our results pertaining to the proxy do not depend on the particular feature of the synapses and are therefore likely to extend to different synapses as well. For instance, it should not matter for the quality of our suggested proxy whether or not the synaptic weights are changing due to plasticity since the weight changes will be reflected in the synaptic currentsterm·dom extracted from the LIFterm·mth network as well. Including slower synapses, such as NMDA synapses in the model setup, will on the other hand affect the LFPterm·dom frequency content, particularly at low frequencies. This effect could be captured by a proxy including NMDA in the sum of synaptic currentsterm·dom with a weight depending on the number and spatial distribution of NMDA synapses. As with the synaptic weight distributions discussed above, the inclusion of NMDA synapses when computing the LFP proxyterm·key presupposes that it is also included in the LIFterm·mth network model (which was beyond the scope of this study). Moreover, we did not model subthreshold active dendritic conductance [ 55 ], nor the active channels behind spiketerm·dom generation. The contributions from the latter is expectedly negligible for at least the low frequencies of the LFPterm·dom [ 56 ] (but see [ 57 – 59 ]), while the effect of the former should be explored in future projects.

The present suggested proxy assumes the LFPterm·dom contribution following spikesterm·dom to be spatiotemporally separable, i.e., factorizable into a product of a function of time with a function of space [ 30 ]. Due to, for example, the intrinsic filtering effect [ 29 , 36 , 60 ] this is not strictly true as the spatial distribution of the transmembrane currentsterm·dom setting up the LFPterm·dom depends to some extent on the frequency. However, if warranted the present proxy can be extended, for example by assuming a more detailed proxy consisting of a sum of such spatiotemporally separable kernels.

Importance of this work for analytical studies of LFPsterm·dom Recently, we presented an analytical method to estimate the LFPterm·dom spectrum from the dynamics of a LIFterm·mth network [ 61 ] using as LFP proxyterm·key the sum of the absolute values of synaptic currentsterm·dom. By fitting a recurrent excitatory-inhibitory LIFterm·mth network model to LFPterm·dom recordings from monkeys presented with visual stimuli, we were able to estimate the LIFterm·mth model that best fitted the observed LFPterm·dom, and to predict at least in part the observed firing rateterm·key and some of the visual features in the receptive field that elicited the observed neural activity. In this recent work [ 61 ], the time evolution of the LFPterm·dom was computed analytically from the LIFterm·mth network as a function of the external input by applying linear response theory to the mean-field approximations of each kind of synaptic currentterm·dom separately and then summing their absolute values over pyramidalterm·key neurons (as in [ 15 ] and in Eq 2 ). In principle, it is possible to extend this analytical calculation by using the more efficient proxy presented here by simply changing the coefficients in the final sum of the synaptic currentsterm·dom. This paves the work for obtaining realistic analytical estimations of LFPsterm·dom from recurrent LIFterm·mth networks. As discussed in [ 35 ], an efficient analytical approach could be at the heart of the development of model-based analysis methods for performing inferential statistics of network models on LFPsterm·dom, analogous to the role played by Dynamic Causal Modelling [ 62 , 63 ] in the analysis of EEG and fMRI recordings.

Outlook Here we studied proxies for the LFPterm·dom produced by a local 3D network, corresponding to a single cortical layer. Experimentally recorded LFPsterm·dom, however, are most likely containing contributions from several layers [ 20 ]. Therefore, a natural extension of this work would be to study the LFPterm·dom generated by several connected 3D networks forming a full cortical columns [ 64 , 65 ] and determine how LFPterm·dom proxies should be designed in this context. Since electrical potentials in the nervous tissue are assumed to add linearly, we expect LFPterm·dom proxies to be constructed in largely the same manner as presented here, by summing synaptic contributions from different cortical layers, possibly with a weighting depending on the recording depths. Constructing the LFPterm·dom signal from a full cortical column model [ 65 ] is the topic for a separate ongoing project [ 66 ].

We expect our proxy to also work well for other brain structures where pyramidalterm·key neurons are elongated and arranged in an almost parallel way, such as the CA1 and CA2 regions of the hippocampus. On the other hand, many subcortical structures have a neuronal architecture so different from the cortex that we that we cannot a priori expect the present rules of LFPterm·dom prediction to be applicable. A possible future line of research will be to apply the combination of LIFterm·mth dynamics and 3D morphology we used in this work to investigate such areas to find a compact way to study the mechanisms generating the LFPterm·dom observed there.

We focused in the present study on the LFPterm·dom signal, but finding good models for relating activity in spiking network models and experimentally measured signals is relevant also for other types of commonly recorded signals such as the EEG, MEG and VSD. Since the biophysical mechanisms generating these signals are in principle known, we believe our framework could be extended to also study other measurement modalities in the future.

Methods Leaky integrate-and-fireterm·mth (LIFterm·mth) recurrent network model We summarize here the structure of the LIFterm·mth network that generated the spiking dynamics. We refer to [ 15 , 46 ] for full details. The network was composed of LIFterm·mth neurons with current-based synapses whose time evolution was modeled as difference between exponentials (see below), fixed threshold, fixed refractory time [ 67 ], and fixed conduction delay of 1 ms. Subthreshold dynamics for each neuron were given by τ m V ˙ m ( t ) = − V m ( t ) + ∑ P S C s y n ( t ) (5) where τ m corresponded to the membrane time constant due to the leak, V m was the membrane potentialterm·key, and PSC were the occurring synaptic events as a function of time t. When the membrane potentialterm·key V m crossed a threshold value of 18 mV above resting potential, a spiketerm·dom occurred, the potential dropped to a reset value of 11 mV above the reset potential and no spiketerm·dom could be emitted for a refractory time of 2 ms.

Post-synaptic currents (PSCs) were determined by the spikesterm·dom emitted by the pre-synaptic neurons in the LIFterm·mth network as well as by the external inputs. The time course of PSCs was described by the difference of two exponentials simulating the opening and closing process of the synapse. The equation can be written with two first order differential equations introducing the auxiliary variable x : τ d s y n P S C ˙ ( t ) = − P S C ( t ) + x ( t ) (6) τ r s y n x ˙ ( t ) = − x ( t ) + τ m ( J s y n ∑ s y n δ ( t − t s y n − τ l ) ) (7) where τ r / dsyn indicate the rise and decay times of the synapses, and J syn indicates the synaptic strength. The latency time of the synapses τ l was set to 1 ms. Compound synaptic currentsterm·dom were the linear sum of contributions induced by single pre-synaptic spikesterm·dom occurring at time t syn . We included two types of synapses: AMPAterm·key and GABAterm·key. Pyramidalterm·key neurons had AMPA-like synapses, and interneurons had GABA-like synapses. Moreover, each neuron received excitatory external drive from (1) a long range cortico-cortical input activating AMPAterm·key synapses identical to those of the recurrent connections and (2) a thalamic input activating AMPAterm·key synapses with timescales and strengths resembling those of thalamocortical synapses. Synaptic parameters such as rise time, decay time, and amplitude depended on the type of synapse and the category of the post-synaptic neuron. All simulationterm·mth parameters were in the range of the values reported in the literature [ 68 – 70 ] and are listed in Table 2 . We verified that modifying these values did not affect the results qualitatively [ 15 , 46 ].

The default network was composed of 4000 pyramidalterm·key neurons and 1000 interneurons ( Fig 1A ). The LIFterm·mth network connectivity was random and sparse, with a directed connection probability of 0.2 between any pair of cells. This resulted in an inhomogeneous connectivity with an average of 1000 pre- and post- synaptic connections for each cell. Each neuron received inhibitory and excitatory inputs from the neurons in the network, and also cortico-cortical and thalamic excitatory drives as described above. The long-range cortico-cortical drive represented the ongoing activity and the global contributions from other areas of cortex. Since ongoing cortical activity has most power for slow frequencies, this external drive was generated by an Ornstein-Uhlenbeck process with a low pass cut-off frequency of 10 Hz and a 0.25 mV standard deviation. Thalamic inputs were time-invariant in this set of simulationsterm·mth. Synapses carrying both types of external inputs were activated by random Poisson spiketerm·dom trains, with time-varying rates identical for all neurons. Details can be found in Table 1 and 2 .

Simulationsterm·mth were computed with time steps of 0.05 ms and lasted 10.1 seconds, with the first 100 ms removed to limit the analysis to the network steady state.

Current based and conductance based LIFterm·mth model source codes are identical to those used in [ 47 ] and are already available on the ModelDB sharing repository ( http://senselab.med.yale.edu/ModelDB/ShowModel.asp?model=152539 ) with accession number 152539.

Morphological model of a simplified layered cortical column In order to compute the transmembrane currentsterm·dom that lead to an LFPterm·dom signal, we constructed 3D morphological neuron models that captured the main morphological features of the cortical network described by point neurons in the LIFterm·mth model. The algorithm used to construct the model morphologies was based on the fact that dendrites connect to their presynaptic partners in a manner minimizing their total length and conduction times from all synapses to the soma [ 71 ]. In such a framework, pyramidalterm·key cell dendrites can be seen as tree structures connecting as directly as possible to axons that are distributed in two distinct layers [ 51 ]. The generation of synthetic trees and subsequent analysis were performed using the TREES toolbox [ 71 , 72 ] [ http://www.treestoolbox.org ]. Two cylinders (250 μm radius and 250 μm height each) were therefore stacked to form a cylindrical column ( Fig 1B ). Somata of all cells were homogeneously distributed in the lower cylinder for both cell types. Axons were distributed isotropically within planes perpendicular to the cortical depth at random depth values. Pyramidalterm·key cells were connected first to the axons in the upper cylinder and then to the axons of the lower cylinder, this resulted in characteristic apical and basal dendritic trees. Stellate cells were only connected to the axons in the lower cylinder. Using 160 axons in each layer and a maximal reach distance of 150 μm for any dendrite to an input axon, resulted in realistic membrane surfaces, cable lengths and branch point number distributions (see S1 Fig ). Diameter taper was selected to equalize synaptic democracy [ 73 ] and yielded good fits to the real counterparts with similar parameters for interneurons and pyramidalterm·key cells. The resulting pyramidalterm·key cell somatic input resistance was about 200 MΩ with specific membrane resistances R m = 20000 Ωcm 2 and axial resistances R a = 150 Ωcm. The stellate cell input resistance was 175 MΩ with specific membrane resistances R m = 10000 Ωcm 2 and axial resistances R a = 150 Ωcm. In order to spatially embed the simplified LIFterm·mth network model, 4000 such pyramidalterm·key cells and 1000 interneurons were generated to populate the simplified columnar architecture. The resulting morphologies were then exported to NEURON [ 74 , 75 ] using the TREES toolbox functions.

In order to continuously alter the “pyramidalness” of cortical neurons as in Fig 8 we simply modulated the distance between the two cylinders corresponding to the two layers. With a distance of 0 μm, a perfect overlap of both cylinders, the resulting shape was symmetric as for the stellate cell. As the distance was increased between the two cylinders, the shape of the cortical cell traversed the shape of layer 2/3 pyramidalterm·key cells (distance of 250 μm), layer 4 pyramidalterm·key cells (distance of 350 μm) and layer 5 pyramidalterm·key cells (distance larger than 500 μm). The corresponding validation of morphological features compared with real dendrite reconstructions as can be observed in S1 Fig .

As a control for use of algorithmically constructed morphologies we derived an alternative model using multiple copies of real reconstructions distributed within the columnar arrangement. We used reconstructions from NeuroMorpho.org [ 76 ], made available by the group of Markram [ 76 ], of both stellate cells and layer 2/3 pyramidalterm·key cells in young rat somatosensory cortex. Since only 4 stellate cells and 36 layer 2/3 pyramidalterm·key cell morphologies were available, we reached the number of 1000 interneurons and 4000 pyramidalterm·key neurons by randomly selecting copies of the smaller sample and distributing them within the simplified columnar geometry. Cell body locations were chosen to preserve a fairly homogeneous distribution of membrane throughout the cylinders. This alternative model was then injected with the same synaptic currentterm·dom stimuli as the original model based on algorithmically developed morphologies, and yielded similar results (compare Fig 6 and S4 Fig ).

Connecting two simulationterm·mth frameworks Spiketerm·dom trains generated by the LIFterm·mth network were used as input in the 3D network model used for LFPterm·dom generation. Each multi-compartmental neuron model in the 3D network was associated with a given point neuron in the LIFterm·mth network. To make sure the total synaptic currentsterm·dom in each cell were identical in the two simulationterm·mth environments, we used the connectivity structure of the LIFterm·mth network to determine the presynaptic LIFterm·mth neurons for each postsynaptic multi-compartmental neuron in the 3D network. We triggered the synaptic currentsterm·dom in the multi-compartmental neurons of the 3D network at the precise times given by the spiketerm·dom trains generated by the presynaptic cells during LIFterm·mth network simulationsterm·mth. Note that we did not take into account synaptic latency time. In the 3D network we associated with each presynaptic cell a single specific synapse in the postsynaptic cell. Synaptic dynamics in the 3D network was identical to the one in the LIFterm·mth network (Eqs 4 and 5 ). In addition we recreated the external inputs (“Thalamic” and “Cortical”, see Fig 1A ) used in the LIFterm·mth network simulationsterm·mth and injected the same patterns of external spiketerm·dom trains in specific AMPAterm·key synapses in the 3D network neurons.

Since the LIFterm·mth neuron model used in the LIFterm·mth network simulationsterm·mth lacked spatial structure, we needed to make additional assumptions regarding the synapse placement when simulating the multi-compartmental neurons in the 3D network. The cylinders that were used to create the morphologies of the multi-compartmental models (see above) were also used to broadly define the synaptic regions. Our default setting was to place GABAterm·key synapses only in the lower cylinder, while AMPAterm·key synapses were placed in both cylinders. We tested also other scenarios in the “Dependency of the LFPterm·dom signal on the distribution of synapses” subsection of Results ( Fig 9 ). We randomly chose the detailed spatial position on the dendritic structure for each synapse, with the probability for a section to be selected being proportional to its membrane area, such that the resulting synaptic density was homogeneous within the selected cylinder Calculation of local field potentialterm·dom (LFPterm·dom) from the morphological 3D network We calculated the model LFPterm·dom signal from the transmembrane currentsterm·dom in the multi-compartmental neuron populations based on volume conduction theory and the line-source approximation implemented in the Python package LFPy ( http://lfpy.github.io/ ) [ 34 ]. We first simulated transmembrane currentsterm·dom resulting from synaptic activity using the NEURON simulationterm·mth environment [ 74 , 75 ] after which extracellular potentials were calculated as a weighted sum of those transmembrane currentsterm·dom [ 31 , 32 , 34 ]. The extracellular potentials were computed for 32 equispaced vertically aligned points in space (simulating a laminar multielectrode), set at 25 μm intervals along the central vertical axis of the 3D network cylinder ( Fig 1B ). For the analysis illustrated in the subsection “Spatial distribution of simulated LFPterm·dom signal” the recording locations were set at different distances from the vertical axis of the 3D network cylinder. To directly match the LIFterm·mth network simulationsterm·mth, morphological neurons used current synapses in the reference case, except in the simulationterm·mth discussed in the subsection “Difference between current-based and conductance-based synapses for the LFPterm·dom signal” were conductance synapses were adopted ( Table 3E ). The calculations of transmembrane currentsterm·dom in the morphological model were performed using passive neuron models with the parameters listed above ( Table 4 ). Following volume conductor theory, the model neurons were assumed to be surrounded by an infinitely sized extracellular medium with conductivity assumed to be real, scalar (the same in all directions) and homogeneous (the same everywhere) with σ = 0.3 S/m [ 77 ]. For further discussion on these assumptions see [ 32 ].

The Python codes we used to generate LFPterm·dom from artificial morphologies injected with LIFterm·mth spiketerm·dom dynamics are available on the LFPy official site ( http://lfpy.github.io/ ).

LFPterm·dom proxies for LIFterm·mth networks We tested several simple models to match the LFPterm·dom simulationterm·mth based on the different variables describing the activity in the LIFterm·mth network: firing rateterm·key, membrane potentialterm·key, AMPAterm·key and GABAterm·key synaptic currentsterm·dom. We considered variables computed over the set of all pyramidalterm·key neurons, of all interneurons or both populations. We considered proxies based on these variables and on the simple sum or the sum of absolute values of synaptic currentsterm·dom as in [ 15 , 46 ]. Then we considered linear combinations of synaptic currentsterm·dom with different time delays. We tested the accuracy of the proxy in describing the time evolution of the LFPterm·dom given by the morphological model by using the mean of squared values of the correlation coefficient R (which is equivalent to the fraction of variance explained). The quality of the proxy was tested separately for each depth. We computed the cross-correlation between the simulated LFPterm·dom signal and the corresponding proxy and we determined the delay as the lag of the cross-correlation peak (see Fig 4 ). For this delay we determined the best linear fit using the Matlab function polyfit for single regressors and the Matlab function regress for regressor combinations. We estimated the quality of the proxy as the squared correlation coefficient between the best fit and the LFPterm·dom. The proxy for each depth is defined by the optimal delays and the coefficients of the different components for regressor combinations. To compare the performance of the different proxies taking into account the different number of free parameters between WS, RWSterm·mth and all the other proxies, we used the Bayesian Information Criterion (BIC, [ 54 , 78 ]) B I C = − 2 l + K log n (8) where l is the optimized loglikelihood function, K the number of estimable parameters and n the sample size. Under the assumption of Gaussian noise, −2 l can be approximated as c o n s t a n t + n log R S S n [ 79 ] where RSS is the sum of the residual squares, so the BIC criterion becomes B I C = n log R S S n + K log n (9) which is the criterion we adopted in the manuscript.

Supporting Information (A) Layer 2/3 pyramidalterm·key cells: Comparison of amount of cellular membrane surface between anatomically reconstructed cells (left) and synthetic morphologies (right). Results shown in units of μm 2 per cell for 25 μm bins of cortical depth. Green lines indicate profiles of individual cells and black lines are average traces. (B) As in panel A, but for stellate cells.

(TIF) Click here for additional data file.

(A) Fraction of LFPterm·dom variance explained by the different proxies at different depths. FR: Firing rateterm·key; Vm: membrane potentialterm·key; AMPAterm·key/GABAterm·key: input current values; ∑I: sum of the two input currents; ∑|I|: sum of the absolute values of the two input currents. Quantities with _exc/_inh subscript are computed only over the set of excitatory/inhibitory neurons, e.g., AMPAterm·key_exc is the average AMPAterm·key current input into excitatory neurons. Quantities with no subscript are computed over all neurons. (B) Optimal time lag for different depths for same proxies as previous panel. (C) Fraction of variance explained by combination of AMPAterm·key and GABAterm·key currents with same coefficients as Eq ( 4 ) and different delays. Optimal (and reference) combination is indicated by an X.

(EPS) Click here for additional data file.

Same as panel 5C but showing for each proxy the BIC value (see Methods ) instead of the fraction of variance explained.

(EPS) Click here for additional data file.

Same as Fig 6B and 6D when using real reconstructed morphologies of cortical stellate and layer 2/3 pyramidalterm·key cells (see Methods ).

(EPS) Click here for additional data file.