| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Department of Mathematical Sciences and Technology and Center for Integrative Genetics, Norwegian University of Life Sciences, Ås, Norway; 2Athinoula A. Martinos Center, Massachusetts General Hospital, Charlestown, Massachusetts; 3Departments of Neurosciences and Radiology, University of California, San Diego, La Jolla, California; and 4Institute for Psychology of the Hungarian Academy of Sciences, Budapest, Hungary
Submitted ; accepted in final form
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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% O2 in air. Ventilation parameters were adjusted to maintain PaCO2 between 35 and 45 mmHg, PaO2 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 (a3,t1). 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.
|
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(zi,tj) and M(zi,tj), respectively. Here zi refers to electrode contact position and tj to time. The time resolution is 0.5 ms.
LPA
Two-dimensional data like our LFP [L(zi,tj)] and MUA [M(zi,tj)] data can be expanded in a series of spatiotemporally separable functions
 | (1) |
In PCA, the expansion functions fn(zi) (called spatial loading) and gn(tj) (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 [fn(zi)] and the temporal scores [gn(tj)] 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) |
Because the MUA data are rectified, we impose the constraint that the spatial functions Mn(zi) should be nonnegative. The interpretation of rn(tj) 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 rn(tj) 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 rn(tj) 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) |
rn)(tj) is the temporal convolution given by
 | (4) |
