msg_pert — MSG Perturbation Analysis Toolbox

Systematic perturbation analysis for Magnetospinography (MSG) forward models: generating anatomically shifted geometry files and perturbed tissue conductivities, then quantifying the sensitivity of BEM, FEM, Biot-Savart, and single-sphere forward models to perturbations in source space, sensor array position, and tissue conductivity. Includes a self-contained simulations/ package that adds realistic sensor noise to evoked responses.

Developed by Maike Schmidt at the Department of Imaging Neuroscience, University College London.

---

Overview

This toolbox answers questions about MSG forward models:

Which forward model is most sensitive to imperfect anatomical knowledge? How large a shift in anatomy still produces consistent source reconstruction? And once realistic sensor noise is added, how much of the field survives?

Starting from an original geometry file produced by msg_coreg, msg_pert generates three sets of perturbed configurations and provides the tools to analyse the resulting forward fields once they have been computed in msg_fwd.

Two modalities are supported side by side: MSG (triaxial magnetometer array) and ESG (tangential/radial surface electrodes). Both are configured in one place and the analysis pipeline loops over them, so nothing has to be edited between modalities. A final combined stage compares the two directly.

Perturbation types

Type	What is perturbed	What stays fixed	Configurations
Source space	Spinal cord mesh, bone mesh, source positions	Torso, heart, lungs	24: 3 bundles × 8 random shifts
Sensor array	Sensor coil / electrode positions	All anatomy meshes	24: 3 bundles × 8 random shifts
Conductivity (BEM only)	Tissue compartment conductivities	Geometry (original)	24: 3 bundles × 8 random scalings

All three use the same bundle structure: three error bundles (small / medium / large) with 8 random realisations each.

Source space shifts model uncertainty in anatomical understanding — the whole cord–bone unit moves together. Bundles are ~2 mm (small), ~5 mm (medium), ~10 mm (large) per axis.

Sensor array shifts model registration error — e.g. the scanner cast placed a few mm off its true position. Same ~2/5/10 mm bundle scales. For ESG (surface electrodes) the shift is applied in X and Y only (Z held at 0).

Conductivity perturbations scale each BEM tissue compartment's conductivity upward by a random factor (bundles: up to +5%, +10%, +50%). Because HBF bakes conductivity into the BEM transfer matrices, each perturbation rebuilds the head model (see msg_fwd/run_conductivity_perturbation.m).

---

Two-phase workflow

Phase 1 — Perturbation generation

  msg_coreg geometry file
         │
         ▼
  pt_generate_source_shifts        →  24 source-shift geometry .mat files
  pt_generate_sensor_shifts        →  24 sensor-shift geometry .mat files
                                       (+ shift vectors to paste into config_pert)
         │
         │  (run forward models in msg_fwd)
         ▼
  msg_fwd: run_bem_leadfields / run_fem_leadfields / Biot-Savart / sphere
           on each shifted geometry (front + back arrays)
  msg_fwd: run_conductivity_perturbation  →  24 BEM conductivity leadfields
         │
         ▼
  leadfield .mat files (per geometry, per method, per array)


Phase 2 — Analysis  (run_perturbation_analysis)

         ├─ pt_load_leadfields      load BEM/FEM/BS/sphere + BEM-cond leadfields
         │                          for all perturbations → leadfields_organised.mat
  leadfields_organised.mat
         │
         ├─ pt_compute_rsq          r² per source/orientation/method (source, sensor, cond)
         ├─ pt_plot_curves          r² vs cord distance (detail, summary, cross-model)
         ├─ pt_plot_heatmaps        pairwise RE and r² heatmaps
         ├─ pt_plot_displacement    displacement / % change vs r²  (cervical + full cord)
         ├─ pt_plot_slope_vs_position   slope of r² change vs cord position
         └─ pt_compute_table        summary tables (.txt and .csv)


Optional — Realistic-measurement simulation  (simulations/, see below)

         run_simulation_analysis    evoked response + sensor noise across
                                    systems (SQUID/OP-MSG/ESG) and geometries

---

