Automatic segmentation of spinal cord lesions in MS: A robust tool for axial T2-weighted MRI scans Open Access

Multiple sclerosis (MS) is a complex autoimmune disease affecting the central nervous system, including the brain and spinal cord (Jakimovski et al., 2024). With over 2.8 million people impacted worldwide, MS is a leading cause of severe neurological disability in young adults (Walton et al., 2020). The disease is characterized by inflammatory demyelinating lesions, which are visible as hyperintensities on T2-weighted (T2w) MRI scans. Consequently, MRI has become the primary tool for diagnosing (Thompson et al., 2017) and monitoring MS (M. A. Rocca et al., 2024). However, most MRI studies focus on the brain and the measurement of atrophy in its anatomical regions (Calabrese et al., 2010; Kidd et al., 1999; Lucchinetti et al., 2011), despite spinal cord lesions being integral to MS diagnosis (Kearney et al., 2015; Thompson et al., 2017) and contributing significantly to disability (Bussas et al., 2022; Lauerer et al., 2024; M. Rocca et al., 2022).

Spinal cord MRI poses unique challenges due to its small size, tubular structure, proximity to moving organs, and the heterogeneous composition of surrounding tissues (Jasperse, 2024). These factors make lesion detection and manual segmentation cumbersome, time-consuming, and prone to variability (Gros et al., 2018). Existing open-source segmentation methods (De Leener et al., 2014, 2015; Gros et al., 2018; Naga Karthik et al., 2024) trained on specific regions of the spine fail to generalize under shifts in the data distribution caused by highly anisotropic axial scans. Furthermore, methods developed for segmenting brain lesions (Gentile et al., 2023; Mendelsohn et al., 2023; Schmidt et al., 2012; Shiee et al., 2010; Wiltgen et al., 2024) do not necessarily generalize to those of the spinal cord. In this study, we adopt an approach that simultaneously segments the spinal cord and lesions, focusing on axial T2w scans. Axial scans are superior for lesion detectability and can enhance diagnostic certainty in differentiating MS from other diseases, especially when high sensitivity for spinal cord lesions is required to confirm findings from sagittal scans or to detect small lesions potentially missed on sagittal sequences (Breckwoldt et al., 2017; Galler et al., 2016; Kerbrat et al., 2020; Weier et al., 2012). Such protocols are increasingly applied in clinical practice (Kearney et al., 2015), although they require extended scanning times and, usually, acquisition in three chunks covering the cervical, upper thoracic, and thoracolumbar part of the spine; for brevity, we will refer to the latter two simply as thoracic and lumbar (part of the spine).

Automated spinal cord lesion segmentation presents several additional challenges. The appearance of lesions along with their ambiguous boundaries due to partial volume effects in anisotropic scans poses a significant challenge to the generalizability of segmentation models. Still, available training data (i.e., images with manually annotated reference standard lesion masks) are sparse, with a particular shortage for the thoraco-lumbar spinal cord. This is critical because, although most lesions develop in the cervical cord, a considerable proportion occurs in the thoracic spinal cord (Eden et al., 2019; Hua et al., 2015; Poulsen et al., 2021). Even in the most caudal part of the spinal cord (i.e., in the medullary conus), lesions can occur and be decisive for the individual patient regarding disability and differential diagnosis (Dubey et al., 2019; Mariano et al., 2021). In addition, it is unclear which data preprocessing¹ and training strategies lead to the most optimal delineation of lesions. For instance, given the similar lesion intensities within axial slices across vertebral levels, does segmentation performance increase through additional training data from other (possibly even distant) spinal levels? We hypothesize that training on whole-spine scans obtained by stitching individual chunks and training a 3D model which leverages additional context surrounding the cord could provide more insights to this question. Further, how does the curvature of the spinal cord impact the segmentation performance? Can straightening, a preprocessing strategy that virtually erects the spinal cord to a vertical column (De Leener, Mangeat, et al., 2017), simplify the lesion segmentation task? Our hypothesis here is that straightening could potentially improve the segmentation performance by removing curvature-related variance in the spinal cord. Lastly, does segmenting both the spinal cord and the lesions simultaneously result in better performing models than those trained to segment only the lesions? Echoing the conclusion of a previous study for segmenting brain tumors (Isensee et al., 2020), we hypothesize that treating the spinal cord and lesion masks as overlapping regions could result in better gradients during training and improve segmentation performance overall.

Lastly, given previous studies have primarily focused on the automatic segmentation of MS lesions in cervical or cervico-thoracic spine (Eden et al., 2019; Gros et al., 2018), we systematically study the impact of the above data preprocessing and training strategies to develop a robust deep learning-based method for the automatic segmentation of MS lesions in axial T2w scans covering the entire spine. Specifically, we:

1. Use a region-based approach to simultaneously segment MS lesions and the spinal cord. We improve upon the performance of the existing tools (Bédard et al., 2025; Gros et al., 2018), particularly in the lower-thoracic and lumbar spines.

2. Evaluate preprocessing methods such as stitching individual chunks versus training on raw chunks and the effects of straightening versus no straightening.

3. Compare (patch-wise) 2D-/3D hybrid and (slice-by-slice) 2D kernels for training spinal cord lesion segmentation models with convolutional neural networks.

4. Investigate generalizability across various sites and imaging protocols for spinal cord lesion segmentation.

MRI data were retrospectively collected from four different sites, encompassing both cross-sectional and longitudinal datasets, with partial and full spinal cord coverage, from a total of 582 patients, following internal review board approval and performed in accordance with the Declaration of Helsinki.

