Abstract

With rapidly advancing multi-electrode recording technology, the local field potential (LFP) has again become a popular measure of neuronal activity in both research and clinical applications. Proper understanding of the LFP requires detailed mathematical modeling incorporating the anatomical and electrophysiological features of neurons near the recording electrode, as well as synaptic inputs from the entire network. Here we propose a hybrid modeling scheme combining efficient point-neuron network models with biophysical principles underlying LFP generation by real neurons. The LFP predictions rely on populations of network-equivalent multicompartment neuron models with layer-specific synaptic connectivity, can be used with an arbitrary number of point-neuron network populations, and allows for a full separation of simulated network dynamics and LFPs. We apply the scheme to a full-scale cortical network model for a ∼1 mm2 patch of primary visual cortex, predict laminar LFPs for different network states, assess the relative LFP contribution from different laminar populations, and investigate effects of input correlations and neuron density on the LFP. The generic nature of the hybrid scheme and its public implementation in

hybridLFPy
form the basis for LFP predictions from other and larger point-neuron network models, as well as extensions of the current application with additional biological detail.

Introduction

The local field potential (LFP), the low-frequency component (500Hz) of the extracellular potential recorded in the brain, is commonly used as a measure of neuronal activity (Buzsáki et al. 2012; Einevoll et al. 2013). The LFP originates from transmembrane currents (Nicholson and Freeman 1975), and at the single-cell level the biophysical origin of such extracellular potentials is well understood (see, e.g., Rall and Shepherd 1968; Holt and Koch 1999; Buzsáki et al. 2012; Einevoll et al. 2013). However, the interpretation of the LFP remains difficult due to the large number of neurons contributing to the recorded signal. In neocortex, for example, the measured LFP is typically generated by thousands or even millions of neurons near the recording electrode (Kajikawa and Schroeder 2011; Lindén et al. 2011; Łe¸ski et al. 2013). Moreover, the LFP reflects also synaptic input generated by remote populations, for example, inputs from other cortical or subcortical areas in addition to local network interactions (Herreras et al. 2015). A thorough theoretical description of the LFP, therefore, needs to account not only for the anatomical and electrophysiological features of neurons in the vicinity of the recording electrode, but also for the entire large-scale neuronal circuitry generating synaptic input to these cells.

Modeling large-scale neural-network dynamics with individual spiking neurons is challenging due to the memory required to represent the large number of synapses. With current technology and using the largest supercomputers available today, simulations of neural networks comprising up to 109 neurons and 1013 synapses (roughly corresponding to the size of a cat brain) are feasible for simplified model neurons (Diesmann 2013; Kunkel et al. 2014). Typically, these simplified models neglect the spatial aspects of neuronal morphologies and describe neurons as points in space (point-neuron models). Despite their simplicity, point-neuron-network models explain a variety of salient features of neural activity observed in vivo, such as spike-train irregularity (Softky and Koch 1993; van Vreeswijk and Sompolinsky 1996; Amit and Brunel 1997; Shadlen and Newsome 1998), membrane-potential fluctuations (Destexhe and Paré 1999), asynchronous firing (Ecker et al. 2010; Renart et al. 2010; Ostojic 2014), correlations in neural activity (Gentet et al. 2010; Okun and Lampl 2008; Helias et al. 2013), self-sustained activity (Ohbayashi et al. 2003; Kriener et al. 2014), and realistic firing rates across laminar cortical populations (Potjans and Diesmann 2014). Point-neuron networks are amenable to mathematical analysis (see, e.g., Brunel 2000; Deco et al. 2008; Tetzlaff et al. 2012; Helias et al. 2013; de Kamps 2013; Schuecker et al. 2015; Bos et al. 2016) and can be efficiently evaluated numerically (Brette et al. 2007; Plesser et al. 2007; Helias et al. 2012; Kunkel et al. 2014). The mechanisms governing networks of biophysically detailed multicompartment model neurons, in contrast, are less accessible to analysis and these models are more prone to overfitting. Existing multicompartment neuron network models accounting for realistic cell morphologies are restricted to sizes of ∼104–105 neurons (Hines et al. 2008; Reimann et al. 2013; Migliore et al. 2014; Markram et al. 2015). Large-scale models are, however, necessary to include contributions to the LFP from distant populations in situations where the spatial reach of the LFP is known to be large (Lindén et al. 2011; Łe¸ski et al. 2013).

Although point-neuron networks capture many features of in vivo spiking activity, they fail to predict extracellular potentials that result from transmembrane currents distributed across the cell surface. According to Kirchhoff's law of current conservation, the sum of all transmembrane currents, including all ionic and capacitive currents, must be zero for each neuron. In a point-neuron model, all transmembrane currents are collapsed in a single point in space. The net transmembrane current, and hence the extracellular potential, therefore vanishes. Only the spatial separation between current sinks and sources leads to a nonzero extracellular potential (Pettersen et al. 2012; Einevoll et al. 2013). A priori, the prediction of extracellular potentials, therefore, requires spatially extended neuron models accounting for the spatial distribution of transmembrane currents, commonly handled using multicompartment neuron models (De Schutter and Van Geit 2009). Note that the principle of current conservation implies a current sum rule for multicompartment neuron models as well: the sum of all single-cell transmembrane currents remains zero, also across neuron populations and the whole column. In several previous studies (Bazhenov et al. 2001; Hill and Tononi 2005; Ursino and La Cara 2006; Mazzoni et al. 2008; 2010; 2011), the activity of point-neuron networks (e.g., population firing rates, synaptic currents and membrane potentials) has nevertheless been used as a proxy for the LFP when comparing with experiments. In a recent study comparing different candidate proxies, it was found that a suitably chosen sum of synaptic currents could provide a good LFP proxy, but only for the case when the LFP is generated from transmembrane currents of a single population of pyramidal neurons (Mazzoni et al. 2015). In cortex, however, several populations in general contribute to the LFP, and there are spatial cancellation effects when positive LFP contributions from one population overlap in space with negative LFP contributions from other populations. This effect cannot be accounted for by a simple LFP proxy.

In this article, we present a hybrid modeling scheme that combines the simplicity and efficiency of point-neuron network models and the biophysical principles underlying LFP generation captured by multicompartment neuron models with anatomically reconstructed morphologies. The scheme allows for arbitrary numbers of LFP-contributing populations, and directly incorporates spatial cancellation effects. Furthermore, the spatially extended LFP-generating neurons assure that effects from intrinsic dendritic filtering of synaptic inputs are included in the predicted LFP (Lindén et al. 2010). The scheme assumes that the spiking activity of the neural network (Fig. 1B) generating the synaptic input reflected in the LFP is well described by a point-neuron network model (Fig. 1A). The network spiking activity serves as synaptic input to a population of mutually unconnected multicompartment model neurons with realistic morphologies positioned in 3-dimensional (3D) space (Fig. 1C) and is thereby translated into a distribution of transmembrane currents and, hence, an LFP (Fig. 1D). Thus each multicompartment model neuron has its equivalent in the point-neuron network and receives input spikes from the same presynaptic neurons as this point-neuron equivalent.

Figure 1.

Overview of the hybrid LFP modeling scheme for a cortical microcircuit model. (A) Sketch of the point-neuron network representing a 1 mm2 patch of early sensory cortex (adapted from Potjans and Diesmann 2014). The network consists of 8 populations of LIF neurons, representing excitatory (E) and inhibitory neurons (I) in cortical layers 2/3, 4, 5, and 6. External input is provided by a population of TC neurons and cortico-cortical afferents. The color coding of neuron populations is used consistently throughout this paper. Red arrows: excitatory connections. Blue arrows: inhibitory connections. See Tables 12, 56 for details on the network model. (B) Spontaneous (t<900ms) and stimulus-evoked spiking activity (synchronous firing of TC neurons at time t = 900 ms, denoted by thin vertical line) generated by the point-neuron network model shown in panel A, sampled from all neurons in each population. Each dot represents the spike time of a particular neuron. (C) Populations of LFP-generating multicompartment model neurons with reconstructed, layer-, and cell-type specific morphologies. Cells are distributed within a cylinder spanning the cortex. Layer boundaries are marked by horizontal black lines (at depths z relative to cortex surface z = 0). Only one representative neuron for each population is shown (see Fig. 4 for a detailed overview of cell types and morphologies). Sketch of a laminar recording electrode (gray) with 16 contacts separated by 100μm (black dots). (D) Depth-resolved LFP traces predicted by the model (cf. Tables 3 and 4). Note that channel 1 is at the pial surface, so that channel 2 corresponds to a cortical depth of 100 μm and so forth.

Figure 1.

Overview of the hybrid LFP modeling scheme for a cortical microcircuit model. (A) Sketch of the point-neuron network representing a 1 mm2 patch of early sensory cortex (adapted from Potjans and Diesmann 2014). The network consists of 8 populations of LIF neurons, representing excitatory (E) and inhibitory neurons (I) in cortical layers 2/3, 4, 5, and 6. External input is provided by a population of TC neurons and cortico-cortical afferents. The color coding of neuron populations is used consistently throughout this paper. Red arrows: excitatory connections. Blue arrows: inhibitory connections. See Tables 12, 56 for details on the network model. (B) Spontaneous (t<900ms) and stimulus-evoked spiking activity (synchronous firing of TC neurons at time t = 900 ms, denoted by thin vertical line) generated by the point-neuron network model shown in panel A, sampled from all neurons in each population. Each dot represents the spike time of a particular neuron. (C) Populations of LFP-generating multicompartment model neurons with reconstructed, layer-, and cell-type specific morphologies. Cells are distributed within a cylinder spanning the cortex. Layer boundaries are marked by horizontal black lines (at depths z relative to cortex surface z = 0). Only one representative neuron for each population is shown (see Fig. 4 for a detailed overview of cell types and morphologies). Sketch of a laminar recording electrode (gray) with 16 contacts separated by 100μm (black dots). (D) Depth-resolved LFP traces predicted by the model (cf. Tables 3 and 4). Note that channel 1 is at the pial surface, so that channel 2 corresponds to a cortical depth of 100 μm and so forth.

In the proposed hybrid modeling scheme, the LFP stems from the presynaptic spiking activity, but does not affect the spike-generation dynamics. Thus, the modeling of the spike trains and the LFP generation are separated so that the effects of the spatial and electrophysiological properties of the postsynaptic (multicompartment) neurons on the LFP can be investigated independently of the spike-generation dynamics. Due to the linearity of Maxwell's equations and volume conduction theory linking transmembrane currents to extracellular potentials (Pettersen et al. 2012; Einevoll et al. 2013), the compound LFP results from the linear superposition of all single-cell LFPs generated by the collection of neurons in the multicompartment model neuron population (Einevoll et al. 2013). Note that this linear superposition principle applies even for nonlinear cell dynamics (e.g., nonlinear synaptic integration, action-potential generation and active conductances) as in Reimann et al. 2013. As ephaptic interactions (Anastassiou et al. 2011) are neglected, the LFP contribution from each multicompartment model neuron can be treated independently from the others. The computational hybrid LFP scheme proposed here exploits the methodological and conceptual advantages due to the independence of the contributions to the LFP from each multicompartment model neuron: the evaluation of the LFPs becomes “embarrassingly parallel” (see Foster 1995) and simulations of the multicompartment model neuron dynamics can be easily distributed in parallel across many compute units (i.e., CPUs). Although tailored toward use on high-performance computing facilities, the hybrid simulation can in principle be run on a single laptop.

The hybrid scheme predicts spatially and temporally resolved neural activity at various scales: spikes, synaptic currents, membrane potentials, current-source densities (CSD, see e.g., Nicholson and Freeman 1975; Pettersen et al. 2006; 2008), and LFPs. It therefore allows for investigation of relationships between different measures of neural activity. Thus, although point-neuron networks until now only have connected to in vivo experiments via measurement of spikes, single-neuron membrane potentials and currents, the present hybrid scheme allows for comparison of model predictions also with measured LFPs (and associated CSDs).

As an illustration, we apply the hybrid scheme to a multilayered point-neuron network model of an early sensory cortical microcircuit (Potjans and Diesmann 2014). We thereby demonstrate how to obtain LFP predictions from point-neuron network models using additional spatial connectivity information from anatomical data (Binzegger et al. 2004; Izhikevich and Edelman 2008). The example illustrates how the hybrid scheme can be used to examine the relation between single-neuron and population signals, that is, spikes and LFPs, the effect of network dynamics on the LFP, and the interpretation of the LFP in terms of underlying laminar neuron populations. We further use the example to demonstrate that synaptic-input correlations result in a nontrivial dependence of the LFP on the neuron density. Correct LFP predictions can therefore only be obtained by accounting for realistic neuron densities.

The network model of Potjans and Diesmann (2014) is chosen here since it has a minimum level of detail in the sense that individual neurons have simplified leaky integrate-and-fire (LIF) dynamics, but still represents a cortical column with full density of neurons and connections. The connectivity in such a full-scale circuit alone suffices to explain realistic firing rates across populations as well as propagation of activity through layers (Potjans and Diesmann 2014). Applicability of the scheme is, however, not restricted to this model as it in principle can be used for all network models generating spikes.

In Methods and Materials, we detail the components of the hybrid scheme and their application to the cortical microcircuit model: the point-neuron-network model, the populations of multicompartment neurons, the synaptic connectivity of the point-neuron network and the multicompartment model neuron populations, and the biophysical forward-modeling scheme of extracellular potentials. We further describe the analysis of the data generated by the simulations, as well as the

hybridLFPy
software implementation. In Results, we apply the hybrid scheme to the cortical microcircuit model of Potjans and Diesmann (2014) and study the effects of network dynamics on the LFP, the contributions of individual cortical subpopulations to the LFP, the role of correlations and neuron density, and how well the LFP can be predicted from population firing rates (rather than from individual spikes). In Discussion, we outline implications of our work, and in particular future applications and extensions of the hybrid LFP modeling scheme.

Methods and Materials: Hybrid LFP Modeling Scheme

Point-Neuron Network Model

The point-neuron network model is a key component of the hybrid scheme. The hybrid scheme enables LFP predictions from network models with an arbitrary number of populations and thus permits application to a large class of networks with arbitrarily complex single-neuron and synapse dynamics. The example network of “spike-generators” used here is, except for some minor adjustments (see below), the multilayered model of a cortical microcircuit published by Potjans and Diesmann (2014). The model is implemented and included in