Companion Repositories

This toolbox is part of the MSG toolbox family and depends on:

msg_coreg — MSG Coregistration Toolbox (produces the input geometry files)
https://github.com/maikeschmidt/msg_coreg

msg_fwd — MSG Forward Modelling Toolbox (runs BEM/FEM/Biot-Savart/sphere on shifted geometries)
https://github.com/maikeschmidt/msg_fwd

Both must be cloned as sibling directories to msg_pert.

---

Directory Structure

msg_pert/
├── pert_path.m                    — path function (locates repository root)
├── pt_add_functions.m             — dependency setup (msg_coreg, msg_fwd, HBF)
├── config_pert.m                  — all paths and perturbation parameters
│
├── pt_generate_source_shifts.m    — PHASE 1: generate 24 source-shift geometries
├── pt_generate_sensor_shifts.m    — PHASE 1: generate 24 sensor-shift geometries
│                                    (conductivity leadfields come from
│                                     msg_fwd/run_conductivity_perturbation.m)
│
├── run_perturbation_analysis.m    — PHASE 2: master script, loops over modalities
├── pt_modality.m                  — holds the active modality (MSG / ESG)
├── pt_run_one_modality.m          — runs the 7 analysis steps for one modality
├── pt_run_compare.m               — runs the combined MSG-vs-ESG comparison
├── pt_load_leadfields.m           — load + organise all leadfields (incl. BEM-cond)
├── pt_compute_rsq.m               — per-source r² (source, sensor, conductivity)
├── pt_plot_curves.m               — r² vs cord distance figures
├── pt_plot_heatmaps.m             — pairwise RE and r² heatmaps
├── pt_plot_displacement.m         — displacement / % change vs r²
├── pt_plot_slope_vs_position.m    — slope of r² change vs cord position
├── pt_compute_table.m             — summary tables (.txt and .csv)
├── pt_compare_perturbations.m     — cross-perturbation and MSG-vs-ESG statistics
├── pt_diagnose_leadfields.m       — inspect a loaded leadfield set when results look wrong
│
├── simulations/                   — self-contained realistic-measurement package
│   ├── run_simulation_analysis.m  — master script (6 steps)
│   ├── config_sim.m               — models, systems, noise floors, geometry list
│   ├── sim_run_geometries.m       — evoked response + noise, looped over geometries
│   ├── sim_plot_original.m        — base noise curves for the original geometry
│   ├── sim_plot_comparison.m      — noise curves overlaid across geometry variants
│   ├── sim_plot_worstcase.m       — systems compared under the largest shift
│   ├── sim_plot_topoplots.m       — perfect-field topoplots per model
│   ├── sim_plot_noise_topoplot.m  — measured (noisy) topoplot at a chosen source
│   ├── functions/                 — sim helpers: sim_load_leadfield,
│   │                                sim_sensor_positions, sim_geom_file,
│   │                                sim_lf_path, sim_evoked_noise_rsq
│   └── README.md
│
└── README.md

---

Requirements

1. MATLAB (R2020a or later recommended)

2. SPM (developmental version)
   https://www.fil.ion.ucl.ac.uk/spm/

3. FieldTrip (bundled with SPM — do not install standalone)

4. Helsinki BEM Framework (HBF) by Matti Stenroos
   Clone into msg_coreg/hbf_lc_p:
   https://github.com/MattiStenroos/hbf_lc_p

5. msg_coreg sibling repository
   https://github.com/maikeschmidt/msg_coreg

6. msg_fwd sibling repository
   https://github.com/maikeschmidt/msg_fwd

---

Getting Started

Step 1: Set up dependencies

pt_add_functions;   % adds msg_pert to path; checks msg_coreg, msg_fwd, HBF

Step 2: Configure

Everything is set in config_pert.m. Paths, sensor description and method availability are declared per modality, in one block each:

base_geom_name = 'original';   % short stem in file names, WITHOUT the
                               % leading 'geometries_' prefix

