## Abstract

We present a new method, laminar population analysis (LPA), for analysis of laminar-electrode (linear multielectrode) data, where physiological constraints are explicitly incorporated in the mathematical model: the high-frequency band [multiunit activity (MUA)] is modeled as a sum over contributions from firing activity of multiple cortical populations, whereas the low-frequency band [local field potential (LFP)] is assumed to reflect the dendritic currents caused by synaptic inputs evoked by this firing. The method is applied to stimulus-averaged laminar-electrode data from barrel cortex of anesthetized rat after single whisker flicks. Two sample data sets, distinguished by stimulus paradigm, type of applied anesthesia, and electrical boundary conditions, are studied in detail. These data sets are well accounted for by a model with four cortical populations: one supragranular, one granular, and two infragranular populations. Population current source densities (CSDs; the CSD signatures after firing in a particular population) provided by LPA are further used to estimate the synaptic connection pattern between the various populations using a new LFP template-fitting technique, where LFP population templates are found by the electrostatic forward solution based on results from compartmental modeling of morphologically reconstructed neurons. Our analysis confirms previous experimental findings regarding the synaptic connections from neurons in the granular layer onto neurons in the supragranular layers and provides predictions about other synaptic connections. Furthermore, the time dependence of the stimulus-evoked population firing activity is predicted, and the temporal ordering of response onset is found to be compatible with earlier findings.

## INTRODUCTION

Extracellular recordings using multielectrode arrays are uniquely suited to bridge the gap between classical single-cell electrophysiology and noninvasive extracranial EEG/MEG recordings. However, the relationship between extracellular field potentials and population firing rates, and the behavior of specific neuronal cell types and circuits, is poorly understood. In this paper we describe a novel framework for estimation of spiking and synaptic activity in identified cortical neuronal populations, based on laminar array recordings. This approach also provides a means to constrain and validate computational models of neuronal cortical circuits (Dayan and Abbott 2001; Koch and Segev 1998).

Technology for large-scale electrical recording of neural activity is rapidly improving, enabling routine measurement of electric potentials from large numbers of closely spaced locations (Buzsaki 2004; Nadasdy et al. 1998). Local field potentials (LFP), i.e., the low-frequency band of extracellularly recorded potentials, are thought to predominantly stem from dendritic processing of synaptic inputs, but a direct interpretation in terms of the underlying neural activity is difficult (Nicholson and Freeman 1975; Pettersen et al. 2006). A standard measurement procedure has been to record the field potential at equidistant, linearly positioned electrode contacts using a laminar electrode vertically penetrating the cortical layers (Mitzdorf 1985; Nadasdy et al. 1998; Rappelsberger et al. 1981; Ulbert et al. 2001). The current source density (CSD) can be estimated from the second spatial derivative of the recorded LFP (Nicholson and Freeman 1975; Rappelsberger et al. 1981) or by using the newly developed inverse CSD (iCSD) method (Pettersen et al. 2006). However, the CSD measures the net transmembrane current provided by all the underlying neural populations and the contributions from individual populations are difficult to assess. In a series of papers, Barth and collaborators pursued such an assessment by the application of principal component analysis (PCA) of the CSD (Barth and Di 1991; Barth et al. 1989, 1990; Di et al. 1990), and in their study of stimulus-evoked data from the rat barrel cortex, they identified putative cortical populations in supra- and infragranular layers of the barrel column (Di et al. 1990).

PCA is one of several statistical methods where functions of two variables (here electrode contact positions and time) are expanded into sums over spatiotemporally separable functions, i.e., functions that can be written as a product of a spatial function and a temporal function. Such an expansion can be done in an infinite number of ways, and additional constraints are needed to make the expansion unique. In PCA, the functions (i.e., components) are constrained to be orthogonal both in space and time. In this paper, we instead use physiological constraints to specify the expansion. As a consequence, the expansion will not only be compatible with the physiology, but each component will also have a clear physiological interpretation.

The basic underlying constraint inherent in the laminar population analysis (LPA) presented here is that the observed LFP is evoked by the firing of action potentials in the modeled neuronal populations. An experimental measure of this population firing is obtained from the high-frequency component of the laminarly recorded extracellular potentials, i.e., the multiunit activity (MUA) (Ulbert et al. 2001). The outcomes of the laminar population analysis are *1*) identification of the relevant laminar cortical populations and their vertical spatial position and extent, *2*) estimates of the firing-rates of these populations, and *3*) estimates of the spatiotemporal LFP signature (and CSD signature) after action potential firing in the individual populations.

Here we apply LPA to stimulus-evoked laminar-electrode data from rat somatosensory (barrel) cortex, a well-studied example of topographic mapping where each of the large facial vibrissa (whiskers) is mapped onto a specific barrel column (Woolsey and Van der Loos 1970). Several sample data sets, distinguished by stimulation paradigm, anesthesia conditions, and electrical boundary conditions at the cortical surface, are considered. Despite the different experimental situations, LPA seems to identify the same set of cortical populations (typically 1 supragranular, 1 granular, and 2 infragranular). Furthermore, the CSD signatures after population firing are generally found to be similar between experiments, in particular for the three topmost populations accounting for most of the observed LFP.

While these population CSDs give more information about the underlying neural circuit activity than traditional CSD analysis, they still do not provide a complete understanding of the synaptic connection pattern between the various populations. To estimate this pattern in this experimental situation, we use a new template-fitting analysis where LFP population templates are fitted to the observed LFP. These LFP population templates are calculated using the electrostatic forward solution with transmembrane currents obtained from compartmental modeling of morphologically reconstructed neurons. The results are found to agree with previous findings about the projection from the granular population onto the layer 2/3 supragranular population (Lübke et al. 2000; Petersen and Sakmann 2001) and provide predictions about other synaptic connections.

Preliminary results from this project were presented earlier in poster format (Einevoll et al. 2004).

## METHODS

### Experimental procedure and initial data processing

All experimental procedures were approved by the Massachusetts General Hospital Subcommittee on Research Animal Care.

Male Sprague-Dawley rats (250–350 g, Taconic) were used in the experiments. Glycopyrrolate (0.5 mg/kg, im) was administered 10 min before the initiation of anesthesia. Two anesthesia regimens were used. *1*) Rats were anesthetized with 1.5% halothane in oxygen for surgery. After surgery, halothane was discontinued, and anesthesia was maintained with 50 mg/kg intravenous bolus of α-chloralose followed by continuous intravenous infusion at 40 mg/kg/h. *2*) Rats were anesthetized with a single dose of ketamine hydrochloride (50 mg/kg, ip), supplemented with sodium pentobarbital (12 mg/kg) during the surgery. During the experiment, animals were supplemented by constant intraperitoneal infusion of 50 mg/kg/h of ketamine hydrochloride. During surgery, a tracheotomy was performed, and cannulas were inserted in the femoral artery and vein. All incisions were infiltrated with 2% lidocaine. After tracheotomy, rats were mechanically ventilated with 30% O_{2} in air. Ventilation parameters were adjusted to maintain Paco_{2} between 35 and 45 mmHg, Pao_{2} between 140 and 180 mmHg, and pH between 7.35 and 7.45. Heart rate, ECG, and blood pressure were monitored continuously. Body temperature was maintained at 37.0 ± 0.5°C with a homeothermic blanket (Harvard Apparatus, Holliston, MA). The animal was fixed in a stereotaxic frame. An area of skull overlying the primary somatosensory cortex was exposed and thinned with a dental burr. The thinned skull was removed, and the dura matter was dissected to expose the cortical surface. A barrier of dental acrylic was built around the border of the exposure and filled with either saline or mineral oil. Mapping with a single metal microelectrode (FHC, 2–5 MΩ) was done to determine the positioning of the laminar electrode array. The optimal position was identified by listening to an audio monitor while stimulating different whiskers. After the mapping procedure, the electrode was withdrawn, and the array was slowly introduced at the same location perpendicular to the cortical surface. Contact 1 was positioned at the cortical surface using visual control. A linear array multielectrode with 23 contacts with diameter 0.04 mm spaced at 0.1 mm was used (Ulbert et al. 2001), and the data from contacts 2–23 were used in the modeling.

Single whiskers were deflected upward by a wire loop coupled to a computer-controlled piezoelectric stimulator. We used a fast, randomized event-related stimulus presentation paradigm. The stimulus sequence was optimized for event-related response-estimation efficiency using the approach described by Dale (1999). The stimulation paradigm consisted of single deflections of varying amplitude with interstimulus interval (ISI) of 1 s.