The longitudinal dataset consisted of 317 patients acquired at a single site, namely the Klinikum Rechts der Isar, Technical University of Munich, Munich, Germany (referred to as the TUM dataset). Each patient underwent two MRI examinations (sessions), with each session containing 3 individual chunks of highly anisotropic axial T2-weighted scans covering the cervical, thoracic, and lumbar spines, respectively. The chunks are numbered 1-3 in the cranio-caudal direction for each spinal region. Participants for this site spanned a large variety of clinical phenotypes and conditions, namely, clinically isolated syndrome (CIS; $n=11$⁠), radiologically isolated syndrome (RIS; $n=1$⁠), primary progressive MS (PPMS; $n=11$⁠), secondary progressive MS (SPMS; $n=14$⁠), relapsing-remitting MS (RRMS; $n=277$⁠), and unknown (n = 3). The dataset comprises patients with spinal cord lesions and those without lesions at one or both longitudinal time points. The mean time interval between baseline and follow-up scan (in years) was $3.1±2.6$⁠, with an average age of $36.1±10.4$ at baseline. Of the 317 patients, the female-/male- ratio was 209/108.

The ground-truth (GT) spinal cord and lesion masks were manually annotated by two experienced raters (SR and RW). Annotations were guided by an iterative pre-labeling scheme, which involved an initial training phase where a segmentation network was trained on a manually labeled subset of the cohort. Intermediate segmentations were obtained from this trained network followed by manual corrections by raters SR and RW. This pre-labeling approach significantly expedited the annotation process, enabling the inclusion of a large number of patients and whole spine scans, which would have been practically infeasible with manual annotation alone.

Starting from the raw chunks of the individual spinal regions in the native space (Chunks Native), Figure 1A shows the dataset characteristics of the different preprocessed variants of the TUM dataset obtained after stitching and straightening the chunks. More details are given in Section 2.4 and Figure 2.

For the cross-sectional datasets, 153 patients were scanned at the NYU Langone Medical Center, New York, USA, 80 subjects were scanned at the Brigham and Women’s Hospital, Harvard Medical School, Boston, USA, and 32 were scanned at the Zuckerberg San Francisco General Hospital, San Francisco, California, USA. The scans from all three sites primarily covered the cervical or cervico-thoracic spines with only either of the chunks for each patient. The lesion masks were annotated manually by raters at the respective sites. The spinal cord masks were generated with sct_deepseg_sc using Spinal Cord Toolbox (De Leener, Lévy, et al., 2017) followed by manual corrections wherever necessary. Clinical data were not available for the MS patients across the sites. Following our experiments with the chunks from the TUM dataset, the scans from these three sites were primarily used to improve and test the segmentation model’s generalization capabilities. Figure 1B shows how the data from different sites are combined for developing the lesion segmentation model. Table 1 shows a detailed overview of the image characteristics and acquisition parameters from all sites in this study.

In clinical practice, full spinal cord scans are impractical due to long acquisition times. An alternative is to stitch individual segments into whole-spine scans. In this paper, we refer to stitching as the process of aligning and merging multiple image segments or chunks, in this case, axial T2-weighted (T2w) scans acquired during the same session, into a single image. It must be noted that chunks may also be registered to the template individually; however, it is preferred to directly register the stitched image as this simplifies the computation of lesion morphometrics by eliminating the need to reconcile overlapping regions or risk double-counting that can occur when multiple chunks are processed separately (Graf, Möller, et al., 2024).

To compare whether the model benefits from the extended context around the spinal cord in stitched whole-scans, we employed the stitching algorithm implemented in (Graf, Platzek, et al., 2024; Lavdas et al., 2019). This algorithm has four main steps: (i) the corner points are computed for each chunk, (ii) using the affine matrix of each chunk, corner points are rotated to find the smallest volume bounding box that fits all the chunks, resulting in optimal rotation and spacing parameters, (iii) the chunks are resampled to the new (common) image space and occupancy maps containing overlapping regions between chunks are stored, and (iv) the resampled chunks are blended using ramp-based interpolation with a weighted averaging function to ensure smooth transitions between chunks. Manual corrections are applied to the stitched mask, if required, to reduce interpolation artifacts. Figure 2A-C shows the result of the stitching algorithm and the subsequent appearance of lesions in the whole-spine scan.

While the stitching procedure described in (Graf, Platzek, et al., 2024) includes a registration-based alignment alternative, we chose not to use it in this study, instead relying on the global alignment of chunks within the scanner coordinate space. This approach could theoretically introduce misalignment due to patient movement, similar to intra-scan motion. However, our quality checks showed that misalignment was very rare. When it did occur, it was typically in non-overlapping regions where empty slices resulted from suboptimal MRI planning. Moreover, given the minimal chunk overlap and the inherent challenges of registering axial T2w spinal cord images—such as high anisotropic resolution, partial volume effects, and the lack of distinct anatomical landmarks for segmentation-based registration—we found that a registration-based approach introduced additional alignment errors and was ultimately inferior in our stitching experiments.

