Corentin Ambroise
email: corentin.ambroise@cea.fr

Antoine Grigis
email: antoine.grigis@cea.fr

Josselin Houenou
email: josselin.houenou@aphp.fr

Vincent Frouin
email: vincent.frouin@cea.fr

INTERPRETABLE AND INTEGRATIVE DEEP LEARNING FOR DISCOVERING BRAIN-BEHAVIOUR ASSOCIATIONS WITH STABILITY ANALYSIS PREPRINT

whether they are published or not. The documents may come

Interpretable and integrative deep learning for discovering brain-behaviour associations with stability analysis

Introduction

Today, psychiatry, in both its diagnostic and therapeutic dimensions, is moving from a paradigm based on the study of syndromes to a new one built on an understanding of their underlying neurobiological mechanisms. Achieving this goal of precision psychiatry is an ongoing process that requires the combination of scientific advances, technological innovations, and a patient-centered approach [Williams, 2022]. Identifying relationships between behaviours and brain measures is a key aspect of this paradigm. However, behaviours are complex and often result from a combination of genetic, environmental, and psychological factors. In addition, many behaviours observed in mental health conditions can be present in multiple disorders, and the same disorder may manifest with different behaviours in different individuals. This complexity challenges the idea that a single behaviour signature can correspond to a sole disorder. The NIMH Research Domain Criteria (RDoC) [Insel et al., 2010] provides recommendations for properly addressing this complexity. In contrast to traditional approaches that aim to find a diagnosis from a specific score or modality (as outlined in the Diagnostic and Statistical Manual of Mental Disorders (DSM) [START_REF] Association | Diagnostic and statistical manual of mental disorders : DSM-5â[END_REF]), the RDoC promotes dimensional and transdiagnostic approaches. These dimensions include genetic, biological, environmental, and lifestyle factors in the research on personalised psychiatry. In recent years, the transdiagnostic literature in psychopathology has developed along these lines [START_REF] Fusar-Poli | Transdiagnostic psychiatry: a systematic review[END_REF]. Transdiagnostic studies propose to examine a general psychopathological factor, the so-called p-factor [START_REF] Caspi | The p factor: One general psychopathology factor in the structure of psychiatric disorders?[END_REF], Caspi et al., 2023]. This p-factor underlies mechanisms common to several psychiatric syndromes and is usually represented as a single dimension calculated on the basis of different symptoms. The goal of such studies is to search for neural correlates of the p-factor using various biological markers, such as imaging or genetics.

For these novel experimental studies in populations, new methods are being developed to integrate multimodal data, including structural or functional characteristics of the brain, tabular data from report forms, genotyping, or lifestyle conditions. At the same time, cohorts of at-risk individuals and numerous clinical initiatives are collecting datasets consistent with the RDoC framework. In imaging genetics, publicly available cohorts show diverse psychological scores as well as imaging phenotypes such as regional cortical thickness or gyrification. Among the currently available integrative methods, some aim to adjust a classifier, a regressor, a clustering to find biomarkers (selection of combined features) or significant associations (univariate approach). In this paper, our focus will derive from the latter type of approach -integrative association study -, specifically aimed at uncovering associations between measurements coming from multiple modalities. Integrative association studies represent a departure from conventional univariate association analyses that consider a single measurement from one modality and find its associations with other modalities. Instead, we combine multivariate tools with association tools in accordance with the principles of precision psychiatry. These tools allow for modeling intra-and inter-modality structures, enabling the estimation of joint relationships between multiple measures across modalities. Integrative association methods face specific challenges, most notably the problem of missing data. This problem does not only affect data in a sporadic manner but can occasionally strikes all measurements of a modality for a patient. Furthermore, research in integrative association is also very active, especially regarding three axes: first, the integrative capacity of the considered multivariate models -whether they are linear or not -, second their interpretability, and third their generalisability or stability. These three points are important in assessing the contribution of these methods to precision psychiatry and are detailed below.

Integrative capacity. Pioneering works on integrative association models are Canonical Correlation Analysis (CCA) [Hotelling, 1936] or Partial Least Squares (PLS) [START_REF] Wold | Pls-regression: a basic tool of chemometrics[END_REF]. We retain these two methods as emblematic of integration analysis. Two essential concerns when using these approaches are the limited capacity of linear framework used, and the need to normalise the different data blocks beforehand. Beyond linearly adjustable covariates such as age or sex, the removal of other known confounding factors, such as imaging acquisition site, is more challenging. These factors introduce variability that may entangle the signal of interest in a manner that depends on the recruitment configuration. The need to address normalisation challenges has drawn attention to artificial neural network-based methods. For example, mutli-view AutoEncoders (AEs) have been shown to embed disentanglement capabilities [START_REF] Lee | Private-Shared Disentangled Multimodal VAE for Learning of Latent Representations[END_REF], which in turn could replace normalisation. In line with the idea proposed by [START_REF] Andrew | Deep canonical correlation analysis[END_REF], we propose to use Deep Learning (DL) for its greater integrative capacity and to leverage this disentanglement to avoid normalisation. Of particular interest, multimodal Variational AE (mVAE) is an extension of traditional VAE, as proposed by [START_REF] Kingma | Auto-encoding variational bayes[END_REF], designed to handle multiple views or modalities. It extends the idea of a VAE to capture the joint distribution of data from different views. The choice of the prior distribution in latent space represents the assumed distribution of the latent variables. Thus, it plays a crucial role in constraining the learning process and influences the structure of the latent space. Classical choices for the prior distribution in an mVAE include standard or multimodal Gaussian distributions [START_REF] Kingma | Auto-encoding variational bayes[END_REF], Suzuki et al., 2017, Wu and Goodman, 2018, Antelmi et al., 2019]. Recent proposals for prior distributions include Product of Experts (PoE), Mixture of Experts (MoE), or a combination of the two (MoPoE) [START_REF] Wu | Multimodal generative models for scalable weaklysupervised learning[END_REF], Shi et al., 2019, Sutter et al., 2021]. The rationale for adopting a mVAE lies in its ability to capture shared latent representations that account for intra-view correlation structure and model the inter-view correlation. However, questions remain regarding some limitations of mVAE in modeling only a shared latent space where information from different views is integrated [START_REF] Daunhawer | On the limitations of multimodal vaes[END_REF]. Recent research advocates exploring alternatives involving view-specific latent spaces [Sutter et al., 2021, Lee and Pavlovic, 2021, Daunhawer et al., 2021, Palumbo et al., 2022]. In this work, we propose to use a MoPoE-VAE as a specific mVAE, which promotes integrative capacity. This model can elaborate modality-specific and shared across modalities latent representations during the training phase. This shared representation conveniently dampens the influence of non-interest factors.

Interpretability. Using an mVAE to perform integrative association analysis comes at the expense of interpretability, a strength of their linear counterparts, such as CCA or PLS. Indeed, these classical methods have been praised for their interpretability. For example, adaptations of CCA and PLS propose to obtain interpretability by imposing latent covariation exploration with various sparsity constraints, including L1 or total variation regularisations [START_REF] Tenenhaus | Variable selection for generalized canonical correlation analysis[END_REF], Cao et al., 2011]. These implementations not only identify associations between different blocks, but also select the variables within each block that contribute to these correlations. This selection step facilitates the translation of these associations into interpretable knowledge. Bringing interpretation capabilities to mVAE, and more generally in DL, remains an opened and significant challenge. In previous work, we developed an interpretability method to explain the information learned by a mVAE. This method is used to evaluate the effect of controlled variations in the input of one view on the output of other views as generated by the mVAE [START_REF] Ambroise | Multi-view variational autoencoders allow for interpretability leveraging digital avatars: Application to the hbn cohort[END_REF]. In the present work, we adopt this novel interpretability framework, which is rooted in what we call a Digital Avatar (DA). This framework is used to study brain-behaviour relationships, allowing us to showcase interpretable associations between features across different views. Harnessing the generative capabilities of a trained mVAE, we generate sets of DAs from left-out subjects. Starting from a subject and varying its behavioural score, the model reconstructs the corresponding brain image for each score value. The analysis performed on the obtained sets of DAs yields what we hereafter call a Digital Avatar Analysis (DAA). It allows the exploration of brain-behaviour relationships through linear model fitting for all pairs of behavioural scores and regional image measures. The interpretability of DAs is enhanced by their exclusive reliance on the joint latent representation. However, DAA results are strongly affected by the epistemic variability [START_REF] Kendall | What uncertainties do we need in bayesian deep learning for computer vision?[END_REF] of the mVAE inherent to training stochasticity and weight initialisation. We mitigate the effects of epistemic variability through ensembling, using multiple trained mVAE models with different initialisation parameters on the same training subjects. This strategy, which acts as a regularisation of the DAA, is parameterised by the number of trained models and hereafter referred to as regularised-DAA (r-DAA). In summary, we expect that an ensemble of DAAs will provide a more comprehensive understanding of brain-behaviour associations.

