Abstract

Background:

Spinal cord atrophy is a key biomarker for tracking disease progression in neurological disorders, including multiple sclerosis, amyotrophic lateral sclerosis, and spinal cord injury. Recent MRI advancements have improved atrophy detection, particularly in the cervical region, facilitating longitudinal studies. However, validating atrophy quantification algorithms remains challenging due to limited ground truth data.

Objective:

This study introduces SynSpine, a workflow for generating synthetic spinal cord MRI data (i.e., digital phantoms) with controlled levels of artificial atrophy. These phantoms support the development and preliminary validation of spinal cord imaging pipelines designed to measure degeneration over time.

Methods:

The workflow consists of two phases: (1) generating synthetic MR images by isolating, extracting and scaling the spinal cord, simulating atrophy on the PAM50 template; (2) performing non-rigid registration to align the synthetic images with the subject’s native space, ensuring accurate anatomical correspondence. A proof-of-concept application utilizing the Active Surface and Reg methods implemented in Jim demonstrated its effectiveness in detecting atrophy across various levels of simulated atrophy and noise.

Results:

SynSpine successfully generates synthetic spinal cord images with varying atrophy levels. Non-rigid registration did not significantly affect atrophy measurements. Atrophy estimation errors, estimated using Active Surface and Reg methods, varied with both simulated atrophy magnitude and noise level, exhibiting region-dependent differences. Increased noise led to higher measurement errors.

Conclusion:

This work presents a novel and modular framework for simulating spinal cord atrophy data using digital phantoms, offering a controlled setting for testing spinal cord analysis pipelines. As the simulated atrophy may over-simplify in vivo conditions, future research will focus on enhancing the realism of the synthetic dataset by simulating additional pathologies, thus improving its application for evaluating spinal cord atrophy in clinical and research contexts.

1 Introduction

Spinal cord (SC) atrophy is an important biomarker for assessing disease progression in various neurological conditions, including multiple sclerosis (MS) (Rocca et al., 2005; Schlaeger et al., 2014; Kearney et al., 2015; Bischof et al., 2022; Tsagkas et al., 2022), amyotrophic lateral sclerosis (ALS) (Cohen-Adad et al., 2013; El Mendili et al., 2019; Wendebourg et al., 2024), and SC injury (SCI) (Lundell et al., 2011; Trolle et al., 2023; Seif et al., 2018, 2019, 2022).

In MS, SC atrophy is particularly prominent in progressive forms of the disease and is considered among the strongest predictor of clinical disability (Tsagkas et al., 2018; Mina et al., 2021). It progresses significantly over time at different rates, with higher rates associated with clinical worsening (Rocca et al., 2019). In addition, SC atrophy seems to be unrelated to cerebral atrophy suggesting partially independent patterns of neurodegeneration in these two compartments (Cohen et al., 2012; Ruggieri et al., 2015). In ALS, a notable reduction in the cross-sectional area (CSA) of the cervical SC has been reported in patients compared to healthy individuals (Barry et al., 2022). The extent of cervical atrophy has been found to progress over time, to be predictive of shorter life expectancy in ALS patients and found to correlate with functional deficits (de Albuquerque et al., 2017; Grolez et al., 2018). In SCI, there is a rapid and progressive degeneration of the SC, which is further aggravated by a neuro-immune response (Azzarito et al., 2020; Van Broeckhoven et al., 2021).

Recent advancements in MRI technology have greatly facilitated acquisition of detailed images of the SC, especially in the cervical region. Accurate measurement of annual atrophy rates as low as a few percent is now possible, even in short-term studies with moderate sample size, through semi-automated analysis on images obtained at C1, C2, or C2–C5 and beyond (Valsasina et al., 2015; Rocca et al., 2019).

