lesiontrack

Longitudinal multiple sclerosis lesion evolution from two or more MRI timepoints: which lesions are new, enlarging, shrinking, stable or resolved, and which are slowly expanding lesions (SELs) by the Jacobian definition of Elliott et al. 2019, with a synthetic backtest gate that has to pass before a real number is reported.

Nothing in this repository is evidence about multiple sclerosis. It is research software. It is not a medical device, it has not been validated against histology or clinical outcome, and it makes no diagnostic claim.

Why this exists

Every MS clinic that pays for follow-up MRI analysis gets lesion-level change and, in trials, SEL counts. In the open there was one SEL implementation, bound to a GUI, needing FreeSurfer and ANTs, taking hours per subject and carrying its authors' own warning not to use it. There is no public SEL ground truth, and the field's thresholds were set by visual assessment. This package is the missing open engine, and the backtest is the missing check.

• one dependency for registration: greedy, installed as a wheel, x86-64 and arm64, CPU only
• minutes per subject pair, not hours
• every threshold in one file with its citation (lesiontrack/config.py)
• every output a plain NIfTI or TSV, on the baseline grid, ready for any viewer

Method

1. Halfway space. Rigid registration of the follow-up T1w to the baseline T1w, then the square root of that transform moves both scans into an unbiased midpoint space (Nakamura et al.; Elliott et al. 2019).
2. Deformable registration driven by T1w and FLAIR jointly (normalised cross correlation, step 0.7, update smoothing 2 voxels), baseline fixed, follow-up moving.
3. Jacobian determinant of the deformation field, converted to local percent volume change per year: (det J - 1) * 100 / years.
4. SEL candidates (Elliott 2019): voxels inside baseline lesions expanding at

   = 12.5 %/year seed 18-connected components; each seed grows by hysteresis into neighbours expanding at >= 4 %/year, competing so distinct expansions stay distinct (marker watershed on the expansion map); candidates under 10 voxels are dropped.

5. Scores (Elliott 2019 as described by Vanden Bulcke et al. 2025): concentricity is the robust slope of mean expansion per 1-voxel shell against depth from the candidate edge; constancy is the normalised residual of a robust linear fit of expansion over the intermediate timepoints (needs three or more scans). Both are z-scored across the cohort; S = z(concentricity) - z(constancy); S >= 0 is a definite SEL. Because the z-scores are cohort-relative, the definite/possible split is a ranking within your cohort, not an absolute threshold; raw values are kept. Candidates with no interior voxel have no concentricity and are marked unscorable.
6. Lesion tracking from the binary masks: connected components at baseline and follow-up are matched on the bipartite overlap graph, so merges and splits are explicit. Volume change is annualised and classified as enlarging or shrinking beyond +/-9 %/year, stable within +/-4 %/year, and trend_up / trend_down between (bands of Vanden Bulcke et al. 2025). A follow-up component with no baseline overlap is new if it clears the 1-voxel-dilated baseline mask and adjacent_fragment if it touches it; a baseline component with no follow-up overlap is resolved.

The backtest gate

lesiontrack backtest takes a real baseline (T1w, FLAIR, lesion mask), picks lesions, expands each by a known volume factor (1.15, 1.25, 1.40) with a smooth radial deformation, perturbs the whole head rigidly, rescales intensities and adds noise, then runs the full pipeline on the pair and reports:

• sensitivity: expanded lesions that receive a SEL candidate
• false positive rate: untouched lesions that receive one
• recovery fraction: robust slope of measured against true expansion (1.0 = exact)
• the noise floor inside untouched lesions and, per candidate, how well its mean rate separates injected from spurious candidates

It is the only place the pipeline is checked against a known answer. Run it on your own data before trusting a count. The numbers for MSLesSeg are in the section below.

Install

pip install lesiontrack            # once released; until then:
pip install git+https://github.com/CedricConday/lesiontrack

Python 3.10+. Wheels for greedy exist for Linux x86-64 and arm64 and macOS.

Use

Write a tab-separated manifest, one row per scan, paths relative to the manifest:

