Year: 2021

Abstract:
A plethora of machine learning methods have been applied to imaging data, enabling the construction of clinically relevant imaging signatures of neurological and neuropsychiatric diseases. Oftentimes, such methods do not explicitly model the heterogeneity of disease effects, or approach it via nonlinear models that are not interpretable. Moreover, unsupervised methods may parse heterogeneity that is driven by nuisance confounding factors that affect global brain structure or function, rather than heterogeneity relevant to a pathology of interest. On the other hand, semi-supervised clustering methods seek to derive a pathology-related dichotomous subtypes, ignoring the truth that disease heterogeneity spatially and temporally extends along a continuum. To address the aforementioned limitations, herein, we propose a novel method, termed Surreal-GAN (Semi-SUpeRvised ReprEsentAtion Learning via GAN). Using cross-sectional imaging data, Surreal-GAN dissects underlying disease-related heterogeneity under the principle of semi-supervised clustering [cluster mappings from the normal control (CN) to patient (PT) domain], proposes a continuously dimensional representation, and infers the disease severity of patients at individual level along each dimension. The model first learns a transformation function from the CN domain to the PT domain with latent variables controlling transformation directions. An inverse mapping function together with regularization on function continuity, pattern orthogonality and monotonicity was also imposed to make sure that the transformation function captures necessarily meaningful imaging patterns with clinical significance. We first validated the model through semi-synthetic experiments, and then demonstrated its potential in capturing biologically plausible imaging patterns in Alzheimer’s disease (AD).

heading: 1 INTRODUCTION
text: Neuroimaging, conventional machine learning, and deep learning have offered unprecedented opportunities to understand the underlying mechanism of brain disorders over the past decades (Davatzikos (2018)), and pave the road towards individualized precision medicine (Bzdok & MeyerLindenberg (2018)). A large body of case-control studies leverage mass univariate group comparisons to derive structural or functional neuroimaging signatures (Habeck et al. (2008), Hampel et al. (2008), Ewers et al. (2011)). However, these studies suffer from under-powered statistical inferences since the heterogeneous nature of neurological diseases usually violate the assumption that each group population are pathologically homogeneous. A developing body of methodologies aim to disentangle this heterogeneity with various machine learning methods. Zhang et al. (2016) used the Bayesian Latent Dirichlet Allocation (LDA) model to identify latent atrophy patterns from voxel-wise gray matter (GM) density maps. However, due to the nature of the LDA model, the density maps need to be first discretized, and the model hypothesized to exclusively capture brain atrophy while ignoring the potential regional enlargement. Young et al. (2018) proposed another method, Sustain, to uncover temporal and phenotypic heterogeneity by inferring both subtypes and stages. However, the model only handles around 10-20 large brain regions without more detailed information. Notably, both methods applied directly in the PT domain, confront main limitations in avoiding potential disease-unrelated brain variations. In con-
trast, semi-supervised clustering methods (Varol et al. (2016), Dong et al. (2015), Wen et al. (2021), Yang et al. (2021)) were proposed to cluster patients via the patterns or transformations between the reference group (CN) and the target patient group (PT). A recent proposed deep learning-based semi-supervised clustering method, termed Smile-GAN (Yang et al. (2021)), achieved better clustering performance by learning multiple mappings from CN to PT with an inverse clustering function for both regularization and clustering. Albeit these methods demonstrated potential in clinical applications, they have a major limitation. In particular, these methods model disease heterogeneity as a dichotomous process and derive a discrete output, i.e., subtype, ignoring the fact that brain disorders progress along a continuum. Moreover, variations among subtypes were contributed by both spatial and temporal differences in disease patterns, thus confounding further clinical analysis.(Fig. 1a) To address the aforementioned limitations, we propose a novel method, Surreal-GAN (SemiSupervised Representation Learning via GAN), for deriving heterogeneity-refined imaging signatures. Surreal-GAN is a significant extension of Smile-GAN under the same principle of semisupervised learning by transforming data from CN domain X to PT domain Y. Herein, the key extension is to represent complex disease-related heterogeneity with low dimensional representations with each dimension indicating the severity of one relatively homogeneous imaging patterns (Fig. 1b). We refer to these low dimensional representations as R-indices (ri, where i represent the ith dimension) of disease pattern. The first key point of the method is modelling disease as a continuous process and learning infinite transformation directions from CN to PT, with each direction capturing a specific combination of patterns and severity. The idea is realized by learning one transformation function which takes both normal data and a continuous latent variable as inputs and output synthesized-PT data whose distribution is indistinguishable from that of real PT data. The second key point is controlling the monotonicity (from light to severe with increasing index value) of disease patterns through one monotonicity loss and double sampling procedure. The third key point is boosting the separation of captured disease patterns through one orthogonality regularization. The fourth key point is introducing a inverse mapping function consisting of a decomposition part and a reconstruction part which further ensure that patterns synthesized through transformation are exactly captured for inferring representations of disease patterns, the R-indices. Besides the above mentioned regularization, function continuity, transformation sparsity and inverse consistency further guide the model to capture meaningful imaging patterns. To support our claims, we first validated Surreal-GAN on semi-synthetic data with known ground truth of disease patterns and severity. We compared performance of Surreal-GAN to NMF (Lee & Seung (1999)), Factor Analysis and LDA (Blei et al. (2001)) models, and subsequently studied the robustness of model under various conditions. Finally, we applied the model to one real dataset for Alzheimer’s disease (AD), and demonstrated significant correlations between the proposed imaging signatures and clincial variables.

