The physics in this tool

These are the models and relations the analysis assumes, with a link to the lines of code that implement each, so that they can be checked directly. Build: 7988056203, built 2026-10-20.

The signal models and geometry

All of these are in one short file, nmrdash/physics.py, which uses only NumPy; every other part of the tool calls these functions, so what you read there is what runs.

Exponential T₂ decay

S(t) = A · exp(−t / T₂)

Transverse relaxation of one population, measured as the decay of the echo amplitudes in a CPMG train.

The code: physics.py: exp_decay, physics.py: decay_factor

Multi-exponential decay

S(t) = Σᵢ Aᵢ · exp(−t / T₂ᵢ) + c

Several populations relaxing independently (slow exchange), plus a baseline c. The multi-exponential fit.

The code: physics.py: multi_exp

Saturation recovery (T₁)

S(t) = Σᵢ Aᵢ · (1 − exp(−t / T₁ᵢ)) + c

Recovery of the longitudinal magnetisation after saturation, for one or two T₁ populations.

The code: physics.py: saturation_recovery, physics.py: multi_saturation, physics.py: recovery_factor

Diffusion attenuation

exp(−b · D)

Attenuation of the signal by free diffusion with coefficient D at diffusion weighting b.

The code: physics.py: diffusion_attenuation

Diffusion weighting b of the stimulated echo

b = (γG)² · τ² · (Δ + 2τ/3)

How strongly each diffusion delay τ weights diffusion, for a stimulated echo with storage time Δ in a static gradient G. Calculated from the gradient setting, Δ and τ in acqu.par (the echo time does not enter).

The code: physics.py: b_value_stimulated_echo, physics.py: geometric_grid, physics.py: gradient_implied_by_kea

T₂–D signal

S(t, b) = Σₖ Aₖ · exp(−t / T₂ₖ) · exp(−b · Dₖ), Aₖ ≥ 0

Several populations, each with its own T₂ and D (the few-population T₂–D model; no offset).

The code: physics.py: t2d_signal, physics.py: t2d_factors

Inversion kernels

K₁D[i, j] = exp(−tᵢ / T₂ⱼ); K_D[i, j] = exp(−bᵢ Dⱼ); K₂D = K₁D ⊗ K_D

The matrices the regularised inversions invert to obtain a T₂ distribution (1D) or a T₂–D map (2D, separable).

The code: physics.py: t2_kernel, physics.py: diffusion_kernel

Range of D the b-values can resolve

1 / b_max ≤ D ≤ 5 / b_min

A rule of thumb for which diffusion coefficients the measured b-values can tell apart; outside it a regularised map is extrapolation.

The code: physics.py: resolvable_diffusion_window

Position along the gradient

z = f / (γ·G); Δz = 1 / (N · dwell · γ·G)

Frequency to position for the static-gradient profile (FT of each echo), and the spatial resolution set by the acquisition window.

The code: physics.py: z_um_from_frequency, physics.py: native_resolution_um

What is not modelled

Restricted diffusion, or a distribution of D within a population (each population has one free-diffusion coefficient).
Non-exponential relaxation within a population, and exchange of magnetisation between populations (the populations are independent).
Finite pulse lengths and imperfect refocusing pulses in the diffusion encoding and in the echo train.
Diffusion during the CPMG echo train in the static gradient. It adds γ²G²D·tₑ²/12 to the decay rate (tₑ is the echo time), so the fitted T₂ is the apparent one. It is negligible for solids and slow diffusers and can matter for free liquids at long echo times (for D = 10⁻⁹ m²/s and G = 23.67 T/m: about 6.5 s⁻¹ at tₑ = 44 µs, about 130 s⁻¹ at 200 µs). The tool does not correct for it; only users who know their D and need the true T₂ should.

How numbers are estimated from the models

Noise, goodness of fit, model selection, intervals and the regularisation of the inversions are in nmrdash/inference.py.

Phasing and integrating the echoes

θ = −arg(Σ all points); decay_k = Σ_points Re(e^{iθ} · echo_k)

One zero-order phase puts the signal in the real channel; each echo is then summed to one number, giving the decay.

The code: inference.py: phase_and_integrate

Noise level from the imaginary channel

σ = 1.4826 · median(|x − median(x)|)

After phasing the imaginary channel holds only noise, so its robust spread measures the noise on the real decay without any model.

The code: inference.py: noise_sigma_mad

Fit quality

RSS = Σ (yᵢ − mᵢ)²; χ²_red = RSS / (σ² · dof)

The residual sum of squares and its comparison with the measured noise (about 1 when the model explains everything the noise allows).

The code: inference.py: sum_squared_residuals, inference.py: reduced_chi2, inference.py: chi2_per_point

Choosing the number of components

AIC = N ln(RSS/N) + 2k; BIC = N ln(RSS/N) + k ln N

Penalised measures for comparing fits with different numbers of components; only differences between models mean anything.

The code: inference.py: information_criteria

Confidence intervals from least squares

estimate ± t(level, dof) · SE

Intervals from the least-squares covariance with a Student-t multiplier (derived quantities such as fractions use the delta method).

The code: inference.py: student_t_critical

Regularisation of the inversions

min ‖K x − y‖² + w² ‖D₂ x‖², x ≥ 0; w = 1/√(pasdata² · pas³ · α)

The smoothness-penalised, non-negative least-squares problem the 1D and 2D inversions solve; larger α means less smoothing.

The code: inference.py: second_derivative_matrix, inference.py: tikhonov_weighting