subject	session	time_years	t1	flair	mask	brainmask
P1	T1	28.09	P1/T1/T1.nii.gz	P1/T1/FLAIR.nii.gz	P1/T1/MASK.nii.gz	P1/T1/brain.nii.gz
P1	T2	28.51	P1/T2/T1.nii.gz	P1/T2/FLAIR.nii.gz	P1/T2/MASK.nii.gz
P1	T3	29.49	P1/T3/T1.nii.gz	P1/T3/FLAIR.nii.gz	P1/T3/MASK.nii.gz

time_years is any consistent clock (age, or years since baseline). Masks come from whatever segmenter you trust (LST-AI, SAMSEG, nnU-Net, manual); the brain mask is optional and only restricts the registration metric.

lesiontrack run manifest.tsv --out derivatives          # every subject
lesiontrack run manifest.tsv --out derivatives --subject P1
lesiontrack cohort --out derivatives                    # z-scores and S across subjects
lesiontrack backtest --t1 T1.nii.gz --flair FLAIR.nii.gz --mask MASK.nii.gz --out backtest

Per subject you get, on the baseline grid in halfway space:

file	content
expansion_<tp>_pct_per_year.nii.gz	Jacobian-derived local volume change
sel_candidates.nii.gz, sel_candidates.tsv	candidate labels; size, peak and mean expansion, concentricity, constancy
baseline_lesion_labels.nii.gz, lesion_labels_<tp>.nii.gz	lesion components
lesion_tracking.tsv	one row per matched group with class, volumes, rate, merge/split flags
new_or_enlarging_<tp>.nii.gz	follow-up lesion voxels outside the dilated baseline
overview_<tp>.png	four axial slices: baseline, warped follow-up, expansion with SEL outlines
summary.json	counts, volumes, intervals, every parameter used

lesiontrack cohort adds sel_cohort.tsv with z_concentricity, z_constancy, score_S and definite_sel.

Results on public data

Backtest on MSLesSeg patient P1 (2026-09-26)

Three synthetic series (seeds 0, 1, 2), each with scans at 0, 0.5 and 1 year, eight of eighteen lesions expanded by nominal factors 1.15 / 1.25 / 1.40 (analytic mean inside the lesion 11 to 25 % because elongated lesions extend past the radial window), rigid perturbation about 1 degree and 1 voxel per scan, intensity rescale 0.9 to 1.1, 2 % noise. Default registration (greedy stationary-velocity, NCC 3x3x3, halfway space).

quantity	value
recovery fraction (measured / true, robust slope)	0.29
recovery by lesion size: >= 500 / 100-500 / < 100 voxels	0.57 / 0.27 / 0.29
candidate sensitivity, Elliott thresholds	0.88
candidate false positive rate, Elliott thresholds	0.47 (1.00 on untouched lesions >= 500 voxels)
definite SEL (S >= 0, z-scored within the run, so about half by construction): sensitivity / FPR	0.42 / 0.30
noise floor: 95th percentile of per-voxel expansion in untouched lesions	5 to 10 %/yr
candidate mean expansion as a classifier, AUC	0.91 to 0.96

Three conclusions the tool would not have earned without the gate:

1. The Jacobian under-reports by about 3x on sub-voxel boundary shifts at 1 mm, and more so for small lesions. SEL rates from any pipeline of this kind are lower bounds.
2. Elliott's candidate rule is noise-limited on confluent lesions. A one-voxel peak over 12.5 %/yr, grown by hysteresis to 4 %/yr, makes a 10-voxel candidate; large lesions contain enough voxels to produce one by chance. Concentricity and constancy, the published filters, were at chance here (AUC 0.54 and 0.40).
3. The candidate's mean rate is the filter that works (AUC 0.90). It is exposed as min_mean_pct_per_year, off by default. Trade-off on this registration (mean over the three seeds):

min mean (%/yr)	sensitivity	false positive rate
0 (Elliott)	0.88	0.47
8	0.67	0.10
10	0.38	0.00

Full record with the registration sweep: docs/BACKTEST_NOTES.md.

MSLesSeg cohort (24 patients, 2026-09-26)

All longitudinal MSLesSeg patients, baseline against their last scan (median interval 1.5 years, range 0.16 to 6.0), default settings. One patient (P49) was refused by the input gate: its follow-up T1 in the released dataset has 0.1 % nonzero voxels.