In the last decade, several research groups have developed different methods for measuring atrophy, utilizing segmentation techniques based on the contrast between the SC tissue and the surrounding cerebrospinal fluid (CSF) in MR images. Among the most popular are a semi-automatic method based on an active surface (AS) model (Horsfield et al., 2010) and two fully-automated methods based on either propagation segmentation (PropSeg) (De Leener et al., 2014) or on convolutional neural networks (Deepseg) (Gros et al., 2019), respectively. Both PropSeg and Deepseg are implemented in the Spinal Cord Toolbox (SCT) (De Leener et al., 2017). The approach of extrapolating longitudinal changes from cross-sectional measurements involves numerically subtracting CSAs obtained separately at two different time points. However, the relatively high measurement noise and low reproducibility associated with specific segmentation-based methods when measuring the small structure of the SC can affect the accuracy and reliability of the results (Prados and Barkhof, 2018). To overcome the obstacles and enhance the assessment of SC atrophy, registration-based techniques have been developed, such as the Generalized Boundary Shift Integral (GBSI) (Prados et al., 2020), Reg (Valsasina et al., 2022) and recently SIENA-SC (Luchetti et al., 2024). Unlike segmentation-based techniques, GBSI involves capturing intensity changes in the cord profile over time. Reg consists in a refinement of the AS method including accurate, slice-wise registration between time points, while SIENA-SC is an adapted version of the original SIENA method (Smith et al., 2002) designed to directly calculate the percentage of SC volume change over time. These methods have shown to result in more reliable and consistent measurements of longitudinal changes compared to “purely” cross-sectional ones (Prados et al., 2020).

In the context of MS, there is compelling evidence suggesting the presence of subclinical disease progression characterized by SC atrophy (Mina et al., 2021), which occurs before the emergence of clinical worsening (Bischof et al., 2022). It is theorized that the loss of SC volume initiates during the early phase of the disease and precedes the clinical signs of progression (Zeydan et al., 2018; Rocca et al., 2019). For this reason, there is a recognized need to develop reliable imaging analysis techniques that can be automated for use in clinical trials. Such advancements would be valuable in assessing the effectiveness of treatments aimed at slowing down the neurodegenerative mechanisms of progression (Ontaneda et al., 2023).

The concept of simulated brain imaging methodologies is well established for verifying and validating novel methods and protocols (Stöcker et al., 2025). While earlier studies used physical phantoms (Amiri et al., 2019), digital phantoms that provide more accurate and more realistic models to simulate MRI studies are increasingly being utilized (Alfano et al., 2011). Recent studies that took advantage of digital phantoms, or simulated brain imaging methodologies, include the study by Prados et al. (2020), which involved an anisotropic axial shrinkage on the cord-straightened image, or the study by Bautin and Cohen-Adad, which used digital phantoms to rescale T1-weighted (T1-w) images isotropically using different scaling factors (Bautin and Cohen-Adad, 2021). So far, these approaches served as valuable resources for evaluating the performance of various algorithms. However, one notable limitation of these assessments is the use of global scaling, which resizes all anatomical structures proportionally. In a realistic scenario, SC volume decreases only in its soft tissues, but not the surrounding bones and muscles. In the case of segmentation-based methods, where each follow-up scan is treated as an independent readout, this scaling approach might be considered acceptable. However, for registration-based methods, which rely heavily on the registration between scans, the global scaling effect is typically canceled out in the processing pipeline.

The primary objectives of this study are twofold. Firstly, we aim to establish a comprehensive workflow for generating synthetic longitudinal SC MRI data overcoming the above-mentioned limitations, thus incorporating an artificial rate of cord atrophy but at the same time being suitable for testing both on segmentation-based and registration-based methods. Secondly, we seek to conduct a proof-of-concept evaluation of two well-established methods for the quantification of changes in the upper cervical cord area over time. By accomplishing these objectives, we aim to enhance our understanding of SC analysis techniques and contribute to the development of improved diagnostic tools. In accordance with these objectives, the present work focuses on demonstrating the feasibility and reproducibility of the proposed workflow rather than capturing population-level variability in SC anatomy. The simulated atrophy provides a controlled framework for testing and validating analysis pipelines.

2 Materials and methods

