1.

Introduction

Distributed learning has become a promising alternative to centralized learning for training machine learning (ML) models in medical image analysis, addressing many patient privacy regulations and overcoming administrative barriers.¹ It offers a practical solution for accessing large and diverse datasets by enabling ML model training in distributed environments. Despite its success, distributed learning does not address the inherent challenge of image acquisition biases across multiple centers,¹ which is especially important when no imaging protocol harmonization is used.

Image acquisition biases can, for example, be caused by differences in imaging acquisition protocols and/or scanner types across multiple centers, potentially leading the ML model to learn patterns unrelated to the main task (i.e., spurious correlations).² Several previous studies have explored the effects of image acquisition biases. For instance, Refs. 3 and 4 demonstrated that combining data acquired with different pulse sequences or scanners in multi-center studies can lead to biases in brain volumes when conducting brain morphology analysis for cross-sectional or longitudinal studies. In addition, variations in cortical thickness measurements were observed due to distinct imaging acquisition protocols, magnetic field strength, and pulse sequence parameters.^5–7 Moreover, Glocker et al.⁸ found that simple intensity-based harmonization techniques cannot eliminate all of a scanner’s encoded effects. Similar findings were also reported by Ref. 9, demonstrating that ML models can accurately identify the site where brain scans were acquired, even after intensity harmonization. Recently, Souza et al.² showed that scanner types can be classified from the internal feature representations of a model trained for Parkinson’s disease (PD) classification in a centralized approach. Acquisition bias has also been reported for other data modalities when working with multi-center datasets such as molecular data¹⁰ and magnetic resonance imaging (MRI)-derived features.^11,12 All of these studies reinforce that ML models can identify and possibly exploit image acquisition biases, such as scanner type, as potential shortcuts for the disease classification task, practically diminishing their clinical utility.¹³

Shortcuts can emerge in ML models for various reasons, as described in more detail by Geirhos et al.¹³ First, models may disproportionately focus on making good predictions for the majority group, leading to potential misclassifications in minority group(s). Second, ML models often operate with minimal effort, meaning that once a model identifies a feature to perform its task, it may rely exclusively on that feature, even if it represents a data artifact. Third, a model might learn features unrelated to the intended task, such as classifying scanner types instead of disease groups. Last, a model may combine multiple features to make a decision. For example, during a disease classification task, the model might incorrectly associate scanner type with disease status, even though the scanner used to acquire the data should not influence the diagnosis.

Hence, recent ML research has aimed to develop methods to harmonize data across different centers to reduce imaging biases and improve ML models’ reliability, generalizability, usability, and validity.¹⁴ Data harmonization is the process of transforming data from different sources into a common format or reference frame to enable meaningful comparisons or integration. Most data harmonization methods in ML that employ MRI volumetric data utilize generative adversarial networks (GANs).^15–20 However, many of these methods have significant limitations, requiring paired data (e.g., data from the same individual using every scanner or image modality to be harmonized across), and large quantities of data from each site. These factors can hinder the implementation of GANs in distributed learning environments, especially in scenarios involving numerous centers, where each center contributes only a few data points for training.

Distributed learning is commonly implemented following the federated learning (FL) paradigm. In FL, all centers train a copy of a model initialized by a server in parallel. After each training round, the server aggregates the parameters from the models trained at each center, updates the global model, and sends it back to the centers for further training. This process is repeated until the aggregated model meets predefined convergence criteria.²¹ Managing centers with limited datasets presents significant challenges in FL setup, including the risk of the local models overfitting and the difficulty of defining an aggregation function that does not marginalize centers with fewer datasets.²² Although data harmonization methods have been proposed specifically for the FL,^23–25 this has been limited to the scenario where centers contribute large amounts of data for training.