In the analysis of brain MS lesions, MRI scans of the brain are typically registered to the MNI anatomical template using rigid registration (Carass et al., 2017; Wiltgen et al., 2024). Following this approach of analyzing the lesions in a template space, one would have to register the individual chunks to, for example, the PAM50 template (De Leener et al., 2018). As it is difficult to obtain the vertebral levels in axial scans with thick sagittal slices, we opted for a simpler approach involving straightening of the spinal cord (De Leener, Mangeat, et al., 2017). Straightening essentially eliminates all curvature-related variance in the spinal cord, offering a simpler alternative to template-based registration and facilitating cohort-level analysis. The straightening procedure using Spinal Cord Toolbox (SCT) (De Leener, Lévy, et al., 2017) is described as follows: Given an axial T2w scan (chunk or stitched), we applied sct_straighten_spinalcord using the spinal cord segmentation mask as reference to obtain the straightened cord along with the warping field (warp_curve2straight). Then, the output warping field was applied using sct_apply_transfo with linear interpolation to straighten both the GT spinal cord and the lesion masks. Lastly, both the straightened GT masks were binarized using a threshold of 0.5. Figure 2D-F shows the straightened individual chunks, straightened whole-spine scans and the resulting appearance of lesions across spinal regions, with a reduced field of view resulting from the straightening algorithm.

The continuing dominance of nnUNet (Isensee et al., 2021, 2024) across several open-source segmentation challenges has shown that a well-tuned convolutional neural network (CNN) architecture is robust and continues to achieve state-of-the-art results over novel transformer-based architectures (Hatamizadeh et al., 2021; Karthik et al., 2024; Ma et al., 2024). Based on this rationale, we used nnUNetV2 (Isensee et al., 2021) as the segmentation framework. Specifically, we used nnUNet’s region-based training strategy, where the model treats the masks of the spinal cord and lesions as overlapping regions and segments them simultaneously. The rationale for going with this approach is two-fold: (i) combining the label of a larger object (i.e., the spinal cord) with a smaller one (i.e., the lesion) might help obtain better gradients for effective learning early on during the training, and (ii) it has been shown that directly optimizing for overlapping regions, forming a hierarchy of sub-structures (e.g., spinal cord $→$ lesions), results in better performing models than optimizing for each of the objects as independent classes (Isensee et al., 2020). In practice, this essentially reframes the two-class segmentation task (with independent spinal cord and lesion masks) into two binary segmentation tasks (where, the first class is now the union of the spinal cord and the lesion masks and the second class is the lesion mask) optimized using the sigmoid activation with the compound Dice and binary cross-entropy loss. While region-based training requires both spinal cord and lesion masks as inputs, the rise of robust spinal cord segmentation tools (Bédard et al., 2025; Gros et al., 2018; Naga Karthik et al., 2024) across various MRI contrasts and pathologies enables obtaining spinal cord segmentation masks without significant manual annotation costs.

We used the 2d and 3d_fullres models in their default configurations as automatically set by nnUNet. The 2D model employs 2D convolutional kernels throughout all stages of the U-Net architecture. However, for 3d_fullres model, the default configurations resulted in a hybrid 2D/3D kernels given the highly anisotropic nature of our dataset. In particular, out of the 6 layers in the encoder, the 1st layer used a $3×3×1$ kernel with stride $1×1×1$ (i.e., no downsampling/pooling), the 2nd and 3rd layers used $3×3×1$ kernels with a stride of $2×2×1$ (i.e. in-plane downsampling only), and the 4th, 5th, and 6th layers used $3×3×3$ kernels with a stride of $2×2×2$⁠. Note that the first 3 layers have singleton dimensions and the downsampling/pooling is only done when the resolution (or, slice thickness) of the out-of-plane dimension lies within a factor of 2 of the in-plane resolution. Therefore, it must be noted that despite being a 3D model (i.e., trained on 3D input patches), the actual configuration consists of a mix of 2D and 3D convolutional kernels.

Default data augmentation transforms applied by nnUNet were used, namely, random rotation, scaling, mirroring, Gaussian noise addition, Gaussian blurring, adjusting image brightness and contrast, low-resolution simulation, and gamma transformation. All scans were converted to RPI orientation and preprocessed using Z-score normalization. For each of the four preprocessed datasets, 2D and 3D models were trained using three-fold cross validation resulting in a total of 24 models. All models were trained for 1000 epochs, with a batch size of 2, using the stochastic gradient descent optimizer with a polynomial learning rate scheduler.

While nnUNet is known for automatically configuring the patch sizes of the images used during training based on the hardware capacity, its default patch sizes could be sub-optimal depending on the size of the input images. Specifically, in the case of the stitched datasets, the automatically-configured patch sizes for 2D/3D-hybrid models were defined to be half the size of the median shape of the images from the training set. This means that the model effectively received only half of the total length of the (stitched) spinal cord in the S-I plane as a 3D input patch during training. To evaluate whether fitting the whole context of the cord within a single input patch would improve the segmentation performance, we also trained a model on a patch size covering the entire S-I plane, meaning that the model has access to the whole context of the spinal cord in the sagittal plane.

All segmentation networks were trained with nnUNet v $2.4.1$ using Python $3.10.13$ and Pytorch $2.2.1$⁠. Training both 2D and 3D models (with hybrid 2D/3D kernels) on any given dataset variant took approximately 1 day on a single NVIDIA A100 GPU with 40 GB memory.

Experiments were designed following two major themes: (i) investigating the effects of data preprocessing (namely, stitching and straightening the chunks) on lesion segmentation using the TUM dataset, and building on this (ii) leveraging data from additional sites to obtain a robust and generalizable segmentation model for MS lesions in axial T2w scans.

Regarding the first theme, we created four versions of the TUM training dataset (see Fig. 1A) to evaluate how the different formats of the input scans affected the downstream lesion segmentation performance.

