Comparative validation of AI and non-AI methods in MRI volumetry to diagnose Parkinsonian syndromes

In this section, we describe the brain MRI data (Section “Data preparation”), FS implementation (Section “Brain structure segmentation: baseline with FS”), and DL method implementation (Section “DL models for brain structure segmentation”) for the volumetric analysis of key brain structures to diagnose neurodegenerative diseases. Figure 1 shows an overview of the study process considering the evaluation and comparisons between FS and DL models (i.e., modified V-Net and UNETR representing CNN and ViT DL architectures, respectively). Supplementary Fig. S1 shows a diagram of the overall performance comparison. We developed DL models with faster processing but similar segmentation performance to FS. The DL models were trained to reproduce and segment the results of FS for each brain structure \(F_i \in [0,1]^{256 \times 256 \times 128}\) as model output \(V_i \in [0,1]^{256 \times 256 \times 128}\) by taking skull-stripped brain image \(I \in \mathbb {R}^{256 \times 256 \times 128}\) as input (\(i \in \{pallidum, \, putamen, \, caudate, \, third \, ventricle, \, midbrain, \, pons \}\)), with resolution (h, w, d) (height \(h=256\), width \(w=256\), depth \(d=128\)). The DL segmentation results for the six brain structures were stored as 3D binary masks (\(F_i\) and \(V_i\) indicate the FS and DL-model masks for brain structure i, respectively), where each mask output contained intensities between 0 and 1 (area outside and inside the target brain structure, respectively). By calculating the absolute volume of each or all the brain structures predicted by FS or DL models, we performed binary classification of PD, MSA-C, MSA-P, PSP, and normal cases, and calculated the area under the curve (AUC) of segmentation.

Ethical approval

All authors of this study confirm that all methods or experiments were performed in accordance with the Declaration of Helsinki and the relevant guidelines and regulations provided by the policies of the Nature Portfolio journals. This study was approved by the Institutional Review Board of the Samsung Medical Center (IRB number: SMC 2021-07-026). The written informed consent of the patients was waived by the Institutional Review Board of Samsung Medical Center because we used deidentified and retrospective data.

Data preparation

Study population and clinical assessments

We retrospectively screened patients from the Neurology Department of Samsung Medical Center between January 2017 and December 2020. Patients diagnosed with PD, probable MSA, or probable PSP were included in this study. The diagnosis for each patient was determined by movement disorder specialists based on the following criteria: PD was determined according to the United Kingdom PD Society Brain Bank criteria³⁹ using [18F] N-(3-fluoropropyl)-2β-carbon ethoxy-3β-(4-iodophenyl) nortropane positron emission tomography, while probable MSA and PSP were diagnosed according to the second consensus diagnosis of MSA⁴⁰ and movement disorder society clinical diagnostic criteria for PSP⁴¹, respectively. MSA cases were further classified as either MSA-P or MSA-C after reaching consensus⁴⁰. Patients with concomitant or structural brain lesions, including stroke and tumors, which may affect brain MRI scans, were excluded from the study. An age-matched healthy elderly population was included as the control group. Demographic information on age, sex, and disease duration until the brain MRI examination was collected, as listed in Table 4. We analyzed the data from 411 individuals and performed threefold cross-validation to train and evaluate the DL models. Each group consisted of 105 healthy controls and 105 PD, 69 PSP, 69 MSA-C, and 63 MSA-P cases.

We applied cross-validation with three outer folds for evaluation to mitigate bias in the validation and test sets and analyze the effect of set composition (combinations of cases in groups). The data were randomly divided into three sections, one for testing and two for training. Each group comprised 35 normal, 35 PD, 23 PSP, 23 MSA-C, and 21 MSA-P cases.

Data acquisition and standardization

Axial brain MRI scans were acquired using a standard protocol for T1-magnetization-prepared rapid acquisition of gradient echo, with repetition/echo time of 11,000/125 ms, inversion time of 2800 ms, field of view of 240 mm, acquisition matrix size of \(320 \times 249\), echo train length of 27, 1 signal average, slice thickness of 5 mm, interslice gap of 1.5 mm, and scanning time of 198 s.

