a1. Continuous LFP forward model.

The extra-cellular potential ϕ (that is, the LFP) is generated by transmembrane currents that set up an electric field E and induce an associated extra-cellular current density J in the tissue volume [12]. The extra-cellular potential is related to J by (1) where σ denotes the conductivity tensor of the tissue [12]. Given that cortical tissue is predominantly organized in the intra- and inter-laminar directions, it is often assumed that σ is a diagonal matrix, when expressed in Cartesian coordinates (x, y, z), where (x, y) denote the intra-laminar (horizontal) and z denotes the inter-laminar (vertical) location. Let us denote the diagonal entries of σ by (σ_x, σ_y, σ_z). From Eq (1) and using Cartesian coordinates, it follows that (2) where we have introduced the current source density C (3)

It is a scalar quantify with the dimension of current per unit-of-volume. Eq (2) is known as the anisotropic Poisson’s equation. If the tissue is homogeneous, that is, if σ do not depend on location, Eq (2) can be solved by applying the coordinate transformation (4) which converts it into the isotropic Poisson’s equation. The latter can be solved explicitly and, after applying the inverse coordinate transformation, yields (5) where the kernel K_σ is given by (6) [12] and where the integral is taken over the tissue volume .

a2. Discretization of the forward model.

For practical use, both sides of the continuous forward model (Eq (5)) need to be discretized. Discretization of the left-hand-side of amounts to calculating the potential ϕ at the electrode tips. Discretization of the right-hand-side of such an equation (Fredholm integral equation of the first kind) is generally done either by numerical integration or by expansion of C using a set of basis functions [32]. Since K_σ is singular for (u, v, w) = (x, y, z), numerical integration becomes problematic when the electrode tip is located in active tissue (that is, at locations for which C ≠ 0). In the case of inter-laminar LFP recordings, the singularity can be dealt with by assuming C to have the form C = C_hC_v, where C_h and C_v are intra- and inter-laminar components of C, respectively [25, 30]. For a given choice of C_h, which corresponds to choosing an a priori intra-laminar source profile, C_h is integrated out of Eq (5) (either analytically or numerically), yielding a one-dimensional continuous forward model with a non-singular kernel that can be discretized either by numerical integration or basis function expansion. Although in this study we also assume that C can be decomposed into an intra- and an inter-laminar component, it might be advantageous for future studies to have a more general way of discretizing Eq (5), which can be done by expanding C using suitable basis functions. Before choosing the basis functions, we briefly outline this approach.

The source space of the continuous forward model is the infinite-dimensional vector space of square-integrable functions C defined on the tissue volume , denoted by . We now choose n linear independent CSDs and restrict the CSDs to the subspace space spanned by ρ₁, ⋯, ρ_n. For every , there are unique coefficients C₁, ⋯, C_n such that (7) and we can therefore identify every with the vector of its expansion coefficients C = (C₁, ⋯, C_n)^t, where t denotes matrix transpose. The source space of the discretized forward model hence is the n-dimensional Euclidean vector space. By substituting Eq (7) into Eq (5), it follows that (8) which reduces the calculation of ϕ(x, y, z) to calculating the potential (at (x, y, z)) generated by the basis functions ρ₁, ⋯, ρ_n.

To obtain a discrete forward model, the potential field ϕ has to be discretized as well. We thus select p locations (u_i, v_i, w_i) (i = 1, ⋯ p) within the tissue volume , which correspond to the locations of the p electrode tips. Denote the corresponding potentials by Φ₁, ⋯, Φ_p. Note that this reduces the data space to the p-dimensional Euclidean space. Define the leadfield matrix by (9) for i = 1, ⋯, p and j = 1, ⋯, n. This gives the discretized forward model (10) where Φ = (Φ₁, ⋯, Φ_p)^t.

Which basis functions are appropriate to represent C? In previous LFP inverse modeling studies, different basis functions have been used, including step functions, balls, splines, Gaussians, and data kernels [25–27, 30]. Although choosing the data kernels is an interesting option because it yields a low-dimensional model space, the data kernels of the general forward model (Eq (6)) are singular, which prohibits their use as basis functions in the current context. Instead, as basis functions we choose homogeneous voxels, whose potential has been explicitly calculated for the cubic and isotropic case [33]. In S1 Text, we generalize this formula to the case of rectangular and anisotropic voxels. These are the indicator functions of voxels within the tissue volume , scaled to have unit norm. The advantage of rectangular monopoles over other basis functions such as Gaussians or balls is that monopolar basis functions are orthonormal, in which case the relationship between the singular value expansion (SVE) of the continuous forward model and the singular value decomposition (SVD) of the discretized model is well-understood [32]. Moreover, their supports constitute a (arbitrary fine-grained) partition of which allows them to efficiently represent arbitrary current distributions. Another advantage is that only the voxel dimensions have to be chosen, while in the case of balls or Gaussians, besides their dimensions, their locations have to be chosen as well, which introduces unnecessary (free) parameters in the discretization of the forward model.