1. chunks native, the original dataset containing the individual chunks of axial T2w scans corresponding to cervical, thoracic, and lumbar regions of the spine,

2. stitched native, consisting of a single scan per patient with the individual cervical, thoracic and lumbar chunks stitched together to form the whole spinal cord,

3. chunks straightened, where the individual chunks in their native space are straightened using the spinal cord segmentation masks and the model is trained in the straightened image space, and,

4. stitched straightened, where the stitched images in their native space are straightened using the whole-spine cord segmentation masks.

We trained 2D and 3D models on each of the four dataset variants, resulting in a total of 8 models to investigate the impact of input spatial context. As mentioned in the previous section, the 3D models used a mix of 2D and 3D convolutional kernels owing to the dataset characteristics. Please refer to Table A1 in the Appendix A for a detailed comparison of training configurations. Only anisotropic, axially acquired T2-weighted scans were used in our experiments. Furthermore, because of the evolution in the appearance of lesions between sessions, the images were treated as independent inputs for training the segmentation model. The TUM dataset was split patient-wise according to 80-20% train/test ratio to ensure that the chunks for a given patient belonged either to the training set or testing set but not to both. For the chunks native dataset, this resulted in a total of 1522 individual chunks from 254 patients for training/validation and 380 chunks from 63 patients for testing. For the stitched native dataset, this translated into 508 stitched whole-spine scans from 254 patients in the training/validation set and 126 scans from 63 patients in the testing set.

As for the second theme focusing on robustness and generalization, all 153 patients’ chunks from the NYU site were added to the training set while keeping the BWH and UCSF sites unseen, in order to test the model’s generalization to out-of-distribution data. Such a training and evaluation approach partly overlaps with a lifelong learning scenario where the segmentation model is continually enriched with data from new sites by adding them to the current training set and testing on data from unseen sites (Naga Karthik et al., 2022). Figure 3 shows the distribution of lesions in terms of the lesion volume across train and test splits for each site.

We compared our models on the TUM test set with three other open-source methods available in SCT: (i) sct_propseg (De Leener et al., 2014), (ii) sc_deepseg_sc (Gros et al., 2018), and (iii) the recently proposed contrast-agnostic segmentation model (Bédard et al., 2025). The test split was done at the patient-level (and not at the chunks-level), ensuring that the axial scans of a particular patient strictly belonged only to the testing set, resulting in a total of 126 patients. The performance of these models was then evaluated independently on each of the three chunks (i.e., cervical, thoracic and lumbar) from the TUM dataset.

We created three test sets: the TUM $(n=63)$ and the chunks from sites, BWH $(n=80)$⁠, and UCSF $(n=32)$ (which were kept unseen during training) to evaluate the model’s performance on out-of-distribution data. First, to determine the dataset variant and training strategy that leads to best performance, 8 models (4 dataset variants, 2 models each after combining the results from 3 folds) were compared. Then, the best model out of the 8 (based on the Dice score), the model trained on chunks from two sites, along with sct_deepseg_lesion (Gros et al., 2018), were evaluated on the TUM, BWH, and UCSF test sets.

While the ANIMA toolbox (Commowick et al., 2021), specifically animaSegPerfAnalyzer, is widely used for the evaluation of brain MS lesion segmentation models (Commowick et al., 2021), it fails to output appropriate metrics in special cases where the ground-truth mask contains no lesions, irrespective of whether the model predicts false positive lesions or not. Skipping such cases in the evaluation of spinal cord MS lesions would result in (biased) higher scores for the metrics by not accounting for the false positive rate of the segmentation models. To address these limitations, we used the open-source MetricsReloaded package (Maier-Hein et al., 2024; Reinke et al., 2023), designed to mitigate several shortcomings of existing open-source toolboxes. In the context of this study, when the GT and the automatic prediction are both empty (i.e., contain no lesions), the voxel-wise Dice is set to 1, indicating the model’s ability to learn correctly. When either of them is non-empty, the respective false-positive and false-negative scores are computed by formal definition. As for metrics, for spinal cord segmentation, we computed the voxel-wise Dice score, normalized surface distance (NSD), and relative volume error (RVE). For lesion segmentation, we reported the voxel-wise Dice score, NSD, and average absolute difference of the lesion count between the number of predicted lesions and GT lesions, and lesion-wise $F1$ score, where a prediction was considered true positive if it showed at least a $10%$ overlap with the GT lesion.

To ensure a fair comparison between the models trained on the four variants of the TUM dataset, we used the following post-processing strategies before computing the metrics:

1. Chunks Native: For each patient, the 3 chunks corresponding to cervical, thoracic, and lumbar spines were stacked together to obtain a single prediction for a given patient (instead of having 3 lesion prediction masks per patient). As the sizes of the chunks differed depending on the spinal region, an appropriate amount of padding was applied to the chunks to enable stacking. Figure 4 illustrates our concept of stacking and algorithm 1 describes it in practice.

2. Stitched Native: No post-processing was done in this case, as we already have a single prediction per patient.

3. Chunks Straightened: To ensure that all predictions are evaluated in the same space, the straightened chunks were brought back to the image space by applying the inverse warping field. This was followed by stacking the chunks (now in the native space) to obtain a single prediction per patient.

4. Stitched Straightened: Similar to the straightening of chunks, the straightened stitched images were transformed back to the native space by applying the inverse warping field.

