2.1.

Data

DWI scans from two sites were combined for joint analysis: biomarkers of cognitive decline among normal individuals: the BIOCARD cohort (BIOCARD)²⁷ and Vanderbilt Memory and Aging Project (VMAP).²⁸ BIOCARD patients were scanned on a 3T Philips Achieva scanner (Eindhoven, The Netherlands) at Johns Hopkins University in Baltimore, Maryland, United States. Diffusion-weighted images were acquired from a spin echo sequence (TR = 7.5 s, TE = 75 ms, $\text{resolution}=0.828\mathrm{mm}\times 0.828\mathrm{mm}\times 2.2\mathrm{mm}$, $b\text{-}\text{values}=0$, $700\mathrm{s}/{\mathrm{mm}}^2$, and number of gradients = 33). VMAP patients were scanned on a Philips 3T Achieva scanner (Best, The Netherlands) at Vanderbilt University in Nashville, Tennessee, United States. Diffusion-weighted images were acquired from a spin echo sequence (TR = 8.9 ms, TE = 4.6 ms, $\text{resolution}=2\mathrm{mm}\times 2\mathrm{mm}\times 2\mathrm{mm}$, $b\text{-}\text{values}=0$, $1000\mathrm{s}/{\mathrm{mm}}^2$, and number of gradients = 32).

Data used in this study are split into training and testing cohorts. The training cohort comprised 347 scans from BIOCARD (site 1) and 324 scans from VMAP (site 2), all free of cognitive impairment. Training data from VMAP have ages $74.1\pm 7.5$ and 114 women. BIOCARD training data are ages $73.1\pm 5.7$ with 205 women. The testing cohort is 77 matched participants (one scan from each participant), free of cognitive impairment, ages $72.9\pm 7.6$, and 57% percent women.

To isolate site-wise differences between these two datasets without traveling subjects,^29,30 we curate a subset of patients that have the same demographics. The matching process finds a 1:1 cross-site mapping of each patient based on demographics. For example, a 59-year-old female from the VMAP cohort is matched with a 59.5-year-old female from the BIOCARD cohort. Matched patients are the same sex and $\pm 1$ year in age. The matching was done using the pyPheWAS maximal group matching tool (version 4.1.1).³¹

2.2.

Diffusion Processing

The proposed model learns from connectome representations of diffusion tractography and is derived from DWI outlined in Sec. 2.1. DWI from all participants were first preprocessed to remove eddy current, motion, and echo-planar imaging distortions prior to any model fitting.³² We used MRTrix³³ to perform tractography over fiber orientation distributions with anatomically constrained tractography framework (seeded on gray matter–white matter interface, allowed backtracking, terminating with five-tissue-type mask, and generated 10 million streamlines³⁴). Afterward, we map the tractogram to a connectome representation using the Desikan–Killany atlas³⁵ with 84 cortical parcellations from Freesurfer.³⁶ We use the Brain Connectivity toolbox (version-2019-03-03) to compute 12 brain network measures⁵ and filter connections less than 0.00001% of the total number of streamlines.³⁴ Network measures computed with this toolbox are the site-biased ground truth used in model training.

We consider modularity, average node betweenness centrality, assortativity, average node participation coefficient, average node clustering, average node strength, average local efficiency, global efficiency, density, rich club coefficient, characteristic path length, and number of edges in the characteristic path length (“edge count”).⁵ Modularity is the quality of division of the network into modules.⁵ Betweenness centrality is the sum of the ratio between the shortest paths that contain the node and the total number of the shortest paths.⁵ Assortativity is the correlation coefficient among the degrees of all nodes on two opposite ends of a connection and reflects network resilience.⁵ Participation coefficient characterizes the diversity of intermodular connections of nodes.⁵ Clustering coefficient is the possibility that any two neighbors of a node are also connected.⁵ Here, node strength is the number of streamlines connecting a node.⁵ Global efficiency is the ability for information to move around the whole network, and local efficiency is the ability of information to move around nodal subsets.⁵ Density is the number of present connections to the total number of connections possible.⁵ Rich club networks have subgroups of central nodes that tend to interact with one another.³⁷ The characteristic path length is the average shortest path, in millimeters, to traverse the network.⁵