a3. Utah forward model for cortical evoked responses.

In the simulations we focus on inverse modeling of LFPs recorded with the Utah intracortical electrode array [34], which is one of the most frequently used arrays and comprises 100 electrodes, arranged in a 10 × 10 array with 400 μm inter-electrode spacing. To fully specify the leadfield matrix G, a tissue volume , together with a discrete sampling have to be chosen. We take to have intra-laminar extent 7.2 × 7.2 mm and inter-laminar extent 3.1 mm, which is about the thickness of the macaque neocortex. Fig 1A and 1B provide illustrations. Note that the intra-laminar extent of is twice that of the Utah array (which equals 3.6 mm), which allows simulating currents outside the array, which is generally the most realistic scenario. Next, the source space has to be discretized. As the inverse methods behave essentially different for source spaces that have a higher resolution than the Utah array (that is, whose voxel-length is larger than 400 μm) and those that have a lower resolution than the Utah array (that is, whose voxel-length equals 400 μm), we considered two different discretization schemes, corresponding to voxel dimensions of 100 × 100 × 100 μm³ and 400 × 400 × 100 μm³, respectively. In both cases, the Utah array was placed at a depth of 1 mm at the intra-laminar center of the tissue volume. Fig 1C provides an illustration. The high-resolution source space thus contains n_v = 32 horizontal cortical slices, each containing voxels and the total number of voxels hence equals [calculate this leadfield again]. The low-resolution source space contains n_v = 31 cortical slices, each containing voxels and the total number of voxels hence equals .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Source space for the Utah array.

A. Schematic drawing of a piece of cortex enclosing the rectangular source space (black rectangle) and the Utah electrode array (red line). B. Close-up of the source space and Utah array. The array is located 1 mm under the pial surface. C. Two discretizations of the (intra-laminar) source space (left panel: high-resolution, right panel: low-resolution). Black dots and red circles denote the centers of source voxels and locations of the recording electrodes, respectively. In the high-resolution source space, each cortical slice contains 61 × 61 voxels with intra-laminar lengths of 100 μm. In the low-resolution source space, each cortical slice contains 18 × 18 voxels with intra-laminar lengths of 400 μm. The high and low resolution source spaces comprise 31 cortical slice, each 100 μm thick.

https://doi.org/10.1371/journal.pone.0187490.g001

An evoked CSD is specified by its value on each of the source voxels. It will be convenient to use its vectorization which is obtained by stacking its columns on top of each other. The evoked LFPs recorded at the p = 100 electrode tips of the Utah array at the k-trial are given by the following forward model: (11) where is the leadfield matrix and where denotes measurement noise which we assume to be normally distributed with expectation zero and covariance matrix (assumed to be the same for each trial) and to be independent across trials. Note that we assume here that the evoked CSD is the same on every trial. We thus adopt the “signal-plus-noise model” for evoked responses, according to which the CSD can be written as the sum of an evoked response (the “signal”) and a term that models the spontaneous background activity (the “noise”). The signal is assumed to be the same for every trial and the noise is assumed to be not-locked to the stimulus and therefore to be independent over trials. The background activity is absorbed into the measurement noise term. Averaging both sides of Eq (11) over trials yields the trial-averaged forward model (12) where V and ξ denote the trial-averaged potentials and measurement noise, respectively. Note that the covariance matrix Σ_ξ of ξ equals divided by the number of trials. In the simulations we can assume assume both to be proportional to the identity matrix since the generally large number of trials allows an accurate sample estimate of Σ_ξ, which can subsequently be used to prewhiten the data [32, 35]. Eq (12) describes the trial-averaged forward model whose inversion is the main aim of this study.

b1. General solution.

Inverse modeling aims at “inverting” the forward model for evoked LFPs (Eq (12)), that is, to reconstruct C from observed V, given the leadfield matrix G. In distributed source modeling, the reconstructed CSD, denoted by , depends linearly on the observed data through an inverse matrix, denoted by G^♯: (13) where G^♯ is given by (14) where Σ_ξ denotes the trial-averaged noise covariance matrix and can be interpreted as an a priori covariance matrix for Vec(C) [31]. Indeed, Eq (14) can be derived by following the standard Bayesian procedure using a (multivariate) normally distributed prior on C. This class of source modeling methods is called “distributed” because the reconstructions are allowed to be distributed across the source-space, in contrast to EEG/MEG dipole localization methods, in which the number of sources is restricted [31]. Although Eq (14) is the general solution to the inverse problem, in order to obtain reasonable reconstructions, an inter-laminar CSD profile needs to be assumed [26, 27]. To incorporate this assumption into the CSD covariance matrix, we denote the inter- and intra-laminar profiles of C by and , respectively (“h” and “v” stand for “horizontal” and “vertical”, respectively.). The above assumption means that C can be decomposed as (15) where C_v assumed to be known. Thus, for each intra-laminar location, the CSD is given by C_v up to a multiplicative constant. Since the vertical CSD profile is assumed to be known, the forward model (Eq (12)) can be reduced to the following horizontal forward model: (16) where V and ξ are as in Eq (12) and is the horizontal leadfield matrix, which, for each intra-laminar location, is obtained by taking the inner-products of the inter-laminar entries of G with C_v. The horizontal forward model makes explicit that only C_h needs to be reconstructed. Note, however, that G_h depends on the choice for C_v so that a different choice leads to a different forward model. The solution to the horizontal forward model is given by the following inverse matrix: (17) where denotes the a priori covariance matrix of C_h. The distributed inverse methods that we consider (MNE, WMNE, LORETA, and LORETA*) differ only in the choice of Σ_h.

