Brain Imaging

Derivation 1: The Hemodynamic Response Function

The BOLD signal $y(t)$ is the convolution of neural activity $n(t)$ with the hemodynamic response function (HRF) $h(t)$:

$$y(t) = \int_0^\infty h(\tau) \cdot n(t - \tau) \, d\tau + \epsilon(t)$$

The canonical HRF is modeled as the difference of two gamma functions:

$$h(t) = \frac{t^{a_1 - 1} e^{-t/b_1}}{b_1^{a_1} \Gamma(a_1)} - c \cdot \frac{t^{a_2 - 1} e^{-t/b_2}}{b_2^{a_2} \Gamma(a_2)}$$

with parameters $a_1 = 6, b_1 = 1, a_2 = 16, b_2 = 1, c = 1/6$ (SPM defaults). The HRF peaks at ~5 s, has an undershoot at ~15 s, and returns to baseline by ~30 s. In the General Linear Model (GLM), the BOLD data is:

$$\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}$$

where $\mathbf{X}$ is the design matrix of convolved stimulus regressors. Statistical testing uses t-contrasts: $t = \mathbf{c}^T\hat{\boldsymbol{\beta}} / \sqrt{\mathbf{c}^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{c} \cdot \hat{\sigma}^2}$, with correction for multiple comparisons across voxels.

Derivation 2: The EEG Forward and Inverse Problem

The forward problem relates source currents $\mathbf{J}(\mathbf{r})$ to scalp potentials$\mathbf{V}$ via the lead field matrix $\mathbf{L}$:

$$\mathbf{V} = \mathbf{L} \mathbf{J} + \boldsymbol{\eta}$$

For a current dipole at position $\mathbf{r}_0$ with moment $\mathbf{q}$ in a spherical head model, the potential at electrode $k$ at position $\mathbf{r}_k$ is:

$$V_k = \frac{1}{4\pi\sigma} \frac{\mathbf{q} \cdot (\mathbf{r}_k - \mathbf{r}_0)}{|\mathbf{r}_k - \mathbf{r}_0|^3}$$

The inverse problem ($\mathbf{J}$ from $\mathbf{V}$) is ill-posed: infinitely many source configurations produce the same scalp distribution. Minimum norm estimation (MNE) regularizes with an L2 penalty:

$$\hat{\mathbf{J}} = \mathbf{L}^T (\mathbf{L}\mathbf{L}^T + \lambda \mathbf{I})^{-1} \mathbf{V}$$

The regularization parameter $\lambda$ balances data fit against source norm. Other approaches include beamforming (LCMV), sLORETA (standardized low-resolution), and sparse Bayesian methods.

Derivation 3: Spectral Analysis and Event-Related Potentials

EEG power spectrum reveals oscillatory brain rhythms. The power spectral density is:

$$S(f) = \lim_{T \to \infty} \frac{1}{T} \left|\int_0^T V(t) e^{-2\pi i f t} dt\right|^2$$

Characteristic bands include delta (0.5–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–100 Hz). Time-frequency analysis using wavelets reveals event-related spectral perturbations (ERSP):

$$\text{ERSP}(f, t) = 10 \log_{10}\left(\frac{|W(f, t)|^2}{\overline{|W(f, t_{\text{base}})|^2}}\right) \text{ dB}$$

where $W(f, t)$ is the wavelet transform and the baseline is the pre-stimulus period. Event-related potentials (ERPs) are obtained by trial-averaging, enhancing phase-locked components while canceling non-phase-locked activity.

Derivation 4: Kinetic Modeling and Binding Potential

The simplified reference tissue model (SRTM) estimates receptor binding without arterial blood sampling. For target region concentration $C_T(t)$ and reference region $C_R(t)$:

$$C_T(t) = R_1 C_R(t) + \left(k_2 - \frac{R_1 k_2}{1 + BP_{\text{ND}}}\right) C_R(t) \otimes e^{-k_2 t / (1 + BP_{\text{ND}})}$$

where $R_1 = K_1 / K_1'$ is the relative delivery ratio,$BP_{\text{ND}} = k_3 / k_4 = f_{\text{ND}} B_{\text{avail}} / K_D$ is the binding potential:

$$BP_{\text{ND}} = \frac{B_{\max}}{K_D} \cdot f_{\text{ND}}$$

where $B_{\max}$ is receptor density, $K_D$ is dissociation constant, and$f_{\text{ND}}$ is the free fraction. Changes in $BP_{\text{ND}}$ reflect changes in receptor availability, enabling studies of neurotransmitter release, receptor occupancy by drugs, and neurodegeneration.

Derivation 5: Diffusion MRI and Tractography

Diffusion-weighted imaging (DWI) measures the Brownian motion of water molecules, which is restricted by axonal membranes. The signal attenuation for a gradient direction$\hat{\mathbf{g}}$ is:

$$S(\hat{\mathbf{g}}) = S_0 \exp\left(-b \, \hat{\mathbf{g}}^T \mathbf{D} \hat{\mathbf{g}}\right)$$

where $\mathbf{D}$ is the diffusion tensor (3x3 symmetric positive-definite matrix) and $b$ is the diffusion weighting factor. The fractional anisotropy (FA) measures directional preference:

$$FA = \sqrt{\frac{3}{2}} \frac{\sqrt{(\lambda_1 - \bar{\lambda})^2 + (\lambda_2 - \bar{\lambda})^2 + (\lambda_3 - \bar{\lambda})^2}}{\sqrt{\lambda_1^2 + \lambda_2^2 + \lambda_3^2}}$$

where $\lambda_1 \geq \lambda_2 \geq \lambda_3$ are the eigenvalues and$\bar{\lambda}$ is their mean. FA ranges from 0 (isotropic, CSF) to 1 (perfectly anisotropic, major fiber tracts). Tractography algorithms follow the principal eigenvector to reconstruct white matter pathways. Crossing fibers require higher-order models (CSD, Q-ball imaging) using spherical deconvolution.