Stability. The stability of results in integrative association studies within a multimodal dataset is a critical consideration to ensure the generalisability of the approach. Traditional approaches, such as sparse CCA, employ cross-validation strategies that often include the definition of a downstream auxiliary predictor for model selection. However, it has been shown that these procedures may not consistently identify stable associations [START_REF] Labus | Multivariate morphological brain signatures predict patients with chronic abdominal pain from healthy control subjects[END_REF], Baldassarre et al., 2017, Ing et al., 2019, Mihalik et al., 2020, Chegraoui et al., 2023]. This is particularly true when correlations between views are low or sample sizes are small [START_REF] Cao | Sparse PLS discriminant analysis: biologically relevant feature selection and graphical displays for multiclass problems[END_REF], Helmer et al., 2023, Yang et al., 2021, Nakua et al., 2023]. Clearly, the challenge of stability remains an open question, especially in the context of the DAA described above. To address this concern, we propose to use a sound framework. Stability selection is an inspiring technique developed in machine learning that efficiently exploits the regularisation mechanisms embedded in classical optimisation problems such as fitting a classifier [START_REF] Meinshausen | Stability selection[END_REF]. This procedure consistently identifies stable structures upon which the classifier is built, which inherently promotes robustness. In this study, we leverage Meinshausen's concepts regarding the stability of selected features and extend their application to the stability of selected associations. In fact, the DAA procedure, which generates the brain-behaviour associations, is not only affected by the epistemic variability, but is also sensitive to aleatoric variability (characteristics of the train / left-out split of the dataset). To reduce it, the stability selection procedure splits the dataset multiple times to generate different training and left-out splits. This allows the generation of stability paths, which identify stable brain-behaviour associations. The overall goal of this approach is to increase the robustness of our findings. This, in turn, contributes to a more reliable exploration of brain-behaviour relationships in the context of multimodal data analysis.

Overall, we present an empirical framework for extracting a consistent set of brain-behaviour associations from a multimodal dataset containing behavioural reports and brain imaging features. This framework is based on the stable selection of r-DAA output associations. In summary, our methological contributions are threefold. First, we propose a r-DAA that uses a weighted ensembling procedure to consolidate interpretations. Second, we introduce a stability selection procedure to generate a set of robust brain-behaviour associations. Finally, we show that our procedure, which is built on mVAEs, inherits from their ability to handle data with missing modality, and alleviates the need for additional normalisation. The proposed framework is applied to a publicly available cohort of at-risk children with notable behavioural symptoms. The cohort consists of structural MRI data and behavioural scores. Our results reveal interesting transdiagnostic brain-behaviour associations common to several psychiatric syndromes, such as Autism Spectrum Disorder (ASD) or Attention Deficit Hyperactivity Disorder (ADHD). These n E DAA outputs are then ensembled with an aggregation and a binary support decision function into n select associations. The r-DAA procedure is repeated for each of the N = 100 different splits of the complete dataset. (e) A stability selection procedure is used to retain associations with stability path Π ≥ π thr high enough. The whole framework outputs stable associations between clinical questionnaires and imaging.

Material and methods

This section introduces the proposed concept of performing stable association discovery in multimodal data using association study and stability selection. The digital avatar-based association study, referred to as DAA, is designed to find associations from a trained mVAE in heterogeneous multimodal datasets with potential missing data. Specifically, a regularised DAA (r-DAA) procedure is employed to mitigate the epistemic uncertainty inherent in training neural networks. Regularisation is achieved by ensembling the results of DAAs to retain meaningful associations. Finally, stability selection is applied to r-DAA by resampling the original dataset. This strategy is employed to handle the aleatoric uncertainty associated with the inherent variability in the data. Associations consistently selected by the r-DAA across these resampling splits are considered stable associations. The proposed approach consists of three nested steps listed from the innermost to the outermost: 1) DAA -training a mVAE and conducting an association study, 2) r-DAA -ensembling results from several DAAs derived on a same left-out set but with different initialisations and training batches, and 3) stability selection -performing several r-DAAs on different left-out sets to select stable associations. These steps are described below and illustrated in Fig. 1, and Supplementary Alg. D.1 and D.2.

Multimodal dataset

This study uses the Healthy Brain Network (HBN) cohort [Alexander et al., 2017]. This dataset is a multi-center, transdisciplinary clinical study. It includes a variety of assessments, including imaging, and a comprehensive set of psychological and clinical assessments to better understand psychiatric disorders. Inclusion criteria are not diagnosis dependent, but rather encompass an at-risk population with notable behavioural symptoms. Specifically, subjects were selected based on the presence of behavioural constructs related to Autism Spectrum Disorder (ASD) or Attention Deficit Hyperactivity Disorder (ADHD). Consensus diagnoses are not available for the majority of subjects enrolled in the HBN cohort. As such, the dataset allows for the study of the different manifestations of psychiatric syndromes within the data. In particular, it provides an opportunity to explore methods for new biomarker discovery and dimensional analysis.

In a previous study, our group identified seven behavioural assessments that capture the most salient dimensions in ASD patients [START_REF] Mihailov | Cortical signatures in behaviorally clustered autistic traits subgroups: a population-based study[END_REF]. Specifically, a quantitative measure of clinical autistic traits was defined as the parent Social Responsiveness Scale (SRS), hyperactivity levels were determined using the hyperactivity subscale within the Strengths and Difficulties Questionnaire (SDQ-ha), the anxiety was measured using the total score from the Screen for Child Anxiety Related Disorders Parent-Report (SCARED), irritability was defined using the total score of the Affective Reactivity Index Parent-Report (ARI), and finally, levels of depression, aggression, and attention problems were determined using subscales of the same names within the Child Behavior Checklist (CBCL-wd, CBCLab, and CBCL-ap, respectively). Filtering out subjects with sporadic missing data in questionnaires results in 2454 individuals.

The MRIs were acquired at four sites. A mobile 1.5T Siemens Avanto on Staten Island, a 3T Siemens Tim Trio at the Rutgers University Brain Imaging Center, and 3T Siemens Prisma at the CitiGroup Cornell Brain Imaging Center and at the CUNY Advanced Science Research Center. Collected T1-weighted images were preprocessed using FreeSurfer [Dale et al., 1999]. All results were manually reviewed in-house. The Euler number is used as a quality metric summarising the topological complexity of the reconstructed cortical surfaces [START_REF] Rosen | Quantitative assessment of structural image quality[END_REF]. Specifically, a single Euler number exclusion threshold of -217 is applied to yield 2042 selected subjects. Finally, cortical measures based on three metrics -cortical thickness, curvature, and area -are averaged in the 148 cortical regions of interest (ROIs) defined by Destrieux's parcellation [Destrieux et al., 2010].