Lesion tracking, all subjects, last follow-up (baseline lesions 573, matched groups): new 216, resolved 250, enlarging 127, shrinking 115, stable 36, trend 45. That churn is real signal, not just segmentation noise: voxels the masks call new go from FLAIR z 0.73 at baseline to 1.49 at follow-up and only 6 % were already lesion-bright; resolved voxels go from 1.45 to 0.87. But masks are segmented per timepoint in MSLesSeg, so "new" and "resolved" also carry every difference between two segmentations.

SEL candidates. 530 candidates in 24 subjects, median 15 per subject (12.5 with mean

= 8 %/yr); the literature reports about 2 to 4 SELs per patient over two years. On this registration the Elliott rule is noise-dominated on real data as it was on the backtest.

The Jacobian does not see what the masks see. Per subject, mask-based total lesion volume change and the Jacobian integrated over baseline lesions correlate at Spearman 0.10 (median magnitudes 21 % and 6 %). Per 1:1 matched lesion (n = 276) the correlation is 0.09; 0.18 for lesions over 100 voxels. The largest mask growers (+100 to +680 %) have Jacobian means near zero. A candidate is no more likely on a mask-enlarging lesion (42 %) than on a shrinking one (45 %).

Read together with the backtest, the picture is consistent: a 1 mm deformable registration recovers a fraction of true tissue displacement, and lesion growth that arrives as a change of intensity at the boundary rather than as displacement is largely invisible to it. Whether a stronger registration (ANTs SyN as in the original studies, or multi-resolution around lesions) changes this is exactly the question the backtest lets anyone answer for their own pipeline: inject, run, report recovery fraction and noise floor. We would like to see those numbers for the pipelines behind published SEL counts.

Cohort tables: derivatives/mslesseg/cohort_summary.tsv, sel_cohort.tsv, check_volume_change_by_lesion.tsv, check_flair_consistency.tsv (scripts/summarize_cohort.py, scripts/cohort_checks.py).

Data used

• MSLesSeg (Guarnera et al. 2025, Sci Data 12:920, CC BY): 25 patients with 2 to 4 timepoints, already in MNI152 1 mm space with a lesion mask per timepoint. scripts/make_manifest_mslesseg.py builds the manifest; intervals are age differences.
• Lesjak et al. 2016 longitudinal database (CC BY): 20 patients, 2 timepoints, expert masks of lesion change. Not yet used; it lacks baseline lesion masks.

Limitations, stated plainly

• The Jacobian sees any local volume change, including registration error at lesion edges, atrophy, and interpolation. The backtest measures the pipeline's own noise floor on untouched lesions; report it with your SEL counts.
• Two-timepoint studies cannot score constancy; S then rests on concentricity alone and is labelled so.
• Thresholds are Elliott's, set by eye in 2019 on 3D scans resampled to 1x1x3 mm. Vanden Bulcke et al. 2025 found only weak lesion-level agreement (Dice about 0.4) between SEL masks from different segmenters or resolutions. Subject-level counts were stable; lesion-level identity was not. Expect the same here.
• Masks are taken as given. A segmenter that hallucinates at follow-up will produce "new" lesions; the tool does not second-guess it.

References

• Elliott C, et al. Slowly expanding/evolving lesions as a magnetic resonance imaging marker of chronic active multiple sclerosis lesions. Mult Scler 2019;25:1915-1925.
• Vanden Bulcke C, et al. Do slowly expanding lesions correspond to chronic active multiple sclerosis lesions? An integrated imaging analysis study. medRxiv 2025, doi:10.64898/2025.12.03.25341192.
• Guarnera F, et al. MSLesSeg: baseline and benchmarking of a new Multiple Sclerosis Lesion Segmentation dataset. Sci Data 2025;12:920.
• Yushkevich PA, et al. Fast automatic segmentation of hippocampal subfields and medial temporal lobe subregions in 3 Tesla and 7 Tesla T2-weighted MRI (greedy). 2016.

License

MIT. Written by Cedric Conday with Claude (Anthropic) as coding partner.