Functional Brain Imaging and Brain–Computer Interfaces

Contents

• Introduction and overview
• 0. Review: vector calculus, linear algebra and the Gaussian linear model
• 1. From neural currents to electromagnetic fields
• 2. The M/EEG forward problem and finite elements
• 3. Preprocessing and denoising M/EEG signals
• 4. The M/EEG inverse problem
• 5. fMRI: the BOLD signal and the general linear model
• 6. Hypothesis testing and multiple comparisons
• 7. Connectivity, networks and covariance learning
• 8. Brain–computer interfaces
• 9. Complements: boundary elements, sparse imaging, random fields, decoding, effective connectivity
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, computer science or physics who wants to turn the electromagnetic or haemodynamic traces of brain activity into maps, statistics and commands, whether for neuroscience, for the clinic or for a brain–computer interface.

Empirical first: Bennett and colleagues (2010) put a dead Atlantic salmon in an MRI scanner, “showed” it photographs of people in social situations, and thresholded the voxelwise statistics at an uncorrected p < 0.001: a few voxels in the fish’s brain came out “active”. Eklund, Nichols and Knutsson (2016) then ran standard group analyses on resting-state data, where no task effect exists, and parametric cluster inference gave familywise false-positive rates far above the nominal 5 %: both failures come from statistics (ch. 6), not from the scanner.

The course is about two non-invasive windows on the working brain. Magneto- and electroencephalography (M/EEG) record, within milliseconds, the fields created by currents in the cortex. Functional MRI records, within millimetres but over seconds, a haemodynamic signal tied to the energy used by brain areas.

The page follows the official syllabus: the physics of the signals, the forward and inverse problems of M/EEG, the statistical analysis of functional images and networks, and brain–computer interfaces (BCIs). It shares the neural modelling of Mathematical Methods for Neuroscience, the regularised inversion of Inverse Problems and Imaging and the clinical statistics of Biostatistics.

What you will be able to do

• Derive the quasistatic equations of bioelectromagnetism and compute the potential and magnetic field of a current dipole.
• Explain why MEG is blind to radial sources in a spherical head, and prove it.
• Write the M/EEG forward problem as a boundary value problem, solve it with finite elements, and build a lead field by reciprocity.
• Denoise M/EEG recordings by filtering, averaging, whitening, signal-space projection and independent component analysis.
• Localise sources with dipole fitting, MUSIC, beamformers and minimum-norm imaging, and validate them with resolution matrices and simulations.
• Model an fMRI time series with the general linear model, test contrasts, and run a group analysis with mixed effects.
• Control false positives in statistical maps with permutation tests, familywise and false discovery rate procedures, and cluster inference.
• Estimate functional connectivity with partial correlations, shrinkage, the graphical lasso and the geometry of covariance matrices, and recognise the effect of volume conduction on M/EEG coherence.
• Build and evaluate a BCI pipeline (common spatial patterns, linear discriminant analysis, Riemannian classifiers) without leakage, and compute its information transfer rate.
• Discuss the role of feedback and user learning in a closed-loop BCI.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states that the course deals with functional brain imaging and its application to brain–computer interfaces, based on two non-invasive techniques that measure brain activity. M/EEG measures the electromagnetic field created by cortical electric currents; functional MRI measures a signal (BOLD) linked to the energy consumption of brain areas.

The description announces a course of 8 sessions. It lists the themes below, grouped in four parts, followed by the heading “Les intervenants” (the speakers). The Université Paris-Saclay course detail adds only the nature of the assessment.

Extensions beyond the syllabus: the review ch. 0; the subtraction source model and discrete reciprocity (ch. 2); the proof that sLORETA has no localisation bias for one source (Theorem 4.8); the Riemannian geometry of covariance matrices (7.5); the complements of ch. 9 and the research directions of ch. 10.

Where this course fits in the programme

Official prerequisites: integral calculus, elements of statistics, numerical computation, linear algebra and optimization. The course is a second-semester Modelling course of 5 ECTS of the Santé track, taught by Bertrand Thirion and Marie-Constance Corsi (INRIA, as listed on master-mva.com).

• Refreshers — Linear Algebra, Analysis and Optimization for the SVD, Tikhonov regularisation and the Lagrange multipliers of ch. 0 and ch. 4; Probability and Statistics for the Gaussian linear model and the tests of ch. 5–ch. 6; Scientific Python, NumPy and PyTorch for the code.
• First-semester courses that prepare it — Mathematical Methods for Neuroscience (where the currents come from); Convex Optimization and Applications in Machine Learning (graphical lasso, proximal methods of section 9.2); Introduction to Statistical Learning (LDA, cross-validation); Computational Statistics (resampling); Machine Learning for Time Series (filtering, spectra); Medical Image Analysis Based on Generative, Geometric and Biophysical Models (MRI, registration, the Santé track).
• Related second-semester courses — Inverse Problems and Imaging (Bayesian regularisation); Large Random Matrices and Applications to Machine Learning (covariance estimation in high dimension); Biostatistics (mixed models, clinical validation); Deep Learning for Medical Imaging; Graphs in Machine Learning (brain networks as graphs); Kernel Methods for Machine Learning; Audio Signal Processing: Time–Frequency Analysis (spectral estimation).
• Alongside and after the course — the Hospital Immersion: Paired Project with a Physician, where clinical EEG and MRI are met in practice, and the Research Internship.

Course outline

• Review: vector calculus, linear algebra and the Gaussian linear model — divergence theorem and Green’s identity, the fundamental solution of the Laplacian, SVD and Tikhonov regularisation, Rayleigh quotients, the distribution of least squares.
• From neural currents to electromagnetic fields — the origin of M/EEG signals, the quasistatic Maxwell equations, the current dipole, Biot–Savart, spherical models and the silence of radial dipoles, source models.
• The M/EEG forward problem and finite elements — the boundary value problem, its weak form, P1 finite elements, singular sources, the lead field and Helmholtz reciprocity.
• Preprocessing and denoising M/EEG signals — filtering and averaging, noise covariance and whitening, signal-space projection, independent component analysis.
• The M/EEG inverse problem — ill-posedness, dipole fitting, MUSIC, beamformers, minimum-norm imaging, dSPM and sLORETA, validation.
• fMRI: the BOLD signal and the general linear model — the haemodynamic response, design matrices, contrasts, prewhitening, mixed-effects group models.
• Hypothesis testing and multiple comparisons — permutation tests, familywise error rate, false discovery rate, cluster-level inference.
• Connectivity, networks and covariance learning — partial correlations, shrinkage, the graphical lasso, the manifold of covariance matrices.
• Brain–computer interfaces — paradigms, information transfer rate, common spatial patterns, linear and Riemannian classifiers, evaluation, feedback.
• Complements — boundary elements and the Geselowitz formula, Maxwell filtering, sparse source imaging, random field theory, decoding and the interpretation of weights, M/EEG coherence and volume conduction, effective connectivity (self-study).
• Research directions and open problems — reproducibility, large cohorts, deep decoding of non-invasive signals, transfer across sessions and subjects, user training (self-study).

Session plan

24 h in class = 8 sessions of 3 h.

The volume is the Université Paris-Saclay programme table (lectures 24 h), and master-mva.com announces a course of 8 sessions, which gives 3 h per session. The split of the content into sessions is a proposal: it follows the four parts of the official syllabus in their order.

Level 1 = foundations, 2 = core models and algorithms, 3 = advanced or research-level material; the level never decreases.

Self-study, outside the 24 h (about 8 h in total): ch. 9 Complements: boundary elements and Maxwell filtering, sparse imaging (Code 9.4), random field theory, decoding, M/EEG coherence, effective connectivity (≈ 5 h, with Ex. 18) · ch. 10 Research directions and open problems (≈ 2 h) · the written-exam proposal of Projects and assessment as a self-test (≈ 1 h).

How to study this module

• Read ch. 0 first, even quickly: Green’s identity, the push-through identity and the distribution of least squares are used without warning in ch. 2, ch. 4 and ch. 5.
• Keep the orders of magnitude of Example 1.3 and Example 1.6 in mind: they explain most modelling choices of the M/EEG chapters.
• Run every code block: they use NumPy, SciPy and scikit-learn only, with fixed seeds and synthetic data, and each runs in a few seconds. Then change one parameter (the noise level, the depth of a source, the number of trials) and predict the output before running again.
• After each session, solve the exercises of the Practice column; the ★★★ exercises are good starting points for a TD or a project.
• For real data, the open-source libraries MNE-Python (M/EEG), Nilearn (fMRI) and pyRiemann (BCI) implement the methods of this page; the projects use them on public data sets.
• Read ch. 9 section by section as the matching chapter is taught: 9.1 after session 3, 9.2 after session 5, 9.3 after session 7, 9.4–9.6 after session 8.

0. Review: vector calculus, linear algebra and the Gaussian linear model

Three M1 toolboxes are used: vector calculus for the physics of ch. 1–ch. 2, regularised least squares for ch. 4, and the Gaussian linear model for ch. 5–ch. 6. This chapter states them with short proofs.

0.1 Vector calculus

0.2 Singular values and regularised least squares

The second form costs one m \times m inversion. With m \approx 300 sensors and p \approx 10^4 sources, this is what makes minimum-norm imaging cheap (ch. 4).

0.3 The Gaussian linear model

1. From neural currents to electromagnetic fields

This chapter follows the signal from the synapse to the sensor and reduces Maxwell’s equations to the two formulas used by M/EEG software: a Poisson equation for the potential and the Biot–Savart law for the field.

1.1 Where the signals come from

When a neurotransmitter opens ion channels on the dendrite of a cortical pyramidal neuron, a current enters the cell, flows along the dendrite, and returns through the extracellular medium. The intracellular part of this loop is the primary current; the extracellular return is the volume current.

Pyramidal cells are aligned perpendicular to the cortical surface, so their primary currents add up when synchronous; a measurable signal needs tens of thousands of cells active together. Brief action potentials contribute little at a distance.

The magnetic field of the brain is about 10^8 times smaller than the Earth’s field of 50 µT; MEG needs superconducting or optically pumped magnetometers in a magnetically shielded room.

1.2 Maxwell’s equations in the quasistatic regime

Brain signals lie below about 1 kHz. Two dimensionless numbers then justify dropping the time derivatives: the ratio of displacement to ohmic currents, \omega\varepsilon/\sigma, and the ratio of the head size to the skin depth \delta = \sqrt{2/(\omega \mu_0 \sigma)}.

The time dependence is now a parameter: at every instant the potential is the static response to the primary currents of that instant. This is why M/EEG models are linear and instantaneous.

1.3 The current dipole and the potential in a homogeneous medium

The potential decreases as the inverse square of the distance, with a positive lobe in the direction of q and a negative lobe behind it.

1.4 The magnetic field: Biot–Savart

Both primary and volume currents contribute. Splitting J = J^p - \sigma \nabla V gives B = B_0 + B_v, where B_0 is the field of the primary current alone.

1.5 Spherical head models

In general the volume currents depend on the whole conductivity profile. A sphere is the exception that makes MEG modelling simple.

MEG is thus blind to radial sources, located at the crests of the gyri, and most sensitive to tangential sources in the walls of the sulci. EEG sees both, which is one reason to combine the two modalities.

import numpy as np

mu0 = 4e-7 * np.pi

def b_primary(r, r0, q):
    """Biot-Savart field of the primary current dipole alone (T)."""
    a = r - r0
    return mu0 / (4 * np.pi) * np.cross(q, a) / np.linalg.norm(a) ** 3

def b_sarvas(r, r0, q):
    """Field outside a spherically symmetric conductor centred at 0 (Sarvas, 1987)."""
    a_vec = r - r0
    a, rn = np.linalg.norm(a_vec), np.linalg.norm(r)
    F = a * (rn * a + rn**2 - r0 @ r)
    gradF = ((a**2 / rn + a_vec @ r / a + 2 * a + 2 * rn) * r
             - (a + 2 * rn + a_vec @ r / a) * r0)
    qxr0 = np.cross(q, r0)
    return mu0 / (4 * np.pi * F**2) * (F * qxr0 - (qxr0 @ r) * gradF)