mods_cfg.msg.geoms_path           = '';   % original geometry .mat
mods_cfg.msg.perturbed_geoms_path = '';   % output for shifted geometry files
mods_cfg.msg.forward_fields_base  = '';   % msg_fwd leadfield output
mods_cfg.msg.save_base_dir        = '';   % figures and tables
mods_cfg.msg.sensor_n_axes        = 3;    % 3 = triaxial MSG
mods_cfg.msg.sensor_is_meg        = true;
mods_cfg.msg.have_bem             = true; % which forward models you computed
% ... and the same block again for mods_cfg.esg (sensor_n_axes = 2)

combined_results_dir = '';                % MSG-vs-ESG comparison output
pert_modalities      = {'msg', 'esg'};    % e.g. {'msg'} to run MSG only

sensor_n_axes and sensor_is_meg are declared, not inferred: an ESG electrode count can also be divisible by 3, so guessing would mis-split the leadfield.

Step 3: Generate source-shift geometries (Phase 1a)

pt_generate_source_shifts;

This produces 24 shifted geometry files (3 bundles × 8 random shifts) in perturbed_geoms_path and prints the shift vectors to paste into config_pert.m plus a filename list for msg_fwd.

Step 4: Generate sensor-shift geometries (Phase 1b)

pt_generate_sensor_shifts;

This produces 24 shifted geometry files and prints:

1. The shift vectors — paste into config_pert.m under sensor_shift_vectors
2. A filename list for msg_fwd

Step 5: Run forward models in msg_fwd

Copy the printed filename lists into msg_fwd and run the forward models you need for both front and back arrays:

• run_bem_leadfields.m — BEM (required)
• run_fem_leadfields.m, Biot-Savart, single sphere — optional extra methods
• run_conductivity_perturbation.m — BEM conductivity leadfields (for the conductivity mode)