We used three stimulus conditions. *1*) In the nine-stimulus condition case, the vertical displacement of the whisker varied from 0.12 (amplitude 1) to 1.08 mm (amplitude 9). Intervening amplitudes had whisker displacements with equal increments on a linear scale. The stimulus rise time (time from onset to maximum displacement) was 30 ms, and the (mean) stimulus angular velocity increased from 76 (amplitude 1) to 660°/s (amplitude 9). *2*) In the second stimulus condition, with altogether 27 different stimuli, three stimulus amplitudes [vertical displacements 0.40 (*a1*), 0.80 (*a2*), and 1.2 mm (*a3*)] each with nine stimulus rise times [20 (*t1*), 30 (*t2*), . ., 100 ms (*t9*)] were applied. The stimulus angular velocity varied between 76 (*a1,t9*) and 1,090°/s (*a*3,*t*1). *3*) In the third stimulus condition, a stimulus rise time of 30 ms was used, with 27 different, linearly incrementing, vertical displacements ≤1.2 mm.

The recorded potential was amplified and filtered into two signals: a low-frequency part (0.1–500 Hz, sampled at 2 kHz with 16 bits) and a high-frequency part (150–5,000 Hz, sampled at 20 kHz with 12 bits; see Ulbert et al. 2001 for details). The low-frequency part is referred to as the LFP. The high-frequency part was further filtered digitally between 750 and 5,000 Hz using a zero phase-shift second-order Butterworth filter and rectified along the time axis to provide the MUA. Finally, the MUA data were decimated by a factor 10 along the time axis to provide the same time resolution (0.5 ms) as the LFP data. Note that the temporal form of the MUA was found to be essentially the same without decimation or when 500 or 1,000 Hz was used as the lower cut-off frequency in the band-pass filtering. Example traces of LFP and MUA (before rectification) for a single trial are shown in Fig. 1.

Stimulus-averaged data were obtained by averaging over all trials for each stimulus type (40 trials for experiments with 27 stimulus conditions, 120 trials for experiment with 9 stimulus conditions). Baseline-subtracted data were used in the modeling, and this was obtained by subtracting the mean signal for a 50- or 100-ms period before stimulus onset for each electrode contact channel separately.

In this paper we consider four experiments. In *experiment 1*, α-chloralose was used as anesthesia, the stimulus-condition paradigm with three amplitudes and nine rise times was applied, and saline covered the exposed cortex around the laminar electrode. In *experiment 2*, ketamine was used as anesthesia, the nine stimulus-condition paradigm was applied, and oil covered the exposed cortex. In *experiment 3*, α-chloralose, saline, and the 27-amplitude stimulus condition was used, and in *experiment 4*, α-chloralose, saline and the three amplitude, nine rise-time stimulus condition was used. *Experiment 1b* corresponds to a separate 27-amplitude experiment performed on the rat in *experiment 1*. Eight of 5,400 trials were removed because of artifacts in the MUA. In *experiments 2–4*, the MUA signal from contact 9 was erratic, and for this channel, the arithmetic mean of the MUA signals from the neighboring contacts (contacts 8 and 10) was used instead.

The stimulus-averaged, baseline-subtracted LFP and MUA data used in the modeling are denoted as *φ*_{L}(*z*_{i},*t*_{j}) and *φ*_{M}(*z*_{i},*t*_{j}), respectively. Here *z*_{i} refers to electrode contact position and *t*_{j} to time. The time resolution is 0.5 ms.

### LPA

Two-dimensional data like our LFP [*φ*_{L}(*z*_{i},*t*_{j})] and MUA [*φ*_{M}(*z*_{i},*t*_{j})] data can be expanded in a series of spatiotemporally separable functions (1)

In PCA, the expansion functions *f*_{n}(*z*_{i}) (called spatial loading) and *g*_{n}(*t*_{j}) (called temporal scores) are chosen so that the first component picks up most of the data variance, the second component most of the remaining variance, and so on. Furthermore, the PCA algorithm forces both the spatial loadings [*f*_{n}(*z*_{i})] and the temporal scores [*g*_{n}(*t*_{j})] to be orthogonal, i.e., , for *n* ≠ *m* (Gershenfeld 1999).

In LPA, the MUA and LFP data are modeled jointly. First, we assume the MUA data to be modeled as a sum over spatiotemporally separable contributions from several neuronal populations (2) where the superscript *m* is added to distinguish model from data. Here *N*_{pop} is the number of populations, *r*_{n}(*t*_{j}) is the population firing rate of neuronal population *n*, and *M*_{n}(*z*_{i}) is the corresponding MUA spatial profile after action potential firing in this population. This spatial profile will depend on the physiological properties of the neurons in the population, the distribution of their spatial positions in the cortical lamina, and the properties of the electrical conductivity in the extracellular electrical medium.

Because the MUA data are rectified, we impose the constraint that the spatial functions *M*_{n}(*z*_{i}) should be nonnegative. The interpretation of *r*_{n}(*t*_{j}) as the population firing rate should imply that these functions should also be nonnegative. However, in our analysis of stimulus-averaged data, we subtract the baseline, i.e., mean MUA before stimulus onset. This means that *r*_{n}(*t*_{j}) should be interpreted as the population firing rate measured relative to baseline and should thus be allowed to have small negative values.

Next, we assume that the LFP data are driven by the population firing activities *r*_{n}(*t*_{j}) seen in the MUA; firing of an action potential of a neuron in population *n* will lead to postsynaptic transmembrane currents (including both the ligand-gated synaptic currents and consequent return currents). These transmembrane currents will in turn contribute to the LFP (Nicholson and Freeman 1975). Consequently, we assume that the LFP data can be decomposed into contributions from each of the neuronal populations, and we suggest the following form (3) Here (*h*_{n} ⊗ *r*_{n})(*t*_{j}) is the temporal convolution given by (4) *L*_{n}(*z*_{i}) represents the spatial profile of the contribution to the LFP data after action potential firing in population *n*. The temporal coupling kernel *h*_{n}(*t*_{k}) accounts for the temporal delay and spread in the generation of the LFP after firing of neurons in population *n,* and from causality, it follows that *h*_{n}(*t* < 0) = 0. The LFP is thought mainly to be caused by synaptic currents and their associated membrane response of the postsynaptic neurons. The form of this temporal transfer function will thus depend both on synaptic delays, synaptic time constants, and physiological properties of the postsynaptic neurons. The task is now to identify the choice of spatial profiles [*M*_{n}(*z*_{i}), *L*_{n}(*z*_{i})], temporal transfer functions *h*_{n}(*t*_{k}), and firing-rate functions *r*_{n}(*t*_{j}) [given the above-mentioned nonnegativity constraint on *M*_{n}(*z*)] to make the deviations between the model expressions φ_{M}^{m}*(z*_{i},*t*_{j}) and φ_{L}^{m}*(z*_{i},*t*_{j}) and the experimental data φ_{M}(*z*_{i},*t*_{j}) and φ_{L}(*z*_{i},*t*_{j}) as small as possible.

##### FORM OF SPATIAL POPULATION PROFILES *L*_{n}(*z*_{i}) and *M*_{n}(*z*_{i}).

The spatial profiles *L*_{n}(*z*_{i}) are allowed to vary freely. However, for the MUA spatial profiles *M*_{n}(*z*_{i}), we assume particular parameterized functional forms. The exact profiles will depend both on the morphology of the neurons in the population and the spatial distribution of the neurons constituting the population and are difficult to know in detail a priori. However, the results of Somogyvari et al. (2005) for the vertical extent of the extracellular potential after an action potential (measured with the same type of laminar electrode) indicate that the MUA signal is mainly caused by sums of action potential contributions from neurons with their somas positioned within a distance 0.1 mm above or below the recording position. With the further assumption that each population is spatially continuous and covers a certain depth range in cortex, this implies that the MUA spatial profiles of the different populations will be largely nonoverlapping and that the individual spatial profiles will be spatially confined (bump-like) with rather sharp interface slopes and limited overlap to the neighboring populations. We thus model *M*_{n}(*z*_{i}) as a trapezoidal symmetric profile centered at the position *z*_{0n}_{.} The profile has a constant height between *z*_{0n} − *a*_{n}/2 and *z*_{0n} + *a*_{n}/2 and decays linearly to zero over a length *b*_{n} on each side. We here constrain the flat maxima of the spatial profiles *M*_{n}(*z*_{i}) to be nonoverlapping, (e.g., *z*_{0n} + *a*_{n}/2 < *z*_{0N} − *a*_{N}/2 if population *N* is positioned immediately below population *n* in cortex). We further constrain *b*_{n} to be smaller than 0.1 mm to enforce a rather sharp interface between neighboring populations.

##### FORM OF TEMPORAL COUPLING KERNEL *h*_{n}(*t*_{k}).

This kernel can be thought of as the time-course of the effect of an action potential on the LFP, and one could in principle derive the functional form from knowledge of the electrophysiological effects of synaptic inputs onto neuronal dendrites. However, as described in the Results, our sample *experiment 2* offers a way to estimate a suitable form directly from the experimental data. The resulting form is found to be (5) This is a (normalized) delayed exponential, i.e., the temporal kernel of the delayed low-pass RC filter, where Θ(*t*) is the Heaviside unit step function [Θ(*t* ≥ 0) = 1, Θ(*t* < 0) = 0], and τ_{n} and Δ_{n} are the time constant and a time-delay parameter, respectively, of population *n*. In this application, we make the additional simplifying assumption that all temporal kernels are identical, i.e., same τ_{n} and Δ_{n} for all populations so that *h*_{n}(*t*_{k}) = *h*(*t*_{k}).