The proposed brain-behaviour integrative analysis considers these two blocks of data. First, the electronic Clinical Record Forms (eCRF) view consists of the p eCRF = 7 behavioural scores. Second, the ROI view is composed of the p ROI = 444 cortical features from the 3 considered metrics across the 148 ROIs. In total, our dataset comprises 2991 subjects. Among them, 1505 have both complete views (referred to as complete dataset in Fig. 1(a) and 1486 have only one of the two views available (referred to as incomplete dataset in Fig. 1(a). Missing views remain a common problem in data integration. The factors contributing to missing data are usually not known in advance. Most traditional data mining and machine learning approaches operate on complete data and fail when data is missing. As fallback strategies, some models rely only on samples with all views available or on an auxiliary inference step that generates missing views. Our goal here is to use a model able to accommodate missing views, in order to use a maximum number of available samples.

Digital avatar analysis

Multimodal deep learning, which does not impose data normalisation and provides integrative capacitive models, is a promising approach to decipher the relationship between clinical dimensions, neuroimaging and a pathology. It has the potential to support the development of personalised medicine. If the considered model has generative capabilities, the exploration of a disorder spectrum can be achieved through the generation of so-called Digital Avatars (DAs). In the proposed work, we employ a mVAE and specifically use a MoPoE-VAE (see Supplementary A for details). We keep L = 301 subjects with complete data for the left-out set and use the remaining subjects, whether having complete or incomplete data, for the training set. The MoPoE-VAE can indeed be trained on data with missing views. Once trained, a MoPoE-VAE cannot directly provide brain-behaviour associations. However, the information learned by the model can be used to discover relevant pairwise associations between views, which we consider in our work to be relevant pairwise associations between each brain imaging feature and each score. Starting from the left-out subjects, we modify each subject's clinical score and generates a set of T = 200 DAs. Thus, as being predicted by the model, each DA bears brain imaging features at a given clinical state. From this set of virtual neuroimaging data, a conventional association study is performed to find brain-behaviour relationships. For each score of the eCRF view, we perform p ROI univariate association studies. In a previous publication, we introduced this DA analysis [START_REF] Ambroise | Multi-view variational autoencoders allow for interpretability leveraging digital avatars: Application to the hbn cohort[END_REF]. Hereafter, we describe in detail how the DAs are obtained and analysed.

For the DAA step, we use the left-out set for which all subjects have complete data. Without loss of generality and for simplicity, we consider the case where the eCRF view has only one score in what follows. Our aim is to create a set of DAs from each of the L left-out subjects. Let s = (s 1 , . . . , s L ) ∈ R L represent the scores from the eCRF view, and m = [m 1 , . . . , m L ] ∈ R L×p ROI represent all cortical measures from the ROI view. Note that this specific case can be generalised to all clinical scores and additional views. The general idea is to modify s into ŝ and observe how these modifications affect the reconstructed cortical measures m. In this way, virtual pairs (ŝ, m) are obtained. A first, although simplistic, approach to perturb s is to sample linearly or randomly between feature percentiles or even to bootstrap values across subjects. Conversely, the proposed approach relies on a simulation scheme to realistically perturb s. The realism of these perturbations stems from the generative aspect of the model. We hypothesise that such an approach, which models the composite variability of subjects measurements, will facilitate the discovery of interesting brain-behaviour relationships. Indeed, the variability in the simulated avatars integrates both subjectspecific and population-level variability. The latter source of variability is captured through the learned variance that is shared across the population. Our goal is to comprehensively assess what the model has learned by introducing subject-level perturbations and investigating their effects on the reconstructed digital avatars. The proposed strategy can be divided into three stages (inference, simulation, and association), which are outlined below and illustrated in Fig. 1(c) and Supplementary Fig. S2.

Inference -estimating likelihood distributions

For a given subject i ∈ {1, . . . , L}, we propose to use p θ (s i |z i ), the likelihood distribution of observations conditioned on the latent variable learned from the data, to sample DA score values. In p θ (s i |z i ), θ represents the weights of the eCRF view decoder, and

z i = (z i eCRF , z i Joint )
is the latent representation of s i , consisting of the eCRF specific and joint latent representations, respectively (see Supplementary A). This likelihood model will tend to sample likely values in the sense of the model. Therefore, provided our model is properly trained, such perturbed clinical scores will at best reflect the training data. The estimation of p θ (s i |z i ) is obtained by drawing D = 1000 realisations of z i ∼ q ϕ (z i |s i , m i ), where ϕ are the weights of the MoPoE-VAE encoders (see Supplementary A). Passing these latent representations to the eCRF decoder provides a good estimate of p θ (s i |z i ) [START_REF] Kingma | Auto-encoding variational bayes[END_REF] (by averaging the D decoded reconstructions). Note that the same strategy can be applied to categorical data by approximating the parameters of a Bernouilli distribution instead of a Gaussian.

Simulation -sampling realistic values

Given a subject i ∈ {1, . . . , L}, T samples are drawn from p θ (s i |z i ), resulting in T perturbations ŝi ∈ R T . Repeating this sampling for all subjects i, the generated perturbed set forms our DA perturbed eCRF observations ŝ = [ŝ 1 , . . . , ŝL ] ∈ R L×T . Importantly, when we generalise to multiple clinical scores, only one score is perturbed at a time. Finally, the perturbed ROI measures are reconstructed using a forward pass by considering ŝ and the corresponding ROI features m as input to the model. This results in a set of T perturbed ROI measures representing our DAs m ∈ R L×T ×p ROI . As a compromise between computational cost and accuracy, we set T = 200.

Association -computing inter-view associations

We search for associations between a clinical score and each ROI measure. We refer to the generated DA score values as ŝ ∈ R L×T , and generated DA ROI measures as m ∈ R L×T ×p ROI . Associations are obtained using hierarchical regression models [START_REF] Bryk | Hierarchical Linear Models: Applications and Data Analysis Methods[END_REF]. Specifically, for each image feature k ∈ {1, . . . , p ROI }, a linear regression of the form ŝi

= c i k + b i k mi k + ϵ i k is fitted for each subject i.
The resulting slopes are averaged over all subjects as

β k = 1 L L i=1 b i k .
Performing the same analysis for all available eCRF scores results in the association vector β ∈ R p , where p = p ROI * p eCRF . This vector encompasses all potential associations between the p eCRF clinical scores of the eCRF view and the p ROI cortical features of the ROI view.

Regularised digital avatar analysis

In our seminal paper [START_REF] Ambroise | Multi-view variational autoencoders allow for interpretability leveraging digital avatars: Application to the hbn cohort[END_REF], we observed that the associations were not stable, even when the model was retrained on the same training set, and evaluated on the same left-out set. This instability is likely due to epistemic variability [Kiureghian andDitlevsen, 2009, Kendall andGal, 2017], and has mainly been studied in supervised settings [Gal andGhahramani, 2016, Lakshminarayanan et al., 2017]. In feature selection using deep learning, studies have shown that ensembling can enhance stability and thus mitigate epistemic variability [START_REF] Gyawali | Ensembling improves stability and power of feature selection for deep learning models[END_REF]. We propose to repeat the previous step n E times: we train a MoPoE-VAE on the training subjects, then perform the DAA on the obtained model using the left-out set. Note that the training set and left-out set remain the same throughout these n E procedures, only the random weights initialisations and batches during training vary. This is illustrated in Fig. 1. By subsequently using ensembling, our objective is to identify stable associations from n E DAAs association matrix β = [β 1 , . . . , β n E ] ∈ R n E ×p . Like classical deep ensembling, our approach involves ensembling candidate associations (i.e., model predictions in supervised settings [START_REF] Lakshminarayanan | Simple and scalable predictive uncertainty estimation using deep ensembles[END_REF]) proposed by the n E DAAs, and regrouped in the β matrix. The proposed ensembling procedure requires the definition of an aggregation function f and a decision function g. The aggregation function summarises the n E associations coefficients. The decision function generates a binary decision support from the aggregated coefficients, which is designed to retain meaningful associations. Formally, let f : R n E ×p → R p be an aggregation function and g : R p → {0, 1} p be a decision function. The composition g • f forms the proposed ensembling, which is outlined below and illustrated in Fig. 1(d) and Supplementary Alg. D.1(b) and D.2(b).

Aggregationf

The role of the aggregation function f is to assign an importance to each association from the association matrix β. Initially, we opt for a classical mean, although alternative functions such as median or maximum could also be relevant. However, the quality of the estimated n E models is not accounted for by any of these functions. In fact, some models may converge to local minima of the optimised loss function, resulting in less relevant models and representations. Such behaviour can produce outlier associations that affect the stability of selected associations [START_REF] Gyawali | Ensembling improves stability and power of feature selection for deep learning models[END_REF]. Therefore, we employ a weighted average where the weights reflect the quality of the estimated models. The weights are determined by the model's joint latent space ability to capture a significant proportion of the eCRF-related variability. Practically, we use the representational similarity analysis (RSA) [Kriegeskorte et al., 2008], which generates correlations between the joint latent space and the eCRF scores (see Supplementary B and Supplementary Alg. D.2 for details).

Decisiong

The purpose of the decision function g is to obtain binary decision support from the aggregated coefficients on a scoremetric basis. This support allows us to choose a subset of the most informative associations for each score-metric pair (s, m), enhancing reproducibility and interpretability. While a simple strategy involves specifying a numeric threshold or defining built-in heuristics to find this threshold, we prefer setting the number of features explicitly rather than using a threshold. This approach provides consistency across score-metric pairs and DAAs. In practice, we choose to select the top n select = 12 associations with the greatest amplitudes.

Stability selection

Stability selection, derived from penalised machine learning [START_REF] Meinshausen | Stability selection[END_REF], is designed to reinforce feature selection through the definition of stability paths. Stability paths represent the probability of selecting each feature when training the same algorithm with some regularisation parameter on different random splits of a dataset. We extend this methodology to our brain-behaviour association study. In our context, the regularisation parameter is the number of models n E considered in the ensembling step. We propose to repeat the r-DAA (see Section 2.3) N times, using different splits of the original dataset, and setting the regularisation parameter n E ∈ {1, . . . , 20}. This approach allows us to model aleatoric uncertainty inherent to population variability [START_REF] Kendall | What uncertainties do we need in bayesian deep learning for computer vision?[END_REF]. The implementation comprises three aspects outlined below and illustrated in Fig. 1(e) and Supplementary Alg. D.1(c) and D.2(c). First, we define a valid splitting strategy. Then, we repeat the r-DAA N times by varying the regularisation parameters n E , enabling us to estimate stability paths. Finally, we define a criterion to assess the stability of the results obtained across the N splits.

Splitting strategy

Our dataset consists of two parts: 1505 subjects with complete data and 1486 subjects with incomplete data (see Section 2.1). It's important to note that subjects with missing views can only be used to train the model. For each split, out of the 2991 subjects available, 2690 are used in the training set. The latter includes the 1486 subjects with incomplete data and 1204 randomly selected subjects with complete data. The remaining L = 301 subjects with complete data form the left-out set (see Fig. 1). To maintain population statistics, including age, sex, and acquisition site distributions, we employ shuffled iterative stratification [Sechidis et al., 2011]. Although using the entire incomplete dataset in each training iteration may appear as a limitation, it is effectively mitigated by employing different shuffled batches at each training epoch. We opt for using N = 100 splits.

Stability paths

Each r-DAA with regularisation parameter n E produces a decision support S j (n E ) ∈ {0, 1} p for a given split j ∈ {1, . . . , N }. The stability paths Π n E ∈ [0, 1] p are the probability for each association to be selected across the N splits. It is obtained by averaging all the calculated binary decision supports as:

Π n E = 1 N N j=1 S j (n E ) (1)

Stability criterion

The stability criterion analyses the stability paths and defines which associations are considered stable. ability of an association happens to be greater than a user-defined threshold π thr ∈ [0, 1] for a specific regularisation parameter n E , then the corresponding association is selected. Finally, we define the set of stable features as:

S stable = {k : max n E ∈{1,...,20} Π n E k ≥ π thr } (2)
More simply, this means that associations with a high probability of selection for any regularisation parameter are retained as stable associations. Conversely, those with a low probability of selection are dropped.

Results

Models inspection

To evaluate the information encoded in the various latent spaces of the trained models, we employ the RSA [Kriegeskorte et al., 2008] tool introduced in Supplementary B. The RSA results in Table 1 depict the relationships between different latent spaces (eCRF, Joint, and ROI) and each eCRF score, alongside other relevant covariates like age, sex, and image acquisition site. In short, we computed each of these Kendall τ correlations for every N = 100 splits and corresponding n E = 20 models. The reported correlations τ are averaged across these N * n E values. See Supplementary B for further details.

Notably, each eCRF score strongly correlates with the eCRF-specific latent space, as highlighted in pink in Table 1. This is noteworthy since the models successfully learned in a single one-dimensional space, capturing a significant amount of variability associated with all eCRF view scores. This outcome is expected, considering the correlations between the eCRF scores (see Supplementary Fig. S1). Additionally, it underscores that this latent space is not informing about age, sex, or image acquisition site. In contrast, the latent space specific to the ROI view does not significantly correlate with the eCRF scores with the only exception of ARI, yet with a small correlation. Note that it significantly correlates with age, sex, and image acquisition site, as shown in the cells highlighted in orange. Finally, we highlight in salmon in Table 1 what the joint latent space has learned. These representations moderately and significantly correlate with each score but not with acquisition site or sex. Moreover, it appears to correlate with age. It seems that the brain-behaviour relationships, modeled in the MoPoE-VAE joint latent space, relies on information related to age.

Stability selection from r-DAA

By increasing the number of models and applying a stability selection procedure, we aim to increase the stability of the sets of brain-behaviour associations supported by our multimodal dataset. In the following, we consider the SRS score associations with the thickness metric. Figure 2 illustrates the associations trajectories taken by the considered 148 ROIs. This illustration is inspired by the figure of merit proposed in the seminal paper [START_REF] Meinshausen | Stability selection[END_REF]. In Figure 2(a), we display trajectories of aggregated β values (i.e., effect sizes), using a mean function for f , against the number of model n E . Some trajectories become more prominent than others as the number of models increases. However, they are hardly distinguishable from other trajectories, which remain densely packed with low values. This general aspect is displayed here for one of the N = 100 stability selection splits. This suggests that we can isolate some, but not all, stable associations by increasing the number of models aggregated in the r-DAA.

Figure 2(b) presents the stability paths Π for each ROI with respect to the number of models n E . This probability is computed as described in Eq. 1 and summarises the N = 100 stability selection splits. Here, the aggregation function f is the mean function. The results highlights that the most stable ROIs, characterised by higher probabilities Π, are easily identifiable. The stable paths reach a plateau with the number of models (i.e., the considered regularisation hyperparameter) between 5 and 10 models, depending on the considered ROI.

Figure 2(c) is similar to Figure 2(b), but the procedure uses the RSA-weighted average as the aggregation function f instead of a simple average (see Supplementary Alg. D.2). The results suggest that stability is only slightly affected by the choice of the aggregation function f . We keep the RSA-weighted average because it may be more robust to outlier models. We apply a threshold of π thr = 0.4 (as illustrated by the dashed red line) to retain stable associations, displayed in colors. These colors correspond to the same ROIs throughout the different plots. Comparing Figure 2(b) or Figure 2(c) with Figure 2(a) reveals that a high coefficient amplitude is not always a good indicator of stability. In fact, the ROIs represented by pink or light purple paths in Figure 2(a) are indistinguishable from other black dotted lines (i.e., unselected ROIs). However, these ROIs convey some of the most stable associations (Π ≃ 0.7 when n E = 20). Note that the variability of the SRS score is not the best captured by our models on our dataset (see Table 1). To test the influence of the selected score and metric on the stability paths, we generate similar figures by examining the association coefficients between the SDQ-ha score and the area metric, as well as the stability paths. Very similar observations can be made, as shown in Supplementary E.

Transdiagnostic-factor spatial support

At the core of our approach is a MoPoE-VAE network designed to construct a latent space consolidating joint information between the entire ROI view and the entire eCRF view, along with ROI-and eCRF-specific latent spaces. The interpretation of joint latent representations amounts to extracting brain-behaviour associations that effectively represent this joint information. This strategy supports the research on the general psychopathology dimension introduced by Caspi [START_REF] Caspi | The p factor: One general psychopathology factor in the structure of psychiatric disorders?[END_REF], Caspi et al., 2023]. In line with these ideas, our goal is to identify the brain regions and metrics that underpin these similarities. These brain-behaviour associations retained by our approach for each score-metric pair (s, m) are listed in see Supplementary F.

In this section, we focus on transdiagnostic associations whose ROIs exhibit co-association with exactly four specific scores. These regions could be considered as components of a transdiagnostic spatial support for mental disorders. We focus on four specific scores of interest: SRS, which assesses social interaction and is correlated with autism disorder; SCARED, which assesses fear and anxiety issues and is correlated with anxiety disorders; SDQ-ha, which is related to hyperactivity disorders; and CBCL-wd, an indicator of depression, a symptom common to most of these pathologies. For each cortical measure, Figure 3 highlights the ROIs found to be associated with these four scores. Among these selected transdiagnostic regions, many belong to the pericalosal and cingulate regions. These regions, whether considered with the area or curvature metrics, are consistently associated with each examined score. This means that the associations between selected ROIs and scores show similar covariations. In the identified pericalosal and cingulate regions, a decrease in both area and curvature metrics is associated with an increase in the SRS, SCARED, or CBCL-wd scores. However, the opposite is observed for the SDQ-ha score (Figure 3(b)-(c)). Looking at the thickness metric and grouping the curvature and area metrics, we find two disjoint sets of regions among the selected transdiagnostic regions. First, the left and right occipital poles for the thickness metric. Second, the cingulate regions for the curvature/area metrics. Only the right cingulate mid / posterior region (R.Ci.Mid.Post.) and left cingulate posterior dorsal gyrus (L.G.Ci.Post.Dors) seem to be related to the four scores for both metrics. Finally, looking at the SDQ-ha score and comparing to the grouped SRS, SCARED and CBCL-wd scores, we find that the associations identified have systematically opposite signs. Overall, the transdiagnostic regions are mostly bilateral and present a spatially smooth pattern of association with the eCRF scores.

Transdiagnostic-factor spatial support in ASD and ADHD

Recall that our method, with the presented setting, is designed to find transdiagnostic associations, in line with the general psychopathology dimension. Among the structural biomarkers involved in ASD (or ADHD), we can study the ones that can attributed to transdiagnostic factors. In this section, we focus on transdiagnostic associations whose ROIs are specifically associated with one symptom/score. We perform this analysis for two scores. First, the SRS score, which is considered as a proxy for ASD [START_REF] Constantino | Validation of a brief quantitative measure of autistic traits: comparison of the social responsiveness scale with the autism diagnostic interview-revised[END_REF]. Although the SRS score alone may not be sufficient, clinicians often use it to determine a patient's status. This score has also a high correlation with diagnosis (unpublished results based on our databases) in datasets adhering to diagnosis-balanced inclusion criteria such as ABIDE I (SRS-1, 0.82), ABIDE II (SRS-1, 0.86 and SRS-2, 0.72), or EU-AIMS (SRS-2, 0.85). With the same general caveats, we use the SDQ-ha score as a proxy for ADHD-related symptoms and diagnosis [Goodman, 1997].

The retained associations can be found in Supplementary F and are summarised below.

Markers related to SRS (ASD).

In Figure 4, we display ROIs associated in thickness and area with the score SRS. Interestingly, the selected ROIs display a rather symmetrical pattern. With the thickness metric, the SRS appears to positively covary with the bilateral occipital poles, bilateral pericallosal region, and a right prefrontal region, while negatively covarying with the bilateral temporal sulci and right postcentral region. With area metrics, the SRS appears to positively covary with the left prefrontal cortex, left subcallosal gyrus and bilateral circular inferior insular sulci, while negatively covarying with the bilateral cingulate and left parieto-temporal cortex.

Discussion

In this paper, we present an interpretability method dedicated to deep learning-based integration models using mVAEs.

We demonstrate its application in the study of the transdiagnostic dimension within the HBN cohort, integrating neuroimaging data and psychological assessments in an at-risk population with notable behavioural symptoms. Our method is endowed with a stability selection procedure to retain associations and we investigate its effects. A sufficient number of models is the only need to achieve stability. Finally, with only a prior on the expected number of associations, our method enables the identification of stable associations between measures of the brain cortical surface and symptom scores.

Deep mVAE model are expressive multi-modal integration tools

We take advantage of the versatile definition of the latent space in the MoPoE-VAE to choose a representation setting that naturally disentangles specific and shared sources of variation between view-specific and joint latent spaces.

Other work already identified the disentanglement achieved with this type of architecture [START_REF] Lee | [END_REF]Pavlovic, 2021, Qiu et al., 2022], but we apply it here to neuroimaging and score questionnaires integration. We also leverage this source separation using our interpretation module based on digital avatars. By mitigating or eliminating confounding effects such as the MRI acquisition site in shared representation, we expect interpretations to be unaffected, without performing any standardisation or harmonisation of the data prior to learning. Moreover, the joint latent spaces that significantly correlate with age highlight how our models handle age information, although it is never explicitly provided. This contrasts with conventional approaches using age residualisation.