In contrast to FL, the traveling model (TM) visits one center at a time in a sequential process.^22,26 Therefore, the model is initialized at a server or the first center, followed by training with the data available in the first center. Subsequently, the updated model is transferred to the next center, and training continues with locally available data. This process is repeated until the final center is reached, completing one training cycle. Multiple cycles can be conducted to enhance the model’s overall performance.²² Despite being less explored than FL, previous research has found that the TM approach to distributed learning is especially effective when centers are only able to provide a few datasets for training.^22,26 Instances in which healthcare facilities possess limited data are commonly found in scenarios involving rare diseases with low prevalence,²⁷ small or rural hospitals admitting a reduced number of patients, and hospitals in low- or middle-income countries due to a lack of sufficient imaging equipment, skilled clinicians, and radiologists.²⁸ The effectiveness of TM arises from the iterative training of a single model, which alleviates the FL challenge of the local model’s overfitting to some extent.

Recently, the TM method has been successfully applied to classify PD using three-dimensional (3D) MRI data acquired across 83 centers worldwide.²⁶ Although the TM demonstrated state-of-the-art performance on this task [achieving an area under the receiver operating characteristic curve (AUROC) of 83%, comparable with the centralized approach of 80%], it did not address or investigate the inherent challenge of image acquisition biases. In this specific example, attempting to directly apply existing data harmonization methods designed for FL to the TM approach becomes impractical as 51.8% (43 out of 83) of the centers in this PD dataset contribute less than 10 data points, with half of them (21 out of 43) contributing less than five data points for training. Unfortunately, restricted data access is a common issue in medical image analysis, which significantly affects the effectiveness of data harmonization methods, especially when insufficient data are available for training. Thus, the current state-of-the-art techniques in this domain may severely counterbalance the advantageous capability of the TM to empower centers to contribute smaller amounts of data effectively in a distributed learning scenario.

An alternative to GAN-based and FL-specific data harmonization approaches is to remove domain-specific information (e.g., scanner or acquisition protocol) from the model’s learned feature representations while trying to retain as much essential disease-related information as possible. For example, Dinsdale et al.²⁹ introduced a data harmonization technique for centralized learning trained on neuroimaging data following this idea. Their work employed an adversarial training setup, with domain-specific information being “unlearned” from the features used by the model for the main task (e.g., brain age prediction). In their setup, the adversarial network comprises an encoder for extracting features from input data, a classification head for the main task, and a classification head for the domain task (i.e., scanner). Despite the effectiveness of their work, certain restrictions, such as ensuring that batches include representation from every scanner and oversampling the smallest dataset, must be adapted for a distributed environment where batches are composed of data from a single center. Adapting this method to distributed learning is challenging because centers do not share information, which limits the unlearning step, particularly in the FL setup, where training occurs in parallel. To address this challenge, Dinsdale et al.²³ proposed an adaptation that tracks site information to generate features on the central server. Although sharing these features with each site enables the unlearning procedure to be performed effectively, it contradicts the principle of distributed learning, which aims to avoid data sharing in any way. The TM training approach not only allows centers with limited local data to participate but also eliminates the need to share information or features during the unlearning step. Thus, investigating the adversarial training setup is theoretically viable for the TM approach, although this has not yet been experimentally verified.

Therefore, this work introduces HarmonyTM, which adapts and extends the harmonization framework proposed by Ref. 29 to the TM approach, which is subsequently tested for PD classification. Overall, we aim to learn a feature representation that minimizes the effect of image acquisition biases (i.e., spurious correlations) while retaining high performance for the main task and enabling centers to contribute very small sample sizes. Our major contributions include: (1) the development of the first data harmonization method for the TM approach and (2) the evaluation of HarmonyTM in reducing the impact of scanner differences while successfully distinguishing between patients with PD and healthy participants.

2.

Material and Methods

In this work, we evaluate HarmonyTM using the most extensive multi-center database for PD classification published to date.

2.3.

Harmonization Strategy