In this section, we outline SynSpine, a workflow for generating synthetic MR images of the spinal cord (i.e., digital phantoms) to simulate different levels of atrophy. Furthermore, we demonstrate the practical application of these digital phantoms by highlighting its potential in evaluating the performance of tools designed for analyzing SC atrophy. The generation of synthetic images was done using MATLAB R2023b (version: 23.2) and Spinal Cord Toolbox (SCT) v6.0 (De Leener et al., 2017). The SynSpine workflow is implemented within the MRTool software suite (version 1.5.0), which is freely available at https://www.nitrc.org/projects/mrtool/. To promote reproducibility and facilitate further research, all synthetic images used in this study are publicly available at https://www.nitrc.org/projects/synspine/. Figure 1 shows an overview of the processing pipeline.

Workflow for the generation of synthetic MR images of the spinal cord. The workflow for generating and registering synthetic spinal cord images consists of two main phases: generation of synthetic MR images in template space (top) and the registration of the images from template space to subject native space (bottom). Phase 1: Firstly, the spinal cord was excised from the image (spinal cord stripping), and the excised region was filled using MATLAB’s regionfill function (canal filling), which interpolates pixel values based on the surrounding CSF. Secondly, the excised spinal cord was in-plane scaled (axial scaling) to simulate different atrophy levels and resliced using trilinear interpolation to minimize artifacts. Finally, the scaled spinal cord was superimposed onto the canal-filled images, creating synthetic templates with varying levels of simulated atrophy (e.g., PAM50 with n% atrophy). Phase 2: The non-rigid registration (W) was computed to ensure an accurate mapping of vertebral levels between the subject and the template, thus allowing a precise correspondence between the two images. Eventually, the inverted deformation (W-1) was applied to all synthetic images, transforming them into the subject’s native space and aligning their anatomical features with the subject’s spinal cord.

2.1 Description of the workflow2.1.1 Generation of synthetic MR images in template space

The initial phase involves a sequence of image manipulation steps applied to the PAM50 SC template, which establishes a standardized reference image that is commonly utilized to facilitate population-based analyses. This multimodal MRI template of the SC and the brainstem is anatomically compatible with the ICBM152 brain template and uses the same coordinate system (De Leener et al., 2018). We used only T1-w images, as these are the most commonly used in clinical and research settings to evaluate tissue structures and determine important cord geometry measures, such as the CSA (Cohen-Adad et al., 2021b).

The procedure began with the dilation of the cord mask, a component of the PAM50 template, on each axial level. This was performed using the imdilate function in MATLAB, which takes advantage of a circular structuring element with a diameter of 4 pixels to achieve the dilation effect. This enlarged mask allowed for the cord to be excised from the image without compromising the cord boundaries. Once the cord was excised (SC stripping), the corresponding area was filled by employing MATLAB’s regionfill function. This function selects pixels of CSF that encircle the excised area and smoothly interpolates them inward, taking cues from the values at the region’s outer boundary. The function executes this by computing the discrete Laplacian over the specified area and solving the Dirichlet problem to seamlessly fill in the space (canal filling). Following the region-filling step, the extracted SC image underwent in-plane scaling to the desired sizes, and then it was resliced using trilinear interpolation to reduce artifacts at the interface between white matter and CSF. In the final step, the rescaled cord images were superimposed onto the image with the filled canal, creating multiple templates with simulated atrophy levels (PAM50 with n% atrophy). For a comprehensive understanding of the described processing steps, please refer to Figure 1 for a detailed depiction.

2.1.2 Registration of synthetic images from template space to subject space

The second part of the workflow involves a series of image segmentation and registration steps utilizing SCT (De Leener et al., 2017) to align the PAM50 template (from which synthetic images with artificial atrophy rates are derived) to a real T1-w image. An initial segmentation delineates the SC in the subject’s image for precise template registration. SCT’s DeepSeg algorithm performs cord detection, centerline computation, image cropping, and final segmentation using a neural network applied to the cropped image (Gros et al., 2019). Next, vertebral labeling is performed (Ullmann et al., 2014) aligning vertebral levels between subject and template images. Two reference points, the first and last vertebral levels, are used to guide this alignment, ensuring accurate anatomical correspondence. Finally, multi-step non-rigid deformation adjusts the subject’s cord shape to match the template (De Leener et al., 2018). The first step handles large deformations; the second refines alignment. This transformation maps synthetic images from template space into the subject’s native space.

2.1.3 Addition of noise