b2. Choices for the a priori covariance matrix.

The minimum norm estimate (MNE) is a popular inverse method within the field of EEG/MEG [31, 36] and corresponds to taking the a priori covariance matrix Σ_h to be proportional to the identity matrix. The a priori covariance matrix thus has the following form: (18) where denotes the a priori source variance and denotes the identity matrix of dimension . This means that the MNE assumes the neural currents at two different intra-laminar locations to be uncorrelated and having equal strength. When applied to EEG/MEG recordings, the MNE is known to overestimate superficial sources, for example those located on gyral crowns, and to underestimate deeper sources, for example those located on sulcal walls and fundi. This undesirable property of the MNE is known as surface bias [31] because MNE reconstructions tent to concentrate all source power at surface locations in the brain. In Section b2 of Materials and Methods we will see that in the case of two-dimensional LFP recordings, surface bias takes the form of overestimating currents in the proximity of the electrode tips.

One way to reduce surface bias is to counterbalance the bias by weighting the a priori source covariance matrix. This gives the weighted minimum norm estimate (WMNE) [31]. The WMNE corresponds to the following a priori choice for the intra-laminar covariance matrix: (19)

Note that in the special case , the WMNE reduces to the MNE. In applications to EEG/MEG data, W is usually taken to be a diagonal matrix, thereby reducing the choice of weights to specifying a weighting vector w (the diagonal of W). Most often the entries of w are chosen to be powers of the Euclidean norms of the corresponding leadfields: (20) where G_•,i denotes the i-th column of G (the i-th leadfield) and where q ≥ 0 is a weighting parameter [36, 37]. Although q is usually set to 0.5, its optimal value depends on several factors and determining its optimal value is an empirical issue [37].

Another way of dealing with surface bias is low resolution electrical tomography (LORETA), which combines the weighting matrix W of WMNE with a constraint on spatial smoothness [31, 38] and has already been successfully applied to three-dimensional LFP recordings [29]. LORETA thus corresponds to the following a priori intra-laminar covariance matrix: (21) where denotes the discrete two-dimensional Laplace operator, which can be written in terms of the (Kronecker) tensor sum ⊕ as (22) where Δ_xx and Δ_yy denote the discrete Laplace operators in the x and y directions, respectively. Since the Laplace operator is a spatial lowpass filter (it performs a local averaging), LORETA biases the reconstructions to spatially smooth ones, which reduces surface bias. W we take identical to the weighting matrix defined in the previous section. Because this choice of W differs from that in the original version of LORETA [38], we also considered a version of LORETA without weighting, denoted as LORETA*, which corresponds to the following prior structure on the intra-laminar covariance matrix: (23)

Table 1 lists the different inverse methods.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Linear distributed methods considered in this study.

First column: abbreviations of the methods’ names: (W)MNE = (weighted) minimum norm estimate, LORETA = low resolution electrical tomography, LORETA* = LORETA without weighting. Second column: corresponding prior structure on the intra-laminar CSD covariance matrix. W denotes a weighting matrix and Δ denotes the two-dimensional discrete Laplacian operator.

https://doi.org/10.1371/journal.pone.0187490.t001

b3. Model tuning.

To actually calculate a CSD reconstruction from observed LFPs, an estimate of the trial-averaged measurement noise covariance matrix Σ_ξ is required. Because Σ_ξ is proportional to the single-trial noise covariance matrix , it suffices to estimate the latter. This is usually done by averaging the sample covariance matrices of the pre-stimulus data over all trials. Since in event-related studies, there typically are a large number of trials, this estimate will be accurate [35]. The multiplicative factor relating Σ_ξ to combines with one due to the a priori source variances, yielding a free parameter λ > 0 in front of Σ_ξ in the general solution (Eq (17)), and hence to an inverse matrix G^♯ that depends on λ. To obtain a CSD reconstruction, λ needs to be chosen, a procedure referred to as model tuning. An appropriate value of this regularization parameter can be derived from the observed LFPs and a robust and fast method to do this is generalized cross-validation (GCV) [39] and has already been shown to work for tuning LORETA reconstructions of three-dimensional CSDs [29]. GCV chooses the value of λ for which the following function g is minimized: (24) where 1_p denotes the p-dimensional identity matrix (p is the number of electrodes), V the observed LFPs, and Tr denotes the matrix trace. GCV is fast and it can be shown to approximate the value of λ obtained by using the computationally expensive leave-one-out cross-validation [39].