Given the difference in the magnitude of $F1$-scores computed on individual and stacked chunks, all metrics described in Section 3 were computed on stacked chunks. Notably, the process of stitching chunks can lead to the interpolation of lesions in overlapping areas, while the process of stacking may result in double-counting these lesions. To assess whether this has a significant impact on our dataset, we analyzed our TUM test set. Among a total of 126 stitched scans (or respectively 378 chunked scans), we observed 450 lesions in the chunked dataset and 448 in the stitched dataset, a difference that can be attributed to a single patient. Given the minimal frequency and negligible impact of this discrepancy, we decided not to account for this effect in our analysis.

Statistical analysis was performed using the SciPy Python library v $1.9.1$ (Virtanen et al., 2019). Data normality was tested using D’Agostino and Pearson’s normality test. Group-wise comparison between the 2D and 3D models for each of the four variants of the TUM dataset (thus resulting in 8 models to compare) was performed using the non-parametric Kruskal-Wallis H-test. We decided on the Kruskal-Wallis H-test because of the non-parametric nature of our data (as concluded from the normality test) and due to the fact that we were comparing more than two groups. Post hoc pairwise test for multiple comparisons was performed using Dunn’s test with holm correction.

This section compares the results of our proposed model for spinal cord segmentation against various existing methods (Section 3.1). Next, we examine how spinal cord straightening affects the performance of the lesion segmentation model (Section 3.2). We then demonstrate the model’s capability to detect lesions of varying sizes and evaluate its performance in detecting lesions across different chunks (Section 3.3). Finally, we identify the best model based on these findings and enhance it further by incorporating additional data from external sites (Section 3.4).

Figure 5 compares the performance of 2D/3D models trained on chunks/stitched images respectively with state-of-the-art methods for spinal cord segmentation evaluated independently across different spinal regions (cervical, thoracic, and lumbar). Interestingly, spinal cord segmentations from sct_propseg and sct_deepseg_sc vary drastically across chunks, with the lowest performance on chunk-3 covering thoracolumbar levels. The contrast-agnostic v2.5² model (Bédard et al., 2025) reduces the gap in segmentation accuracy across chunks with a notable improvement in the lumbar spine segmentation. However, the models trained with TUM data (chunk-wise or stitched) perform significantly better than other methods for our test set.

In this section, we compare the lesion segmentation performance across the different variants of the TUM dataset described in Section 2.4. Figure 6 shows multiple raincloud plots comparing the lesion segmentation performance of models trained on individual axial slices (2D) and three-dimensional image patches (3D). Training a 2D model resulted in a higher Dice and $F1$-score compared to the 3D models in all datasets except for the stitched native variant.

As mentioned in Section 2.3, along with default patch sizes defined by nnUNet, we also experimented with modified patch sizes covering the entire length of the cord within a single input patch to see if extended patch sizes would improve segmentation performance. Interestingly, we found no statistically significant differences between the 3D models trained with default patch sizes and the models trained with a full S-I coverage. Hence, we proceeded with the default patch sizes for all models.

Considering the role of spinal cord curvature in lesion segmentation, we observed that the models trained on straightened cords performed slightly worse than the models trained in the native space of the input scans. Specifically, in 2D, a median Dice score of 0.72 (chunks native) vs. 0.69 (chunks straightened) and, in 3D, a median Dice of 0.72 (stitched native) vs. 0.53 (stitched straightened). This is also apparent from the wider distribution of the Dice scores of predictions (i.e. scatter points) in Figure 6A. In a between-dataset pairwise comparison of Dice scores, the 2D chunks native model significantly outperformed 2D stitched native$(**P<.01)$⁠, 3D chunks straightened$(***P<.001)$⁠, and 2D/3D stitched straightened models $(***P<.001)$⁠. Likewise, for lesion-wise $F1$-scores, the 2D chunks native model showed significant differences between 3D chunks straightened $(***P<.001)$, and 3D stitched straightened models $(***P<.001)$⁠.

Considering the median Dice and $F1$-scores, the Chunks Native 2D model ranks higher than the straightened counterparts and is on par with the models trained on stitched scans. In other words, training directly on chunks is a simple strategy that does not perform worse than the stitched models without the need for any stitching procedure. As this presents a simple and scalable solution, we proceeded with the 2D model and present the results from our subsequent experiments evaluating lesion detection in this section.

Figure 7 shows the lesion detection rates across various categories of lesions sizes. The performance of the model increases with increasing lesion sizes even detecting small lesions between 10-50 $mm3$ at a 60% rate. We observed that lesion detection becomes challenging as lesions get even smaller.

To evaluate whether our model is consistent in detecting lesions occurring at various spinal cord regions, we plotted the distribution of lesion sizes and the corresponding lesion detection rate across spinal cord chunks in Figure 8. As expected, we observed a decreasing trend in the frequency of lesions going from the cervical to thoracolumbar spine. Our model achieved similar lesion detection rates in cervical and cervicothoracic spines (left and middle), with a slight decrease in the thoracolumbar spine (right).

We have now established that the model trained on chunks performs at least on par as the stitched models without any explicit preprocessing (namely, stitching and straightening) with similar lesion detection rates across chunks. Following the experimental design, we added axial T2w chunks from an external site (NYU) and evaluate whether the performance of the model improves after training on two sites (TUM and NYU). We tested the model in two ways: (i) on an in-distribution test set (i.e., TUM) to evaluate whether the model’s performance has drifted after the addition of an external site, and (ii) on unseen (out-of-distribution) data from two additional sites (BWH and UCSF), acquired using different scanners and protocols unseen during training.