We included six brain structures that are involved in Parkinsonian syndromes in the gray matter, namely, the midbrain, pons, putamen, pallidum, caudate, and third ventricle. These areas are reported to have the highest sensitivity and specificity for differentiating Parkinsonian syndromes^13,16. The MRI scans were resized to \(256 \times 256 \times 128\) (i.e., number of slices in the coronal/sagittal/axial planes) to segment each structure.

The FS accepts Digital Imaging and Communications in Medicine (DICOM) or Neuroimaging Informatics Technology Initiative (NIfTI) files as inputs. DICOM is a compelling and flexible but complex format that provides interoperability between several hardware and software tools. Given its complexity, DICOM format was converted to NIfTI format⁴². NIfTI is a more straightforward format than DICOM and preserves the essential metadata. In addition, it maintains the volume as a single file and uses raw data after a simple header, and NIfTI files can be loaded and processed faster than DICOM files for whole brain images. Therefore, we converted files in the brain MRI DICOM format into files in the NIfTI format using MRIcroGL.

Brain structure segmentation: baseline with FS

The extraction of brain structures obtained using atlas-based automated segmentation are necessary for training and validation before establishing an automated DL segmentation model. In this study, we used these results as DL ground-truth labels and evaluated the validity of DL model for generating the same label. As a representative technology for atlas-based automated segmentation (see details in Supplementary Section A.2), we selected FS (version 7.2), which is publicly available for neuroscience research and provides high segmentation performance^{18,19,20,21,43,44}. Additionally, FS no longer supports CUDA, thus unable to calculate the time using GPU.

To segment and extract the six brain structures using FS, it sequentially executes the recon-all⁴⁵ pipeline and Brainstem Substructure pipeline⁴⁶. We used both pipelines because the recon-all pipeline does not support segmentation of brainstem structures (e.g., pons and midbrain). However, because the Brainstem Substructure pipeline receives pre-processed inputs from the recon-all pipeline, both pipelines should be executed. Therefore, the extraction of the six brain structures through FS can be divided into MRI scan pre-processing in the recon-all pipeline and the remaining segmentation of the recon-all pipeline along with segmentation in the Brainstem Substructure pipeline. These processes are explained in Section “Data preparation” and Section “Brain structure segmentation: baseline with FS”.

MRI scan pre-processing for FS: motion correction and skull removal

The MRI scan pre-processing in the recon-all pipeline of FS mainly consists of (1) motion correction, (2) normalization, and (3) skull stripping. Motion correction is conducted before averaging when various source volumes are used, compensating for small motion variations between volumes. FS constructs cortical surface models and the boundary between white matter and cortical gray matter to automatically match the brain images of patients, using software¹⁷. In addition, intensity normalization is applied to the original volume. However, adjusting for intensity fluctuations may hinder intensity-based segmentation. Instead, we scale the intensities of all voxels to the mean value (110) of white matter.

After correcting for motions and normalizing the data, FS removes the skull and provides the skull-stripped brain MRI scan. Removing intracranial brain cavities (e.g., skin, fat, muscle, neck, and eyeballs) may reduce human rater variability⁴⁷ and promote automated brain image segmentation and improve analysis quality. Therefore, brain MRI scans should be pre-processed to isolate the brain from extracranial or nonbrain tissues in a process known as skull stripping⁴⁸. FS developers devised and applied in-house automated skull-stripping algorithms to isolate intracranial cavities by default.

In this study, the steps of brain MRI scan pre-processing (i.e., skull stripping with motion correction and normalization of a brain MRI scan) took approximately 20 min. We converted the final skull-stripped images to NIfTI files with size of \(256 \times 256 \times 128\), while the original brain MRI scan had a size of \(256\times 256 \times 256\), which was adjusted for efficient comparison with the DL models.

FS for brain structure segmentation