##### NUMERICAL OPTIMIZATION PROCEDURE.

In the first step of the LPA procedure, we fit the MUA data to the expression in *Eq. 2*. Each spatial MUA profile *M*_{n}(*z*) is described by three parameters (*z*_{0n}, *a*_{n}, *b*_{n}). Thus in the numerical optimization scheme 3 × *N*_{pop} parameters are varied to minimize the deviation *e*_{M} between the experimental MUA data and the model expression. Here *e*_{M} is the relative mean square error (6) where *N*_{ch} = 22 is the number of data channels (electrode contacts), and *N*_{t} is the number of data time-points. Only the first 100 ms of data after stimulus onset is considered, i.e., 201 data time points for each stimulus condition.

In the second step of the procedure, we fit the LFP data given the fitted firing rates *r*_{n}(*t*_{j}) from the MUA fit. The spatial LFP profiles *L*_{n}(*z*_{i}) are fitted nonparametrically, so the parameter fit involves only the two model parameters of *h*(*t*_{j}), the time constant τ, and the delay Δ. The optimized solution minimizes the relative mean square error, defined in analogy to *Eq. 6*. The details of the numerical optimization procedure are given in appendix 1.

### CSD estimation

The recently developed iCSD method was used to estimate the population CSD *C*_{n}(z_{i}) from the population LFP profiles *L*_{n}(*z*_{i}) (Pettersen et al. 2006). This method is based on the explicit inversion of the electrostatic solutions and can incorporate finite lateral extents of the underlying CSD and discontinuities in the extracellular conductivity at the cortical surface. Here we used the δ-source iCSD method, where the CSD sources are assumed positioned at the electrode contacts and distributed in discs of infinitesimal thickness with lateral extent corresponding to the barrel diameter, combined with a three-point Hamming spatial filter (Pettersen et al. 2006). For the case where the discs have infinite lateral extent and the extracellular conductivity is homogeneous, the δ-source iCSD method was shown in Pettersen et al. (2006) to correspond to the standard CSD method involving the double spatial derivative (Nicholson and Freeman 1975).

### Template-fitting analysis of LFP

More specific information than the CSD can be obtained by taking advantage of the fact that the forward electrostatic solutions, i.e., the LFP generated by synaptic inputs onto a particular dendritic structure (of a morphologically reconstructed neuron), can be evaluated: The transmembrane currents, i.e., CSD, after onset of synaptic input follow immediately from compartmental modeling of neurons, and given the extracellular conductivity, the LFP at any spatial position can be calculated (Gold et al. 2006; Holt and Koch 1999). This approach can be generalized to obtain corresponding LFP signatures for synaptic activation of a population of neurons. The observed LFP can be decomposed as a linear sum over such LFP population templates, and information about the synaptic connection pattern can be extracted from the fitted weights of the linear sum. Here we use this approach, denoted LFP template-fitting analysis, to obtain information about the synaptic connections between the laminar populations identified by LPA. For a detailed description, see appendix 2.

## RESULTS

### Experimental data

In Fig. 2, we show the stimulus-averaged LFP and MUA data for *experiment 1* (α-chloralose/saline/27 stimulus conditions) from the time of stimulus onset until 100 ms after onset. Results for 9 of the 27 stimulus conditions are shown: three different stimulus rise times (*t9,t5,t1*) for all three amplitudes (*a1–a3*). The shortest rise time (*t1*) gives the largest responses, whereas the longest (*t9*) gives the smallest responses. From Fig. 2, we observe that the stimulus rise time has a stronger influence on the response than does the stimulus amplitude. The corresponding stimulus-averaged LFP and MUA data for *experiment 2* (ketamine/oil/9 stimulus conditions) are shown in Fig. 3. Here, results for all nine amplitude stimulus conditions are shown (*a1–a9*).

Comparison of Figs. 2 and 3 reveals a significant qualitative difference between the LFP data for *experiment 1* compared with *experiment 2*: the LFP of *experiment 1* is generally negative, whereas the LFP of *experiment 2* has a clear positivity immediately below the cortical surface. This presumably simply reflects the difference in electrical boundary conditions at the cortical surface. Saline has a higher electrical conductivity than cortical tissue, and the saline covering the cortical surface around the inserted laminar electrode in *experiment 1* thus acts to short circuit the electrical potential toward the zero value at the distant reference electrode. With the oil cover in *experiment 2*, a generic surface positivity caused by, say, a source-sink pair below the cortical surface may instead be enhanced (because oil has a much lower electrical conductivity than cortical tissue). For further elucidation of the effects on the LFP from various surface discontinuities in the electrical conductivity, see Pettersen et al. (2006) (in particular their Fig. 4).

For both experiments 1 and 2, the earliest MUA signals appear ∼0.6–0.8 mm below the cortical surface. This MUA signal is expected to be dominated by neurons with a vertical displacement of less than ∼0.1 mm from the electrode contact position (Somogyvári et al. 2005). Thus this earliest MUA signal is interpreted as being largely caused by action potentials from layer 4 neurons, i.e., granular cells. This interpretation is supported by previous single-unit studies with single flick stimuli by Armstrong-James et al. (1992): they found that the first action potential firing is from layer 4 neurons positioned at a cortical depth between 0.5 and 0.75 mm.

Comparison of Figs. 2 and 3 also reveals that the maximum LFP and MUA signals are larger in *experiment 1* than in *experiment 2*. This may reflect differences in electrode placement and type and depth of anesthesia. We further see that the MUA in *experiment 1* has a longer duration and extends somewhat deeper down in cortex than in *experiment 2*.

### Principal component analysis

Di et al. (1990) used PCA in the interpretation of laminar-electrode LFP data from the barrel cortex of anesthetized rats using a similar stimulus paradigm. From stimulus-averaged LFP data, they derived an estimate of the CSD by performing a smoothed discrete spatial second derivative (Nicholson and Freeman 1975). They found that two principal components (PCs) accounted for essentially all (98%) of the observed variance in CSD data obtained by averaging results from five animals. Although the direct interpretation of the PCs is difficult, PCA is a useful first analysis step for elucidating the effective dimensionality of the data, providing a hint at the model complexity needed and warranted in the further analysis.

In Fig. 4, we show the results from applying PCA on the stimulus-averaged LFP and MUA data in *experiment 1* including all 27 stimulus conditions. For the LFP data, the first PC accounts for 90.6% of the variance and the second PC for 7.5%, and together they thus account for ∼98% of the data variance. For the MUA, the contributions from the first two PCs are 92.0% and 2.6%, respectively. To facilitate a comparison with the CSD analysis of Di et al. (1990), we also did PCA of CSD data obtained by using a smoothed double spatial derivative (cf. Ulbert et al. 2001) on the LFP data. Here the two first PCs were found to account for 70.8% and 18.3% of the data variance, respectively. The spatial loadings and temporal scores of the first two PCs resulting from the application of PCA on the LFP, MUA, and CSD data are shown in Fig. 4.

In Fig. 5, we show the corresponding results from applying PCA on the data in *experiment 2* including data for all nine stimulus conditions. Here the first two PCs account for 88.0% and 8.7% of the variance for the LFP data, 87.5% and 6.0% for the MUA data, and 89.7% and 7.2% for the CSD data. Comparison of the spatial loadings and temporal scores of the LFP and MUA data with the corresponding results for *experiment 1* in Fig. 4 show many similar features. The main differences in the spatial loadings occur high up in the cortex, consistent with the different electrical boundary conditions at the cortical surface. The differences in temporal scores are somewhat larger. For example, the first PC for the LFP in *experiment 1* has a shorter time-course than in *experiment 2*, whereas the situation is somewhat reversed for the MUA. For the CSD, the spatial loadings of the first PC have a similar shape for experiments 1 and 2, but the dominant sink extends deeper into cortex for *experiment 1*. For the second CSD PC, the difference is even larger, a difference also reflected in the difference in variance accounted for by this PC (18.3% for *experiment 1* vs. 7.2% for *experiment 2*).

A comparison with the first two CSD PCs obtained by Di et al. (1990) shows a qualitative resemblance with the corresponding CSD components from our data, in particular with *experiment 2* where ketamine was used. Di et al. (1990) used a combination of ketamine, xylazine, and atropine sulfate as anesthesia. The application of PCA to our sample laminar data shows, as in Di et al. (1990), that a few components can account for most of the data variance. Thus the data seem to be low-dimensional, suggesting that a laminar population model with a modest number of populations should be able to account for the data.

### LPA

##### ESTIMATION OF TEMPORAL COUPLING KERNEL.