Table 2 and Figure 9 compare three models qualitatively and quantitatively: (i) sct_deepseg_lesion, (ii) the model trained on TUM dataset’s chunks (i.e. single site), and (iii) the model trained on two sites. As data from BWH and UCSF are acquired at different sites and with acquisition parameters, the results on these sites show the model’s robustness to out-of-distribution data. Both variants of the models trained on chunks performed significantly better than sct_deepseg_lesion (Gros et al., 2018), which tends to miss lesions (shown with yellow arrows) resulting in higher false negatives. Importantly, adding data from the NYU site did not degrade the model’s performance on the original TUM test set; rather, it showed a slight improvement across all metrics. As for the evaluation metrics on unseen sites, the model trained on two sites performed better than the single site model, although no statistically significant difference was found. A comparison with seg_ms_lesion, a contemporary work accessible in SCT, is presented in Appendix A.2.

In this study, we developed an automatic tool for the joint segmentation of the spinal cord and MS lesions in axial T2w MRI scans using a region-based training approach. Starting with a comprehensive comparison between chunks, stitched, and straightened variants of our dataset, we showed that training a 2D model directly on the raw chunks in the native space achieved the best performance, presenting a simple and scalable solution. With respect to spinal cord segmentation, our model improved over the existing methods on the lumbar spine, whereas, for lesions, the proposed model detected up to $90%$ of lesions larger than 50 $mm3$⁠, while struggling to detect small lesions reliably. Given that MS lesions manifest along the entire length of the spinal cord, our model could detect lesions across the entire spinal axis, with slightly lower detection capabilities in the thoracolumbar regions. We summarize our key insights and perspectives on spinal cord and MS lesion segmentation in this section.

All benchmark methods (namely, sct_propseg, sct_deepseg_sc and contrast_agnostic) showed a decrease in performance towards the lower thoracic and lumbar spinal cord. This is not surprising as these models were trained on datasets that predominantly comprised cervical and cervicothoracic chunks. In contrast, given that the proposed model was trained exclusively on chunks covering the whole spine, it achieved consistent segmentations throughout various regions of the cord. Particularly, as the current state-of-the-art (Bédard et al., 2025) struggles to segment the thoracolumbar spinal cord, integrating our model into SCT will provide a novel solution for the segmentation of thoracolumbar spine axial T2w scans. Another likely reason for the better performance of our model is that the benchmark methods were not trained on the TUM data, leading to an expected drop in performance when evaluated on out-of-distribution data.

Straightening the spinal cord negatively impacted the lesion segmentation performance. At first glance, this may seem surprising, given that straightening eliminates curvature-related morphological variations across individuals. However, the straightening procedure introduces heavy deformations in the vertebrae and the discs surrounding the spinal cord, possibly blurring anatomically meaningful information; in addition, small hyperintensities occurring at the concave side of the cord shrink, while lesions on the convex side get enlarged. We hypothesize that these non-linear deformations negatively impacted the performance of the models. Furthermore, resampling involved in the straightening procedure could degrade the signal in the images caused due to interpolation, further affecting the performance. Additionally, for models trained in the native space, lesion segmentation can be performed using the raw chunks as input, whereas segmenting lesions on straightened scans requires spinal cord masks a priori at inference.

In addition to straightening the spinal cord, we could have registered the images to an anatomical template (e.g., PAM50 (De Leener et al., 2018)) and subsequently train a segmentation model in the template space, to further reduce the inter-patient variability. Similar approaches involving the MNI template have been reported in brain MS lesion segmentation studies (Carass et al., 2017; Wiltgen et al., 2024). However, registration to the PAM50 template requires robust identification of the vertebral levels, which was difficult to achieve here due to thick slices in the sagittal plane. While several deep learning-based methods exist for vertebral labeling in CT scans (thanks to the open-source, labeled datasets (Sekuboyina et al., 2021)), it remains a challenge in MRI scans. Existing studies are limited in their approaches concerning fixed fields-of-view, spinal region-specific models (i.e., cervical/lumbar) and do not provide vertebral labels for whole-spine images (Azad et al., 2021; Moller et al., 2024). On the other hand, TotalSpineSeg (Warszawer et al., 2024), a recently-introduced model, obtains robust vertebral labels on a wide variety of contrasts. Using vertebral labels from TotalSpineSeg and registering to the PAM50 template, future works could explore lesion segmentation directly in the template space, thus simplifying the tasks of lesion mapping and computing lesion morphometrics.

Recent studies highlight the significant impact of MS lesions at cervical, cervicothoracic, and thoracolumbar levels (Bussas et al., 2022; Eden et al., 2019; Hua et al., 2015; Kearney et al., 2015; Poulsen et al., 2021) and recommend extended axial coverage of the spinal cord (Breckwoldt et al., 2017; Galler et al., 2016). However, as acquiring full spinal cord scans suffers from long acquisition times, stitching independently acquired segments offers a simpler method to obtain whole-spine scans. We observed that training models on whole-spine scans did not result in statistically significant performance improvements over training on individual chunks, nor did our additional experiments with extended patch sizes covering the entire spinal cord (in the S-I plane). This suggests two things: (i) the model learns to be sensitive to hyperintensities in the cord without requiring any additional 3D patch-wise or length-wise contextual information, (ii) on a practical note, one can simplify the lesion segmentation task on axial T2w scans by training on raw chunks (i.e., without stitching). Notably, this approach is more scalable as many acquisition protocols do not cover the whole spinal cord.