Our framework based on the MoPoE-VAE can conveniently handle incomplete data, a common issue in multi-view integration. The requirement of complete data often significantly reduces the number of subjects that can be used and impairs the statistical power of most studies. This resilience made it possible for us to use an openly available cohort with minimal missing data control and include nearly all available subjects. This has provided us with a substantial sample size collected from multiple acquisition sites. This characteristic should improve reproducibility and replicability when employing such an approach in multi-modal neuroimaging studies, aligning with state-of-theart guidelines [START_REF] Klapwijk | Opportunities for increased reproducibility and replicability of developmental neuroimaging[END_REF].

We also show that MoPoE-VAE networks, when equipped with a stable interpretation module, can provide associations between variables of different views. This stable interpretation module can be configured with only one a priori parameter which is the expected number of associations (n select = 12 here).

Deep mVAE models for transdiagnostic studies

Transdiagnostic researches [START_REF] Fusar-Poli | Transdiagnostic psychiatry: a systematic review[END_REF] focus a common risk factor of pathology in psychiatry, often referred to as p-factor [START_REF] Caspi | The p factor: One general psychopathology factor in the structure of psychiatric disorders?[END_REF], computed using various clinical assessments and summarised in a unique dimension. These approaches address the fact that many psychiatric illnesses share common symptoms and comorbidities, and that the heterogeneity of the population cohorts makes it difficult to obtain reproducible results based on a single diagnosis. Studies hypothesise a common etiology, that would explain some underlying mechanisms common to multiple psychiatric disorders. They usually express this global pathological variability as p factor, and then inspect its association with brain imaging to find potential biomarkers for these common mechanisms.