MR images can be degraded by thermal noise which negatively affects image processing task, such as registration and segmentation (Gudbjartsson and Patz, 1995). To investigate its impact on the workflow above, increasing levels of synthetic noise were added to the original data to simulate sequentially lower signal-to-noise ratio (SNR) conditions. While magnitude images in MRI are assumed to have Gaussian noise, the non-linear mapping from real and imaginary images creates a Rician noise distribution. It is also assumed that Gaussian noise in these images has equal standard deviation. For the purpose of this simulation, images have been produced with escalating degrees of noise (percentage relative to the median intensity of the SC) according to the following formula:

With NI and I as the image with added noise and the noise-free image, respectively. n is a tridimensional array of random numbers drawn from the normal distribution with mean equal to zero and standard deviation: σ_l = n_l⋅m_sc. n_l refers to the specific noise level and m_sc refers to the median SC intensity (extracted using the PAM50 cord mask). Figure 2 shows the synthetic PAM50 template image with 1, 2, 4, and 8% added Rician noise.

Simulated images with different levels of Rician noise on an axial slice from PAM50 template. From top to bottom: 1, 2, 4, and 8% added noise.

2.2 Method validation and statistics

In order to evaluate this novel workflow, a total of six subjects were randomly chosen from the open-access quantitative MRI dataset spine generic (Cohen-Adad et al., 2021b). This dataset consists of participants from various centers who underwent scanning following a predefined protocol (Cohen-Adad et al., 2021a). The data was collected, organized, and analyzed using a well-documented procedure, which can be found at https://spine-generic.rtfd.io. Out of all the sequences available in the dataset, the defaced 3D sagittal T1-w image was utilized for the purpose of validation. Table 1 provides detailed information regarding the demographics and scanner parameters for the selected subjects, and Figure 3 shows and example of the workflow output obtained for the six randomly selected.

Subject IDSexAgeHeight (cm)Weight (kg)InstitutionScannerReceiver coilsub-fslAchieva04F2616863Santa Lucia Foundation IRCCS, Rome, ItalyPhilips AchievaNAsub-mgh04M2818379Athinoula A. Martinos Center for Biomedical Imaging, Charlestown, MA, USASiemens Skyra64ch + spinesub-mni01M3717476McConnell Brain Imaging Centre, Montreal Neurological Institute, CanadaSiemens Prisma-fit64ch + spinesub-perform05F2816754Concordia University, Perform Center, Montreal, CanadaGE MR750NAsub-sherbrooke02F3216751Center d’imagerie moléculaire de Sherbrooke, Sherbrooke, CanadaPhilips Ingenia16ch head/neck + 12ch posteriorsub-stanford02F3916354Richard M. Lucas Center, Stanford University School of Medicine, Stanford, CA, USAGE MR75016ch neurovascular

Example of the workflow output obtained for the six randomly selected subjects from the spine generic dataset. In each panel, the subject real MRI image is shown on the left and the respective synthetic images with different degree of simulated atrophy (in percentage) on the right.

To evaluate the overall workflow, we conducted a visual assessment of the impact of artificial scaling on intensity distributions within the SC mask provided with the PAM50 template. This was achieved by plotting the pairwise distributions between the intensity of the atrophy-free image and each image exhibiting incremental levels of atrophy. To evaluate the accuracy of reinserting rescaled cord images onto the canal-filled PAM50 template, as well as the effectiveness of trilinear interpolation in minimizing artifacts at the CSF/SC interface, we performed a quantitative analysis of the partial volume effect at the cord boundary. For each simulated atrophy level, a boundary ring mask was generated by computing the difference between dilated and eroded versions of the corresponding stripped SC mask. Dilation and erosion were performed using MATLAB’s imdilate and imerode functions, respectively, with a circular structuring element of 1-pixel diameter. This process produced a 2-pixel–wide ring in the axial plane, capturing the CSF/SC transition zone specific to each rescaled cord image. The ring mask for each simulated atrophy level was then applied to extract intensity values from the corresponding canal-filled PAM50 template. These values were used to quantify changes at the boundary region across simulated atrophy levels. Independent-sample t-tests were conducted to compare the original (non-atrophied) PAM50 template with each rescaled version, providing a statistical assessment of potential differences in boundary intensity resulting from the reinsertion process. To quantify the magnitude of these differences, Cohen’s d was calculated as an estimate of effect size.