heading: 2 METHODS
text: The schematic diagram of Surreal-GAN is shown in Fig. 1d. The model is applied on regional volume data (volumes of different graymatter (GM) and whitematter (WM) regions) derived from structural MRI. The essential idea is to learn one transformation function f : X ∗ Z → Y which transforms CN data x to different synthesized PT data y′ = f(x, z), with latent variable z specifying distinct combinations of patterns and severity. Here, the Latent (LAT) domain, Z = {z : 0 ≤ zi ≤ 1 ∀1 ≤ i ≤ M}, is a class of vectors with dimension M (predefined number of patterns). We denote four different data distributions as x ∼ pcn(x), y ∼ ppt(y), y′ ∼ psyn(y′) and z ∼ plat(z), respectively, where z ∼ plat(z) is sampled from a multivariate uniform distribution U [0, 1]M , rather than a categorical distribution as in Smile-GAN. A continuous z variable lay the foundation for learning continuous representations. In addition to the transformation function, an adversarial discriminator D is introduced to distinguish between real PT data y and synthesized PT data y′, thereby ensuring that the synthesized PT data are indistinguishable from real PT data (i.e. minimizing the distance between real PT data distribution and synthesized PT data distribution). A continuous latent variable does not necessarily lead to desired properties of CN to PT transformation (e.g. separation/monotonicity of captured patterns). That is, there are a number of functions potentially achieving equality in distributions, thus making it hard to guarantee that the transformation function learned by the model is closely related to the underlying pathology progression process. Moreover, during the training procedure, the transformation function backboned by the neural network tends to ignore the latent variable z. Even when the Latent variable z is emphasized
through regularization, different components of z variable do not naturally lead to separate disease patterns and are not positively correlated with pattern severity. Therefore, with the assumption that there is one true underlying function for real PT data y = h(x, z), Surreal-GAN aims to boost the transformation function f to approximate the true underlying function h, while having the latent variable, z, contributing to recognizable imaging patterns rather than being random ”activators” of the GAN. For this purpose, we constrain the function class of f via different types of regularization as follows: 1) encourage sparse transformations; 2) enforce Lipschitz continuity of functions; 3) introduce one inverse function g : Y → Z serving to decompose and reconstruct; 4) boost orthogonality among synthesized/captured patterns; 5) enforce monotonicity of synthesized patterns and their positive correlations with components of z. The training procedure of Surreal-GAN is thus a adversarial learning process constrained by these regularization.

heading: 2.1 ADVERSARIAL LOSS
text: The adversarial loss (Goodfellow et al. (2014)) is used for iterative training of discriminator D and the transformation function f , which can be written as:
LGAN(D, f) = Ey∼ppt(y)[log(D(y))] + Ex∼pcn(x),z∼Plat(z)[log(1−D(f(x, z)))] (1) = Ey∼ppt(y)[log(D(y))] + Ey′∼psyn(y′)[log(1−D(y ′))] (2)
The transformation function f attempts to transform CN data to synthesized PT data so that they follow similar distributions as real PT data. The discriminator, providing a probability that y comes from the real data rather than transformation function, tries to identify the synthesized PT data and distinguish it from real PT data. Therefore, the discriminator D attempts to maximize the
adversarial loss function while the transformation function f attempts to minimize against it. The training process can be denoted as:
min f max D
LGAN(D, f) = Ey∼ppt(y)[log(D(y))] + Ey′∼psyn(y′)[log(1−D(y ′))] (3)

heading: 2.2 REGULARIZATION

heading: 2.2.1 SPARSE TRANSFORMATION
text: We assumed that disease process will not change brain anatomy dramatically and primarily only affect certain regions throughout most of disease stages, which means that the true underlying transformation only moderately changes some regions while keeping the rest mildly affected with most of z variable values. Therefore, to encourage sparsity and control the distance of transformations, we defined one change loss to be the l1 distance between the synthesized PT data and the original CN data:
Lchange(f) = Ex∼pcn(x),z∼plat(z)[||f(x, z)− x||1] (4)

heading: 2.2.2 LIPSCHITIZ CONTINUITY AND INVERSE MAPPING
text: First, by enforcing the transformation function f to be K1-Lipschitz continuous, we can derive that, for fixed latent variable z = a and ∀x1,x2 ∈ X, ||f(x1,a)− f(x2,a)||2 ≤ K1||x1−x2||2. Thus, by controlling the constant K1, we constrained the transformation function to preserve the original distances among CN data instead of scattering them into the PT domain if a same combination of disease patterns is imposed to them. Second, to prevent latent variable z from being ignored, we introduced an inverse mapping function from PT domain to Latent domain, g : Y → Z. By Lemma1 (section A.1), with function g being K2-Lipschitz continuous and d(·, ·) being any distance metric satisfying triangular inequality, we can derive that, ∀z1, z2 ∼ plat(z), z1 6= z2, and x̄ ∼ pcn(x), d(f(x̄, z1), f(x̄, z2)) is lower bounded by (d(z1,z2)K2 − 1 K2
(d(g(f(x̄, z1)), z1) + d(g(f(x̄, z2)), z2))). Therefore, by minimizing the distance between sampled latent variable z and reconstructed latent variable g(f(x, z)), we can control differences between synthesized patterns to be positively related to the distances between latent variables, and thus to be non-trivial (i.e., same CN data is transformed to significantly different PT data with different z variables). Therefore, we define a reconstruction loss as the l2 distance between g(f(x, z)) and z:
Lrecons(f, g) = Ex∼pcn(x),z∼plat(z)[||g(f(x, z))− z||2] (5)