Deep learning MoPoE-VAE models are trained from data in an at-risk population cohort in which subjects are assessed with questionnaires, expressing symptoms of several psychiatric syndromes (SRS, SDQ-ha, etc). This MoPoE-VAE training is organised to integrate all available data to create distinct representations of shared and specific information. This shared representation contains multivariate variability linking multiple symptom scores to imaging. Our stable interpretation pipeline inspects this joint latent space to produce transdiagnostic associations between behaviour and cortical measures. By displaying regions linked with multiple symptom scores, we capture regions often reported as associated with transdiagnostic factors in the literature.

In particular, results from our transdiagnostic perspective in section 3.3 highlight fairly symetrical regions, mostly associated in area and curvature in cingulate regions and in thickness in the occipital poles. Studies in functional MRI already identified cingulo-opercular network (CON) [START_REF] Sheffield | Transdiagnostic Associations Between Functional Brain Network Integrity and Cognition[END_REF] and default mode network (DMN), and particularly the anterior cingulate areas [START_REF] Gong | Network-level dysconnectivity in drug-naïve first-episode psychosis: Dissociating transdiagnostic and diagnosis-specific alterations[END_REF], MacNamara et al., 2017] as related to transdiagnostic factors. Other fMRI studies identified anterior / middle cingulate and occipital areas, among other, associated with general p-factor [START_REF] Elliott | A connectome-wide functional signature of transdiagnostic risk for mental illness[END_REF], as well as in studies of transdiagnostic population with other factors [START_REF] Feldker | Transdiagnostic brain responses to disorder-related threat across four psychiatric disorders[END_REF], Tong et al., 2022], and in particular with ASD and ADHD [START_REF] Bush | Anterior cingulate cortex dysfunction in attention-deficit/hyperactivity disorder revealed by fmri and the counting stroop[END_REF], D'Cruz et al., 2016, Lukito et al., 2020].