The per-geometry leadfield .mat files are the input to Phase 2. (msg_pert organises them itself via pt_load_leadfields in the next step — you do not run msg_fwd's own load_and_organise_leadfields here.)

Step 6: Run perturbation analysis (Phase 2)

Set the have_<method> flags for each modality in config_pert.m to match the forward models you actually computed in msg_fwd:

mods_cfg.msg.have_bem      = true;    % BEM via Helsinki BEM Framework
mods_cfg.msg.have_fem      = false;   % FEM via DUNEuro
mods_cfg.msg.have_bslaw    = false;   % Biot-Savart (infinite space)
mods_cfg.msg.have_sphere   = false;   % Single sphere (Sarvas analytical)
mods_cfg.msg.have_bem_cond = false;   % BEM with perturbed conductivities

Then run the full pipeline. It loops over every modality in pert_modalities, then runs the combined MSG-vs-ESG comparison:

run_perturbation_analysis;

To run one modality's steps standalone, choose it first, then call the steps directly:

pt_modality('set', 'msg');   % or 'esg'

pt_load_leadfields;          % load and organise leadfields (run first)
pt_compute_rsq;              % r² for source, sensor, conductivity
pt_plot_curves;              % r² vs cord distance
pt_plot_heatmaps;           % pairwise RE / r² heatmaps
pt_plot_displacement;       % displacement / % change vs r²
pt_plot_slope_vs_position;  % slope of r² vs cord position
pt_compute_table;           % summary tables

---

Script Reference

Script	Phase	Description
pert_path	—	Returns the absolute path to the msg_pert root directory
pt_add_functions	—	Adds msg_pert/functions, msg_coreg, msg_fwd, HBF, and FieldTrip wrappers to the MATLAB path
config_pert	—	Shared configuration: paths, source/sensor/conductivity parameters, naming, plot styling
pt_generate_source_shifts	1	Generate 24 geometry files for 3 bundles × 8 random source-space shifts (~2/5/10 mm)
pt_generate_sensor_shifts	1	Generate 24 geometry files for 3 bundles × 8 random sensor-array shifts
run_perturbation_analysis	2	Master script: loops over pert_modalities, then runs the combined comparison
pt_modality	2	Get/set the active modality; survives the clearvars in each sub-script
pt_run_one_modality	2	Runs the 7 analysis steps for whichever modality is active
pt_run_compare	2	Runs the combined MSG-vs-ESG comparison in an isolated workspace
pt_load_leadfields	2	Load and organise BEM/FEM/BS/sphere + BEM-conductivity leadfields; saves leadfields_organised.mat
pt_compute_rsq	2	Compute per-source r² for source, sensor, and conductivity perturbations vs the original
pt_plot_curves	2	r² vs cord distance figures (detail, summary, cross-model) for all three modes
pt_plot_heatmaps	2	Pairwise RE and r² heatmaps (within- and cross-method) for all three modes
pt_plot_displacement	2	Displacement (mm) or % conductivity change vs r² (individual: cervical; combined + trend table: full cord)
pt_plot_slope_vs_position	2	Slope of r² change vs cord position, from the displacement trend tables
pt_compute_table	2	Write median r², min r², and first-drop thresholds as .txt and .csv
pt_compare_perturbations	2	Source-vs-sensor and MSG-vs-ESG comparisons: permutation tests (aggregate and per cord position, cluster-corrected) plus Wilcoxon for reference
pt_diagnose_leadfields	2	Inspect a loaded leadfield set — sizes, scales, missing entries — when results look wrong

---

Realistic-measurement simulation (simulations/)

The perturbation pipeline above asks how the noise-free forward field changes when the model is wrong. The self-contained simulations/ package asks the complementary question: given the forward field, how much does realistic sensor noise degrade a measured evoked response, and how does that differ between sensor systems and between geometry variants?

It simulates a Gaussian-windowed evoked burst at every source on the cord, projects it through a chosen leadfield, adds trial-averaged white sensor noise across a sweep of levels, and scores r² against the noise-free field. Three systems are compared — SQUID MSG, OP-MSG, and ESG — each with its own published noise floor and measurement bandwidth. Because MSG (fT) and ESG (µV) noise floors are not comparable in absolute terms, everything is expressed as a multiple of each system's own baseline.

Run simulations/run_simulation_analysis (configure simulations/config_sim.m first). Outputs: perfect-field topoplots, per-geometry noise curves, cross-variant comparison curves (mean across cord + IQR band), a worst-case system comparison (largest shift, scored against the original field), and measured noisy topoplots. See simulations/README.md for full detail.

---

Metrics

Metrics are computed by msg_fwd's shared lf_metrics implementation, so msg_pert numbers and msg_fwd numbers are directly comparable. The unperturbed geometry is always the reference.

r² (squared Pearson correlation)

Computed per source position by comparing the full leadfield vector of the shifted model against the unshifted original:

r² = (Pearson r)^2 between shifted and original leadfield at each source
r² = 1.0 — identical leadfields (no effect of perturbation)
r² = 0.0 — no correlation

Computed separately per dipole orientation (VD / RC / LR) and sensor axis. Edges (first and last source) are excluded.

Threshold conventions:

• r² < 0.99 — first position where the perturbation has a measurable effect
• r² < 0.95 — first position where the effect is practically significant

RE (relative error) is also reported, in the pairwise heatmaps and summary tables, under the same reference-normalised definition msg_fwd uses: ‖L_shifted − L_original‖₂ / ‖L_original‖₂ × 100. See msg_fwd/INTERPRETATION.md for how to read RE against r².

---

Coordinate Convention

All geometries follow the msg_coreg scanner-cast frame (mm):

Axis	Direction
X	Left → Right
Y	Posterior → Anterior (Rostral → Caudal along cord)
Z	Inferior → Superior (Ventral → Dorsal for cord cross-section)

Source shifts are applied in this frame: a +Y shift moves the cord anterior relative to the sensor array; a +Z shift moves the cord superior.

---

Citation

If you use this toolbox, please cite it along with the companion toolboxes you used:

msg_coreg: https://github.com/maikeschmidt/msg_coreg
msg_fwd: https://github.com/maikeschmidt/msg_fwd

---

Contact

For questions, issues, or contributions, open an issue or pull request on GitHub.
Contact: maike.schmidt.23@ucl.ac.uk