Subsequently, to quantify the influence of non-rigid registration on the simulated atrophy as defined by the PAM50 template, we computed the voxel-wise differences between the intensity distributions extracted from the PAM50 SC mask on the atrophy-free image and the corresponding images with varying degrees of atrophy. This analysis was performed for both the PAM50 template and the simulated images of each subject. A two-sample Kolmogorov-Smirnov (KS) test was performed to compare the cumulative differential distributions derived from the PAM50 template with those from the subject-specific images featuring incremental levels of simulate atrophy. The KS is a nonparametric test that compares the cumulative distributions of two independent samples (PAM50 vs. synthetic subject). The test statistics D represents the maximum difference between the two cumulative distributions and it ranges between 0 and 1, where 0 occurs if the two distributions are identical, and 1 if the two are completely distinct. To assess the impact of the anatomy, the analysis was performed for each subject at three cervical levels: (i) cord segment between upper C1 extremity and C5/C6 intervertebral disk (“C1C5”); (ii) cord segment between upper C1 extremity till C2/C3 intervertebral disk (“C1C2”); and (iii) cord segment between C2/C3 and C5/C6 intervertebral disks (“C2C5”). In order to have accurate definitions of these vertebral landmarks, spinal level anchoring was defined at the mid-disc level between vertebrae using SCT labeling tools.

2.3 Proof of concept

The potential of using synthetic images to assess the effectiveness of methods for measuring SC atrophy was investigated. Specifically, the Active Surface (AS) and Reg methods within the Jim software (Xinapse Systems, Colchester, United Kingdom)1 were chosen for this assessment. The AS algorithm applies a semiautomatic active surface model, which is based on SC surface parametrization, yielding reproducible measurements of cord cross-sectional areas (Horsfield et al., 2010). In contrast, the Reg method extends the AS approach by incorporating accurate slice-wise registration between time points, enabling improved assessment of SC atrophy in longitudinal studies (Valsasina et al., 2022). The performance of both methods was then tested for increasing levels of cord atrophy (0.5, 1, 1.5, 2, 3, 4, 5, 6 and 10%) and noise (i.e., 1, 2, 4, 8%, and no noise). CSA was calculated on synthetic images in two different cervical segments (C1C2 and C2C5, as defined above). Cord atrophy is calculated by comparing the numerical differences in CSA between two time points: one with a specified simulated atrophy rate and the other serving as a reference with no atrophy applied. Notably, atrophy and noise were simulated across the entire SC in a single procedure, rather than dividing the process into two separate parts for C1C2 and C2C5. The performance of both methods was evaluated using both simulated images based on the PAM50 template and synthetic images derived from six subjects extracted from the spine generic dataset. For the PAM50 analysis, the error was determined as the difference between the measured and simulated atrophy values for each combination of noise and atrophy. In the case of the spine generic dataset, the error was assessed using the Root Mean Squared Error (RMSE) across the six subjects. This approach quantifies the overall deviation, with lower RMSE values indicating higher accuracy. In addition, the association between AS and Reg was evaluated using Pearson’s correlation coefficient (r) for each subject, separately for the cervical segments C1C2 and C2C5, and across all noise levels.

3 Results3.1 Assessing effectiveness of the workflow to generate synthetic SC atrophy images

To evaluate whether SynSpine was able to generate images with simulated levels of atrophy, the distribution of intensities of the original SC binary mask (without atrophy, 0%) was compared to the corresponding intensities of the SC with different levels of induced artificial atrophy (0.5, 1, 1.5, 2, 3, 4, 5, 6 and 10%) in the PAM50 template (Figure 4a), and each of the 6 subjects (results for a representative subject [sub-stanford02] are shown in Figure 4b; for all remaining subjects, please refer to Supplementary Figure 1). In both cases, the varying degrees of induced atrophy can be observed, as demonstrated by the deviation from the diagonal line. This deviation is proportional to the degree of induced atrophy.