Alongside fMRI studies, many studies report structural alterations associations with transdiagnostic factors, such as changes in cortical thickness or grey matter volume of the cingulate and occipital cortices, as part of more global patterns [START_REF] Clementz | Identification of distinct psychosis biotypes using brain-based biomarkers[END_REF], Yin et al., 2022, Parkes et al., 2021]. Such findings were highlighted as well in specific ASD [START_REF] Oblak | Decreased gabab receptors in the cingulate cortex and fusiform gyrus in autism[END_REF], Chien et al., 2021, Mihailov et al., 2020, Ecker et al., 2022] or ADHD [START_REF] Amico | Anterior cingulate cortex gray matter abnormalities in adults with attention deficit hyperactivity disorder: A voxel-based morphometry study[END_REF], He et al., 2015, Bayard et al., 2020] studies, as well as research studying these diseases jointly [START_REF] Rommelse | Structural brain imaging correlates of asd and adhd across the lifespan: a hypothesis-generating review on developmental asd-adhd subtypes[END_REF], Lukito et al., 2020]. Moreover, occipital lobe grey matter reduction was specifically identified as linked with increases in p-factor [START_REF] Romer | Structural alterations within cerebellar circuitry are associated with general liability for common mental disorders[END_REF], Romer et al., 2021]. Finally, cingulate related white matter tracks fractional anisotropy extracted from diffusion MRI were found associated with transdiagnostic factors as well [START_REF] Stefanik | Brain-behavior participant similarity networks among youth and emerging adults with schizophrenia spectrum, autism spectrum, or bipolar disorder and matched controls[END_REF].

Cingulate regions are well-known to be implicated in cognition, emotion processing [START_REF] Bush | Cognitive and emotional influences in anterior cingulate cortex[END_REF], Vogt et al., 1992, Denson et al., 2009], while occipital regions are responsible for primary visual processing. Since these functions are often altered when expressing psychiatric symptoms, their implication in a common mechanism PREPRINT to multiple syndromes is very likely. Our results support the hypothesis that cingulate and occipital regions and their related functional or structural networks are important in general psychopathology and warrant further investigation.

Perspectives

The method described in this paper focuses on an integrative multi-view approach based on deep learning and featuring interpretation capacity. We discussed the main benefits of its use and we list below some perspectives related to some limitations of the proposed approach, while highlighting more general remarks regarding the integration of the ever growing volume and variety of available data in neuroscience.

Effect size. In our transdiagnostic study, the stable associations found are characterised by rather small effect sizes. This underscores the difficulty of exploring the common processes that hypothetically contribute to the etiology of multiple psychopathologies. While individual effect sizes may be modest, the cumulative evidence supports the existence of transdiagnostic factors. One could expect an increase of the observed effect sizes by considering functional MRI (fMRI) data instead of structural MRI. Association studies using fMRI offer numerous advantages, including direct measurement of brain activity and localisation of function. Current research with fMRI shows strong associations compared to structural MRI [Oblong et al., 2023] and would be worth investigating using our tool.

Data harmonisation. Several recent works have shown that machine learning models are strongly biased by the MRI acquisition site and do not generalise well to new MRI images from sites that have never been seen before [START_REF] Glocker | Machine learning with multi-site imaging data: An empirical study on the impact of scanner effects[END_REF], Wachinger et al., 2021]. This problem is due to differences in scanner manufacturers, specifications, settings, and hardware. While traditional residualisation technique applied to remove the site effect marginally improves the performance of machine learning models, it does not bring any improvement for deep learning models [START_REF] Dufumier | Exploring the potential of representation and transfer learning for anatomical neuroimaging: application to psychiatry[END_REF]. In this work, we show that the influence of non-interest factors, in particular the site effect, can be effectively eliminated by disentangling modality-specific and joint latent representations. Nevertheless data harmonisation remains an ongoing research area [START_REF] Dinsdale | Deep learning-based unlearning of dataset bias for mri harmonisation and confound removal[END_REF], Bashyam et al., 2022]. For example, the OpenBHB challenge [START_REF] Dufumier | Openbhb: a large-scale multi-site brain mri data-set for age prediction and debiasing[END_REF] on brain age prediction with site effect removal could bring new harmonisation techniques developed by the community.

Views with different internal correlation structure. As mentioned above, a wide variety of assessments are available in large population cohorts. Most often, each type of data is exploited separately for regression, classification and segmentation. Deep learning has been extensively used on imaging data and has demonstrated considerable benefits over traditional methods for some tasks, such as the segmentation of medical anomalies and anatomical structures [START_REF] Menze | The multimodal brain tumor image segmentation benchmark (brats)[END_REF], Henschel et al., 2020]. But in the medical field, many assessments, such as imaging or genomic data, can be more naturally represented as graphs due to their underlying biological properties (e.g. functional / structural connectomes or protein interactions). DL specific operators [START_REF] Bronstein | Geometric deep learning: Going beyond euclidean data[END_REF], Ghosal et al., 2022] have been developed to take into account such correlation structure when learning from such data. In the current study, we integrate cortical measures on ROIs as tabular data. It would be highly relevant to return to brain surface data and model their actual correlation structure using dedicated convolution operators [START_REF] Zhao | Spherical u-net on cortical surfaces: Methods and applications[END_REF]. The same applies to other available modalities such as genomic data.

Dimensional discrepancy of the views. The demonstrated integrative capacity of mVAEs comes from modelling modality-specific and joint latent spaces. Handling views with very different number of variables is an ongoing research question. For example, genotyping data may have millions of variables. Training a mVAE in such setting is challenging. Intuitively, learned modality specific latent spaces will have different sizes, somehow proportional to their input size, but this is not so clear for the joint representations. Pioneering work defines a sparse multichannel VAE [Antelmi et al., 2019]. It leverages a variational dropout regularisation [START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF] that identifies an optimal number of joint latent dimensions.

Disentanglement. In the development of digital avatars, transparency of the trained NN models is a key feature. Importantly, the chosen mVAE model must learn disentangled representations that separate modality specific and joint variability. Indeed, confounding factor like site effect are usually specific to one modality, and such disentanglement would ensure the joint latent is devoid of such unwanted effects. Recent works propose contrastive VAE which uses deep encoders to capture higher-level semantics [START_REF] Aglinskas | Contrastive machine learning reveals the structure of neuroanatomical variation within autism[END_REF], as compared to its former linear conterpart contrastive PCA [START_REF] Abid | Exploring patterns enriched in a dataset with contrastive principal component analysis[END_REF]. From two encoders, they typically structure the learned latent space into two parts containing the background (i.e., common to the studied population) and salient (i.e., specific to a pathology) variabilities. Followup works show that classical VAE losses alone can not effectively separate joint and salient variability, and that further constraints and regularisations are needed to satisfy the assumptions of the generative process and to promote disentanglement in the latent space [START_REF] Abid | Contrastive variational autoencoder enhances salient features[END_REF], Choudhuri et al., 2020, Weinberger et al., 2022, Zou et al., 2022]. In the case of mVAEs, the same observations hold. Modeling modality-specific and joint latent spaces somehow separates joint from specific variability [START_REF] Lee | Private-Shared Disentangled Multimodal VAE for Learning of Latent Representations[END_REF], as observed in this study. However, to ensure disentanglement, additional regularisations during training are needed [START_REF] Daunhawer | Self-supervised Disentanglement of Modality-Specific and Shared Factors Improves Multimodal Generative Models[END_REF].

Replication. The design of a replication study will contribute to the generalisability of our findings. In particular, it would be interesting to apply our approach to cohorts such as the Dunedin Longitudinal Study [START_REF] Caspi | The p factor: One general psychopathology factor in the structure of psychiatric disorders?[END_REF], Romer et al., 2021], the ABCD Study [Casey et al., 2018, Karcher andBarch, 2021] or the Duke Neurogenic Study [START_REF] Romer | Structural alterations within cerebellar circuitry are associated with general liability for common mental disorders[END_REF], Elliott et al., 2018]. These studies were not explicitly focused on transdiagnostic research, but their design, including comprehensive assessments with a multidisciplinary context, has provided valuable insights into the understanding of psychiatric disorders from a transdiagnostic perspective.

A Multimodal Variational AutoEncoders

This appendix presents details on the multiview variational autoencoder (mVAE) models, used to perform the Digital Avatar Analysis, once they are trained.

Multimodal (or multiview) VAEs are extensions of the well-established VAE [START_REF] Kingma | Auto-encoding variational bayes[END_REF]] that allow joint integration and reconstruction of two or more views. These two VAE models are themselves derived from AEs which attempt to find a compressed representation that contains as much information as possible about the input data. They are based on the assumption that the input data can be generated from a lower dimensional space [Kramer, 1991]. AEs have been extensively studied [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF] and used in various applications such as medical image segmentation [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. In the variational version of anAE, the main idea is to learn a probability distribution over the low dimensional latent space. Finally, multimodal VAEs allow the construction of view-specific and shared latent spaces via the different modeling schemes. This is of great interest in medicine, for example to separate common sources of variability from disease-specific one [START_REF] Wu | Multimodal generative models for scalable weaklysupervised learning[END_REF], Shi et al., 2019, Sutter et al., 2021].