We numerically determined the minimum of the function g by evaluating it in the values λ = 10⁻²⁰, 10⁻¹⁹, ⋯, 10⁵. It is known that GCV generally works best when the measurement noise is uncorrelated, that is (Σ_ξ = 1_p). In all simulations, we will assume this to be the case and note that the case of correlated noise can be reduced to this case by pre-whitening the data [32]. For evoked responses, this can be done since the large number of trials allow Σ_ξ to be estimated accurately.

In Section b2 of Materials and Methods, we assess the bias β of the difference inverse methods, which, for the general forward model Eq (12), is given by (25) where C is the true CSD and R_λ = G^♯G denotes the resolution matrix associated with the inverse matrix G^♯. Generally, for Tikhonov estimators, ||β|| increases as a function of λ, which reflects the trade-off between bias and uncertainty. In particular, in the absence of measurement noise (λ = 0), the bias is minimal, and it is referred to as the projection bias. In Section b2 of Materials and Methods, we evaluate the projection bias of the different inverse methods. Due to the finite accuracy of numerically computing the inverse operator, we do not set λ to zero, but to an extremely small value (λ = 10⁻³⁰).

b4. Measuring performance.

When applying inverse methods to event-related EEG/MEG data, one is often interested in estimating the locations of relatively localized sources. For this reason, the performance of distributed inverse methods is often characterized using resolution matrices [35]. Although resolution matrices contain information on the accuracy of the reconstructions for localized sources, for more extended sources they are less informative [40]. In the case of multi-electrode LFP data, the dipole approximation which is central to inverse modeling of EEG/MEG data is no longer valid and one generally deals with extended source distributions, rather than with discrete point sources. For this reason, we chose to measure the performance of the different inverse methods using the relative mean squared error, which is defined below.

Let be a simulated CSD and let be its reconstruction obtained by applying one of the inverse methods to the horizontal forward model (Eq (16)). We compare and C_h by calculating the relative means squared error (rMSE): (26) where the factor 10² is included to express the rMSE as a percentage (instead of a fraction). There are two complications in using this performance measure. First, because C_v will generally not be equal to the true inter-laminar profile, C_h and do not have the same scale and therefore cannot be compared directly. Therefore, before calculating the rMSE, is scaled by a constant α so as to minimize , which is achieved for . Second, since the horizontal extent of the source space is larger than that of the electrode array, direct calculation of the rMSE will yield relatively high values. This is because the CSD lateral of the electrode array is generally not well reconstructed due to the high dependency of the corresponding leadfields (see S1 Fig). We therefore calculated the rMSE by first restricting C_h and to intra-laminar locations that are covered by the array. The resulting values hence measure the ability of the inverse methods to reconstruct that part of the CSD that is covered by the electrode array. In the rest of the paper, we simply refer to the rMSE as the (reconstruction) error.

c1. Simulated data.

To test the different inverse methods, we simulate LFPs by applying the trial-averaged forward model to simulated event-related CSDs. Because the CSDs are assumed to be decomposable into an intra- and inter-laminar profile, simulation of an evoked CSD at a given latency amounts to specifying these profiles (C_h and C_v). The intra-laminar profiles of experimental event-related responses are diverse and depend on several factors, including the species under consideration, recording region, response latency, state of the preparation, stimulus properties (frequency, duration, strength, etc.), and spontaneous background activity [2, 41–46]. Evoked responses can be confined to single cortical columns (< 500 μm), as is the case for rat barrel cortex after weak stimulation of individual whiskers [43] or be complex and distributed (up to several millimeters in extent) containing depolarizing (source) as well as hyperpolarizing (sink) components, as is the case for odor-induced responses in salamander olfactory bulb [41]. The inter-laminar profiles of evoked responses display similar diversity and stimulus dependencies [1, 14, 47]. Furthermore, while the CSD is commonly assumed to be balanced, recent studies have reported observing unbalanced and monopolar current sources, probably arising through ionic diffusion processes [22, 23].

To be relevant to a broad range of experimental recordings, therefore, we tested the inverse methods on both simple as well as complex and distributed evoked responses. Specifically, the complex-valued intra-laminar current profile C_h is modeled as a superposition of N two-dimensional Gaussian densities: (27) where (x_n, y_n) denotes the intra-laminar location of the n-th density, ϕ_n denotes a random phase (which gives rise to depolarizing and hyperpolarizing components), and γ_h denotes the spatial width of the densities, which is used to control the spatial scale of C_h. We set N = 100 and selected the locations randomly from the intra-laminar source space. There is no compelling reason to set N = 100 other than that it yields reasonably looking source distributions. The simulation parameters should hence not be taken as claims about the true number of generators of experimental evoked responses; the merely serve to generate a large number of test currents to evaluate the imaging methods. In addition to these responses, we simulated localized responses by setting N = 1, ϕ = 0, and (x₁, y₁) = (3.4, 3.4) (that is, in the center of the electrode array).