Characterization of the intensity distribution for different levels of simulated atrophy in comparison to 0% atrophy for PAM50 (a) and a representative subject from the spine generic dataset (b). The intensity distributions refer to the cervical segment from C1 to C5.

3.2 Validation of SC rescaling and reinsertion: preservation of boundary integrity and intensity profiles

To evaluate the impact of cord rescaling and reinsertion on the CSF/SC boundary, intensity values were extracted within the 2-pixel–wide boundary ring mask applied to the canal-filled PAM50 templates. Visual inspection of axial slices confirmed that the mask effectively captured the transition zone, providing a sensitive region for partial volume analysis. Quantitative analysis of these intensity values revealed no significant deviations between the original (non-atrophied) template and the rescaled versions across all simulated atrophy levels. Independent-sample t-tests showed no significant differences in boundary intensity between the original template and any rescaled version (p > 0.05), indicating that reinsertion did not introduce detectable artifacts or distortions at the CSF/SC interface (Figure 5). Trilinear interpolation during reinsertion maintained smooth intensity transitions, minimizing partial volume effects and preserving anatomical realism at the cord boundaries. Assessment of Cohen’s d further confirmed the absence of meaningful differences in intensity profiles. These results validate that the simulation and reinsertion process preserves both the structural integrity and intensity characteristics of the SC in the PAM50 template across varying levels of simulated atrophy.

Validation of cord rescaling and reinsertion in the PAM50 template: preservation of CSF/SC boundary integrity and intensity profiles. Intensity values were extracted within a 2-pixel–wide boundary ring mask applied to the canal-filled PAM50 templates. This mask effectively captured the CSF/SC transition zone, providing a sensitive region for partial volume assessment. Boxplots show the distribution of boundary intensity values across increasing levels of simulated atrophy (0–10%). Independent-sample t-tests revealed no significant differences (p > 0.05) in boundary intensity between the original (non-atrophied) and rescaled templates at any atrophy level. Corresponding Cohen’s d values indicated negligible effect sizes (d ≤ 0.05), confirming the absence of meaningful deviations in intensity profiles. These findings demonstrate that the rescaling and reinsertion process preserves both the anatomical and intensity continuity of the CSF/SC interface, validating the robustness of the simulation approach.

3.3 Impact of non-rigid registration

Additionally, to evaluate the impact of the non-rigid registration used to transform the artificial image from PAM50 template space to the original subject space (as described in Figure 1), the intensity differences between the SC without atrophy and the cord with a specific atrophy rate were computed for both the PAM50 template and the non-linearly registered PAM50 template in subject space. Table 2 shows the two-sample KS test statistic (D) for different simulated atrophy values computed for each subject at three cervical levels (C1C5, C1C2, and C2C5). The D value reported in the table represents the maximum difference between the two cumulative distributions. No significant differences (p> 0.05) were observed at any cervical level for any subject (p-values are reported in Supplementary Table 1).

CervicalSubjectTwo-sample Kolmogorov-Smirnov test statistics (D) by percent simulated atrophylevelID0.5%1%1.5%2%3%4%5%6%10%C1C5sub-fslAchieva040.0130.0140.0120.0120.0140.0140.0140.0130.016sub-mgh040.0140.0150.0190.0210.0230.0210.0220.0220.021sub-mni010.020.0190.0210.020.0210.0220.0210.0210.021sub-perform050.0140.0210.0220.0240.0240.0230.0250.0240.024sub-sherbrooke020.0070.0130.0120.0130.0130.0130.0130.0120.012sub-stanford020.0070.0190.020.0210.0230.0210.0220.0230.022C1C2sub-fslAchieva040.0330.0420.0380.0420.0440.0470.0460.0480.048sub-mgh040.0140.0210.0180.0210.0210.0230.0210.0220.021sub-mni010.0290.040.0350.0390.040.0380.0410.040.038sub-perform050.020.0270.0250.0240.0250.0240.0230.0250.021sub-sherbrooke020.0140.0220.0220.0250.0220.0250.0250.0240.025sub-stanford020.0250.0240.030.0280.0290.0270.0280.0290.027C2C5sub-fslAchieva040.0120.0190.0220.0230.0190.0210.0210.0210.021sub-mgh040.0110.0170.0180.0210.0190.020.02