The selected model Let x = {x 1 , . . . , x K } be an observation of K views, where each x k can have different dimensions J k . The general form of the latent of the variable is {z, z 1 , . . . , z k , . . . , z K } with z denoting the d-dimensional latent variable shared by all x k , and z k the view-specific d k -dimensional latent variables. We assume the following generative process for the observation set:

z ∼ p(z) z k ∼ p(z k ) x k ∼ p θ (x k |z, z k )
where p(z), p(z k ) are prior distributions for the latent variables and p θ (x k |z, z k ) is a likelihood distribution of the observations conditioned on the latent variable. The generative process p θ (x) is approximated by its Evidence Lower Bound (ELBO) . The bottleneck of this approximation problem comes to computing the posterior p θ (z, z k |x). A first usual assumption about p θ (z|x) considers that it contains all the shared information between the views, and we can write p θ (z, z k |x) = p θ (z|x) × p θ (z k |x k ). Other assumptions are needed since these posteriors are not analytically tractable. The posterior distributions are classically approximated by parametrized distributions q ϕ (z|x) (the so-called joint variational posterior) and q ϕ (z k |x k ) respectively.

State-of-the-art models solve this problem with different ways of modeling the joint variational posterior distribution q ϕ (z|x). These approaches compose q ϕ (z|x) using the Gaussian experts q ϕ (z|x k ) unimodal posteriors, which may be combined in various ways. In Product of Experts (PoE) approach [START_REF] Wu | Multimodal generative models for scalable weaklysupervised learning[END_REF], q ϕ (z|x) is defined as a product of individual Gaussian experts, assuming the conditional independence of the experts, hence remains Gaussian. The product operator used in PoE may dampen the specific variability of a given block described by gaussian unimodal posterior, and other work proposes to model [Shi et al., 2019] q ϕ (z|x) using a Mixture of Experts (MoE), assuming each unimodal experts contributes equally to the joint posterior. MoPoE approaches defines q ϕ (z|x) as a Mixture of Product of Experts. Each product of experts is computed over the different subsets of views [Sutter et al., 2021]. We use the latter strategy, referred to as MoPoE, as it is a generalization of the two previous approaches. Each strategy has different ways of handling missing views, that may differ between train and test time. MoPoE simply uses the subsets of available views for each observation. As regards to the unimodal posteriors q ϕ (z k |x k ), we classicaly use Gaussian distributions.

The MoPoE optimizes a generalized multimodal ELBO objective for learning view-specific and a joint distribution of multiple views x with potential missing data, defined as:

C Digital Avatar Analysis

D Algorithm pseudocodes

We display below the pseudocode of the two procedures developed for our association discovery pipeline, relying on r-DAA and stability selection. x j train , x j left-out ← X[j] ▷ a random split 7:

for n E ∈ {1, . . . , 20} do 8:

Fit a M oP oE j n E model on x j train using init random weights θ n E 9:

Compute the DAA with x j left-out association weights β j n E ∈ R p Compute aggregation scores f ((β j 1 , . . . , β j n E )) ∈ R p 16:

A j n E ← f ((β j 1 , . . . , β j n E ))

17:

Compute decision support g(A j n E ) ∈ {0, 1} p 18:

S j (n E ) ← g(A j n E )

19:

end for 20:

Compute the selection probability of each feature Π n E ∈ R p 21:

Π x j train , x j left-out ← X[j] ▷ a random split 7:

for n E ∈ {1, . . . , 20} do 8:

Fit a M oP oE j n E model on x j train using init random weights θ n E 9:

Compute the DAA with x j left-out association weights β j n E ∈ R p 10:

Compute the left-out subjects pairwise joint latent representation dissimilarity matrix C j n E ∈ R 301×301 Compute aggregation scores f ((β j 1 , . . . , β j n E ), (τ j 1 , . . . , τ j n E )) ∈ R p 25:

A j n E ← f ((β j 1 , . . . , β j n E ), (τ j 1 , . . . , τ j n E )) =

Figure 1 :

 1 Figure 1: Illustration of the interpretation framework based on DAA. In (a) the HBN cohort with two modalities: imaging data (metric × ROI) and clinical questionnaires (eCRF). For each split (orange bounding box) of the dataset, complete data are split into training (plain orange removing the dotted part) and left-out (dotted orange part) subjects. Each r-DAA (blue bounding box) procedure uses one split of the dataset. In the inner dotted blue bounding box, we repeat the following n E times : (b) a MoPoE-VAE is trained using all the subjects with one modality and the training subjects with both. (c) Left-out subjects are used for DAA using this trained model, outputting coefficients β. (d)These n E DAA outputs are then ensembled with an aggregation and a binary support decision function into n select associations. The r-DAA procedure is repeated for each of the N = 100 different splits of the complete dataset. (e) A stability selection procedure is used to retain associations with stability path Π ≥ π thr high enough. The whole framework outputs stable associations between clinical questionnaires and imaging.

  Figure2: Investigation of the ROIs associated with SRS in thickness. Each line corresponds to a ROI, dotted black ones are not selected as below the threshold π thr = 0.4 when using as aggregation function a mean weighted by RSA correlations. The ROIs thus selected are colored and the colors are consistent across the plots. The red horizontal dotted line highlights the threshold π thr = 0.4. (a) Mean coefficients aggregated across models for a given split over the N = 100 splits, plotted against the number of models. (b) ROIs stability paths Π when using an uniformly weighted mean to aggregate the coefficients output from DAA on each model in the r-DAA, against the number of models used. (c) Same as the last one, except that the mean was weighted using RSA correlations. The latter strategy is used to select the ROIs associated in thickness with the SRS score, with Π > π thr = 0.4.

Figure 3 :

 3 Figure 3: Display of ROIs associated with selected eCRF scores for a transdiagnostic assessment, e.g. SRS, SDQ-ha, SCARED and CBCL-wd, for each cortical measure. (a) ROIs color-coded, corresponding to color in the following polar plots. Lateral and medial view are displayed for each hemisphere (left hemisphere on the left hand side and right hemisphere on the right). The polar plots indicate the sign of each association (negative when the bar is close to the center, postive when it is close to the edge of the circle). Only retained associations are displayed. Each polar plot correspond to a metric : (b) displays assocations in thickness, (c) associations in curvature and (d) in area. L: left, R: right, Ci: cingulate, Oc: occipital, Te: temporal, S: sulcus, G: gyrus, Mid: middle, Post: posterior, Dors: dorsal, med: medial, Inf: inferior, Marg: marginal, Ant: anterior, PeriCal: pericallosal, PostCe: postcentral, Sup: superior, Li: lingual, Pl: plan, Pola: polar.

Figure 4 :

 4 Figure 4: Display of ROIs associated with the SRS score.

Figure 5 :

 5 Figure 5: Display of ROIs associated with the SDQ-ha hyperactivity score.

LFigure S1 :

 S1 Figure S1: Correlations between the 7 eCRf scores.

Figure S2 :

 S2 Figure S2: Illustration of the interpretation framework based on DAA in a clinical cohort setting with two modalities:imaging data and clinical questionnaires. First, the inference flow estimates output distributions via sampling in the latent space. Then, the simulation flow generates realistic perturbed samples of the view we want to study against others (here, the questionnaires) and infer digital avatars through the model. Finally, meaningful inter-view associations are inspected using hierarchical linear regressions.