The inter-laminar current profile C_v is modeled as a superposition of two one-dimensional Gaussian densities of equal amplitude and opposite sign and centered at depths z₀ ± L/2: (28)

This parametrization corresponds to a dipolar profile with poles located at depths z₀ ± L/2 thus having length L. We set L = 0.8 mm and treat z₀ as a free parameter (see below). Furthermore, γ_v models the width of the poles and is kept fixed at γ_v = L/3.

The trial-averaged evoked LFPs are subsequently calculated by applying the leadfield matrix (Section Utah forward model for cortical evoked responses) and subsequently adding measurement noise. We take the covariance matrix Σ_ξ of the trial-averaged measurement noise to be a diagonal matrix with variance . Because the variance of the measurement noise is inversely proportional to the number of trials, which varies from study to study, we consider different noise-levels. Specifically, following [29], we set (29) which expresses the variance of the noise as a percentage β of the (sample) variance of the noise-free vector of recorded potentials V, where β ranges from 1 to 20 in steps of 2.

In all simulations except the high-resolution simulations in Section b2 of Materials and Methods, we vary two key parameters: the width of the intra-laminar current profile (γ_h) and the depth of the current generator (z₀). Specifically, each of these two parameters takes on two different values (see Table 2). This means that when assessing the effect of errors in the inter-laminar current profile, four sets of simulations are performed for each of the five noise-levels. For each of the four combinations, the reconstruction errors of every inverse method are averaged over 500 independent realizations. To be able to accurately compare the results obtained for different parameters and different types of errors in the a priori inter-laminar current profile, we used the same 500 realizations of all random variables appearing in the simulations (the intra-laminar locations, initial phases, and measurement noise). Fig 2A shows that inter-laminar current profiles of the superficial and deep generators and Fig 2B shows two realizations of their intra-laminar profiles.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Free parameters and their values.

Listed are the free parameters, their symbols, units, and the two values they take on (Value 1 and Value 2).

https://doi.org/10.1371/journal.pone.0187490.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Simulation of evoked current source densities.

A. Inter-laminar current source density (CSD) component modeled by superficial (left) or deep (right) dipolar profiles. Deep and superficial are with respect to the electrode array, which is located at a depth of 1 mm (red line). Zero mm corresponds to the pial surface. B. Two (independent) realizations of simulated intra-laminar CSD components corresponding to local (left column) and global (right column) activations. Red and blue correspond to current sources (depolarization) and current sinks (hyperpolarization), respectively. Black dots denote the electrode locations (400 μm inter-electrode spacing).

https://doi.org/10.1371/journal.pone.0187490.g002

c2. Experimental data.

Electrophysiological data were collected from primary visual cortex of a macaque monkey (Macaca mulatta) that was performing a motion discrimination task. All experimental procedures were in accordance with the animal welfare guidelines of EU directive 2010/63/EU. Ethical review and permission for this work was granted (F149/05) by the regional board Regierungspräsidium Darmstadt. The monkey was group housed with other macaques in facilities of the Ernst StrÜngmann Institute for Neuroscience in accordance with German and EU regulations. The facility provides an enriched environment including toys, wood, natural daylight access and exceeds the size requirements of EU regulations. The monkey received unrestricted access to food and fluids for the duration of the study. On training and recording days, fluid access was controlled contingent on performance during the task. The monkey was bred and purchased from Health Protection Agency, Salisbury, UK. At the end of the study all implants for electrophysiological data collection were removed and the monkey continued to live in his group.

For array and headpost implantation surgeries, anesthesia was induced with a Ketamin/Dexmeditomidin injection and maintained with volatile Isoflurane. Pain was managed with Remifentanil. The data analyzed here is from the trial period where the monkey only had to keep fixation. Black and white square-wave gratings (80% contrast, horizontal orientation, static, 2° diameter, 1.25 cycles/°, surrounded by a 2 px wide annulus with 2.6° diameter) were presented on a 24” Samsung 2233RZ screen at 120 Hz refresh rate at a viewing distance of 86 cm with gray background. The monkey had to keep fixation within 0.6° radius on a small fixation spot in the center of the screen. Stimuli were presented at 4.5° of eccentricity. Local field potentials (LFPs) were acquired from a 64 multi-electrode grid (“Utah” array, 8 × 8 layout, 400 μm inter-electrode spacing, 1 mm electrode length) using a CerePlex E headstage and connected to a Cerebus System (Blackrock Inc.) at 30 kHz sampling rate. The reference electrode was a small wire (∼ 2 mm) reaching out of the array. For the analysis, 535 number of trials were available.

b1. The Utah leadfield matrix.

The Utah leadfield matrix G—which maps cortical currents to extra-cellular potentials—contains the electrical properties and geometry of the tissue and recording set-up, and as such, sets limits on what can be recovered by the different inverse methods. Therefore, before considering the inverse methods, it will be instructive to have a closer look at the leadfield matrix itself and especially the high-resolution leadfield matrix. Although the a priori choice of an inter-laminar current profile reduces the inverse problem to a two-dimensional (intra-laminar) problem by collapsing the inter-laminar dimension of G. in this section we consider the full three-dimensional leadfield matrix. For one thing, considering the inter-laminar dimension of G aids in understanding how inter-laminar volume-conduction effects “contaminate” cortical LFPs [48]. We analyze G through the sensitivities of its constituent leadfields.