r0 = np.array([0.0, 0.0, 0.07])           # source 7 cm from the centre
r = np.array([0.03, 0.01, 0.105])         # sensor about 4 cm away
rhat = r / np.linalg.norm(r)
q_tan = np.array([10e-9, 0.0, 0.0])       # 10 nAm tangential dipole
q_rad = np.array([0.0, 0.0, 10e-9])       # 10 nAm radial dipole

for name, q in [("tangential", q_tan), ("radial", q_rad)]:
    B, B0 = b_sarvas(r, r0, q), b_primary(r, r0, q)
    print(f"{name:10s} |B| = {np.linalg.norm(B)*1e15:7.1f} fT, |B0| = "
          f"{np.linalg.norm(B0)*1e15:7.1f} fT, radial parts: "
          f"{B @ rhat*1e15:7.2f} vs {B0 @ rhat*1e15:7.2f} fT")

1.6 Source models

Column j of G is the topography that a unit source at position j would produce on the sensors. Computing it is the forward problem of ch. 2; inverting it is the inverse problem of ch. 4.

2. The M/EEG forward problem and finite elements

Given a source, what do the sensors of a realistic head record? This forward problem is a boundary value problem; this chapter solves it with finite elements and uses reciprocity to get a whole lead field from one solve per electrode.

2.1 The boundary value problem

Typical values are about 0.33 S/m for the scalp and the grey matter, 1.79 S/m for the cerebrospinal fluid, and between 0.004 and 0.02 S/m for the skull; the white matter is anisotropic. The skull value is the main source of modelling error in EEG.

The potential is only defined up to a constant: EEG always measures differences, with respect to a reference electrode or to the average of all electrodes.

2.2 Weak formulation

A point dipole is not in L^2: for it, \ell(w) = q \cdot \nabla w(r_0), which is not continuous on H^1. Two classical remedies follow. The direct approach evaluates this formula on piecewise-linear test functions anyway. The subtraction approach removes the singularity analytically (Proposition 2.6).

2.3 Galerkin approximation with P1 elements

For a smooth solution, interpolation gives an error of order h in the energy norm and h^2 in L^2 (Ern and Guermond, 2004). A dipole solution is not smooth near r_0, which spoils this rate unless the singularity is removed.

2.4 Lead field and reciprocity

Computing G column by column takes one solve per source, that is 10^4 solves for a cortical source space. Reciprocity reverses the roles of sources and sensors.

One solve per electrode (relative to the reference) thus gives every row of G. The discrete version is exact: since K is symmetric, (e_a - e_b)^\top K^{+} b = b^\top K^{+} (e_a - e_b), which is the identity of Theorem 2.8 for the FEM solution.

import numpy as np
from scipy.spatial import Delaunay
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve

def disk_mesh(n_rings):
    """Nodes on concentric rings of the unit disk; the last ring is the scalp."""
    pts = [np.zeros((1, 2))]
    for k in range(1, n_rings + 1):
        t = 2 * np.pi * np.arange(6 * k) / (6 * k)
        pts.append(k / n_rings * np.c_[np.cos(t), np.sin(t)])
    p = np.vstack(pts)
    return p, Delaunay(p).simplices, np.arange(len(p) - 6 * n_rings, len(p))

def stiffness(p, tri, sigma=1.0):
    """P1 stiffness matrix and the gradients of the hat functions per triangle."""
    M = np.concatenate([np.ones((len(tri), 3, 1)), p[tri]], axis=2)
    G = np.linalg.inv(M)[:, 1:, :].transpose(0, 2, 1)       # (T, 3, 2)
    area = np.abs(np.linalg.det(M)) / 2
    Ke = sigma * area[:, None, None] * G @ G.transpose(0, 2, 1)
    I, J = np.repeat(tri, 3, axis=1).ravel(), np.tile(tri, (1, 3)).ravel()
    return coo_matrix((Ke.ravel(), (I, J)), shape=(len(p),) * 2).tocsr(), G

def solve_neumann(K, b):
    """Ground node 0 (V is defined up to a constant), then remove the mean."""
    K = K.tolil(); K[0, :] = 0; K[0, 0] = 1; b = b.copy(); b[0] = 0
    v = spsolve(K.tocsr(), b)
    return v - v.mean()

def v_inf(x, r0, q):                       # 2-D infinite-medium dipole, sigma = 1
    d = x - r0
    return (d @ q) / (2 * np.pi * np.sum(d**2, axis=1))

def grad_v_inf(x, r0, q):
    d = x - r0; d2 = np.sum(d**2, axis=1)[:, None]
    return (q / d2 - 2 * (d @ q)[:, None] * d / d2**2) / (2 * np.pi)

def rhs_direct(p, tri, G, r0, q):          # b_i = q . grad phi_i(r0)
    t = Delaunay(p).find_simplex(r0)
    b = np.zeros(len(p)); b[tri[t]] = G[t] @ q
    return b

def rhs_subtraction(p, bd, r0, q):         # b_i = - int dV_inf/dn phi_i ds
    b, (gx, gw) = np.zeros(len(p)), np.polynomial.legendre.leggauss(3)
    s, w = (gx + 1) / 2, gw / 2
    for i, j in zip(bd, np.roll(bd, -1)):
        e = p[j] - p[i]; L = np.linalg.norm(e); n = np.array([e[1], -e[0]]) / L
        flux = grad_v_inf(p[i] + s[:, None] * e, r0, q) @ n
        b[i] -= L * np.sum(w * flux * (1 - s)); b[j] -= L * np.sum(w * flux * s)
    return b

def exact(x, r0, q, h=1e-6):               # image solution in the unit disk
    def N(s):
        return (np.log(np.linalg.norm(x - s, axis=1))
                + np.log(np.linalg.norm(x - s / (s @ s), axis=1))) / (2 * np.pi)
    v = -sum(qk * (N(r0 + h * e) - N(r0 - h * e)) / (2 * h)
             for qk, e in zip(q, np.eye(2)))
    return v - v.mean()

r0, q = np.array([0.3, 0.2]), np.array([1.0, 0.5])
for n in [10, 20, 40, 80]:
    p, tri, bd = disk_mesh(n)
    K, G = stiffness(p, tri)
    va = exact(p[bd], r0, q)
    v1 = solve_neumann(K, rhs_direct(p, tri, G, r0, q))
    v2 = solve_neumann(K, rhs_subtraction(p, bd, r0, q)) + v_inf(p, r0, q)
    err = [np.linalg.norm(v[bd] - v[bd].mean() - va) / np.linalg.norm(va) for v in (v1, v2)]
    print(f"rings {n:2d}, h = {1/n:.4f}: direct {err[0]:.1e}, subtraction {err[1]:.1e}")