The deep learning architecture used in this work is based on the state-of-the-art simple fully convolutional network (SFCN),⁴⁵ which achieved high performance for PD classification using multi-center T1-weighted MRI scans in centralized and TM approaches.^26,46 The model’s encoder comprises seven blocks: The initial six blocks mirror the structure of the original SFCN model. These include five blocks featuring a 3D convolutional layer with $3\times 3\times 3$ kernel filters, batch normalization, $2\times 2\times 2$ max pooling, and ReLU activation, whereas one block incorporates a 3D convolutional layer with $1\times 1\times 1$ kernel filters, batch normalization, and ReLU activation. The last block is tailored for our specific task and includes a 3D average pooling layer, a dropout layer with a 0.2 rate, and a flattening layer with 768 features (neurons). The classification head utilizes the encoder output and consists of a single dense layer using a sigmoid activation function for the binary output, distinguishing between patients with PD and healthy participants. Meanwhile, the domain head is a single layer employing a softmax activation function for the multiclass output, categorizing 23 different scanner types.

Before removing domain-specific details, it is necessary to pre-train the network components with the TM approach. Therefore, the encoder and disease classification head are initially trained until convergence. Following this, the encoder is frozen, and the scanner classification head is trained until convergence. Finally, utilizing these pre-trained models, the scanner harmonization procedure is implemented in three steps as follows. It is important to highlight that these steps are performed for each batch, such that the three training steps occur at each center before transferring the model to the next center (see Fig. 1).

• 1. Optimize the encoder and disease classification head for the PD classification task.

• 2. Optimize the scanner classification head for identifying scanners from the feature representation of the frozen encoder that is trained in step 1.

• 3. Optimize the encoder by employing an adversarial confusion loss to eliminate scanner-specific information. This loss guides the scanner classification head output toward chance-level performance. In essence, chance-level means that the model would make predictions purely by random guessing, eliminating any shortcuts related to scanners that might be exploited for disease classification.

Fig. 1

HarmonyTM steps to remove scanner-specific information from the feature representation of the PD classification model.

The scanner harmonization process involves the application of four distinct losses. In the first step, a PD classification loss [Eq. (1)] is utilized to minimize the binary cross entropy. This consists of adjusting the model’s parameters to diminish the difference between its predictions and the ground truth labels (PD versus HP). In the second step, a scanner classification loss [Eq. (2)] is employed with a similar objective as in step 1. However, categorical cross-entropy is minimized in this instance, given that the prediction involves identifying scanner types and only the domain head is optimized in this step. An adversarial confusion loss [Eq. (3)] is computed in the third step, introducing a counteractive objective to step 2. Step 2 strives for the model to recognize scanner types precisely, whereas step 3 aims to counteract this ability by removing information related to scanner types from the feature representations, where $N$ represents the batch size, which varies depending on the amount of data available at each center. In the end, the total loss [Eq. (4)] is computed as the sum of the three losses described. This cumulative loss is subsequently used to optimize the encoder. Compared with the implementation proposed in Ref. 29, two important adaptations were made. (1) Introducing $N$ to manage varying batch sizes during training. In our case, $N$ accounts for the variation in batch size due to differing amounts of data across centers. (2) Eliminating the requirement for oversampling, which was used in the original method²⁹ to ensure that every scanner type was represented in each batch. Their centralized approach allowed for control over such a batch composition, ensuring the representation of every scanner type and oversampling underrepresented types to achieve balance. However, this is not feasible in the TM approach, where each center only has access to its own data. The harmonization process was iterated through 30 cycles, taking $\sim 2.5\mathrm{h}$ to complete. The code is available in a GitHub repository available at: https://github.com/RaissaSouza/scanner-harmonization:

Eq. (1)

$$Lpd(y,\hat{y})=-\frac{1}{N}\sum _{n=1}^N{y}_n\times \log (y-\widehat{y}),$$

Eq. (2)

$$Lsc(y,\widehat{y})=-\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^M{y}_{i,j}\times \log ({\widehat{y}}_{i,j}),$$

Eq. (3)

$$Lconf(\widehat{y})=-\sum _{n=1}^N\frac{1}{N\times 23}\log ({\widehat{y}}_n),$$

Eq. (4)

$$Ltotal=Lpd+Lsc+Lconf.$$

3.

Results