In LPA, the functional form of the temporal coupling kernel *h*_{n}(*t*_{j}) used in *Eq. 3* must be specified. Data from *experiment 2* offer a way to estimate it; close inspection of the MUA contour plots in Fig. 3*B* shows that the MUA response for the smallest stimulus (*a1*) is not only weaker and delayed compared with the higher-stimulus responses, it also seems to be spatially more confined. A detailed comparison with the corresponding LFP response in Fig. 3*A* further shows that the onset of the LFP response is delayed compared with the MUA onset. The spatial position of the MUA response (maximum response at ∼0.7 mm cortical depth) supports an interpretation that this response stems from layer 4 cells, and this suggests that the subsequent LFP is caused by synaptic inputs from these layer 4 cells. The absence of any significant MUA response after this LFP indicates that these synaptic inputs are not sufficiently strong to evoke significant firing of the postsynaptic cells for this particularly weak stimulus (Brecht et al. 2003). The lack of any observed LFP before the onset of MUA in layer 4 indicates that the magnitude and spatial distribution of the thalamic synaptic input (as well as the morphological properties of the postsynaptic cells) are such that their contribution to the present LFP is small.

Our hypothesis is thus that for the stimulus condition *a1* firing in the layer 4 population dominates the MUA and that this firing drives the subsequently observed LFP. The time-courses of these can be estimated by performing a PCA decomposition keeping only the first and thus dominant PCs. The temporal scores of these PCs, denoted *g*_{M1}(*t*_{j}) for the MUA and *g*_{L}_{1}(*t*_{j}) for the LFP, can be used to obtain a direct estimate of the temporal coupling kernel *h*(*t*_{j}). Figure 6 shows that the temporal coupling kernel for this situation is excellently modeled by the delayed exponential kernel in *Eq. 5* with the parameters τ = 16.4 ms and Δ = 0.6 ms. We also note that the peak in the temporal score of the LFP PC *g*_{L}_{1}(*t*_{j}) occurs about 5 ms after the peak in corresponding MUA temporal score *g*_{M}_{1}(*t*_{j}). These results seem to support the fundamental assumption of the model that the presently observed LFP predominantly results from MUA and that their coupling can be described using the simple delayed exponential temporal coupling kernel *h*_{n}(*t*_{j}) given in *Eq. 5*. We will therefore use this in the further analysis of the data from all stimuli.

##### APPLICATION OF LPA TO SAMPLE DATA.

In Fig. 7 *C*, we show the results (for 9 of the 27 stimulus conditions) of fitting the MUA model (*Eq. 2*) to data from *experiment 1* assuming four cortical populations, i.e., *N*_{pop} = 4. The fitting error (*Eq. 6*) is *e*_{M} = 0.055, i.e., the model accounts for almost 95% of the MUA data variance. A visual comparison with the experimental data in Fig. 2*B* shows the overall good fit, and this is further shown in Fig. 7*D* showing the difference between the experimental data and the fitted model results. The corresponding fitted MUA spatial profiles *M*_{n}(*z*_{i})*, n* = 1, 2, 3, 4, and the corresponding population firing rates *r*_{n}(*t*_{j}) are shown in Fig. 8, *A* and *C*. In the following, we number the four identified populations from the top: population 1 (black) is found to be centered at *z*_{01} = 0.423 mm cortical depth, population 2 (blue) is centered at *z*_{02} = 0.714 mm, population 3 (red) is centered at *z*_{03} = 1.155 mm, and population 4 (green) is centered at *z*_{04} = 1.796 mm. Note that we have chosen the maximum values of all spatial profiles *M*_{n}(*z*_{i}) to be unity.

In Fig. 7*A* we show results for the fitted LFP model (*Eq. 3*) for *experiment 1* based on the fitted firing rates *r*_{n}(*t*_{j}) obtained from the MUA model fit in Fig. 7*C*. The fitting error *e*_{L} (defined in analogy with *Eq. 6*) is 0.10, i.e., the model accounts for 90% of the LFP data variance). Figure 7*B* shows the difference between the experimental LFP data and the fitted model results. Overall, the deviation is small, but it is more prominent for the stimuli providing the largest responses (e.g., *a3,t1*). For these stimuli, the actual LFP seems to increase slightly more rapidly than the model predicts. The fitted model parameters for the temporal coupling kernel (*Eq. 5*) are τ = 13.4 ms and Δ = 2.0 ms.

The fitted LFP spatial profiles *L*_{n}(*z*_{i}), *n* = 1, 2, 3, 4, are shown in Fig. 8*B*, and in Fig. 8*D* we show the time-courses of the postsynaptic LFP activity, (*h* ⊗ *r*_{n})(*t*_{j}), for the four populations. These time-courses are obtained by convolving the population firing rates *r*_{n}(*t*_{j}) in Fig. 8*C* with the delayed, exponentially decaying kernel in *Eq. 5* using the fitted parameter values for τ and Δ. This convolution effectively corresponds to a low-pass filtering of the population firing rates, and comparison of (*h* ⊗ *r*_{n})(*t*_{j}) in Fig. 8*D* with *r*_{n}(*t*_{j}) in Fig. 8*C* thus shows much smoother time-courses. The contribution to the LFP from each population *n* is determined by (*h* ⊗ *r*_{n})(*t*_{j}) multiplied by the corresponding *L*_{n}(*z*_{i}). From Fig. 8, *B* and *D,* we thus see that populations 2 (blue) and 3 (red) give large contributions to the LFP, whereas population 4 (green) generally provides a meager contribution.

In Fig. 9, *C* and *A,* we show results for the fitted MUA and LFP models, respectively, for *experiment 2*, assuming four cortical populations. The same procedure as for the model fitting to *experiment 1* was used, and the fitting errors were found to be *e*_{M} = 0.066 and *e*_{L} = 0.087, respectively. The corresponding difference between the experimental data and the fitted model results are shown in Fig. 9, *D* and *B*. Overall, a very good agreement is observed. The corresponding fitted spatial profiles for the MUA and LFP [*M*_{n}(*z*_{i}) and *L*_{n}(*z*_{i}), *n* = 1, 2, 3, 4] are shown in Fig. 10, *A* and *B*, whereas the fitted population firing rates *r*_{n}(*t*_{j}) and the postsynaptic time-courses (*h* ⊗ *r*_{n})(*t*_{j}) are shown in Fig. 10, *C* and *D*, respectively.

The spatial MUA profiles predicted from *experiment 1* (Fig. 8*A*) and *experiment 2* (Fig. 10*A*) are seen to be very similar. The main difference is the shape of the profile for the top population: for *experiment 1*, *M*_{1}(*z*_{i}) is seen to decay to zero at the top, whereas for *experiment 2*, *M*_{1}(*z*_{i}) is seen to extend to the topmost channel. This is consistent with the predicted effect of the differences in electrical boundary conditions, i.e., saline covering the cortical surface in *experiment 1* and oil in *experiment 2*.

With the exception of population 4 (green), the predicted LFP spatial profiles *L*_{n}(*z*_{i}) for the two experiments are seen to share most of the qualitative features (cf. Figs. 8*B* and 10*B*). The main difference occurs at the top where the influence on the shape of the LFP profiles from the different electrical boundary conditions at the cortical surface is strongest. The fitted model parameters for the temporal coupling kernel are τ = 21.2 ms and Δ = 1.2 ms.

### Population CSD

The interpretation of the predicted LFP spatial profiles *L*_{n}(*z*_{i}) (shown in Figs. 8*B* and 10*B*) is less obvious than for the MUA spatial profiles *M*_{n}(*z*_{i}). These LFP profiles represent the spatial form of the LFP after firing of action potentials of neurons in population *n* and will consist of a sum over contributions to the LFP from the CSD generated by the synaptic projections from population *n* onto all other populations. By analogy with previous CSD analyses, we can thus separately estimate the CSD resulting from activity in each population, i.e., the population CSD denoted *C*_{n}(*z*_{i}). This measure differs from the traditional estimation of CSD that gives the total CSD across all populations.

The estimation of the population CSD *C*_{n}(*z*_{i}) from the LFP profiles *L*_{n}(*z*_{i}) raises the same methodological issues as the estimation of the total CSD from laminar LFP recordings. In general, the LFP stemming from a particular CSD distribution will depend on both the spatial extent of the CSD and the extracellular conductivity (Pettersen et al. 2006). Here we use the iCSD method (Pettersen et al. 2006) to incorporate these effects in the CSD estimation. The effective lateral sizes of the CSDs [*C*_{n}(*z*_{i})] underlying the LFP profiles [*L*_{n}(*z*_{i})] are not known, however. We therefore consider two different situations: *1*) infinite lateral extent of the CSD (as assumed in the standard CSD method, see Pettersen et al. 2006) and *2*) cylindrical CSD distribution of diameter 0.5 mm centered at the axis of the laminar electrode, with constant CSD in the horizontal direction. In addition to the common choice of assuming a constant electrical conductivity σ, we also consider other electrical boundary conditions at the cortical surface. For *experiment 1* with highly conductive saline added at the cortical surface, we also evaluate CSD estimates for the situation where σ/σ_{0} = 0, i.e., the conductivity σ_{0} above the cortical surface (*z* < 0) is infinitely large. For *experiment 2* with poorly conducting oil added at the cortical surface, we also evaluate CSD estimates for the situation where σ_{0} = 0, i.e., the conductivity above the cortical surface is zero. Although none of these situations can be completely realistic, qualitative features of the population CSD [*C*_{n}(*z*_{i})] present in all sets of estimates are more likely to be genuine than features present only for a particular choice of CSD diameter and/or electrical boundary condition.