The sensitivity of the k-th leadfield of G, that is, its k-th column, is defined to be its Euclidean norm. It is a measure for the strength with which the corresponding monopolar current of unit amplitude contributes to the array potentials. Because the columns of the LFP leadfield matrix are the leadfields of monopolar currents, all entries are positive. This is in contrast to MEG/EEG leadfields, which contain the sensor-projections of dipolar currents and hence generally contain both positive and negative entries. Fig 4A shows the sensitivity profiles of G along three horizontal slices through the modeled tissue volume. The left panel shows the sensitivity profile at the depth of the array (1 mm). It shows that currents at the electrode tips (the black dots) contribute the strongest to the recorded LFPs and that sensitivity drops sharply lateral from the array. The middle and right panels shows that with increasing depth (relative to the electrode array) the sensitivity profile becomes more homogeneous (middle and right panels) and decreases steadily. Fig 4B shows the sensitivity profiles for three vertical slices through the modeled tissue. The fast decrease in sensitivity with depth is clear from the left and middle panels, while the right panel shows that the sensitivity for locations lateral to the electrode array is relatively homogeneous.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Sensitivity of the Utah leadfield matrix.

A. Sensitivity profile of the leadfields along three intra-laminar (horizontal) slices at different cortical depths (left: 1 mm, middle: 1.4 mm, right: 1.8 mm). The Utah array is located at a depth of 1 mm. B. Sensitivity profile of the leadfields along three vertical slices at different lateral locations (left: through the center of the array and on the electrode line, middle: through the center of the array and off the electrode line, right: 200 μm lateral to the array). In all panels, the same color-scaling has been applied so that the sensitivities can be directly compared. Green and red correspond to low and high values, respectively. Black dots in A and B denote the electrodes of the 10 × 10 Utah electrode array (400 μm inter-electrode spacing).

https://doi.org/10.1371/journal.pone.0187490.g004

Since the sensitivity of the LFP leadfields decreases fast for locations outside the electrode array (in both intra- and inter-laminar directions), it might seem as if volume-conduction does not pose difficulties in the interpretation of LFP data. However, the contribution of a current at a particular location is obtained by weighting the locations’ sensitivity with the strength of the current at that location and this is why LFPs can contain contributions from neural activity several millimeters or even centimeters away. In [1, 48], for example, visually and auditory evoked potentials are shown to contain contributions of different belts of primary sensory cortices. It also explains the often observed discrepancies between simultaneously recorded voltage-sensitive dye (VSD) signals—which are not contaminated by volume-conduction—and LFP recordings from the same region [13, 49].

b2. Proximity bias.

The MNE inverse method solves a penalized least squares problem in which the penalty is proportional to the power (that is, the squared Euclidean norm) of the current distribution [31]. When applied to EEG and MEG data, this has the undesirable consequence that the power of reconstructed current distributions tends to concentrate in locations that are “electrically close” to the EEG electrodes/MEG sensors. More precisely, current power is over- and under-estimated at locations whose corresponding leadfields have high and low sensitivity, respectively. Since electric and magnetic fields attenuate with distance, EEG and MEG leadfield sensitivities correlate with physical proximity. The consequence is that current power will be overestimated in cortical gyri—which are close to the electrodes/sensors—and underestimated in cortical sulci. For EEG and MEG data, this “surface bias” in MNE reconstructions can be reduced by weighting the leadfield matrix with the inverse leadfield sensitivities (yielding the WMNE), by spatial smoothing (yielding LORETA), or by appropriate normalization of the MNE reconstructions (yielding dSPM and sLORETA) [31].

Although in the case of intra-cortical LFPs there is no surface bias, leadfield sensitivities can still be analyzed and it is important, therefore, to consider the presence of a more general bias, which we will refer to as proximity bias. Although the sensitivities of the LFP leadfield matrix are all finite—due to the use of voxels instead of point monopoles—in the previous section we have seen that the leadfields at the electrode locations are particularly sensitive (Fig 4A, left panel). This suggests that MNE reconstructions might be biased towards electrode locations and that one of the other imaging methods might perform better. Fig 4A also showed that proximity bias is largest in the intra-laminar plane containing the electrode array. To assess the presence of a proximity bias, we therefore simulated CSDs confined to a single high-resolution intra-laminar slice containing the electrode array.