The results of this work show that the PD classification performance of the TM not only remains stable but also improves after removing scanner-specific information from the feature representation (see Table 2). As expected, the scanner classification performance drops across all metrics after harmonization. Moreover, it can be observed that the TM approach is less prone to encoding scanner information (53% accuracy) before harmonization when compared with the centralized approach (65% accuracy). Nevertheless, in both cases, the harmonization method reduces scanner classification abilities to 30% accuracy and sensitivity, 96% specificity, 12% precision, and 16% $F1$-score. Figure 2 supports these findings as illustrated by the confusion matrices of the centralized and TM approaches before and after harmonization. The matrices display the counts of predictions for each pair of actual and predicted classes, where each row represents the actual scanner types, and each column represents the predicted scanner types. In Fig. 2(b), it is evident that after scanner harmonization, the majority of datasets are mistakenly classified as Siemens Skyra and Siemens Trio Tim. Together, these two scanners constitute the largest portion (37.5%) of the test set (Fig. S1 in the Supplementary Material shows the proportion of each scanner type in the test set). Furthermore, our results demonstrate that the TM approach leads to the greatest improvements in PD classification after scanner harmonization, achieving 76% accuracy, 82% AUROC, 83% sensitivity, and 75% $F1$-score compared with the centralized approach. None of the models showed improvements in specificity and precision after harmonization, although their performance remained comparable.

Table 2

Parkinson’s disease and scanner classification performances before and after scanner harmonization.

Fig. 2

Confusion matrices of the scanner predictions for the centralized and TM approaches. Panel (a) illustrates classification before harmonization, whereas panel (b) presents them after harmonization.

Figures 3 and 4 display the distributions of the first two modes of PCA applied to the feature representations of the models, with color-coded by scanner type in Fig. 3 and color-coded by disease class in Fig. 4. In Fig. 3, it can be seen that before scanner harmonization, the distribution of scanner types within the feature representation is more distinct and less overlapping compared with the corresponding results after harmonization. In contrast, Fig. 4 shows that the distributions of patients with PD and healthy participants are more similar before harmonization but become more distinguishable after harmonization. Logistic regression analysis reveals that, before scanner harmonization, PCA mode 1 encodes information for disease classification, achieving accuracies of 58% and 53% for the centralized and TM approaches, respectively. After harmonization, these percentages increase to 71% and 75%, respectively. Moreover, logistic regression analysis revealed that some information is encoded through PCA modes 1 and 2 for scanner classification, yielding accuracies of 38% and 35% before harmonization, which subsequently decreases to 28% after harmonization. The complete logistic regression analysis results are provided in Tables S1 and S2 in the Supplementary Material.

Fig. 3

2D scatter plots of the first two PCA modes, considering the scanner type. Panel (a) illustrates the scatterplots before harmonization, whereas panel (b) presents the scatterplots after harmonization.

Fig. 4

2D scatter plots of the first two PCA modes, considering the presence of disease. Panel (a) illustrates the scatterplots before harmonization, whereas panel (b) presents the scatterplots after harmonization.

4.

Discussion

In this work, we introduced HarmonyTM, the first data harmonization method specifically designed and evaluated for the TM approach. The results showed that HarmonyTM is effective in creating a feature representation that reduces image acquisition biases (i.e., spurious correlations) and enhances disease-related information, achieving the highest classification performance for PD after scanner harmonization. We specifically focused this first feasibility analysis on scanners because of the well-known bias effects.⁹ Although our study demonstrated the effectiveness of HarmonyTM in harmonizing imaging data from different scanners using T1-weighted MRI scans, the method can be readily adapted to address other potential spurious correlations and different imaging modalities. However, further experimental validation is needed to confirm its broader applicability.