In Fig. 11, we show the resulting estimates for the population CSDs for experiments 1 and 2. A first feature to note is that effects of different electrical boundary conditions in the estimation procedure are mainly seen at the electrode contacts positioned 0.1 and 0.2 mm below the cortical surface. The effects of lateral spatial confinement of the CSD sources are seen at all electrode positions, however.

Experiments 1 and 2 were done under different conditions (anesthesia, electrical boundary conditions, stimuli). Nevertheless, a strong qualitative similarity is observed between the predicted population CSD profiles for the three topmost populations (black, blue, red). With the assumption of infinite lateral extent of the CSD (as in the standard CSD method), the synaptic projections from population 2 (blue) are seen to generate a dominant sink centered at ∼0.5–0.6 mm depth, with a dominant source above. For population 3 (red), this pattern is reversed. In addition, population 3 seems to generate a sink at ∼1 mm. For population 1 (black), both experiments predict a dominant sink at ∼0.3 mm, but the source below is predicted to extend deeper into cortex in *experiment 2*. Furthermore, in *experiment 2*, a source is predicted at the top electrode position (0.1 mm), whereas for *experiment 1*, such a source is only seen if a constant conductivity is assumed (solid black). The CSD profiles for the deepest population, population 4 (green), are seen to be quite different in the two experiments, in particular in the upper 1 mm of cortex, but this population also accounts for a rather small part of the total LFP, especially for *experiment 2*.

Most of the same qualitative features of the population CSDs are observed when a population diameter of *d* = 0.5 mm is assumed. Even though the sinks and sources generally are predicted to be somewhat wider and smoother compared with the infinite-diameter assumption, the main predictions described above seem to be robust.

In Fig. 12 *B*, we show the population CSD profiles, assuming infinite population diameters and homogeneous electrical conductivity, for the three topmost populations for all our five data sets from four rats. As observed, the salient features of these population CSDs seem to be quite robust across the experiments, in particular for the two topmost populations. We note, however, that for two of the experiments (3 and 4), the population CSD profiles seem to have been shifted toward lower electrode contacts. Because the same is observed for the MUA profiles in Fig. 12*A*, it seems that the laminar-electrode position was shifted in the vertical direction (relative to the cortical layers) compared with the other experiments.

### Estimation of synaptic connection pattern from LFP template-fitting analysis

The population LFPs [*L*_{n}(*z*); and CSDs, *C*_{n}(*z*)] estimated in the LPA procedure will reflect the LFP (CSD) generated by synaptic inputs from a particular presynaptic population, but it does not specify the postsynaptic neuronal populations receiving these synaptic projections. Such information can be obtained by taking advantage of the fact that the forward electrostatic solution, i.e., the LFP generated by synaptic activation onto a population of similar neurons, can be evaluated. We here apply an LFP template-fitting analysis (see methods and appendix 2) to estimate the synaptic connection pattern between the neuronal populations for our two sample data sets, experiments 1 and 2.

##### FORWARD CALCULATION OF CSD PROFILES.

In the *left columns* of Fig. 13, *A* and *B*, we show the steady-state CSD profiles of a layer 5 and a layer 3 cortical neuron (morphology taken from Mainen and Sejnowski 1996) after injection of a sink-like steady-state current (corresponding to excitatory inputs) into the dendritic trees at different vertical positions measured from the soma. Current has been injected into all dendritic branches crossing a horizontal plane at the particular vertical position, and the total injected current reflects the number and circumference of the intersections with the dendritic branches. From Fig. 13, we see that for the layer 5 cell, excitatory input onto the lower basal dendrites gives a deep sink accompanied by a spatially more distributed source from the return current around the soma and at the top of the apical dendrite. Conversely, excitatory input at the apical dendrite gives a high sink accompanied by a spatially distributed return-current source with a notable contribution from the soma area. A similar, but spatially more compressed pattern is seen in the CSD profile for the layer 3 neuron.

##### FORWARD CALCULATION OF LFP POPULATION PROFILES.

The corresponding steady-state LFP patterns at the center axis of cylindrical populations (with a diameter of 1 mm) of the layer 5 and layer 3 neurons are shown in the *right columns* of Fig. 13, *A* and *B*. In these examples, the CSD (taken from the calculated CSD profiles for the individual neurons) has been assumed to be homogeneously distributed in the horizontal directions inside the cylindrical column. The population LFP profiles are seen to have simple dipolar patterns: For the layer 5 cell, excitatory input at the basal dendrites is seen to give a rather standardized dipolar LFP profile with a low negativity, whereas excitatory input at the apical dendrite gives a similar dipolar with the signs reversed. The same pattern, although with a smaller dipole, is seen for the layer 3 cell. This implies that synaptic activation of populations of pyramidal neurons yields LFPs of limited variability: they can be simply described by a dipole (of either sign) whose spatial size is essentially determined by the vertical extent of the neuron. Note that while a cylindrical column of diameter 1 mm was assumed in the example in Fig. 13, similar dipolar patterns were found for a wide range of population diameters. Inhibitory synaptic inputs correspond to current sources instead of current sinks, but the dipolar patterns above prevail with reversed signs.

##### ESTIMATION OF LFP POPULATION TEMPLATES.

The simple but characteristic dipolar LFP patterns observed in Fig. 13, *A* and *B,* for synaptic activation for these two example pyramidal neurons suggests that the population LFPs [*L*_{n}(*z*)] can be expanded as a weighted sum over such LFP patterns characteristic for the neuronal populations in barrel cortex. Furthermore, the standardized dipolar pattern forms suggest that the LFP from such synaptic activation can be well described by the first PC (found in a singular-value decomposition) of the LFP patterns in Fig. 13.

As detailed in the Discussion, the LPA results suggests the following identification of populations: population 1 (black in Fig. 8*A*) corresponds to a population of layer 2/3 pyramidal cells, population 2 (blue) corresponds to a population of predominantly stellate layer 4 cells, population 3 (red) corresponds to a population of layer 5 pyramidal cells, and population 4 (green) correspond to a population of layer 6 (or deep layer 5) pyramidal cells. Contributions from projections onto the (predominantly) stellate cells (population 2) are assumed negligible because of the closed-field structure of the stellate cells in this layer. As seen in Figs. 2 and 3, the electrical boundary conditions significantly affect the LFP at the topmost electrode and is thus included in the forward LFP modeling.

As described in appendix 2, the characteristic spatial LFP patterns, called LFP templates, for synaptic activation of each of the three pyramidal populations (L2/3, L5, and L6) were obtained in a numerical optimization procedure identifying the population diameters and electrical conductivity-jump parameter *W*_{im} that gives the templates providing the best fit to the total LFP. The soma positions *z*_{n0} were kept fixed at the values given by the MUA fit, whereas the population diameters and the electrical conductivity jump at the cortical surface were allowed to vary freely. The LFP profiles resulting from minimizing the value of *e*_{L} for the two experiments are shown in Fig. 13*C*. For *experiment 1*, the temporal sum over the three optimized LFP templates was found to account for 98.6% of the LFP data, whereas for *experiment 2*, 96.8% of the LFP data were accounted for. Thus our LFP templates are able to account for most of the LFP data.

##### TEMPLATE-FITTING ANALYSIS.

As described in appendix 2 the population LFPs [*L*_{n}(*z*)] can now be projected onto the three fitted LFP templates, produced by activation of the different cell populations: the magnitudes of the fitted weight coefficients measure the proportion of the population LFPs caused by projections from particular presynaptic and onto particular postsynaptic populations, whereas the sign of the fitted weight coefficients tells about the spatial positioning of the synaptic projections (i.e., basal or apical). For example, for both sample experiments, the layer 4 population is found to have a significant positive weight of the template corresponding to the layer 2/3 population (seen in Fig. 13*C*). This means that the layer 4 population is predicted to have a significant excitatory synaptic projection onto the lower basal dendrites of the layer 2/3 population. Likewise, the analysis predicts the layer 4 population to have a weaker excitatory synaptic projection onto the basal dendrites of the layer 5 population (alternatively, an inhibitory projection onto the apical dendrites of the same population). If we assume the dominant synaptic projections to be excitatory, some further predictions from the analysis of our two experiments are as follows: the layer 2/3 population projects to its own basal dendrites, to the apical dendrites of the layer 5 population, and (weakly) to the apical dendrites of the layer 6 population. The layer 5 population projects to its own basal dendrites and to the apical dendrites of the layer 2/3 population. The predicted projections from the layer 6 population are, however, found be conflicting in the two experiments.