After pre-processing (Section “Brain structure segmentation: baseline with FS”), FS segments the six brain structures by applying the remaining processes of the recon-all pipeline and the complete Brainstem Substructure pipeline. After skull stripping, registration-based segmentation proceeds as follows. FS determines and refines the white and gray matter interfaces for both hemispheres. Then, FS searches for the edge of the gray matter, which represents the pial surface. With pial surfaces, FS expands and inflates sulci banks and gyri ridges. Subsequently, it extends again into a sphere and parcellates the cortex. After applying these processes, FS segments the brain. The recon-all pipeline encompasses some brain structures (i.e., putamen, caudate, pallidum, and third ventricle), while the Brainstem Substructure pipeline segments the midbrain and pons.

In this study, the final segmentation result was assessed with the same input size of \(256 \times 256 \times 128\). The original size of the segmentation result was \(256 \times 256 \times 256\), but it was adjusted to \(256 \times 256 \times 128\) for comparison with the DL models. In addition, we replaced FS with a DL model applied to the skull-stripped MRI scan (i.e., pre-processing result of the recon-all pipeline) to perform segmentation. For the replacement, we evaluated whether the DL analysis is faster than FS analysis and whether the segmentation result of DL is sufficiently reproducible compared with that of FS. The differences between FS and DL segmentation are illustrated in Fig. 2.

DL models for brain structure segmentation

In this study, we used DL models and FS to segment the same skull-stripped images (i.e., images pre-processed by the FS recon-all pipeline, as described in Section “Brain structure segmentation: baseline with FS”). The original size of the skull-stripped image generated by FS was \(256 \times 256 \times 256\), which was adjusted to \(256 \times 256 \times 128\) for DL segmentation owing to the limited GPU memory. We evaluated and compared the performance and analysis time of the DL models by replacing the segmentation process of FS after skull stripping with DL. FS may be inefficient because it segments the entire brain image, requiring many hours of processing. In fact, FS takes at least 4.5 h to segment the six brain structures considered in this study because it requires atlas-based registration to transform the coordinates of the entire MRI scan to segment specific brain structures. Consequently, FS cannot notably reduce the processing time even if only six brain structures were to be segmented. On the other hand, we verified that DL segmentation (e.g., using V-Net or UNETR) takes less than 1 min to 18 min per case to segment the six target brain structures. As DL models do not require complex registration, unlike non-artificial-intelligence methods (e.g., FS), they can substantially increase the processing efficiency. The implementation details of the DL models are described herein. As DL models, we adopted the CNN-based V-Net²⁹ and ViT-based UNETR³⁰ using the segmentation results provided by FS as labels (Section “Brain structure segmentation: baseline with FS”). The two models were trained to reproduce FS segmentation.

CNN-based V-Net

Architecture of CNN-based 3D segmentation using V-Net. ResBlock, MaxPooling, and UpConvolution were used to reduce the depth, height, and width. The output shown in the figure is the segmentation of pallidum. (Conv convolution layer, BN batch normalization).

V-Net has been used to segment an entire volume after training an end-to-end CNN on MRI volumes for revealing the prostate^29,49,50 The architecture of V-Net is V-shaped, where the left part of the network is a compression path, whereas the right part decompresses the features until the original input size is recovered. The left part of the network is separated into stages that operate at varying resolutions.

In this study, one to three convolutional layers were used in each step. A residual function was learned at each level. The input of the residual part was used in the convolutional layers and nonlinear operations. This output was added to the last convolutional layer of the stage. The rectified linear unit (ReLU) was used as the nonlinear activation function. Convolutions were applied throughout the compression path. The right part of the network learned a residual function similar to that of the left part. V-Net has shown promising segmentation results, and using this model in our application improved performance. The model was adjusted according to the available memory. The proposed architecture is illustrated in Fig. 3. The left part used a residual block (ResBlock) and maximum pooling (MaxPooling). ResBlock was applied to all the blocks with an input size of \(256 \times 256 \times 128\). On the other hand, 3D MaxPooling reduced the depth, height, and width of the feature maps to reduce their resolution. The right part also used ResBlock but replaced MaxPooling with UpConvolution, which consisted of 3D upsampling, batch normalization, ReLU activation, and convolutional layers (\(5 \times 5\times 5\) filter, same padding, and stride of 1). Upsampling increased the resolution of the feature maps, and batch normalization improved convergence throughout the network⁵¹.

ViT-based UNETR