It is crucial to highlight that the TM approach itself was found less prone to learning spurious correlations than the centralized approach even before data harmonization. This resistance may be attributed to the fact that each training data batch originates from a single center in the TM setup. By avoiding the inclusion of data from multiple centers in a single batch, the likelihood of the model learning image acquisition biases associated with different factors (e.g., scanner types) across centers is reduced. Conversely, the TM approach’s use of single-center batches may explain why it achieved less balanced metrics compared with the centralized approach after harmonization. The centralized approach, which includes data from multiple centers in each batch, increases the likelihood of training in batches where both classes are represented. This contrasts with the TM approach, where training batches from centers such as OASIS, SALD, and ADNI consist only of healthy participants, and centers such as UOA, RUH, and some from PPMI include only PD cases.

Our results indicate that before scanner harmonization, the difference in the distribution of patients with PD and healthy participants within the feature representation of the last encoder’s layer was poorly noticeable in the centralized as well as TM approach. Conversely, the distributions of scanner types exhibited more noticeable differences, allowing for some degree of classification, with a higher level of accuracy of the logistic regression model from the two PCA modes computed from the model trained using the centralized approach compared with the model trained using the TM approach. These observations are consistent with previous research,² which demonstrated that the feature representation of a PD classifier trained with the centralized approach can be used directly for classifying scanner types, achieving even better performance in scanner type classification compared with PD classification. Following scanner harmonization, a shift occurred: the distinctions between patients with PD and healthy participants became more pronounced, whereas the differences in scanner-type distributions became less prominent. This pattern aligns with the results reported by Dinsdale et al.,²⁹ who initially introduced this method for the centralized training approach. Although their research implies that removing image acquisition biases could lead to a loss of disease-related information when the entire dataset is used for harmonization, we chose to apply harmonization to the full training dataset. As a result, we did not impose restrictions based on data availability, which could substantially reduce center participation in real-life scenarios. Despite these concerns, we found that HarmonyTM benefits from using the full dataset for harmonization. Rather than negatively impacting performance, this approach led to improved disease classification outcomes.

It is essential to highlight some of the limitations of this work. First, we utilized a single established PD classifier model, which was trained using a realistic multi-site database. However, only one specific but widely used image modality (T1-weighted MRI) was investigated in this work for harmonization. Therefore, the good performance of HarmonyTM in other scenarios involving different deep learning architectures, other medical imaging modalities, and additional medical imaging tasks has yet to be demonstrated. It is crucial to recognize that biases must be first identified in the data, as it is impossible to unlearn unknown biases. It is noteworthy, however, that our work utilized the largest PD multi-center database, encompassing various scanner types and small datasets from many centers, compared with the datasets used in many other multi-center analyses, potentially enhancing the generalizability of results. Moreover, our investigation was limited to a single neuroimaging modality—T1-weighted MRI data. Although HarmonyTM can be applied to any two-dimensional (2D) or 3D data, the results may vary depending on the degree of spurious correlations present in the specific image modality. Furthermore, although we used a computer equipped with an NVIDIA GeForce RTX 3090 GPU, we believe that any center with comparable resources should be capable of training HarmonyTM. Last, it is important to note that we used only one training epoch at each center per cycle to minimize the risk of overfitting. However, we did not assess the performance at each center to ensure that overfitting did not occur. Nevertheless, additional strategies, such as regularization and the use of an external validation set, could be implemented to further analyze and prevent overfitting. Future work could address several limitations identified in this study. These include exploring HarmonyTM’s effectiveness across different deep learning architectures, medical imaging modalities, and diagnostic tasks to assess its generalizability. In addition, further investigation into strategies for mitigating overfitting risks and determining the minimal computational resources required would enhance its viability for under-resourced centers.

5.

Conclusion

This work introduced HarmonyTM, a method specifically designed for harmonizing 3D MRI data in the context of the TM approach. HarmonyTM tackles the issue of image acquisition biases across different centers, a common challenge in systems that learn from data distributed across multiple locations. To the best of our knowledge, this is the first work implementing a data harmonization method for the TM approach. Our findings demonstrate the effectiveness of HarmonyTM in generating features with reduced influence from image acquisition biases, such as scanner types, while not only maintaining but also improving performance in classifying PD. In addition, our results emphasize that the TM approach has inherent resistance to learning image acquisition biases. This aspect is crucial for developing clinically useful deep learning models with broad applicability.