title: A Joint Network Optimization Framework to Predict Clinical Severity from Resting State Functional MRI Data
authors: Niharika Shimona D'Souza Mary Beth Nebel Nicholas Wymbs Stewart H. Mostofsky Archana Venkataraman
text: 
  
  
  
  
  
  
  
  wse]N.S. D'Souzacor1
  Shimona.Niharika.Dsouza@jhu.edu
  KKI]M.B. Nebel
  KKI]N. Wymbs
  KKI,Neu,Ped]S.H. Mostofsky
  wse]A. Venkataraman
  
  
  [cor1]Corresponding author
  
  [wse]Department of Electrical and Computer Engineering, Johns Hopkins University, USA
  [KKI]Center for Neurodevelopmental & Imaging Research, Kennedy Krieger Institute
  [Neu]Department of Neurology, Johns Hopkins School of Medicine, USA
  [Ped]Department of Pediatrics, Johns Hopkins School of Medicine, USA
  
  
  We propose a novel optimization framework to predict clinical severity from resting state fMRI (rs-fMRI) data. Our model consists of two coupled terms. The first term decomposes the correlation matrices into a sparse set of representative subnetworks that define a network manifold. These subnetworks are modeled as rank-one outer-products which correspond to the elemental patterns of co-activation across the brain; the subnetworks are combined via patient-specific non-negative coefficients. The second term is a linear regression model that uses the patient-specific coefficients to predict a measure of clinical severity. We validate our framework on two separate datasets in a ten fold cross validation setting. The first is a cohort of fifty-eight patients diagnosed with Autism Spectrum Disorder (ASD). The second dataset consists of sixty three patients from a publicly available ASD database. Our method outperforms standard semi-supervised frameworks, which employ conventional graph theoretic and statistical representation learning techniques to relate the rs-fMRI correlations to behavior. In contrast, our joint network optimization framework exploits the structure of the rs-fMRI correlation matrices to simultaneously capture group level effects and patient heterogeneity. Finally, we demonstrate that our proposed framework robustly identifies clinically relevant networks characteristic of ASD. 
  
  
  
  Matrix Factorization, Dictionary Learning, Functional Magnetic Resonance Imaging, Clinical Severity
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Resting State fMRI (rs-fMRI) is a non-invasive neuroimaging modality that captures steady-state patterns of co-activation in the brain in the absence of a task paradigm. It is believed that these correlation patterns reflect the intrinsic communication between brain regions [<cit.>]. Consequently, rs-fMRI has become ubiquitous in the characterization of neuropsychiatric disorders such as Autism Spectrum Disorder (ASD) [<cit.>], Attention Deficit Hyperactivity Disorder (ADHD) [<cit.>], and schizophrenia [<cit.>]. Traditional rs-fMRI analysis has concentrated on comparing the statistics of the rs-fMRI data, or variations in these statistics, across individuals or between different cohorts. For example, statistical differences in rs-fMRI features between a patient cohort and neurotypical controls have been considered as biomarkers of a particular disorder. However, the high dimensionality of rs-fMRI data, along with the considerable inter-patient variability, make it extremely difficult to reliably predict clinical manifestations on a patient-specific basis. 
  There has been considerable work in developing statistical methods to analyze rs-fMRI data. A large number of these studies build on standard multivariate [<cit.>] or random effects models [<cit.>] to capture population level differences in functional connectivity. Although these studies identify functional connections affected by the disease, they often fail to generalize on a patient-specific level. Additionally, these techniques do not  adequately characterize distributed impairments across multiple brain systems, which is crucial for understanding the complex pathologies associated with neuropsychiatric disorders [<cit.>]. This limitation has warranted the development of network-based models to study the inter and intra-subject variation across populations.
  Network-based rs-fMRI studies typically group voxels in the brain into regions of interest (ROIs) using a standard anatomical or functional atlas. Further, the synchrony between the average regional time courses is summarized using a similarity matrix, which is the input for further analyses. This extraction procedure is demonstrated in Fig. <ref>. The works of [<cit.>] use graph theoretic notions of connectivity based on aggregate network measures, such as node degree, betweenness centrality, and eigenvector centrality to study the functional organization of the brain. These measures are extremely useful to compactly summarize the connectivity information onto a restricted set of nodes which map to brain regions. A more global network property is small-worldedness [<cit.>], which describes an architecture of sparsely connected clusters of nodes. Changes in small-worldedness have been implicated in many neurological disorders [<cit.>,<cit.>]. These characterizations are quite successful at capturing global connectivity information, but often fail to illuminate the underlying etiological mechanisms. 
  To address the limitations of aggregate graph theoretic notions, recent focus has shifted towards mechanistic network models, which incorporate hierarchy onto existing graph connectivity notions. Community detection techniques are a class of population-level models which are used to identify highly interconnected subgraphs within a larger network. These techniques have become popular for understanding the organization of complex systems like the brain network architecture [<cit.>]. An application of this approach to identify regions having abnormal connectivity in schizophrenia patients can be found in [<cit.>]. Similarly, Bayesian community detection algorithms developed in [<cit.>] have provided valuable insights in characterizing the social and communicative deficits associated with autism. An alternative network topology is the hub-spoke model, which targets regions associated with a large number of altered rs-fMRI connections [<cit.>, <cit.>, <cit.>]. However, the above methods focus on group characterizations, and even studies that consider patient variability [<cit.>] have little generalization power on new subjects.  
  Machine learning techniques cast the neuroimaging prediction problem as a two stage procedure. Essentially, the first step is a feature selection or a representation learning stage, while the second stage uses the output of the first to predict the subject characteristics. A simple representation learning framework entails a careful sub-selection of specialized biomarkers [<cit.>]. On a whole brain level, data-driven approaches treat the patient connectivity information as a feature map and estimate lower dimensional projections, typically through PCA, kernel-PCA [<cit.>] or ICA [<cit.>]. From here, the most popular classifier (i.e. a stage two algorithm) is a Support Vector Machine (SVM) [<cit.>], which optimizes the decision boundary between patients and neurotypical controls [<cit.>]. SVMs have also been shown to identify disease sub-types [<cit.>] from the lower dimensional features with high accuracy. Along similar lines, the work of [<cit.>] proposes an ensemble learning based encoder-decoder model, that is able to provide competitive performance on large fMRI datasets at discriminative tasks. 
  While this two stage pipeline has been successful in the classification realm, characterizing finer-grained measures of clinical severity in the fMRI literature has been restricted to associative analysis, as opposed to an actual prediction on unseen data. For example, the work of [<cit.>] identifies key visual and motor ICA components, which are then used to compute a visuo-motor measure that is significantly correlated with social-communicative and motor deficit measures in ASD. In the context of a continuous value prediction, [<cit.>] develops a modified random forest regression algorithm for stacked multi-output score estimation from multiple ROI-voxel correlation maps. They demonstrate that it outperforms single score prediction. This strategy, however, does not permit a straightforward interpretation of the co-activation patterns explaining an individual severity score. Rather, it identifies regions that explain the complete set of scores jointly. Finally, deep learning methods have become popular for several neuroimaging data analysis. These models have the ability to efficiently learn complex abstractions of the input data without requiring careful feature engineering. As a result, they have been quite successful in a number of case/control classification tasks [<cit.>]. However, a downside to these models is the requirement of large amounts of training data for adequate generalization, which is rarely the case with clinical neuroimaging. Consequently, there has been limited success in predicting behavior from rs-fMRI data using neural networks. In summary, the unification of rs-fMRI and behavioral severity prediction, remains an open challenge.
  Dictionary learning [<cit.>] methods move away from the pipelined representations, and have recently gained traction due to their ability to simultaneously model both group level and patient specific information. The work of [<cit.>] proposed a correlation matrix decomposition strategy, in which, multiple rank one outer products capture an underlying generative basis. The sparse basis representation identifies meaningful co-activation patterns common to all the patients, while patient-specific coefficients combine the subnetworks and model the individual variability in the dataset. An extension of their work [<cit.>] looks at classification of young adults versus children, again, by the addition of an SVM like hinge loss. Our work builds on this representation by using the discriminative nature of these coefficients to predict their clinical severity via a linear regression penalty. This Joint Network Optimization (JNO) framework combines both a generative and discriminative term, as opposed to a pipelined hyperparameter search. The generalizability of the model is indicated by the regression performance on unseen data, instead of the correlation fit as used in [<cit.>]. This refinement demonstrates the potential of our JNO framework in identifying patient-predictive biomarkers of a given disorder. 
  We validate our framework on an rs-fMRI study of Autism Spectrum Disorder (ASD). Patients with ASD are known to manifest a wide spectrum of impairments in terms of social reciprocity, communicative functioning, and repetitive/restrictive behaviours [<cit.>]. This variation is typically quantified by a clinical severity measure obtained from an expert assessment. We find that our method outperforms several graph theoretic and machine learning feature representation techniques in predicting these severity scores on unseen data. Additionally, our model automatically extracts key networks commonly associated with altered functioning in the ASD literature. Finally, we quantify the merit of our joint objective function by comparing the combined predictive performance with that of a similarly defined two-stage decomposition and regression. We demonstrate that the joint objective bridges the gap in the two representative views of the data, thus aiding the resting state ASD characterization.
  
  A preliminary version of our work appeared in [<cit.>]. Here, we provide a detailed analysis of our model. We demonstrate the predictive performance of our algorithm on three different clinical severity measures, which capture varied social, behavioral and cognitive deficits associated with ASD. We evaluate our model on two clinical datasets to demonstrate the reproducibility of the method. We identify resting state networks explaining the different behavioral manifestations, and accordingly discuss the robustness of our brain basis characterization across behavioral measures. Lastly, we study hyperparameter sensitivity, generalizability in a test-retest setting and include mitigation strategies to improve the robustness of the framework. 
  
  
  
  § MATERIALS AND METHODS 
  
  
  
  
   §.§ A Joint Model for Connectomics and Clinical Severity
  
  Fig. <ref> presents a graphical overview of our model. The two inputs to our model are the rs-fMRI similarity matrices (upper left) and the scalar clinical severity scores for each patient (lower right). As mentioned earlier, Fig. <ref> illustrates the construction of the similarity matrix from the data. These matrices quantify the Pearson's Correlation Coefficient between the average time courses for each region of interest (ROI). The clinical scores are obtained from an expert evaluation and quantify the severity of the symptoms for the individual.
  Notice that the correlation matrices in Fig. <ref> have a dual representation. The generative part of the model is indicated in the purple box. Here, we decompose the correlation matrix into a basis term and a patient coefficient term. The columns of the basis capture ROI co-activation patterns common to the entire cohort, while the coefficients differ across patients and quantify the strength of each basis column in the matrix representation. The green box indicates the discriminative part of the model. Here, we leverage the information from the patient-specific coefficients to estimate a given measure of clinical severity via a linear regression model for each individual.
  
  
  
    
  Rs-fMRI Data Representation We define Γ_n∈ℛ^P × P as the correlation matrix for patient n, where P is the number of regions given by the parcellation. As seen in Fig. <ref>, we model Γ_n using a group average basis representation and a patient-specific network strength term. The matrix 𝐁∈ℛ^P × K is a concatenation of  K elemental bases vectors 𝐛_k∈ℛ^P × 1, i.e. 𝐁 := [𝐛_1  𝐛_2   ...   𝐛_K], where K ≪ P. These bases capture steady state patterns of co-activation across regions in the brain. While the bases are common to all patients in the cohort, the combination of these subnetworks is unique to each patient and is captured by the non-negative coefficients 𝐜_nk. We include a non-negativity constraint 𝐜_nk≥ 0 on the coefficients to preserve the positive semi-definite structure of the correlation matrices {Γ_n}. Our complete rs-fMRI data representation is:
  
      Γ_n≈∑_k𝐜_nk𝐛_k𝐛_k^T    s.t.     𝐜_nk≥ 0
  
  As seen in Eq. (<ref>), we model the heterogeneity in the cohort using a patient specific term in the form of 𝐜_n: = [𝐜_n1   ...   𝐜_nK]^T∈ℛ^K × 1. Taking diag(𝐜_n) to be a diagonal matrix with the K patient coefficients on the diagonal and off-diagonal terms set to zero, Eq. (<ref>) can be re-written in matrix form as follows:
  
      Γ_n≈𝐁diag(𝐜_n)𝐁^T    s.t.     𝐜_nk≥ 0
  
  Overall, this formulation strategically reduces the high dimensionality of the data, while providing a patient level description of the correlation matrices.
  
  
  
    
  Modeling Behavioral Scores As shown in the green box of Fig. <ref>, the patient coefficients {𝐜_nk} from the representation term, are used to model the clinical severity score 𝐲_n using a linear regression vector 𝐰∈ℛ^K × 1
  
      𝐲_n≈𝐜_n^T𝐰
  
  Concatenating the vectors 𝐜_n into a matrix 𝐂 := [𝐜_1   ...   𝐜_N] ∈ℛ^K × N, and the severity scores into a vector 𝐲 := [𝐲_1   ...   𝐲_N]^T∈ℛ^N × 1, Eq. (<ref>) can be equivalently represented in matrix form:
  
      𝐲≈𝐂^T𝐰
  
  
  
    
  Joint Objective for Representation and Prediction
  We combine the two contrasting viewpoints described above into a joint optimization function by summing the contributions of Eq. (<ref>) and Eq. (<ref>) below:
  
      𝒥(𝐁,𝐂,𝐰) = ∑_n||Γ_n - 𝐁𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T||_F^2
       + γ||𝐲 - 𝐂^T𝐰||^2_2  s.t.   𝐜_nk≥ 0 ,
  
  Here, ∑_n||Γ_n - 𝐁diag(𝐜_n)𝐁^T||^2_F is the total error in the representation of the N patient correlation matrices, and
  ||𝐲 - 𝐂^T𝐰||^2_2 is the prediction error for the behavioral data. Finally, γ is the trade-off between the rs-fMRI data-representation and score prediction terms. 
  
  
  
    
  Regularization Penalties Since we wish to capture a compact, yet clinically informative subnetwork representations, we add an ℓ_1 penalty to encourage sparsity in 𝐁. Intuitively, this regularizer will sub-select a small number of nonzero entries in 𝐁 that explain the data. From an optimization perspective, notice that scaled solution pairs {𝐁,𝐂} and {α𝐁, 1/α^2𝐂}, as well as {𝐂,𝐰} and {β𝐂, 1/β𝐰} give rise to equivalent data representations. As a result, we introduce a quadratic penalty on 𝐂 to act as a bound constraint. Similarly, we add an ℓ_2 regularization term to the regression vector 𝐰 analogous to ridge regression. Mathematically, the three regularizers can be written as:
  
      λ_1||𝐁||_1 + λ_2||𝐂||^2_F + λ_3||𝐰||^2_2
  
  The penalty terms in Eq. (<ref>) are added to the main objective in Eq. (<ref>). The final joint objective is as follows:
  
      𝒥(𝐁,𝐂,𝐰) = ∑_n||Γ_n - 𝐁𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T||_F^2
       + γ||𝐲 - 𝐂^T𝐰||^2_2
          + λ_1||𝐁||_1
       + λ_2||𝐂||^2_F + λ_3||𝐰||^2_2   s.t.   𝐜_nk≥ 0 ,
  
  The parameter λ_1 controls the number of nonzero elements in 𝐁 by scaling the contribution of the ℓ_1 penalty. Similarly, λ_2 and λ_3 relate to element wise bounds on the entries in 𝐂 and 𝐰 since they scale the contribution of their respective ℓ_2 norms. 
  
  
   §.§ Optimization Algorithm
  
  
  
  We employ an alternating minimization technique in order to infer the set of latent variables {𝐁,𝐂,𝐰}. Here, we optimize the JNO objective function from Eq. (<ref>) for each output variable, while holding the estimates of the other unknowns constant.
  Proximal gradient descent [<cit.>] is an attractive algorithm to handle the non-differentiable sparsity penalty on  𝐁 in Eq. (<ref>), when the supporting terms in the variable of interest are convex. However, from Eq. (<ref>), we see that the Frobenius norm terms expand to a biquadratic  representation in 𝐁, which is non-convex. We circumvent this problem by introducing N constraints of the form 𝐃_n = 𝐁diag(𝐜_n). We enforce these constraints using the Augmented Lagrangian [<cit.>], denoting the set of Lagrangian matrices by {Λ_n}. The modified objective function in Eq. (<ref>) takes the form:
  
      𝒥(𝐁, 𝐂,𝐰,𝐃_n,Λ_n) = ∑_n||Γ_n-𝐃_n𝐁^T||_F^2
       + γ||𝐲-𝐂^T𝐰||_2^2 + ∑_n[Λ_n^T(𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n))]
       + ∑_n1/2||𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n)||_F^2 + λ_1||𝐁||_1
       + λ_2||𝐂||^2_F + λ_3||𝐰||^2_2      s.t.   𝐜_nk≥ 0
  
  Such that [𝐌] is the trace operator, which sums the diagonal elements of the argument matrix 𝐌. The additional Frobenius norm terms ||𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n)||_F^2 act as regularizers for the trace constraints. Observe that Eq. (<ref>) is now convex in both 𝐁 and the set {𝐃_n}, which allows us to optimize them via standard procedures.
  Fig. <ref> provides an overview of the alternating minimization strategy employed. Each individual block in our optimization is described below. We refer the interested reader to Appendix A, which systematically delineates the supporting calculations from this section. 
  
  
    
  Proximal Gradient Descent on 𝐁
  Given the fixed learning rate parameter t, the proximal update for 𝐁 is:
  
      𝐁^k+1 = 𝐬𝐠𝐧(𝐋).^*(𝐦𝐚𝐱(|𝐋|-t,0)) 
          s.t.      𝐋 = 𝐁^k - (t/λ_1)∂𝒥/∂𝐁
  
  Here, -∂𝒥/∂𝐁 is a descent direction for the 𝐁 update, and t controls the magnitude of the step we take in this direction. In practice, we fix t at 10^-4 for stable convergence. The derivative of 𝒥 with respect to 𝐁, is computed as:
  
      ∂𝒥/∂𝐁 = ∑_n[2[𝐁𝐃_n^T𝐃_n-Γ_n𝐃_n]-𝐃_ndiag(𝐜_n)]
      
       +∑_n[𝐁diag(𝐜_n)^2-Λ_ndiag(𝐜_n)]
  
  At a high level, Eq. (<ref>) performs an iterative shrinkage thresholding operation to handle the non-smoothness of the ||𝐁||_1 using a locally smooth quadratic model.
  
  
    
  Optimizing C using Quadratic Programming
  The objective is quadratic in 𝐂 when 𝐁 and 𝐰 are held constant. Furthermore, the 𝐝𝐢𝐚𝐠(𝐜_n) term decouples the updates for 𝐜_n across patients. We use N quadratic solvers of the form given below to estimate the vectors {𝐜_n} : 
  
      1/2𝐜_n^T𝐇_n𝐜_n + 𝐟_n^T𝐜_n  s.t.   𝐀_n𝐜_n≤𝐛_n
  
  The objective and constraint matrices for our quadratic programming solvers are given by:
  
      𝐇_n = ℐ_K∘ (𝐁^T𝐁) + 2γ𝐰𝐰^T+ 2λ_2ℐ_K     
          𝐟_n = -2[ℐ_K∘(𝐃_n^T+Λ_n^T)𝐁]1  -2γ y_n𝐰;   
      𝐀_n = -ℐ_K   𝐛_n = 0
  
  Here, we use ∘ to denote the Hadamard product between two matrices and 1 to denote a vector of all ones. This strategy helps us find the globally optimal solutions for 𝐜_n. The constraint matrices 𝐀_n and 𝐛_n project the solutions onto the K dimensional space of positive reals.
  
  
    
  Closed Form Update for 𝐰 The global minimizer of 𝐰 can be computed by setting the gradient of Eq. (<ref>) equal to zero. Thus, the closed form solution for 𝐰 is given by: 
  
      𝐰 = (𝐂𝐂^T+ λ_3/γℐ_K)^-1(𝐂𝐲)
  
  This is analogous to a regularized linear regression , i.e. ridge regression update for 𝐰.
  
  
    
  Optimizing the Constraint Variables 𝐃_n and Λ_n
  Similar to the case of 𝐰, each of the primal variables {𝐃_n} has a closed form solution given by: 
  
      𝐃_n = (𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T+ 2Γ_n𝐁 - Λ_n)(ℐ_K+2𝐁^T𝐁)^-1
  
  We update the dual variables {Λ_n} via gradient ascent: 
  
      Λ_n^k+1 = Λ_n^k + η_k(𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n))
  
  The updates for 𝐃_n and Λ_n ensure that the proximal constraints are satisfied with increasing certainty at each iteration. The learning rate parameter η_k for the gradient ascent step  of the augmented Lagrangian is chosen to guarantee sufficient decrease for every iteration of alternating minimization. In practice, we initialize this value to 10^-3, and scale it by 0.5 at each iteration.
  
  
   §.§ Prediction on an unseen patient:
  
  
  
  In order to estimate the coefficients 𝐜̂ for a new patient, we re-solve the quadratic program in Eq. (<ref>) using the {𝐁^*,𝐰^*} computed from the training data via the procedure outlined in Section <ref>. We explicitly set the contribution from the data term in Eq. (<ref>) to 0, since the corresponding value of ŷ is unknown for the new patient. We also implicitly assume that the conditions for the proximal operator hold, i.e. the constraint 𝐃̂ = 𝐁^*𝐝𝐢𝐚𝐠(ĉ) is exactly satisfied. The estimation of the unseen patient's coefficients are thus mathematically formulated as follows:
  
      ĉ = _𝐜||Γ_n -  𝐁𝐝𝐢𝐚𝐠(𝐜)𝐁^T||_F^2 + λ_2||𝐜||^2_2
         s.t.   𝐜_k≥ 0
  
  Once again, Eq. (<ref>) can be formulated as a quadratic program. The parameters from Eq. (<ref>) correspond to:
  
      𝐇_n = 2(𝐁^T𝐁)∘(𝐁^T𝐁) + 2λ_2ℐ_K      
      𝐟_n = -2[ℐ_K∘(𝐁^TΓ_n𝐁)]1;             
              𝐀_n = -ℐ_K       𝐛_n = 0
  
  The estimate for the behavioral score for the test patient is given by the vector product ŷ= ĉ^T𝐰^*.
  
  
   §.§ Baseline Comparison Techniques
  
  
  We evaluate the performance of our method against a set of well established statistical, graph theoretic, and data-driven frameworks that have been used to provide rich feature representations. Fig. <ref> describes a general two stage pipeline for our task. The first stage is a representation learning step used for feature extraction. Stage 2 is a regression model to map the learned features to behavioral data. We evaluate our method against several choices of linear and non-linear algorithms for Stage 1. These are combined with a regularized linear regression in Stage 2, similar to our method. Additionally, we evaluate the performance obtained by omitting a Stage 1 and training a deep neural network end-to-end on the input correlation features. Lastly, we demonstrate the advantage provided by combining the neuroimaging and behavioral representations in the JNO framework. For this, we present a comparison where the feature learning and prediction stages are decoupled, similar to the baselines.
  
  
    §.§.§ Machine Learning Approach (PCA):
  
  
  We start with the P × P correlation matrix Γ_n for each patient. Since this matrix is symmetric, we have M =P × (P-1)/2 distinct rs-fMRI correlation pairs between various communicating sub-regions. Accordingly, the features from every individual are composed into a descriptor matrix 𝐗∈ℛ^M × N. We further concentrate these feature into a small number of representative bases. The basis extraction procedure in Stage 1 corresponds to a linear mapping in the original correlation space via a Principal Component Analysis (PCA). In Stage 2, we construct a regularized linear regression (ridge regression) on the projected features to predict the clinical severity. PCA projects the observations onto a set of uncorrelated principal component basis by means of an orthogonal linear decomposition. Mathematically, PCA poses the following dimensionality reduction problem:
  
      ℱ(𝐔,𝐙,μ) = _μ,𝐔,𝐘||𝐗-μ1^T-𝐔𝐙||^2_F
          s.t.   𝐔^T𝐔 = ℐ_d,   𝐙1 =0
  
  Here, 𝐔∈ℛ^N × d is the d dimensional subspace basis which best approximates the information from 𝐗 in the Frobenius norm sense, computed by calculating the eigenvectors of the sample covariance matrix 𝐗𝐗^T. Consequently, 𝐙∈ℛ^d × N is a compact d dimensional representation of 𝐗, where d ≪ M. 1 is a d dimensional vector of ones. The constraint 𝐙1 =0 centers the representation 𝐙.
  
  
    §.§.§ Statistical Approach (ICA)
  
  Here, we use Independent Component Analysis (ICA) as the Stage 1 algorithm combined with ridge regression. ICA operates on the raw time series data to extract representative spatial patterns that explain rs-fMRI connectivity. ICA has become ubiquitous for identifying group level as well as individual-specific connectivity signatures. It decomposes a multivariate signal into `independent' non-Gaussian components based on the statistics of the data. Mathematically, ICA models the components {𝐲_k} of the observed signal 𝐲 = [ 𝐲_1, … ,𝐲_m] as a sum of n independent components 𝐒 = [𝐬_1, … ,𝐬_n] combined via the mixing matrix 𝐀 = [ 𝐚_1, … ,𝐚_n]
  
      𝐲 = ∑_i=1^n𝐬_i𝐚_i    i.e.  𝐘 = 𝐀𝐒
  
  𝐬 can be recovered by multiplying the observed signals 𝐘 with the inverse of the mixing matrix 𝐖 =𝐀^-1. We adaptively estimate both the mixing matrix 𝐀 and the components 𝐬 by setting up a cost function that maximizes the non-gaussianity of 𝐬_i = 𝐰_i^T𝐲 or minimizes the mutual information.
  Group ICA extends this algorithm to a multi-subject analysis for extracting independent spatial patterns common across patients, but combined via individual time courses. We use the GIFT [<cit.>] software in order to perform Group-ICA to derive independent spatial maps for each patient. The correlation values between the identified components are fed to the regression model.
  
  
  
    §.§.§ Graph Theoretic Approach (Node Degree):
  
  
  Each correlation matrix Γ_n can be thresholded and considered a graph adjacency matrix, which we denote by Ψ∈ℛ^P × P. The element Ψ_ij gives the strength of association between two communicating sub-regions i and j. The underlying graph topology can be summarized using node/edge based importance measures [<cit.> <cit.>]. Again, we use a regularized linear regression technique to estimate the severity score from the reduced representation. This treatment closely parallels the machine learning approach, as we can view the graph measures as a dimensionality reduction. We compute Node Degree (D_N) from the adjacency graph followed by a ridge regression on the features.
  Given the adjacency matrix Ψ, the degree of region v is equal to the number of edges incident on v, with loops counted twice. Mathematically, the degree 𝐃_N(v) is computed as follows:
  
      𝐃_N(v) = ∑_j ≠ v1(Ψ_jv>0)
  
  where, 1(.) is the indicator function, which takes the value 1 if the condition is satified, and 0 otherwise. This metric captures the importance of each node in explaining the graph, which in our case, corresponds to the average connectivity strength of each region in the brain.
  
  
  
    §.§.§ A Neural Network Approach:
  
  
  Recently, there has been an upsurge in using neural networks to investigate neuroimaging correlates of developmental disorders [<cit.>]. Here, we test the efficacy of a simple Artificial Neural Network (ANN) for predicting the severity score from the correlation feature matrix 𝐗 defined above. The network architecture encodes a series of non-linear transformations of the input correlations to approximate the severity score. Recall that the size of the input is dependent on our choice of parcellation, which could be of considerable width (of the order of ≈5000 connections for P=100). After evaluating several architectures, we employ a two hidden layer network with widths 8000 and 10 respectively. We use a Rectified Linear Unit (ReLU) non-linearity after the first hidden layer and a Tanh non-linearity after the second hidden layer. We used the ADAM optimizer with an initial learning rate of 10^-4, scaled by 0.9 per 10 epochs, and a momentum of 0.9 to train the network.
  
  
  
  § EXPERIMENTS:
  
  
  
   §.§ Validation on Synthetic Data
  
  
  As a sanity check, we first sample data from the generative model in Eq. (<ref>) and use the optimization outlined in Section <ref> to estimate the unknowns {𝐁,𝐂,𝐰}. This procedure helps us analyze the performance of the algorithm under different noise scenarios. The inputs to our model are the correlation matrices {Γ_n} and the clinical scores {𝐲_n}. We note that the model gives a complete description of each Γ_n in terms of the basis vectors {𝐛_k} and the patient coefficients {𝐜_n}. Since the data representation terms for each patient are coupled solely through the basis representation, the coefficient descriptors are independent of each other. In a similar observation, each score 𝐲_n is explained by the corresponding 𝐜_n, independent of the remaining subjects, when we fix the regression vector 𝐰. We use this information to describe the observed data {Γ_n,𝐲_n} using a generative model with the likelihood model based on the hidden variables {𝐁,𝐂,𝐰}.
  
  Notice that, when treated as a Bayesian log-likelihood (i.e. taking a negative exponent of the objective), the ℓ_2 norms in Eq. (<ref>) translate into Gaussian distributions, and the ℓ_1 norm is equivalent to a Laplacian prior. The corresponding graphical model is shown in Fig. <ref>. The observed variables are indicated by the shaded circles. The white circles contain the hidden variables. The distribution parameters for the hidden variables are indicated in the corresponding rectangle pointing to the variable. The Laplacian parameter σ_B controls the overlap in the patterns of sparsity in 𝐁, which relates to λ_1. 𝐂 and 𝐰 are described by Gaussians with means zero (i.e. ℓ_2 norm offset). The variances σ^2_𝐂 and σ^2_𝐰 are related to the penalty parameters λ_2 and λ_3 respectively. The non-negativity constraint on 𝐜_n is handled by folding (i.e. taking the absolute value of) the normal distribution to restrict the 𝐜_n values to be positive reals. The observed variable {𝐲_n}, translates to a Gaussian with mean μ_𝐲_n=𝐜_n^T𝐰, and variance parameters σ_𝐲_n. This is again folded to reflect positive values of 𝐲_n. The correlation matrices {Γ_n} are drawn from a Gaussian distribution with mean μ_Γ_n= 𝐁𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T (which is positive semi-definite by construction) and variance σ_Γ_n.
  There are two sources of noise for the observed variables, which include the error in the correlation matrices Γ_n, and the error in the severity scores 𝐲_n. These scenarios can be directly related to controlling the variance parameters σ_Γ_n and σ_𝐲_n respectively. Additionally, we are interested in the performance of the algorithm under varying levels of overlap in the sparsity patterns in 𝐁.  
  We evaluate the performance using an average inner-product measure of similarity S between each recovered network, 𝐛̂_k, and its corresponding best matched generating network, 𝐛_k, both normalized to unit norm, i.e.:
  
      S =1/K∑_k|𝐛_k^T𝐛̂_k|/||𝐛_k||_2||𝐛̂_k||_2.
  
  Fig. <ref> depicts the performance of the algorithm in these three cases. The x-axis corresponds to increasing the levels of noise, while the y-axis indicates the similarity metric S computed for the particular setting. In the leftmost plot, an x-axis value close to 0 indicates high percentage of sparsity in 𝐁, while increasing values correspond to denser basis matrices. Throughout this experiment, the values of the other free parameters in the generative model were held constant. The middle plot evaluates subnetwork recovery when the noise in the scores, i.e. σ_𝐲_n is increased. The x-axis reports normalised values of σ_𝐲_n while the remaining free parameters were held constant. Similarly, the rightmost plot in Fig. <ref> indicates performance under varying noise in the correlation matrices Γ_n. Again, normalized σ_Γ_n values are reported on the x-axis. All numerical results have been aggregated over 100 independent trials.
  As expected, increasing the noise in the correlation matrices and scores worsens the recovery performance of the algorithm. This is indicated by the decay in the similarity measure with increasing noise parameters as well as an increase in the corresponding variance. Additionally, the algorithm performs better when there is lesser overlap in the columns of 𝐁, i.e. when the generating basis is sparse. However, we observe that our algorithm is robust in the noise regime estimated from the real-world rs-fMRI data (0.01-0.2) and recovered sparsity levels (0.1-0.4). In addition, we identify the stable parameter settings for the algorithm which guide our real world experiments.
  
  
  
   §.§ A Population Study of Autism :
  
  We evaluate the efficacy of our JNO framework on two separate cohorts. Our first dataset consists of 58 children with high functioning ASD acquired at the Kennedy Krieger Institute in Baltimore, USA. We refer to this as the KKI dataset. The age of the subjects ranged from 10.06 ± 1.26 and the IQ as 110±14.03. The second cohort is a subset of the publicly available Autism Brain Imaging Data Exchange (ABIDE I) [<cit.>] acquired at the New York University (NYU) site consisting of 63 patients. Social and communicative deficits in autism are believed to arise from abnormal interactions between functionally linked regions in the brain [<cit.>]. Therefore, identifying long-range correlative patterns in the brain directly linked to clinical severity is an important stepping stone to understanding and quantifying the neural underpinnings of the disorder. 
  
  
    
  Neuroimaging Data  
  For the KKI dataset, rs-fMRI scans were acquired on a Phillips 3T Achieva scanner using a single shot, partially parallel gradient-recalled EPI sequence with TR/TE = 2500/30ms, flip angle 70, res =3.05 × 3.15 × 3mm, having 128 or 156 time samples. The children were instructed to relax with eyes open and focus on a central cross-hair while remaining still for the duration of the scan.
  Slice time correction, rigid body realignment, and normalization to the EPI version of the MNI template was performed as a part of pre-processing using SPM [<cit.>]. The time courses were temporally detrended in order to remove gradual trends in the data. From here, we used a CompCorr [<cit.>] strategy for the estimation and removal of spatially coherent noise from the white matter and ventricles, along with the linearly detrended versions of the six rigid body realignment parameters and their first derivatives. We performed a spatial smoothing with a 6mm FWHM Gaussian kernel followed by a temporal filtering using a 0.01-0.1Hz pass band. Lastly, we removed spikes from the data via tools from the AFNI package [<cit.>] as an alternative to motion scrubbing.
  For the NYU cohort, the rs-fMRI data was pre-processed using the configurable pipeline for the analysis of connectomes, that has been integrated with ABIDE. The pre-processing steps involved are skull-stripping, global mean intensity normalization, spatial normalization to the MNI template, nuisance regression and CompCorr, followed by bandpass filtering, but without global signal regression. 
  This work relies on the Automatic Anatomical Labelling (AAL) atlas [<cit.>], which defines 116 cortical, subcortical and cerebellar regions. We compute a 116 × 116 correlation matrix for each patient based on the Pearson's Correlation Coefficient between the average time series for these regions. Empirically, we observed a consistent noise component with nearly constant contribution from all brain regions and low predictive power for both datasets. Therefore, we subtracted out the first eigenvector contribution from each of the correlation matrices and used the residuals as the inputs {Γ^t_n} to the algorithm and the baselines.
  
  
    
  Behavioral Data We analyzed three independent measures of clinical severity. These include:
  
  
      
    1 Autism Diagnostic Observation Schedule, Version 2 (ADOS-2) total raw score
      
    2 Social Responsiveness Scale (SRS) total raw score
      
    3 Praxis total percent correct score
   
  ADOS and SRS are standard assessments and are available for both the KKI and NYU datasets, while the Praxis has been collected for the KKI dataset only. 
  The ADOS consists of different sub-scores which  quantify the social and communicative deficits of the patient along with the restrictive/repetitive behaviors [<cit.>]. The test is administered by a trained clinician who evaluates the child against a set of guidelines. The total score is computed by adding the individual sub-scores. The dynamic range for ADOS is between 0 and 30, with higher score indicating greater impairment.
  The SRS scale characterizes the social responsiveness of an individual [<cit.>]. Typically, a parent/care-giver or teacher completes out a standardized questionnaire assessing various aspects of the child's behavior. SRS reporting tends to be more variable across patients compared to ADOS, since the responses are heavily biased by the parent/teacher attitudes. The SRS dynamic range is between 70-200 for ASD patients, with higher values corresponding to higher severity in terms of social responsiveness.
  Praxis was assessed using the Florida Apraxia Battery (modified for children) [<cit.>], which assesses ability to perform skilled motor gestures to command, to imitation, and with actual tool use. Several studies (<cit.>, <cit.>, <cit.>, <cit.> etc) reveal that children with ASD show marked impairments in Praxis i.e., developmental dyspraxia, and that impaired Praxis correlates with impairments in core autism social-communicative and behavioral features. Performance is videotaped and later scored by two trained research-reliable raters, with total percent correctly performed gestures as the dependent variable of interest. Scores therefore range from 0-100, with higher scores indicating better Praxis performance. This measure was available for 52 of the total patients in the KKI dataset. 
  
  
   §.§ Evaluating Predictive Performance
  
  We wish to compare the performance of our JNO framework against a wide class of algorithms described in Section <ref>. In this paper, we use regularized linear regression (i.e. ridge regression) for all baselines, except the ANN. Our Supplementary Results document includes further comparison with a non-linear (Random Forest) regression model. We emphasize that the baseline performance is nearly identical for both the linear and the non-linear models.
  
  We characterize the performance of each method using a 10 fold cross validation strategy as illustrated in Fig. <ref>. For a given parameter setting, we first split the data set into 10 training and test folds. For each of the folds, we train the models on a 90 percent training set split of the data. We report the score prediction on the held out 10 percent, which constitutes the testing set for that fold. Note that each datapoint is a part of the test set in exactly one of the 10 folds.
  We report two quantitative measures of performance. Median Absolute Error (MAE) quantifies the absolute distance between the measured and predicted scores: 
  
      MAE = median(|ŷ-𝐲|),
  
  where the median is computed across the set of patients. We report MAE values with the standard deviation of the error. Lower MAE indicates better testing performance.
  Normalized Mutual Information (NMI) assesses the similarity in the distribution of the predicted and observed score distributions across test patients. NMI is computed as follows:
  
      NMI(𝐲,ŷ) = H(𝐲) + H(ŷ) - H(𝐲,ŷ)/min{H(𝐲),H(ŷ) }
  
  where H(𝐲) denotes the entropy of 𝐲 and H(𝐲,ŷ) is the joint entropy between 𝐲 and ŷ. NMI ranges from 0-1 with higher values indicating a better agreement between predicted and measured score distributions, and thus characterizing improved performance.
  
  
  
  Along with MAE and NMI, we perform a statistical test on the error distribution to evaluate the performance gain of our framework. We first calculate the Cumulative Distribution Function (CDF) of errors and then report the Kolmogorov-Smirnov statistic on the CDF. This statistic helps us quantify the differences in the error distributions of each baseline compared to our algorithm. A Kolmogorov-Smirnov statistic lower than α = 0.05 is widely accepted in literature as statistically significant. We indicate comparisons which fall within this threshold in bold and near misses using an underline.
  
  
  
  
  
    §.§.§ Parameter Settings
  
  
  Our method has five user-specified parameters {γ,λ_1,λ_2,λ_3,K}. Recall that K is the number of basis networks, γ is the penalty tradeoff between the representation and regression terms, λ_1 is the sparsity penalty, while λ_2 and λ_3 are the regularization penalties on the coefficients 𝐂 and regression weights 𝐰 respectively.
  
  We use the knee point of the eigenspectrum of the correlation matrices Γ_n to select the number of bases (K=8 for both datasets). For reference, we have included the scree-plots in Fig. 8. Empirically, the JNO model is insensitive to the choice of λ_3 and γ, so we fix both at one. Effectively, we are left with two free parameters, which we optimize by performing a bivariate grid search. We note that the generalization accuracy is dependent on the dynamic range of the scores and is sensitive to λ_1 and λ_2. As described in Section 3.6, we have identified a stable range of operation across a single order of magnitude for these parameters. Based on the cross validation results, we finally use the following settings in our experiments: For the KKI dataset, {λ_2=0.2,λ_1=30} for ADOS, {λ_2=0.9,λ_1=50} for SRS and {λ_2=0.6,λ_1=20} for Praxis; for the NYU dataset, {λ_2=0.1,λ_1=20} for ADOS, {λ_2=0.9,λ_1=40} for SRS. 
  
  To provide a fair comparison with our JNO framework, we use a joint grid search on the Stage 1 hyperparameters and the Stage 2 ridge penalty to optimize these values for every baseline method. Again, we report the best performance in a ten fold cross validation setting.
  
  We select 10 PCA components for the KKI dataset, and 15 for the NYU dataset. For ICA, we obtained good performance for 35 spatial maps obtained from GIFT [<cit.>]. For the graph theoretic baseline, we threshold the correlation matrices {Γ_n} at 0.2 to obtain valid adjacency matrices {Ψ_n}. In conjunction with these, the ridge penalty parameter was swept across four orders of magnitude. 
  
  
  As described in Section <ref>, we fixed the network architecture to be a two hidden layer network with widths 8000 and 10 respectively, having a Rectified Linear Unit (ReLU) non-linearity after the first hidden layer and a Tanh non-linearity after the second hidden layer. We use a standard weight decay regularizer, with the regularization parameter varied over three orders of magnitude for each baseline comparison. We trained the network using the ADAM optimizer with an initial learning rate of 10^-4, scaled by 0.9 per 10 epochs, and a momentum of 0.9.
  
  Finally, we include the performance upon decoupling the ridge regression and the matrix decomposition in Eq. (<ref>) as a sanity check. This is akin to the two stage treatment in the baselines where the two terms are not explicitly coupled as in the JNO objective.
  
  
  
   §.§ Performance on Real World Data
  
  
  
  
  
  
  
  
   
  
   
  Figs. <ref>-<ref> compare the performance of our method against the baselines described in Section <ref> for the prediction of ADOS, SRS and Praxis respectively for the KKI dataset. Similarly, Fig. <ref> and Fig. <ref> illustrate the performance comparison on the NYU dataset for ADOS and SRS respectively. We plot the score predicted by the algorithm on the y-axis against the measured ground truth score on the x-axis. The bold x=y line indicates ideal performance. The red points correspond to training data, while the green points represent the held out testing data for all the folds in the cross validation. Our method is indicated at the top left corner in each plot. We observe that, although the training performance of the baselines is good (i.e. the red points follow the x=y line), the JNO achieves the best training performance in all cases. Furthermore, we notice that all the two stage baseline testing performances track the mean value of the held out data (indicated by the black horizontal line). Our method clearly outperforms the baselines and is able to capture a trend in the data, beyond a mean value estimation in case of both datasets for all scores. This can be observed by the spread of the green points about the x=y line in the case of the JNO method. Through our experiments, we noticed that the testing performance of the ANN is dependent on the choice of architecture. For example, the architecture chosen in Section <ref> performs well on predicting ADOS for the KKI dataset, but performs poorly on all other comparisons. Our empirical evaluations could not identify a single architecture that performed well in all cases, like our JNO framework. The failure of the two stage decomposition in the bottom right comparison figures strengthens our hypothesis that a joint modeling of the neuroimaging and behavioral data is necessary in the context of generalization onto unseen data.
  The lackluster generalization performance of the baselines is testament to the difficulty of the task at hand. The number of connections or features available to us are of the order of a 6670 dimensional vector representation for 58 or 63 patients. Both the machine learning and graph theoretic techniques we selected for a comparison are well known in literature for being able to robustly provide compact characterizations for high dimensional datasets. However, we see that PCA and ICA are unable to estimate a reliable projection of the data that is particularly indicative of clinical severity. Similarly, the node degree measure heavily rely on being able to accurately identify informative network topologies from the observed correlation matrices. However, its aggregate nature captures general trends and is not successful in characterizing subtle patient level differences. The failure of the decoupled matrix factorization and ridge regression makes a strong case for including the regression term as a part of our JNO objective. The basis directions obtained in this case are not indicative of clinical severity, due to which the regression performance suffers. Despite sweeping parameters across several orders of magnitude, we observe that the baselines are only good at capturing group level information, as is indicated by the training fit. However, they fail to characterize patient level differences for an unseen subject and simply predict the mean of the given cohort.
  On the other hand, the generalization power of the ANN is contingent on the model order choice. This is demonstrated by its inability to perform well on comparisons outside of ADOS for the KKI dataset. Said another way, we have to change the network architecture for different severity measures across datasets. This is a major computational disadvantage when compared with our method.
  A key difference between the JNO framework and the baselines is that we utilize the structure of the correlation matrices to guide the predictive model. In essence, we optimize for the tradeoff between the neuroimaging and behavioral data representations jointly, instead of posing it as a two stage problem. The matrix decomposition we employ explicitly models the group information through the basis, and the patient differences through the coefficients. The limited number of basis elements we employ to decompose the data provides us with compact representations which explain the connectivity information well. The regularization terms and constraints ensure that the problem is well posed, while providing clinically meaningful and informative representations about the data. We also quantify the performance indicated in these figures in Tables <ref> (KKI dataset) and Table <ref> (NYU dataset) based on the validation metrics mentioned earlier.
  
  
   §.§ Subnetwork Identification
  
  
  
  
  
  
  
  
  
  
  
  Figs. <ref>-<ref> illustrate the subnetworks in 𝐁, as trained on the ADOS, SRS and Praxis in the KKI dataset, respectively. Since each column of the basis corresponds to a set of co-activated subregions, we plot the normalized values stored in these columns onto the corresponding AAL ROIs. The colorbar indicates subnetwork contribution to the AAL regions. Regions colored as negative values are anticorrelated with regions storing positive ones. We rank the 8 subnetworks obtained from SRS and Praxis according to their overlap with the subnetworks from ADOS. As seen from these figures, corresponding subnetworks show considerable overlap in regional co-activation patterns. The individual variations can arise from the fundamental differences in the behavioral traits that each score is trying to capture.
  From a clinical standpoint, Subnetwork 7 includes competing i.e. anticorrelated contributions from regions of the default mode network (DMN) and somatomotor network (SMN). Abnormal connectivity within the DMN and SMN has been previously reported in ASD [<cit.>]. Subnetwork 5 comprises of competing contributions from SMN regions. Additionally, it includes higher order visual processing areas in the occipital and temporal lobes, which is consistent with behavioral reports of reduced visual-motor integration in the ASD literature [<cit.>]. Subnetwork 1 has competing  from prefrontal and subcortical contributions, mainly the thalamus, amygdala and hippocampus. The thalamus is responsible for relaying sensory and motor signals to the cerebral cortex in the brain. The hippocampus is known to play a crucial role in the consolidation of long and short term memory, along with spatial memory to aid navigation. Altered memory functioning has been shown to manifest in children diagnosed with ASD [<cit.>]. Along with the amygdala, which is known to be associated with emotional responses, these areas may be crucial for social-emotional regulation in ASD.  Finally, Subnetwork 2 is comprised of competing contributions from the central executive control network and the insula, which is thought to be critical for switching between self-referential and goal-directed behavior [<cit.>]. 
  
  
    §.§.§ Robustness in Subnetwork Recovery
  
  Notice that we estimate a different basis matrix 𝐁 for each cross validation fold. Therefore, one important property to verify is that these subnetworks are similar across different cohorts of the data.
  We observed an average similarity of 0.79  0.06 for the ADOS networks, 0.86  0.04 for the SRS networks, and 0.76  0.06 for the Praxis networks across their cross validation runs. Additionally, upon a cross comparison between the ADOS and SRS networks, we obtained an average similarity of 0.82  0.07. Similarly, the overlap between ADOS and Praxis is 0.79  0.04, and between SRS and Praxis is 0.77  0.06. For a convenient visual inspection, we have arranged the networks in Fig. <ref> (SRS) and Fig. <ref> (Praxis) in the order of their inner product similarity with the ADOS networks in Fig. <ref>. This finding strengthens the hypothesis that our model is successful at capturing the stable underlying mechanisms which explain the different sets of deficits of the disorder.
  
  
  
    §.§.§ Comparing Subnetwork Representations
   In this section, we compare the subnetworks identified by  the JNO to the representations learned by the baseline methods. Recall that we have used a regularized linear regression as the Stage 2 predictor for the baselines. Therefore, we can probe the learned regression weights to characterize the baseline network representations.
  
  Degree centrality looks at the relative importance of each brain region or `node' to the overall representation. To visualize the pattern identified by the degree centrality + ridge regression baseline, we display the regression weights on the brain surface plots in Fig. <ref>, normalized to unit norm. The colorbar indicates the strength of co-activation. Regions storing negative values are anticorrelated with regions storing positive weights. We again observe patterns from the DMN in the subnetwork plot. Note that the DMN was also a key connectivity pattern identified by the JNO. However, several other subnetworks identified by the JNO do not figure in this representation.
  
  On the other hand, for the PCA + ridge regression baseline, the regression weights inform us of the relative importance of the principal components in prediction. Since the features fed into PCA are the M = (P×(P-1))/2 correlation values, we are left with a 6670 dimensional edge connectivity representation for the AAL per component. We first examine the absolute value of the regression weights learned, and then display the connectivity in the top 2 basis components in Fig. <ref>. We render this connectivity measure using the BrainNet Viewer [<cit.>] software. For clarity, we have chosen to display the top 5 percent of the connections obtained. The solid edges signify retained connections, while the blue spheres correspond to nodes of the AAL regions. The colorbar to the right indicates the strength of the connections. We notice that the components consist of several crossing connections spread across different regions of the brain. As compared to our model, which pinpoints key subnetworks already known to be associated with ASD, the representation obtained is not immediately interpretable.
  
  In the ICA + ridge regression baseline, the input to the regression model are the correlation values between the components identified by ICA. After the model is fit, we sort the input correlations based on the learned regression weights. This helps us identify the features important for prediction. In Fig. <ref> , we display the spatial maps of the top 2 connections identified by the algorithm. We again, observe patterns from the DMN and visual areas. However, it fails to capture several other subnetwork patterns that the JNO identifies as important for ASD.
  
  For the ANN, we use the weight matrix learned at the input layer to inform us of the subnetwork connectivity. Recall that this matrix is of dimension M × D, where M = (P × (P-1))/2 = 6670 for the AAL atlas. For our application, D is of width 8000. We first take the absolute values of these weights, and then normalize the columns of this matrix to unit norm. We then average across the rows to obtain a single 6670 dimensional edge-edge connectivity vector. Again, we use the BrainNet connectivity plots to display this information in Fig. <ref>. We have chosen to display the top 1 percent of the connections obtained. The solid edges signify retained connections, while the blue spheres correspond to nodes of the AAL regions. The colorbar to the right indicates the strength of the connections. We observe several overlapping connectivity patterns spread across the entire brain despite applying a stringent threshold. Additionally, the narrow range of values indicates that the ANN assigns nearly equal weight to all connections on an average. Similar to the PCA baseline, this representation is unable to capture interpretable connectivity patterns which explain behavior.
  
  Finally, we examine the representation learned by performing the matrix decomposition and prediction separately, i.e. the decoupled case. Note that the learned basis matrix 𝐁 follows the same interpretation as that of the JNO. We display the corresponding co-activation patterns in Fig. <ref>. Again, the colorbar indicates the strength of activation of the AAL ROIs. Negative regions are anticorrelated with the positive regions. For convenience, we have ordered the 8 subnetworks according to their similarity with the ADOS subnetworks identified in Fig. <ref>. Since we use the same matrix decompotion as the JNO, we observe several similarities in the learned representations. We also notice subtle differences in the patterns on account of the coupling with the predictive term in the JNO. We conjecture that these learned differences are what gives the JNO the leverage to generalize to unseen data.
  
  
  
  
   §.§ Evaluating Model Generalizability
   
  
  
  We have shown both predictive power and interpretabilitiy of our model thus far. Furthermore, characterizing model generalizability is important for future application of our framework. Here, we first examine the sensitivity of our prediction results with respect to the model hyperparameters. Then, we show robustness of our formulation to two common obstacles in generalizabilitiy for rs-fMRI analysis, namely choice of parcellation scheme and test-retest reliability. Lastly, we show mitigation strategies to handle hyperparameter sensitivity that make our framework more robust. 
  
  
  
    §.§.§ Hyperparameter Sensitivity
  
  As initially described in Section <ref>, our JNO framework is insensitive to the regression tradeoff γ and ridge penalty λ_3. We also have a natural way to set the number of subnetworks K. However, we observe that our JNO framework is fairly sensitive to the sparsity on 𝐁 and the ridge penalty on the coefficients 𝐜_n, i.e. λ_1 and λ_2 respectively. Fig. <ref> represents the MAE recovery performance of the algorithm for varying settings of λ_1 and λ_2, holding the remaining parameter settings constant when evaluated on the KKI dataset. The red plots in each case indicate the performance of the JNO framework. The green plots allude to our mitigation strategy which we will present later in Subsection <ref>. The 𝐱-axis denotes the parameter value, while the 𝐲-axis quantifies the MAE from cross validation. Observe that the best λ_1 and λ_2 settings for the individual scores are different, i.e ADOS-{λ_1=30, λ_2=0.2}, SRS-{λ_1=50, λ_2=0.9}, and Praxis-{λ_1=20, λ_2=0.6}. Additionally, the kinks in the plots (shown by the black arrow) also indicate that small changes in the sparsity and coefficient regularization lead to a dramatic change in performance, i.e. the operating points for these two parameters are narrow. We suspect that the hyperparameter differences can be partially attributed to the different dynamic range of each clinical score. Specifically, these differences impact the tradeoff between the representation learning and prediction terms in the JNO optimization. This in turn affects the generalization performance at a particular hyperparameter setting. These observations further illustrate the difficulty of the problem we are trying to address. Keeping these subtleties in mind, the next subsection focuses exclusively on studying the generalization performance of the JNO and the baselines at the optimal hyperparameter settings we obtained in Section <ref>.
  
  
  
  
  
    §.§.§ Examining Test-Retest Performance:
  
  Based on the observations above, we design two experiments to characterize the generalizability of the JNO with regards to its free parameters, i.e. λ_1 and λ_2, in a test-retest setting. The first of these experiments investigates the impact of varying the rs-fMRI parcellation scheme, which not only changes the dimensionality of the input correlation matrices, but also changes the distribution of correlation measures. The second experiment is a cross site comparison between KKI and NYU where the same parcellation (i.e. AAL) is maintained. However, site differences like acquisition protocol and pre-processing, impact the distribution of the input correlations. Overall, these two experiments validate generalizabilitiy and robustness of our model, as we show good performance under sub-optimal parameter settings.
  
  
  
  
  
  
  
  For the first experiment, we compute correlation matrices for the KKI rs-fMRI dataset using the Brainetome-246 atlas. Note that the Brainetome [<cit.>] provides a much finer ROI resolution as compared to the AAL-116 atlas used previously. We then predict the clinical scores from the new correlation matrices using the AAL hyperparameter settings. This is repeated for all the baselines along with the JNO. The plots in Figs. <ref>-<ref> illustrate the prediction of ADOS, SRS and Praxis respectively by the JNO. In each figure, the left plots denote the Brainetome performance using the hyperparameter selection outlined in Section <ref>, while the right plots indicate the performance after transferring the hyperparameter settings from the AAL experiments. Notice that the performance is slightly worse than the case where they were explicitly set by cross validation (left plots). This is expected since we have not explicitly optimized for the new data characteristics. Table <ref> provides the corresponding quantitative metrics. We emphasize that the performance trends are very similar to those in Section <ref>, namely, we outperform nearly all the baselines, while retaining the interpretability of the subnetwork decomposition. Since the ROI resolution differs here, we can treat this comparison as a surrogate for out of sample generalization. Additionally, this experiment suggests that the JNO is agnostic to the parcellation scheme used. We have included additional baseline comparisons in the supplementary results document. We believe that the joint optimization helps us balance the tradeoff between the representation learning and prediction terms. Said another way, the matrix factorization term regularizes the problem and helps provide stability while handling the changes in the input data distribution. As a result, we are able to transfer the learned hyperparameter settings, yet outperform the baselines at severity prediction.
  
  
  In our second experiment, we again use the hyperparameter settings learned from the KKI dataset, but this time, predict clinical scores from the NYU rs-fMRI correlation matrices. The parcellation scheme i.e. AAL, is maintained for both datasets. Figs. <ref> and <ref> compare these predictions (right sub-plots) with those from Section <ref> (left sub-plots) for ADOS and SRS respectively. Table <ref> delineates the quantitative comparisons against the baselines. Again, we observe that the JNO still outperforms the baselines, though there is a slight reduction in quantitative performance. This reduction is expected, especially given the differences in scanning protocols across sites. As in the previous experiment, we are able to use the transferred parameter settings, and still optimize for a representation which is predictive of clinical severity. The cross dataset comparison is a good indicator of the robustness of the JNO to its hyperparameters. The observations from both experiments suggest that the JNO has the potential to capture robust and interpretable phenomenon related to the neurological disorder of interest. 
  Lastly, both the test-retest experiments study the overall generalizability of the method by examining a form of out of sample prediction, either by varying the parcellation scheme, or via a cross dataset comparison. This is very similar to the principles of nested cross validation, which aims to reduce bias in parameter selection.
  
  
  
  
  
  
  
    §.§.§ Mitigating Parameter Sensitivity
  
  
  Finally, we propose two main modifications to tackle the observed hyperparameter sensitivity in λ_1 and λ_2. Given that the dynamic ranges of the scores are quite different and potentially impact generalization, our first mitigation strategy is to rescale the measures to a fixed interval. Since ADOS is the most widely accepted observational measure of clinical autism severity, we have scaled and offset the remaining scores to have a range of 0–30 (similar to ADOS).  To mitigate the narrow `operating point', we include an extra template average correlation term in Eq. (<ref>). We now model the residual outer-product terms as deviations around a mean template correlation matrix 𝐁_avg. The rationale behind this additional term is that it encourage sparsity in the basis matrix along with the explicit ℓ_1 penalty. The modified objective is as follows:
  
      𝒥(𝐁,𝐁_avg,𝐂,𝐰) = ∑_n||Γ_n - 𝐁_avg -𝐁𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T||_F^2
       + γ||𝐲 - 𝐂^T𝐰||^2_2 + λ_1||𝐁||_1
       + λ_2||𝐂||^2_F + λ_3||𝐰||^2_2  s.t.   𝐜_nk≥ 0 ,
  
  
  
  Notice that 𝐁_avg has a closed form update, which does not add much computational overhead. The updates for the remaining variables follow the same procedure as described in Section. <ref>, except that the term, {Γ_n} is replaced with {Γ_n-𝐁_avg} in every update.
  The green plots in Fig. <ref> illustrate the cross validated performance of the modified JNO framework from Eq. (<ref>). The operating point {λ_1,λ_2} for the modified framework is fairly consistent across the scores. Moreover, the green plots exhibit a larger stable range (highlighted in yellow) compared to the red plots. Accordingly, we identify the settings {λ_1 = 10–30, λ_2 = 0.08–0.6} as the operating range for the modified JNO objective, which is roughly an order of magnitude larger than the original formulation and does not exhibit any kinks. Fig. <ref> and Fig. <ref> illustrate the best generalization performance for SRS and Praxis using the two algorithms.
  
  Notice that the modified JNO has a slight tradeoff in regression performance at the expense of the gain in parameter stability. We highlight the importance of this exploration, as future applications of our work include applying our method to rs-fMRI and severity scores from a variety of neurological disorders. Our modified formulation provides additional flexibility in this sense, and extends the overall generalizability of our model.
  
  
  
  
  § DISCUSSION
  
  Our JNO model cleverly exploits the structure intrinsic to rs-fMRI correlation matrices through an outer product representation. The regression term further guides the basis decomposition to explain the group level and patient specific information. The compactness of our representation serves as a dimensionality reduction step that is related to the clinical score of interest, unlike the pipelined treatment commonly found in the literature. As seen from the results, our JNO framework outperforms a wide range of well established baselines from the machine learning and graph theoretic methods ubiquitous in fMRI analysis on two separate real world datasets.
  We conjecture that the baseline techniques fail to extract representative patterns from the correlation data, and learn only the group level representation for the cohort. Consequently, they overfit the training set, despite sweeping the parameters across several orders of magnitude. Any patient level symptomatic and connectivity level differences are lost due to the restrictive pipelined procedure and the group level confounds. 
  Our Joint Network Optimization Framework is agnostic to the choice of parcellation scheme. This was demonstrated by our additional experiments on the KKI dataset, where we chose the 246 region Brainnetome parcellation to extract correlation matrices (Section <ref>). We emphasize that our framework makes minimal assumptions on the data. Provided we have access to a valid behavioral and network similarity measure, this analysis can be easily adapted to other neurological disorders and even predictive network models outside the medical realm. This greatly broadens the scope of the method to numerous potential applications. 
  Finally, notice that the training examples (red points) in Figs. <ref>-<ref> follow the 𝐱=𝐲 line nearly perfectly. Here, we explain this (potentially misleading) phenomenon in terms of the parametrizatization of our joint objective in Eq. (<ref>).
  Recall that Section <ref> describes the procedure for calculating the coefficients for an unseen patient 𝐜̅_n from the training solution set {𝐁^*,𝐰^*}. Recall that we explicitly set the contribution from the data term in Eq. (<ref>) to 0. Since the patient is not a part of the training set, the corresponding value of ŷ is unknown. In contrast, the training performance is computed based on the estimated coefficients 𝐜_n, which have access to the severity scores. Here, we examine the effect of removing the severity information when calculating the coefficients for the training patients. In other words, we estimate the corresponding severity 𝐲 excluding the ridge regression term. Accordingly, Fig. <ref> highlights the differences in training fit with and without this term is not included in estimating 𝐜_n. Notice that in the latter, the training accuracy has the same distribution as the testing points in Figs. <ref>-<ref>. Taken together, we conclude that, the linear predictive term overparamterizes the search space of solutions for 𝐜_n to yield a near perfect fit. We use this observation to emphasize that the subnetworks and regression model learned by our JNO framework are capturing the underlying data distribution and not simply `overfitting' the training data.
  
  
  
  § CONCLUSION
  
  
  We present an elegant matrix decomposition strategy to combine neuroimaging and behavioral data information. As opposed to generic prediction frameworks, our model directly captures key representative information from the correlation matrices. In contrast, conventional analysis methods dramatically fall short of unifying the two data viewpoints reliably enough to implicate predictive functional patterns in the brain. Our joint optimization framework robustly identifies brain networks characterizing ASD and provides a key link to quantifying and interpreting the clinical spectrum of manifestation of the disorder across patients. Moreover, our evaluation on two separate real world dataset supports the reproducibility of the framework.
  We are working on a multi-score extension which can incorporate data from different behavioral domians. In the future, we will explore extensions of this model that learn a patient versus controls distinction in addition to predicting symptom severity. A potential extension of this model includes replacing the linear regression term in Eq. (<ref>) with its non-linear counterpart, thus providing us with the flexibility to model more complex decision functions which can better map the behavioral space.
  Another avenue for refinement is to incorporate structural connectivity information in the form of anatomical priors. Typically, structural modalities such as Diffusion Tensor Imaging (DTI) are used to define and track existing anatomical pathways in the brain. Incorporating this information into the network optimization model could be an important step towards unifying anatomical, functional and behavioural domains to better understand altered brain functioning in the context of neurological disorders such as Autism, ADHD, and Schizophrenia.
  Our experiment on the KKI dataset using the Brainetome parcellation supports the observation that the method can be tuned with different connectivity and behavioral information, which is an added flexibility. Thus, it could be used to characterize the efficacy of behavioral therapies for neuropsychiatric disorders, as well as, for the development of patient-specific disease biomarkers. We believe that the all of benefits offered by our JNO framework could make it an important diagnostic tool for personalized medicine in the future. 
  
  
    
  Acknowledgements
  This  work  was  supported  by the National Science Foundation CRCNS award 1822575 and CAREER award 1845430, the National  Institute  of Mental Health (R01 MH085328-09, R01 MH078160-07, K01 MH109766 and R01 MH106564), the National Institute of Neurological Disorders and Stroke (R01NS048527-08), and the Autism Speaks foundation.
  
  
  § APPENDIX A
  
  
  In this section, we derive the alternating minimization updates from Section <ref>:
  
  
    
  Optimizing 𝐁 via Proximal Gradient Descent We first write out the optimization problem with respect to 𝐁 when the estimates of {𝐂,𝐰} are held constant:
  
      𝐁^k+1 = _𝐁λ_1||𝐁||_1+∑_n||Γ_n-𝐃_n𝐁^T||_F^2
      
      + ∑_n[Λ_n^T(𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n))]
       + ∑_n1/2||𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n)||_F^2        
      𝐁^k+1 = _𝐁||𝐁^k||_1+ 1/λ_1𝒢(𝐁^k)
  
  We see that the proximal gradient iteration is the solution to the following fixed point problem:
  
      0∈1/λ_1∂𝒢/∂𝐁 + ∂(||𝐁||_1)
  
  The derivative of 𝒢 with respect to 𝐁, is computed as:
  
      ∂𝒢/∂𝐁 = ∑_n[2[[𝐁𝐃_n^T-Γ_n]𝐃_n]-𝐃_ndiag(𝐜_n)] 
       +∑_n[𝐁diag(𝐜_n)^2-Λ_ndiag(𝐜_n)]
  
  Given the fixed learning rate parameter t, the proximal update for 𝐁 is easily computed as:
  
      𝐁^k+1 = 𝐬𝐠𝐧(𝐗)∘(𝐦𝐚𝐱(|𝐗|-t,0)) 
      𝐗 = 𝐁^k - (t/λ_1)∂𝒢/∂𝐁
  
  This step first estimates a locally smooth quadratic model at each iterate 𝐁^k and applies a step of iterative shrinkage thresholding to the compute the local solution of 𝐁. The resulting iterative algorithm is computationally efficient compared to the counterpart sub-gradient based descent methods and arrives at a good local solution for an appropriate choice of the learning rate.
  
  
    
  Optimizing 𝐂 using Quadratic Programming
  The objective is quadratic in 𝐂 when 𝐁, and 𝐰 are held constant. Furthermore, the 𝐝𝐢𝐚𝐠(𝐜_n) term decouples the updates for 𝐜_n across patients. Each 𝐜_n is the solution to the a separate optimization problem of the following form:
  
      𝐜^k+1_n = _𝐜_n∈ℛ^K+[Λ^T_n (𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜^k_n))] +λ_2||𝐜^k_n||^2_2
       + 1/2||𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜^k_n)||^2_F  + γ((𝐜_n^k)^T𝐰-𝐲_n)^2
  
  Hence, we use N quadratic programs (QP) of the form below to solve for the vectors {𝐜_n} : 
  
      1/2𝐜_n^T𝐇_n𝐜_n + 𝐟_n^T𝐜_n  s.t.   𝐀_n𝐜_n≤𝐛_n
  
  The QP parameters for our problem are given by:
  
      𝐇_n = ℐ_K∘ (𝐁^T𝐁) + 2γ𝐰𝐰^T+ 2λ_2ℐ_K     
           𝐟_n = -2[ℐ_K∘(𝐃_n^T+Λ_n^T)𝐁]1  -2γ y_n𝐰;   
      𝐀_n = -ℐ_K   𝐛_n = 0
  
  The non-negativity constraint requires us to project the quadratic programming solution to the space of positive reals in K dimensions for each 𝐜_n through 𝐀_n and 𝐛_n. Since the Hessians {𝐇_n} for our problem are positive definite, there exist polynomial time algorithms for solving the bound constrained QPs to the global optimum value. The decoupling of the {𝐜_n} allows us to solve for each coefficient vector in parallel.
  
  
    
  Closed Form Update for 𝐰 The global minimizer of 𝐰 is computed at the first order stationary point of the convex objective, which is:
  
      𝒥(𝐰) = λ_3||𝐰||^2_2 + γ||𝐂^T𝐰-𝐲||^2_2  
                 ∂𝒥/∂𝐰 = 0 = 2 λ_3𝐰 + 2 γ (𝐂𝐂^T𝐰 - 𝐂𝐲)
  
  The closed form update can be expressed as:
  
      𝐰 = (𝐂𝐂^T+ λ_3/γℐ_K)^-1(𝐂𝐲)
  
  Thus, the ratio λ_3/γ acts as a regularizer for the matrix inversion in our estimate, ensuring that the update for 𝐰 is well defined at each iterate. 
  
  
    
  Optimizing the Constraint Variables 𝐃_n and Λ_n
  A closed form solution for the primal variables {𝐃_n} can be obtained by setting their first derivatives to zero: 
  
      ∂𝒥/∂𝐃_n = 0 = 𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T + 2Γ_n𝐁
       - Λ_n - 𝐃_n  -2 𝐃_n𝐁^T𝐁
  
  
      𝐃_n = (𝐝𝐢𝐚𝐠(𝐜_n)𝐁^T+ 2Γ_n𝐁 - Λ_n)(ℐ_K+2𝐁^T𝐁)^-1
  
  The gradient ascent update on  {Λ_n} is as follows:  
  
      ∂𝒥/∂Λ_n = 𝐃_n-𝐁𝐝𝐢𝐚𝐠(𝐜_n)
      
      Λ_n^k+1 = Λ_n^k + η_k∂𝒥/∂Λ_n
  
  Similar to the case of the coefficients 𝐜_n, each of the N pairs of updates {𝐃_n,Λ_n} are decoupled from each other, and can be solved in parallel. Overall, the sets of Λ_n gradient ascent updates ensure that the respective set of constraints 𝐃_n= 𝐁diag(𝐜_n) is satisfied with increasing certainty at each iteration. The Augmented Lagrangian construct ||𝐃_n-𝐁diag(𝐜_n)||^2_F prevents trivial Lagrangian Λ_n solutions. 
  
  
  
  