Fig 5A shows the reconstruction errors averaged over 100 independently generated evoked responses with intra-laminar spatial width equal to 0.45 mm and measurement noise set to zero. Note that MNE has relatively large errors compared to WMNE, LORETA, and LORETA*. These errors are largely due to proximity bias as becomes clear when inspecting the reconstructions (Fig 5C) of a single simulated CSD (Fig 5B). Note that MNE overestimates the currents at the electrodes locations and (slightly) underestimates the currents at locations in between the electrodes. In other words: MNE tends to current power in the electrode locations, that is, it suffers from proximity bias. Although weighting of the leadfields (WMNE) yields smaller errors (Fig 5A), it does not completely remove the bias. Following [37], we have tested a range of values for the weighting parameter p, but this did not substantially reduce the bias. Fig 5A also shows, however, that WMNE is more effective in reducing proximity bias than normalization of MNE reconstructions (dSPM and sLORETA), which tends to overcompensate, leading to underestimation of current power at the electrode locations (see Fig 5C). We did not consider alternative normalization matrices and it therefore remains an open question to what extent dSPM and sLORETA can be adjusted to more effectively reduce proximity bias. It seems that LORETA* reconstructions are free of proximity bias, which is due to its spatial smoothing and absence of weighing (as is done in LORETA).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Proximity bias in high-resolution LFP imaging.

A. Reconstruction errors averaged over 100 independently generated evoked responses for each of the inverse methods. Intra-laminar spatial width of the responses was set to 0.45 mm, measurement noise was set to zero, and the responses were confined to the intra-laminar slice containing the electrode array. B. Single realization of a simulated evoked response (showing only the part that is covered by the electrode array). C. Reconstructions of the response in B. using the different inverse methods (MNE, WMNE, LORETA*, LORETA, dSPM, and sLORETA).

https://doi.org/10.1371/journal.pone.0187490.g005

In the remainder of the manuscript, we restrict to imaging using the low-resolution source-space for which proximity bias is absent. We do not further consider dSPM and sLORETA as they practically yielded the same results as MNE.

c1. Performance of the imaging methods.

In this section we assess the performance of the inverse methods in the absence of errors in the a priori inter-laminar current profile. In other words, we assume that the current profiles are known a priori. The average reconstruction errors over 500 realizations are displayed in Fig 6 (solid lines). Before we consider performance differences between the methods, we have a look at the properties of the current reconstructions that are common to all methods. First, irrespective of the type of activation (local/global and superficial/deep), reconstruction errors increase with increasing noise-level. MNE and WMNE reconstruction errors for global activations seem to be an exception, but in this case, the average reconstruction errors are less accurate because of their large variance over realizations. Such a decrease in performance with increasing noise-level is to be expected because higher noise-levels lead to stronger regularization, which increases the bias in the reconstructions. Second, for noisy data, the errors vary more across measurements. This is because higher noise levels make it more difficult to select an appropriate value for the noise regularization parameter. This fact is relevant for experimental data, because in addition to yielding high errors on average (that is, over a large number of data-sets), high noise-levels make the reconstructions more uncertain. We also note that although the inverse methods considered in this study are linear, this is only true when the value of the noise regularization parameter is set. Indeed, when model tuning is viewed as part of the inverse method, the methods are non-linear. Third, deep generators are more difficult to reconstruct than superficial ones, especially for local activations. The reason for this is that the forward model is a spatial low-pass filter: spatial detail (short-wavelength activity) in the intra-laminar current profiles is suppressed and hence is harder to recover, especially in the presence of measurement noise. This is why deep generators are difficult to reconstruct especially for local activations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Performance of the imaging methods.

A. Mean reconstruction errors for the four inverse methods (MNE (blue), WMNE (red), LORETA (green), and LORETA* (black)) as a function of noise-level and for each of the four combinations of simulated currents (superficial/deep and local/global). Noise-levels are 1, 5, 10, 15, and 20%. The mean errors were obtained by averaging over 500 realizations. B. Same format as in A. but displaying the error standard deviations instead of their means. In A and B, the solid lines correspond to the case of no mismatch in the a priori inter-laminar current profiles. The dashed lines correspond to the case of a mismatch in the a priori inter-laminar current profiles (see text).

https://doi.org/10.1371/journal.pone.0187490.g006

Concerning the performance differences between the methods, we make two remarks. First, MNE and WMNE perform roughly equal and the same holds for LORETA and LORETA*. This means that weighting of the leadfields only has a small effect on the resulting reconstructions. To explain this, recall that for the low-resolution source space, each voxel that is covered by the electrode plane corresponds to an electrode: their are no voxels in-between electrodes (see Fig 1C). Therefore, the leadfields in the electrode plane have similar norms so that weighting by (a power of) their norms effectively doesn’t make a difference. For high-resolution imaging, leadfield weighting does make a difference (see Section b1 of Results). Second, irrespective of the noise-level and generator depth; for local activations, the performance of MNE (and WMNE) and LORETA (LORETA*) is similar, while for global activations, LORETA performs much better than MNE. This is to be expected since global currents better agree with the a priori assumptions of LORETA (spatial smoothness). Corresponding results for unbalanced sources are shown in S2 Fig.