heading: 2.2.3 PATTERN DECOMPOSITION
text: Nevertheless, a simple function g directly reconstructing sampled z from synthesized PT y′ do not reconstruct each component zi completely based on changes led by zi in transformation. We are more interested in pattern representations of real PT data than generating synthesized PT, and inverse function g will be used for inferring R-indices as introduced in section 2.4. Thus, it is very important for the inverse function to accurately capture patterns synthesized by intensively regularized transformation function f . For this purpose, we further separated the function g : Y → Z into g1 : Y → RM∗S (where M is number of patterns and S is number of input ROIs) and g2 : RS → R. Decomposer g1 serves to reconstruct changes synthesized by each component zi in transformation process: qi = f(x,ai) − x, where ai is a vector with ith component aii = zi and aij = 0, ∀i 6= j. Let q̂f(x,z) = [qT1 , · · · ,qTM ]T be concatenation of all synthesized changes qi, we defined the decomposition loss as:
Ldecom(f, g1) = Ex∼pcn(x),z∼plat(z)[||g1(f(x, z))− q̂f(x,z)||2] (6)
g2 serves to further reconstruct each component of sampled z variable from g1(f(x, z)) independently, so that g(x) = l̂g2(g1(f(x,z)) = [g2(g1(f(x, z))0:S), · · · , g2(g1(f(x, z))S∗(M−1):S∗M )]T and the reconstruction loss can be rewritten as
Lrecons(f, g) = Lrecons(f, g1, g2) = Ex∼pcn(x),z∼plat(z)[||l̂g2(g1(f(x,z)) − z||2] (7)
Function g here can be also considered as the approximation of expectation of posterior distribution. In this sense, minimization of l2 distances can also be interpreted as maximization of mutual information between latent variable z and synthesized PT data y′ = f(x, z) as explained in Remark1 in Appendix (Chen et al. (2016)). From this perspective, the transformation function is forced to best utilize information from latent variable z, while also keeping mutual information between synthesized PT data y′ and original CN data x by controlling transformation distances.

heading: 2.2.4 ORTHOGONALITY OF PATTERNS
text: Inclusion of latent variable z in the transformation function does not guarantee that different components of the latent variable, zi, contributes to relatively different patterns in synthesized PT data. Instead, during training procedure, different components, zi, are more likely to synthesize patterns in same regions and lead to accumulated severity. Therefore, we define another orthogonal loss to encourage separation among patterns synthesized by different components. With changes led by each component, qi defined in section 2.2.3, we constructed a matrix Af(x,z) with the ith column Af(x,z):,i = |qi|. To encourage separation of changes led by different components, we boosted the matrix Af(x,z) to be relatively orthogonal by minimizing the following orthogonality loss:
Lortho(f) = Ex∼pcn(x),z∼plat(z)[||A T f(x,z)Af(x,z) − I||F ] (8)

heading: 2.2.5 MONOTONICITY AND POSITIVE CORRELATION
text: Besides separation of synthesized patterns, a positive correlation between values of each component zi (from 0 to 1) and severity of synthesized patterns can not be acquired spontaneously either. Regarding severity of the disease pattern, we assumed that, as a pattern is becoming more severe, absolute values of changes in regional volumes can only increase monotonically or remain constant. In another word, there is no oscillation in regional volumes as disease pattern is becoming more severe, while more regions can be included gradually with absolute changes switched from 0 to a positive value. To satisfy this requirement, we sampled another latent variable z′ ∼ psev(z′|z), conditioned on previous sampled z variable, such that z′i ≥ zi, ∀1 ≤ i ≤ M . With these double sampled latent variables, we defined the monotonicity loss to be:
Lmono(f) = Ex∼pcn(x),z∼plat(z),z′∼psev(z′|z)[||max(|f(x, z)− x| − |f(x, z ′)− x|,0)||2] (9)
Minimization of this monotonicity loss penalize the case that volume changes induced by z is larger than that induced by z′, while the other direction is permitted. To further ensure a positive correlation (small z lead to mild patterns; large z lead to severe patterns), we introduced a cn loss, which serves to more strictly constrain the distance between synthesized PT data y′ and CN data x, when latent variable z is close to 0 vector. By letting pcn(z) = U(0, 0.05)M to be a multivariate uniform distribution, the cn loss is defined as:
Lcn(f) = Ex∼pcn(x),zcn∼pcn(z)[||f(x, z cn)− x||1] (10)

heading: 2.3 FULL OBJECTIVE
text: With the adversarial loss and all other loss functions defined for regularization, we can write the full objective as:
L(D, f, g1, g2) = LGAN(D, f) + γLchange(f) + κLdecom(f, g1) + ζLrecon(f, g1, g2)
+ λLortho(f) + µLmono(f) + ηLcn(f) (11)
with γ, κ, ζ, λ, µ and η be parameters controlling the relative importance of different loss functions. Through training process, we want to derive parametrized functions f ,g1 and g2 satisfying:
f, g1, g2 = arg min f,g max D L(D, f, g1, g2) (12)

heading: 2.4 REPRESENTATION-INDICES INFERENCE
text: With the assumption that psyn(f(x, z)) ≈ ppt(y) and f satisfying all constraints, we consider learned transformation function f to be a good approximation of the true underlying function h, such
that f(x, z) ≈ h(x, σ(z)), where σ ∈ Ω and Ω is a class of permutation functions which change orders of elements in vector z. Since the order of indices in derived representation is not important and we can always reorder them to find the best matching, we simply rewrite it as f(x, z) ≈ h(x, z) without loss of generality. For any real PT data, ȳ = h(x̄, r) ∼ ppt(y), we can derive that the ground truth R-indices r ≈ g(f(x̄, r)) ≈ g(h(x̄, r)) ≈ g(ȳ). Therefore, reconstruction function g can be directly applied to infer R-indices of disease patterns, r, of any unseen PT data.

heading: 3 IMPLEMENTATION DETAILS

heading: 3.1 MODEL ARCHITECTURE
text: The general architecture of networks can be understood from Fig. 1c. The transformation function f utilizes a encoding-decoding structure. Latent variable z and CN data x are first mapped to two vectors with dimension 34 and the element-wise multiplication of them are then decoded to construct the synthesized PT data with dimension 139. The inverse function g first decodes real/synthesized PT data y/y′ to larger vector with dimension 139*M and then maps them into a vector with dimension M. The discriminator D utilizes the encoding structure which maps y/y′ to a vector with dimension 2. More details of model architectures can be found in Appendix Table 2 and 3.

heading: 3.2 LIPSCHITZ CONTINUITY
text: In implementation, we performed weight clipping to ensure Lipschitz continuity of transformation function f and inverse function g (Arjovsky et al. (2017)). With Θf and Θg denoting the weight spaces of function f and g, the compactness of them imply the Lispchitz continuity of two functions, with Lipschitz constants K1 and K2 only depends on Θf and Θg . The compactness of Θf and Θg was guaranteed by clapping weight spaces into two closed box, Θf = [−cf , cf ]d and Θg = [−cg, cg]d. In implementation, both cf and cg were set to be 0.5 and these bounds can be further relaxed. Weight clipping was claimed not to be the most satisfactory method to ensure Lipschitz continuity because of difficulty in choosing the appropriate clipping bounds. However, in our case, we performed weight clipping for function f and g rather than discriminator D for a different purpose and weight clipping contributed to good performances for deriving representations.

heading: 3.3 TRAINING DETAILS
text: Regarding optimization procedure, ADAM optimizer(Kingma & Ba (2014)) was used with a learning rate (lr) 4 ∗ 10−5 for Discriminator and 2 ∗ 10−4 for transformation function f and clustering function g. β1 and β2 are set to be 0.5 and 0.999 respectively. For hyper-parameters, we set γ = 6, κ = 80, ζ = 80, µ = 500, η = 6. However, performance of model on different tasks are robust to varying hyper-parameters as shown in section A.7. Parameter λ determines the relative importance of orthogonality loss and controls the degree of overlapping among synthesized patterns. Without prior knowledge on overlapping degree of underlying disease patterns, we set λ to different values in the following experiments and used agreements among multiple trained models to determine the optimal one. Moreover, for all experiments, the batch size was set to be 1/8 of the PT data sample sizes. The model was trained for at least 100000 epochs and saved until the reconstruction loss and the monotonicity loss are smaller than 0.003 and 6 ∗ 10−4 respectively. Higher lr might improve convergence speed in some cases. Detailed training procedure is revealed by Appendix Algorithm1.

heading: 4 EXPERIMENTS

heading: 4.1 DATA PREPROCESSING
text: For semi-synthetic and real data experiments, we used T1-weighted (T1w) MR imaging (MRI) and cognition from an international multicentric consortium for cognitive impairment and dementia (anonymous for review process). All baseline MRIs were corrected for intensity inhomogeneities (Sled et al. (1998)). Then, 139 regions of interest (ROIs) of brain volumes were extracted using a multi-atlas label fusion method (Doshi et al. (2015)). The ROIs were then harmonized (Pomponio et al. (2020)) to remove site effect and utilized as input features for Surreal-GAN. (Fig. 6)

heading: 4.2 SEMI-SYNTHETIC EXPERIMENTS
text: Semi-synthetic data construction To set the ground truth of disease patterns and severity a prior, we generated semi-synthetic data by imposing simulated disease-related patterns to real CN data, thereby retaining the variation from the real CN data. Specifically, we split 1392 cognitively normal (age < 70 years) into one real CN group with 492 subjects and one Pseudo-PT group with 900 subjects. For the ith Pseudo-PT subject, we sampled a three-dimensional vector as pattern severity from a multivariate uniform distribution, si ∼ U [0, 1]3. Three types of atrophy patterns were introduced to 900 Pseudo-PT subjects. Different regions were included in different patterns as shown in Appendix Table 4. Five different datasets were constructed for testing the robustness of model under different conditions: (1) Basic dataset: For the ith Pseudo-PT subject and the jth ROI included in pattern k, we decreased the volume vij depending on sampled severity: vij = vij−vij ∗sik ∗N(1, 0.05)∗0.3. Each pattern consists of changes in 14 ROIs and has four of them shared with other patterns; (2) Large overlapping dataset: Same as the Basic dataset, except that each pattern shares 8 ROIs with others; (3) Scarce regions dataset: Each pattern consists of changes in only 4 ROIs with other parts the same as the Basic dataset; (4)Within-pattern noise dataset: Regional changes incorporated in the same pattern have larger within-pattern variances with vij = vij − vij ∗ sik ∗ N(1,0.2) ∗ 0.3; (5) Mild atrophy dataset: Compared to the basic case, smaller changes were imposed by letting vij = vij − vij ∗ sik ∗N(1, 0.05) ∗ 0.2. Evaluation and Agreement Metric In our experiments, Concordance Index (c-index) was used as the evaluation metric. Specifically, for each pattern, we calculated a c-index between the inferred R-index and the ground truth. The average of M derived C-indices, refereed as pattern-c-index, was used for model evaluations. Moreover, to quantify the pair-wise agreement between two trained models, we proposed another Pattern-agr-index which equals pattern-c-index between R-indices derived by two different models. Experiment Procedure We first validated the relationship between model agreements and their accuracy. On each dataset, we ran the model 10 times with parameter λ = 0.1, 0.2, 0.4, 0.6, 0.8 respectively. With each λ, the average pair-wise pattern-agr-index among 10 models was compared with their average pattern-c-indices. Further, the λ value leading to the best agreement was considered as the best choice, and results from 10 corresponding models were reported as model performances for each task. Moreover, we fixed the optimal λ value for each task and changed other parameters to 0 to understand contributions from each regularization term (section A.8). Lastly, we compared our model with five other models. NMF and Factor Analysis (FA) are well-known for deriving lower-dimensional representation of complex dataset. The variant opNMF (Sotiras et al. (2014)) has shown promise in parsing complex brain imaging data with extra orthogonal constraint. Discriminant-NMF (DNMF) (Lee et al. (2012)) was proposed to learn representations which best classify subgroups. LDA, though known for extracting topics from documents, have been recently applied to uncover heterogeneity of AD disease from neuroimaging data (Zhang et al. (2016)). Data preprocessing for these models were introduced in section A.9.

heading: 4.3 REAL DATA EXPERIMENTS
text: Data Selection We applied Surreal-GAN to an AD dataset. For AD, we defined the CN group (N=850) to be subjects with Mini-mental state examination (MMSE) scores above 29, and the PT group (N=2204) as subjects diagnosed as mild cognitive impairment (MCI) or AD at baseline. Experiments Procedure We estimated covariates effects (age&gender) in the CN group using a linear regression model. Then, all CN and PT subjects were residualized to rule out the covariate effects. The adjusted features were then standardized with respect to CN group. Without ground truth, we selected both the optimal number of patterns, M , and λ parameter by measuring agreement among repetitively trained models. For each combination of M = 2, 3, 4 and λ = 0.1, 0.2, 0.4, 0.6, 0.8, we ran the model 10 times respectively. The M and λ values leading to the highest agreement was considered optimal. Among the 10 corresponding models, the one having the highest mean pair-wise agreements with the other models was used to derive R-indices for all PT subjects. To visualize patterns corresponding to the ith dimension of R-indices, ri, we grouped subjects into three different subgroups: ri < 0.4, 0.4 < ri < 0.7 and ri > 0.7, with rj < 0.4 for all j 6= i. Voxel-wise group comparisons were performed between each subgroup and the CN group via AFNI-3dttest (R.W. (1996)) and GM-tissue map (Davatzikos et al. (2002)) which encodes, for each tissue type, the volumetric changes observed during the registration. Pattern rep-
resentations, R-indices, of PT subjects were further compared with other clinical variables to test their significance. Details of clinical variables and statistical tests can be found in section A.6.

heading: 5 RESULTS

heading: 5.1 RESULTS ON SEMI-SYNTHETIC DATA
text: We first validated that agreements among multiple trained models (Fig. 2b) are good indicators of model accuracy (Fig. 2a) for parameter selections. Besides, from Fig. 2a, we proved the model’s ability in capturing ground truth of patterns and severity under different cases including pattern overlapping and large within-pattern noise. However, the model experienced declining performances when scarce disease-related regions were included in each pattern and showed even worse performances when very mild atrophies are led by disease process. This is expected since, for the mild atrophy case, the simulated atrophy rate is 0-20% of the regional volume which is even smaller than the standard deviation of most regional volumes among the normal population. Finally, SurrealGAN significantly outperformed all five methods on five data sets we constructed (Table 1). However, none of these methods optimally utilize the CN data as a reference distribution, which is one of important limitations preventing them from precisely learning disease-specific pattern representations rather than capturing disease-unrelated confounding variations. The family of DNMF methods capture discriminant representation by minimizing and maximizing within class and among-class variations respectively, thus did not capture heterogeneity within the PT class accurately either. Designed for parsing disease heterogeneity, semi-supervised methods focused on clustering patients into hard categorical subtypes (Varol et al. (2016), Dong et al. (2015), Wen et al. (2021), Yang et al. (2021)). With different end points, these methods (including Smile-GAN) can not be directly
compared with Surreal-GAN. To the best of our knowledge, our model is the first semi-supervised approach for learning disease-related representations, using a healthy reference group as a means for comparison. Therefore, to better understand the model performance, we leveraged the simulated ground truth information (not available in real experiments), derived ”practical upper bounds” for compared unsupervised models, and also computed the gap between Surreal-GAN and supervised methods. (section A.10 and A.11)

heading: 5.2 RESULTS ON REAL DATA
text: On the AD dataset, models show highest agreement when the number of patterns M = 2. From voxel-based comparison results (Fig. 3), we can clearly visualize two patterns related to AD. Each dimension of R-indices, ri, shows a clear positive correlation with the corresponding pattern severity. Pattern1, correlated with r1, shows diffuse cortical atrophies (e.g., bilateral insula, orbitofrontal and frontal poles), while pattern2 reveals focal atrophy in medial temporal lobes. Reproducibility of these two patterns were validated in section A.12. Furthermore, the 2D R-indices were effective and interpretable signatures for disease diagnosis and prognosis, achieving similar performance in AD classification compared to original 139 ROIs (section A.13). Continuity and monotonicity of inferred R-indices further enable us to easily test correlation among them and other clinical variables (Table 5). Both r1 and r2 are highly correlated with CSF-Abeta value and memory dysfunction measured by ADNI-MEM . However, r1 shows stronger correlations with white matter lesion volumes, presence of hypertension, executive function and language functioning, while r2 is more correlated with CSF-Tau, CSF-pTau and presence of APOE-E4 alleles.

heading: 6 CONCLUSION
text: In this study, we have proposed a novel method, Surreal-GAN, for learning representations of underlying disease-related imaging patterns. This model has overcome limitations in previous published semi-supervised clustering methods and shown great performance on semi-synthetic data sets. On AD dataset, the model parse heterogeneous PT data into two concise and clinically informative indices with preserved discriminant power, further proving its huge potential in uncovering disease heterogeneity from MRI data.

heading: A APPENDIX

heading: A.1 LEMMA1
text: Let f : X ∗ Z → Y and g : Y → X be the transformation and inverse mapping function respectively. If g is K-Lipschitz continuous and d(·, ·) is any distance metric satisfying triangular inequality, then for ∀z1, z2 ∈ Z, and ∀x ∈ X, d(f(x, z1), f(x, z2)) is lower bounded by d(z1,z2)
K − 1 K (d(g(f(x, z1)), z1) + d(g(f(x, z2)), z2))
Proof: By triangular inequality, we can derive that:
d(f(x, z1), f(x, z2)) ≥ 1
K d(g(f(x, z1)), g(f(x, z2))) (13)
≥ 1 K (d(z1, z2)− d(g(f(x, z1)), z1)− d(g(f(x, z2)), z2)) (14)
= d(z1, z2) K − 1 K (d(g(f(x, z1)), z1) + d(g(f(x, z2)), z2)) (15)

heading: A.2 REMARK1
text: By minimizing the reconstruction loss Lrecons defined in section 2.2.2 and 2.2.3, we are maximizing a lower bound of the mutual information between latent variable z and synthesized PT data y′ = f(x, z). Considering Q(z|y′) to be an approximation of distribution P (z|f(x, z)) and assuming that it follows a Gaussian distribution N(g(y′),Σ), mutual information denoted by I and entropy denoted by H , we can derive that: I(z, f(x, z)) = H(z)−H(z|f(x, z)) (16)
= Ey′∼psyn(y′)[Ez′∼p(z|y′)[logP (z ′|y′)]] +H(z) (17) = Ey′∼psyn(y′)[DKL(P (·|y ′)||Q(·|y′)) + Ez′∼p(z|y′)[logQ(z′|y′)]] +H(z) (18) ≥ Ey′∼psyn(y′)[Ez′∼p(z|y′)[logQ(z ′|y′)]] (19) = Ez∼plat(z),y′∼psyn(f(x,z))[logQ(z|y ′)] (20)
= Ez∼plat(z),x∼pcn(x)[[logQ(z|f(x, z))] (21)
≈ 1 nm n∑ i=1 m∑ j=1 (− ln |Σ| 2 − 1 2 (zi − f(xj , zi))TΣ−1(zi − f(xj , zi)) + constant)
(22) (21) and (22) follows the Lemma 5.1 in Info-GAN (Chen et al. (2016)) and law of the unconscious statistician (Lotus) Theorem respectively. Therefore, we derive that the mutual information between the latent variable z and synthesized PT data y′ is bounded below by 1nm ∑n i=1 ∑m j=1(− ln |Σ| 2 − 1 2 (zi− f(xj , zi))
TΣ−1(zi− f(xj , zi)) + constant). With a further assumption that Σ is an identity matrix, maximization of this lower bound is equivalent to minimizing the reconstruction loss Lrecons.

heading: A.3 NETWORK ARCHITECTURE
text: Table 2 and Table 3 display network architectures of transformation function f , discriminatorD and reconstruction function g (consisting of g1 and g2).

heading: A.4 ALGORITHM
text: Detailed training procedure of Surreal-GAN is disclosed by Algorithm 1.
Algorithm 1: Surreal-GAN training procedure. lc represents cross entropy loss and ei represents a one hot vector with 1 at ith component. q̂f(x,z), l̂g2(g1(f(x,z)) and Af(x,z) are defined in section 2.2.3 and 2.2.4. while not meeting stopping criteria or reaching max epoch do
for all batches {xi}mi=1,{yi}mi=1 do sample m latent vectors {zi}mi=1 from multivariate uniform distribution with zi ∼ U [0, 1]M
Update weights of discriminator D: Use ADAM to update θD with gradient: ∇θD 1m ∑m i=i[(lc(D(yi), e1) + lc(D(f(xi, zi), e0)))]
sample m latent vectors {z′i}mi=1 conditioned on previously sampled {zi}mi=1, with each z′ij ∼ U(zij , 1] M sample another m latent vectors {zcni }mi=1 with each zcni ∼ U(0, 0.05)M
Update weights of transformation function f : Use ADAM to update θf with gradient: ∇θf 1m ∑m i=i[lc(D(f(xi, zi), e1))) + γ(||f(xi, zi)− xi||1) +
κ||g1(f(xi, zi))− q̂f(xi,zi)||2 + ζ||l̂g2(g1(f(xi,zi)) − zi||2 + λ||ATf(xi,zi)Af(xi,zi) − I||F + µ||max(|f(xi, zi)− xi| − |f(xi, z′i)− xi|,0)||2 + η||f(xi, zcni )− xi||1]
Update weights of function g1: Use ADAM to update θg1 with gradient: ∇θg1 1 m ∑m i=i[||g1(f(xi, zi))− q̂f(xi,zi)||2]
Update weights of function g2: Use ADAM to update θg2 with gradient: ∇θg2 1 m ∑m i=i[||lg2(g1(f(xi,zi)) − zi||2
end end

heading: A.5 REGIONS WITH SIMULATED ATROPHY IN SEMI-SYNTHETIC DATASETS
text: Regions of interest (ROI) included in three patterns (P1-P3) in different semi-synthetic datasets are revealed by Table 4. Small overlapping case indicates the distribution of ROIs in three different dataset as described in section 4.2: (1) Basic dataset; (4)Within-pattern noise dataset and (5)Mild atrophy dataset. Each ROI is a combination of left and right parts of the region, so each check sign represents inclusion of 2 out of 139 regions.

heading: A.6 ASSOCIATIONS WITH CLINICAL VARIABLES
text: Each dimension of R-indices, ri, is compared with different clinical variables for AD disease. For each clinical variable, we only selected out samples who have it available for statistical analyses. ADNI composite scores, including ADNI-EF, ADNI-MEM and ADNI-LAN, measure subjects’ functions in executive functions, memory and language respectively. They are only available among subjects from Alzheimer’s Disease Neuroimaging Initiative (ADNI) study. Digit Symbol Substitution Test (DSST) measures a range of cognition operations including motor speed, attention, and visuoperceptual. For all mentioned clinical scores, a higher value indicates a better performance. CSFAbeta, CSF-Tau and CSF-PTau are three very important hallmarks of AD disease corresponding to amyloid and tau hypotheses. Since measures of these three values have great variance among different studies, only ADNI samples with CSF measures were included for analysis. For dichotomous variables (Apoe-E4 Alleles Carrier & Hypertension), point biserial correlation with corresponding p value was derived between the variable and each pattern dimension, ri. Pearson correlations were derived for other continuous variables.

heading: A.7 HYPER-PARAMETER SELECTION
text: We trained the model on the basic case and the mild changes case (the easy and hard case) with varying hyper-parameters to test the model robustness to hyper-parameter selections. Specifically, we changed each parameter from 50% to 150% of the preset value (introduced in section 3.3) while keeping other parameters fixed. With each set of parameters, we trained the model 10 times and reported average pattern-c-index and average pattern-agr-index as described in section 4.2. From Fig. 4a, we can see that model performances are mostly robust to different choices. For the mild changes case, models, with different η values, do show higher variations in performances. However,
pattern-agr-indices still have good indication of the best choice, and this metric can always be used for parameter selections on different datasets.

heading: A.8 CONTRIBUTION FROM REGULARIZATION
text: Similar to section A.7, we set different parameters to 0 to understand contributions from corresponding regularization terms. Specifically, we set κ = 0 (no decomposition loss), λ = 0 (no orthogonality loss), µ = 0 (no monotonicity loss), ζ = 0 (no cn loss), µ/ζ = 0 (no constraint on monotonicity and positive correlation), κ/λ/µ/ζ = 0 (no extra regularization compared to Smile-GAN framework) respectively with other parameters fixed. With different sets of parameter values, we trained the model 10 times and report average pattern-c-indices. From Fig.4c, we can clearly understand the importance of different regularization terms and visualize how additional regularization terms lead to improved representation performances compared to a pseudo comparable ”Smile-GAN” (Using same regularization as Smile-GAN but with continuous latent variable).

heading: A.9 DATA PREPROCESSING FOR COMPARED MODELS
text: Z-scores of ROIs with respect to the CN group were directly used as input features of FA model. However, for NMF, OpNMF and LDA model, we need to borrow a prior knowledge that atrophy (volume decreasing) is more common in dementia-related brain diseases, and thus, reversed the sign of z-scores, kept positive z-scores, and set negative values to 0. For LDA model, we further discretized the value by first multiplying z-scores by 10 and then mapping them to the closest integer as described in Zhang et al. (2016). For the DNMF model, both CN and PT data and their corresponding diagnosis labels are provided as input, with ROI data first standardized to lie in range 0-1.

heading: A.10 ”PRACTICAL UPPER BOUNDS” OF UNSUPERVISED METHODS
text: As discussed in section 5.1, unsupervised methods were not designed to best utilize the CN data as a reference distribution. Thus, we derived ”practical upper bounds” for compared models on different tasks, with the assumption that these methods have some components capturing the ground truth while having other components representing non-disease related variations. On each semisynthetic test, we ran each method with number of components, M , ranging from 3-10. For M > 3, we intensively searched for the combination of 3 components that best match the ground truth and calculated the pattern c-index. The best result among M = 3 − 10 is considered as an ”practical upper bound” for the model. As shown in Table 6, upper bounds of NMF methods are comparable to results from Surreal-GAN, while other methods still show much worse performances. Note that these are ”practical upper bounds” derived with known ground truth. In real cases, we are not interested in a set patterns that includes the true disease effects but want to find disease-related patterns only. However, it will be very hard to know what the best number of components is and which components are truly disease-related with these compared models. In contrast, Surreal-GAN can directly and accurately output disease-specific representations without indices capturing nondisease related variations, which is one of the key advantages of the proposed method.

heading: A.11 COMPARISONS WITH SUPERVISED METHODS
text: We further compared the model with linear regression (LR) and support vector regression (SVR) with RBF kernel to understand the gap between our model and the most optimal model performances in semi-synthetic tests. A five fold cross validation were run with three different models. Average pattern-c-indices on both train and test sets are reported. We can first observe that Surreal-GAN, as an semi-supervised methods, shows almost equal performances on train and test set, proving a good ability in generalization. On five different tasks (table7), there is a 0.04-0.07 pattern-c-indices gap between Surreal-GAN and supervised model performances. This gap is considered reasonable but there is still space for model improvements.

heading: A.12 GENERALIZATION AND REPRODUCIBILITY ON REAL DATASET
text: We randomly half split the real data (described in section 4.3) in train and test set, with both CN and PT data equally distributed in two sets. Following same experiment procedure described in section 4.3, we first trained the model on the training set, and then applied the model to derive 2-
dimensional R-indices for all CN and PT data in the test set. Voxel-wise group comparisons were conducted between subgroups and CN in the test set for visualization. From Fig. 5, we can observe that R-indices are correlated with almost same patterns shown in Fig. 3, though much smaller sample size led to relatively lower statistical significance.

heading: A.13 LOW DIMENSIONAL REPRESENTATION FOR AD DIAGNOSIS
text: We used R-indices of test data introduced in section A.12, and compared their ability to classify CN and PT with that of 139 ROI volumes. Support vector machine (SVM) was selected as the classification model. 20% leave-out experiments were performed 100 times with derived R-indices and 139 ROI as input features respectively. In each run, we randomly sampled 80% of data as training set and applied trained models on the rest 20% to derive the AUC metric. Grid search were used for parameter selection. With 2D R-indices as features, SVM model was able to reach an average AUC, 0.895 ± 0.011, which is only mildly lower than AUC derived with 139 ROIs as features, 0.929±0.015. Therefore, we can conclude that, 2D R-indices, kept almost all discriminant information among PT data, while parsing the heterogeneity in a clear, concise and interpretable way.

Review summaries:
In this paper, the authors pointed out the lack of a good method to predict disease progression in a continuum based on representations explicitly modeling specific disease effects. To tackle this problem, they proposed a novel method called Surreal-GAN, a tailored version of generative adversarial networks, to learn separate representations of neurological and neuropsychiatric diseases in an unsupervised manner. The method frames the diseased brains as multiple disease-related features at various severity imposed on normal brains, and by learning these features and severities, it can generate a representation with improved reliability and explanability.

In this paper, the authors present a method for learning a representation of disease-related image patterns from regional volume information generated from structural MRI images. The method is based on a GAN approach, and different regularization methods on the latent representation are presented to improve the learning of disease-related representations. There are four contributions presented in the paper: 1. Modelling a disease as a continuous process 2. Controlling the monotonicity of the disease pattern 3. Separation of captured disease pattern 4. Inverse Mapping to ensure that disease patterns synthesized are captured  The evaluation of the method is performed on semi-synthetic data sets and an actual Alzheimer’s disease data set. The results show that the method can identify two clinically informative patterns. 

Surreal-GAN aims to create fake pathological data with a latent variable and a healthy input and includes an inverse function that predicts the latent variable from a fake/real pathological data. However, unlike its predecessor model Smile-GAN, this adds 5 regularizations (L1 loss between generated output and input, L2 loss between the latent variable and the reverse generated latent variable, decomposing the latent variable, making the latent variables orthogonal to each other, and adding positive correlation to disease pattern and the latent variable) to make the latent variable more interpretable and continuous. This would increase the interpretability since instead of clustering the real data it could be used as understanding how severe the patient’s disease progress is. The results show that the model has a higher c-index than NMF, LDA and FA model predictions. It also shows that the best number of patterns for the model to predict is two, and it indeed has positive correlation to the severity of the disease.