2.3.

Model Architecture

We implemented a variational autoencoder (VAE) with site-conditional restrictions on the latent space, $z$. The proposed model architecture is outlined in Fig. 1. First, the encoder processes the input data using two layers of linear layers and two layers of non-linear activation (ReLU). Following encoding, the latent space is reparametrized by sampling from the learned distribution. Next, the decoder takes this compact representation and reconstructs it back into the original data format while accounting for site-specific characteristics. In addition to recreating the original connectome, the model is trained to predict other useful information, such as the patient’s sex, age, and brain connectivity features. These predictions help ensure that the latent space is organized in a way that captures meaningful patterns and relationships in the data. During test time, we set the decoding conditions to project network measures and the reconstructed connectome to a site domain of the user’s choice.

We chose this architecture to learn site-minimal, bio-maximal representations of the connectome. We accomplish this goal by maximizing mutual information ($I(x,y)$) of $z$ and biological attributes (let $a$ be patient age and $s$ be patient sex) and $z$ with network connectivity properties ($B$). Simultaneously, we minimize mutual information of $z$ with site (let $c$ be the site encoding). We translate these relationships to loss functions as follows: to maximize mutual information, we minimize loss [mean squared error (MSE) or binary cross entropy]. To minimize mutual information between $z$ and site, we minimize MSE loss of the conditional reconstruction [the first term in Eq. (1)] and maximize KL divergence [the second term in Eq. (1)].

The second term in this bound is difficult to compute directly but, as shown in Moyer et al.,²⁰ can be approximated using the standard normal Gaussian in place of induced marginal $q(z)$. This fits nicely with existing VAE literature,³⁸ is computationally tractable, and, as we show, performs well empirically for removing site information. We use Eq. (1) to replace $I(z,c)$ in Eq. (2).

Therefore, the overall loss function, ${\ell }_{\text{total}}$, is the sum of five sub-component losses: connectome reconstruction (${\ell }_{\mathrm{recon}}$), mean squared site-conditional prediction error for brain network measures for BIOCARD (${\ell }_B(c)$, $c=1$) and VMAP (${\ell }_B(c)$, $c=2$), mean squared age prediction error (${\ell }_{\text{age}}$), and binary cross entropy loss of sex predictions (${\ell }_{\text{sex}}$). The final sub-component loss is KL divergence, which is a key component in our learning as it limits site information in $z$.

Toward that end, our component loss function can be rewritten as mutual information terms (up to constant entropic terms):

2.4.

Co-Learning and Individual Learning Schemes

We explore two learning schemes for optimizing the latent space. In the first scheme, the model learns to predict one brain network measure from the site-invariant latent space (“individual learning”). We then train 12 separate models, one for each brain network measure. In the second scheme, the model learns to predict all 12 brain network measures at the same time from the site-invariant latent space (“co-learning”).

2.5.

Bootstrapping

To assess training stability, we bootstrapped training data. At each bootstrap iteration, 80% of training data are randomly sampled without replacement. The same testing data are used to evaluate each iteration.

2.6.

Evaluating Harmonization Performance

We measure successful harmonization by evaluating differences in the connectomes and network measures with respect to their site variable. Although we do not use traveling participants in this study, we do match patients across sites based on their age, sex, and cognitive status. We expect network measures computed for this matched dataset to have similar distributions. We measure the disparity between sites with Cohen’s $D$ and the Mann–Whitney $U$-test of medians (with $p\text{-}\text{value}<0.05$ as significant). Cohen’s $D$ is a standardized effect size that reports the difference between two means compared with the overall data variability.³⁹ Successful harmonization would reduce large (Cohen’s $D>0.5$) differences to small (Cohen’s $D<0.2$) and remove significant differences in the medians among sites.³⁹

The second benchmark of harmonization performance is preserving biological variation. Previous studies^40–42 demonstrated that connectomes contain relevant signals related to patient age and sex. During model training, the latent space is shaped by five tasks. Two of which, sex and age prediction, are included to encourage the model to retain relevant biological covariates.