## DISCUSSION

### Laminar population identification

Inspection of the predicted spatial and temporal MUA profiles (Figs. 8 and 10) and population CSDs in Fig. 11 suggests the following interpretation: population 1 (black) positioned above ∼0.5 mm cortical depth corresponds to a population of layer 2/3 pyramidal cells, population 2 (blue) centered at ∼0.7 mm corresponds to a population of predominantly stellate layer 4 cells, population 3 (red) centered at ∼1.1–1.2 mm corresponds to a population of layer 5 pyramidal cells, and population 4 (green) positioned below ∼1.5 mm corresponds to a deep infragranular population of layer 6 (or deep layer 5) pyramidal cells. In the following, we will for clarity use this suggested identification when we denote the populations. This population identification is supported by several observations:

##### CORTICAL DEPTH POSITIONS OF ESTIMATED POPULATION BORDERS.

The depth positions extracted by LPA are compatible with previously suggested borders between cortical layers in barrel cortex of similarly sized rats (Armstrong-James et al. 1992; Di et al. 1990). In Armstrong-James et al. (1992), the border between layers 3 and 4 were estimated to be around 0.5 mm, the border between layers 4 and 5 around 0.75–0.8 mm, and the border between layers 5 and 6 at 1.5 mm.

##### TEMPORAL ORDERING OF POPULATION FIRING.

Armstrong-James et al. (1992) measured the response latency after single whisker flicks. In the principal barrel, they found the following approximate sequence of first firing of action potentials: *1*) layer 4, *2*) layer 5b, *3*) layer 3 (slightly later), *4*) layer 5a and 2, and *5*) layer 6. The temporal resolution of our predicted firing rates *r*_{n}(*t*) is lower, but for *experiment 1* (Fig. 8*B*), it is clearly seen that the layer 4 population initiates firing before the layer 2/3 and layer 5 populations. For *experiment 2*, a clear temporal ordering of the initiation of firing cannot be discerned, but the dominance of firing of the layer 4 population for the weakest stimuli (in particular *a1*) is compatible with the established notion of layer 4 being the dominant sensory input layer acting as the main gateway to cortex (Pinto et al. 2003; Swadlow et al. 2002).

##### POPULATION CSD AND ESTIMATED SYNAPTIC CONNECTION PATTERN.

The suggested population identifications are also supported by the observed CSD population profiles and the estimated synaptic connection pattern from the LFP template-fitting analysis: the layer 4 population CSD profiles shown in Figs. 11 and 12 are seen to have large sinks centered at ∼0.5–0.6 mm cortical depth with a large source above, and in the LFP template-fitting analysis the population was predicted to have a significant excitatory synaptic projection onto the lower basal dendrites of the layer 2/3 population. This is in accordance with the known axonal projections from the excitatory layer 4 barrel neurons. Lübke et al. (2000) and Petersen and Sakmann (2001) observed these neurons to have significant axonal projections to pyramidal neurons in layer 2/3 of the same barrel column. Furthermore, the excitatory contacts onto the layer 2/3 neurons have been found to predominantly be positioned at their basal dendrites (Feldmeyer et al. 2002; Lübke et al. 2003). Such a connection pattern would imply a sink around the position of the basal dendrites of the layer 2/3 neurons with a source above caused by the return currents, in agreement with the abovementioned observations for the layer 4 population CSD profile.

The intralaminar connections to other layer 4 neurons should in principle also affect the CSD profile of the layer 4 population. However, most of the layer 4 cells appear to have a stellate shape and thus have a closed-field dendritic structure (Feldmeyer and Sakmann 2000; Feldmeyer et al. 2002; Lübke et al. 2000; Petersen and Sakmann 2001). The net CSD may thus be small, and consequently, the effect of the projection on the observed LFP was neglected in the LFP template-fitting scheme.

The CSD profile for the output of the layer 2/3 population is typically seen in Figs. 11 and 12 to have a dominant sink positioned higher up in cortex (centered at ∼0.3–0.4 mm depth) compared with the CSD profile for the layer 4 population. The LFP template-fitting analysis suggests that this sink is at least partially caused by excitatory projections predominantly onto the basal dendrites of other layer 2/3 neurons, a suggestion supported by the observed strong overlap between the dendritic trees and axonal projections of these neurons (Petersen and Sakmann 2001). This interpretation of the dominant sink implies a superficial source reflecting the return currents in the apical dendrites. Such a source is observed for *experiment 2* in Fig. 11 both for small (*d* = 0.5 mm) and (infinitely) large CSD diameters. In *experiment 1*, such a superficial source is only observed for large CSD diameters when homogeneous electrical boundary conditions (solid black curve in Fig. 11, *top left*) are assumed.

While our LPA approach also provides other predictions for the effective synaptic connection between the identified laminar populations, there is currently little well-established knowledge available for comparison for other parts of the population microcircuit (although the existence of an ubiquity of intracortical connections have been established in, e.g., dual intracellular recordings of single neurons, see Thomson and Bannister 2003). We will discuss these predictions next.

### Model predictions

##### SYNAPTIC CONNECTIONS BETWEEN POPULATIONS.

In addition to the predicted projection from layer 2/3 neurons onto other layer 2/3 neurons, the layer 2/3 pyramidal neurons also project onto layer 5 pyramidal neurons (Thomson et al. 2002). The deep sources (positioned ∼0.5–1 mm below the dominant sink) typically seen in the layer 2/3 CSD profiles in Figs. 11 and 12 may thus stem from a return current from the basal dendrites of layer 5 neurons after excitatory input on their apical dendrites. (The sink from this excitatory input would then contribute to the dominant sink above.) Such a synaptic connection is also predicted by the LFP template-matching analysis.

The CSD profiles from the layer 5 population typically has a dominant source between ∼0.3 and ∼0.8 mm depth with sinks above and below. Layer 5 pyramidal axons arborize most densely within layer 5, but can also send projections to all other layers (Thomson and Bannister 2003). Synaptic contacts onto the basal dendrites of layer 5 neurons from other layer 5 neurons have been shown (Feldmeyer and Sakmann 2000; Thomson and Bannister 2003). This suggests the interpretation that the sinks seen between ∼0.9 and ∼1.5 mm can have contributions from the projections from the population of layer 5 neurons onto the basal dendrites of other layer 5 neurons. This would be accompanied by a more superficial source in the CSD profile, as seen in Figs. 11 and 12. This synaptic projection is also predicted by the LFP template-matching analysis. Furthermore, this scheme predicts an excitatory synaptic projection from the layer 5 population onto the tufted, apical dendrites of the layer 2/3 population (alternatively an inhibitory projection onto the basal dendrites, or both apical excitation and basal inhibition).

The LFP and CSD profiles from the layer 6 population are quite different in the experiments. Furthermore, this population generally accounts for a minor part of the LFP. An attempt at an interpretation of this profile thus does not seem warranted.

##### POPULATION FIRING RATES.

By visual inspection of Fig. 8*C* showing the estimated temporal population firing rates for *experiment 1*, we see that the layer 4 population tends to have the earliest firing, whereas the layer 5 population tends to terminate last. We assessed the temporal ordering of this firing more quantitatively by evaluating the time delay corresponding to the maximum peak of the cross-correlation function (Dayan and Abbott 2001) between the population firing rates. These cross-correlation functions were evaluated over the entire 100-ms window and included all stimulus conditions. For *experiment 1*, these cross-correlation evaluations showed that the layer 4 population firing on average preceded the layer 2/3 population firing by 1 ms and the layer 5 and layer 6 population firing by ∼2.5 ms. For *experiment 2*, a clear temporal ordering is hard to discern by visual inspection of Fig. 10*C*, and the cross-correlation between the layer 4 and layer 2/3 population firing indicated on average no time delay. However, the layer 5 population was found to fire on average 2 ms after the two topmost populations.

The notion of a sequential activation of cortical populations after thalamic input (as hinted in particular by the temporal order of initiation of firing for *experiment 1*) is supported by the stimulus dependence of the observed firing rates in *experiment 2*. As seen in Fig. 10*C*, the weakest stimulus (*a1*) mainly activates the layer 4 population, suggesting the interpretation that the population firing in this case is too weak to inflict a significant response in the supragranular and infragranular populations, i.e., “ignite” the full cortical microcircuit.

### Model assessment

##### DEPENDENCE ON NUMBER OF POPULATIONS.