As described in Section b1 of Materials and Methods, the different imaging methods require the specification of an a priori inter-laminar current profile. Since generally we cannot expect this profile to be accurately known, it is important to consider the methods’ performance in case of a mismatch between the true and a priori profiles. We conducted the same sets of simulations as above, but with the difference that in case of a superfical/deep generator, we a priori assume a deep/superficial generator (see Table 2). The results are shown in Fig 6 (dashed lines). We make several remarks. First, the mismatch affects LORETA/LORETA* differently than MNE/WMNE in that for LORETA/LORETA*, the changes in error means and standard deviations are independent of the noise-level, while for MNE/WMNE they depend on the noise-level. Second, for all methods, the increases in error means are larger for superficial than for deep generators. Third, while the mismatch lead to increases in the error variances for superficial generators, it generally leads to decreases for deep generators. The latter reflects a more stable selection of the regularization parameter λ. The cause of this is that if a deep generator is (erroneously) assumed to be superficial, its intra-laminar spatial extent will be overestimated, in order to fit the observed data. This is due to the fact that the LFP forward model acts like a spatial lowpass filter. Especially for local activations and noisy measurements, larger (assumed) sources are easier to distinguish from spatially uncorrelated noise, which is what (generalized) cross-validation tries to do when selecting an appropriate value for the regularization parameter. This also explains the increased error variance for superficial generators because (erroneously) assuming superfical generators to be deep, leads to underestimation of their intra-laminar spatial extent, which renders model tuning more difficult.

c2. The effect of electrode-montage.

In the simulations above we have used the forward model given by Eq (12), which describes the relation between the CSD and the ensuing (trial-averaged) absolute electrical potential at the recording electrodes, that is, the potentials referenced to an electrode located at infinity. In practice, however, the reference electrode must be in contact with the preparation because electric potentials are measured indirectly via the current between the recording and reference electrode. For in vivo cortical LFPs, the reference electrode is often located at the surface of the contra-lateral cortex or on the skull. Because existing LFP volume-conductor models are local [9, 24] and therefore cannot simulate LFPs that are referenced to a distant electrode, the LFPs have to be re-referenced, prior to the estimation of the CSD. A common misconception in the field of EEG imaging is that switching to a different electrode montage (that is, re-referencing) does not influence the source reconstructions and this issue is given surprisingly little attention in the literature. In fact, most EEG (and ECoG) inverse modeling studies assess the performance of inverse methods on absolute potentials only [31] and whose practical relevance therefore, is limited. As far as we know, reference issues have not been discussed in the literature on LFP inverse modeling and existing studies have focused exclusively on absolute potentials [25–28, 30].

It is straightforward to show that re-referencing changes the LFP forward model and hence the CSD reconstructions. Consider again the horizontal forward model for absolute potentials (Eq (16)): (30) where , , and where denotes Gaussian measurement noise with expectation zero and covariance matrix Σ_ξ. Furthermore, let for certain q ≤ p denote a montage transformation that transforms the absolute potentials V to the re-referenced potentials V_M: (31)

The forward model for the re-referenced potentials is given by (32) where denotes the re-referenced horizonal leadfield matrix and (33) denotes the re-referenced measurement noise, which has covariance matrix MΣ_ξ M^t. The re-referenced horizontal forward model shows that re-referencing has two effects: It changes the leadfield matrix and it changes the covariance structure of the measurement noise. Because we have assumed the data to be prewhitened, MΣ_ξ M^t is a scaled copy of the identity matrix. Moreover, since we will focus on the average-reference montage, which corresponds to taking , where denotes the vector containing all ones, MM^t is close to the identity matrix and hence MΣ_ξ M^t approximately equals Σ_ξ. We therefore focus on the effect that re-referencing has on the leadfield matrix.

To assess the effect of re-referencing on the quality of the CSD reconstructions, we carried out the standard set of simulations, but used average-referenced potentials instead of absolute potentials. Fig 7 shows the difference in mean reconstruction errors with those obtained by using single-wire potentials. Note that for local activations, changing from absolute potential to average-reference potentials has practically no effect on the mean reconstruction errors and this holds for all methods. In contrast, for global activations, changing to the average-reference montage drastically increased the mean errors. Also note that these effects are independent of the noise-level. It turns out that the large errors for global activations can largely be accounted for by errors in the off-set of the reconstructions. This we checked by recalculating the errors modulo an additive constant: the resulting errors were practically the same as those obtained from the absolute potentials. The errors in off-set arise because imaging average-reference potentials yields approximately balanced (that is, sum-zero) intra-laminar current profiles. Since localized activations are approximately balanced, error don’t increase much. Global activations, however, tend to be unbalanced (see Fig 2B), giving rise to larger errors when using the average-reference montage.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Effects of re-referencing the data.

Differences in mean reconstruction errors obtained using the average-reference montage and the single-wire montage for the four inverse methods (MNE (blue), WMNE (red), LORETA (green), and LORETA* (black)) as a function of noise-level and for each of the four combinations of simulated currents (superficial/deep and local/global). Noise-levels are 1, 5, 10, 15, and 20%. The mean errors were obtained by averaging over 500 realizations.

https://doi.org/10.1371/journal.pone.0187490.g007