The stitching procedure requires resampling both individual chunks and their ground truth masks to a common resolution. While the additional interpolation may affect boundary-sensitive metrics such as the Dice score, our post-stitching quality control ensured that small lesions or lesions in the overlapping regions between chunks are preserved, thus not affecting lesion sensitivity or false positive rates.

Spinal cord lesions, which often appear as elongated, blob-like structures, frequently span across multiple axial slices and vertebral levels (see Fig. 2). Existing studies (Gros et al., 2018; Walsh et al., 2024) have only reported results from exclusively 3D models, motivating the comparison of 2D and 3D models in this study. Given the highly anisotropic axial scans, the 3D models were trained with a combination of 2D kernels with singleton dimensions and fully 3D kernels at different stages of the encoder to balance the in-plane and out-of-plane dimensions. Our experiments revealed that 2D-only models often outperformed the hybrid 2D/3D models, or at least performed on par. This suggests that the additional spatial context provided by 3D patches did not add significant value. This is evident in Table 2, where sct_deepseg_lesion (Gros et al., 2018), a 3D model trained on patch sizes of $48×48×48$⁠, performed significantly worse than our proposed 2D models, even though it was trained on a subset of axial scans. In contrast, pure 2D models effectively mitigate the destabilizing effect of partial volume by training slice-by-slice on high-resolution axial slices and any variations in image intensities across vertebral levels may have a less significant impact on 2D kernels than on 3D architectures.

Consistent with previous studies (Bussas et al., 2022; Waldman et al., 2024), we observed a progressively decreasing trend in the occurrence of lesions from the cervical to the thoracolumbar spine within the TUM cohort. Notably, our analysis revealed that the lesion detection rate remains consistent across different spinal segments. This suggests that detection performance is influenced more by lesion size than by spinal level, highlighting the need to address limitations in segmenting smaller lesions by incorporating reliable training samples for categories that the model currently struggles to segment accurately.

In particular, our results demonstrate that the model struggles to reliably segment spinal cord lesions smaller than 10 $mm3$⁠. Existing studies (Walsh et al., 2023) also reported high inter-rater variability for small lesions among expert raters. More importantly, our approach and model improve substantially with increasing lesion size, achieving approximately 60% accuracy for lesions between 10 and 50 $mm3$⁠. Although these lesions are relatively small, they are clinically significant (Bussas et al., 2022). Recent studies have shown that using synthetically generated lesions as an additional data augmentation method improved the performance over traditional data augmentation strategies (Basaran et al., 2022; Zhang et al., 2021). Using such methods targeting the synthetic generation of small lesions could help improve our model’s performance across small lesion sizes.

In MS, distinguishing patients with and without spinal cord lesions is crucial (Lauerer et al., 2024; M. A. Rocca et al., 2024). Ideally, the training data should contain a small proportion of patients with no lesions (or even healthy controls) and the evaluation metrics should account for the empty ground-truth masks. These considerations help the model in distinguishing between MS lesions and common spinal cord artifacts (appearing as abnormal hyperintensities). The ANIMA toolbox, commonly used in lesion segmentation challenges (Commowick et al., 2016, 2021), skips empty lesion masks, preventing proper estimation of false positive rates. In contrast, MetricsReloaded (Maier-Hein et al., 2024) assigns a Dice score of 1 when both GT and predictions are empty and 0 otherwise, ensuring empty cases are accounted for. While this might skew the Dice score, it is justified given the importance to correctly predict for empty masks. Thus, we used MetricsReloaded in this work and introduced an evaluation approach that stacks individual chunks before computing metrics for a fair comparison between models trained on chunks vs. whole-spine scans.

An important aspect to consider is the use of different (pre-)labeling strategies. The TUM masks were pre-segmented with manual labels and an iterative approach employing nnUNetV2, while, in the data from NYU, UCSF, and BWH, lesions were labeled entirely manually by different raters and spinal cord masks were prelabeled with SCT (Gros et al., 2018) and subsequently corrected. This variability contributed positively to the generalizability of our proposed models.

Directly adding raw chunks from the NYU site highlighted the simplicity of training on chunks. The two-site model outperformed the single-site model on the TUM test set and two unseen external datasets, benefiting from the greater diversity in scanners and acquisition protocols. This improved robustness makes it more suitable for integration into SCT (De Leener, Lévy, et al., 2017).

Our work has several limitations: (1) we note that only axial T2w scans were considered in this work. Although our training and test datasets include axially acquired scans with a wide range of resolutions, we did not evaluate our model on high-resolution axial or isotropic scans falling outside the resolution range of our datasets, (2) the proposed model is constrained to the segmentation of lesions in T2w scans. Incorporating a variety of orientations and contrasts in the training set would improve its performance on several contrasts (Bédard et al., 2025). However, this generalization may come at the cost of reduced accuracy for specific contrasts and orientations. Furthermore, clinical protocols frequently include complementary scans in various orientations (Clara Weyer, 2021; Wattjes et al., 2021), such as axial and sagittal T2w, which offer additional information that our current models are unable to use for its predictions. Leveraging multiple contrasts or views is technically challenging, as current deep learning models typically require registered images to be aligned and stacked as multiple channels using the same field of view. Enhancing models to effectively leverage such complementary information, for example, via super-resolving scans first, is an ongoing area of research (McGinnis et al., 2023). (3) We note that nnUNet’s automatic configuration of patches, kernels, and other training hyper-parameters may still be suboptimal depending on the task. Although we experimented with minor adjustments to the initial parameters (detailed in 2.3), further optimization is still an open area for investigation. Our comparison of training plans for the same model across different dataset configurations revealed that nnUNetV2 adapts settings such as patch size and stride based on the data. While this demonstrates nnUNetV2’s ability to tailor its plans to the data, these adjustments can also affect key variables, such as lesion sensitivity, which we aim to optimize. Given the vast number of possible configurations, a comprehensive ablation study would be computationally intensive and beyond the scope of this study. We acknowledge this limitation and emphasize the need for further exploration in future work. (4) The performance of our model is limited by the inter-rater variability in the GT masks, as they were obtained independently at each site with different annotation strategies. Understanding the effect of inter-rater agreement in MS lesion segmentation is an active area of research (Walsh et al., 2023), our study does not explore this in greater detail. Lastly, the model was developed using data from a private cohort. While this prevents future reproducibility of our study, the model is made publicly-available as part of SCT, thus establishing a baseline for future studies in MS lesion segmentation in axial T2w scans.