NEST
(http://www.nest-simulator.org, Eppler et al. (2015)) and was recently made freely available (http://www.opensourcebrain.org/projects/potjansdiesmann2014).

The network model describes 1 mm2 of primary sensory cortex and consists of 4 layers with one excitatory (E) and one inhibitory (I) neuron population each, as illustrated in Fig. 1A. The network receives modulated thalamic input in addition to stationary external input. While the neuron (LIF) and synapse (static, exponential-current-based) model are intentionally left simple, the focus of this network implementation is on the complex connectivity which integrates multiple sources of anatomical and electrophysiological data (Potjans and Diesmann 2014). Apart from the layer identity, the model does not explicitly account for cell positions. For the full network description, see Tables 1, 2, and 5. The microcircuit model reproduced experimentally observed distributions of firing rates across populations and propagation of activity across layers (Potjans and Diesmann 2014). It thus forms a suitable starting point for LFP predictions in a cortical column.

Table 1

Description of point-neuron network for the cortical microcircuit model (continued in Table 2) following the guidelines of Nordlie et al. (2009)

Model summary 
Structure Multilayered excitatory-inhibitory (E–I) network 
Populations 8 cortical in 4 layers, 1 thalamic (TC) 
Connectivity Random, independent, population-specific, fixed number of connections 
External input Cortico-cortical: constant current with population-specific strength 
Neuron model Cortex: LIF; TC: point process 
Synapse model Exponential postsynaptic currents, static weights, population-specific weight distributions 
Measurements Spike activity, input currents, membrane potential of each neuron 
Network model 
Connectivity Connection probability CYX (X,Y{L2/3,L4,L5,L6}×{E,I}TC, CYX=0 for Y=TC
 •Fixed number of synapses KYX between populations X and Y 
 •Binomial in-/out-degrees 
Input Cortico-cortical direct current IYext 
Neuron model 
Cortex 
Type Leaky integrate-and-fire neuron (LIF) 
Description Dynamics of membrane potential Vi(t) (neuron i[1,N]): 
  •Spike emission at times til with Vi(tli)θ 
  •Subthreshold dynamics: τmVi˙=Vi+RmIi(t) if l:t(tli,tli+τref] 
  •Reset + refractoriness: Vi(t)=Vreset if l:t(tli,tli+τref] 
 Exact integration with temporal resolution dt (Rotter and Diesmann 1999
 Uniform distribution of membrane potentials at t = 0 
Thalamus 
Type •DC current for constant background input 
 •Nonstationary Poisson process for modulation 
Description DC current included in external DC input 
 
 Types of thalamic input modulation: 
 •Spontaneous activity: no modulation in activation of thalamic neurons 
 •Thalamic pulses: fixed-interval coherent activation of all NTC thalamic neurons 
 •AC modulation: Poisson spike trains with sinusoidally modulated rate profile (discretized with time resolution dt):  
(18)
νth(t)=ν¯TC+ΔνTCsin(2πtfTC)
 
Model summary 
Structure Multilayered excitatory-inhibitory (E–I) network 
Populations 8 cortical in 4 layers, 1 thalamic (TC) 
Connectivity Random, independent, population-specific, fixed number of connections 
External input Cortico-cortical: constant current with population-specific strength 
Neuron model Cortex: LIF; TC: point process 
Synapse model Exponential postsynaptic currents, static weights, population-specific weight distributions 
Measurements Spike activity, input currents, membrane potential of each neuron 
Network model 
Connectivity Connection probability CYX (X,Y{L2/3,L4,L5,L6}×{E,I}TC, CYX=0 for Y=TC
 •Fixed number of synapses KYX between populations X and Y 
 •Binomial in-/out-degrees 
Input Cortico-cortical direct current IYext 
Neuron model 
Cortex 
Type Leaky integrate-and-fire neuron (LIF) 
Description Dynamics of membrane potential Vi(t) (neuron i[1,N]): 
  •Spike emission at times til with Vi(tli)θ 
  •Subthreshold dynamics: τmVi˙=Vi+RmIi(t) if l:t(tli,tli+τref] 
  •Reset + refractoriness: Vi(t)=Vreset if l:t(tli,tli+τref] 
 Exact integration with temporal resolution dt (Rotter and Diesmann 1999
 Uniform distribution of membrane potentials at t = 0 
Thalamus 
Type •DC current for constant background input 
 •Nonstationary Poisson process for modulation 
Description DC current included in external DC input 
 
 Types of thalamic input modulation: 
 •Spontaneous activity: no modulation in activation of thalamic neurons 
 •Thalamic pulses: fixed-interval coherent activation of all NTC thalamic neurons 
 •AC modulation: Poisson spike trains with sinusoidally modulated rate profile (discretized with time resolution dt):  
(18)
νth(t)=ν¯TC+ΔνTCsin(2πtfTC)
 
Table 2

Description of point-neuron network for the cortical microcircuit model (continuation of Table 1)

Synapse model 
Type Exponential postsynaptic currents, static weights 
Description Input current of neuron j of synapses formed with presynaptic neurons i
  
Ij(t)=iJjilexp((ttlidi)/τs)H(ttlidi)+Ijext
 
  •Static synaptic weights Jji=sgn(X)JYX (iX, jY); sgn(X)=1 for X{L2/3E,L4E,L5E,L6E,TC},1 otherwise 
  •Absolute weights JYX drawn from lognormal distribution  
(19)
p(JYX)=12πσYXJYXexp((lnJYXμYX)22σYX2)
or normal distribution  
(20)
p(JYX)=12πσYXexp((JYXμYX)22σYX2)
with μYX=gYXJ and σYX=σJ,relμYX 
  •Delays di=dX (iX) drawn from (left-clipped) Gaussian distribution  
(21)
p(dX)=12πσXexp((dXμX)22σX2)
with mean μX=dE,dI for X exc., inh., standard deviation σX=σd,relμX and dX[dt,) 
 H(t)=1fort0, and 0 elsewhere. 
  •External DC input Ijext=IYext=IextkYext (jY
Synapse model 
Type Exponential postsynaptic currents, static weights 
Description Input current of neuron j of synapses formed with presynaptic neurons i
  
Ij(t)=iJjilexp((ttlidi)/τs)H(ttlidi)+Ijext
 
  •Static synaptic weights Jji=sgn(X)JYX (iX, jY); sgn(X)=1 for X{L2/3E,L4E,L5E,L6E,TC},1 otherwise 
  •Absolute weights JYX drawn from lognormal distribution  
(19)
p(JYX)=12πσYXJYXexp((lnJYXμYX)22σYX2)
or normal distribution  
(20)
p(JYX)=12πσYXexp((JYXμYX)22σYX2)
with μYX=gYXJ and σYX=σJ,relμYX 
  •Delays di=dX (iX) drawn from (left-clipped) Gaussian distribution  
(21)
p(dX)=12πσXexp((dXμX)22σX2)
with mean μX=dE,dI for X exc., inh., standard deviation σX=σd,relμX and dX[dt,) 
 H(t)=1fort0, and 0 elsewhere. 
  •External DC input Ijext=IYext=IextkYext (jY
Table 3

Description of multicompartment-neuron populations for the cortical microcircuit model (continued in Table 4)

Model summary 
Topology Cortical column under 1mm2 of cortical surface 
Populations 8 excitatory and 8 inhibitory cell types 
Input Spiking activity of thalamic and cortical populations as modeled by point-neuron network 
Neuron model Multicompartment, passive cable formalism 
Synapse model Exponential postsynaptic current, static weights 
Measurements Current source density (CSD), local field potential (LFP) 
Topology 
Type Cylindrical volume with layer-specific distribution of cell types and synapses 
Description Cylinder radius r 
 Laminar, defining upper/lower boundaries of layers 1, 2/3, 4, 5, 6 
Populations 
Type Each cell type y assigned to population Y, yY 
Description Populations Y{L2/3,L4,L5,L6}×{E,I} (population size NY, cell types yY
 (e.g., L4E={p4,ss4(L23),ss4(L4)}, cf., Fig. 2). 
 Cell types y
 •Size Ny=FyYNY, Fy is the occurrence of cell type y in the full model 
 •Morphology My 
 •Extrapolated according to spatial connectivity data (Table 7
 Somatic placement, population Y
 •Random soma placement in cylindrical volumes with radius r, thickness h 
 •Volumes centered between boundaries of layers 2/3–6 
Morphologies 
Type 3D histological reconstructions from slice preparations (see Jacobs et al. 2009; De Schutter and Van Geit 2009) of cat visual and somatosensory cortices 
Description One morphology My per cell type: 
 •Excitatory and inhibitory cells in layers 2/3–6 
 •For all cells jy:Mj = My 
 •For some cell types y,y:My=My (limited availability) 
 Orientations: 
 •Pyramidal cells: apical dendrites oriented along depth axis with random depth-axis rotation 
 •Interneurons, stellate cells: random rotation around all axes 
 Corrections: 
 •Apical dendrites of pyramidal cells elongated to accommodate spatial connectivity 
 •Axons removed if present 
 Reconstructed morphologies (cf., Fig. 2): 
 •Cat visual cortex (Kisvárday and Eysel 1992; Mainen and Sejnowski 1996; Contreras et al. 1997; Stepanyants et al. 2008
 •Cat somatosensory cortex from NeuroMorpho.org (Contreras et al. 1997; Ascoli et al. 2007). 
Model summary 
Topology Cortical column under 1mm2 of cortical surface 
Populations 8 excitatory and 8 inhibitory cell types 
Input Spiking activity of thalamic and cortical populations as modeled by point-neuron network 
Neuron model Multicompartment, passive cable formalism 
Synapse model Exponential postsynaptic current, static weights 
Measurements Current source density (CSD), local field potential (LFP) 
Topology 
Type Cylindrical volume with layer-specific distribution of cell types and synapses 
Description Cylinder radius r 
 Laminar, defining upper/lower boundaries of layers 1, 2/3, 4, 5, 6 
Populations 
Type Each cell type y assigned to population Y, yY 
Description Populations Y{L2/3,L4,L5,L6}×{E,I} (population size NY, cell types yY
 (e.g., L4E={p4,ss4(L23),ss4(L4)}, cf., Fig. 2). 
 Cell types y
 •Size Ny=FyYNY, Fy is the occurrence of cell type y in the full model 
 •Morphology My 
 •Extrapolated according to spatial connectivity data (Table 7
 Somatic placement, population Y
 •Random soma placement in cylindrical volumes with radius r, thickness h 
 •Volumes centered between boundaries of layers 2/3–6 
Morphologies 
Type 3D histological reconstructions from slice preparations (see Jacobs et al. 2009; De Schutter and Van Geit 2009) of cat visual and somatosensory cortices 
Description One morphology My per cell type: 
 •Excitatory and inhibitory cells in layers 2/3–6 
 •For all cells jy:Mj = My 
 •For some cell types y,y:My=My (limited availability) 
 Orientations: 
 •Pyramidal cells: apical dendrites oriented along depth axis with random depth-axis rotation 
 •Interneurons, stellate cells: random rotation around all axes 
 Corrections: 
 •Apical dendrites of pyramidal cells elongated to accommodate spatial connectivity 
 •Axons removed if present 
 Reconstructed morphologies (cf., Fig. 2): 
 •Cat visual cortex (Kisvárday and Eysel 1992; Mainen and Sejnowski 1996; Contreras et al. 1997; Stepanyants et al. 2008
 •Cat somatosensory cortex from NeuroMorpho.org (Contreras et al. 1997; Ascoli et al. 2007). 
Table 4

Description of multicompartment-neuron populations for the cortical microcircuit model (continuation of Table 3)

Neuron models 
Type Passive, multicompartment, reconstructed morphologies 
Description Compartment n membrane potential Vmjn of cell j having length ljn, diameter djn and surface area Ajn:  
(22)
CmjndVmjndt=k=1mIajknGLjn(VmjnEL)iIjin,
 
(23)
Cmjn=cmAjn,
 
(24)
Iajkn=Gajkn(VmjkVmjn),
 
(25)
Gajkn=π(djk2+djn2)/(4ra(ljk+ljn)),
 
(26)
GLjn=Ajn/rm,
 
(27)
Imjn=CmjndVmjndt+GLjn(VmjnEL)+iIjin.
Cmjn is compartment capacitance, GLjn its passive leak conductance, EL the passive leak reversal potential, Iajkn axial current between compartment n and neighboring compartment k (out of m compartments), Gajkn axial conductance between n and k, Ijin synaptic currents, and Imjn transmembrane current of compartment n. For specific parameter values, see Table 6. Membrane potentials and transmembrane currents are computed using
NEURON
through
LFPy
(Carnevale and Hines 2006; Lindén et al. 2014), assuming the extracellular potential to be zero everywhere on the outside of the neuron, that is, an infinite extracellular conductivity. 
Synapse model 
Type Exponential postsynaptic current, static weights 
Description Neuron j input current of synapse formed with presynaptic neuron i:  
(28)
Iji(t)=Iji,maxlexp((ttlidi)/τs)H(ttlidi),
 
(29)
Iji,max=CmμYXofpointneuronnetwork,
 
(30)
H(t)=1fort0,and0elsewhere.
•Static synaptic weights Jji=μYX (jY, iX) (see Table 2
 •Delays di from Gaussian distribution with mean dX (iX), relative standard deviation σd,rel 
 •Synapse activation times: network spike trains plus delay 
 •No cortico-cortical connections: Iext=0 (cf., Table 5
Input 
Type Spike times til of spiking neuron network (including thalamic input spikes), no cortico-cortical input 
Description Synapse placement, postsynaptic cell jy,yY (see Methods): 
 •Number of synapses from presynaptic population X in layer L: kyXL (Eq. 9
 •Compartment specificity of connections: Ajn/nLAjn, compartment nL 
 •Synapse locations within layers are chosen randomly among dendritic compartments only 
Measurements 
Type Local field potential (LFP) and current source density (CSD) 
Description Laminar multi-electrode, see parameter values in Table 6
 •Axis perpendicular to pial surface 
 ncontacts: number of contacts 
 hcontacts: intercontact distance 
 rcontact: contact surface radius 
Neuron models 
Type Passive, multicompartment, reconstructed morphologies 
Description Compartment n membrane potential Vmjn of cell j having length ljn, diameter djn and surface area Ajn:  
(22)
CmjndVmjndt=k=1mIajknGLjn(VmjnEL)iIjin,
 
(23)
Cmjn=cmAjn,
 
(24)
Iajkn=Gajkn(VmjkVmjn),
 
(25)
Gajkn=π(djk2+djn2)/(4ra(ljk+ljn)),
 
(26)
GLjn=Ajn/rm,
 
(27)
Imjn=CmjndVmjndt+GLjn(VmjnEL)+iIjin.
Cmjn is compartment capacitance, GLjn its passive leak conductance, EL the passive leak reversal potential, Iajkn axial current between compartment n and neighboring compartment k (out of m compartments), Gajkn axial conductance between n and k, Ijin synaptic currents, and Imjn transmembrane current of compartment n. For specific parameter values, see Table 6. Membrane potentials and transmembrane currents are computed using
NEURON
through
LFPy
(Carnevale and Hines 2006; Lindén et al. 2014), assuming the extracellular potential to be zero everywhere on the outside of the neuron, that is, an infinite extracellular conductivity. 
Synapse model 
Type Exponential postsynaptic current, static weights 
Description Neuron j input current of synapse formed with presynaptic neuron i:  
(28)
Iji(t)=Iji,maxlexp((ttlidi)/τs)H(ttlidi),
 
(29)
Iji,max=CmμYXofpointneuronnetwork,
 
(30)
H(t)=1fort0,and0elsewhere.
•Static synaptic weights Jji=μYX (jY, iX) (see Table 2
 •Delays di from Gaussian distribution with mean dX (iX), relative standard deviation σd,rel 
 •Synapse activation times: network spike trains plus delay 
 •No cortico-cortical connections: Iext=0 (cf., Table 5
Input 
Type Spike times til of spiking neuron network (including thalamic input spikes), no cortico-cortical input 
Description Synapse placement, postsynaptic cell jy,yY (see Methods): 
 •Number of synapses from presynaptic population X in layer L: kyXL (Eq. 9
 •Compartment specificity of connections: Ajn/nLAjn, compartment nL 
 •Synapse locations within layers are chosen randomly among dendritic compartments only 
Measurements 
Type Local field potential (LFP) and current source density (CSD) 
Description Laminar multi-electrode, see parameter values in Table 6
 •Axis perpendicular to pial surface 
 ncontacts: number of contacts 
 hcontacts: intercontact distance 
 rcontact: contact surface radius 
Table 5

Parameters of the cortical microcircuit model

Global simulation parameters 
Symbol Value Description 
T 5200 ms Simulation duration 
dt 0.1 ms Temporal resolution 
Point-neuron network 
Populations and external input 
Symbol Value Description 
X L23E L23I L4E L4I L5E L5I L6E L6I TC Name 
NX 20,683 5834 21,915 5479 4850 1065 14,395 2948 902 Size 
kXext 1600 1500 2100 1900 2000 1900 2900 2100  Ext. in-degree per neuron 
Iext τsynνbgJ,νbg=8Hz       DC ampl. per ext. input 
Connectivity 
CYX from
  L23E L23I L4E L4I L5E L5I L6E L6I TC 
toY L23E 0.101 0.169 0.044 0.082 0.032 0.0 0.008 0.0 0.0 
 L23I 0.135 0.137 0.032 0.052 0.075 0.0 0.004 0.0 0.0 
 L4E 0.008 0.006 0.050 0.135 0.007 0.0003 0.045 0.0 0.0983 
 L4I 0.069 0.003 0.079 0.160 0.003 0.0 0.106 0.0 0.0619 
 L5E 0.100 0.062 0.051 0.006 0.083 0.373 0.020 0.0 0.0 
 L5I 0.055 0.027 0.026 0.002 0.060 0.316 0.009 0.0 0.0 
 L6E 0.016 0.007 0.021 0.017 0.057 0.020 0.040 0.225 0.0512 
 L6I 0.036 0.001 0.003 0.001 0.028 0.008 0.066 0.144 0.0196 
Connection parameters 
Symbol Value Description 
J 87.81 pA Reference synaptic strength. All synapse weights are measured in units of J
σJ,rel  Relative width of synaptic strength distribution 
 •for lognormal distribution 
 0.1 •for Gaussian distribution 
gYX  Relative synaptic strength: 
 X{TC,L23E,L4E,L5E,L6E}, 
 4 X{L23I,L4I,L5I,L6I}, except for: 
 (X,Y)=(L4E,L23E) 
 4.5 (X,Y)=(L4I,L4E) 
dE 1.5 ms Mean excitatory spike transmission delay 
dI 0.75 ms Mean inhibitory spike transmission delay 
σd,rel 0.5 Relative width (stdev/mean) of transmission delay distributions 
Neuron model 
Symbol Value Description 
Rm 40 MΩ Membrane resistance 
Cm 250 pF Membrane capacitance 
τm RmCm (10 ms) Membrane time constant 
EL −65 mV Resting potential 
θ −50 mV Fixed firing threshold 
Vm(t=0) [65,50]mV Uniformly distributed initial membrane potential 
Vreset EL Reset potential 
τref 2 ms Absolute refractory period 
τsyn 0.5 ms Postsynaptic current time constant 
Thalamocortical (TC) input 
Symbol Value Description 
ν¯TC 30 s−1 Mean firing rate per TC neuron 
ΔνTC 30 s−1 Firing-rate modulation amplitude per TC neuron 
fTC 15 Hz Frequency of sinusoidal firing-rate modulation 
Global simulation parameters 
Symbol Value Description 
T 5200 ms Simulation duration 
dt 0.1 ms Temporal resolution 
Point-neuron network 
Populations and external input 
Symbol Value Description 
X L23E L23I L4E L4I L5E L5I L6E L6I TC Name 
NX 20,683 5834 21,915 5479 4850 1065 14,395 2948 902 Size 
kXext 1600 1500 2100 1900 2000 1900 2900 2100  Ext. in-degree per neuron 
Iext τsynνbgJ,νbg=8Hz       DC ampl. per ext. input 
Connectivity 
CYX from
  L23E L23I L4E L4I L5E L5I L6E L6I TC 
toY L23E 0.101 0.169 0.044 0.082 0.032 0.0 0.008 0.0 0.0 
 L23I 0.135 0.137 0.032 0.052 0.075 0.0 0.004 0.0 0.0 
 L4E 0.008 0.006 0.050 0.135 0.007 0.0003 0.045 0.0 0.0983 
 L4I 0.069 0.003 0.079 0.160 0.003 0.0 0.106 0.0 0.0619 
 L5E 0.100 0.062 0.051 0.006 0.083 0.373 0.020 0.0 0.0 
 L5I 0.055 0.027 0.026 0.002 0.060 0.316 0.009 0.0 0.0 
 L6E 0.016 0.007 0.021 0.017 0.057 0.020 0.040 0.225 0.0512 
 L6I 0.036 0.001 0.003 0.001 0.028 0.008 0.066 0.144 0.0196 
Connection parameters 
Symbol Value Description 
J 87.81 pA Reference synaptic strength. All synapse weights are measured in units of J
σJ,rel  Relative width of synaptic strength distribution 
 •for lognormal distribution 
 0.1 •for Gaussian distribution 
gYX  Relative synaptic strength: 
 X{TC,L23E,L4E,L5E,L6E}, 
 4 X{L23I,L4I,L5I,L6I}, except for: 
 (X,Y)=(L4E,L23E) 
 4.5 (X,Y)=(L4I,L4E) 
dE 1.5 ms Mean excitatory spike transmission delay 
dI 0.75 ms Mean inhibitory spike transmission delay 
σd,rel 0.5 Relative width (stdev/mean) of transmission delay distributions 
Neuron model 
Symbol Value Description 
Rm 40 MΩ Membrane resistance 
Cm 250 pF Membrane capacitance 
τm RmCm (10 ms) Membrane time constant 
EL −65 mV Resting potential 
θ −50 mV Fixed firing threshold 
Vm(t=0) [65,50]mV Uniformly distributed initial membrane potential 
Vreset EL Reset potential 
τref 2 ms Absolute refractory period 
τsyn 0.5 ms Postsynaptic current time constant 
Thalamocortical (TC) input 
Symbol Value Description 
ν¯TC 30 s−1 Mean firing rate per TC neuron 
ΔνTC 30 s−1 Firing-rate modulation amplitude per TC neuron 
fTC 15 Hz Frequency of sinusoidal firing-rate modulation 

The stationary thalamic Poisson input and cortico-cortical input to the microcircuit present in the original model of Potjans and Diesmann (2014) are here replaced by DC currents. DC input slightly increases the degree of synchrony (see e.g., Brunel (2000)), but retains network dynamics and firing rate distributions across populations as in Potjans and Diesmann (2014).

The network of Potjans and Diesmann (2014) shows slightly synchronous behavior due to the E–I network of layer 4 being close to the synchronous irregular (SI) regime (Brunel 2000; Bos et al. 2016). In order to reduce synchrony, we here increased the average synaptic weight from neurons in population L4I (inhibitory) to L4E (excitatory) neurons by 12.5%, resulting in attenuated oscillations in layer 4. Taking advantage of the fact that point-neuron networks are amenable for theoretical analysis, we derived modified weights based on predictions from dynamical mean-field theory applied to the microcircuit model (Bos et al. 2016). Moreover, we found that high-frequency network oscillations seen for Gaussian synaptic weight distributions are reduced when using lognormally distributed synaptic weights (Sarid 2007; Iyer et al. 2013; Teramae and Fukai 2014). This made the dynamics more similar to experimental observations (Song et al. 2005; Buzsáki and Mizuseki 2014), and we thus also chose this for our network. Henceforth, we refer to our modified network as the “reference network”. Modulated activity of each thalamo-cortical (TC) neuron in the external thalamic population was modeled as synchronous spikes or as independent non-stationary Poisson processes with sinusoidally oscillating rate profiles (cf. Eq. (18)).

Populations of Multicompartment Model Neurons

Cancellation effects from positive and negative contributions to extracellular potentials and effects of intrinsic dendritic filtering can only be captured with spatially extended multicompartment neuron models (Einevoll et al. 2013). In the hybrid scheme, extracellular potentials are estimated from the spiking activity in the point-neuron network through synaptic activation of populations of multicompartment model neurons (“LFP generators”). In principle, these mutually unconnected model neurons mirror their network counterparts and receive inputs from exactly the same point neurons.

In addition to the description of the point-neuron network model, different types of spatial information are thus needed to predict LFPs. For one, detailed dendritic morphologies are required for each individual network population (Fig. 2). Furthermore, the positions of neurons and synaptic connections must also be specified, as well as the separation of network populations into morphologically distinct cell types (Figures 2 and 4).

Figure 2.

Cell types and morphologies of the multicompartment-neuron populations. The 8 cortical populations Y of size NY in the microcircuit network model are represented by 16 subpopulations of cell type y with detailed morphologies My (Binzegger et al. 2004, Izhikevich and Edelman 2008). Neuron reconstructions are obtained from cat visual cortex and cat somatosensory cortex (source: NeuroMorpho.org (Ascoli et al. 2007), Contreras et al. (1997), Mainen and Sejnowski (1996), Kisvárday and Eysel (1992), Stepanyants et al. (2008) cf. Table 7). Each morphology My is here shown in relation to the layer boundaries (horizontal lines). Colors distinguish between network populations as in Figure 1. The number of compartments (ncomp), frequencies of occurrence (Fy), relative occurrence (FyY), and cell count (Ny) are given for each cell type yY.

Figure 2.

Cell types and morphologies of the multicompartment-neuron populations. The 8 cortical populations Y of size NY in the microcircuit network model are represented by 16 subpopulations of cell type y with detailed morphologies My (Binzegger et al. 2004, Izhikevich and Edelman 2008). Neuron reconstructions are obtained from cat visual cortex and cat somatosensory cortex (source: NeuroMorpho.org (Ascoli et al. 2007), Contreras et al. (1997), Mainen and Sejnowski (1996), Kisvárday and Eysel (1992), Stepanyants et al. (2008) cf. Table 7). Each morphology My is here shown in relation to the layer boundaries (horizontal lines). Colors distinguish between network populations as in Figure 1. The number of compartments (ncomp), frequencies of occurrence (Fy), relative occurrence (FyY), and cell count (Ny) are given for each cell type yY.

Availability of detailed cell-type specific connectivity of neural circuits, especially including information about synapse positions, is limited due to the substantial experimental effort involved. However, several ongoing large-scale neuroscience projects (Kandel et al. 2013) address this issue, and detailed connectomes are beginning to become publicly available (Jiang et al. 2015; Reimann et al. 2015; Markram et al. 2015). In the present example application, we used the connection probabilities as given by Izhikevich and Edelman (2008) derived from Binzegger et al. (2004). Note that the point-neuron network connectivity was partially derived from the same data (Potjans and Diesmann 2014). Quantitative data were provided for the number of connections in 5 cortical layers (layer 1 (L1), layers 2 and 3 grouped into a joint layer 2/3 (L2/3), and layers 4 (L4), 5 (L5), and 6 (L6)) between 17 cortical cell types, cortico-cortical connections from other areas, and 2 TC relay cell types. We follow the nomenclature of Izhikevich and Edelman (2008), where y=p23 denotes pyramidal cell types in layer 2/3, y=b23 and y=nb23 basket interneurons and non-basket interneurons within the same layer, y=ss4(L23) spiny stellate cells in layer 4 with targets mainly within layer 2/3, y=p4 layer 4 pyramidal cells and so forth. Out of the 17 covered intracortical cell types only the y=nb1 cell type is not associated with any point-neuron network population in our scheme. To account for the lack of layer 1 neurons in our model, we renormalized the connection probabilities for the remaining 16 cortical cell types including the 2 TC relay-cell types, such that the occurrences Fy of all cell types y summed to 100% as given in Table 8.

Furthermore, we assumed that the excitatory point-neuron network populations within one layer are composed of pyramidal cells and spiny stellate cells if both are present in the layer, and that inhibitory network populations encompass both types of interneurons. This results in the grouping of cell types y into postsynaptic populations Y illustrated in Figure 2. The neuron count Ny of each cell type is then trivially computed from the frequency of occurrence Fy as given in Table 8 and Figure 2.

Inclusion of cell-type and layer-specific connections in the present hybrid scheme has some implications for how we proceed with setting up equivalent populations consisting of morphologically detailed model neurons. Different cell types belonging to a particular population may have different spatial distributions of synapses, or the populations may consist of different morphological classes of neurons (Kisvárday and Eysel 1992; Nowak 2003; Stepanyants et al. 2008). An example is layer 4 in which spiny stellate cells lack apical dendrites, while pyramidal cells have apical dendrites extending into layer 1. To incorporate some of this morphological diversity, we considered altogether 16 cell types for the 8 cortical network populations.

For each of the 16 cell types, we acquired representative morphological reconstructions of predominantly cat visual cortex neurons from several sources (Contreras et al. 1997; Kisvárday and Eysel 1992; Mainen and Sejnowski 1996; Stepanyants et al. 2008) (cf. Fig. 2, Table 7). Morphology files were obtained either from NeuroMorpho.org (Ascoli et al. 2007) or through personal communication with the authors. Constrained by layer boundary depths (Stepanyants et al. 2008) and laminar connectivities (see below) we applied an intermediate preprocessing step to our pyramidal cell morphologies. Assuming that the soma compartments of each cell type were centered in their corresponding layer, and noting that the layer-specific connectivity (cf. Table 8) implies connections to layer 1, we stretched the apical dendrites along the axis perpendicular to the cortical surface such that they reached the pial surface. The only exception was the p6(L4) morphology, which we extended to reach the center of layer 2/3 in accordance with Stepanyants et al. (2008) and the observation that Table 8 predicts zero connections within layer 1 and very few connections in layer 2/3 to the p6(L4) cell type. Due to lack of available morphologies of sufficient reconstruction quality, certain cell types were represented by the same neuron morphology. Interneuron types and spiny stellate cells in a given layer shared morphologies, the same interneuron morphology was reused in both layers 5 and 6, and finally the p5(L23) and p6(L56) cell morphologies were similar except for the stretching of the apical dendrites.

Table 7

Morphology types and file names used for each cell type in the model (p — pyramidal cell, ss — spiny stellate, i — interneuron)

Morphology files 
Cell type y Morphology My File Source Online source 
p23 p23 oi24rpy1.hoc (Kisvárday and Eysel 1992#NMO_00851 (#NMO_10045) 
b23 i23 oi38lbc1.hoc (Stepanyants et al. 2008— 
nb23 i23 oi38lbc1.hoc (Stepanyants et al. 2008— 
p4 p4 oi53rpy1.hoc (Kisvárday and Eysel 1992#NMO_00855 (#NMO_10040) 
ss4(L23) ss4 j7_L4ste.hoc (Mainen and Sejnowski 1996#MDB_2488 (#NMO_00905) 
ss4(L4) ss4 j7_L4ste.hoc (Mainen and Sejnowski 1996#MDB_2488 (#NMO_00905) 
b4 i4 oi26rbc1.hoc (Stepanyants et al. 2008— 
nb4 i4 oi26rbc1.hoc (Stepanyants et al. 2008— 
p5(L23) p5v1 oi15rpy4.hoc (Kisvárday and Eysel 1992#NMO_00850 (#NMO_10046) 
p5(L56) p5v2 j4a.hoc (Mainen and Sejnowski 1996#MDB_2488 
b5 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
nb5 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
p6(L4) p6 51–2a.CN.hoc (Contreras et al. 1997#NMO_00879 
p6(L56) p5v1 oi15rpy4.hoc (Kisvárday and Eysel 1992#NMO_00850 (#NMO_10046) 
b6 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
nb6 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
Morphology files 
Cell type y Morphology My File Source Online source 
p23 p23 oi24rpy1.hoc (Kisvárday and Eysel 1992#NMO_00851 (#NMO_10045) 
b23 i23 oi38lbc1.hoc (Stepanyants et al. 2008— 
nb23 i23 oi38lbc1.hoc (Stepanyants et al. 2008— 
p4 p4 oi53rpy1.hoc (Kisvárday and Eysel 1992#NMO_00855 (#NMO_10040) 
ss4(L23) ss4 j7_L4ste.hoc (Mainen and Sejnowski 1996#MDB_2488 (#NMO_00905) 
ss4(L4) ss4 j7_L4ste.hoc (Mainen and Sejnowski 1996#MDB_2488 (#NMO_00905) 
b4 i4 oi26rbc1.hoc (Stepanyants et al. 2008— 
nb4 i4 oi26rbc1.hoc (Stepanyants et al. 2008— 
p5(L23) p5v1 oi15rpy4.hoc (Kisvárday and Eysel 1992#NMO_00850 (#NMO_10046) 
p5(L56) p5v2 j4a.hoc (Mainen and Sejnowski 1996#MDB_2488 
b5 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
nb5 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
p6(L4) p6 51–2a.CN.hoc (Contreras et al. 1997#NMO_00879 
p6(L56) p5v1 oi15rpy4.hoc (Kisvárday and Eysel 1992#NMO_00850 (#NMO_10046) 
b6 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 
nb6 i5 oi15rbc1.hoc (Stepanyants et al. 2008— 

Online source numbers #NMO_* refer to NeuroMorpho.org identifiers, #MDB_* refer to ModelDB identifiers.

Table 8

Spatial distribution of cell types and synapses of the cortical microcircuit model adapted from Binzegger et al. (2004) and Izhikevich and Edelmann (2008). Note that occurrences are renormalized between cell types within layers 2/3, 4, 5 and 6 so that yFy=100%

Postsynaptic Layer Occurr. Num. Presynaptic populations X and cell types x 
pop. cell type   syn. L23E L23I L4E L4I L5E L5I L6E L6I TC 
     p23 b23 nb23 ss4(L4) ss4(L23) p4 b4 nb4 p5(L23) p5(L56) b5 nb5 p6(L4) p6(L56) b6 nb6 TCs TCn 
Y y L Fy (%) kyL pyxL (%) 
L23E p23 2/3 26.7 5800 59.9 9.1 4.4 0.6 6.9 7.7 — 0.8 7.4 — — — 2.3 — — 0.8 — — 
   1306 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L23I b23 2/3 3.2 3854 51.6 10.6 3.4 0.5 5.8 6.6 — 0.8 6.3 — — — 2.1 — — 0.7 — 0.5 
 nb23 2/3 4.3 3307 48.6 11.4 3.3 0.5 5.5 6.2 — 0.8 5.9 — — — 1.8 — — 0.6 — 0.7 
L4E ss4(L4) 9.4 5792 2.7 0.2 0.6 11.9 3.7 4.1 7.1 0.8 0.1 — — 32.7 — — 5.8 1.7 1.3 
 ss4(L23) 9.4 4989 5.6 0.4 0.8 11.3 3.8 4.3 7.2 2.1 1.1 0.1 — — 31.1 — — 5.5 1.7 1.3 
 p4 9.4 5031 4.3 0.2 0.6 11.5 3.6 4.2 7.2 2.1 1.2 0.1 — — 31.4 0.1 — 5.9 1.7 1.3 
  2/3  866 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   806 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L4I b4 5.5 3230 5.8 0.5 0.8 11 3.8 4.2 8.4 2.4 1.1 — — — 30.3 — — 5.4 1.6 1.2 
 nb4 1.5 3688 2.7 0.2 0.6 11.7 3.6 8.2 2.3 0.8 0.1 — — 32.2 — — 5.7 1.7 1.3 
L5E p5(L23) 4.8 4316 45.9 1.8 0.3 3.3 7.5 — 0.9 11.7 0.8 1.1 2.3 2.1 — 11.5 0.1 0.4 
   283 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  412 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   185 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
 p5(L56) 1.3 5101 44.3 1.7 0.2 3.2 7.3 — 0.8 11.3 1.2 0.8 1.1 2.3 2.5 0.3 11.3 0.2 0.5 
   949 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  1367 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   5658 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L5I b5 0.6 2981 45.5 2.3 0.2 3.3 7.5 — 1.1 11.6 0.9 1.3 2.3 — 11.4 0.1 0.4 
 nb5 0.8 2981 45.5 2.3 0.2 3.3 7.5 — 1.1 11.6 0.9 1.3 2.3 — 11.4 0.1 0.4 
L6E p6(L4) 14.0 3261 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.3 1.2 13.2 7.7 7.7 0.6 2.9 
   1066 46.8 0.8 0.3 3.4 2.1 7.7 — 0.6 11.9 0.6 0.8 2.3 2.1 — 11.7 0.1 0.4 
   1915 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  121 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
 p6(L56) 4.6 5573 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.3 1.2 13.2 7.8 7.8 0.6 2.9 
   257 46.8 0.8 0.3 3.4 2.1 7.7 — 0.6 11.9 0.6 0.8 2.3 2.1 — 11.7 0.1 0.4 
   243 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  286 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   62 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L6I b6 2.0 3230 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.4 1.2 13.2 7.7 7.7 0.6 2.9 
 nb6 2.0 3230 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.4 1.2 13.2 7.7 7.7 0.6 2.9 
Postsynaptic Layer Occurr. Num. Presynaptic populations X and cell types x 
pop. cell type   syn. L23E L23I L4E L4I L5E L5I L6E L6I TC 
     p23 b23 nb23 ss4(L4) ss4(L23) p4 b4 nb4 p5(L23) p5(L56) b5 nb5 p6(L4) p6(L56) b6 nb6 TCs TCn 
Y y L Fy (%) kyL pyxL (%) 
L23E p23 2/3 26.7 5800 59.9 9.1 4.4 0.6 6.9 7.7 — 0.8 7.4 — — — 2.3 — — 0.8 — — 
   1306 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L23I b23 2/3 3.2 3854 51.6 10.6 3.4 0.5 5.8 6.6 — 0.8 6.3 — — — 2.1 — — 0.7 — 0.5 
 nb23 2/3 4.3 3307 48.6 11.4 3.3 0.5 5.5 6.2 — 0.8 5.9 — — — 1.8 — — 0.6 — 0.7 
L4E ss4(L4) 9.4 5792 2.7 0.2 0.6 11.9 3.7 4.1 7.1 0.8 0.1 — — 32.7 — — 5.8 1.7 1.3 
 ss4(L23) 9.4 4989 5.6 0.4 0.8 11.3 3.8 4.3 7.2 2.1 1.1 0.1 — — 31.1 — — 5.5 1.7 1.3 
 p4 9.4 5031 4.3 0.2 0.6 11.5 3.6 4.2 7.2 2.1 1.2 0.1 — — 31.4 0.1 — 5.9 1.7 1.3 
  2/3  866 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   806 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L4I b4 5.5 3230 5.8 0.5 0.8 11 3.8 4.2 8.4 2.4 1.1 — — — 30.3 — — 5.4 1.6 1.2 
 nb4 1.5 3688 2.7 0.2 0.6 11.7 3.6 8.2 2.3 0.8 0.1 — — 32.2 — — 5.7 1.7 1.3 
L5E p5(L23) 4.8 4316 45.9 1.8 0.3 3.3 7.5 — 0.9 11.7 0.8 1.1 2.3 2.1 — 11.5 0.1 0.4 
   283 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  412 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   185 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
 p5(L56) 1.3 5101 44.3 1.7 0.2 3.2 7.3 — 0.8 11.3 1.2 0.8 1.1 2.3 2.5 0.3 11.3 0.2 0.5 
   949 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  1367 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   5658 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L5I b5 0.6 2981 45.5 2.3 0.2 3.3 7.5 — 1.1 11.6 0.9 1.3 2.3 — 11.4 0.1 0.4 
 nb5 0.8 2981 45.5 2.3 0.2 3.3 7.5 — 1.1 11.6 0.9 1.3 2.3 — 11.4 0.1 0.4 
L6E p6(L4) 14.0 3261 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.3 1.2 13.2 7.7 7.7 0.6 2.9 
   1066 46.8 0.8 0.3 3.4 2.1 7.7 — 0.6 11.9 0.6 0.8 2.3 2.1 — 11.7 0.1 0.4 
   1915 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  121 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
 p6(L56) 4.6 5573 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.3 1.2 13.2 7.8 7.8 0.6 2.9 
   257 46.8 0.8 0.3 3.4 2.1 7.7 — 0.6 11.9 0.6 0.8 2.3 2.1 — 11.7 0.1 0.4 
   243 2.8 0.1 0.7 12.2 3.8 4.2 5.2 1.5 0.8 0.1 — — 33.7 — — 5.9 1.8 1.4 
  2/3  286 63.1 5.1 4.1 0.6 7.2 8.1 — 0.6 7.8 — — — 2.5 — — 0.8 — — 
   62 6.3 0.1 1.1 — — 0.1 — — 0.1 — — — — — — — — 4.1 
L6I b6 2.0 3230 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.4 1.2 13.2 7.7 7.7 0.6 2.9 
 nb6 2.0 3230 2.5 0.1 0.1 0.7 0.9 1.3 — 0.1 0.1 4.9 — 0.4 1.2 13.2 7.7 7.7 0.6 2.9 

Preserving the laminar cell density under 1 mm2 surface area of the point-neuron network model, we created for each postsynaptic cell type y model populations where somas were assigned random locations in 3D within cylindrical slabs with radius r=564μm and thickness h=50μm, each centered in their respective layer (illustrated in Fig. 1D, see also Fig. 6E). Regardless of the vertical offset of the soma of pyramidal cells, postsynaptic target dendrites were therefore present within the 80μm thick (Stepanyants et al. 2008) uppermost layer 1 except for cell type p6(L4). For simplicity, each cell type was represented by a single reconstructed morphology in the present application. The full specification of the populations is given in Table 3.

Each neuron is modeled using the multicompartmental, passive cable formalism (Rall 1964; 2009; De Schutter and Van Geit 2009), describing the changes in membrane voltage and the associated transmembrane currents throughout all parts of the neuron geometry (Table 4). We used (non-plastic) exponential current-based synapses as in the point-neuron network model (Table 2). Synapse locations were randomly assigned onto cell compartments assuming a probability proportional to the compartment's surface area divided by the total surface area of the same cell within the target layer. Tables 46 summarize parameters relevant for the synapse models, synapse locations, and passive parameters of the multicompartment models.

Table 6

Parameters of the multicompartment model neuron populations and calculations of extracellular potentials

Multicompartment model neurons 
Symbol Value Description 
cm 1.0 μFcm−2 Membrane capacity 
rm τm/cm Membrane resistivity 
ra 150 Ωcm Axial resistivity 
EL EL Passive leak reversal potential 
Vinit EL Membrane potential at t = 0 ms 
λf 100 Hz Frequency of AC length constant 
λd 0.1 Factor for
d
_
lambda
rule (Hines and Carnevale 2001
σe 0.3Sm1 Extracellular conductivity 
r 10002/πμm Population radius 
h 50 μm Soma layer thickness 
ncontact 16 Number of electrode contacts 
helec 100 μm Laminar-electrode intercontact distance 
rcontact 7.5 μm Electrode contact-point radius 
Multicompartment model neurons 
Symbol Value Description 
cm 1.0 μFcm−2 Membrane capacity 
rm τm/cm Membrane resistivity 
ra 150 Ωcm Axial resistivity 
EL EL Passive leak reversal potential 
Vinit EL Membrane potential at t = 0 ms 
λf 100 Hz Frequency of AC length constant 
λd 0.1 Factor for
d
_
lambda
rule (Hines and Carnevale 2001
σe 0.3Sm1 Extracellular conductivity 
r 10002/πμm Population radius 
h 50 μm Soma layer thickness 
ncontact 16 Number of electrode contacts 
helec 100 μm Laminar-electrode intercontact distance 
rcontact 7.5 μm Electrode contact-point radius 

Values for τmandEL are inherited from network parameters in Table 5.

Spatial Synaptic Connectivity

A full description of the connectivity in networks of multicompartment model neurons requires a 3D representation, for example in the form of sparse NX×NY×ncomp matrices of synaptic weights and spike-transmission delays between presynaptic neurons i[1,NX] and compartments n[1,ncomp] of postsynaptic cells j[1,NY]. Here, NX and NY denote the number of presynaptic and postsynaptic neurons in populations X and Y, respectively, and ncomp the number of compartments of the postsynaptic cell.

In point-neuron networks, in contrast, connectivity is by definition only 2-dimensional (2D) as the cell morphology is collapsed into a single point and, consequently, the specificity of synapse locations on the postsynaptic morphology is ignored. In the proposed hybrid modeling scheme, the connectivity within the point-neuron network is consistent with the connectivity between point neurons and multicompartment model neurons. Ideally, each multicompartment model neuron has its equivalent in the point-neuron network and receives inputs from exactly the same presynaptic sources as its point-neuron counterpart. Synapses should be positioned on the dendritic tree according to anatomical data, and synaptic weights and time constants should be adapted such that the somatic membrane potential or somatic current match the point-neuron counterparts. Such mapping between point neurons and passive multicompartment neurons is feasible (Koch and Poggio 1985; Wybo et al. 2013; 2015).

In the current application of the hybrid scheme to the cortical microcircuit model, we make the simplest approximation to the mapping problem and fixed the current amplitudes Iji,max and synaptic time constants as in the network model, with compartment specificity of connections dependent on compartment surface area (see Table 4). We further preserve only the statistics of connections (average number of inputs, distribution of spike-transmission delays) for each pair of presynaptic and postsynaptic neuron populations, exploiting that connections between network populations are drawn randomly with fixed probabilities. Finally, we simplify the positioning of synapses to a layer specificity of connections. The activation times of each synapse are then given by the spike train of a randomly drawn point neuron in the network model, with random delays consistent with the delay distribution in the network (Tables 2 and 4).

We first show how to derive a 2D point-neuron connectivity from a given 3D multicompartment-neuron connectivity and describe the case where the complexity of the point-neuron network is further reduced by pooling cell types. Then we describe the opposite procedure, connecting an existing (published) point-neuron network with a predefined 2D connectivity to a population of multicompartment model neurons such that the resulting 3D connectivity is as consistent as possible with anatomical data sets accounting for the compartment specificity of connections (e.g., the layer specificity of connections as in the anatomical data published by Binzegger et al. 2004). The procedures outlined below allow a reduction of complexity within the point-neuron network while accounting for the full diversity in cell types and synapse locations for multicompartment-neuron populations which is essential for predicting extracellular potentials (Fig. 3).

Figure 3.

Example LFP responses from single-synapse activations of layer 4 neurons. (A) Illustration of the nontrivial relationship between apical synaptic input (red circle) onto a reconstructed morphology (black) of a pyramidal cell in layer 4 and the corresponding extracellular potential. The exponential synaptic input current Ii,j(t) (upper inset) results in deflections in the extracellular potential ϕ(r,t) here shown as time courses at 2 locations in proximity to the input site and the basal dendrites (green and blue circles, respectively; lower inset). The color-coded isolines show the magnitude of the scalar extracellular potential at t = 2 ms (vertical black line in insets) in the vicinity of the cell. (B) Same as in panel A, however, with the synaptic input current relocated to a basal dendrite, resulting in an extracellular potential with a different spatiotemporal signature less dependent on the geometry of the apical dendritic tree. At the location denoted by the blue circle, the extracellular potential changes sign with time due to interactions between signal propagation in the passive model neuron and volume conduction. (C) Same as panels B and C for a spiny stellate cell in layer 4 receiving an excitatory synaptic input on a basal dendrite.

Figure 3.

Example LFP responses from single-synapse activations of layer 4 neurons. (A) Illustration of the nontrivial relationship between apical synaptic input (red circle) onto a reconstructed morphology (black) of a pyramidal cell in layer 4 and the corresponding extracellular potential. The exponential synaptic input current Ii,j(t) (upper inset) results in deflections in the extracellular potential ϕ(r,t) here shown as time courses at 2 locations in proximity to the input site and the basal dendrites (green and blue circles, respectively; lower inset). The color-coded isolines show the magnitude of the scalar extracellular potential at t = 2 ms (vertical black line in insets) in the vicinity of the cell. (B) Same as in panel A, however, with the synaptic input current relocated to a basal dendrite, resulting in an extracellular potential with a different spatiotemporal signature less dependent on the geometry of the apical dendritic tree. At the location denoted by the blue circle, the extracellular potential changes sign with time due to interactions between signal propagation in the passive model neuron and volume conduction. (C) Same as panels B and C for a spiny stellate cell in layer 4 receiving an excitatory synaptic input on a basal dendrite.

Figure 4.

Constructing spatial synaptic connectivity for the cortical microcircuit model. (A) Illustration of pooling of presynaptic cell types. Presynaptic populations X in the point-neuron model (left box; here X= L4E) consist of multiple cell types x (here x{p4,ss4(L4),ss4(L23)}). The layer-specific number of synapses kyXL (dash-dotted lines) formed between one cell of postsynaptic cell type y (right part of panel A: morphology projected onto cortical layers 1–6; here y = p5(L56)) and a presynaptic population X is given by the sum of all individual cell-type resolved synapse counts kyxL (dotted or dash-dotted lines). (B) Bi-directional cell- and layer-specific pooling and dispersing of synapses between presynaptic and postsynaptic cell types. Postsynaptic populations Y (right box; here Y = L5E) in the point-neuron model consist of multiple cell types y (here y{p5(L56),p5(L23)}). A given presynaptic population X (left box; here X = L4E) containing cell types x (here x{p4,ss4(L4),ss4(L23)}) forms cell-type and layer-specific connections within Y (black connection tree). For the number of synapses KyXL between population X and cells of type y in layer L (right-most branching of connection tree) the synapse count KYX between all cells in X and Y can be obtained by pooling all synapses onto cell types yY and input layers L. Conversely, for a given total number of synapses KYX between all cells in X and Y, the number of synapses KyXL onto a specific cell type y and layer L can, as described by Equation (9), be obtained by calculating the cell-type and layer specificity of connections TyX and LyXL (see Fig. 5) from anatomical data (Table 7).

Figure 4.

Constructing spatial synaptic connectivity for the cortical microcircuit model. (A) Illustration of pooling of presynaptic cell types. Presynaptic populations X in the point-neuron model (left box; here X= L4E) consist of multiple cell types x (here x{p4,ss4(L4),ss4(L23)}). The layer-specific number of synapses kyXL (dash-dotted lines) formed between one cell of postsynaptic cell type y (right part of panel A: morphology projected onto cortical layers 1–6; here y = p5(L56)) and a presynaptic population X is given by the sum of all individual cell-type resolved synapse counts kyxL (dotted or dash-dotted lines). (B) Bi-directional cell- and layer-specific pooling and dispersing of synapses between presynaptic and postsynaptic cell types. Postsynaptic populations Y (right box; here Y = L5E) in the point-neuron model consist of multiple cell types y (here y{p5(L56),p5(L23)}). A given presynaptic population X (left box; here X = L4E) containing cell types x (here x{p4,ss4(L4),ss4(L23)}) forms cell-type and layer-specific connections within Y (black connection tree). For the number of synapses KyXL between population X and cells of type y in layer L (right-most branching of connection tree) the synapse count KYX between all cells in X and Y can be obtained by pooling all synapses onto cell types yY and input layers L. Conversely, for a given total number of synapses KYX between all cells in X and Y, the number of synapses KyXL onto a specific cell type y and layer L can, as described by Equation (9), be obtained by calculating the cell-type and layer specificity of connections TyX and LyXL (see Fig. 5) from anatomical data (Table 7).

Construction of Point-Neuron Network Connectivity

For our example point-neuron network model, the cortical microcircuit model by Potjans and Diesmann (2014), the connectivity is to a large extent based on anatomical data from cat visual cortex (Binzegger et al. 2004; Izhikevich and Edelman 2008) (cf. Table 8).

From Table 8 we obtain 1) the number Ny of neurons belonging to cell type y, 2) the average total number kyL of synapses on all compartments in layer L (input layer) of a single postsynaptic neuron of type y, and 3) the fraction pyxL of the kyL synapses formed with presynaptic neurons of cell type x. The quantity  

(1)
kyxL=pyxLkyL
defines the number of synapses between all presynaptic cells of type x and a single postsynaptic cell of type y in input layer L (cf. network connectivity in Izhikevich and Edelman 2008). The number of synapses between all neurons in x and all neurons in y, irrespective of the input layer L, is given by  
(2)
Kyx=NyLkyxL.

The number Kyx of connections in combination with a chosen connectivity model (e.g., random graphs with binomially distributed (Erdős and Rényi 1959) or fixed in-/out-degree (Newman 2003) or random graphs with defined higher order statistics (Song et al. 2005; Zhao et al. 2011)) is sufficient for setting up the point-neuron network. Assuming independently drawn synapses (allowing multiple connections between neurons), the probability Cyx of at least one connection between a neuron of type x and a neuron of type y can be obtained from Kyx as (Potjans and Diesmann 2014) 

(3)
Cyx=1(11NxNy)Kyx.

In our case, the point-neuron microcircuit model consists of excitatory and inhibitory populations X,Y (see Tables 1 and 2) pooling different presynaptic and postsynaptic cell types xX and yY (cf. Fig. 4). Given a single multicompartment model neuron of type y, we compute the number kyXL of incoming connections (in-degree) from cell types x in each presynaptic population X in a given layer L by pooling all connections as illustrated in Figure 4A as  

(4)
kyXL=xXkyxL.

The total number of connections onto postsynaptic cells y from cells in X is then  

(5)
KyXL=NykyXL.
The layer-specific connection probability CyXL (Fig. 5B) can be derived from Equation (5) analogous to Equation (3) for a presynaptic population size NX (here, NX=xXNx).

Figure 5.

Connectivity of the cortical microcircuit model. (A) Connection probability CYX between presynaptic population X and postsynaptic population Y of the cortical microcircuit model by Potjans and Diesmann (2014) given in Table 5. Zero values are shown as gray here and in subsequent panels. (B) Layer- and cell-type specific connectivity map CyXL, where X, y, and L denote presynaptic populations, postsynaptic cell types, and the synapse location (layer), respectively. This map is computed from the connectivity of the point-neuron network (panel A), cell-type (panel C), and layer specificity (panel D) of connections. (C) Cell-type specificity TyX of connections quantified as the fraction of synapses between presynaptic and postsynaptic populations X and Y formed with a specific postsynaptic cell type y. (D) Layer specificity LyXL of connections denoting the fraction of synapses between population X and cell type y formed in a particular layer L. Both TyX and LyXL in panels C and D, respectively, are calculated from anatomical data (Binzegger et al. 2004, Izhikevich and Edelman 2008), cf. Table 8.

Figure 5.

Connectivity of the cortical microcircuit model. (A) Connection probability CYX between presynaptic population X and postsynaptic population Y of the cortical microcircuit model by Potjans and Diesmann (2014) given in Table 5. Zero values are shown as gray here and in subsequent panels. (B) Layer- and cell-type specific connectivity map CyXL, where X, y, and L denote presynaptic populations, postsynaptic cell types, and the synapse location (layer), respectively. This map is computed from the connectivity of the point-neuron network (panel A), cell-type (panel C), and layer specificity (panel D) of connections. (C) Cell-type specificity TyX of connections quantified as the fraction of synapses between presynaptic and postsynaptic populations X and Y formed with a specific postsynaptic cell type y. (D) Layer specificity LyXL of connections denoting the fraction of synapses between population X and cell type y formed in a particular layer L. Both TyX and LyXL in panels C and D, respectively, are calculated from anatomical data (Binzegger et al. 2004, Izhikevich and Edelman 2008), cf. Table 8.

In order to obtain the connectivity within the point-neuron network, that is, between populations X and Y, we also need to pool over all synapses of input layers L and cell types y within the postsynaptic population Y (dashed/dotted lines in Fig. 4B). Thus  

(6)
KYX=yYKyX=yYLKyXL,
which yields the connectivity of the simplified network structure CYX (cf. Eq. 3, Fig. 5A).

Figure 6.

Overview of output signals obtained from application of the hybrid scheme to a cortical microcircuit (spontaneous activity). Point-neuron network: (A) Spiking activity. Each dot represents the spike time of a point neuron (color coding as in Figure 1). (B) Population-averaged firing rates for each population. (C) Population-averaged somatic input currents (red: excitatory, blue: inhibitory, black lines: total). (D) Population-averaged somatic voltages. Averaged somatic input currents and voltages are obtained from 100 neurons in each population. Multicompartment model neurons: (E) Somas of excitatory (triangles) and inhibitory (stars) multicompartment cells and layer boundaries (gray/black ellipses). Illustration of a laminar electrode (gray) with 16 recording channels (black circles). (F) Depth-resolved CSD obtained from summed transmembrane currents in cylindrical volumes centered at each contact. (G) Depth-resolved LFP calculated at each electrode contact from transmembrane currents of all neurons in the column. Channel 1 is at pial surface, channel 2 at 100μm depth, etc.

Figure 6.

Overview of output signals obtained from application of the hybrid scheme to a cortical microcircuit (spontaneous activity). Point-neuron network: (A) Spiking activity. Each dot represents the spike time of a point neuron (color coding as in Figure 1). (B) Population-averaged firing rates for each population. (C) Population-averaged somatic input currents (red: excitatory, blue: inhibitory, black lines: total). (D) Population-averaged somatic voltages. Averaged somatic input currents and voltages are obtained from 100 neurons in each population. Multicompartment model neurons: (E) Somas of excitatory (triangles) and inhibitory (stars) multicompartment cells and layer boundaries (gray/black ellipses). Illustration of a laminar electrode (gray) with 16 recording channels (black circles). (F) Depth-resolved CSD obtained from summed transmembrane currents in cylindrical volumes centered at each contact. (G) Depth-resolved LFP calculated at each electrode contact from transmembrane currents of all neurons in the column. Channel 1 is at pial surface, channel 2 at 100μm depth, etc.

From Pooled to Specific Network Connectivity

In case of an already existing point-neuron network model, the reverse task of creating a spatial connectivity CyXL from a given point-neuron network connectivity CYX is necessary. This inverse procedure compared with pooling over cell types and input layers entails introducing the cell-type specificity  

(7)
TyX=KyXKYX,
which describes the fraction of synapses between populations X and Y that are formed with a specific postsynaptic cell type y (Fig. 5C), and the layer specificity of connections  
(8)
LyXL=KyXLKyX,
denoting the fraction of synapses between population X and all cells of cell type y formed in a particular layer L (Fig. 5D). The product TyXLyXL defines the probability of a synapse between populations X and Y formed with a specific postsynaptic cell type y in a particular layer L (Fig. 4B). Thus, if KYX is given, the total number of connections in layer L onto postsynaptic cells y from cells in X is  
(9)
KyXL=KYXTyXLyXL.

If KYX is constructed from the same data as TyX and LyXL, Equations (79) are fully consistent. However, KYX can also be computed from any given point-neuron network connectivity CYX. This is particularly relevant for the network connectivity CYX (Fig. 5A) of Potjans and Diesmann (2014) that includes additional data sets for which spatial information on synapse locations is not available. Here the number of synapses kyXL from population X established in layer L on each multicompartment model neuron of type y is obtained from Equation (9) as kyXL=KyXL/Ny.

Forward Modeling of Extracellular Potentials

The LFP signal reflects transmembrane currents weighted according to the distance from the source to the measurement location (Einevoll et al. 2013), and here we compute the LFP from the model neurons using a now well-established forward-modeling scheme combining multicompartment neuron modeling and electrostatic (volume-conduction) theory (Holt and Koch 1999; Gold et al. 2006; Pettersen et al. 2008; Lindén et al. 2010; 2011; Reimann et al. 2013; Lindén et al. 2014; Tomsett et al. 2014).

Each morphology was spatially discretized into compartments using the

d
_
lambda
rule (Hines and Carnevale 2001) with electrotonic length constants computed at f = 100 Hz.

In this forward modeling scheme, localized synaptic activation of a morphologically detailed neuron results in spatially distributed transmembrane currents across the neuronal membrane as calculated using standard cable theory, see, for example, De Schutter and Van Geit (2009). The extracellular potentials, including the LFP, are in turn given as a weighted sum of transmembrane currents as described by volume conduction theory (Holt and Koch 1999; Einevoll et al. 2013).

The cable-equation description is summarized in box D in Table 4. Equation (27) relates synaptic-input currents Ijin onto compartment n in a neuron j from presynaptic neurons i, the membrane voltages Vmjn and transmembrane currents Imjn, and is derived from the assumption (Kirchoff's law) of current balance in the intracellular node of the equivalent electrical circuit of a cylindrical compartment n with m neighboring compartments.

We use the standard convention that a positive membrane current is a positive current from the intracellular to the extracellular space across the membrane. Imjn is assumed to be homogeneously distributed across the outer surface of the cylindrical compartment, and the calculation assumes that the electrical potential on the outside boundary of the membrane is zero at all times. Hence, there are no mutual interactions (i.e., ephaptic coupling (Anastassiou et al. 2011)) between the extracellular potential estimated using volume conduction theory, the transmembrane currents, and intracellular potentials.

The associated extracellular potential resulting from the transmembrane currents is calculated based on volume conduction theory (Nunez and Srinivasan 2006; Einevoll et al. 2013). In the present application where the signal frequencies are well below 1000 Hz, this calculation is simplified by applying the quasistatic approximation to Maxwell's equations, that is, terms with time derivatives of the electrical and magnetic fields are omitted, cf. (Hämäläinen et al. 1993). Furthermore, we assume the extracellular medium to be linear, isotropic, homogeneous, and ohmic (Pettersen et al. 2012; Einevoll et al. 2013) and represented by a scalar extracellular conductivity σe.

Given a time-varying point current source with magnitude I(t) at position r, the scalar extracellular potential ϕ(r,t) at position r and time t is then given by (Nunez and Srinivasan 2006; Lindén et al. 2014) 

(10)
ϕ(r,t)=I(t)4πσerr.

Contributions to the extracellular potential from multiple current sources, that is, transmembrane currents of all individual compartments n from all cells j in a population of N cells sum linearly. In accordance with the assumed homogeneous current distribution along each cylindrical compartment, the line-source approximation is used for dendritic compartments (Holt and Koch 1999). The line-source forward-modeling formula is obtained by integrating Equation (10) along the cylindrical axis of each compartment n, and summing the contributions from all ncomp compartments (Holt and Koch 1999; Pettersen et al. 2008; Lindén et al. 2014):  

(11)
ϕ(r,t)=j=1Nn=1ncompImjn(t)4πσe1rrjndrjn.

Presently, we approximate the thick soma compartments as spherical current sources, and thus combine the point-source equation (Eq. 10) with the line-source formula (Eq. 11) for dendrite compartments, obtaining (Lindén et al. 2014):  

(12)
ϕ(r,t)=j=1N14πσe(Imj,soma(t)rrj,soma+n=2ncompImjn(t)rrjndrjn)=j=1N14πσe(Imj,soma(t)rrj,soma+n=2ncompImjn(t)Δsjnlnhjn2+rjn2hjnljn2+rjn2ljn).

Here, Δsjn denotes compartment length, rjn the perpendicular distance from the electrode point contact to the axis of the line compartment, hjn the longitudinal distance measured from the start of the compartment, and ljn=Δsjn+hjn the longitudinal distance from the other end of the compartment. If the distance between electrode contacts and dendritic current sources becomes smaller than the radius of the dendritic segment, an unphysical singularity in our extracellular potential may occur. In these cases, singularities are avoided by setting rrj,soma or rjn equal to the compartment radius.

Electrode contacts of real recording devices have finite spatial extent and are not point contacts as assumed above. However, the recorded signal can be well approximated as the mean of the potential averaged across the uninsulated surface (Robinson 1968; Nelson et al. 2008; Nelson and Pouget 2010; Ness et al. 2015), at least for current sources positioned further away than an electrode radius or so (Ness et al. 2015). Here we employed the disc-electrode approximation to the potential (Camuñas Mesa and Quiroga 2013; Lindén et al. 2014; Ness et al. 2015):  

(13)
ϕdisc(r,t)=1ASSϕ(r,t)d2r1mh=1mϕ(rh,t).

We further considered circular electrode contacts with a radius of rcontact=7.5μm, and we averaged the point-contact potential in Equation (12) over m = 50 random locations rh across the contact surface S, AS being the surface area. The chosen locations were distributed with uniform probability on circular discs representing each contact surface, with surface vectors oriented perpendicular to the electrode axis (Lindén et al. 2014). Calculations of extracellular potentials were facilitated by

LFPy
(http://LFPy.github.io, Lindén et al. 2014), in which
NEURON
simulation software is used for calculations of transmembrane currents (i.e., solving Eq. (27)) (Hines and Carnevale 2001; Hines et al. 2009; Carnevale and Hines 2006).

Data Analysis and Software

Model Measurements

The main simulation output of the hybrid scheme consists of spike trains of each neuron in the point-neuron network, “ground-truth” CSD and LFP of each neuron in the morphologically detailed postsynaptic model populations (see Tables 3 and 4). Here the term “ground-truth” refers to the fact that the CSD is computed from transmembrane currents rather than estimated from the LFP.

As the transmembrane current of each compartment (Eq. 27) is known at each simulation time step, we follow the procedure of Pettersen et al. (2008) to compute the “ground-truth” CSD in addition to the LFP. From N model neurons with ncomp compartments having membrane currents Imjn(t) and lengths Δsjn, we calculate the CSD signals ρ(r,t) inside cylinder elements Vρ(r) centered around each electrode contact as  

(14)
ρ(r,t)=1πr2helecj=1Nn=1ncompImjn(t)Δsjn,inside(r)Δsjn.
Δsjn,inside(r) denotes the length of the line source contained within Vρ(r). In contrast to Pettersen et al. (2008), we do not apply a spatial filter to the CSD. The volumes have radii equal to the population radius r and heights equal to the electrode separation helec (cf. Table 6).

In the present example application, the extracellular potential is computed at locations corresponding to a laminar multi-electrode array with 16 recording electrodes with an interelectrode distance of helec=100μm, positioned at the cylindrical axis of the model column with the topmost contact at the pial surface (cf. Fig. 1D, see Table 6 for details). Each electrode contact is set to have a radius of 7.5μm (cf. Eq. (13)). In the network, we also record membrane voltages and input currents from a subset of cells in each of the eight cortical populations (see Table 9).

Table 9

Measurements and derived signals obtained from the hybrid scheme for LFP simulations

Measurements and derived signals 
Measurements 
Symbol Description Number of recorded units 
Iiex(t) Excitatory synaptic input current of neuron i 100 per population X 
Iiin(t) Inhibitory synaptic input current of neuron i 100 per population X 
Vi(tSomatic voltage of neuron i 100 per population X 
si(tSpike train of neuron i NX neurons 
ϕi(r,t) Single-cell LFP generated by neuron i NX neurons 
ρi(r,t) Single-cell CSD generated by neuron i NX neurons 
Derived signals 
Symbol Definition Description 
Ii(tIiex(t)+Iiin(t) Total synaptic input current of neuron i 
I¯Xex(t) 1naviXnavIiex(t) Average excitatory synaptic input current of population X (nav=100
I¯Xin(t) 1naviXnavIiin(t) Average inhibitory synaptic input current of population X (nav=100
I¯X(t) 1naviXnavIi(t) Average total synaptic input current of population X (nav=100
V¯X(t) 1naviXnavVi(t) Average membrane voltage of population X (nav=100
νX(t) nXs(t)tbin Instantaneous population (X) firing rate, with nXs(t) being the number of spikes in (t,t+tbin) of all cells in population X, tbin=1ms 
ν¯X(t) νX(t)NX Average instantaneous firing rate of population X 
ϕX(r,t) iXϕi(r,t) Population LFP of population X 
ρX(r,t) iXρi(r,t) Population CSD of population X 
ϕ(r,t) XϕX(r,t) Compound LFP of all cells 
ρ(r,t) XρX(r,t) Compound CSD of all cells 
Rescaled signals 
ϕγ(r,t) XiXXϕi(r,t) Compound LFP signal from subset of neurons iXX with NX=γNX, γ[0,1] 
ϕγξ(r,t) ξϕγ(r,t)(31) Compound LFP signal from subset of neurons iXX with NX=γNX, γ[0,1], rescaled by factor ξ 
Measurements and derived signals 
Measurements 
Symbol Description Number of recorded units 
Iiex(t) Excitatory synaptic input current of neuron i 100 per population X 
Iiin(t) Inhibitory synaptic input current of neuron i 100 per population X 
Vi(tSomatic voltage of neuron i 100 per population X 
si(tSpike train of neuron i NX neurons 
ϕi(r,t) Single-cell LFP generated by neuron i NX neurons 
ρi(r,t) Single-cell CSD generated by neuron i NX neurons 
Derived signals 
Symbol Definition Description 
Ii(tIiex(t)+Iiin(t) Total synaptic input current of neuron i 
I¯Xex(t) 1naviXnavIiex(t) Average excitatory synaptic input current of population X (nav=100
I¯Xin(t) 1naviXnavIiin(t) Average inhibitory synaptic input current of population X (nav=100
I¯X(t) 1naviXnavIi(t) Average total synaptic input current of population X (nav=100
V¯X(t) 1naviXnavVi(t) Average membrane voltage of population X (nav=100
νX(t) nXs(t)tbin Instantaneous population (X) firing rate, with nXs(t) being the number of spikes in (t,t+tbin) of all cells in population X, tbin=1ms 
ν¯X(t) νX(t)NX Average instantaneous firing rate of population X 
ϕX(r,t) iXϕi(r,t) Population LFP of population X 
ρX(r,t) iXρi(r,t) Population CSD of population X 
ϕ(r,t) XϕX(r,t) Compound LFP of all cells 
ρ(r,t) XρX(r,t) Compound CSD of all cells 
Rescaled signals 
ϕγ(r,t) XiXXϕi(r,t) Compound LFP signal from subset of neurons iXX with NX=γNX, γ[0,1] 
ϕγξ(r,t) ξϕγ(r,t)(31) Compound LFP signal from subset of neurons iXX with NX=γNX, γ[0,1], rescaled by factor ξ 

We ran our simulations for a total duration of T=5200 ms using a temporal resolution of dt=0.1 ms (cf. Tables 5 and 11). However, LFP and CSD signals were resampled prior to file storage to a temporal resolution of dtψ=1ms by 1) applying a fourth-order Chebyshev type I low-pass filter with critical frequency fc=400Hz and 0.05 dB ripple in the passband using a forward–backward linear filter operation, and 2) then selecting every 10th time sample.

Post-processing and Data Analysis

As the contributions to the CSD and LFP of the different cells sum linearly (cf. Eq. 12), we compute population-resolved signals as the sum over contributions from all cells in a population, and the full compound signals as the sum over all population signals (see Table 9). In Results, we derive a rescaled “low-density predictor” ϕγξ(r,t) of the LFP from random subsets of neurons in all populations. Thereby, we make a downscaled LFP-generating model setup with the same column volume, but with neuron density reduced to a factor γ(0,1) of the original density, while preserving the in-degrees, that is, the number of synaptic connections onto individual neurons. The LFP from the downscaled setup is multiplied by an overall scaling factor ξ chosen to roughly preserve the LFP from the full-scale model. For analysis and plotting, the initial 200 ms of results after simulation onset was removed, and the signal mean was subtracted from LFP and CSD traces emulating DC filtering during experimental data acquisition. The contribution from each population to the overall LFP signal was assessed by computing and comparing the signal variances (see Table 10).

Table 10

Data analysis of model output signals

Data analysis 
Symbol Definition Description 
ψ,ψ ψ,ψ{ϕi,ϕX,ϕ,ρi,ρX,ρ,νX} Signal (LFPs, CSDs, firing rates) 
μψ(r) dtTh=1T/dtψ(r,hdt) Temporal mean of signal ψ(r,t) 
covψψ(r) dtTh=1T/dtψ(r,hdt)ψ(r,hdt)μψ(r)μψ(r) Temporal covariance of signals ψ(r,t),ψ(r,t) 
σψ2(r) covψψ(r) Temporal variance of signal ψ(r,t) 
ccψψ(r)  
(32)
covψψ(r)/σψ2(r)σψ2(r)
 
Zero time-lag correlation coefficient of signals ψ(r,t),ψ(r,t) 
Cψψ(r,f) F[ψ](r,f)*F[ψ](r,f) (implemented using Welch's method) Pairwise cross-spectral density of signals ψ(r,t),ψ(r,t), F[ψ]:Fourier transform of ψ 
Pψ(r,f) Cψψ(r,f) PSD of signal ψ(r,t) 
Pϕ¯(r,f)  
(33)
1Ni=1NPϕi(r,f)
 
Average single-cell LFP power spectrum 
Cϕ¯(r,f)  
(34)
1N(N1)i=1Nj=1jiNCϕiϕj(r,f)
 
Average cross-spectrum between single-cell LFPs 
κϕ¯(r,f)  
(35)
Cϕ¯(r,f)/Pϕ¯(r,f)
 
Average LFP coherence between cells 
Pϕ(r,f)  
(36)
NPϕ¯(r,f)+N(N1)Cϕ¯(r,f)
 
PSD of the compound LFP 
Pϕ0(r,f)  
(37)
NPϕ¯(r,f)
 
PSD of the compound LFP signal omitting pairwise cross-correlations 
Data analysis 
Symbol Definition Description 
ψ,ψ ψ,ψ{ϕi,ϕX,ϕ,ρi,ρX,ρ,νX} Signal (LFPs, CSDs, firing rates) 
μψ(r) dtTh=1T/dtψ(r,hdt) Temporal mean of signal ψ(r,t) 
covψψ(r) dtTh=1T/dtψ(r,hdt)ψ(r,hdt)μψ(r)μψ(r) Temporal covariance of signals ψ(r,t),ψ(r,t) 
σψ2(r) covψψ(r) Temporal variance of signal ψ(r,t) 
ccψψ(r)  
(32)
covψψ(r)/σψ2(r)σψ2(r)
 
Zero time-lag correlation coefficient of signals ψ(r,t),ψ(r,t) 
Cψψ(r,f) F[ψ](r,f)*F[ψ](r,f) (implemented using Welch's method) Pairwise cross-spectral density of signals ψ(r,t),ψ(r,t), F[ψ]:Fourier transform of ψ 
Pψ(r,f) Cψψ(r,f) PSD of signal ψ(r,t) 
Pϕ¯(r,f)  
(33)
1Ni=1NPϕi(r,f)
 
Average single-cell LFP power spectrum 
Cϕ¯(r,f)  
(34)
1N(N1)i=1Nj=1jiNCϕiϕj(r,f)
 
Average cross-spectrum between single-cell LFPs 
κϕ¯(r,f)  
(35)
Cϕ¯(r,f)/Pϕ¯(r,f)
 
Average LFP coherence between cells 
Pϕ(r,f)  
(36)
NPϕ¯(r,f)+N(N1)Cϕ¯(r,f)
 
PSD of the compound LFP 
Pϕ0(r,f)  
(37)
NPϕ¯(r,f)
 
PSD of the compound LFP signal omitting pairwise cross-correlations 

Cross-correlations between single-cell LFPs ϕi(r,t) were quantified by analyzing the power spectrum of the compound extracellular signals. Given that the compound LFP/CSD is a linear superposition of single-cell LFP/CSD contributions, the power spectrum of the compound signal Pϕ(r,f) can be obtained as the sum of all single-cell power spectra Pϕi(r,f) and all pairwise cross-spectra Cϕiϕj(r,f), or equivalently from the average single-cell power spectrum Pϕ¯(r,f), the average pairwise cross-spectrum Cϕ¯(r,f), and the total cell count N, as shown in Table 10. Note that Cϕ¯(r,f) and hence the average pairwise single-cell LFP coherence κϕ¯(r,f) (Eq. 35) are real, while Cϕiϕj(r,f) is complex. Note that this definition allows for negative values of κϕ¯(r,f). In the sum over all i and j in Equation (34), the imaginary parts of Cϕiϕj(r,f) and Cϕjϕi(r,f) cancel because Cϕiϕj(r,f)=Cϕjϕi(r,f)* (∗ denotes the complex conjugate). The power spectrum Pϕγξ(r,f) of the compound signals of the “downscaled” network (i.e., low-density LFP predictor, see above) is given by the reduced cell count γN, the average single-cell power spectrum Pϕγξ¯(r,f), and the average pairwise single-cell cross-spectrum Cϕγξ¯(r,f) calculated from that subset of neurons (see Table 10). These single-cell averages are the same as the respective single-cell averages of the full-scale model setup, apart from variability due to subsampling.

Throughout this paper, signal power spectra are estimated using Welch's average periodogram method (Welch 1967) (with the

matplotlib.mlab.psd
implementation in
Python
, see Table 11). Temporal cross-correlations are quantified as the zero time-lag correlation coefficient (Eq. 32 in Table 10). Distributions of chance correlation coefficients were determined using surrogate data that preserved the overall spectra. We generated this data by transforming one input vector to the Fourier domain, assigned a random phase [0,2π) to each frequency f>0, and transformed the data back to the time domain prior to computing the correlation coefficients (Schreiber and Schmitz 2000). We show 1% significance levels obtained by computing the 99th percentile of chance correlation coefficients across 1000 trials. Spike-trigger-averaged LFP (stLFP) signals are computed as the cross-covariance between the time-resolved population spike rate νX(t) and the compound LFP ϕ(r,t), divided by the total number of spikes (i.e., 0TνX(t)dt).

Table 11

Post-processing and data analysis parameters

Post-processing 
Parameters 
Symbol Value Description 
Ttrans 200 ms Start-up transient 
dtψ 1 ms Signal resolution 
 
Power spectral density settings 
Method 
plt.mlab.psd
Welch's average periodogram 
Symbol Value Description 
Tψ 5000 ms Signal length 
NFFT 256 Number of data points used in each block for the FFT 
Fs 1 kHz Sampling rate 
noverlap 128 Number of overlapping data points between blocks 
window 
plt.mlab.window
_
hanning
Window filter (*
plt
denotes
matplotlib
Post-processing 
Parameters 
Symbol Value Description 
Ttrans 200 ms Start-up transient 
dtψ 1 ms Signal resolution 
 
Power spectral density settings 
Method 
plt.mlab.psd
Welch's average periodogram 
Symbol Value Description 
Tψ 5000 ms Signal length 
NFFT 256 Number of data points used in each block for the FFT 
Fs 1 kHz Sampling rate 
noverlap 128 Number of overlapping data points between blocks 
window 
plt.mlab.window
_
hanning
Window filter (*
plt
denotes
matplotlib

The
hybridLFPy Python
Package

To facilitate usage of the hybrid scheme by other users, a novel

Python
software package,
hybridLFPy
, has been made publicly available under the General Public License version 3 (GPLv3,http://www.gnu.org/licenses/gpl-3.0.html) on GitHub (http://github.com/INM-6/hybridLFPy). Compatibility with a host of different machine architectures and operating systems (∗nix, OSX, Windows) is ensured with the freely available, object-oriented programming language
Python
(http://www.python.org).
Python
adds tremendous flexibility in terms of interfacing a large number of packages and libraries, for example, for performing numerical analysis and data visualization, such as
numpy
(http://www.numpy.org) and
matplotlib
(http://www.matplotlib.org), while several other neural simulation softwares also come with their own
Python
interfaces, such as
NEST
(http://www.nest-initiative.org, Eppler et al. 2008) and
NEURON
(http://www.neuron.yale.edu, Hines et al. 2009).

The source code release of

hybridLFPy
provides a set of classes implementing the hybrid scheme, as well as example network simulation codes implemented with
NEST
(Eppler et al. 2015) (at present a simplified two-population network (Brunel 2000) and the full cortical microcircuit model of Potjans and Diesmann (2014) adapted from public codes). The class
hybridLFPy.CachedNetwork
uses an efficient
sqlite3
database implementation for reading in all point-neuron network spike events and interfacing network spike events with the main simulation in which the LFP and CSD are calculated. The class
hybridLFPy.Population
defines populations of multicompartment model neurons representing each cell type, assigns synapse locations across the laminae, selects spike trains for each synapse location from the appropriate presynaptic population, and calculates the LFP and CSD. The single-cell calculations are handled using
LFPy
(http://LFPy.github.io) (Lindén et al. 2014), which builds on
NEURON
(Carnevale and Hines 2006; Hines et al. 2009). As there are no mutual interactions between the multicompartment model neurons in the calculation of LFPs, these calculations remain embarrassingly parallel operations. Finally, the class
hybridLFPy.PostProcess
constructs the full compound signals in terms of LFPs and CSDs created by multiple instances of the
Population
class, and performs the main analysis steps as described above. Reproducible simulation and data analysis are assured by tracking code revisions using
git
and by fixing random number generation seeds and the number of parallel processes. Further documentation and information on installing and using
hybridLFPy
is provided online, see http://github.com/INM-6/hybridLFPy.

The present implementation of

hybridLFPy
(v.0.1.3) and corresponding simulations were made possible by
Open MPI
(v.1.6.2),
HDF5
(v.1.8.13),
sqlite3
(v.3.6.20),
Python
(v.2.7.3) with modules
Cython
(v.0.23dev),
NeuroTools
(v.0.2.0dev),
SpikeSort
(v.0.13),
h5py
(v.2.5.0a0),
ipython
(v.0.13),
matplotlib
(v.1.5.x),
mpi4py
(v.1.3),
numpy
(v.1.10.0.dev-c63e1f4),
pysqlite
(v.2.6.3) and
scipy
(v.0.17.0.dev0–357a1a0). Point-neuron network simulations were performed using
NEST
(v.2.8.0 ff71a29), and simulations of multicompartment model neurons using
NEURON
(v.7.4 1186:541994f8f27f) through
LFPy
(dev. v.d393d7). All software was compiled using
GCC
(v.4.4.6). Simulations were performed in parallel (256 threads) on the Stallo high-performance computing facilities (NOTUR, the Norwegian Metacenter for Computational Science) consisting of 2.6 GHz Intel E5–2670 CPUs running the Rocks Cluster Distribution (Linux) operating system (v.6.0).

Results: LFP Generated by a Cortical Microcircuit

To illustrate the application of the hybrid scheme for predictions of LFPs from point-neuron networks, we here present results for a modified version of the point-neuron network model of Potjans and Diesmann (2014) with ∼78, 000 neurons mimicking a 1 mm2 patch of cat primary visual cortex (see Methods). The microcircuit model has realistic cell density and deliberately neglects many biological details on the single-cell level, focusing on the effect of the connectivity on the local network dynamics in such circuits. Despite this simplicity, the model displayed firing rates across populations in agreement with experimental observations (Potjans and Diesmann 2014), as well as propagation of spiking activity across layers. Likewise, the microcircuit model in conjunction with simplified, passive multicompartment populations is used here to study the effect of the (spatial) connectivity on the laminar pattern of spontaneous and stimulus-evoked CSD and LFP signals. For this, CSD and LFP signals for a laminar multi-electrode recording at different cortical depths are computed. The large network size of the model is further used to illustrate the effect of correlations and neuron density on CSD and LFP predictions.

Spontaneous vs. Stimulus-Evoked LFP

We first consider the LFP generated by spontaneous network activity. The output of our hybrid scheme covers various scales and measurement modalities, from spikes of each neuron (Fig. 6A), population-averaged firing rates (Fig. 6B), excitatory, inhibitory, and total synaptic-input currents (Fig. 6C), and membrane voltages (Fig. 6D), to the compound CSD and LFP stemming from all populations of different cell types (Fig. 6EG). For spontaneous activity, that is, no modulated thalamic input (cf. Table 1), we observe asynchronous irregular spiking in all populations (Fig. 6A) and firing rates similar to the original model (Potjans and Diesmann 2014) (Fig. 6B). In particular, layer 2/3 exhibits low firing rates, and in general inhibitory neurons fire with higher rates than excitatory neurons of the respective layers. The network is in a balanced regime (Brunel 2000), reflected by the substantial cancellation of population-averaged excitatory and inhibitory input currents (Fig. 6C). The population-averaged membrane potential fluctuates below the fixed firing threshold θ=50mV down to 80mV (Fig. 6D).

The corresponding compound CSD and LFP signals with contributions from all cortical populations are shown in Figure 6F and G, respectively. As expected for spontaneous cortical network activity, the LFP signal amplitudes are small, 0.1mV (intriguingly close to what has been seen in experiments, see, e.g., Fig. 1 in Maier et al. 2010 for macaque visual cortex, Fig. 7 in Hagen et al. 2015 for mouse visual cortex), and exhibit strong cross-channel covariance and coherence in line with experimental observations, for example Einevoll et al. 2007; Maier et al. 2010; Riehle et al. 2013; Hagen et al. 2015. In the model, these correlations stem from dendritic cable properties and volume conduction effects (Pettersen et al. 2008; Lindén et al. 2011; Łe¸ski et al. 2013). The correlations across channels are, as expected, generally less visible in the more localized CSD signal where volume conduction effects are absent (Nicholson and Freeman 1975; Pettersen et al. 2006; 2008).

LFP and associated CSD studies have commonly been used to investigate stimulus-evoked responses in sensory cortices, see, for example, Mitzdorf and Singer (1979); Mitzdorf (1985); Di et al. (1990); Swadlow et al. (2002); Einevoll et al. (2007); Sakata and Harris (2009); Szymanski et al. (2009); Maier et al. (2010); Jin et al. (2011); Szymanski et al. (2011), as well as Einevoll et al. (2013) and references therein. To model the situation with a sharp onset of a visual stimulus (or direct electrical stimulation of the TC pathway (Mitzdorf and Singer 1979)), we drive the network with a short thalamic pulse mimicking a volley of incoming spikes onto primary visual cortex from the visual thalamus (lateral geniculate nucleus). The activation targets populations in layers 4 and 6 (see Fig. 1A or Fig. 5A,B) and propagates in the network to populations in layers 2/3 and 5 (Fig. 7A,B). At the level of spiking activity, the results match the behavior of the original model (Potjans and Diesmann 2014) and even agree qualitatively with experimental findings in rodents from stimulus-evoked activation of auditory cortex (Sakata and Harris 2009) and somatosensory cortex (Armstrong-James et al. 1992; Einevoll et al. 2007; Reyes-Puerta et al. 2015). This points to a general biological plausibility of this generic network model, based largely on data from cat V1.

Figure 7.

Network activity following transient activation of TC afferents. (A) Raster plot of spiking activity before and after δ-shaped thalamic stimulus presented at t = 900 ms (vertical black line in panels A, C, and D). (B) Population-averaged firing rate histogram for each population (color coding as in panel A). (C) Depth-resolved compound CSD of all populations (shown both in color and by the black traces). (D) Depth-resolved compound LFP (shown both in color and by the black traces) at each electrode channel as generated by all populations. Channel 1 is at pial surface, channel 2 at 100 μm depth, etc.

Figure 7.

Network activity following transient activation of TC afferents. (A) Raster plot of spiking activity before and after δ-shaped thalamic stimulus presented at t = 900 ms (vertical black line in panels A, C, and D). (B) Population-averaged firing rate histogram for each population (color coding as in panel A). (C) Depth-resolved compound CSD of all populations (shown both in color and by the black traces). (D) Depth-resolved compound LFP (shown both in color and by the black traces) at each electrode channel as generated by all populations. Channel 1 is at pial surface, channel 2 at 100 μm depth, etc.

The corresponding CSD and LFP profiles across depth associated with this spiking activity are determined by the synapse locations and dendritic filtering by cell-type specific morphologies (see Figs 3 and 5D). For the CSD (Fig. 7C), a complex alternating spatiotemporal pattern of current sinks (negative CSD) and sources (positive CSD) is observed. Due to volume conduction, this detailed spatial pattern is largely smeared out in the LFP profile (Fig. 7D), which displays a strong positivity across the middle layer around t = 910 ms. We also note the different spatiotemporal profiles of the spiking activity (Fig. 7A,B) compared with laminar CSD and LFP profiles. Not only do the CSD and LFP signals typically fade out 5–10 ms later than the spiking, the spatial profiles are also very different. For example, the LFP signal is very weak between channels 11 and 12 (i.e., between 1100 and 1200μm depth), even if the firing rate of layer 5 positioned between these channels is very high. We also note that the predicted LFP magnitudes are an order of magnitude larger compared with the LFP predicted for spontaneous activity, that is, ∼ 1 mV for stimulus-evoked versus ∼ 0.1 mV for spontaneous activity.

Although the present example has not been tuned to address specific experiments, we nevertheless observe that the model predictions display several features seen in experiments. For example, stimulus-evoked LFP amplitudes on the order of 1 mV are similar to the maximal amplitudes (∼ 1–3 mV) observed in cat V1 following electric stimulation of TC axons (optical radiation) (Mitzdorf 1985) and ventroposterior lateral pathway axons projecting to rat somatosensory cortex (Castro-Alamancos and Connors 1996, Fig. 3), and in rat somatosensory (barrel) cortex following whisker flicks (Di et al. 1990; Einevoll et al. 2007; Reyes-Puerta et al. 2015). Also some qualitative features of the spatiotemporal LFP and CSD patterns from the cat visual cortex experiments of Mitzdorf (1985) can be recognized. One example is the early CSD sink around layer 4 (i.e., at channel 6 close to the boundary between layers 2/3 and 4), another is the large positivity in LFP extending through most of the layers following the initial response to the thalamic input volley.

The observed time courses of the LFP and CSD in our simulations with δ-pulse thalamic activations are not as directly comparable with experimental stimulus-evoked activity, where the input is temporally filtered by several cell populations before reaching thalamus on the way to cortex. However, we note that very swift stimulus-evoked responses lasting not much longer than the ∼ 10 ms response volleys seen in our simulations also are observed in the somatosensory system (Di et al. 1990; Einevoll et al. 2007; Reyes-Puerta et al. 2015).

Effect of Network Dynamics on LFP

The spiking activity in the microcircuit model is highly sensitive to modification of intrinsic model parameters and external input (Bos et al. 2016). In general LFPs reflect synaptic input both from local and distant neurons (Herreras et al. 2015), and also depend on network state (see, e.g., Kelly et al. 2010; Gawne 2010). With the hybrid scheme, we first illustrate as an example the dependence of the spontaneous LFP, that is, the LFP without thalamic input, on local network dynamics as determined by intrinsic network parameters. In particular, we compare our reference network model with the original model proposed by Potjans and Diesmann (2014) where the strength of connections from L4I to L4E neurons is weaker and the synaptic weights are drawn from a Gaussian distribution. We next investigate the LFP when the network is stimulated by sinusoidally modulated thalamic input.

For spontaneous activity (Fig. 8AE), the spiking (Fig. 8A) is asynchronous irregular in all populations. The firing-rate power spectra (Fig. 8B) vary from relatively flat to more band-pass-like with a maximum power at ∼ 80 Hz. The suppressed power at lower frequencies arises from active decorrelation due to inhibitory feedback (Tetzlaff et al. 2012). In combination with the low-pass filtering involved in the generation of LFPs from spiking activity (Lindén et al. 2010; Łe¸ski et al. 2013), the firing-rate spectra translate into an LFP power spectrum that, depending on recording depth, has either low-pass or band-pass filter characteristics (Fig. 8D). The notably sharper attenuation of the LFP power spectra compared with the firing-rate power spectra above 100Hz is expectedly due to the intrinsic dendritic filtering effect (Lindén et al. 2010). As the effect of this filtering depends on the position of the electrode compared with the neuronal morphology (Lindén et al. 2010), the result is a variable LFP power spectral density (PSD) profile across cortical depths, even in the absence of any structured external input (Fig. 8E).

Figure 8.

Effect of network dynamics on LFP. Comparison of two different thalamic input scenarios and two different networks. Top: Reference network, spontaneous activity. Center: Reference network, oscillatory thalamic activation. Bottom: Original model by Potjans and Diesmann (2014), spontaneous activity. (A,F,K) Population-resolved spiking activity. (B,G,L) Population-averaged firing rate spectra. (C,H,M) Depth-resolved LFP. (D,I,N) LFP power spectra in layer 1 and at typical somatic depths of network populations. (E,J,O) LFP power spectra across all channels. Channel 1 is at pial surface, channel 2 at 100 μm depth, etc.

Figure 8.

Effect of network dynamics on LFP. Comparison of two different thalamic input scenarios and two different networks. Top: Reference network, spontaneous activity. Center: Reference network, oscillatory thalamic activation. Bottom: Original model by Potjans and Diesmann (2014), spontaneous activity. (A,F,K) Population-resolved spiking activity. (B,G,L) Population-averaged firing rate spectra. (C,H,M) Depth-resolved LFP. (D,I,N) LFP power spectra in layer 1 and at typical somatic depths of network populations. (E,J,O) LFP power spectra across all channels. Channel 1 is at pial surface, channel 2 at 100 μm depth, etc.

There is experimental evidence that cortical microcircuits receive oscillating input at various frequencies from remote areas and subcortical structures (Bastos et al. 2015, van Kerkoerle et al. 2014; Ito et al. 2014). To illustrate the effect of an oscillatory external input in our model, we model thalamic input as independent realizations of non-stationary Poisson processes with a sinusoidal rate profile (see Table 5, Fig. 8FJ). The spiking activity of all populations as well as the CSD and LFP across depth are sensitive to thalamic input, showing that the network response goes beyond the layers receiving thalamic input, that is, layers 4 and 6. The stimulus-evoked spiking activity follows the 15 Hz modulated rate of the thalamic input (Fig. 8F). The 15 Hz oscillation is reflected in the firing-rate spectrum as a peak at ∼15 Hz (Fig. 8G) and is robustly transferred to the LFP (Fig. 8HJ). The LFP oscillation strength varies with depth and is greatest in channels 1 and 2 and channels 8–14, while the oscillation is barely seen in channels 3–6. Populations L4E/I and L6E/I receive thalamic inputs around the depths of channels 8–10 and channels 14 and 15, and these channels are strongly affected by the stimulus. However, recurrent connections between populations, dendritic propagation of currents, and volume conduction produce strong LFP oscillations also in other channels.

The LFP amplitudes are not only influenced by the temporal structure of the external input, but also by synaptic weights in two ways: first via their influence on the spiking dynamics in the point-network simulation and second via the influence on the size of the synaptic currents setting up the transmembrane currents in the multicompartment models in the LFP-computing step. In order to illustrate the weight dependence, we compare the LFP under spontaneous activity in our reference network model (Fig. 8AE) with the corresponding spontaneous LFP in the original model by Potjans and Diesmann (2014) (Fig. 8K–O). The lower inhibition from population L4I onto L4E and the narrow weight distribution compared with our model gives a higher degree of spike synchrony in all populations. Although our model exhibits asynchronous irregular activity (Fig. 8A), the original model is closer to an SI regime (Brunel 2000) (Fig. 8K), resulting in high-frequency oscillations ∼80 Hz and in the 300–400 Hz band (Fig. 8L), both associated with delay loops in the multilayered network (Bos et al. 2016). The 80 Hz oscillation also appears in the LFP and its corresponding power spectra (Fig. 8M–O), but the magnitude of the peak is not constant across depth. The lowest magnitudes are located in the vicinity of layers 2/3 and 5 (in channels 3 and 4 and channels 11–13).

Contributions from Individual Populations to CSD and LFP

The direct interpretation of CSD and LFP signals in terms of the underlying activity of different populations or input pathways is inherently ambiguous and thus difficult: for example, a CSD sink observed in cortical layer 2/3 can alternatively stem from excitatory synaptic inputs to the basal dendrites of layer 2/3 cells, similar inputs into the apical dendrites of layer 5 cells, or even return currents from appropriately placed inhibitory inputs onto the same cells (Lindén et al. 2010, Einevoll et al. 2013). Several schemes for decomposition of CSD and LFP data into contributions from cortical populations have thus been proposed: principal component analysis (Di et al. 1990), laminar population analysis (LPA, Einevoll et al. 2007), and independent component analysis (Łe¸ski et al. 2010; Makarov et al. 2010; Głąbska et al. 2014; Herreras et al. 2015). In our model, we have the benefit of having the contributions from the various populations, connections, or different synapse types to the CSD and LFP signals directly accessible.

Here, we focus on the LFP and CSD contributions from individual populations and different synapse types. To quantify the contributions from each postsynaptic population, that is, the CSD and LFP stemming from the transmembrane currents of a population, we simply summed all single-cell CSD and LFP contributions from all neurons in the population. Results for spontaneous network activity are shown in Figure 9. As different cell types are assigned appropriate morphologies and cortical depths based on available anatomical data, a trivial consequence is that the neurons in each population make their main contributions to the CSD and LFP at depths spanned by their dendrites (Fig. 9B,C for post-synaptic populations L23E and L6E, respectively). For example, L6E neurons have a high density of afferent synapses on apical dendrites in layer 4, and as a consequence, the amplitude of CSDs and LFPs generated by population L6E (Fig. 9C) is large in the vicinity of layer 4 and not just in layer 6.

Figure 9.

Composition of CSD and LFP during spontaneous activity. (A) Representative morphologies of each population Y illustrating dendritic extent. (B) LFP (black traces) and CSD (color plot) produced by the superficial population L23E for spontaneous activity in the reference network. (C) Similar to panel B for population L6E (summing over contributions of y{p6(L4),p6(L56)}). (D) CSD variance as function of depth for each individual subpopulation (colored lines) and for the full compound signal (thick black line). (E) Same as in panel D, but for LFPs. Variances σ2<107mV2 not shown. (F) Compound LFP (black traces) and CSD (color plot) resulting from only excitatory input to the LFP-generating multicompartment model neurons. (G) Conversely, LFP (black traces) and CSD (color plot) resulting from only inhibitory input to the neurons. (H) Full compound LFP (black traces) and CSD (color plot) resulting from both excitatory and inhibitory synaptic currents. (I) Compound CSD variance as a function of depth with all synapses intact (thick black line), or having only excitatory (red) or inhibitory synapse input (blue). (J) Same as in panel I, but for the LFP signal.

Figure 9.

Composition of CSD and LFP during spontaneous activity. (A) Representative morphologies of each population Y illustrating dendritic extent. (B) LFP (black traces) and CSD (color plot) produced by the superficial population L23E for spontaneous activity in the reference network. (C) Similar to panel B for population L6E (summing over contributions of y{p6(L4),p6(L56)}). (D) CSD variance as function of depth for each individual subpopulation (colored lines) and for the full compound signal (thick black line). (E) Same as in panel D, but for LFPs. Variances σ2<107mV2 not shown. (F) Compound LFP (black traces) and CSD (color plot) resulting from only excitatory input to the LFP-generating multicompartment model neurons. (G) Conversely, LFP (black traces) and CSD (color plot) resulting from only inhibitory input to the neurons. (H) Full compound LFP (black traces) and CSD (color plot) resulting from both excitatory and inhibitory synaptic currents. (I) Compound CSD variance as a function of depth with all synapses intact (thick black line), or having only excitatory (red) or inhibitory synapse input (blue). (J) Same as in panel I, but for the LFP signal.

To quantify the relative contributions from the various populations, we show in Figure 9D,E the CSD and LFP variances across time (corresponding to PSDs summed over all nonzero frequencies) for all depths (Lindén et al. 2011, Łe¸ski et al. 2013). For the compound CSD, only the L23E neurons contribute substantially in the superficial channels (channels 1–4) (Fig. 9D). At deeper contacts, the main contributing population is L6E. L4E and L5E populations make sizable contributions only in the channels closest to their somatic location reflecting that the net associated return currents of their distributed synaptic inputs are largely restricted to somatic regions (Lindén et al. 2010). The bulk of the variance of the CSD in layer 5 arises in about equal parts from the L6E and relatively sparse L5E populations. The magnitude of the depth-resolved CSD variances of each layer's inhibitory population (L23I, L4I, ...) is consistently one order of magnitude or more smaller than that of the corresponding excitatory cell populations. Moreover, they span comparatively small depth ranges as determined by the maximum extent of the dendrites.

The depth-resolved LFP variance (Fig. 9E) has similar features as the CSD variance, as expected from their common biophysical origin. However, volume conduction has some qualitative effects, such as the reduced relative contribution from the L6E population in layer 2/3. As correlated synaptic inputs are known to amplify and increase the spread of the LFP generated from a cortical population (Lindén et al. 2011, Łe¸ski et al. 2013), this may reflect that the synaptic inputs to the L23E population are more correlated than to the L6E population.

Overall, we conclude that for the spontaneous activity in our network, the L23E and L6E populations dominate the compound LFP and CSD with only smaller contributions from the other excitatory populations. The signal variances from the inhibitory populations are typically much smaller than the contributions from the excitatory populations, suggesting that they can be safely neglected, in line with recent findings of Mazzoni et al. (2015).

Although transmembrane currents in inhibitory neurons provide little of the observed CSD and LFP, the inhibitory synaptic inputs onto excitatory neurons provide a substantial contribution. In the present network, inhibitory synaptic currents have a 4-fold larger amplitude compared with most excitatory synapses (Table 5), inhibitory neurons have higher overall firing rates compared with excitatory neurons (Potjans and Diesmann 2014), and inhibition specifically targets soma-proximal sections (Markram et al. 2004). Since our LFP-generating model is linear (passive cable formalism and linear synapse model), we can decompose the compound signal into contributions from each synapse type. For example, selective removal of either inhibitory or excitatory synaptic currents in the CSD and LFP modeling (Fig. 9F,G) shows that the CSD and LFP signals (Fig. 9H) are dominated by inhibitory synaptic currents and their associated return currents (Fig. 9I,J) when the network operates in the spontaneous asynchronous irregular firing regime. Furthermore, at most depths the variance of the signal arising from inhibitory input exceeds the compound variance, implying that the inhibitory component is generally negatively correlated with the excitatory component (Fig. 9I,J). Visual inspection of Figure 9FH also reveals that the inhibitory dominance appears particularly strong at high frequencies, in accordance with the firing-rate PSDs in Figure 8B showing less power of inhibitory spiking, and thus inhibitory synaptic-input currents, at low frequencies.

The relative contribution from excitation and inhibition to the CSD and LFP depends on thalamic input and network state, however. For oscillatory thalamic input, the CSDs and LFPs from excitatory (Fig. 10A) or inhibitory synapses (Fig. 10B) alone show much stronger oscillations than the compound signals (Fig. 10C). This follows from the observations that the contributions to the CSD and LFP from excitatory and inhibitory synapses are anticorrelated, which, in turn, is a consequence of the dynamical balance between excitation and inhibition in asynchronous-irregular states of balanced random networks (Hertz 2010; Renart et al. 2010; Tetzlaff et al. 2012).

Figure 10.

Decomposition of CSD and LFP into contributions due to excitatory and inhibitory inputs for thalamic activation. (AE) Oscillatory thalamic activation (f = 15 Hz). (FJ) Transient thalamic activations at t=900+n·1000ms for n{0,1,2,3,4}. Same row-wise figure arrangement as in Figure 9F–J.

Figure 10.

Decomposition of CSD and LFP into contributions due to excitatory and inhibitory inputs for thalamic activation. (AE) Oscillatory thalamic activation (f = 15 Hz). (FJ) Transient thalamic activations at t=900+n·1000ms for n{0,1,2,3,4}. Same row-wise figure arrangement as in Figure 9F–J.

For transient thalamic activation, the same cancellation of excitation and inhibition can be observed (Fig. 10FH), although it is more pronounced at some depths (layers 5 and 6, i.e., channels 11–13) than at others, depending on whether excitation and inhibition are in phase or not (Fig. 10I,J). For such strong and transient thalamic input, the network activity is briefly imbalanced as inhibition cannot keep up with thalamic excitation on the short-time scales. Since thalamic input is most prominent in layer 4 (channels 7–10), fewer cancellation effects are present here: in channel 9 in Figure 10I the total CSD variance is, for example, seen to be larger than the individual contributions from inhibitory and excitatory synapses.

Effect of Input Correlations

Synaptic inputs to two neighboring cells are typically correlated because 1) they receive, to some extent, inputs from the same presynaptic sources (“shared-input correlation”), and 2) the spike trains of the presynaptic neurons may be correlated (“spike-train correlation”). The net synaptic-input correlation is determined by the interplay between these two contributions, shared-input correlations and spike-train correlations (Renart et al. 2010; Tetzlaff et al. 2012). As the LFP is largely generated by synaptic inputs, synaptic-input correlations result in correlated single-cell LFP contributions ϕi(r,t) (for details, see Lindén et al. 2011; Łe¸ski et al. 2013). As outlined in the following, these single-cell LFP correlations play a dominating role for the spectrum of the compound LFP.

The power spectrum Pϕ(r,f) (Eq. 36) of the compound LFP ϕ(r,t)=i=1Nϕi(r,t) of a population of N neurons is given by  

(15)
Pϕ(r,f)=NPϕ¯(r,f)+N(N1)Pϕ¯(r,f)κϕ¯(r,f).
Here Pϕ¯(r,f) is the average single-cell LFP power spectrum (Eq. 33) and κϕ¯(r,f) the average pairwise single-cell LFP coherence (Eq. 35), a measure for cross-correlations, across all cells. Note that, while the first term in Equation (15) scales linearly with the number of neurons N, the second term is proportional to N(N1)N2 for large N. Hence, for large N, even small cross-correlations may dominate the spectrum of the compound LFP. Here, we investigate this situation by calculating the power spectrum Pϕ0(r,f) of the compound LFP under the assumption of zero cross-correlation (where it simply reduces to a sum over single-cell spectra Pϕi(r,f)), and compare to the true spectrum Pϕ(r,f). The ratio between these quantities is given by  
(16)
Pϕ(r,f)Pϕ0(r,f)=1+(N1)κϕ¯(r,f).
With weak or no cross-correlations, that is, (N1)κϕ¯(r,f)1, the ratio approaches unity, and the power of the compound LFP is essentially the sum of the power of the single-cell LFPs. For Nκϕ¯(r,f)1, that is, in the correlation-dominated regime, this ratio is instead proportional to the number of neurons N. Note also that anti-correlated signals (κϕ¯(r,f)<0) may lead to a ratio Pϕ(r,f)/Pϕ0(r,f)<1.

In the example application, both for spontaneous (Fig. 11A,B) and for evoked activity (Fig. 11C,D) the compound power spectra Pϕ(r,f) are systematically (across channels and frequencies) larger than Pϕ0(r,f), demonstrating the importance of cross-correlations in the present network. Depending on the recording depth and frequency, the ratio varies from ∼ 1 to 103 (see Fig. 11B,D). For spontaneous activity (see Fig. 8AE), the largest effects of cross-correlations are typically found at higher frequencies (Fig. 11A,B). At low frequencies, cross-correlations are suppressed by inhibitory feedback (cf. Fig. 8B, Tetzlaff et al. (2012)). The thalamic sinusoidally modulated input to the network (Fig. 8FJ) synchronizes single-cell CSDs and LFPs at the stimulus frequency and gives a large boost of the LFP power at this frequency (see peak in Fig. 11C,D). Close inspection reveals that there is also in fact a boost of the power at ∼ 80 Hz, but much less so than at ∼ 15 Hz. Note that the external activation hardly affects the single-cell spectra (see red curves in Fig. 11A,C). LFP synchronization is thus mainly encoded in the phase of ϕi(r,t).

Figure 11.

Effect of single-cell LFP cross-correlations on compound-LFP power spectra during spontaneous activity (A,B) and for oscillatory thalamic input (C,D). A,C Compound-LFP power spectra Pϕ(r,f) (black traces) and compound spectra Pϕ0(r,f) (red traces) obtained when omitting cross-correlations between single-cell LFPs (red traces; computed for 10% of the cells and multiplied by a factor 10) at recording channels corresponding approximately to the centers of layers 1, 2/3, 4, 5, and 6. B,D Depth and frequency-resolved ratio Pϕ(r,f)/Pϕ0(r,f) of LFP power spectra, cf. Equation (16).

Figure 11.

Effect of single-cell LFP cross-correlations on compound-LFP power spectra during spontaneous activity (A,B) and for oscillatory thalamic input (C,D). A,C Compound-LFP power spectra Pϕ(r,f) (black traces) and compound spectra Pϕ0(r,f) (red traces) obtained when omitting cross-correlations between single-cell LFPs (red traces; computed for 10% of the cells and multiplied by a factor 10) at recording channels corresponding approximately to the centers of layers 1, 2/3, 4, 5, and 6. B,D Depth and frequency-resolved ratio Pϕ(r,f)/Pϕ0(r,f) of LFP power spectra, cf. Equation (16).

In conclusion, cross-correlations between single-cell LFP contributions play a pivotal role in shaping the compound LFP spectra (similar for CSD spectra, results not shown). To account for the dominant features of the LFP (CSD) in such models, it is therefore essential to include the main factors determining the synaptic-input correlations, that is, realistic correlations in presynaptic spike trains and shared-input structure. The findings presented here for the cortical microcircuit model hold in general in the presence of correlated activity. Only the details of the spectra depend on the specific underlying network dynamics.

Network Downscaling

Due to the computational cost associated with modeling LFPs, it would be desirable to downsize the postsynaptic populations of multicompartment model neurons to a fraction γN (γ(0,1)) while leaving the point-neuron network at full size (N) and at the same time preserve the in-degrees, that is, the number of inputs, of each postsynaptic cell. The power spectrum of the full-scale LFP can indeed be estimated from the population-averaged single-cell power spectra Pϕγ¯(r,f)Pϕ¯(r,f) and coherences κϕγ¯(r,f)κϕ¯(r,f) computed for downsized networks by means of Equation (15). These quantities are preserved except for deviations due to smaller sampling size γN (Eq. 36 in Table 10). However, due to lack of phase information in the power spectra, one cannot estimate the LFP time course.

One could attempt to obtain a time-course estimate ϕγξ(r,t), that is, a “low-density LFP prediction”, of the full-scale signal ϕ(r,t) by upscaling single-cell LFPs ϕi(r,t) computed in the downsized setup by a scalar factor ξ (cf. Eq. 31). Such a naive upscaling can grossly recover the amplitude of the full-scale LFP ϕ(r,t), but it still only partially reconstructs its detailed time course (Fig. 12A,E). Also, this approach does not generally give accurate power spectra as the two terms in Equation (15) scale differently with ξ: The rescaling introduces a prefactor ξ2 in the population-averaged single-cell power spectra Pϕγξ¯(r,f)ξ2Pϕ¯(r,f), while the coherences κϕγξ¯(r,f)κϕ¯(r,f) are unchanged. Thus, the compound spectra Pϕ(r,f) and Pϕγξ(r,t) of the full-size LFP and the low-density LFP predictor, respectively, differ. Their ratio  

(17)
Pϕ(r,f)Pϕγξ(r,f)=1+(N1)κϕ¯(r,f)γξ2+γξ2(γN1)κϕ¯(r,f)={1/(γξ2)forκϕ¯(r,f)=01/(γ2ξ2)forκϕ¯(r,f)=1
demonstrates that in the general case there is no scaling factor ξ which allows for the recovery of the full-size compound LFP power, that is, makes the ratio in Equation (17) equal to one for all spatial positions r and frequencies f. This can only be done in the special case where κϕ¯(r,f) is a constant c. Here the two extreme cases correspond to no correlation (κϕ¯(r,f)=0 with ξ=1/γ) and full correlation between all single-cell signals (κϕ¯(r,f)=1 with ξ=1/γ).

Figure 12.

Prediction of LFPs from downsized networks. Top row: Spontaneous activity. Bottom row: Oscillatory thalamic activation. (A,E) Full-scale LFP traces ϕ(r,t) (black) and low-density predictors ϕγξ(r,t) (red) obtained from a fraction γ=0.1 of neurons in all populations and upscaling by a factor ξ=γ1/2. (B,F) Correlation coefficients (gray bars) between full-scale LFP and low-density predictor shown in panels A and E, respectively. The dashed lines denote 1% significance levels obtained after computing the chance correlation coefficients for 1000 trials. (C,G) Power spectra Pϕ(r,f) and Pϕγξ(r,f) of full-scale LFPs (black) and low-density predictors with γ=0.1 and ξ=γ1/2 (red) or ξ=γ1 (gray). (D,H) Ratio Pϕ(r,f)/Pϕγξ(r,f) between power spectra of full-scale LFP and low-density predictor with γ=0.1 and ξ=γ1/2 (cf. Eq. 17).

Figure 12.

Prediction of LFPs from downsized networks. Top row: Spontaneous activity. Bottom row: Oscillatory thalamic activation. (A,E) Full-scale LFP traces ϕ(r,t) (black) and low-density predictors ϕγξ(r,t) (red) obtained from a fraction γ=0.1 of neurons in all populations and upscaling by a factor ξ=γ1/2. (B,F) Correlation coefficients (gray bars) between full-scale LFP and low-density predictor shown in panels A and E, respectively. The dashed lines denote 1% significance levels obtained after computing the chance correlation coefficients for 1000 trials. (C,G) Power spectra Pϕ(r,f) and Pϕγξ(r,f) of full-scale LFPs (black) and low-density predictors with γ=0.1 and ξ=γ1/2 (red) or ξ=γ1 (gray). (D,H) Ratio Pϕ(r,f)/Pϕγξ(r,f) between power spectra of full-scale LFP and low-density predictor with γ=0.1 and ξ=γ1/2 (cf. Eq. 17).

The substantial scaling effects observed for our microcircuit model in the asynchronous state (Fig. 12AD) suggest that correlations cannot be neglected even when modeling the LFP for spontaneous network activity. Choosing the scaling factor ξ=1/γ corresponding to κϕ¯(r,f)=0 (red lines in Fig. 12A,C) leads to a severe underestimation of the full-size compound power spectrum (Fig. 12C,D). Even though the correlation (i.e., Pearson's correlation coefficients) between the full-size LFP signals and low-density LFP predictions are quite high (Fig. 12B), the power ratios (Fig. 12D) reveal that the rescaled signals are systematically wrong in frequency bands where single-cell LFPs are most strongly correlated (i.e., the frequencies for which the compound spectra are much larger than the predictions when omitting cross-correlations (Fig. 11A,B)). Assuming the full-correlation scaling factor ξ=1/γ, however, typically overestimates the full-size compound power spectrum (cf. gray spectra in Fig. 12C), particularly at low frequencies.

The results for the sinusoidally stimulated network (Fig. 12EH) are quite similar to the spontaneous-activity results, except around the stimulation frequency 15 Hz where the modulated input leads to strongly correlated single-cell LFP contributions and a strong boost of the compound LFP. The approximate downscaling procedure assuming the full-correlation scaling factor ξ=1/γ thus essentially agrees with the full-size compound spectrum for 15 Hz (while giving a strong overestimation for frequencies other than ∼ 15 and ∼ 80 Hz).

We note in passing that in contrast to the power spectra, the computation of the spike-trigger-averaged LFP (stLFP) (Swadlow et al. 2002; Nauhaus et al. 2009; Jin et al. 2011; Denker et al. 2011) in downsized networks do not have a principled problem due to cross-correlations between single-cell LFPs. As stLFPs are linearly dependent on the single-neuron LFP contributions, the only principled problem with downsizing is increased noise in the estimates due to sampling over fewer postsynaptic neurons.

LFP Prediction from Population Firing Rates

An important question in systems neuroscience is to what extent the dynamics of networks of thousands or millions of neurons can be described by much simpler mathematical descriptions in terms of neural populations (Deco et al. 2008; Blomquist et al. 2009). Likewise, we here ask the question of whether LFPs can be predicted from knowledge of the population firing rate (Einevoll et al. 2007; Moran 2008; Einevoll et al. 2013). The hybrid scheme is excellently suited for testing and development of simplified numerical schemes for LFP prediction as the ground truth, that is, the LFP from the full network, is available as benchmarking data.

The use of current-based synapses and passive dendrites in the present application of the hybrid scheme renders synaptic events independent of each other in the LFP prediction. This inherent linearity results in a unique spatio-temporal relation HXi(r,τ) for τ[,] between a spike event of a point neuron i in population X and its contribution to the compound LFP ϕ(r,t) from all its postsynaptic multicompartment model neurons. In this scheme, the link is causal, that is, the spikes drive the LFP, so that HXi(r,τ)=0 for τ<0 (as in LPA (Einevoll et al. 2007)). HXi(r,τ) encompasses connectivity, spike transmission delays and all postsynaptic responses including effects of synaptic-input currents and passive return currents. With a linear, current-based model such as our example cortical column, it is in principle possible by linear superposition to fully reconstruct the compound LFP if HXi(r,τ) and the spike times til are known for all neurons i in each population of the network.

It is, however, in the case of large networks impractical to assess each HXi(r,τ), as the LFP response needs to be determined for every neuron separately. In contrast, a large reduction in dimensionality can be achieved by determining the population-averaged LFP responses H¯X(r,τ) of a spike within each population X. We thereby ignore heterogeneity in kernels HXi(r,τ) due to the variability in the connections from neurons in population X. An approximate compound LFP ϕ*(r,t) based on population firing rates (Einevoll et al. 2007) can be computed from these extracted population kernels by means of the convolution ϕ*(r,t)=X(νX*H¯X)(r,t), where νX(t) are the instantaneous population firing rates.

Here we estimate the population LFP kernels by computing the response to synchronous activation of all neurons in a population (Fig. 13). The spatio-temporal kernels H¯X(r,τ) are extracted from time slices [tX20ms,tX+20ms] of the compound LFP response ϕ(r,t)/NX, where NX is the number of neurons in a presynaptic population. The procedure results in unique kernels H¯X(r,τ) for each excitatory and inhibitory population in the network (Fig. 13A).

Figure 13.

Linear prediction of LFPs from population firing rates. (A) LFP responses H¯X(r,τ) (kernels) to simultaneous firing of all neurons in a single presynaptic population X (see subpanel titles) at time τ = 0 ms, normalized by size NX of the presynaptic population (red/blue: responses to firing of excitatory/inhibitory presynaptic populations). (B) Spike-trigger-averaged LFPs (stLFP) triggered on spikes of L5E neurons during spontaneous activity (left) and oscillatory thalamic network activation (right), averaged across all L5E spikes (T = 5000 ms simulation time). (C,F) LFP traces of the full model (black) compared with predictions (red) obtained from superposition of linear convolutions of population firing rates νX with LFP kernels H¯X(r,τ) shown in panel A. (D,G) Correlation coefficients (gray bars) between LFPs and population-rate predictors shown in panels C and F. The dashed lines denote 1% significance levels obtained after computing the chance correlation coefficients for 1000 trials. (E,H) Power spectra of LFPs (black) and the population-rate predictors (red) for different recording channels. Panels C–E and F–H show results for spontaneous activity and oscillatory thalamic activation, respectively.

Figure 13.

Linear prediction of LFPs from population firing rates. (A) LFP responses H¯X(r,τ) (kernels) to simultaneous firing of all neurons in a single presynaptic population X (see subpanel titles) at time τ = 0 ms, normalized by size NX of the presynaptic population (red/blue: responses to firing of excitatory/inhibitory presynaptic populations). (B) Spike-trigger-averaged LFPs (stLFP) triggered on spikes of L5E neurons during spontaneous activity (left) and oscillatory thalamic network activation (right), averaged across all L5E spikes (T = 5000 ms simulation time). (C,F) LFP traces of the full model (black) compared with predictions (red) obtained from superposition of linear convolutions of population firing rates νX with LFP kernels H¯X(r,τ) shown in panel A. (D,G) Correlation coefficients (gray bars) between LFPs and population-rate predictors shown in panels C and F. The dashed lines denote 1% significance levels obtained after computing the chance correlation coefficients for 1000 trials. (E,H) Power spectra of LFPs (black) and the population-rate predictors (red) for different recording channels. Panels C–E and F–H show results for spontaneous activity and oscillatory thalamic activation, respectively.

In the example application, the population kernels H¯X(r,τ) differ significantly between populations. Excitatory spike events result in prominent LFP negativities at depths where most connections are made, such as in layer 4 (channels 7–10) for TC connections (H¯TC(r,τ), column 1 in Fig. 13A, cf. Fig. 5C). In contrast, spikes of inhibitory point neurons on average produce prominent LFP positivities in their corresponding layer, such as in layer 2/3 (channels 3–6) for population L23I (column 2 in Fig. 13A). In all cases, the signatures of opposite-sign return currents and also other, weakly connected populations are seen across depth.

As seen in Figure 13C,F the population-rate predictions are in good qualitative agreement with the ground-truth LFP for both spontaneous and sinusoidally modulated network activity with correlation coefficients (cc) between 0.48 and 0.94 for spontaneous activity (Fig. 13D) and 0.52 and 0.98 for thalamically evoked oscillations (Fig. 13G). Overall, the population-rate predictions appear to be best for the lower frequencies (Fig. 13E,H), while the inherent variability in the individual point-neuron kernels HXi(r,τ) (which is not accounted for in the population approximation) has a larger effect on the higher frequencies. This can be understood on biophysical grounds, as the lower frequencies are expected to mainly reflect the gross anatomical features of the postsynaptic populations and their presynaptic connections patterns where the individual variability plays a lesser role (Lindén et al. 2010; Pettersen et al. 2012). The correlation coefficients and power spectra of the population-rate prediction thus show that the population-rate LFP predictor is more accurate than the low-density LFP predictors (Fig. 12C,F) in case of substantial downscaling.

Although the spike-trigger-averaged LFP (stLFP) (Swadlow et al. 2002; Nauhaus et al. 2009; Jin et al. 2011; Denker et al. 2011), calculated as the cross-covariance between the population spike rate νX(t) and the compound LFP ϕ(r,t) divided by the total number of spikes (i.e., 0TνX(t)dt), is related to our LFP population kernels, it measures very different aspects of cortical dynamics. The population kernels H¯X(r,τ) are causal and independent of effects of spike-train correlations. The stLFP, however, is non-causal and strongly depends on spike-train correlations, and thus also network state (Einevoll et al. 2013). The stLFP is thus not only very different from H¯X(r,τ), it also varies strongly between the spontaneous and sinusoidally modulated network state (see example for L5E neurons in Fig. 13B).

Discussion

We have here described a hybrid modeling scheme for computing the LFP incorporating both large-scale neural-network dynamics and the biophysics underlying LFP generation on the single-neuron level. The hybrid modeling scheme was illustrated with a full-scale network model of a cortical column in early sensory cortex (Potjans and Diesmann 2014), and the impact of individual populations, network dynamics, and cell density on the mesoscopic LFP signal was investigated.

The Hybrid LFP Modeling Scheme

The hybrid scheme combines the simplicity and efficiency of point-neuron network models with the biophysics-based modeling of LFP by means of multicompartment model neurons with detailed dendritic morphologies. The neuronal network dynamics are governed by the point-neuron network model independent of LFP predictions. The spikes of the point-neuron network are distributed to the synapses of the multicompartment model neurons with realistic cell-type and layer-specific connectivity. Synapse activation results in spatially distributed transmembrane currents, which are mapped to an LFP signal according to well-established volume-conduction theory.

A main motivation for developing the hybrid LFP modeling scheme was to obtain the ability to compute LFPs for a key class of network models that are amenable to mathematical analysis and can provide intuitive understanding of emerging network dynamics, namely point-neuron models. Similar to networks of anatomically and biophysically detailed neuron models, point-neuron networks can generate realistic spiking activity. In addition, the hybrid modeling scheme brings a substantial computational advantage: with present-day computing and software technologies, point-neuron networks with ∼ 100, 000 neurons can be modeled with laptop computers, and networks comprising millions of point neurons can be routinely simulated on high-performance compute facilities (Helias et al. 2012; Kunkel et al. 2014). Until now, the largest simulation of LFPs based on networks of multicompartmental neuron models with reconstructed morphologies, in contrast, comprised 12, 000 neurons and was done on a Blue Gene/P supercomputer with 4096 CPUs (Reimann et al. 2013). The linearity of electromagnetic theory allows for the implementation of the LFP hybrid modeling scheme as an “embarrassingly” parallel operation (Foster 1995). Therefore, the results for the cortical microcircuit application with ∼ 78,000 neurons were obtained with only 256 CPUs, and could even be acquired with much smaller computing architectures.

A full implementation of the hybrid scheme is provided by the freely available

Python
module
hybridLFPy
(http://github.com/INM-6/hybridLFPy). Our model implementation in
hybridLFPy
relies on the publicly available
NEURON
software (Carnevale and Hines 2006) as a simulation backend
Python
with
LFPy
(Lindén et al. 2014) for calculating single-cell LFP contributions. This ensures flexibility and compatibility with a large library of existing neuron models, with or without active channels and with morphologies of arbitrary levels of detail, obtained from ModelDB (Hines et al. 2004), NeuroMorpho.org (Ascoli et al. 2007), or other resources. While we did use
NEST
(Eppler et al. 2015) for simulating our reference network, the
hybridLFPy
module can be used in combination with any other neural-network simulation software.

The present hybrid LFP scheme involves several assumptions with respect to 1) the generation of realistic spiking activity, 2) forward modeling of extracellular potentials, and 3) the combined use of point-neuron networks and multicompartment modeling. In the following, we review the main assumptions and discuss potential extensions:

1. Spike-train generation by point-neuron networks: Although highly simplified, single-compartment models of individual neurons (point-neuron models) can mimic realistic spiking for a variety of cell types (Izhikevich 2003; Kobayashi et al. 2009; Yamauchi et al. 2011) and can make accurate predictions of single-cell firing responses under in-vivo like conditions (Jolivet et al. 2008; Gerstner and Naud 2009). Moreover, networks of point neurons can reproduce a number of activity features observed in vivo, such as spike-train irregularity (Softky and Koch 1993; van Vreeswijk and Sompolinsky 1996; Amit and Brunel 1997; Shadlen and Newsome 1998), membrane-potential fluctuations (Destexhe and Paré 1999), asynchronous firing (Ecker et al. 2010; Renart et al. 2010; Ostojic 2014), correlations in neural activity (Gentet et al. 2010; Okun and Lampl 2008; Helias et al. 2013), self-sustained activity (Ohbayashi et al. 2003; Kriener et al. 2014), and realistic firing rates across laminar cortical populations (Potjans and Diesmann 2014). Recently Rössert et al. (2016) demonstrated that the multicompartment-neuron network of Markram et al. (2015) could be reduced to a point-neuron network using an automated step-by-step workflow without significantly affecting the spiking statistics in the network. Nevertheless, under circumstances where simple point neurons may fail to produce plausible spike patterns for a particular input, for example, in the presence of highly nonlinear shaping of postsynaptic responses or dendritic computation (reviewed, e.g., by London and Häusser 2005), more sophisticated point neurons (Breuer et al. 2014) may be required. However, the present hybrid scheme makes no assumptions on generators of spiking activity. In principle, network spikes could, as an example, be replaced by statistical models of spike generation (Lindén et al. 2011; Łe¸ski et al. 2013), multicompartment neuron spikes or even experimentally measured spiking activity.

2. Biophysical forward modeling of LFPs: The biophysical forward model described by Equation (10), implemented in

LFPy
(Lindén et al. 2014), underlies the presently used computational scheme for LFPs of point-neuron networks. This forward model is based on well-established volume conductor theory (Rall and Shepherd 1968; Holt and Koch 1999) and assumes an infinite, isotropic (same in all directions), homogeneous (same in all positions), and ohmic (frequency-independent) extracellular medium represented by a scalar conductivity σe. However, one could generalize the forward model in a straightforward manner to account for anisotropy (Nicholson and Freeman 1975; Logothetis et al. 2007; Goto et al. 2010), or jumps in conductivities at tissue interfaces, for example, using the so-called method of images (Pettersen et al. 2006; Gold et al. 2006; Hagen et al. 2015; Ness et al. 2015). For even more complicated geometrical spatial variations of the conductivity, the forward modeling problem can always be solved by means of finite element modeling (Ness et al. 2015). Recent experiments have only found a small frequency dependence of the extracellular conductivity σe at LFP frequencies (f500Hz, Logothetis et al. 2007; Wagner et al. 2014), but see Gabriel et al. (1996; 2009), Bédard and Destexhe (2009). In any case the forward model could still be applied with frequency-dependent conductivity by means of Fourier decomposition where each frequency component of the LFP signal is considered separately. For more information on possible generalizations of the biophysical forward-modeling scheme, see Pettersen et al. (2012). Finally, we assumed the so-called disc-electrode approximation and averaged the computed LFP signal across the electrode surface (Moulin et al. 2008; Lindén et al. 2014; Ness et al. 2015). Although electrode impedance will affect the measurement, it appears that confounding effects from this can easily be avoided with present-day LFP recording techniques (Nelson and Pouget 2010), hence few compelling reasons exist to incorporate additional temporal filters.

3. Combined use of point-neuron and multicompartmental models: The key approximation in the hybrid LFP scheme comes from the combined use of point-neuron (single-compartment) and multicompartmental neuron models. The multicompartment neurons are mutually unconnected, have no outgoing (efferent) connections, and are solely used to compute the LFP. Furthermore, due to dendritic filtering, the somatic postsynaptic potentials in the multicompartment model neurons are not identical to those of their point-neuron counterparts. This inconsistency could, at least partially, be resolved by adjusting the amplitudes and temporal shapes of the synaptic currents in either the multicompartment neurons or the point neurons (Koch and Poggio 1985; Wybo et al. 2013; 2015).

Applications of the Hybrid LFP Scheme

The hybrid scheme is not limited to the example point-neuron network model and the particular multicompartment neuron models chosen here. It can be applied to networks 1) of arbitrary topology (graph structure, distance dependencies, and dimensions), 2) with any number of populations, 3) with arbitrarily complex point-neuron (e.g., LIF, Izhikevich, MAT, and Hodgkin-Huxley) and synapse dynamics (e.g., current-based, conductance-based, static, and plastic), and 4) any level of biophysical detail in multicompartment neuron models (e.g., morphologies and active channels).

For illustration, we used the hybrid scheme to compute LFPs along a virtual laminar multi-electrode from activity in a multilayered spiking point-neuron network, modeling signal processing in a patch of primary visual cortex. The network consisted of ∼ 78, 000 neurons organized in 4 layers, each with an excitatory and an inhibitory population, representing a cortical patch of ∼ 1 mm2 (Potjans and Diesmann 2014). Altogether 16 different cell types and 10 different morphologically reconstructed neurons were used in the LFP calculation. This cortical microcircuit model is well suited for the illustration of the hybrid scheme due to its 1) minimum level of detail in single-neuron dynamics of both point neurons (LIF) and multicompartment neurons (passive membranes), 2) realistic neuron density allowing investigation of the effects of correlations and scaling of network size, and 3) its spatial organization of multiple populations across cortical layers which yields cancellation effects not captured by LFP proxies such as in Mazzoni et al. (2015).

Even though the example application was based on a generic network model biased toward cat visual cortex and not tuned to address specific experiments, its spiking activity nevertheless matched experimental findings on distributions of firing rates, asynchronous irregular activity and propagation of evoked activity (Potjans and Diesmann 2014). We even observed the predicted LFP to be in qualitative accordance with LFP measurements in primary sensory cortices from a variety of animal species and sensory modalities in terms of 1) LFP amplitude, for both spontaneous (0.11mV, see Maier et al. 2010; Hagen et al. 2015; Reyes-Puerta et al. 2016), electrically stimulated activity (13mV, see Mitzdorf 1985; Castro-Alamancos and Connors 1996) and 2) stimulus-evoked spatiotemporal LFP and CSD patterns (Di et al. 1990; Einevoll et al. 2007; Reyes-Puerta et al. 2015). The fact that our spontaneous and evoked LFP amplitudes are comparable to the observed amplitudes supports the overall biological plausibility of the hybrid LFP scheme.

For the present example, the LFP was dominated by synaptic inputs, and their associated return currents, on excitatory neurons, in particular onto pyramidal cells in layers 2/3 and 6. Furthermore, contributions from inhibitory synaptic inputs typically dominated the contributions from the excitatory inputs, particularly for LFPs stemming from spontaneous network activity. Although the main point of employing the present example was to illustrate the use of the hybrid LFP scheme and not to make predictions for specific neural systems, we note in passing that a dominance of inhibitory synaptic inputs appears to be in agreement with LFPs generated in the CA3 region of hippocampus as observed in an in vitro setting (Bazelot et al. 2010) and in cortex of both humans and macaque in vivo (Telenczuk et al. 2016). In accordance with a previous study (Lindén et al. 2011), we found that correlations in synaptic input play a major role in determining the CSD and LFP stemming from the network activity, for both spontaneous and stimulus-evoked activity. We further showed that due to inevitable correlations between synaptic-input currents the main features of the LFP can only be correctly predicted by a full-scale model. As our ambition is to compute LFPs also for extended point-neuron networks with millions of neurons covering, for example, entire cortical areas, we finally demonstrated how the hybrid modeling scheme can predict LFPs from population firing rates rather than from spikes of individual neurons (Einevoll et al. 2007).

The present microcircuit model application involving simplified, passive multicompartment populations was here used to study the effect of the spatial connectivity on the laminar pattern of spontaneous and stimulus-evoked CSD and LFP signals. The application thus represented a minimal approach incorporating spatial features in LFP predictions of multilaminar point-neuron networks. However, several of the simplifying model assumptions made in the present example application can straightforwardly be generalized. In particular, such generalizations concern 1) the synaptic connectivity between point neurons and the equivalent multicompartment neurons, 2) the absence of active conductances, 3) the positioning of the cells, 4) the reconstructed morphologies, 5) the representation of external inputs, and 6) the fact that the model only encompassed the local circuitry of a ∼ 1 mm2 patch of cortex.

(1) Synaptic connectivity: While the population-specific connection probabilities, delay distributions, synapse time constant and mean synaptic weights were identical for the connections in the point-neuron network and those between point neurons and LFP-generating multicompartment neurons, the exact realizations of the two types of connectivities were different. In contrast to the point-neuron network, each cell of a particular type had a fixed in-degree, that is, a fixed number of synaptic inputs, a fixed synaptic current amplitude and a layer-specific synapse positioning that depended on compartment surface area, in the LFP modeling step. The use of the hybrid LFP scheme, however, is not restricted to these or any other specific assumptions about the synaptic connectivity patterns. One could, for example, gather all point-neuron network connections and corresponding weights and delays for use with the LFP-generating multicompartment neuron populations. However, for large networks this would require additional computing and memory resources as the number of recorded connection weights and delays grow proportionally to N2 in a network of N neurons with fixed connection probabilities. The present layer-specific synapse positioning could also be specialized further, for example, to account for dendrite-specific connections that may be relevant for LFP predictions in particular if nonlinear synapse and dendrite models are used.

(2) Active conductances: In the present study, we neither included the active channels underlying spike generation, nor active dendritic conductances in the multicompartment neuron models (Remme and Rinzel 2011). Experiments suggest that the contribution to the LFP from the former is small in stimulus-evoked recordings from sensory cortex, at least for the low frequencies of the LFP (Pettersen et al. 2008), but see Ray and Maunsell (2011). In any case the contribution of the spikes to the extracellular potentials, including the LFP, could be included in the present scheme by giving each spike produced in the point-neuron network simulations a cell-type-specific spatiotemporal signature in the computation of the extracellular potential (e.g., as calculated in Holt and Koch 1999; Hagen et al. 2015). Substantial effects of active dendritic conductances on the LFP were observed in a two-layered model by Reimann et al. (2013), but this should be further explored. Recent modeling results (Ness et al. 2016) suggest that for the purposes of LFP prediction, active dendritic conductances can, at least for subthreshold potentials, be effectively described by means of “quasi-active linearized” theory (Sabah and Leibovic 1969; Koch 1984; Remme and Rinzel 2011; Ness et al. 2016). This simplifies the LFP modeling substantially and makes the computation similar in complexity to the present case with purely passive membranes.

(3) Soma positioning: For model conciseness, the somas of all neurons belonging to the same specific cortical network populations were set to have the same range across cortical depth, see Figure 6E. By instead assuming a biologically more plausible distribution of soma depths, the CSD profiles are expected to be spatially smoothed compared with the profiles seen for the present distributions, for example, in Figure 6F. Also the LFP profiles will be affected, but to a lesser degree since the LFP profiles already are spatially smoothed due to volume conduction effects.

(4) Reconstructed morphologies: We further chose to rely on a small number of highly detailed reconstructed dendritic morphologies from experimental preparations, and partly reused morphologies across neural populations. Obviously, larger sets of distinct reconstructed dendritic morphologies can be used in future studies as they become available. The effect of dendritic morphologies on generated LFP can be further assessed by the use of stylized (Tomsett et al. 2014; Głąbska et al. 2014) or artificially grown morphologies resembling real neurons (Cuntz et al. 2010; 2011; Torben-Nielsen and De Schutter 2014; Mazzoni et al. 2015). While our morphologies were used as is, effects on the predicted LFP from varying morphologies could be quantified using our present measures of LFP power spectra and variance across depth. Such measures could also be extended to include effects of lateral displacement of the recording electrode (Lindén et al. 2011; Einevoll et al. 2013; Łe¸ski et al. 2013; Tomsett et al. 2014), but also cross-channel coherences and stLFPs should be considered as well.

(5) External input: In our point-neuron network modified from Potjans and Diesmann (2014), we assumed the depolarizing input from surrounding cortex and remote areas to be represented by deterministic, fixed-amplitude input currents (DC inputs) rather than independent Poisson spike trains. We thus avoided the generation and storage of O(105) high-rate uncorrelated Poisson spike trains and distributing the corresponding spike events onto the multicompartment model neurons, but this can be introduced in future applications.

(6) Scale: Our network model represents only an isolated cortical column under 1 mm2 of pial surface (Potjans and Diesmann 2014), but work is underway to extend the network to a larger scale, for example, by incorporating additional cortical areas (Schmidt et al. 2016) and extending the network size in the lateral directions (Senk et al. 2015). In addition to allowing for predictions of LFPs at several lateral positions as measured by multishank electrodes, the anticipated outcome is also an altered network dynamics and consequently altered LFPs. The present columnar model lacks structured input from other parts of cortex which corresponds to ∼ 50% of all excitatory synapses (see Potjans and Diesmann (2014) and references therein). For example, accounting for different cortical areas and their interactions would allow for the detailed investigation of pathway-specific LFP contributions as recently reviewed in hippocampus by Herreras et al. (2015).

Effects of Correlations and Network Size

Synaptic inputs to neurons are typically correlated, and as shown in this article and in earlier studies (Lindén et al. 2011; Łe¸ski et al. 2013), these correlations have a major impact on the properties of the generated LFP (and corresponding CSD). There are different contributions to these input correlations, shared presynaptic neurons and correlations in the presynaptic spiking activity, see, for example, Renart et al. (2010), Tetzlaff et al. (2012), and correct LFP predictions require that both effects are properly taken into account. The generation of presynaptic spike trains with realistic correlation structure requires networks of realistic size (van Albada et al. 2015). The contribution from shared presynaptic input can only be accounted for by using realistic statistics of inputs to the LFP-generating multicompartment model neurons, that is, realistic synaptic connection probabilities resulting in realistic statistics of shared inputs. If these two requirements are fulfilled, the properties of single-cell LFPs as well as the correlations between pairs of single-cell LFPs, will be correctly accounted for.

Given synaptic inputs with realistic statistics from sufficiently large point-neuron networks, do we actually need to represent the full population of LFP generating neurons in order to predict a realistic compound LFP? Or can we alternatively get a good estimate of the compound population-LFP signal from a downscaled population of neuronal LFP generators if we know the correct single cell and pairwise LFP statistics, thereby reducing computational costs? As shown in this article, the answer is negative: in the presence of (even tiny) synaptic-input correlations, a realistic compound LFP can only be generated by multicompartmental neuron populations with realistic size and cell densities.

The hybrid modeling scheme allows one to account for both, that is, realistic sizes of networks to generate spike trains with correct correlation structure and realistic sizes and cell densities of the LFP-generating multicompartmental neuron populations. This can be achieved since point-neuron networks can be simulated very efficiently, and the multicompartmental neurons are independent and can be simulated serially (or in an embarrassingly parallel manner).

Outlook

While we here have focused on the computation of the LFP based on output from spiking point-neuron networks, a similar hybrid approach could be used when the network dynamics is rather modeled in terms of firing rates or even neural fields (Deco et al. 2008). For our example case of a 4-layered cortical network with an excitatory and an inhibitory population in each layer, the present scheme could be adapted directly by replacing the set of spike trains for each population with the corresponding population firing rates (Schuecker et al. 2015) in the LFP-generating step. However, the feasibility and prediction accuracy of such a scheme would have to be investigated in detail.

Another natural development would be to consider other measurement modalities. The present LFP scheme already incorporates the prediction of ECoG (electrocorticography) signals, that is, the electrical signals recorded at the cortical surface, although the LFP forward-modeling scheme may have to be adjusted to account for the discontinuity in electrical conductivity at the cortical surface (Nunez and Srinivasan 2006). One straightforward method to account for discontinuities in our predictions is the so-called method of images, as previously used with our present forward modeling scheme (Hagen et al. 2015; Ness et al. 2015). An extension to EEG (electroencephalography) and MEG (magnetoencephalography) would in principle also be straightforward as the key variable linking single-neuron activity and the measured signal is the single-neuron current dipole moment (Hämäläinen et al. 1993; Nunez and Srinivasan 2006). This dipole moment can be computed from multicompartmental neuron models when the transmembrane or axial currents are known (Lindén et al. 2010; Ahlfors and Wreh 2015). Given the magnitude and orientation of the current dipole moments for all contributing neurons, the EEG and MEG signal can be computed by a linear superposition of single-cell contributions given an appropriate extracellular volume conductor model for the EEG signal. EEG dipole-source localization techniques and EEG forward models have sometimes assumed simplified spherical head models composed of N concentric shells, each with different conductivities (Nunez and Srinivasan 2006). Such models could be used also in the present hybrid scheme as well as realistic head models (Kaiboriboon et al. 2012). Note that MEG signal prediction from the current dipole moment is simplified due to its independence on conductivity, and that the variation in magnetic permeability is insignificant across soft tissue, bone, and air (Hämäläinen et al. 1993). Another measurement that could be modeled is voltage-sensitive dye imaging (VSDi) where the signal largely reflects average membrane potentials of dendrites close to the cortical surface. The spatial profile of the weights in the averaging procedure of the VSDi forward model will be determined by the spatial distribution of dye and the propagation of light in the neural tissue (Chemla and Chavane 2010a, 2010b; Tian et al. 2011).

Even though the LFP has been measured for more than half a century, the interpretation of the recorded data has so far largely been qualitative (Einevoll et al. 2013). The biophysical origin of the signal on the single-cell level appears well understood (Rall and Shepherd 1968; Holt and Koch 1999; Gold et al. 2006; Linden 2010; Buzsáki et al. 2012; Pettersen 2012; Einevoll et al. 2013), and several modeling studies have explored the link between neuron and network activity (Pettersen et al. 2008; Lindén et al. 2010; 2011; Łe¸ski et al. 2013; Reimann et al. 2013; Tomsett et al. 2014; Głąbska et al. 2014). However, the computation of LFPs from network activity has until now been too cumbersome and computer intensive to allow for practical exploration of the links between different types of network dynamics and the resulting LFP. Thus, a validation of network models against measured LFP data has essentially been absent. With the present hybrid LFP scheme, accompanied by the release of the simulation tool

hybridLFPy
, we believe that a significant step has been taken toward the goal of making combined modeling and measurement of the LFP signal a practical research tool for probing neural circuit activity.

Funding

European Union Seventh Framework Programme (FP7/2007–2013); (604102) Human Brain Project, HBP; (269912) BrainScaleS, the Helmholtz Association through the Helmholtz Portfolio Theme “Supercomputing and Modeling for the Human Brain” (SMHB), Juelich Aachen Research Alliance (JARA), the German Research Foundation (DFG; grant DI 1721/3-1 [KFO219-TP9]), the Danish Council for Independent Research and FP7 Marie Curie Actions COFUND (grant id: DFF-1330–00226); the Research Council of Norway (NFR, through ISP, NOTUR -NN4661 K).

Notes

We thank Zoltán F. Kisvárday and Armen Stepanyants for useful discussions in the early phase of this study and providing us with experimentally reconstructed morphologies of cat visual cortex neurons. Conflict of Interest: None declared.

References

Ahlfors
SP
,
Wreh
C
II
2015
.
Modeling the effect of dendritic input location on MEG and EEG source dipoles
.
Med Biol Eng Comput
 .
53
:
879
887
.
Amit
DJ
,
Brunel
N
1997
.
Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex
.
Cereb Cortex
 .
7
:
237
252
.
Anastassiou
CA
,
Perin
R
,
Markram
H
,
Koch
C
2011
.
Ephaptic coupling of cortical neurons
.
Nat Neurosci
 .
14
:
217
223
.
Armstrong-James
M
,
Fox
K
,
Das-Gupta
A
1992
.
Flow of excitation within rat barrel cortex on striking a single vibrissa
.
J Neurophysiol
 .
68
:
1345
1358
.
Ascoli
GA
,
Donohue
DE
,
Halavi
M
2007
.
Neuromorpho.org: a central resource for neuronal morphologies
.
J Neurosci
 .
27
:
9247
9251
.
Bastos
AM
,
Vezoli
J
,
Bosman
CA
,
Schoffelen
J-M
,
Oostenveld
R
,
Dowdall
JR
,
De Weerd
P
,
Kennedy
H
,
Fries
P
2015
.
Visual areas exert feedforward and feedback influences through distinct frequency channels
.
Neuron
 .
85
:
390
401
.
Bazelot
M
,
Dinocourt
C
,
Cohen
I
,
Miles
R
2010
.
Unitary inhibitory field potentials in the CA3 region of rat hippocampus
.
J Physiol
 .
588
:
2077
2090
.
Bazhenov
M
,
Stopfer
M
,
Rabinovich
M
,
Huerta
R
,
Abarbanel
HD
,
Sejnowski
TJ
,
Laurent
G
2001
.
Model of transient oscillatory synchronization in the locust antennal lobe
.
Neuron
 .
30
:
553
567
.
Bédard
C
,
Destexhe
A
2009
.
Macroscopic models of local field potentials and the apparent 1/f noise in brain activity
.
Biophys J
 .
96
:
2589
2603
.
Binzegger
T
,
Douglas
RJ
,
Martin
KAC
2004
.
A quantitative map of the circuit of cat primary visual cortex
.
J Neurosci
 .
24
:
8441
8453
.
Blomquist
P
,
Devor
A
,
Indahl
UG
,
Ulbert
I
,
Einevoll
GT
,
Dale
AM
2009
.
Estimation of thalamocortical and intracortical network models from joint thalamic single-electrode and cortical laminar-electrode recordings in the rat barrel system
.
PLoS Comput Biol
 .
5
:
e1000328
.
Bos
H
,
Diesmann
M
,
Helias
M
2016
. Identifying anatomical origins of coexisting oscillations in the cortical microcircuit. PLoS Comput Biol, accepted.
Brette
R
,
Rudolph
M
,
Carnevale
T
,
Hines
M
,
Beeman
D
,
Bower
JM
,
Diesmann
M
,
Morrison
A
,
Goodman
PH
,
Harris
FC
Jr
, et al
.
2007
.
Simulation of networks of spiking neurons: a review of tools and strategies
.
J Comput Neurosci
 .
23
:
349
398
.
Breuer
D
,
Timme
M
,
Memmesheimer
R-M
2014
.
Statistical physics of neural systems with nonadditive dendritic coupling
.
Physical Review X
 .
4
:
011053
.
Brunel
N
2000
.
Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons
.
J Comput Neurosci
 .
8
:
183
208
.
Buzsáki
G
,
Anastassiou
CA
,
Koch
C
2012
.
The origin of extracellular fields and currents - EEG, ECoG, LFP and spikes
.
Nat Rev Neurosci
 .
13
:
407
420
.
Buzsáki
G
,
Mizuseki
K
2014
.
The log-dynamic brain: how skewed distributions affect network operations
.
Nat Rev Neurosci
 .
15
:
264
278
.
Camuñas Mesa
LA
,
Quiroga
RQ
2013
.
A detailed and fast model of extracellular recordings
.
Neural Comput
 .
25
:
1191
1212
.
Carnevale
NT
,
Hines
ML
2006
.
The NEURON book
 .
Cambridge, UK
:
Cambridge University Press
.
Castro-Alamancos
MA
,
Connors
BW
1996
.
Spatiotemporal properties of short-term plasticity sensorimotor thalamocortical pathways of the rat
.
J Neurosci
 .
16
:
2767
2779
.
Chemla
S
,
Chavane
F
2010
a.
A biophysical cortical column model to study the multi-component origin of the VSDI signal
.
Neuroimage
 .
53
:
420
438
.
Chemla
S
,
Chavane
F
2010
b.
Voltage-sensitive dye imaging: technique review and models
.
J Physiol Paris
 .
104
:
40
50
.
Contreras
D
,
Destexhe
A
,
Steriade
M
1997
.
Intracellular and computational characterization of the intracortical inhibitory control of synchronized thalamic inputs in vivo
.
J Neurophysiol
 .
78
:
335
350
.
Cuntz
H
,
Forstner
F
,
Borst
A
,
Häusser
M
2010
.
One rule to grow them all: a general theory of neuronal branching and its practical application
.
PLoS Comput Biol
 .
6
:
e1000877
.
Cuntz
H
,
Forstner
F
,
Borst
A
,
Häusser
M
2011
.
The TREES toolbox - probing the basis of axonal and dendritic branching
.
Neuroinformatics
 .
9
:
91
96
.
de Kamps
M
2013
. A generic approach to solving jump diffusion equations with applications to neural populations. ArXiv e-prints. 1309.1654v2 [q-bio.NC].
De Schutter
E
,
Van Geit
W
2009
. Modeling complex neurons. In:
De Schutter
E
editor.
Computational modeling methods for neuroscientists
 . 1st ed.
MIT Press
. p.
260
283
, Chapter 11.
Deco
G
,
Jirsa
VK
,
Robinson
PA
,
Breakspear
M
,
Friston
K
2008
.
The dynamic brain: from spiking neurons to neural masses and cortical fields
.
PLoS Comput Biol
 .
4
:
e1000092
.
Denker
M
,
Roux
S
,
Lindén
H
,
Diesmann
M
,
Riehle
A
,
Grün
S.
2011
.
The local field potential reflects surplus spike synchrony
.
Cereb Cortex
 .
21
:
2681
2695
.
Destexhe
A
,
Paré
D.
1999
.
Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo
.
J Neurophysiol
 .
81
:
1531
1547
.
Di
S
,
Baumgartner
C
,
Barth
DS.
1990
.
Laminar analysis of extracellular field potentials in rat vibrissa/barrel cortex
.
J Neurophysiol
 .
63
:
832
840
.
Diesmann
M.
2013
. The road to brain-scale simulations on K. BioSupercomputing Newsletter. 8.
Ecker
AS
,
Berens
P
,
Keliris
GA
,
Bethge
M
,
Logothetis
NK
,
Tolias
AS.
2010
.
Decorrelated neuronal firing in cortical microcircuits
.
Science
 .
327
:
584
587
.
Einevoll
GT
,
Kayser
C
,
Logothetis
NK
,
Panzeri
S.
2013
.
Modelling and analysis of local field potentials for studying the function of cortical circuits
.
Nat Rev Neurosci
 .
14
:
770
785
.
Einevoll
GT
,
Lindén
H
,
Tetzlaff
T
,
Łe¸ski
S
,
Pettersen
KH.
2013
.Local field potential: biophysical origin and analysis. In:
Quiroga
RQ
,
Panzeri
S
editors.
Principles of neural coding
 .
CRC Press
. p.
37
60
, Chapter 3.
Einevoll
GT
,
Pettersen
KH
,
Devor
A
,
Ulbert
I
,
Halgren
E
,
Dale
AM.
2007
.
Laminar population analysis: estimating firing rates and evoked synaptic activity from multielectrode recordings in rat barrel cortex
.
J Neurophysiol
 .
97
:
2174
2190
.
Eppler
JM
,
Helias
M
,
Müller
E
,
Diesmann
M
,
Gewaltig
M-O.
2008
.
PyNEST: a convenient interface to the NEST simulator
.
Front Neuroinform
 .
2
:
1
12
.
Eppler
JM
,
Pauli
R
,
Peyser
A
,
Ippen
T
,
Morrison
A
,
Senk
J
,
Schenck
W
,
Bos
H
,
Helias
M
,
Schmidt
M
, et al
.
2015
. Nest 2.8.0.
Erdös
P
,
Rényi
A.
1959
.
On random graphs
.
Publications Mathematicae
 .
6
:
290
297
.
Foster
I.
1995
. Designing and building parallel programs: concepts and tools for parallel software engineering. In:
Reading
 .
Mass: Addison-Wesley Longman Publishing Co., Inc
.
Gabriel
C
,
Peyman
A
,
Grant
EH.
2009
.
Electrical conductivity of tissue at frequencies below 1 MHz
.
Phys Med Biol
 .
54
:
4863
4878
.
Gabriel
S
,
Lau
RW
,
Gabriel
C.
1996
.
The dielectric properties of biological tissues: II. measurements in the frequency range 10 Hz to 20 GHz
.
Phys Med Biol
 .
41
:
2251
2269
.
Gawne
TJ.
2010
.
The local and non-local components of the local field potential in awake primate visual cortex
.
J Comput Neurosci
 .
29
:
615
623
.
Gentet
LJ
,
Avermann
M
,
Matyas
F
,
Staiger
JF
,
Petersen
CCH.
2010
.
Membrane potential dynamics of GABAergic neurons in the barrel cortex of behaving mice
.
Neuron
 .
65
:
422
435
.
Gerstner
W
,
Naud
R.
2009
.
How good are neuron models
.
Science
 .
326
:
379
380
.
Głabska
H
,
Potworowski
J
,
Łeski
S
,
Wójcik
DK.
2014
.
Independent components of neural activity carry information on individual populations
.
PLoS One
 .
9
:
e105071
.
Gold
C
,
Henze
DA
,
Koch
C
,
Buzsáki
G.
2006
.
On the origin of the extracellular action potential waveform: a modeling study
.
J Neurophysiol
 .
95
:
3113
3128
.
Goto
T
,
Hatanaka
R
,
Ogawa
T
,
Sumiyoshi
A
,
Riera
J
,
Kawashima
R.
2010
.
An evaluation of the conductivity profile in the somatosensory barrel cortex of Wistar rats
.
J Neurophysiol
 .
104
:
3388
412
.
Hagen
E
,
Ness
TV
,
Khosrowshahi
A
,
Sørensen
C
,
Fyhn
M
,
Hafting
T
,
Franke
F
,
Einevoll
GT.
2015
.
ViSAPy: a python tool for biophysics-based generation of virtual spiking activity for evaluation of spike-sorting algorithms
.
J Neurosci Methods
 .
245
:
182
204
.
Hämäläinen
M
,
Haari
R
,
Ilmoniemi
RJ
,
Knuutila
J
,
Lounasmaa
OV.
1993
.
Magnetoencephalography - theory, instrumentation, and application to noninvasive studies of the working human brain
.
Rev Mod Phys
 .
65
:
413
496
.
Helias
M
,
Kunkel
S
,
Masumoto
G
,
Igarashi
J
,
Eppler
JM
,
Ishii
S
,
Fukai
T
,
Morrison
A
,
Diesmann
M.
2012
.
Supercomputers ready for use as discovery machines for neuroscience
.
Front Neuroinform
 .
6
:
1
12
.
Helias
M
,
Tetzlaff
T
,
Diesmann
M.
2013
.
Echoes in correlated neural systems
.
New J Phys
 .
15
:
023002
.
Herreras
O
,
Makarova
J
,
Makarov
V.
2015
.
New uses of LFPs: pathway-specific threads obtained through spatial discrimination
.
Neuroscience
 .
310
:
486
503
.
Hertz
J.
2010
.
Cross-correlations in high-conductance states of a model cortical network
.
Neural Comput
 .
22
:
427
447
.
Hill
S
,
Tononi
G.
2005
.
Modeling sleep and wakefulness in the thalamocortical system
.
J Neurophysiol
 .
93
:
1671
1698
.
Hines
ML
,
Carnevale
NT.
2001
.
Neuron: a tool for neuroscientists
.
Neuroscientist
 .
7
:
123
135
.
Hines
ML
,
Davison
AP
,
Muller
E.
2009
.
Neuron and Python
.
Front Neuroinform
 .
3
:
1
12
.
Hines
ML
,
Eichner
H
,
Schürmann
F.
2008
.
Neuron splitting in compute-bound parallel network simulations enables runtime scaling with twice as many processors
.
J Comput Neurosci
 .
25
:
203
210
.
Hines
ML
,
Morse
T
,
Migliore
M
,
Carnevale
NT
,
Shepherd
GM.
2004
.
ModelDB: a database to support computational neuroscience
.
J Comput Neurosci
 .
17
:
7
11
.
Holt
GR
,
Koch
C.
1999
.
Electrical interactions via the extracellular potential near cell bodies
.
J Comput Neurosci
 .
6
:
169
184
.
Ito
J
,
Roy
S
,
Liu
Y
,
Cao
Y
,
Fletcher
M
,
Lu
L
,
Boughter
J
,
Grün
S
,
Heck
D.
2014
.
Whisker barrel cortex delta oscillations and gamma power in the awake mouse are linked to respiration
.
Nat Commun
 .
5
:
1
10
.
Iyer
R
,
Menon
V
,
Buice
M
,
Koch
C
,
Mihalas
S.
2013
.
The influence of synaptic weight distribution on neuronal population dynamics
.
PLoS Comput Biol
 .
9
:
e1003248
.
Izhikevich
EM.
2003
.
Simple model of spiking neurons
.
IEEE Trans Neural Netw
 .
14
:
1569
1572
.
Izhikevich
EM
,
Edelman
GM.
2008
.
Large-scale model of mammalian thalamocortical systems
.
Proc Natl Acad Sci
 .
105
:
3593
3598
.
Jacobs
G
,
Claiborne
B
,
Harris
K.
2009
. Reconstruction of neuronal morphology. In:
De Schutter
E
editor.
Computational modeling methods for neuroscientists
 . 1 ed.
MIT Press
. p.
187
210
, Chapter 8.
Jiang
X
,
Shen
S
,
Cadwell
CR
,
Berens
P
,
Sinz
F
,
Ecker
AS
,
Patel
S
,
Tolias
AS.
2015
.
Principles of connectivity among morphologically defined cell types in adult neocortex
.
Science
 .
350
:
aac9462
.
Jin
J
,
Wang
Y
,
Swadlow
HA
,
Alonso
JM
2011
.
Population receptive fields of ON and OFF thalamic inputs to an orientation column in visual cortex
.
Nat Neurosci
 .
14
:
232
238
.
Jolivet
R
,
Schürmann
F
,
Berger
TK
,
Naud
R
,
Gerstner
W
,
Roth
A.
2008
.
The quantitative single-neuron modeling competition
.
Biol Cybern
 .
99
:
417
426
.
Kaiboriboon
K
,
Lüders
HO
,
Hamaneh
M
,
Turnbull
J
,
Lhatoo
SD.
2012
.
EEG source imaging in epilepsy–practicalities and pitfalls
.
Nat Rev Neurol
 .
8
:
498
507
.
Kajikawa
Y
,
Schroeder
CE.
2011
.
How local is the local field potential
.
Neuron
 .
72
:
847
858
.
Kandel
ER
,
Markram
H
,
Matthews
PM
,
Yuste
R
,
Koch
C.
2013
.
Neuroscience thinks big (and collaboratively)
.
Nat Rev Neurosci
 .
14
:
659
664
.
Kelly
RC
,
Smith
MA
,
Kass
RE
,
Lee
TS.
2010
.
Local field potentials indicate network state and account for neuronal response variability
.
J Comput Neurosci
 .
29
:
567
579
.
Kisvárday
ZF
,
Eysel
UT.
1992
.
Cellular organization of reciprocal patchy networks in layer III of cat visual cortex (area 17)
.
Neuroscience
 .
46
:
275
286
.
Kobayashi
R
,
Tsubo
Y
,
Shinomoto
S.
2009
.
Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold
.
Front Comput Neurosci
 .
3
:
1
11
.
Koch
C.
1984
.
Cable theory in neurons with active, linearized membranes
.
Biol Cybern
 .
50
:
15
33
.
Koch
C
,
Poggio
T.
1985
.
A simple algorithm for solving the cable equation in dendritic trees of arbitrary geometry
.
J Neurosci Methods
 .
12
:
303
315
.
Kriener
B
,
Enger
H
,
Tetzlaff
T
,
Plesser
HE
,
Gewaltig
M-O
,
Einevoll
GT.
2014
.
Dynamics of self-sustained asynchronous-irregular activity in random networks of spiking neurons with strong synapses
.
Front Comput Neurosci
 .
8
:
1
18
.
Kunkel
S
,
Schmidt
M
,
Eppler
JM
,
Plesser
HE
,
Masumoto
G
,
Igarashi
J
,
Ishii
S
,
Fukai
T
,
Morrison
A
,
Diesmann
M
, et al
.
2014
.
Spiking network simulation code for petascale computers
.
Front Neuroinform
 .
8
:
1
23
.
Łe¸ski
S
,
Kublik
E
,
Świejkowski
DA
,
Wróbel
A
,
Wójcik
DK.
2010
.
Extracting functional components of neural dynamics with independent component analysis and inverse current source density
.
J Comput Neurosci
 .
29
:
459
473
.
Łe¸ski
S
,
Lindén
H
,
Tetzlaff
T
,
Pettersen
KH
,
Einevoll
GT.
2013
.
Frequency dependence of signal power and spatial reach of the local field potential
.
PLoS Comput Biol
 .
9
:
e1003137
.
Lindén
H
,
Hagen
E
,
Łe¸ski
S
,
Norheim
ES
,
Pettersen
KH
,
Einevoll
GT.
2014
.
LFPy: a tool for biophysical simulation of extracellular potentials generated by detailed model neurons
.
Front Neuroinform
 .
7
:
1
15
.
Lindén
H
,
Pettersen
KH
,
Einevoll
GT.
2010
.
Intrinsic dendritic filtering gives low-pass power spectra of local field potentials
.
J Comput Neurosci
 .
29
:
423
444
.
Lindén
H
,
Tetzlaff
T
,
Potjans
TC
,
Pettersen
KH
,
Grün
S
,
Diesmann
M
,
Einevoll
GT.
2011
.
Modeling the spatial reach of the LFP
.
Neuron
 .
72
:
859
872
.
Logothetis
NK
,
Kayser
C
,
Oeltermann
A.
2007
.
In vivo measurement of cortical impedance spectrum in monkeys: implications for signal propagation
.
Neuron
 .
55
:
809
823
.
London
M
,
Häusser
M.
2005
.
Dendritic computation
.
Annu Rev Neurosci
 .
28
:
503
532
.
Maier
A
,
Adams
GK
,
Aura
C
,
Leopold
DA.
2010
.
Distinct superficial and deep laminar domains of activity in the visual cortex during rest and stimulation
.
Front Syst Neurosci
 .
4
:
1
11
.
Mainen
ZF
,
Sejnowski
TJ.
1996
.
Influence of dendritic structure on firing pattern in model neocortical neurons
.
Nature
 .
382
:
363
366
.
Makarov
VA
,
Makarova
J
,
Herreras
O.
2010
.
Disentanglement of local field potential sources by independent component analysis
.
J Comput Neurosci
 .
29
:
445
457
.
Markram
H
,
Muller
E
,
Ramaswamy
S
,
Reimann
MW
,
Abdellah
M
,
Sanchez
CA
,
Ailamaki
A
,
Alonso-Nanclares
L
,
Antille
N
,
Arsever
S
, et al
.
2015
.
Reconstruction and simulation of neocortical microcircuitry
.
Cell
 .
163
:
456
492
.
Markram
H
,
Toledo-Rodriguez
M
,
Wang
Y
,
Gupta
A
,
Silberberg
G
,
Wu
C.
2004
.
Interneurons of the neocortical inhibitory system
.
Nat Rev Neurosci
 .
5
:
793
807
.
Mazzoni
A
,
Brunel
N
,
Cavallari
S
,
Logothetis
NK
,
Panzeri
S.
2011
.
Cortical dynamics during naturalistic sensory stimulations: experiments and models
.
J Physiol Paris
 .
105
:
2
15
.
Mazzoni
A
,
Lindén
H
,
Cuntz
H
,
Lansner
A
,
Panzeri
S
,
Einevoll
GT.
2015
.
Computing the local field potential (LFP) from integrate-and-fire network models
.
PLoS Comput Biol
 .
11
:
e1004584
.
Mazzoni
A
,
Panzeri
S
,
Logothetis
NK
,
Brunel
N.
2008
.
Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons
.
PLoS Comput Biol
 .
4
:
e1000239
.
Mazzoni
A
,
Whittingstall
K
,
Brunel
N
,
Logothetis
NK
,
Panzeri
S.
2010
.
Understanding the relationships between spike rate and delta/gamma frequency bands of LFPs and EEGs using a local cortical network model
.
Neuroimage
 .
52
:
956
972
.
Migliore
M
,
Cavarretta
F
,
Hines
ML
,
Shepherd
GM.
2014
.
Distributed organization of a brain microcircuit analyzed by three-dimensional modeling: the olfactory bulb
.
Front Comput Neurosci
 .
8
:
1
14
.
Mitzdorf
U.
1985
.
Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena
.
Physiol Rev
 .
65
:
37
100
.
Mitzdorf
U
,
Singer
W.
1979
.
Excitatory synaptic ensemble properties in the visual cortex of the macaque monkey: a current source density analysis of electrically evoked potentials
.
J Comp Neurol
 .
187
:
71
83
.
Moran
RJ
,
Stephan
KE
,
Kiebel
SJ
,
Rombach
N
,
O'Connor
WT
,
Murphy
KJ
,
Reilly
RB
,
Friston
KJ
2008
.
Bayesian estimation of synaptic physiology from the spectral responses of neural masses
.
Neuroimage
 .
42
:
272
284
.
Moulin
C
,
Glière
A
,
Barbier
D
,
Joucla
S
,
Yvert
B
,
Mailley
P
,
Guillemaud
R.
2008
.
A new 3-D finite-element model based on thin-film approximation for microelectrode array recording of extracellular action potential
.
IEEE Tran. Bio-Medical Engg
 .
55
:
683
92
.
Nauhaus
I
,
Busse
L
,
Carandini
M
,
Ringach
DL.
2009
.
Stimulus contrast modulates functional connectivity in visual cortex
.
Nat Neurosci
 .
12
:
70
76
.
Nelson
MJ
,
Pouget
P.
2010
.
Do electrode properties create a problem in interpreting local field potential recordings
.
J Neurophysiol
 .
103
:
2315
2317
.
Nelson
MJ
,
Pouget
P
,
Nilsen
EA
,
Patten
CD
,
Schall
JD.
2008
.
Review of signal distortion through metal microelectrode recording circuits and filters
.
J Neurosci Methods
 .
169
:
141
157
.
Ness
TV
,
Potworowski
J
,
Łe¸ski
S
,
Chintaluri
C
,
Głąbska
H
,
Wójcik
DK
,
Einevoll
GT.
2015
.
Modelling and analysis of electrical potentials recorded in multielectrode arrays (MEAs)
.
Neuroinformatics
 .
13
(
4
):
403
426
.
Ness
TV
,
Remme
MWH
,
Einevoll
GT.
2016
.
Active subthreshold dendritic conductances shape the local field potential
.
J Physiol
 .
594
(
13
):
3809
3825
Newman
ME.
2003
.
The structure and function of complex networks
.
SIAM review
 .
45
:
167
256
.
Nicholson
C
,
Freeman
JA.
1975
.
Theory of current source-density analysis and determination of conductivity tensor for anuran cerebellum
.
J Neurophysiol
 .
38
:
356
368
.
Nordlie
E
,
Gewaltig
M-O
,
Plesser
HE.
2009
.
Towards reproducible descriptions of neuronal network models
.
PLoS Comput Biol
 .
5
:
e1000456
.
Nowak
LG
,
Azouz
R
,
Sanchez-Vives
MV
,
Gray
CM
,
McCormick
DA.
2003
.
Electrophysiological classes of cat primary visual cortical neurons in vivo as revealed by quantitative analyses
.
J Neurophysiol
 .
89
:
1541
1566
.
Nunez
PL
,
Srinivasan
R.
2006
.
Electric fields of the brain, The neurophysics of EEG
 . 2nd ed.
Oxford, UK
:
Oxford University Press, Inc
.
Ohbayashi
M
,
Ohki
K
,
Miyashita
Y.
2003
.
Conversion of working memory to motor sequence in the monkey premotor cortex
.
Science
 .
301
:
233
236
.
Okun
M
,
Lampl
I.
2008
.
Instantaneous correlation of excitation and inhibition during ongoing and sensory-evoked activities
.
Nat Neurosci
 .
11
:
535
537
.
Ostojic
S.
2014
.
Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons
.
Nat Neurosci
 .
17
:
594
600
.
Pettersen
KH
,
Devor
A
,
Ulbert
I
,
Dale
AM
,
Einevoll
GT.
2006
.
Current-source density estimation based on inversion of electrostatic forward solution: effects of finite extent of neuronal activity and conductivity discontinuities
.
J Neurosci Methods
 .
154
:
116
133
.
Pettersen
KH
,
Hagen
E
,
Einevoll
GT.
2008
.
Estimation of population firing rates and current source densities from laminar electrode recordings
.
J Comput Neurosci
 .
24
:
291
313
.
Pettersen
KH
,
Lindén
H
,
Dale
AM
,
Einevoll
GT.
2012
.Extracellular spikes and CSD. In:
Brette
R
,
Destexhe
A
editors.
Handbook of neural activity measurement
 .
Cambridge, UK
:
Cambridge University Press
. p.
92
135
.
Plesser
HE
,
Eppler
JM
,
Morrison
A
,
Diesmann
M
,
Gewaltig
M-O.
2007
. Efficient parallel simulation of large-scale neuronal networks on clusters of multiprocessor computers. In:
Kermarrec
A-M
,
Bougé
L
,
Priol
T
editors.
Euro-Par 2007: Parallel processing of lecture notes in computer science
 .
Vol. 4641
:
Springer-Verlag
. p.
672
681
.
Potjans
TC
,
Diesmann
M.
2014
.
The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model
.
Cereb Cortex
 .
24
:
785
806
.
Rall
W.
1964
.Theoretical significance of dendritic trees for neuronal input-output relations. In:
Reiss
RF
editor.
Neural theory and modeling
 .
Redwood City, CA
:
Stanford University Press
.
Rall
W.
2009
. Rall model. Scholarpedia. 4:1369. revision 91692.
Rall
W
,
Shepherd
GM.
1968
.
Theoretical reconstruction of field potentials and dendrodendritic synaptic interactions in olfactory bulb
.
J Neurophysiol
 .
31
:
884
915
.
Ray
S
,
Maunsell
JHR.
2011
.
Different origins of gamma rhythm and high-gamma activity in macaque visual cortex
.
PLoS Biology
 .
9
:
e10001610
.
Reimann
MW
,
Anastassiou
CA
,
Perin
R
,
Hill
SL
,
Markram
H
,
Koch
C.
2013
.
A biophysically detailed model of neocortical local field potentials predicts the critical role of active membrane currents
.
Neuron
 .
79
:
375
390
.
Reimann
MW
,
King
JG
,
Muller
EB
,
Ramaswamy
S
,
Markram
H.
2015
.
An algorithm to predict the connectome of neural microcircuits
.
Front Comput Neurosci
 .
9
:
1
18
.
Remme
MWH
,
Rinzel
J.
2011
.
Role of active dendritic conductances in subthreshold input integration
.
J Comput Neurosci
 .
31
:
13
30
.
Renart
A
,
De La Rocha
J
,
Bartho
P
,
Hollender
L
,
Parga
N
,
Reyes
A
,
Harris
KD.
2010
.
The asynchronous state in cortical circuits
.
Science
 .
327
:
587
590
.
Reyes-Puerta
V
,
Sun
J-J
,
Kim
S
,
Kilb
W
,
Luhmann
HJ.
2015
.
Laminar and columnar structure of sensory-evoked multineuronal spike sequences in adult rat barrel cortex in vivo
.
Cereb Cortex
 .
25
:
2001
2021
.
Reyes-Puerta
V
,
Yang
J-W
,
Siwek
ME
,
Kilb
W
,
Sun
J-J
,
Luhmann
HJ.
2016
.
Propagation of spontaneous slow-wave activity across columns and layers of the adult rat barrel cortex in vivo
.
Brain Struct Funct
 .
1
21
.
Riehle
A
,
Wirtssohn
S
,
Grün
S
,
Brochier
T.
2013
.
Mapping the spatio-temporal structure of motor cortical LFP and spiking activities during reach-to-grasp movements
.
Front Neural Circuits
 .
7
:
1
15
.
Robinson
DA.
1968
.
The electrical properties of metal microelectrodes
.
Proceedings of the IEEE
 .
56
:
1065
1071
.
Rössert
C
,
Pozzorini
C
,
Chindemi
G
,
Davison
AP
,
Eroe
C
,
King
J
,
Newton
TH
,
Nolte
M
,
Ramaswamy
S
,
Reimann
MW
, et al
.
2016
. Automated point-neuron simplification of data-driven microcircuit models. ArXiv e-prints. 1604.00087 [q-bio.NC].
Rotter
S
,
Diesmann
M.
1999
.
Exact digital simulation of time-invariant linear systems with applications to neuronal modeling
.
Biol Cybern
 .
81
:
381
402
.
Sabah
NH
,
Leibovic
KN.
1969
.
Subthreshold oscillatory responses of the Hodgkin-Huxley cable model for the squid giant axon
.
Biophys J
 .
9
:
1206
1222
.
Sakata
S
,
Harris
KD.
2009
.
Laminar structure of spontaneous and sensory-evoked population activity in auditory cortex
.
Neuron
 .
64
:
404
418
.
Sarid
L
,
Bruno
R
,
Sakmann
B
,
Segev
I
,
Feldmeyer
D.
2007
.
Modeling a layer 4-to-layer 2/3 module of a single column in rat neocortex: interweaving in vitro and in vivo experimental observations
.
Proc Natl Acad Sci
 .
104
:
16353
16358
.
Schmidt
M
,
Bakker
R
,
Shen
K
,
Bezgin
G
,
Hilgetag
C-C
,
Diesmann
M
,
van Albada
SJ.
2016
.
Full-density multi-scale account of structure and dynamics of macaque visual cortex
ArXiv e-prints. 1511.09364v4 [q-bio.NC].
Schreiber
T
,
Schmitz
A.
2000
.
Surrogate time series
.
Phys D
 .
142
:
346
382
.
Schuecker
J
,
Diesmann
M
,
Helias
M.
2015
.
Modulated escape from a metastable state driven by colored noise
.
Physical Review E
 .
92
:
052119
.
Senk
J
,
Hagen
E
,
van Albada
SJ
,
Diesmann
M.
2015
. From randomly connected to spatially organized multi-layered cortical network models. In: Proceedings of the 11th Göttingen Meeting of the German Neuroscience Society. pp.
1126
1127
.
Shadlen
MN
,
Newsome
WT.
1998
.
The variable discharge of cortical neurons: implications for connectivity, computation, and information coding
.
J Neurosci
 .
18
:
3870
3896
.
Softky
WR
,
Koch
C.
1993
.
The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs
.
J Neurosci
 .
13
:
334
350
.
Song
S
,
Sjöström
PJ
,
Reigl
M
,
Nelson
S
,
Chklovskii
DB.
2005
.
Highly nonrandom features of synaptic connectivity in local cortical circuits
.
PLoS biology
 .
3
:
e350
.
Stepanyants
A
,
Hirsch
JA
,
Martinez
LM
,
Kisvárday
ZF
,
Ferecskó
AS
,
Chklovskii
DB.
2008
.
Local potential connectivity in cat primary visual cortex
.
Cereb Cortex
 .
18
:
13
28
.
Swadlow
HA
,
Gusev
AG
,
Bezdudnaya
T.
2002
.
Activation of a cortical column by a thalamocortical impulse
.
J Neurosci
 .
22
:
7766
7773
.
Szymanski
FD
,
Garcia-Lazaro
JA
,
Schnupp
JWH.
2009
.
Current source density profiles of stimulus-specific adaptation in rat auditory cortex
.
J Neurophysiol
 .
102
:
1483
1490
.
Szymanski
FD
,
Rabinowitz
NC
,
Magri
C
,
Panzeri
S
,
Schnupp
JWH.
2011
.
The laminar and temporal structure of stimulus information in the phase of field potentials of auditory cortex
.
J Neurosci
 .
31
:
15787
15801
.
Telenczuk
B
,
Dehghani
N
,
Le Van Quyen
M
,
Cash
SS
,
Halgren
E
,
Hatsopoulos
NG
,
Destexhe
A.
2016
. Local field potentials primarily reflect inhibitory neuron activity in human and monkey cortex. bioRxiv.
Teramae
J-N
,
Fukai
T.
2014
.
Computational implications of lognormally distributed synaptic weights
.
Proc. IEEE
 .
102
:
500
512
.
Tetzlaff
T
,
Helias
M
,
Einevoll
GT
,
D