In the present application of LPA, we assumed four cortical population, i.e., *N*_{pop} = 4. With one less population, i.e., *N*_{pop} = 3, for the data of *experiment 1*, the LPA optimization routine essentially merged the layer 4 and layer 5 populations in Fig. 8*A* into a single population with a concomitant significant increase in the MUA fitting error *e*_{M} (from 0.055 to ∼0.083). With one more population (*N*_{pop} = 5), a much smaller decrease in *e*_{M} was obtained (from 0.055 to ∼0.050). In this case, the layer 4 population in Fig. 8*A* was split into two “subpopulations” with the deeper population generally firing earlier than the more superficial. The subsequent fitting of the LFP also supported the choice *N*_{pop} = 4: for four and five populations, the minimum error *e*_{L} was found to be ∼0.10 and fitted parameter values for the time constant τ and the delay Δ were found to be very similar. For *N*_{pop} = 3, however, *e*_{L} was found to be somewhat higher (0.107), and both fitted parameter values were found to be significantly smaller (τ ∼ 11.1 ms, Δ ∼ 0.8 ms) than for *N*_{pop} = 4 and 5. For *experiment 2*, the choice *N*_{pop} = 3 essentially eliminated the layer 6 population seen in Fig. 10*A*. This was accompanied by an increase of *e*_{M} from 0.066 to ∼0.091 and *e*_{L} from 0.087 to ∼0.095. As for *experiment 1*, increasing the number of populations to *N*_{pop} = 5 gave a smaller decrease of *e*_{M}, from 0.066 to ∼0.058, and essentially split the layer 4 population in Fig. 10*A* into two subpopulations. However, *e*_{L} decreased from 0.087 to ∼0.081. While *N*_{pop} = 4 and 5 gave about the same fitted values for τ (∼21.2–21.7 ms) and Δ (∼1.1–1.2 ms), *N*_{pop} = 3 gave somewhat lower values (τ ∼ 18.9 ms, Δ ∼ 0.9 ms).

##### LFP TEMPLATE-FITTING ANALYSIS.

The LFP template-fitting analysis gives more precise predictions about the synaptic connections in the population circuit than a CSD analysis, but the method relies on more assumptions about the underlying neurons. In particular, the population templates must be constructed based on sample neurons and assumptions regarding the variation of the dendritic structure in a population, parameter values of the membrane and intracellular resistances, and lateral and vertical distributions of soma positions. However, as seen in Fig. 13, the population LFP profiles for the layer 3 and layer 5 pyramidal neurons are found to be very similar: they both have a characteristic dipolar form whose size is essentially determined by the vertical extent of the dendritic trees. The calculations presented above assumed that input current density on the membrane was constant across the dendritic tree. This assumption was probed by calculating LFP population templates assuming the input current at each crossed dendritic branch to be equal (and not proportional to the circumference of the intersection with the horizontal plane). Although some minor quantitative changes in the CSD profiles were found, the LFP templates and consequently the predicted synaptic projection pattern were found to be essentially unaltered. We thus do not expect a sensitive dependence of the LFP profiles on dendritic-structure details or the synaptic weight distribution: given a pyramidal structure, the main determining factor seems to be the vertical extent of the dendrites.

The LFP template-fitting analysis involves more assumptions than the traditional CSD analysis, and the predictions might a priori seem more uncertain. However, the qualitatively similar predictions obtained in the two sample experiments for the synaptic connections between three upper laminar populations (layer 2/3, layer 4, layer 5) are encouraging. We thus believe that the approach may be a useful alternative in future analysis of LFP data and not only in the context of LPA.

### Implications of work

LPA allows for identification of laminar populations, estimation of their firing rates, and population CSDs. In combination with LFP template-fitting analysis, these population CSDs can be used to estimate the functional synaptic coupling patterns between the neuronal populations. Such patterns have previously only been measured by dual or triple whole cell patch-clamp recordings (Thomson et al. 2002; Thomson and Bannister 2003), sometimes in combination with caged glutamate photolysis (Kötter et al. 2005; Schubert et al. 2001). These studies are typically done on rat or cat in in vitro preparations. However, laminar-electrode recordings have been done in awake macaques (Schroeder et al. 1998) and even humans (Halgren et al. 2006; Ulbert et al. 2001, 2004; Wang et al. 2005). Thus LPA, combined with template-fitting analysis, offers opportunities for deeper insights into the functional synaptic connections even for these systems.

## APPENDIX 1: NUMERICAL PROCEDURE IN LAMINAR POPULATION ANALYSIS

In the numerical optimization of the LPA, a matrix formalism is convenient, and in this formalism, *Eqs. 2* and *3* can be reformulated as (7) (8)

Here **φ**_{M}^{m} and **φ**_{L}^{m} are *N*_{ch} × *N*_{t} matrices where *N*_{ch} is the number of laminar electrode channels, and *N*_{t} is the number of data time-points for each stimulus condition (which in our case is 201) multiplied by the number of stimulus conditions. For experiments 1, 3, 4, and 1b, we have *N*_{t} = 27 × 201 = 5,427, whereas for *experiment 2*, we have *N*_{t} = 9 × 201 = 1,809. Furthermore, **M**_{n} and **L**_{n} are *N*_{ch}-dimensional column vectors, and the **r**_{n} and **R**_{n} are *N*_{t}-dimensional row vectors. The row vectors **R**_{n} represent the temporal convolution of the firing rates (*h* ⊗ *r*_{n})(*t*_{j}) where the temporal coupling kernel is given by *Eq. 5*. The convolution with the delayed exponential kemel was approximated using the FILTER command in MATLAB.

These equations can be written in a more compact form as (9) (10) where **M** and **L** are *N*_{ch} × *N*_{pop} matrices given by **M** = (**M**_{1} **M**_{2} … **M**_{Npop}) and **L** = (**L**_{1} **L**_{2} … **L**_{Npop}), and **r** and **R** are *N*_{t} × *N*_{pop} matrices given by **r** = (**r**_{1} **r**_{2} … **r**_{Npop}) and **R** = (**R**_{1} **R**_{2} … **R**_{Npop}).

The matrix **M** describing the MUA spatial profiles is fully specified by the 3 × *N*_{pop} parameters *z*_{0n}, *a*_{n}, and *b*_{n} (*n* = 1, 2, …, *N*_{pop}). Given these parameters, we can thus find the best model estimate for the firing rates directly by performing a pseudoinverse (Gershenfeld 1999) (11) where the dagger denotes pseudoinverse. The estimated MUA data are given by (12) and the relative mean square error *e*_{M} between the best model estimate **φ**_{M,est}^{m} and the data **φ**_{M} can be evaluated from *Eq. 6.*

In the first numerical optimization procedure (fitting the MUA data), the parameters *z*_{0n}, *a*_{n}, and *b*_{n} (*n* = 1, 2, …, *N*_{pop}) are given initial values, and at each cycle, given random increments or decrements. If the new parameter values give a lower error *e*_{M}, these parameters are kept and used as starting point for the next cycle. Typically, this optimization runs for a few thousand cycles. In addition, the optimization routine is run numerous times with different starting values for *z*_{0n}, *a*_{n}, and *b*_{n} to reduce the risk for probing only a local minimum in parameter space. The parameter values giving the overall lowest error *e*_{M} are chosen.

If **r**_{est} and the parameters τ and Δ (describing the temporal coupling kernel) are given, the corresponding value **R**_{est} can be directly evaluated. A model estimate for the spatial LFP profiles **L** can also be found from a pseudoinverse (13) and the corresponding model estimate for the LFP data are given by (14) The relative mean square error *e*_{L} between the best model estimate **φ**_{L,est}^{m} and the data **φ**_{L} can now be evaluated from an expression analogous to *Eq. 6*.

In the second numerical optimization procedure (fitting the LFP data), the parameters τ and Δ are given initial values, and a similar procedure as described above for the MUA fit is used to find the parameter values giving the lowest error *e*_{L}.

## APPENDIX 2: LFP TEMPLATES AND TEMPLATE-FITTING ANALYSIS

The procedure for obtaining LFP templates involves three steps: *1*) evaluation of net transmembrane currents as a function of vertical distance from soma (CSD profiles) accompanying synaptic activation of individual pyramidal neurons, *2*) construction of corresponding LFP profiles for a laminar population of similar pyramidal neurons receiving the same synaptic inputs, and *3*) determination of appropriate LFP population templates for a particular laminar-electrode experiment by fitting LFP profiles to the LFP data. The population LFPs were projected onto these fitted templates to estimate the weight and sign of synaptic projection between the populations.

### Evaluation of CSD profiles from synaptic activation of pyramidal neurons