Architecture of ViT-based UNETR directly connected to a CNN-based decoder via skip connections at different resolutions for segmentation. (Deconv deconvolution layer, Conv convolution layer, BN batch normalization, MLP multilayer perceptron).

UNETR³⁰ is a transformer architecture for 3D medical-image segmentation. There is a study that used UNETR as brain tumor segmentation⁵², but no study was held for brain parts segmentation. It uses a transformer as the encoder to learn the sequence representations of the input volume and capture global multi-scale information while adopting U-shaped architectures for the encoder and decoder. The proposed architecture is illustrated in Fig. 4. UNETR followed a contracting–expanding path with an encoder comprising a stack of transformers connected to a decoder through skip connections. The encoder directly used 3D patches and was connected to a CNN-based decoder via a skip connection. A 3D input volume was split into homogeneous nonoverlapping patches and projected onto a subspace using a linear layer. Position embedding was applied to the sequence and then used as input to the transformer. The encoded representations at different levels in the transformer were retrieved and sent to a decoder via skip connections to obtain the segmentation results.

Implementation details of DL models: training and inference

For the DL models, the input comprised a brain mask and the corresponding patient’s segmented brain structures in the MRI scans, which were merged into an array of dimension \(256 \times 256 \times 128\). The ground truth of each brain structure was segmented using FS. For evaluation, threefold cross-validation of the test data was applied to calculate the Dice score and Dice loss. We implemented V-Net in TensorFlow and Keras and trained it for 100 epochs. For UNETR, PyTorch and MONAI⁵³ were applied, and the model was trained for 20,000 iterations. Both models used Python language and were trained using an NVIDIA Tesla V100 DGXS GPU with a batch size of 1 and an initial learning rate of 0.0001. For CPU, Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20Ghz was used.

We evaluated the accuracy of the evaluated models using the Dice score by comparing the expected segmentation with V-Net (or UNETR) and FS outputs. The Dice score measures the overlap between the reference and predicted segmentation masks. A Dice score of 1 indicates perfect spatial correspondence between the two binary pictures, whereas a score of 0 indicates no correlation. We used the Dice loss to determine the performance of the three outer cross-validations on their test sets for the corresponding structures. If \(F_i\) and \(V_i\) are the ground-truth mask and its prediction for each brain structure, respectively (i.e., FS segmentation mask \(F_i\) and its DL prediction mask \(V_i\), respectively, as shown in Fig. 1), the Dice score⁵⁴ for each brain structure \(i \in\) \(\{\)pallidum, putamen, caudate, third ventricle, midbrain, pons\(\}\) is derived as

where \(\circ {}\) denotes the Hadamard product (i.e., component-wise multiplication) and \(||\cdot ||_{1}\) is the L1-norm (i.e., sum of absolute values of all components). Moreover, we measured the segmentation time for evaluation.

Statistical analysis for binary classification of cases

We obtained the absolute volumes from the six segmented brain structures (i.e., pons, putamen, pallidum, midbrain, caudate, and third ventricle) predicted by the DL models (i.e., CNN-based V-Net or ViT-based UNETR) or FS. Based on the absolute volume of the individual brain structures, we calculated the AUC of the binary classification of diseases, normal vs. P-plus, normal vs. PD, and PD vs. P-plus cases. The AUC was computed based on the receiver operating characteristic curve produced by the correlation between the predicted absolute volume of each brain structure and each case.

Disease binary classification was conducted using the six segmented brain structures individually or collectively. For individual analysis, the AUC was derived through thresholding-based binary classification by obtaining the absolute volume of the individual structures. For a comprehensive analysis of all structures, we additionally considered an ML classification algorithm to perform disease binary classification with the six volumes as inputs. For the classification algorithm, binomial logistic regression (LR) and extreme gradient boosting (XGBoost) were used. LR is a statistical model widely used in ML classification^55,56,57. XGBoost is a well-established method that produces advanced results among gradient-boosting-based techniques⁵⁸ (e.g., XGBoost successfully won 17 out of the 29 ML tasks posted on Kaggle by 2015⁵⁹). In both methods, we evaluated the AUC obtained by the DL model and FS through threefold cross-validation.