Surface Laplacian and CSD

Why this matters

The surface Laplacian is the cheapest useful answer to L5.1. It needs no MRI, no head model beyond a sphere, no source space and no inverse operator — only the electrode positions you already have — and it removes the reference from the problem entirely while sharpening the map. It is also the transform most often applied without stating what it throws away, which is a great deal.

Concepts

What the transform is

The surface Laplacian is the second spatial derivative of the scalp potential, taken across the scalp surface. Under a set of assumptions about the head — a locally homogeneous, locally spherical volume beneath the electrodes — it is proportional to the current density entering and leaving the skull beneath each electrode, which is why the result is usually called current source density (CSD).

Two things follow immediately from “second spatial derivative”, and they are the whole lesson.

It is reference-free. Re-referencing subtracts the same number r(t) from every channel at each time point (L2.3). A constant has zero spatial derivative, so it has zero Laplacian: the transform of the re-referenced data is identical to the transform of the original. Not approximately — the reference term is annihilated exactly. The widget shows this directly: switch the reference with the CSD panel open and the CSD map does not move.

It is a spatial high-pass filter. Differentiating twice multiplies each spatial frequency by roughly the square of that frequency, so broad, smooth patterns are suppressed and sharp, local ones are kept or amplified. Everything CSD does well and everything it does badly comes out of that one sentence.

How it is actually computed

You cannot difference a surface you have sampled at 30 points, so the standard implementation fits a smooth interpolating surface first and differentiates that. The usual family is spherical splines: fit a function on a sphere that passes near the measured values while minimising a bending-energy penalty, then take the Laplacian of the fitted function analytically. This is the same interpolation machinery as the topographies of L3.4 and the bad-channel interpolation of L2.2, used for a different purpose.

Two parameters control it, and both should be reported:

• Stiffness / spline order (m, stiffness in MNE, often 3 or 4). Higher values give a smoother fitted surface and therefore a less aggressive sharpening.
• Regularization (lambda2). How far the fitted surface is allowed to deviate from the measured values. Larger values smooth more and are the defence against a noisy channel dominating its neighbourhood.

The asset behind this lesson’s widget was computed with MNE’s compute_current_source_density at its defaults — sphere: auto, lambda2: 1e-05, stiffness: 4, n_legendre_terms: 50 — on the average-referenced grand average. TODO(confirm): the spherical-spline formulation is conventionally attributed to Perrin and colleagues (1989); that citation is not in this site’s reading list and has not been checked against the paper here. (Kayser, 2015) is the reading-list entry for the method’s use in practice.

What it costs

The units change, and they are not comparable with anything referential. CSD is in µV/m² (or, depending on the scaling convention, A/m³). A CSD amplitude cannot be compared with a µV amplitude, cannot be averaged with one, and cannot be quoted against a published referential value. The shipped ERP CORE grand average makes the point bluntly: the target-minus-standard difference at Pz over 300–600 ms measures 3.505 µV under an average reference and 2405.4 µV/m² under CSD. Those are not two estimates of one quantity; they are two different quantities, and the second number’s size means nothing on its own.

Deep and broad sources are attenuated. A deep generator projects a wide, smooth pattern to the scalp (L0.1, and the depth-versus-extent figure in pf-deep-source-claims). A spatial high-pass filter suppresses exactly that. So a component that is genuinely generated by a deep or spatially extended population can shrink or vanish under CSD — not because it was volume-conducted noise, but because CSD cannot see it. “The effect survived CSD” is evidence of a superficial, focal generator; “the effect disappeared under CSD” is not evidence that the effect was spurious.

Noise is amplified. High spatial frequencies include a single channel’s own noise, so one bad or noisy electrode contaminates its neighbourhood, and it does so after the transform in a way that is harder to spot than before it. CSD demands honest bad-channel handling (L2.2) and makes rejection thresholds matter more.