a, b = bd[0], bd[len(bd) // 4]             # reciprocity on the finest mesh
inj = np.zeros(len(p)); inj[a], inj[b] = 1.0, -1.0
U = solve_neumann(K, inj)                  # unit current injected at a, out at b
t = Delaunay(p).find_simplex(r0)
print(f"V(a) - V(b) = {v1[a] - v1[b]:.6f}, q . grad U(r0) = {q @ (G[t].T @ U[tri[t]]):.6f}")

Boundary element methods, which mesh only the interfaces, are the main alternative for nested homogeneous compartments (section 9.1). Finite elements handle anisotropy, holes in the skull and arbitrary geometry, at the price of a volume mesh built from a segmented MRI.

3. Preprocessing and denoising M/EEG signals

A raw recording mixes brain activity of interest with eye blinks, heartbeats, muscle activity, line noise and ongoing background activity, often larger than the signal sought. This chapter gives the linear tools that separate them, and states the assumptions each one makes.

3.1 Signal model and sources of noise

3.2 Filtering and averaging

A band-pass filter keeps the band of interest (1–40 Hz for evoked responses, 8–30 Hz for sensorimotor rhythms) and a notch removes the line frequency. Offline, zero-phase filtering avoids delays; online, as in a BCI, filters must be causal.

3.3 Noise covariance and whitening

Whitening puts channels of different units (EEG in volts, magnetometers and gradiometers in teslas and teslas per metre) on a common scale. Least squares on whitened data is generalised least squares, the efficient estimator under correlated noise (Proposition 5.3 states the same fact for fMRI).

3.4 Signal-space projection

Two consequences matter in practice. The lead field must be projected as well (PG), otherwise source estimates are biased. And the component of every brain topography along U_a is lost: a source whose topography resembles a blink is attenuated with it.

3.5 Independent component analysis

After PCA whitening, z = W x has identity covariance and z = Q s with Q orthogonal, since Q Q^\top = \mathrm{Cov}(z) = I. It remains to find a rotation, which non-Gaussianity determines.

FastICA (Hyvärinen and Oja, 2000) maximises an approximation of negentropy with a nonlinearity g and the fixed-point update, followed by normalisation, w \leftarrow \mathbb{E}[z\, g(w^\top z)] - \mathbb{E}[g'(w^\top z)]\, w . Components that correlate with an EOG or ECG channel, or look like blinks, are zeroed before back-projection.

import numpy as np
from sklearn.decomposition import FastICA
from scipy.stats import kurtosis

rng = np.random.default_rng(0)
fs, T = 250, 20                                   # 250 Hz, 20 s
t = np.arange(fs * T) / fs
alpha = np.sin(2 * np.pi * 10 * t) * (1 + 0.5 * np.sin(2 * np.pi * 0.2 * t))
beta = rng.standard_normal(t.size)                # broadband background source
blink = np.zeros_like(t)
for t0 in rng.uniform(0.5, T - 0.5, size=12):     # 12 blinks, about 0.3 s each
    blink += 8 * np.exp(-0.5 * ((t - t0) / 0.06) ** 2)
S = np.vstack([alpha, beta, blink])               # (3 sources, time)

n_ch = 8
A = rng.standard_normal((n_ch, 3))
A[:, 2] = np.abs(A[:, 2]) * np.linspace(2, 0.2, n_ch)   # blink: frontal channels
X_brain = A[:, :2] @ S[:2]
X = X_brain + np.outer(A[:, 2], blink) + 0.1 * rng.standard_normal((n_ch, t.size))

def report(name, Xc):
    D = Xc - X_brain
    D -= D.mean(axis=1, keepdims=True)            # offsets are removed by high-pass filtering
    err = np.linalg.norm(D) / np.linalg.norm(X_brain)
    print(f"{name:22s} relative error to the artifact-free signal: {err:.3f}")

report("raw", X)

# Signal-space projection: first principal direction of the blink periods
mask = blink > 1.0
U, _, _ = np.linalg.svd(X[:, mask] - X[:, mask].mean(1, keepdims=True), full_matrices=False)
u = U[:, :1]
P = np.eye(n_ch) - u @ u.T                        # projector onto the complement
report("SSP (1 component)", P @ X)

# ICA: the blink component is the one with the largest kurtosis
ica = FastICA(n_components=3, whiten="unit-variance", random_state=0, max_iter=1000)
Sh = ica.fit_transform(X.T)                       # (time, components)
k = np.argmax(kurtosis(Sh, axis=0))
Sh[:, k] = 0
report("ICA (blink removed)", ica.inverse_transform(Sh).T)
print("kurtosis of the components:", np.round(kurtosis(ica.transform(X.T), axis=0), 1))
print("angle between blink topography and SSP vector (deg):",
      round(np.degrees(np.arccos(abs(u[:, 0] @ A[:, 2]) / np.linalg.norm(A[:, 2]))), 2))

4. The M/EEG inverse problem

Given the measurements and the lead field, where are the active sources? The problem is ill-posed, and every method adds an assumption: few dipoles, uncorrelated sources, or a prior on the amplitudes. This chapter covers MUSIC, beamforming and imaging, and their validation.

4.1 Ill-posedness

Helmholtz (1853) already noted that currents inside a conductor cannot be recovered from external measurements; radial dipoles in MEG (Theorem 1.9) are an explicit silent source. The singular values of a realistic lead field decay over several orders of magnitude, so regularisation is unavoidable.

4.2 Dipole fitting

The reduced cost has local minima and is minimised from a grid of starting points. The goodness of fit 1 - \|P^\perp M\|^2/\|M\|^2 is the explained variance. Dipole fitting suits focal responses, such as early auditory components.

4.3 MUSIC

Assume k sources with fixed positions, topographies G_S = [g_1, \dots, g_k] and time courses of covariance C_Q, in white noise (after whitening). The data covariance is R = G_S C_Q G_S^\top + \sigma^2 I_m .

In practice R is estimated and c peaks near, not exactly at, the sources. RAP-MUSIC (Mosher and Leahy, 1999) projects out each source found before scanning again, which avoids picking secondary peaks of the same source. If two sources are perfectly correlated, C_Q is singular and the theorem fails.

4.4 Beamforming

The noise term 1/\alpha grows with depth, because deep sources have small topographies. The neural activity index divides P(r) by 1/(g^\top C_N^{-1} g) to remove it. Beamformers assume uncorrelated sources: two perfectly correlated sources cancel each other (Exercise 10).

4.5 Distributed imaging: minimum-norm estimates

With R = I, the estimate favours superficial sources, whose topographies have a large norm. Depth weighting uses R = \mathrm{diag}(\|g_j\|^{-2\gamma}) with \gamma between 0.5 and 0.8. The regularisation \lambda is often set to 1/\mathrm{SNR}^2 on whitened data.

The theorem says nothing about several simultaneous sources or about noise, where sLORETA can be biased; the unnormalised MNE is biased even in this ideal case. Sparse priors (minimum current or mixed-norm estimates) are treated in section 9.2.

4.6 Validation

For MNE with R = I, \mathcal R = G^\top T G is symmetric: point-spread and cross-talk functions coincide. Validation also uses simulations at realistic SNR, phantoms, known functional anatomy, intracranial recordings or fMRI of the same task, and test–retest reliability.

import numpy as np
rng = np.random.default_rng(1)

def fib_sphere(n, radius, zmin=0.0):
    """Quasi-uniform points on the upper part (z > zmin) of a sphere."""
    k = np.arange(n) + 0.5
    z = 1 - k / n * (1 - zmin)
    phi = np.pi * (1 + 5**0.5) * k
    rho = np.sqrt(1 - z**2)
    return radius * np.c_[rho * np.cos(phi), rho * np.sin(phi), z]

# Sources: radial dipoles on two shells (superficial 7 cm, deep 5 cm); 60 EEG electrodes at 9 cm
src = np.vstack([fib_sphere(300, 0.07, 0.2), fib_sphere(150, 0.05, 0.2)])
ori = src / np.linalg.norm(src, axis=1, keepdims=True)
elec = fib_sphere(60, 0.09, 0.0)
sigma = 0.33
d = elec[:, None, :] - src[None, :, :]                       # (60, 450, 3)
G = np.einsum("ijk,jk->ij", d, ori) / (4 * np.pi * sigma * np.linalg.norm(d, axis=2) ** 3)
H = np.eye(60) - 1 / 60                                      # average reference
G = H @ G                                                     # V/(A m)

def simulate(k, snr=1.0, n_t=200):
    s = 10e-9 * np.sin(2 * np.pi * 10 * np.arange(n_t) / 1000)   # 10 nAm, 10 Hz
    M0 = np.outer(G[:, k], s)
    noise_sd = np.linalg.norm(M0) / np.sqrt(M0.size) / snr
    return M0 + H @ (noise_sd * rng.standard_normal(M0.shape)), noise_sd**2 * H

def localise(M, Cn, lam=1 / 9):
    Cn_r = Cn + 1e-6 * np.trace(Cn) / 60 * np.eye(60)        # regularised noise covariance
    T = np.linalg.inv(G @ G.T + lam * np.trace(G @ G.T) / np.trace(Cn_r) * Cn_r)
    W = G.T @ T                                               # minimum-norm operator
    J = W @ M
    p_mne = np.sum(J**2, axis=1)
    p_dspm = p_mne / np.einsum("ij,jk,ik->i", W, Cn_r, W)     # noise-normalised
    p_slor = p_mne / np.einsum("ij,ji->i", W, G)              # diag of the resolution matrix
    C = M @ M.T / M.shape[1] + 0.05 * np.trace(M @ M.T) / M.shape[1] / 60 * np.eye(60)
    Ci = np.linalg.inv(C)
    p_lcmv = 1 / np.einsum("ij,jk,ki->i", G.T, Ci, G)         # unit-gain beamformer power
    p_lcmv /= np.einsum("ij,jk,ki->i", G.T, np.linalg.inv(Cn_r), G) ** -1  # neural activity index
    U = np.linalg.svd(M, full_matrices=False)[0][:, :1]       # signal subspace, rank 1
    p_music = np.sum((U.T @ G) ** 2, axis=0) / np.sum(G**2, axis=0)
    return {name: np.argmax(p) for name, p in
            [("MNE", p_mne), ("dSPM", p_dspm), ("sLORETA", p_slor),
             ("LCMV", p_lcmv), ("MUSIC", p_music)]}

for depth, ks in [("superficial", range(0, 300, 30)), ("deep", range(300, 450, 15))]:
    errs = {}
    for k in ks:
        M, Cn = simulate(k)
        for name, j in localise(M, Cn).items():
            errs.setdefault(name, []).append(1000 * np.linalg.norm(src[j] - src[k]))
    print(f"{depth:11s} mean localisation error (mm): " +
          ", ".join(f"{n} {np.mean(e):5.1f}" for n, e in errs.items()))

5. fMRI: the BOLD signal and the general linear model

In which of the hundred thousand voxels of an fMRI run does the signal follow the task, and is the effect present in the population, not only in the subjects scanned? The general linear model and its two-level extension answer.

5.1 The BOLD signal

Deoxyhaemoglobin is paramagnetic and oxyhaemoglobin is not. When a brain area becomes active, blood flow increases more than oxygen consumption, so the local concentration of deoxyhaemoglobin falls. The field inhomogeneities it creates decrease, the transverse relaxation time T_2^* lengthens, and the MR signal rises by 1–5 % at 3 T: the BOLD contrast (Ogawa et al., 1990).

The response is slow: it peaks about 5 s after a brief stimulus, then undershoots. It correlates best with local field potentials, hence with synaptic activity (Logothetis et al., 2001). A volume of 2–3 mm voxels is acquired every 1–2 s (the repetition time, TR), corrected for motion, registered to a template and often smoothed.

5.2 The linear model of the response

Linearity holds well for stimuli a few seconds apart; closer stimuli saturate, which nonlinear models such as the balloon model describe (Buxton et al., 1998). Adding the time derivative of h to the design absorbs small shifts of latency.

5.3 Estimation and inference

fMRI noise is autocorrelated; an AR(1) model with coefficient \rho has a bidiagonal whitener, with first row \sqrt{1-\rho^2}, then rows (-\rho, 1). Ignoring the autocorrelation inflates false positives (Code 5.6). The efficiency of a design for c is 1/\big(c^\top (X^\top X)^{-1} c\big): blocks of 15–20 s maximise it; jittered events estimate the response shape better.

5.4 Mixed-effects models for group analysis

This summary-statistics approach (Holmes and Friston, 1998) is robust when within-subject variances are similar. In general, the linear mixed model y = X\beta + Z b + \varepsilon, with b \sim \mathcal{N}(0, \Gamma), has covariance V = Z \Gamma Z^\top + \sigma^2 I; its fixed effects follow by GLS (Proposition 5.3), its variance components by restricted maximum likelihood.

import numpy as np
from scipy.stats import gamma, t as student

rng = np.random.default_rng(2)
TR, n = 2.0, 200                                   # 200 scans, 6 min 40 s

def hrf(tt):                                       # canonical double-gamma (SPM shape)
    return gamma.pdf(tt, 6) - gamma.pdf(tt, 16) / 6

dt = 0.1
tt = np.arange(0, 32, dt)
u = (np.arange(0, n * TR, dt) % 40 < 20).astype(float)   # 20 s on / 20 s off
reg = np.convolve(u, hrf(tt))[: u.size] * dt
reg = reg[:: int(TR / dt)]                         # sample at the scan times
k = np.arange(1, 5)                                # discrete cosine drifts, period > 128 s
time = np.arange(n) * TR
drift = np.cos(np.pi * np.outer(time / (n * TR), k))
X = np.column_stack([reg, drift, np.ones(n)])
c = np.r_[1.0, np.zeros(X.shape[1] - 1)]           # contrast: effect of the task

def ar1_noise(rho, sd, size):
    e = np.zeros(size)
    e[0] = rng.normal(0, sd / np.sqrt(1 - rho**2))
    for i in range(1, size):
        e[i] = rho * e[i - 1] + rng.normal(0, sd)
    return e

def glm_t(y, X, rho=0.0):
    """OLS (rho = 0) or GLS after AR(1) prewhitening with a known rho."""
    W = np.eye(len(y)) - rho * np.eye(len(y), k=-1)
    W[0, 0] = np.sqrt(1 - rho**2)
    Xw, yw = W @ X, W @ y
    beta, *_ = np.linalg.lstsq(Xw, yw, rcond=None)
    res = yw - Xw @ beta
    dof = len(y) - X.shape[1]
    s2 = res @ res / dof
    var_c = s2 * c @ np.linalg.inv(Xw.T @ Xw) @ c
    return c @ beta, (c @ beta) / np.sqrt(var_c), dof

# One voxel: 1% signal change on a baseline of 100, AR(1) noise with rho = 0.4
y = 100 + 1.0 * reg / reg.max() + 0.5 * drift[:, 0] + ar1_noise(0.4, 0.8, n)
for name, rho in [("OLS", 0.0), ("GLS, AR(1)", 0.4)]:
    b, tval, dof = glm_t(y, X, rho)
    print(f"{name:11s} beta = {b:.3f}, t = {tval:.2f}, p = {student.sf(tval, dof):.1e}")

# False-positive rate of OLS under AR(1) noise and no activation (nominal 5 %)
fp = [np.mean([glm_t(100 + ar1_noise(0.4, 0.8, n), X, r)[1] > student.isf(0.05, n - X.shape[1])
               for _ in range(500)]) for r in (0.0, 0.4)]
print(f"null voxels declared active at p < 0.05: OLS {fp[0]:.3f}, GLS {fp[1]:.3f}")

# Group: 12 subjects, true mean effect 0.4, between-subject sd 0.5
betas, var_b = [], []
for s in range(12):
    effect = 0.4 + 0.5 * rng.standard_normal()
    ys = 100 + effect * reg / reg.max() + ar1_noise(0.4, 0.8, n)
    b, tval, dof = glm_t(ys, X, 0.4)
    betas.append(b); var_b.append((b / tval) ** 2)
betas, var_b = np.array(betas), np.array(var_b)
t_fixed = betas.sum() / np.sqrt(var_b.sum())       # ignores between-subject variance
t_random = betas.mean() / (betas.std(ddof=1) / np.sqrt(12))
print(f"group: fixed-effects t = {t_fixed:.2f}, random-effects t = {t_random:.2f} "
      f"(p = {student.sf(t_random, 11):.3f}, 11 dof)")

6. Hypothesis testing and multiple comparisons

A statistical map tests one hypothesis per voxel, or per sensor and time point: tens of thousands of tests at once. This chapter makes precise which error rate is controlled, and gives tests that remain valid without distributional assumptions.

6.1 Mass-univariate testing

At an uncorrected threshold of p < 0.001, a map of 50 000 voxels under the global null contains 50 false positives on average: the dead salmon of the introduction.

6.2 Permutation tests

For a group analysis, the subject effect maps are symmetric around zero under H_0, so flipping their signs leaves the distribution invariant: 2^S transformations (Nichols and Holmes, 2002). When B is too large, a random sample of B - 1 transformations plus the identity keeps the test valid.

6.3 Controlling the familywise error rate

The maximum statistic adapts to the correlation between tests: on a smooth map it gives a lower threshold than Bonferroni. Control extends to any configuration of true nulls under the subset pivotality condition (Westfall and Young, 1993). Random field theory is the parametric alternative (section 9.3).

6.4 The false discovery rate

The bound also holds under positive regression dependence (Benjamini and Yekutieli, 2001), a reasonable model for smooth maps. FDR is an average: on a given map the proportion of false discoveries can exceed q.

6.5 Cluster-level inference

A significant cluster means that at least one of its voxels is active. The result depends on u, which threshold-free cluster enhancement avoids (Smith and Nichols, 2009). Eklund et al. (2016) found parametric cluster inference invalid at liberal thresholds, and permutation valid. For M/EEG, clusters span sensors, times and frequencies (Maris and Oostenveld, 2007).

import numpy as np
from scipy.ndimage import gaussian_filter
from scipy.stats import t as student

rng = np.random.default_rng(4)
n_sub, shape = 16, (40, 40)                        # 16 subjects, 1600 "voxels"
truth = np.zeros(shape, bool); truth[10:18, 10:18] = True   # 64 active voxels
maps = np.array([gaussian_filter(rng.standard_normal(shape), 1.5) for _ in range(n_sub)])
maps /= maps.std()
maps[:, truth] += 1.0                              # effect of 1 noise sd

def tmap(Y):                                       # one-sample t over the first axis
    return Y.mean(0) / (Y.std(0, ddof=1) / np.sqrt(len(Y)))

t_obs = tmap(maps).ravel()
dof, V = n_sub - 1, t_obs.size

# Sign-flipping permutation distribution of the maximum statistic
flat = maps.reshape(n_sub, -1)
max_null = np.array([tmap(rng.choice([-1, 1], size=(n_sub, 1)) * flat).max()
                     for _ in range(2000)])
thr_maxT = np.quantile(max_null, 0.95)

def bh(p, q=0.05):                                 # Benjamini-Hochberg step-up
    order = np.argsort(p)
    ok = p[order] <= q * np.arange(1, p.size + 1) / p.size
    k = np.max(np.nonzero(ok)[0]) + 1 if ok.any() else 0
    out = np.zeros(p.size, bool); out[order[:k]] = True
    return out

p = student.sf(t_obs, dof)
rules = {"uncorrected p<0.05": p < 0.05,
         "uncorrected p<0.001": p < 0.001,
         "Bonferroni 0.05": p < 0.05 / V,
         f"max-T permutation (t>{thr_maxT:.2f})": t_obs > thr_maxT,
         "BH-FDR q=0.05": bh(p)}
for name, rej in rules.items():
    tp, fp = np.sum(rej & truth.ravel()), np.sum(rej & ~truth.ravel())
    print(f"{name:32s} discoveries {rej.sum():4d}  false {fp:3d}  true {tp:2d}/64")
print(f"Bonferroni t threshold: {student.isf(0.05 / V, dof):.2f}")

7. Connectivity, networks and covariance learning

Connectivity analysis asks which brain regions interact. With a few hundred time points and up to hundreds of regions, the main difficulty is statistical: estimating a large covariance matrix and its inverse.

7.1 Functional and effective connectivity

The usual pipeline divides the brain into p regions with an atlas, extracts one time series per region, X \in \mathbb{R}^{n \times p}, and estimates a p \times p matrix, read as a weighted graph whose nodes are the regions.

7.2 Partial correlation and Gaussian graphical models

The sparsity pattern of K is the graph of the Gaussian graphical model. In a chain X_1 \to X_2 \to X_3, X_1 and X_3 are correlated, but their partial correlation is zero: partial correlations separate direct from indirect links (Smith et al., 2011).

7.3 Shrinkage estimation of covariance

The sample covariance S is singular when n \le p and ill-conditioned when n is a few times p: its large eigenvalues are biased upwards and its small ones downwards.

Ledoit and Wolf (2004) estimate a and b consistently from the data, which yields a well-conditioned estimate for any n. Shrinkage is also the standard regulariser of noise covariances (ch. 3) and of LDA (ch. 8).

7.4 Sparse inverse covariance: the graphical lasso

Decreasing \lambda from this value adds edges one by one, which traces a path of graphs. For a group of subjects, a penalty shared across subjects (group sparsity) estimates a common support with subject-specific values (Varoquaux et al., 2010).

import numpy as np
from sklearn.covariance import LedoitWolf, GraphicalLassoCV, empirical_covariance

rng = np.random.default_rng(5)
p, n = 30, 60                                      # 30 regions, 60 time points
K = np.eye(p)                                      # sparse precision: a ring + 5 chords
for i in range(p):
    K[i, (i + 1) % p] = K[(i + 1) % p, i] = -0.4
for i, j in [(0, 15), (3, 20), (5, 12), (8, 25), (10, 27)]:
    K[i, j] = K[j, i] = -0.3
Sigma = np.linalg.inv(K)
X = rng.multivariate_normal(np.zeros(p), Sigma, size=n)
edges = np.abs(np.triu(K, 1)) > 0

def partial_corr(P):
    d = 1 / np.sqrt(np.diag(P))
    return -P * np.outer(d, d)

S = empirical_covariance(X)
lw = LedoitWolf().fit(X)
gl = GraphicalLassoCV(cv=5).fit(X)
print(f"condition number: sample {np.linalg.cond(S):.0f}, "
      f"Ledoit-Wolf {np.linalg.cond(lw.covariance_):.0f} (shrinkage {lw.shrinkage_:.2f})")
for name, P in [("inverse of S", np.linalg.inv(S)),
                ("inverse of Ledoit-Wolf", lw.precision_),
                ("graphical lasso", gl.precision_)]:
    err = np.linalg.norm(P - K) / np.linalg.norm(K)
    pc = np.abs(np.triu(partial_corr(P), 1))
    top = pc >= np.sort(pc[np.triu_indices(p, 1)])[-edges.sum()]   # keep as many as true edges
    print(f"{name:23s} rel. error {err:.2f}, true edges among the top-{edges.sum()}: "
          f"{np.sum(top & edges)}/{edges.sum()}")
nz = np.abs(np.triu(gl.precision_, 1)) > 1e-8
print(f"graphical lasso (alpha = {gl.alpha_:.3f}): {nz.sum()} edges selected, "
      f"{np.sum(nz & edges)} true, {np.sum(nz & ~edges)} spurious")

7.5 Covariance matrices as points of a manifold

Distances between covariances are thus unchanged by any invertible mixing of the channels. Tangent-space embeddings turn connectomes into vectors for standard classifiers and predict phenotypes from resting-state fMRI well (Dadi et al., 2019). The same geometry drives the Riemannian classifiers of ch. 8.

8. Brain–computer interfaces

A brain–computer interface decodes, trial after trial, a few bits of intention from noisy signals and shows the result to the user, who learns in turn. This chapter assembles the previous tools into such a loop.

8.1 Communicating via brain activity

Non-invasive BCIs mostly use EEG and three paradigms:

• Motor imagery: imagining a movement of the right hand reduces the power of the mu (8–12 Hz) and beta (13–30 Hz) rhythms over the left sensorimotor cortex, an event-related desynchronisation (Pfurtscheller and Lopes da Silva, 1999).
• P300: a rare, attended stimulus evokes a positive deflection about 300 ms later; the P300 speller flashes the rows and columns of a 6 \times 6 letter matrix (Farwell and Donchin, 1988).
• SSVEP: targets flickering at distinct frequencies evoke responses at those frequencies over the visual cortex.

8.2 Feature extraction

For evoked potentials, features are the decimated amplitudes of a post-stimulus window, possibly after spatial filtering. For rhythms, they are band powers after a spatial filter that the following theorem optimises.

The features are the log-variances of the 2 or 3 filtered signals at each end of the spectrum (Ramoser et al., 2000). CSP overfits with few trials and is sensitive to outliers. The Riemannian alternative uses the trial covariances themselves (Definition 7.7).

8.3 Classification

For evoked-potential features of dimension in the hundreds, \Sigma is replaced by its Ledoit–Wolf estimate (Proposition 7.3; Blankertz et al., 2011). The minimum distance to the mean classifier assigns a trial covariance to the class whose geometric mean is closest (Definition 7.7); it needs no spatial filter (Barachant et al., 2012).

8.4 Model selection and evaluation

Band, time window, number of filters and shrinkage are hyperparameters, chosen by nested cross-validation. Every data-dependent step, CSP and ICA included, must be fitted inside the training folds; fitted on all trials, CSP separates even random labels. Folds should respect blocks and sessions, since EEG is non-stationary.

Cross-validated accuracies come from overlapping training sets, so a permutation test of the whole pipeline (Theorem 6.2) is the safer assessment (Combrisson and Jerbi, 2015).

import numpy as np
from scipy.linalg import eigh
from scipy.signal import butter, sosfiltfilt
from scipy.stats import binom
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(6)
fs, n_ch, n_tr, n_s = 128, 16, 80, 2 * 128        # 80 trials of 2 s, 16 channels
sos = butter(4, [8, 13], btype="band", fs=fs, output="sos")

def mu_rhythm(size):                              # band-limited (8-13 Hz) source activity
    return sosfiltfilt(sos, rng.standard_normal(size), axis=-1)

A = rng.standard_normal((n_ch, n_ch))             # mixing: sources -> channels
y = np.repeat([0, 1], n_tr // 2)                  # 0 = left hand, 1 = right hand imagery
X = np.empty((n_tr, n_ch, n_s))
for i, c in enumerate(y):
    S = mu_rhythm((n_ch, n_s))
    S[0] *= 0.5 if c == 1 else 1.0                # ERD of source 0 during right-hand imagery
    S[1] *= 0.5 if c == 0 else 1.0                # ERD of source 1 during left-hand imagery
    X[i] = A @ S + 0.5 * rng.standard_normal((n_ch, n_s))
X = sosfiltfilt(sos, X, axis=-1)                  # band-pass, as in any MI pipeline

class CSP(BaseEstimator, TransformerMixin):
    """Common spatial patterns: generalised eigenvectors of (C1, C0 + C1)."""
    def __init__(self, n_filters=4):
        self.n_filters = n_filters
    def fit(self, X, y):
        C = [np.mean([x @ x.T / x.shape[1] for x in X[y == k]], axis=0) for k in (0, 1)]
        vals, W = eigh(C[1], C[0] + C[1])         # ascending eigenvalues in (0, 1)
        idx = np.r_[np.arange(self.n_filters // 2), np.arange(-(self.n_filters // 2), 0)]
        self.W_ = W[:, idx]
        return self
    def transform(self, X):
        Z = np.einsum("ck,ncs->nks", self.W_, X)
        return np.log(np.var(Z, axis=2))          # log-variance features

class LogVar(BaseEstimator, TransformerMixin):    # baseline: log-variance of raw channels
    def fit(self, X, y): return self
    def transform(self, X): return np.log(np.var(X, axis=2))

cv = StratifiedKFold(5, shuffle=True, random_state=0)
for name, pipe in [("log-variance + LDA", make_pipeline(LogVar(), LDA())),
                   ("CSP(4) + LDA", make_pipeline(CSP(4), LDA()))]:
    acc = cross_val_score(pipe, X, y, cv=cv)
    print(f"{name:20s} accuracy {acc.mean():.3f} +/- {acc.std():.3f}")

# Leakage: CSP fitted on all trials with random labels, then cross-validated LDA
y_rand = rng.permutation(y)
Z = CSP(4).fit(X, y_rand).transform(X)
print(f"random labels, CSP outside CV: {cross_val_score(LDA(), Z, y_rand, cv=cv).mean():.3f}")
print(f"random labels, CSP inside CV:  "
      f"{cross_val_score(make_pipeline(CSP(4), LDA()), X, y_rand, cv=cv).mean():.3f}")
k = binom.isf(0.05, n_tr, 0.5) + 1                # smallest k with P(Bin >= k) <= 0.05
print(f"chance level at p = 0.05 with {n_tr} trials: {k / n_tr:.3f}")

8.5 The role of feedback

A BCI is a closed loop: the user sees the decoded output and adapts to it. Controlling sensorimotor rhythms is a skill learned through feedback, so two learners interact, the user and the classifier.

Feedback changes the signals: features recorded during calibration without feedback differ from those of online use. Between 15 and 30 % of users do not reach useful control (Vidaurre and Blankertz, 2010), and standard training protocols have been criticised (Jeunet et al., 2016). Neurofeedback uses the same loop for therapy.

Running averages thus track slow drifts without labels (Vidaurre et al., 2011). Adapting faster than the user learns can destabilise the loop.

9. Complements: boundary elements, sparse imaging, random fields, decoding, effective connectivity

This self-study chapter completes the course chapters with six topics that the official syllabus names or implies: the surface-integral form of the forward problem and Maxwell filtering, sparse source imaging, the parametric theory of thresholds, the interpretation of decoders, M/EEG coherence under volume conduction, and directed connectivity.

9.1 Boundary elements, the Geselowitz formula and Maxwell filtering

The MEG forward problem therefore only needs V on the interfaces. Green’s representation gives, for the potential itself, a second-kind integral equation on the interfaces (Geselowitz, 1967). Discretising it on triangulated surfaces is the boundary element method; the symmetric formulation of Kybic et al. (2005), implemented in OpenMEEG (Gramfort et al., 2010), is accurate even for the thin skull compartment.

9.2 Sparse source imaging with mixed norms

Proximal gradient descent alternates a gradient step on the data term with the block soft-thresholding of item 1; FISTA accelerates it to a rate of 1/k^2 in objective (Beck and Teboulle, 2009).

import numpy as np
rng = np.random.default_rng(7)

def fib_sphere(n, radius, zmin=0.0):
    k = np.arange(n) + 0.5
    z = 1 - k / n * (1 - zmin)
    phi = np.pi * (1 + 5**0.5) * k
    return radius * np.c_[np.sqrt(1 - z**2) * np.cos(phi), np.sqrt(1 - z**2) * np.sin(phi), z]

src, elec = fib_sphere(400, 0.07, 0.2), fib_sphere(60, 0.09)
ori = src / np.linalg.norm(src, axis=1, keepdims=True)
d = elec[:, None] - src[None]
G = np.einsum("ijk,jk->ij", d, ori) / (4 * np.pi * 0.33 * np.linalg.norm(d, axis=2) ** 3)
G = (np.eye(60) - 1 / 60) @ G
G /= np.linalg.norm(G, axis=0)                     # normalised columns

t = np.arange(100) / 500
active = [40, 250]                                  # two sources, different time courses
J_true = np.zeros((400, 100))
J_true[active[0]] = np.sin(2 * np.pi * 10 * t)
J_true[active[1]] = np.exp(-((t - 0.1) / 0.03) ** 2)
M = G @ J_true
M += 0.1 * np.linalg.norm(M) / np.sqrt(M.size) * rng.standard_normal(M.shape)

def mxne(M, G, lam, n_iter=500):
    """FISTA for min 1/2 ||M - G J||_F^2 + lam * sum_j ||J_j.||_2."""
    L = np.linalg.norm(G, 2) ** 2                   # Lipschitz constant of the gradient
    J = Z = np.zeros((G.shape[1], M.shape[1])); tk = 1.0
    for _ in range(n_iter):
        Y = Z - G.T @ (G @ Z - M) / L
        norms = np.linalg.norm(Y, axis=1, keepdims=True)
        J_new = np.maximum(0, 1 - lam / L / np.maximum(norms, 1e-12)) * Y   # block soft-threshold
        tk_new = (1 + np.sqrt(1 + 4 * tk**2)) / 2
        Z = J_new + (tk - 1) / tk_new * (J_new - J)
        J, tk = J_new, tk_new
    return J

lam_max = np.max(np.linalg.norm(G.T @ M, axis=1))   # above it, J = 0
J = mxne(M, G, 0.3 * lam_max)
support = np.nonzero(np.linalg.norm(J, axis=1) > 1e-8)[0]
print("lambda_max =", round(lam_max, 2), "| selected sources:", support.tolist())
for k in active:
    j = support[np.argmin(np.linalg.norm(src[support] - src[k], axis=1))]
    print(f"true source {k}: nearest selected {j}, distance "
          f"{1000 * np.linalg.norm(src[j] - src[k]):.1f} mm")
J_mne = G.T @ np.linalg.solve(G @ G.T + 0.1 * np.eye(60), M)
energy = np.linalg.norm(J_mne, axis=1) ** 2
print(f"MNE: sources holding 90% of the energy: "
      f"{np.searchsorted(np.cumsum(np.sort(energy)[::-1]) / energy.sum(), 0.9) + 1}")

9.3 Random field theory

9.4 Decoding and the interpretation of linear weights

Decoding predicts the stimulus or the mental state from a distributed pattern of activity (Haxby et al., 2001), over the whole brain or in moving spheres called searchlights (Kriegeskorte, Goebel and Bandettini, 2006). Linear models with \ell_2 or sparse penalties dominate. Their accuracy must be estimated by cross-validation across runs or subjects, whose variance is large on typical sample sizes (Varoquaux et al., 2017).

9.5 M/EEG connectivity and volume conduction

Coherence between sensors, or between sources estimated by a linear inverse (ch. 4), is inflated by the spatial spread of the fields; the imaginary part of coherency and related lagged measures are insensitive to it, at the price of missing true zero-lag interactions.

9.6 Effective connectivity: Granger causality and dynamic causal models

On fMRI, the haemodynamic delay varies between regions by more than the neural delays, which makes lag-based methods unreliable on BOLD signals (Smith et al., 2011). On M/EEG, Granger measures are computed per frequency and must be protected against volume conduction (section 9.5).

10. Research directions and open problems

This self-study chapter surveys active questions around the course, with entry points into the literature.

10.1 Reproducibility and the size of studies

• Analytic flexibility. Botvinik-Nezer et al. (2020) gave the same fMRI data set to 70 teams to test nine hypotheses; the teams’ conclusions differed substantially, mostly because of analysis choices. Standardised preprocessing such as fMRIPrep (Esteban et al., 2019) and the BIDS data standard (Gorgolewski et al., 2016) reduce this variability without removing it.
• Sample size. Marek et al. (2022) showed that associations between brain measures and individual traits are small and need thousands of participants to replicate. Open question: which designs and which within-subject measures give reliable individual-level biomarkers with fewer subjects?

10.2 Prediction from large cohorts

• Connectomes. Dadi et al. (2019) benchmarked pipelines that predict phenotypes from resting-state fMRI and found that tangent-space covariance embeddings (section 7.5) with linear models perform well across cohorts.
• Without source modelling. Sabbagh et al. (2019) showed that Riemannian regression on sensor covariances predicts targets from M/EEG without inverse modelling; Engemann et al. (2020) combined MEG and MRI to improve the prediction of brain age. Open question: how to transfer such models across scanners, sites and populations.

10.3 Deep learning for non-invasive signals

• Architectures. Convolutional networks designed for EEG decode motor imagery and evoked responses (Schirrmeister et al., 2017; Lawhern et al., 2018), with gains over CSP and Riemannian pipelines that depend strongly on the amount of data.
• Decoding language. Défossez et al. (2023) trained a model with a contrastive objective to identify, from a few seconds of MEG or EEG, which speech segment a participant was hearing.
• Open problems. Generalisation across subjects and sessions, self-supervised pretraining on large unlabelled recordings, and interpretation of the learned filters in the sense of Proposition 9.7.

10.4 Transfer learning and calibration-free BCIs

• Calibration sessions are tiring for patients. Re-centring each session’s covariances at their own mean before tangent-space classification aligns sessions and subjects (Zanini et al., 2018).
• The MOABB benchmark (Jayaram and Barachant, 2018) compares BCI pipelines on many public data sets. Open question: which representation of the signals transfers best between users with different anatomies and strategies?

10.5 The user in the loop

• A sizeable fraction of users fail to control a BCI (section 8.5). Corsi et al. (2020) related training performance to the functional network organisation of the brain measured with M/EEG; Jeunet et al. (2016) questioned the instructional design of standard protocols.
• Open questions: which feedback (modality, timing, adaptation rate) speeds up learning, and how to model the joint dynamics of user and classifier.

10.6 Invasive interfaces as a benchmark

• Intracortical and electrocorticographic interfaces reach far higher rates: decoding attempted handwriting (Willett et al., 2021) or speech (Moses et al., 2021) in paralysed participants, or restoring walking through a brain–spine interface (Lorach et al., 2023).
• These results set a target for non-invasive BCIs and raise questions shared by both: long-term stability, adaptation and the protection of neural data.

10.7 Further open problems

• Individual head models: estimating skull conductivity in vivo, which dominates EEG forward errors (ch. 2).
• Uncertainty quantification for source estimates, beyond point estimates such as those of ch. 4.
• Fusion of M/EEG and fMRI, which have complementary resolutions but different physiological origins.
• High-field (7 T) fMRI resolving cortical layers, which calls for new models of the BOLD response across depth.

Exercises

1. ★ Orders of magnitude. A dipole of 20 nA·m sits in a medium of conductivity 0.33 S/m. (a) Compute the potential 3 cm away along its axis (Proposition 1.5). (b) Compute the primary magnetic field 3 cm away, perpendicular to the dipole (Proposition 1.8). (c) Compare it with the Earth’s field of 50 µT.

   Solution: (a) and (b) give V = \frac{2 \times 10^{-8}}{4\pi \times 0.33 \times 0.03^2} \approx 5.4 \ \mu\text{V}, \qquad B_0 = \frac{10^{-7} \times 2 \times 10^{-8}}{0.03^2} \approx 2.2 \text{ pT}.

   3. The ratio is \frac{2.2 \times 10^{-12}}{5 \times 10^{-5}} \approx 4 \times 10^{-8}, which is why MEG needs a shielded room and gradiometers.

2. ★ How many trials? A P300 of 5 µV is recorded in background EEG of standard deviation 20 µV per trial. (a) How many trials give an SNR of 2 after averaging? (b) The response is a Gaussian bump of standard deviation 50 ms and its latency jitters with a Gaussian of standard deviation 50 ms. By what factor is the averaged peak reduced, and how many trials are then needed?

   Solution: (a) 5\sqrt N/20 \ge 2 gives N \ge 64 (Proposition 3.2). (b) The average is the bump convolved with the jitter density: a Gaussian of standard deviation \sqrt{50^2 + 50^2} \approx 71 ms, with the same area, so its peak is multiplied by 50/71 \approx 0.71. The condition becomes 0.71 \times 5\sqrt N / 20 \ge 2, so N \ge 128 .

3. ★ A minimum-norm estimate by hand. Two sensors and three sources have the lead field G = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}, with C_N = I, R = I and \lambda = 1. (a) Give a silent source. (b) Compute W (Theorem 4.6). (c) Compute the MNE and sLORETA maps for the data (1, 1) and (1, 0), and check Theorem 4.8.

   Solution: (a) G(1, 1, -1)^\top = 0. (b) GG^\top + I = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}, so T = \frac18 \begin{pmatrix} 3 & -1 \\ -1 & 3 \end{pmatrix} and W = G^\top T. The diagonal of WG is (3/8, 3/8, 1/2).

   3. For M = (1,1), the third source: \hat J = (0.25, 0.25, 0.5) and sLORETA (0.167, 0.167, 0.5). For M = (1, 0), the first source: \hat J = (0.375, -0.125, 0.25) and sLORETA (0.375, 0.042, 0.125). Both maps peak at the true source, with amplitudes shrunk by the regularisation.

4. ★ A contrast by hand. Six scans follow a block design with centred regressor x = (1, 1, 1, -1, -1, -1) and a constant; the data are y = (3, 4, 5, 1, 2, 3). Estimate \beta, \hat\sigma^2, the t statistic of the task effect and its one-sided p-value (Theorem 5.2).

   Solution: the columns are orthogonal, so \hat\beta_{\text{task}} = x^\top y / 6 = 1 and \hat\beta_0 = \bar y = 3. The residuals are (-1, 0, 1, -1, 0, 1), so \hat\sigma^2 = 4/4 = 1 with n - p = 4 degrees of freedom. The standard error is \sqrt{1/6} = 0.408, hence t = 2.45, \qquad p = 0.035 \ (t_4) .

5. ★ Bonferroni and Benjamini–Hochberg. Ten tests give the p-values 0.0001, 0.0004, 0.0019, 0.0095, 0.020, 0.028, 0.030, 0.034, 0.046, 0.32. Which are rejected by Bonferroni at FWER 0.05 and by Benjamini–Hochberg at FDR 0.05?

   Solution: Bonferroni compares with 0.005 and rejects the first three. Benjamini–Hochberg compares p_{(i)} with 0.005\, i. The largest such i is 8, since 0.034 is below 0.040 while 0.046 exceeds 0.045: it rejects the eight smallest (Theorem 6.4).

6. ★ Partial correlations of a chain. Let X_1 = Z_1, X_2 = X_1 + Z_2, X_3 = X_2 + Z_3 with independent standard Gaussian variables Z_i. Compute the covariance and the precision matrices, then the correlation and the partial correlation of X_1 and X_3.

   Solution: one checks that \Sigma K = I for \Sigma = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 2 & 2 \\ 1 & 2 & 3 \end{pmatrix}, \qquad K = \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{pmatrix} .

   The correlation of X_1 and X_3 is 1/\sqrt 3 \approx 0.58, but K_{13} = 0, so their partial correlation is 0 (Theorem 7.2): the link goes through X_2. The partial correlations of the two edges are 0.5 and 1/\sqrt 2 \approx 0.71.

7. ★ Information transfer rates. (a) Compute the ITR of a 4-target SSVEP interface with 95 % accuracy and one selection every 2 s. (b) What is the maximal ITR of a two-class BCI with one trial every 4 s? (c) Which accuracy does it need to reach 10 bits/min?

   Solution: (a) B = 2 + 0.95\log_2 0.95 + 0.05\log_2(0.05/3) = 1.63 bits, that is 49 bits/min (Proposition 8.2). (b) At most 1 bit per trial, so 15 bits/min. (c) 10 bits/min means 2/3 bit per trial: 1 - H(P) \ge 2/3, with H the binary entropy, which requires P \ge 0.94 .

8. ★★ MEG sees the tangential part. In a spherically symmetric conductor, write q = q_r \hat r_0 + q_t with q_t \perp r_0. (a) Show that the external field depends only on q_t. (b) A dipole is tilted by 30° from the radial direction. Compare its field with that of a tangential dipole of the same moment, and check with Code 1.11.

   Solution: (a) B is linear in q (Proposition 1.10), and the radial part contributes nothing (Theorem 1.9, item 3), so B(q) = B(q_t). (b) |q_t| = |q| \sin 30° = |q|/2, and q_t has the direction of the tangential dipole, so the field is exactly half of it at every sensor. In Code 1.11, q = 10e-9 * np.array([np.sin(np.pi/6), 0, np.cos(np.pi/6)]) gives |B| = 78.4 fT, half of 156.9 fT.

9. ★★ Reciprocity in an infinite medium. In an infinite homogeneous conductor, a current I enters at point a and leaves at point b. (a) Give the potential U_{ab}. (b) Check Theorem 2.8 for a dipole at r_0, that is, show with Proposition 1.5 that \frac{q \cdot \nabla U_{ab}(r_0)}{I} = V_\infty(a) - V_\infty(b) .

   Solution: (a) Each point source of current I produces I/(4\pi\sigma|r - a|), since the outgoing flux through a sphere is I. Hence U_{ab}(r) = \frac{I}{4\pi\sigma}\Big(\frac{1}{|r - a|} - \frac{1}{|r - b|}\Big) .

   2. Its gradient at r_0 is \nabla U_{ab}(r_0) = \frac{I}{4\pi\sigma}\Big(\frac{a - r_0}{|a - r_0|^3} - \frac{b - r_0}{|b - r_0|^3}\Big), and its dot product with q/I is V_\infty(a) - V_\infty(b).

10. ★★ Correlated sources defeat the beamformer. Two sources have orthonormal topographies g_1, g_2 and power P, in white noise \sigma^2 I. Compute the LCMV output power at g_1 (Theorem 4.4) (a) for uncorrelated sources and (b) for identical time courses. Let \sigma^2 \to 0.

    Solution: (a) C = P g_1 g_1^\top + P g_2 g_2^\top + \sigma^2 I, so g_1^\top C^{-1} g_1 = 1/(P + \sigma^2) and the output power is P + \sigma^2, as in Proposition 4.5.

    2. With u = (g_1 + g_2)/\sqrt 2, C = 2P\, u u^\top + \sigma^2 I, and Sherman–Morrison gives C^{-1} = \frac{1}{\sigma^{2}}\Big(I - \frac{2P}{\sigma^2 + 2P}\, u u^\top\Big), \qquad g_1^\top C^{-1} g_1 = \frac{\sigma^2 + P}{\sigma^2(\sigma^2 + 2P)} . The output power \frac{\sigma^2(\sigma^2 + 2P)}{\sigma^2 + P} \xrightarrow[\sigma \to 0]{} 0 vanishes: the filter uses source 2 to cancel source 1.

11. ★★ Weighted versus unweighted group means. Three subjects give \hat\beta = (1.0, 0.8, 0.2) with first-level variances v = (0.01, 0.01, 0.25), and \tau^2 = 0.04. Compute the weighted estimate of Proposition 5.5 and its standard error, then the plain mean and its standard error under the same model.

    Solution: the weights are 1/(0.05, 0.05, 0.29) = (20, 20, 3.45), of sum 43.45. The weighted mean is 36.69/43.45 = 0.84, with standard error 1/\sqrt{43.45} = 0.15. The plain mean is 0.67, with variance (0.05 + 0.05 + 0.29)/9 = 0.043, that is a standard error of 0.21. Down-weighting the noisy third subject gains precision; the summary-statistics test is exact only when the v_s are equal.

12. ★★ Sign flipping needs symmetry (Python). For 10 subjects with one voxel each, simulate 4000 null data sets and compute the exact sign-flipping p-value of the one-sample t statistic (all 2^{10} flips, Theorem 6.2), in both tails. Compare symmetric noise (Student t_3) with centred exponential noise, of mean zero but skewed.

    Solution:

    import numpy as np
    from itertools import product
    rng = np.random.default_rng(8)
    S = 10
    flips = np.array(list(product([-1, 1], repeat=S)))           # all 1024 sign flips
    
    def tstat(Y):
        return Y.mean(-1) / (Y.std(-1, ddof=1) / np.sqrt(S))
    
    def pvalues(x):                                               # upper- and lower-tail p-values
        t0, tb = tstat(x), tstat(flips * x)
        return np.mean(tb >= t0), np.mean(tb <= t0)
    
    for name, draw in [("Student t(3)", lambda: rng.standard_t(3, S)),
                       ("centred exponential", lambda: rng.exponential(1.0, S) - 1.0)]:
        p = np.array([pvalues(draw()) for _ in range(4000)])
        print(f"{name:20s} P(p <= 0.05): upper tail {np.mean(p[:, 0] <= 0.05):.3f}, "
              f"lower tail {np.mean(p[:, 1] <= 0.05):.3f}")

    Output: with Student noise, the rejection rates at 0.05 are 0.048 (upper tail) and 0.057 (lower tail). The exact rate is 51/1024 \approx 0.050 in each tail, and the Monte Carlo standard error is about 0.0035, so both lie within about two standard errors of it. With skewed noise they are 0.016 and 0.128: flipping signs is not a symmetry of a skewed law, so the exactness of Theorem 6.2 is lost.

13. ★★ The first edge of the graphical lasso. For two variables with S = \begin{pmatrix} 1 & r \\ r & 1 \end{pmatrix} and 0 < \lambda < |r|, solve the graphical lasso with Proposition 7.5. Show that the estimated partial correlation is the soft-thresholded correlation.

    Solution: write W = \hat K^{-1}. The diagonal condition gives W_{ii} = S_{ii} = 1. In dimension 2, \hat K_{12} has the sign opposite to W_{12}, so W_{12} = r - \lambda\, \mathrm{sign}(r).

    For a 2 \times 2 matrix, the partial correlation read from the precision equals the correlation of its inverse, here W_{12}, so \hat\rho_{12} = \mathrm{sign}(r)\,(|r| - \lambda), and it is 0 for \lambda \ge |r|.

14. ★★ CSP with two channels. The class covariances are C_1 = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix} and C_0 = \begin{pmatrix} 3 & -1 \\ -1 & 3 \end{pmatrix}. Find the CSP filters and eigenvalues (Theorem 8.4), and the class variances along each filter.

    Solution: C_0 + C_1 = 6I, so the problem is the eigenproblem of C_1/6: eigenvalues 2/3 and 1/3, with eigenvectors (1, 1) and (1, -1), normalised to (1, \pm 1)/\sqrt{12}.

    Along (1, 1)/\sqrt{12}, class 1 has variance 2/3 and class 0 has 1/3; along (1, -1)/\sqrt{12} the roles swap. The log-variances of the two filtered signals separate the classes.

15. ★★ Leakage and the chance level (Python). In Code 8.7, set n_tr = 40. For 100 random labellings, compute the cross-validated accuracy (a) with CSP fitted on all trials and (b) with CSP inside the pipeline. Compare the 95th percentile of (b) with the binomial threshold of Proposition 8.6.

    Hints and expected result: with CSP outside the cross-validation, the mean accuracy on random labels is about 0.88, and the 95th percentile 0.95. With CSP inside, the mean is 0.48, slightly below 0.5 as is usual for cross-validation under the null, and the 95th percentile is 0.63, below the binomial threshold of 0.65. The true labels give 0.95.

16. ★★★ Two sources and sLORETA (Python). Using the model of Code 4.10 without noise, activate two superficial sources with the same time course and equal amplitudes. Locate the peak of the sLORETA map as a function of the distance between the sources, and explain why Theorem 4.8 does not extend.

    Hints and expected result: take \lambda small (for instance 10^{-3}\, \mathrm{tr}(GG^\top)/60) and draw 300 random pairs. For sources less than 4 cm apart, the peak falls between them in 35–55 % of the pairs, depending on the draw (mean error 4–6 mm); beyond 4 cm it falls on one of them, and the weaker lobe of the other may not be a local maximum.

    With data g_{k_1} + g_{k_2}, the map of Theorem 4.8 becomes s_j \propto \frac{\langle g_j, g_{k_1} + g_{k_2}\rangle_T^2}{\langle g_j, g_j\rangle_T}, maximised by the topography best aligned with the sum, which need not be a source.

17. ★★★ A drifting BCI (Python). Simulate an LDA calibrated on 60 trials, then used online on 400 trials whose features drift slowly. Compare the static classifier with the unsupervised bias adaptation of Proposition 8.8, using a running mean with forgetting factor 0.05. Then discuss how a learning user, whose class separation grows with feedback, interacts with the adaptation rate.

    Solution:

    import numpy as np
    rng = np.random.default_rng(9)
    d, n_cal, n_on = 4, 60, 400
    m0, m1 = np.zeros(d), np.r_[1.2, 0.6, 0, 0]
    cov = np.eye(d)
    def draw(n, shift):
        y = rng.integers(0, 2, n)
        X = np.where(y[:, None] == 1, m1, m0) + rng.standard_normal((n, d)) + shift
        return X, y
    Xc, yc = draw(n_cal, 0.0)
    mu0, mu1 = Xc[yc == 0].mean(0), Xc[yc == 1].mean(0)
    S = np.cov(np.vstack([Xc[yc == 0] - mu0, Xc[yc == 1] - mu1]).T)
    w = np.linalg.solve(S, mu1 - mu0); b = -w @ (mu0 + mu1) / 2
    acc = {"static": [], "bias-adapted": []}
    mbar, eta = (mu0 + mu1) / 2, 0.05
    for t in range(n_on):
        shift = np.r_[3.0, -2.0, 1.0, 0.0] * t / n_on        # slow global drift
        x, y = draw(1, shift); x, y = x[0], y[0]
        acc["static"].append(((w @ x + b) > 0) == y)
        acc["bias-adapted"].append(((w @ x - w @ mbar) > 0) == y)
        mbar = (1 - eta) * mbar + eta * x                     # unsupervised running mean
    for k, v in acc.items():
        v = np.array(v); print(k, "first 100 trials %.3f, last 100 trials %.3f" % (v[:100].mean(), v[-100:].mean()))

    Output: over the first 100 online trials, the two classifiers reach 0.73 and 0.74; over the last 100, the static one falls to 0.64 while the adapted one keeps 0.72, close to the Bayes accuracy \Phi(\Delta/2) \approx 0.75, where \Delta = \|m_1 - m_0\| \approx 1.34 is the distance between the class means (identity covariance). If the user learns, the class means move apart: a fast adaptation of the bias follows the drift but also follows short-term fluctuations, while re-estimating w requires labels or a model of the user. The adaptation rate trades tracking for stability.

18. ★★★ Filters versus patterns (Python). For the CSP filters fitted in Code 8.7, compute the patterns of Proposition 9.7, using the average trial covariance as \Sigma_x. Compare filters and patterns with the true mixing columns A[:, 0] and A[:, 1] of the two sources that carry the class information.

    Hints and expected result: with W the four CSP filters, the patterns are \Sigma_x W (W^\top \Sigma_x W)^{-1}. The best absolute cosine between a pattern and each true mixing column is about 0.97–0.98, while the best cosine for the filters themselves is about 0.5. Filters cancel the other sources and the noise; patterns show where the discriminative sources project, which is what a physiological interpretation needs.

Projects and assessment

The official assessment is two TDs (travaux dirigés, guided exercise work: “L’évaluation du cours est constituée de 2 TD”, master-mva.com); the resit is a critical commentary of a research paper, and the Université Paris-Saclay course detail gives “Evaluation Continue non Intégrale”. No weighting is displayed.

Everything below is a proposal. Each solo subject is sized like one graded TD; the group subjects suit a team TD or a final project. Most subjects study a named paper, which also trains the resit exercise: a commentary that states the question, the method, the evidence, the limits and what would be needed to reproduce it.

Solo projects

1. [Solo · ≈ 20 h] When is sLORETA unbiased? A paper study of Pascual-Marqui (2002). The question: how far does the zero-localisation-error property of Theorem 4.8 survive noise, depth and several sources?
   • Implement: MNE, dSPM and sLORETA from scratch, checked against MNE-Python (Gramfort et al., 2013), on a three-shell spherical EEG model and a volume source space built with MNE-Python, so that no download is needed.
   • Probe: localisation error and spatial dispersion (Definition 4.9) as functions of SNR, depth, \lambda and the number of simultaneous sources (Exercise 16).
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 1, 2, 4.
   • Success: zero error for single noise-free sources reproduced to machine precision; error curves with 95 % intervals over at least 200 simulated sources per condition.
2. [Solo · ≈ 20 h] Riemannian geometry against CSP for motor imagery. A paper study of Barachant, Bonnet, Congedo and Jutten (2012). Does the minimum distance to the Riemannian mean match or beat CSP with LDA, within and across sessions?
   • Data: BCI Competition IV data set 2a (Brunner et al., 2008; 9 subjects, 4 motor-imagery classes, 22 EEG channels, 2 sessions), available through MOABB (Jayaram and Barachant, 2018).
   • Implement: CSP + LDA (Code 8.7), the geometric mean by fixed-point iteration, MDM and tangent space + logistic regression (Definition 7.7); compare with pyRiemann.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters: ch. 7, 8.
   • Success: per-subject accuracies within session (nested cross-validation) and across sessions, with permutation chance levels and a paired test across subjects.
3. [Solo · ≈ 15 h] Cluster failure at small scale. A paper study of Eklund, Nichols and Knutsson (2016). Do parametric and permutation inference control the familywise error on real null data?
   • Data: the development fMRI data set of Richardson et al. (2018), available through Nilearn (Abraham et al., 2014), analysed with a fake block design whose phase is drawn at random for each subject, so that no effect of the movie is shared across subjects.
   • Implement: first-level GLMs (Code 5.6), random groups of 20 subjects, then voxelwise random-field thresholds (section 9.3), sign-flipping maximum statistic and cluster inference at two cluster-forming thresholds.
   • Deliverables: code, a 4-page report, a 10-minute talk. Chapters: ch. 5, 6, 9.
   • Success: estimated FWER of each method over at least 200 random groups, with binomial confidence intervals, and a reasoned comparison with the paper’s findings.

Group projects

1. [Group of 3 · ≈ 30 h each] From raw MEG/EEG to validated source estimates. On the MNE-Python “sample” data set (auditory and visual stimuli; Gramfort et al., 2013), how consistent are the inverse methods of ch. 4, and how much does each preprocessing step matter?
   • Implement: filtering, SSP or ICA against ocular and cardiac artifacts (ch. 3), epoching and averaging; a BEM forward model from the provided surfaces; shrunk noise covariance; MNE, dSPM, sLORETA, LCMV and a dipole fit of the auditory response.
   • Split: preprocessing and denoising; forward and inverse models; validation by resolution matrices and simulations.
   • Deliverables: repository, an 8–10 page report, a 15-minute defence. Chapters: ch. 1–4, 9.
   • Success: auditory responses in the superior temporal cortex and visual responses in the occipital cortex for each method, with simulation-based error bars and the effect of removing each preprocessing step.
2. [Group of 3 · ≈ 30 h each] Connectome-based prediction: children against adults. A paper study of Dadi et al. (2019). Which connectivity estimator predicts age group best from fMRI time series?
   • Data: the development fMRI data set (Richardson et al., 2018; 155 participants, 122 children and 33 adults), through Nilearn.
   • Implement: region time series with two atlases; correlation, partial correlation with Ledoit–Wolf and the graphical lasso, tangent space at the geometric mean (ch. 7); linear classifiers; nested cross-validation; permutation test of the accuracy.
   • Split: time series and confounds; estimators and tangent space; evaluation and statistics.
   • Deliverables: repository, an 8-page report, a 15-minute defence. Chapters: ch. 6, 7, 9.
   • Success: balanced accuracy with confidence intervals for each pipeline, and a test of the paper’s conclusion that tangent-space embeddings perform best.
3. [Group of 2–3 · ≈ 30 h each] Group inference and reproducibility on the Localizer data. How do detection and reproducibility depend on the group model, the correction and the number of subjects?
   • Data: contrast and t maps of the Brainomics/Localizer data set (Pinel et al., 2007; Papadopoulos Orfanos et al., 2017), available through Nilearn; first-level variances follow from (\text{effect}/t)^2.
   • Implement: summary-statistics and mixed-effects group models (Proposition 5.5), sign-flipping maximum statistic, cluster inference and Benjamini–Hochberg (ch. 6); a subsampling study from 10 to 90 subjects, in the spirit of Thirion et al. (2007).
   • Split: group models; inference; reproducibility study.
   • Deliverables: repository, an 8-page report, a 15-minute defence. Chapters: ch. 5, 6.
   • Success: detection counts per method, a check of FWER on sign-randomised nulls, and Dice overlap of thresholded maps between disjoint subgroups as a function of group size.

Written exam (proposal)

Two hours, one handwritten A4 sheet and a calculator allowed; four questions, 20 points.

1. Question (5 points). (a) State the quasistatic equation for the potential and justify the two approximations with orders of magnitude. (b) Prove that a radial dipole in a spherically symmetric conductor produces no magnetic field outside it; you may admit that volume currents do not contribute to the radial component of B.

   Expected:

   • 1. \nabla\cdot(\sigma\nabla V) = \nabla \cdot J^p, from \nabla \times E = 0 and \nabla \cdot J = 0; \omega\varepsilon/\sigma \approx 10^{-3} and a skin depth of tens of metres below 1 kHz (2 pts).
   • 2. When q \parallel r_0 (3 pts), r \cdot B = r \cdot B_0 = -\frac{\mu_0}{4\pi}\, \frac{r \cdot (q \times r_0)}{|r - r_0|^3} = 0 ; outside, B = -\mu_0 \nabla U and U is the integral of the radial component along rays, so U = 0 and B = 0.

2. Question (5 points). Two sensors and three sources have G = \begin{pmatrix} 2 & 1 & 0 \\ 0 & 1 & 2 \end{pmatrix}, with C_N = I. The data are M = (1, 1). (a) Compute the MNE with R = I and \lambda = 1. (b) Compute the sLORETA map. (c) With the data covariance C = g_2 g_2^\top + I, compute the LCMV output power and the neural activity index at each source.

   Expected:

   • 1. The inverse operator uses (1.5 pts) T = \frac{1}{35}\begin{pmatrix} 6 & -1 \\ -1 & 6 \end{pmatrix}, \qquad \hat J = (2/7, 2/7, 2/7): MNE cannot choose, because source 2 has the smallest topography.
   • 2. (WG)_{jj} = (24/35, 10/35, 24/35), so sLORETA gives (0.12, 0.29, 0.12) and peaks at source 2 (1.5 pts).
   • 3. With (2 pts) C^{-1} = \frac13 \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}, the output powers are (0.375, 1.5, 0.375), where 1.5 = 1 + 1/\|g_2\|^2 (Proposition 4.5); the neural activity index is (1.5, 3, 1.5).

3. Question (5 points). (a) Explain the difference between fixed-effects and random-effects group analyses and when the one-sample t test on subject estimates is exact. (b) Ten voxels have p-values 0.00002, 0.0003, 0.0012, 0.004, 0.006, 0.011, 0.04, 0.09, 0.3, 0.7: which are rejected by Bonferroni (FWER 0.05) and by Benjamini–Hochberg (FDR 0.05)? (c) Prove that Bonferroni controls the FWER, and say why the maximum statistic is less conservative on a smooth map.

   Expected:

   • 1. Fixed effects ignore \tau^2 and concern the subjects scanned; random effects use v_s + \tau^2 and concern the population; exact with equal first-level variances (Proposition 5.5) (1.5 pts).
   • 2. Bonferroni (threshold 0.005): the first four. BH: the largest i with p_{(i)} \le 0.005\, i is 6 (0.011 \le 0.030, while 0.04 > 0.035): the first six (2 pts).
   • 3. Union bound over the true nulls (Proposition 6.3); correlated tests have a maximum whose null quantile is below the Bonferroni threshold (1.5 pts).

4. Question (5 points). (a) Show that CSP filters solve a generalised eigenproblem. (b) Give the filters and eigenvalues for C_1 = \mathrm{diag}(4, 1) and C_0 = \mathrm{diag}(1, 4). (c) A two-class BCI reaches 70 % on 50 test trials, one trial every 5 s. Is it above chance at the 5 % level, and what is its information transfer rate?

   Expected:

   • 1. The variance ratio is a generalised Rayleigh quotient; its stationary points satisfy C_1 w = \mu (C_0 + C_1) w (Theorem 8.4, Proposition 0.6) (2 pts).
   • 2. C_0 + C_1 = 5I: filters e_1 and e_2 (normalised by 1/\sqrt 5), with eigenvalues 0.8 and 0.2 (1 pt).
   • 3. The binomial threshold for 50 trials is 32 correct (64 %), so 35 correct is significant. B = 1 + 0.7\log_2 0.7 + 0.3 \log_2 0.3 = 0.12 bit per trial, that is 1.4 bits/min (2 pts).

Grading rubric (proposal)

References

Official references (master-mva.com)

• The course handout (“Le polycopié du cours”), distributed by the instructors.
• R. A. Poldrack, J. A. Mumford, T. E. Nichols, Handbook of Functional MRI Data Analysis, Cambridge University Press, 2011.
• R. Brette, A. Destexhe (eds.), Handbook of Neural Activity Measurement, Cambridge University Press, 2012.

Textbooks

• P. L. Nunez, R. Srinivasan, Electric Fields of the Brain: The Neurophysics of EEG, 2nd ed., Oxford University Press, 2006.
• A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004.
• A. Hyvärinen, J. Karhunen, E. Oja, Independent Component Analysis, Wiley, 2001.
• P. H. Westfall, S. S. Young, Resampling-Based Multiple Testing, Wiley, 1993.
• R. J. Adler, The Geometry of Random Fields, Wiley, 1981.
• J. R. Wolpaw, E. W. Wolpaw (eds.), Brain–Computer Interfaces: Principles and Practice, Oxford University Press, 2012.
• M. Clerc, L. Bougrain, F. Lotte (eds.), Brain–Computer Interfaces 1 and 2, ISTE–Wiley, 2016.

Physics, forward and inverse problems (ch. 1–4, 9)

• H. von Helmholtz, “Ueber einige Gesetze der Vertheilung elektrischer Ströme in körperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche”, Annalen der Physik und Chemie, 1853.
• D. B. Geselowitz, “On bioelectric potentials in an inhomogeneous volume conductor”, Biophysical Journal 7, 1967.
• J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem”, Physics in Medicine and Biology 32, 1987.
• M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Knuutila, O. V. Lounasmaa, “Magnetoencephalography — theory, instrumentation, and applications to noninvasive studies of the working human brain”, Reviews of Modern Physics 65, 1993.
• J. C. Mosher, R. M. Leahy, P. S. Lewis, “EEG and MEG: forward solutions for inverse methods”, IEEE Transactions on Biomedical Engineering 46, 1999.
• J. Kybic, M. Clerc, T. Abboud, O. Faugeras, R. Keriven, T. Papadopoulo, “A common formalism for the integral formulations of the forward EEG problem”, IEEE Transactions on Medical Imaging 24, 2005.
• A. Gramfort, T. Papadopoulo, E. Olivi, M. Clerc, “OpenMEEG: opensource software for quasistatic bioelectromagnetics”, BioMedical Engineering OnLine 9, 2010.
• M. A. Uusitalo, R. J. Ilmoniemi, “Signal-space projection method for separating MEG or EEG into components”, Medical and Biological Engineering and Computing 35, 1997.
• S. Taulu, M. Kajola, “Presentation of electromagnetic multichannel data: the signal space separation method”, Journal of Applied Physics 97, 2005.
• P. Comon, “Independent component analysis, a new concept?”, Signal Processing 36, 1994.
• N. Delfosse, P. Loubaton, “Adaptive blind separation of independent sources: a deflation approach”, Signal Processing 45, 1995.
• A. Hyvärinen, E. Oja, “Independent component analysis: algorithms and applications”, Neural Networks 13, 2000.
• R. O. Schmidt, “Multiple emitter location and signal parameter estimation”, IEEE Transactions on Antennas and Propagation 34, 1986.
• J. C. Mosher, P. S. Lewis, R. M. Leahy, “Multiple dipole modeling and localization from spatio-temporal MEG data”, IEEE Transactions on Biomedical Engineering 39, 1992.
• J. C. Mosher, R. M. Leahy, “Source localization using recursively applied and projected (RAP) MUSIC”, IEEE Transactions on Signal Processing 47, 1999.
• B. D. Van Veen, W. van Drongelen, M. Yuchtman, A. Suzuki, “Localization of brain electrical activity via linearly constrained minimum variance spatial filtering”, IEEE Transactions on Biomedical Engineering 44, 1997.
• M. S. Hämäläinen, R. J. Ilmoniemi, “Interpreting magnetic fields of the brain: minimum norm estimates”, Medical and Biological Engineering and Computing 32, 1994.
• A. M. Dale et al., “Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity”, Neuron 26, 2000.
• R. D. Pascual-Marqui, “Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details”, Methods and Findings in Experimental and Clinical Pharmacology 24 (Suppl. D), 2002.
• K. Uutela, M. Hämäläinen, E. Somersalo, “Visualization of magnetoencephalographic data using minimum current estimates”, NeuroImage 10, 1999.
• A. Gramfort, M. Kowalski, M. Hämäläinen, “Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods”, Physics in Medicine and Biology 57, 2012.
• A. Beck, M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”, SIAM Journal on Imaging Sciences 2, 2009.
• A. Gramfort et al., “MEG and EEG data analysis with MNE-Python”, Frontiers in Neuroscience 7, 2013.

fMRI and statistics (ch. 5–6, 9)

• S. Ogawa, T. M. Lee, A. R. Kay, D. W. Tank, “Brain magnetic resonance imaging with contrast dependent on blood oxygenation”, PNAS 87, 1990.
• N. K. Logothetis, J. Pauls, M. Augath, T. Trinath, A. Oeltermann, “Neurophysiological investigation of the basis of the fMRI signal”, Nature 412, 2001.
• R. B. Buxton, E. C. Wong, L. R. Frank, “Dynamics of blood flow and oxygenation changes during brain activation: the balloon model”, Magnetic Resonance in Medicine 39, 1998.
• K. J. Friston et al., “Statistical parametric maps in functional imaging: a general linear approach”, Human Brain Mapping 2, 1994.
• A. P. Holmes, K. J. Friston, “Generalisability, random effects and population inference”, NeuroImage 7 (S754), 1998.
• B. Thirion et al., “Analysis of a large fMRI cohort: statistical and methodological issues for group analyses”, NeuroImage 35, 2007.
• T. E. Nichols, A. P. Holmes, “Nonparametric permutation tests for functional neuroimaging: a primer with examples”, Human Brain Mapping 15, 2002.
• Y. Benjamini, Y. Hochberg, “Controlling the false discovery rate: a practical and powerful approach to multiple testing”, Journal of the Royal Statistical Society B 57, 1995.
• Y. Benjamini, D. Yekutieli, “The control of the false discovery rate in multiple testing under dependency”, Annals of Statistics 29, 2001.
• K. J. Worsley et al., “A unified statistical approach for determining significant signals in images of cerebral activation”, Human Brain Mapping 4, 1996.
• S. M. Smith, T. E. Nichols, “Threshold-free cluster enhancement: addressing problems of smoothing, threshold dependence and localisation in cluster inference”, NeuroImage 44, 2009.
• E. Maris, R. Oostenveld, “Nonparametric statistical testing of EEG- and MEG-data”, Journal of Neuroscience Methods 164, 2007.
• C. M. Bennett, A. A. Baird, M. B. Miller, G. L. Wolford, “Neural correlates of interspecies perspective taking in the post-mortem Atlantic salmon: an argument for proper multiple comparisons correction”, Journal of Serendipitous and Unexpected Results 1, 2010.
• A. Eklund, T. E. Nichols, H. Knutsson, “Cluster failure: why fMRI inferences for spatial extent have inflated false-positive rates”, PNAS 113, 2016.
• J. V. Haxby et al., “Distributed and overlapping representations of faces and objects in ventral temporal cortex”, Science 293, 2001.
• N. Kriegeskorte, R. Goebel, P. Bandettini, “Information-based functional brain mapping”, PNAS 103, 2006.
• S. Haufe et al., “On the interpretation of weight vectors of linear models in multivariate neuroimaging”, NeuroImage 87, 2014.
• G. Varoquaux et al., “Assessing and tuning brain decoders: cross-validation, caveats, and guidelines”, NeuroImage 145, 2017.

Connectivity (ch. 7, 9)

• O. Ledoit, M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices”, Journal of Multivariate Analysis 88, 2004.
• J. Friedman, T. Hastie, R. Tibshirani, “Sparse inverse covariance estimation with the graphical lasso”, Biostatistics 9, 2008.
• G. Varoquaux, A. Gramfort, J.-B. Poline, B. Thirion, “Brain covariance selection: better individual functional connectivity models using population prior”, NeurIPS, 2010.
• S. M. Smith et al., “Network modelling methods for FMRI”, NeuroImage 54, 2011.
• K. Dadi et al., “Benchmarking functional connectome-based predictive models for resting-state fMRI”, NeuroImage 192, 2019.
• J.-P. Lachaux, E. Rodriguez, J. Martinerie, F. J. Varela, “Measuring phase synchrony in brain signals”, Human Brain Mapping 8, 1999.
• G. Nolte et al., “Identifying true brain interaction from EEG data using the imaginary part of coherency”, Clinical Neurophysiology 115, 2004.
• J. Geweke, “Measurement of linear dependence and feedback between multiple time series”, Journal of the American Statistical Association 77, 1982.
• K. J. Friston, L. Harrison, W. Penny, “Dynamic causal modelling”, NeuroImage 19, 2003.

Brain–computer interfaces (ch. 8)

• L. A. Farwell, E. Donchin, “Talking off the top of your head: toward a mental prosthesis utilizing event-related brain potentials”, Electroencephalography and Clinical Neurophysiology 70, 1988.
• G. Pfurtscheller, F. H. Lopes da Silva, “Event-related EEG/MEG synchronization and desynchronization: basic principles”, Clinical Neurophysiology 110, 1999.
• H. Ramoser, J. Müller-Gerking, G. Pfurtscheller, “Optimal spatial filtering of single trial EEG during imagined hand movement”, IEEE Transactions on Rehabilitation Engineering 8, 2000.
• J. R. Wolpaw, N. Birbaumer, D. J. McFarland, G. Pfurtscheller, T. M. Vaughan, “Brain–computer interfaces for communication and control”, Clinical Neurophysiology 113, 2002.
• C. Brunner, R. Leeb, G. R. Müller-Putz, A. Schlögl, G. Pfurtscheller, “BCI Competition 2008 — Graz data set A”, Graz University of Technology, 2008.
• C. Vidaurre, B. Blankertz, “Towards a cure for BCI illiteracy”, Brain Topography 23, 2010.
• B. Blankertz, S. Lemm, M. Treder, S. Haufe, K.-R. Müller, “Single-trial analysis and classification of ERP components — a tutorial”, NeuroImage 56, 2011.
• C. Vidaurre, M. Kawanabe, P. von Bünau, B. Blankertz, K.-R. Müller, “Toward unsupervised adaptation of LDA for brain–computer interfaces”, IEEE Transactions on Biomedical Engineering 58, 2011.
• A. Barachant, S. Bonnet, M. Congedo, C. Jutten, “Multiclass brain–computer interface classification by Riemannian geometry”, IEEE Transactions on Biomedical Engineering 59, 2012.
• E. Combrisson, K. Jerbi, “Exceeding chance level by chance: the caveat of theoretical chance levels in brain signal classification and statistical assessment of decoding accuracy”, Journal of Neuroscience Methods 250, 2015.
• C. Jeunet, E. Jahanpour, F. Lotte, “Why standard brain–computer interface (BCI) training protocols should be changed: an experimental study”, Journal of Neural Engineering 13, 2016.
• F. Lotte et al., “A review of classification algorithms for EEG-based brain–computer interfaces: a 10 year update”, Journal of Neural Engineering 15, 2018.
• V. Jayaram, A. Barachant, “MOABB: trustworthy algorithm benchmarking for BCIs”, Journal of Neural Engineering 15, 2018.

Data sets and research directions (projects, ch. 10)

• A. Abraham et al., “Machine learning for neuroimaging with scikit-learn”, Frontiers in Neuroinformatics 8, 2014.
• P. Pinel et al., “Fast reproducible identification and large-scale databasing of individual functional cognitive networks”, BMC Neuroscience 8, 2007.
• D. Papadopoulos Orfanos et al., “The Brainomics/Localizer database”, NeuroImage 144, 2017.
• H. Richardson, G. Lisandrelli, A. Riobueno-Naylor, R. Saxe, “Development of the social brain from age three to twelve years”, Nature Communications 9, 2018.
• K. J. Gorgolewski et al., “The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments”, Scientific Data 3, 2016.
• O. Esteban et al., “fMRIPrep: a robust preprocessing pipeline for functional MRI”, Nature Methods 16, 2019.
• R. Botvinik-Nezer et al., “Variability in the analysis of a single neuroimaging dataset by many teams”, Nature 582, 2020.
• S. Marek et al., “Reproducible brain-wide association studies require thousands of individuals”, Nature 603, 2022.
• D. Sabbagh, P. Ablin, G. Varoquaux, A. Gramfort, D. A. Engemann, “Manifold-regression to predict from MEG/EEG brain signals without source modeling”, NeurIPS, 2019.
• D. A. Engemann et al., “Combining magnetoencephalography with magnetic resonance imaging enhances learning of surrogate-biomarkers”, eLife 9, 2020.
• R. T. Schirrmeister et al., “Deep learning with convolutional neural networks for EEG decoding and visualization”, Human Brain Mapping 38, 2017.
• V. J. Lawhern et al., “EEGNet: a compact convolutional neural network for EEG-based brain–computer interfaces”, Journal of Neural Engineering 15, 2018.
• A. Défossez, C. Caucheteux, J. Rapin, O. Kabeli, J.-R. King, “Decoding speech perception from non-invasive brain recordings”, Nature Machine Intelligence 5, 2023.
• P. Zanini, M. Congedo, C. Jutten, S. Said, Y. Berthoumieu, “Transfer learning: a Riemannian geometry framework with applications to brain–computer interfaces”, IEEE Transactions on Biomedical Engineering 65, 2018.
• M.-C. Corsi et al., “Functional disconnection of associative cortical areas predicts performance during BCI training”, NeuroImage 209, 2020.
• F. R. Willett, D. T. Avansino, L. R. Hochberg, J. M. Henderson, K. V. Shenoy, “High-performance brain-to-text communication via handwriting”, Nature 593, 2021.
• D. A. Moses et al., “Neuroprosthesis for decoding speech in a paralyzed person with anarthria”, New England Journal of Medicine 385, 2021.
• H. Lorach et al., “Walking naturally after spinal cord injury using a brain–spine interface”, Nature 618, 2023.