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

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

Two deformable engines

The rigid step, the halfway frame and the reslicing are greedy in both cases; --engine picks what runs on the halfway-space images. greedy (default) needs nothing extra. ants uses SyN through antspyx, cross-correlation on T1 and FLAIR jointly, and is how the original SEL studies registered:

pip install 'lesiontrack[ants]'
lesiontrack run manifest.tsv --out derivatives_ants --engine ants
lesiontrack backtest --t1 T1.nii.gz --flair FLAIR.nii.gz --mask MASK.nii.gz --out backtest_ants --engine ants

Both engines write the same files, so the backtest compares them on identical inputs. The Jacobian determinant is computed by lesiontrack from the warp, in the image's voxel frame, for either engine; see "The Jacobian is ours" below for why.

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). Numbers re-measured 2026-09-27 with the corrected Jacobian on the same registrations (the 2026-09-26 values, computed through greedy's -rj, are kept in docs/BACKTEST_NOTES.md).

quantity	value
recovery fraction (measured / true, robust slope)	1.16
recovery by lesion size: >= 500 / 100-500 / < 100 voxels (median)	1.00 / 1.22 / 1.21
candidate sensitivity, Elliott thresholds	1.00
candidate false positive rate, Elliott thresholds	0.47 (1.00 on untouched lesions >= 500 voxels, 0.07 under 100)
definite SEL (S >= 0, z-scored within the run): sensitivity / FPR	0.67 / 0.37
noise floor: median peak per-voxel expansion in untouched lesions	10.8 %/yr
candidate mean expansion as a classifier, AUC	0.99 to 1.00

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

1. The Jacobian recovers the injected expansion in full, slightly over on small lesions (the radial window's falloff is counted in the mean). The earlier "under-reports by 3x" was greedy's -rj on LAS images, see "The Jacobian is ours" below.
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; every untouched lesion over 500 voxels produced one. Concentricity and constancy, the published filters, were at chance on the 2026-09-26 run and are not re-evaluated here.
3. The candidate's mean rate is the filter that works (AUC 0.99). 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)	1.00	0.47
8	1.00	0.20
10	1.00	0.07
12.5	0.92	0.00

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

MSLesSeg cohort (24 patients, 2026-09-26, Jacobian re-measured 2026-09-27)

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. 484 candidates in 24 subjects, median 16 per subject (12.5 with mean

= 8 %/yr), 97 definite; 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, and min_mean_pct_per_year around 10 is where the backtest puts the floor.

The Jacobian and the masks agree in rank, not in size. Per subject, mask-based total lesion volume change and the Jacobian integrated over baseline lesions correlate at Spearman 0.41 (medians +12.5 % by mask, -6.3 % by Jacobian, brain-wide mean -1.5 %). Per 1:1 matched lesion (n = 276) the correlation is 0.22; 0.25 for lesions over 100 voxels. Masks are segmented per timepoint in MSLesSeg, so their change carries every difference between two segmentations; the Jacobian sees only displacement of the surrounding tissue, and its median inside persisting lesions is negative.

Whether ANTs SyN (--engine ants) changes any of 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.

The Jacobian is ours

Until 0.2.0 the Jacobian came from greedy's -rj in the reslice command. That option differentiates the composed warp while it is still in physical (LPS, millimetre) units, along voxel index axes, without the image direction matrix, so it is only right for images whose voxel axes are LPS-aligned at 1 mm. On a phantom with an identity (RAS) affine, a sphere growing from radius 8 to 10 voxels (volume ratio 1.95) came back as det 0.68; with a LAS affine, as used by the MSLesSeg scans, as 1.18. For a radial stretch 1+k the true determinant is (1+k)^3; each flipped axis turns one factor into (1-k), so RAS (two flips) reports (1-k)^2 (1+k), which is below 1 for any expansion, and LAS (one flip) reports (1-k)(1+k)^2, roughly a third of the true volume change for small k. greedy's standalone -jac converts the warp first and is correct; the reslice path does not (found 2026-09-27, draft report in docs/upstream/greedy_rj_direction.md).

lesiontrack now computes det(I + grad u) itself from the warp file, after taking the vectors LPS -> RAS -> voxel units through the inverse affine, and tests/test_jacobian.py checks it under four orientations and anisotropic spacing. Registrations cached under an older record are not redone; their Jacobian is recomputed from the stored warp.

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.