Based on the assumption that the observed LFP is mainly caused by synaptic activation of pyramidal neurons, we constructed CSD profiles for the reconstructed layer 3 and layer 5 pyramidal neurons used and made by publicly available by Mainen and Sejnowski (1996). These were reconstructed neurons from cat visual cortex and not rat barrel cortex. However, we found that our results depended only on the gross features of the dendritic morphology (mainly the vertical extent), and we would not expect qualitatively different results if barrel neurons were used instead. The scripts for the simulations of Mainen and Sejnowski, based on the simulator tool NEURON (http://neuron.duke.edu/), were downloaded from http://www.cnl.salk.edu/∼zach/methods.html. All active conductances, as well as the axon, were removed from their model neurons, making the models fully linear. For our application, we needed the steady-state CSD profile, i.e., the profile after the initial transient after onset of synaptic current, and the resistance parameters determining this profile were chosen to be the same as in Mainen and Sejnowski (1996), i.e., *r*_{m} = 30,000 Ω cm^{2} (membrane resistance) and *r*_{a} = 150 Ω cm (axial resistance). The model CSD profiles were obtained in NEURON by injecting a constant current into each dendritic branch crossing an imaginary horizontal plane positioned at a particular height, *z*_{s}, above the soma. The magnitude of the injected current at each branch was chosen to be proportional to the circumference of the intersection (in the horizontal plane) of the dendritic branch. The steady-state transmembrane currents in the soma and all dendritic compartments were recorded. The procedure was repeated for a set of equidistant vertical positions for the input-current planes (with a distance of 0.01 mm) covering the entire vertical extent of the neuron's dendritic structure. From these recorded transmembrane currents and the known three-dimensional (3D) structure of the neuron, the cortical depth variation (*z* variation) of the transmembrane current, *c*_{k}(*z,z*_{s}), i.e., CSD profile, was calculated. This transmembrane current depends both on the vertical position of the synaptic input (*z*_{s}) and the type of neuron (*k*). To remove possible artifacts in *c*_{k}(*z,z*_{s}) caused by the compartmentalization of the neuron, we incorporated a Gaussian distribution of soma depth position with a SD of 0.02 mm. The current representing the synaptic input itself was in the final CSD profile assumed to be Gaussian distributed around the center position with a SD of 0.071 mm.

Example results of *c*_{k}(*z,z*_{s}) for a layer 5 and a layer 3 pyramidal neuron are shown in left column of Fig. 13, *A* and *B*. Note that the CSD profiles depict the response to positive synaptic current onto the dendritic branches, corresponding to excitatory synaptic inputs. Because of the linearity of the forward model, the response to negative synaptic current (corresponding to inhibitory synaptic inputs) would have exactly the same CSD profile, but with opposite sign.

### Evaluation of LFP profiles from synaptic activation of pyramidal cell populations

From these CSD profiles, we calculated the corresponding LFP profiles at the center axis of a cylindrical population of diameter *d* of similar neurons assuming in-plane homogeneity of the CSD inside the discs (of diameter *d*). For a homogeneous extracellular medium with conductivity σ, the potential φ(*z*) at the center axis at position *z* from a homogeneous CSD disc positioned at the position *z′* is given by (Nicholson and Llinas 1971; Pettersen et al. 2006) (15) where *c*_{p} is the planar CSD. The electrical boundary conditions inferred on the electrostatic problem caused by a jump in the conductivity at the cortical surface were included by means of the method of images. With the cortical surface defined to be at *z* = 0 with σ = σ_{0} above cortex (*z* < 0), the expression in *Eq. 15* generalizes to (Nicholson and Llinas 1971) (16) where the parameter *W*_{im} = (σ/σ_{0} − 1)/(σ/σ_{0} + 1) gives the weight of the image current-source potential. For a perfectly conducting media above cortex (σ_{0} → ∞) *W*_{im} = −1, whereas for a perfectly insulating media (σ_{0} → 0) *W*_{im} = 1.

Given the model CSD profile *c*_{k}(*z,z*_{syn}), the depth position of the soma (measured from the cortical surface) of the neuron *z*_{k0}, the diameters of the population activities d_{k}, and the ratio between the extracellular conductivities inside and above cortex (*W*_{im}), the position dependence of the LFP profile at the center axis *l*_{k}(*z,z*_{syn};*z*_{k0}) can be uniquely calculated following the above recipe. Example results for our layer 5 and layer 3 pyramidal neurons are shown in the *left columns* of Fig. 13, *A* and *B*, respectively.

### Determination of LFP templates from experimental LFP data

Example results of the LFP profile at the center axis *l*_{k}(*z,z*_{syn};*z*_{k0}), shown in the *right columns* of Fig. 13, *A* and *B*, show a quite standardized “dipolar” form of the LFP profiles, regardless of the position of the synaptic input. The main difference between a synaptic input onto the tufted apical dendrites and the proximal, basal dendrites is the strength and sign of this “dipolar” LFP profile. In the following analysis, we thus represent *l*_{k}(*z,z*_{syn};*z*_{k0}) by its first PC *l*_{k}^{PC1}(*z*;*z*_{k0}; containing the dipolar contribution) obtained by a singular-value decomposition. For the example LFP profiles shown in Fig. 13 for our layer 5 and layer 3 populations, these first PCs *l*_{k}^{PC1} are found to account for 79% and 86%, respectively, of the total LFP profiles *l*_{k}.

We now model the LFP data, φ_{L}(*z,t*), by a weighted temporal sum over these first principal components (17) To incorporate the effect of differences in the detailed dendritic structure of the pyramidal cells constituting the population and the spatial spread of their soma positions within the layer, we spatially filtered the resulting PC *l*_{k}^{PC1}(*z*;*z*_{k0}) with a symmetric, roughly triangular filter with a half-width of ∼0.3 mm.

We included three postsynaptic populations: a layer 2/3 population based on the layer-3 neuron morphology in Fig. 13, a layer 5 population based on the layer-5 neuron morphology in Fig. 13, and a layer 6 population also based on this layer-5 neuron morphology.

The spatial form of the LFP profiles and their first PCs will depend on the image current-source weight *W*_{im}, the population diameters *d*_{k}, and the depth positions of the neurons in each population (*z*_{k0}). In the fitting of the model expression for φ_{L}^{m}(*z,t*) in *Eq. 17* to experimental data φ_{L}(*z,t*), we kept *z*_{k0} fixed at values determined from the MUA fit in LPA for the layer 2/3, layer 5, and layer 6 populations. *W*_{im} and the population diameters *d*_{k} were varied to minimize the deviation between the model and the experimental data.

For the numerical optimization, a matrix formalism is convenient here, and in this formalism, *Eq. 17* can be written as (18) where **φ**_{L}^{m} is a *N*_{ch}×*N*_{t} matrix, where *N*_{t} is the number of stimulus conditions multiplied with the number of data time points for each stimulus condition (which in our case is 201). Furthermore, **l**_{k} is an *N*_{ch}-dimensional column vector, and **a**_{k} is an *N*_{t}-dimensional row vector. *Equation 18* can be written in a more compact form as (19) Here, **l**^{PC1} is an *N*_{ch}×*N*_{k} matrix given by **l**^{PC1}=(**l**_{1}^{PC1} **l**_{2}^{PC1} … **l**_{Nk}^{PC1}), and **a** is an *N*_{t} × *N*_{k} matrix given by **a** = (**a**_{1} **a**_{2} … **a**_{Nk}).

The matrix **l**^{PC1} describing the spatial profiles of the *N*_{k} postsynaptic neuron populations are fully specified by the calculated CSD depth profile *c*_{k}(*z,z*_{syn}), the population center position *z*_{k0}, the population diameter *d*_{k}, and *W*_{im}. Given these parameters, we can calculate the matrix **l**^{PC1} and find the best model estimate for the time vector matrix **a** directly by performing a pseudoinverse (Gershenfeld 1999) (20) The estimated LFP data are given by (21) and the relative mean square error *e*_{L} between the best model estimate **φ**_{L,est}^{m} and the data **φ**_{L} can be evaluated by (22)

In the numerical optimization procedure, the population depths *z*_{0k} was kept fixed, whereas the values of the population diameters *d*_{k} (*n* = 1, 2, …, *N*_{k}) and *W*_{im} were optimized to give a minimum value of *e*_{L} using a similar stochastic search method as described for the LPA procedure in appendix 1.

Example results from LFP template fits are shown in Fig. 13*C*. Note that these templates are derived from injecting positive synaptic currents, corresponding to excitatory synaptic inputs, into the dendritic branches. By injecting negative synaptic current (corresponding to inhibitory synaptic inputs) instead, the LFP templates would have exactly the same form, but with opposite sign.

### Projection of populations LFPs onto LFP templates

The population LFPs, *L*_{n}(*z*_{i}), obtained from LPA were expanded into the fitted LFP templates (23) or in matrix form (24) Estimates for the connection matrix **w**, i.e., the weight of the synaptic projections, could now be found directly from a pseudoinverse (25)

## GRANTS

This work was supported by the National Institute of Health Grants R01 EB-00790, NS-051188, and NS-1874 and by the Research Council of Norway.

## Acknowledgments

We thank M. Kermit for help with the data processing and S. S. Cash and P. Blomquist for a critical reading of the manuscript.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 by the American Physiological Society