D. 1

 1 Finding associations with a non-weighted r-DAA Algorithm 1 r-DAA using a simple average function for aggregation Input R: a set of ROIs M : a set of metrics S: a set of eCRF scores X: a set with N random train (2690) / left-out (301) splits1: p ← |R| * |M | * |S| ▷ num of features 2: n E ∈ {1, . . . , 20}▷ num of models 3: N ← 100 ▷ num of splits 4: π thr ← 0.4 ▷ stability path thr 5: for j ∈ {1, . . . , N } do 6:

  for 12: f ← mean function ▷ aggregation func 13: for n E ∈ {1, . . . , 20} do 14: for j ∈ {1, . . . , N } do 15:

  associations with a weighted r-DAA Algorithm 2 r-DAA using a RSA weighted average function for aggregation Input R: a set of ROIs M : a set of metrics S: a set of eCRF scores X: a set with N random train (2690)/ left-out (301) splits 1: p ← |R| * |M | * |S| ▷ num of features 2: n E ∈ {1, . . . , 20} ▷ num of models 3: N ← 100 ▷ num of splits 4: π thr ← 0.4 ▷ stability path thr 5: for j ∈ [1, N ] do 6:

  FigureS3: Investigation of the ROIs associated with SDQ in area. Each line corresponds to a ROI, dotted black ones are not selected as below the threshold π thr = 0.4 when using as ensembling function a mean weighted by model ratings. The ROIs selected are colored and the color is consistent across the plots. The red horizontal dotted line highlights the threshold π thr = 0.4. (a) Aggregated coefficients with a simple mean for a given split over the N = 100 splits, plotted against the number of models aggregated with the mean. (b) Feature wise stability Π when using an uniformly weighted mean to aggregate the coefficients output from r-DAA for each model, against the number of models used. (c) Same as the last one, except that the mean was weighted using model ratings. The latter strategy is used to select the ROIs associated in area with the SDQ score, with Π > π thr = 0.4.

Table 1 :

 1 If the prob-Representational similarity analysis results between the different latent spaces (eCRF, Joint, and ROI) and the clinical eCRF scores, and some covariates of interest : the imaging acquisition site (Site), the patient age (Age) and sex (Sex). τ is the corresponding average Kendall τ across models and splits. We computed the p-value associated with this statistic. Values in bold indicate their significativity, i.e. median corrected p-value < 0.01, see Supplementary B for details.

		eCRF	Joint	ROI
	Score	τ (↑)	τ (↑)	τ (↑)
	SRS	0.302 0.018 -0.003
	SCARED	0.101 0.032 -0.004
	ARI	0.256	0.05	0.017
	SDQ ha	0.326 0.058	0.005
	CBCL ab	0.406	0.03	0.008
	CBCL ap	0.443 0.026 -0.002
	CBCL wd 0.152	0.06	0.008
	Site	0.004	0.011	0.156
	Age	-0.005 0.05	0.10
	Sex	0.014	0.003	0.067

  n E ← 1

		N	N j=1 S j (n E )
	22: end for 23: S stable ← {k :	max n E ∈{1,...,20}	Π n E k > π thr }	▷ stable associations
	Output S stable

  E ← [distance(z i , z k ) for (i, k) ∈ {1, . . . , 301} 2 ] ∈ R 301×301Compute the Kendall τ j n E (s) ∈ R between upper triangular parts of C j n E and S j

	11: n 12: C j for s ∈ S do
	13:	R 301×301	Compute the left-out subjects pairwise eCRF score dissimilarity matrix S j n E ∈
	14: 16:	end for	S j n n E
	17: 18:	Compute average Kendall τ j n E ∈ R across the eCRF scores τ j n E ← 1 |S| s∈S τ j n E (s)
	19:	end for
	20: end for
	21: f ← weighted mean function	▷ aggregation func
	22: for n E ∈ {1, . . . , 20} do
	23:	for j ∈ {1, . . . , N } do
	24:

E ← [distance(s i , s k ) for (i, k) ∈ {1, . . . , 301} 2 ] ∈ R 301×301 15:

  Compute the selection probability of each feature Π n E ∈ R p

	26: 27: 28:	Compute decision support g(A j n E ) ∈ {0, 1} p S j (n E ) ← g(A j n E ) end for	1 n E	n E i=1 τ j i β j i
	29:
	30:	Π n E ← 1 N	N j=1 S j (n E )
	31: end for 32: S stable ← {k :	max n E ∈{1,...,20}	Π n E k > π thr }	▷ stable associations
		Output S stable

Data and Code availability

The code to reproduce the experiments and results is made available here. The data is accessible at https://fcon_ 1000.projects.nitrc.org/ .

Author Contributions

C.A. designed the experiments, C.A. conducted the experiments, V.F. et A.G. provided critical feedback, V.F. supervised the project, A.G. pre-processed the data, C.A, V.F. and A.G. redacted the manuscript, all authors reviewed the manuscript.

There are 2 K different subsets contained in the powerset P(x), and qϕ (z|x p ) = x k ∈xp q ϕ (z|x k ) is the product of Gaussian expert posteriors q ϕ (z|x k ), defined for the subset of views in x p . These formulae show our specific implementation of the MoPoE which handles multi-views by modelling viewspecific latent space posterior distributions, and a joint posterior distribution shared between views as a mixture of the products of their marginal posteriors (experts). This implementation differs from the one benchmarked in [START_REF] Daunhawer | On the limitations of multimodal vaes[END_REF][START_REF] Anonymous | Benchmarking multimodal variational autoencoders: Cdsprites+ dataset and toolkit[END_REF] which consider only one joint latent (MoPoE with its initial formulation [Sutter et al., 2021]). These benchmarks also favor the generative abilities when we essentially consider its representation space in our intrepretation framework.

Model specification and training

In our work, the encoders with parameters ϕ are each defined as multilayer perceptrons (MLPs), each with one hidden layer of 256 units and a ReLU activation. The decoders with parameters θ are linear (a single fully connected layer). Note that every decoder has a learnable variance parameter for each reconstructed feature that does not depend on the input observation. This allows the decoders to learn a population-level variability. To handle inputs of different sizes, the dimensions of the view-specific latent spaces were set individually. Based on a previous study [START_REF] Ambroise | Multi-modal latent variable model could help individuals stratification: application to hbn cohort[END_REF], we chose d eCRF = 1 and d ROI = 20. Defining the dimension of the shared latent space jointly with the dimensions (J 1 , . . . , J K ) of the input views seems to be a reasonable criterion, we use the following rule:

To spread the variability of the observations across the view-specific and shared latent representations, we choose d = 3. All data are z-scored before being fed into the network. The model is trained with an Adam optimization, a initial learning rate of 2 × 10 -3 and a batch size of 256.

B Representational Similarity Analysis

Assessing model quality in an unsupervised setting remains an open issue. To investigate the learned latent representations, we derive a Representational Similarity Analysis (RSA) [Kriegeskorte et al., 2008] between these representations and some measures on subjects (e.g. clinical scores or other covariates). We compute the subject-pairwise dissimilarity matrices in the latent space (modality-specific and joint) using the euclidean distance. We derive the same subject-pairwise dissimilarity matrices for the target measures. Finally, the Kendall rank correlation coefficient (Kendall τ ) enables the comparison of these dissimilarity matrices emphasizing the captured information. RSA is used twice in the present work. First we use it as a component of our framework to weight the contribution of a model when aggregating its extracted associations in the regularized digital avatar analysis (r-DAA) (see Section 2.3 and Supplementary alg. D.2). Second we consider it to evaluate the information contained in the different latent representations (see Section 3.1).

Model quality scoring used in r-DAA

We hypothesize that the latent representations of a good model to understand eCRF symptom scores should contain some amount of symptoms-related variability. We evaluate each trained MoPoE-VAE using RSA output Kendall τ between the model's joint latent space and each eCRF score (τ s ), defined as τ = (τ 1 , . . . , τ |S| ) ∈ [-1, 1] |S| . For each of the n E models, we use it's average Kendall τ = 1 |S| |S| s=1 τ s across eCRF scores as it quality score, used to weight the aggregation function in the ensembling of the n E corresponding DAAs (see section 2.3). This is detailed in the Supplementary alg D.2 as well.

Assessing the latent representation spaces with RSA We report in Table 1 aggregated results from RSAs applied to each split and corresponding models. We report the average N × n E Kendall τ between each different measure and latent space. In the table we report the significance of these correlations. We compute it by considering the p-values of the Kendall τ statistics (corresponding to the rejection of the null hypothesis). We correct them for multiple testing using the Hommel [HOMMEL, 1988] correction for multiple dependent tests and multiple independent test using Bonferroni. Then we consider the median of the corrected p-values. We chose the value of 1% for the significance threshold.