We developed an automatic tool for the joint segmentation of spinal cord and MS lesions in axial T2w MRI scans covering the entire spine. We evaluated both patch-wise 3D and slice-by-slice 2D training strategies applied to individual chunks, stitched whole-spine scans, and straightened spinal cords to identify the combination that yields the best lesion segmentation performance. Our experiments demonstrated that slice-by-slice 2D training on raw chunks achieved the highest segmentation accuracy. This approach is simple and scalable, while not involving stitching and straightening. To highlight the practicality of our findings, we further improved the proposed model by training and evaluating on chunks from external sites as they are easily available and need not cover the entire spine. The model performed well on unseen data and slightly improved after extending the training data, suggesting good generalizability. To facilitate broader use, we open-sourced the code and integrated the model into SCT (v7.0 and above).

The data used is private and could be made available upon reasonable request. To facilitate reproducibility and open science principles, all codes, including scripts for preprocessing, training, and generating plots are open-source and can be found at https://github.com/ivadomed/model-seg-ms-axial-t2w. Additionally, the segmentation model is accessible via SCT (v7.0 and above) using sct_deepseg lesion_ms_axial_t2 -i <path-to-image.nii.gz>.

E.N.K.: data curation, formal analysis, investigation, methodology, visualization, and writing (original draft, review & editing). J.M.: data curation, formal analysis, investigation, methodology, visualization, and writing (original draft, review & editing). R.W.: data curation, segmentation. S.R.: data curation, segmentation. R.G.: methodology, writing (review & editing). J.V.: methodology, writing (review & editing). M.L.: data curation, writing (review & editing). P.L.B: data curation, methodology, and writing (review & editing). J.T., R.B., S.T, T.S., A.B., C.Z., and B.H.: data curation. D.R.: conceptualization, supervision, writing (review & editing). B.W.: conceptualization, data curation, methodology, supervision, and writing (review & editing). J.S.K.: conceptualization, data curation, funding acquisition, investigation, methodology, supervision, and writing (review & editing). J.C.A.: conceptualization, data curation, funding acquisition, investigation, methodology, supervision, and writing (review & editing). M.M.: conceptualization, data curation, funding acquisition, investigation, methodology, supervision, and writing (review & editing).

E.N.K. is supported by the Fonds de Recherche du Québec Nature and Technologie (FRQNT) Doctoral Training Scholarship and DAAD (German Academic Exchange Service) Short-term Research Grant. J.M., M.M., and J.S.K. are supported by Bavarian State Ministry for Science and Art (Collaborative Bilateral Research Program Bavaria – Québec: AI in medicine, grant F.4-V0134.K5.1/86/34). R.G. and J.S.K. are supported by European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (101045128-iBack-epic-ERC2021-COG). J.V. received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101107932. J.C.A. is supported by the Canada Research Chair in Quantitative Magnetic Resonance Imaging [CRC-2020-00179], the Canadian Institute of Health Research [PJT-190258], the Canada Foundation for Innovation [32454, 34824], the Fonds de Recherche du Québec - Santé [322736, 324636], the Natural Sciences and Engineering Research Council of Canada [RGPIN-2019-07244], the Canada First Research Excellence Fund (IVADO and TransMedTech), the Courtois NeuroMod project, the Quebec BioImaging Network [5886, 35450], INSPIRED (Spinal Research, UK; Wings for Life, Austria; Craig H. Neilsen Foundation, USA), Mila - Tech Transfer Funding Program. The authors thank Digital Research Alliance of Canada for the computed resources used in this work.

E.N.K., J.M., R.W., S.R., R.G., J.V., P.L.B., M.L., J.T., S.T., T.S., A.B., C.Z., B.H., J.S.K., D.R., J.C.A., and M.M. have no known competing financial interests to declare. R.B. has received speaking honoraria from EMD Serono and Sanofi and research support from Bristol Myers Squibb, EMD Serono, and Novartis. B.W. has received speaker honoraria from Philips and Novartis (unrelated to this study).

We thank Mathieu Guay-Paquet and Joshua Newton for their assistance with dataset management and for their contributions to implementing the sct_qc tool. We extend our gratitude to Erbil’s for their consistent quality of vegan döners, making the completion of this project feasible.

Table A1 shows the training configurations 2D and 3D models on each dataset variant of the TUM dataset. Note that in all 3D models, the convolutional kernels contain a mix of 3D kernels with singleton dimensions (effectively being applied as 2D kernels on the in-plane dimensions) and full $3×3×3$-shaped 3D kernels.