Mm(zi,tj) and Lm(zi,tj) and the experimental data M(zi,tj) and L(zi,tj) as small as possible.
FORM OF SPATIAL POPULATION PROFILES Ln(zi) and Mn(zi).
The spatial profiles Ln(zi) are allowed to vary freely. However, for the MUA spatial profiles Mn(zi), 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 Mn(zi) as a trapezoidal symmetric profile centered at the position z0n. The profile has a constant height between z0n – an/2 and z0n + an/2 and decays linearly to zero over a length bn on each side. We here constrain the flat maxima of the spatial profiles Mn(zi) to be nonoverlapping, (e.g., z0n + an/2 < z0N – aN/2 if population N is positioned immediately below population n in cortex). We further constrain bn to be smaller than 0.1 mm to enforce a rather sharp interface between neighboring populations.
FORM OF TEMPORAL COUPLING KERNEL hn(tk).
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) |
(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 hn(tk) = h(tk).
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 Mn(z) is described by three parameters (z0n, an, bn). Thus in the numerical optimization scheme 3 x Npop parameters are varied to minimize the deviation eM between the experimental MUA data and the model expression. Here eM is the relative mean square error
 | (6) |
In the second step of the procedure, we fit the LFP data given the fitted firing rates rn(tj) from the MUA fit. The spatial LFP profiles Ln(zi) are fitted nonparametrically, so the parameter fit involves only the two model parameters of h(tj), 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 Cn(zi) from the population LFP profiles Ln(zi) (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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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).
|
|
(in particular their Fig. 4).
|
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).
|
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 hn(tj) 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. 3B 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. 3A 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 gM1(tj) for the MUA and gL1(tj) for the LFP, can be used to obtain a direct estimate of the temporal coupling kernel h(tj). 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 gL1(tj) occurs about 5 ms after the peak in corresponding MUA temporal score gM1(tj). 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 hn(tj) given in Eq. 5. We will therefore use this in the further analysis of the data from all stimuli.
|
|
|
= 13.4 ms and = 2.0 ms.
The fitted LFP spatial profiles Ln(zi), n = 1, 2, 3, 4, are shown in Fig. 8B, and in Fig. 8D we show the time-courses of the postsynaptic LFP activity, (h rn)(tj), for the four populations. These time-courses are obtained by convolving the population firing rates rn(tj) in Fig. 8C 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 rn)(tj) in Fig. 8D with rn(tj) in Fig. 8C thus shows much smoother time-courses. The contribution to the LFP from each population n is determined by (h rn)(tj) multiplied by the corresponding Ln(zi). 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 eM = 0.066 and eL = 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 [Mn(zi) and Ln(zi), n = 1, 2, 3, 4] are shown in Fig. 10, A and B, whereas the fitted population firing rates rn(tj) and the postsynaptic time-courses (h rn)(tj) are shown in Fig. 10, C and D, respectively.
|
|
With the exception of population 4 (green), the predicted LFP spatial profiles Ln(zi) for the two experiments are seen to share most of the qualitative features (cf. Figs. 8B and 10B). 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 Ln(zi) (shown in Figs. 8B and 10B) is less obvious than for the MUA spatial profiles Mn(zi). 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 Cn(zi). This measure differs from the traditional estimation of CSD that gives the total CSD across all populations.
The estimation of the population CSD Cn(zi) from the LFP profiles Ln(zi) 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 [Cn(zi)] underlying the LFP profiles [Ln(zi)] 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 [Cn(zi)] 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.
|
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. 12A, it seems that the laminar-electrode position was shifted in the vertical direction (relative to the cortical layers) compared with the other experiments.
|
The population LFPs [Ln(z); and CSDs, Cn(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.
|
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 [Ln(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. 8A) 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 Wim that gives the templates providing the best fit to the total LFP. The soma positions zn0 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 eL for the two experiments are shown in Fig. 13C. 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 [Ln(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. 13C). 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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 rn(t) is lower, but for experiment 1 (Fig. 8B), 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. 8C 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. 10C, 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. 10C, 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., Npop = 4. With one less population, i.e., Npop = 3, for the data of experiment 1, the LPA optimization routine essentially merged the layer 4 and layer 5 populations in Fig. 8A into a single population with a concomitant significant increase in the MUA fitting error eM (from 0.055 to 0.083). With one more population (Npop = 5), a much smaller decrease in eM was obtained (from 0.055 to 0.050). In this case, the layer 4 population in Fig. 8A 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 Npop = 4: for four and five populations, the minimum error eL was found to be 0.10 and fitted parameter values for the time constant and the delay were found to be very similar. For Npop = 3, however, eL 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 Npop = 4 and 5. For experiment 2, the choice Npop = 3 essentially eliminated the layer 6 population seen in Fig. 10A. This was accompanied by an increase of eM from 0.066 to 0.091 and eL from 0.087 to 0.095. As for experiment 1, increasing the number of populations to Npop = 5 gave a smaller decrease of eM, from 0.066 to 0.058, and essentially split the layer 4 population in Fig. 10A into two subpopulations. However, eL decreased from 0.087 to 0.081. While Npop = 4 and 5 gave about the same fitted values for (21.2–21.7 ms) and (1.1–1.2 ms), Npop = 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
 | (7) |
 | (8) |
Here Mm and Lm are Nch x Nt matrices where Nch is the number of laminar electrode channels, and Nt 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 Nt = 27 x 201 = 5,427, whereas for experiment 2, we have Nt = 9 x 201 = 1,809. Furthermore, Mn and Ln are Nch-dimensional column vectors, and the rn and Rn are Nt-dimensional row vectors. The row vectors Rn represent the temporal convolution of the firing rates (h rn)(tj) 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) |
The matrix M describing the MUA spatial profiles is fully specified by the 3 x Npop parameters z0n, an, and bn (n = 1, 2, ..., Npop). Given these parameters, we can thus find the best model estimate for the firing rates directly by performing a pseudoinverse (Gershenfeld 1999)
 | (11) |
 | (12) |
M,estm and the data M can be evaluated from Eq. 6. In the first numerical optimization procedure (fitting the MUA data), the parameters z0n, an, and bn (n = 1, 2, ..., Npop) are given initial values, and at each cycle, given random increments or decrements. If the new parameter values give a lower error eM, 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 z0n, an, and bn to reduce the risk for probing only a local minimum in parameter space. The parameter values giving the overall lowest error eM are chosen.
If rest and the parameters and (describing the temporal coupling kernel) are given, the corresponding value Rest can be directly evaluated. A model estimate for the spatial LFP profiles L can also be found from a pseudoinverse
 | (13) |
 | (14) |
L,estm 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 eL.
|
APPENDIX 2: LFP TEMPLATES AND TEMPLATE-FITTING ANALYSIS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
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., rm = 30,000 cm2 (membrane resistance) and ra = 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, zs, 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, ck(z,zs), i.e., CSD profile, was calculated. This transmembrane current depends both on the vertical position of the synaptic input (zs) and the type of neuron (k). To remove possible artifacts in ck(z,zs) 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 ck(z,zs) 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) |
= 0 above cortex (z < 0), the expression in Eq. 15 generalizes to (Nicholson and Llinas 1971)
 | (16) |
/0 – 1)/(/0 + 1) gives the weight of the image current-source potential. For a perfectly conducting media above cortex (0 
) Wim = –1, whereas for a perfectly insulating media (0 0) Wim = 1. Given the model CSD profile ck(z,zsyn), the depth position of the soma (measured from the cortical surface) of the neuron zk0, the diameters of the population activities dk, and the ratio between the extracellular conductivities inside and above cortex (Wim), the position dependence of the LFP profile at the center axis lk(z,zsyn;zk0) 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 lk(z,zsyn;zk0), 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 lk(z,zsyn;zk0) by its first PC lkPC1(z;zk0; 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 lkPC1 are found to account for 79% and 86%, respectively, of the total LFP profiles lk.
We now model the LFP data, L(z,t), by a weighted temporal sum over these first principal components
 | (17) |
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 Wim, the population diameters dk, and the depth positions of the neurons in each population (zk0). In the fitting of the model expression for Lm(z,t) in Eq. 17 to experimental data L(z,t), we kept zk0 fixed at values determined from the MUA fit in LPA for the layer 2/3, layer 5, and layer 6 populations. Wim and the population diameters dk 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) |
Lm is a NchxNt matrix, where Nt 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, lk is an Nch-dimensional column vector, and ak is an Nt-dimensional row vector. Equation 18 can be written in a more compact form as
 | (19) |
The matrix lPC1 describing the spatial profiles of the Nk postsynaptic neuron populations are fully specified by the calculated CSD depth profile ck(z,zsyn), the population center position zk0, the population diameter dk, and Wim. Given these parameters, we can calculate the matrix lPC1 and find the best model estimate for the time vector matrix a directly by performing a pseudoinverse (Gershenfeld 1999)
 | (20) |
 | (21) |
L,estm and the data L can be evaluated by
 | (22) |
In the numerical optimization procedure, the population depths z0k was kept fixed, whereas the values of the population diameters dk (n = 1, 2, ..., Nk) and Wim were optimized to give a minimum value of eL 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. 13C. 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, Ln(zi), obtained from LPA were expanded into the fitted LFP templates
 | (23) |
 | (24) |
 | (25) |
|
GRANTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: G. T. Einevoll, Dept. of Mathematical Sciences and Technology, Norwegian Univ. of Life Sciences, PO Box 5003, N-1432 Ås, Norway (E-mail: Gaute.Einevoll{at}umb.no)
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX 1: NUMERICAL PROCEDURE...APPENDIX 2: LFP TEMPLATES...GRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Barth DS, Baumgartner C, Di S. Laminar interactions in rat motor cortex during cyclical excitability changes of the penicillin focus. Brain Res 508: 105–117, 1990.[CrossRef][Web of Science][Medline]
Barth DS, Di S. Laminar excitability cycles in neocortex. J Neurophysiol 65: 891–898, 1991.
Barth DS, Di S, Baumgartner C. Laminar cortical interactions during epileptic spikes studied with principal component analysis and physiological modeling. Brain Res 484: 13–35, 1989.[CrossRef][Web of Science][Medline]
Brecht M, Roth A, Sakmann B. Dynamic receptive fields of reconstructed pyramidal cells in layers 3 and 2 of rat somatosensory cortex. J Physiol 533: 243–265, 2003.
Buzsaki G. Large-scale recording of neuronal ensembles. Nature Neurosci 7: 446–451, 2004.[CrossRef][Web of Science][Medline]
Dale AM. Optimal experimental design of event-related fMRI. Hum Brain Mapp 8: 109–114, 1999.[CrossRef][Web of Science][Medline]
Dayan P, Abbott LF.. Theoretical Neuroscience. Cambridge, MA: MIT Press, 2001.
Di S, Baumgartner C, Barth DS. Laminar analysis of extracellular field potentials in rat vibrissa/barrel cortex, J Neurophysiol 63: 832–840, 1990.
Einevoll GT, Devor A, Ulbert I, Kermit M, Halgren E, Dale AM. Physiologically constrained population modelling of cortical activity in rat barrel cortex measured with laminar electrodes. Soc Neurosci Abstr 921.13, 2004.
Feldmeyer D, Lübke J, Silver RA, Sakmann B. Synaptic connections between layer 4 spiny neurone—layer 2/3 pyramidal cells in juvenile rat barrel cortex: physiology and anatomy of interlaminar signalling within a cortical column, J Physiol 538: 803–822, 2002.
Feldmeyer D, Sakmann B. Synaptic efficacy and reliability of excitatory connections between the principal neurons of the input (layer 4) and output layer (layer 5) of the neocortex. J Physiol 525: 31–39, 2000.
Gershenfeld N. The Nature of Mathematical Modeling. New York: Cambridge, 1999.
Gold C, Henze DA, Koch C, Buzsaki G. On the origin of the extracellular action potential waveform: a modeling study. J Neurophysiol 95: 3113–3128, 2006.
Halgren E, Wang C, Schomer DL, Knake S, Marinkovic K, Wu J, Ulbert I. Processing stages underlying word recognition in the anteroventral temporal lobe. Neuroimage, 30: 1401–1413, 2006.[CrossRef][Web of Science][Medline]
Holt GR, Koch C. Electrical interactions via the extracellular potential near cell bodies. J Comput Neurosci 6: 169–184, 1999.[CrossRef][Web of Science][Medline]
Koch C, Segev I. Methods in Neuronal Modeling. Cambridge, MA: MIT Press, 1998.
Kötter R, Schubert D, Dyhrfjeld-Johnsen J, Luhmann HJ, Staiger JF. Optical release of caged glutamate for stimulation of neurons in the in vitro slice preparation. J Biomed Opt 10: 11003, 2005.[CrossRef][Medline]
Lübke J, Egger V, Sakmann B, Feldmeyer D. Columnar organization of dendrites and axons of single and synaptically coupled excitatory spiny neurons in layer 4 of the rat barrel cortex. J Neurosci 20: 5300–5311, 2000.
Lübke J, Roth A, Feldmeyer D, Sakmann B. Morphometric analysis of the columnar innervation domain of neurons connecting layer 4 and layer 2/3 of juvenile rat barrel cortex. Cereb Cort 13: 1051–1063, 2003.
Mainen ZF, Sejnowski TJ. Influence of dendritic structure on firing pattern in model neocortical neurons. Nature 382: 363–366, 1996.[CrossRef][Medline]
Mitzdorf U. Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena. Physiol Rev 65: 37–100, 1985.
Nadasdy Z, Csicsvari J, Penttonen M, Hetke J, Wise K, Buzsaki G. Extracellular recording and analysis of neuronal activity: from single cells to ensembles. In: Neuronal Ensembles, edited by Eichenbaum H, Davis JL. New York: Wiley, 1998, p. 17–55.
Nicholson C, Freeman JA. Theory of current source-density analysis and determination of conductivity tensor for anuran cerebellum. J Neurophysiol 38: 356–368, 1975.
Nicholson C, Llinas R. Field potentials in the alligator cerebellum and theory of their relationship to Purkinje cell dendritic spikes. J Neurophysiol 34: 509–531, 1971.
Petersen CCH, Sakmann B. Functionally independent columns of rat somatosensory barrel cortex revealed with voltage-sensitive dye imaging. J Neurosci 21: 8435–8446, 2001.
Pettersen KH, Devor A, Ulbert I, Dale AM, Einevoll GT. 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, 2006.[CrossRef][Web of Science][Medline]
Pinto DJ, Hartings JA, Brumberg JC, Simons DJ. Cortical damping: analysis of thalamocortical response transformations in rodent barrel cortex. Cereb Cort 13: 33–44, 2003.
Rappelsberger P, Pockberger H, Petsche H. Current source density analysis: methods and application to simultaneously recorded field potentials of the rabbit's visual cortex. Pfluegers 389: 159–170, 1981.[CrossRef]
Schroeder CE, Mehta AD, Givre SJ. A spatiotemporal profile of visual system activation revealed by current source density analysis in the awake macaque. Cereb Cort 8: 575–592, 1998.
Schubert D, Staiger JF, Cho N, Kötter R, Zilles K, Luhmann HJ. Layer-specific intracolumnar and transcolumnar functional connectivity of layer V pyramidal cells in rat barrel cortex. J Neurosci 21: 3580–3592, 2001.
Somogyvári Z, Zalányi L, Ulbert I, Érdi P. Model-based source localization of extracellular action potentials. J Neurosci Methods 147: 126–137, 2005.[CrossRef][Web of Science][Medline]
Swadlow HA, Gusev AG, Bezdudnaya T. Activation of a cortical column by a thalamocortical impulse. J Neurosci 22: 7766–7773, 2002.
Thomson AM, Bannister AP. Interlaminar connections in the neocortex. Cereb Cort 13: 5–14, 2003.
Thomson AM, West DC, Wang Y, Bannister AP. Synaptic connections and small circuits involving excitatory and inhibitory neurons in layers 2–5 of adult rat and cat neocortex: triple intracellular recordings and biocytin labelling in vitro. Cereb Cort 12: 936–953, 2002.
Ulbert I, Halgren E, Heit G, Karmos G. Multiple microelectrode-recording system for human intracortical applications. J Neurosci Methods. 106: 69–79, 2001.[CrossRef][Web of Science][Medline]
Ulbert I, Heit G, Madsen J, Karmos G, Halgren E. Laminar analysis of human neocortical interictal spike generation and propagation: current source density and multiunit analysis in vivo. Epilepsia 45(4 Suppl): 48–56, 2004.
Wang C, Ulbert I, Ives JR, Schomer DL, Blume H, Marinkovic K, Halgren E. Responses of human anterior cingulate cortex microdomains to error detection, conflict monitoring, stimulus-response mapping, familiarity, and orienting. J Neurosci 25: 604–613, 2005.
Woolsey TA, Van der Loos H. The structural organization of layer IV in the somatosensory region (SI) of mouse cerebral cortex. Brain Res 17: 205–242, 1970.[CrossRef][Web of Science][Medline]
This article has been cited by other articles:
|
|
 |
|
  A. Devor, P. Tian, N. Nishimura, I. C. Teng, E. M. C. Hillman, S. N. Narayanan, I. Ulbert, D. A. Boas, D. Kleinfeld, and A. M. Dale Suppressed Neuronal Activity and Concurrent Arteriolar Vasoconstriction May Explain Negative Blood Oxygenation Level-Dependent Signal J. Neurosci., April 18, 2007; 27(16): 4452 - 4459. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