Accurate electrode positions matter more than usual. A derivative is sensitive to the geometry it is taken over. Template positions are an approximation; ds-lemon ships digitized positions for a large number of subjects and is the dataset to use when that difference is the question.

Coverage and density set the spatial scale. The transform can only resolve spatial frequencies the montage samples. On 19 or 30 electrodes the “sharpening” is limited, and near the edge of the montage the fitted surface is extrapolating (L3.4), so the Laplacian there is a property of the spline, not of the head. Report the electrode count and draw the positions.

Using it for connectivity

CSD is a genuine partial answer to L5.1 and it is worth being precise about how partial.

Because the transform is local, each transformed channel is dominated by sources under and immediately around it, so the broad common component that made every sensor pair coherent is strongly suppressed. Connectivity matrices computed on CSD-transformed data are much less uniform than the referential ones, and the reference-dependence problem of L5.2 disappears entirely, which removes one whole family of alternative explanations.

But it does not remove volume conduction. It suppresses the broad part of it. Two nearby electrodes still share sources after the transform — more so, because the sharpening is not a deconvolution — so short-range coherence remains inflated, and the measures of L5.2 are still the right ones to reach for. CSD and wPLI compose; neither makes the other unnecessary. TODO(confirm): quantitative statements about how much CSD reduces spurious zero-lag coupling, and at what spatial scale, should be verified against a primary source before publication.

The data behind this lesson

• ds-erpcore P3, CC BY 4.0, open access: Biosemi ActiveTwo, 30 EEG + 3 EOG electrodes in a 10-20 placement scheme, 1024 Hz, CMS reference, 60 Hz mains, no software filters, 40 participants per paradigm. The widget’s arrays are grand averages over 20 subjects, epoched, baselined over −0.2 to 0 s, re-referenced and resampled to 256 Hz. TODO(confirm): the author mirrors the ERP CORE entry into the catalogue registry and signs off the dataset page; the shipped sidecars also record a licence conflict in the source that the author reconciles.
• The shipped references and their Pz target-minus-standard values over 300–600 ms: original CMS 0.266 µV, Cz 1.676 µV, average 3.505 µV, REST 5.147 µV, linked mastoids 5.6 µV, CSD 2405.4 µV/m². TODO(confirm): the ERP CORE 30-channel montage carries no mastoid electrode, so P9 and P10 — the nearest inferior-posterior pair — stand in for the linked-mastoid reference, and every label that shows it should say so.
• The spread across the five referential numbers is worth more than the CSD number: the same measurement on the same data ranges from 0.27 to 5.6 µV depending only on where the reference sat. That is the sense in which “reference-free” is a substantive property and not a technicality.
• ds-lemon supplies digitized electrode positions for the position-accuracy comparison. TODO(confirm): ds-lemon’s exact licence terms are recorded as unverified in the site’s catalogue.

Explore

What to look for

• Change the reference and watch the CSD map stay still. The referential map swings between red and blue as the zero-crossing moves; the CSD map is pixel-identical. That is reference independence, seen rather than asserted.
• Compare the two maps at the P3’s peak. The CSD map is tighter and has more sign reversals over the same scalp — the same data through a spatial high-pass filter.
• Read the two amplitude axes. µV on one, µV/m² on the other, differing by three orders of magnitude. Convince yourself there is no scale factor that would make them comparable.
• Move the time cursor into a latency where the referential map is broad and smooth. The CSD map is nearly flat there: a broad distribution is what the transform is built to remove, whether it came from mixing or from a genuinely extended generator.
• Look at the outermost electrodes. The spline is extrapolating beyond them, and a derivative of an extrapolation is the least trustworthy thing on the page.

Practice

The notebook applies compute_current_source_density to ds-erpcore and to ds-lemon, showing before-and-after topographies and before-and-after wPLI matrices, and varying the stiffness and regularization parameters so that their effect on both is visible rather than left at defaults.

Exercises