Exploring brain connectivity changes in major depressive disorder using functional-structural data fusion: A CAN-BIND-1 study.

There is a growing interest in examining the wealth of data generated by fusing functional and structural imaging information sources. These approaches may have clinical utility in identifying disruptions in the brain networks that underlie major depressive disorder (MDD). We combined an existing software toolbox with a mathematically dense statistical method to produce a novel processing pipeline for the fast and easy implementation of data fusion analysis (FATCAT-awFC). The novel FATCAT-awFC pipeline was then utilized to identify connectivity (conventional functional, conventional structural and anatomically weighted functional connectivy) changes in MDD patients compared to healthy comparison participants (HC). Data were acquired from the Canadian Biomarker Integration Network for Depression (CAN-BIND-1) study. Large-scale resting-state networks were assessed. We found statistically significant anatomically-weighted functional connectivity (awFC) group differences in the default mode network and the ventral attention network, with a modest effect size (d < 0.4). Functional and structural connectivity seemed to overlap in significance between one region-pair within the default mode network. By combining structural and functional data, awFC served to heighten or reduce the magnitude of connectivity differences in various regions distinguishing MDD from HC. This method can help us more fully understand the interconnected nature of structural and functional connectivity as it relates to depression.

INTRODUCTION

Increasing interest in brain connectivity patterns in illness and in health has given rise to the development of multimodal imaging analysis approaches, utilized to combine functional magnetic resonance imaging (fMRI) data with diffusion tensor imaging (DTI) data (Reddi, 2017; Zhu et al., 2014). Multimodal imaging analysis methods aim to capture the complex interactions between structural and functional connectivity in brain networks and provide new insights into brain connectivity. Complex and heterogeneous disorders such as major depressive disorder (MDD) can benefit from multimodal imaging analysis, which offers a better understanding of the joint structural and functional changes in human brain connectivity patterns (Zhu et al., 2014).

MDD is associated with both structural and functional abnormalities between brain regions within a number of resting‐state networks (Coloigner et al., 2019; Kaiser, Andrews‐Hanna, Wager, & Pizzagalli, 2015). Functional connectivity analyses use resting state fMRI (rsfMRI) to identify synchronous inter‐regional temporal correlations in blood oxygen level dependent signals (Biswal, Van Kylen, & Hyde, 1997). Common approaches to identify these functionally connected brain networks utilize independent component analysis (ICA) or Pearson‐r correlations to isolate and index the connectivity between regions in anatomically segregated brain networks (Yoo et al., 2018). Group ICA has been widely used for multi‐participant studies to identify a set of commonly replicable, temporally synchronized resting‐state networks (RSNs) (Beckmann, DeLuca, Devlin, & Smith, 2005; Damoiseaux et al., 2006; Esposito et al., 2005). Typically, these networks include the visual, somatomotor, default mode (DMN), frontoparietal (FPN), dorsal attention (DAN), ventral attention (VAN), and limbic networks (LIM) (Yeo et al., 2011). From diffusion tensor imaging (DTI), a wide range of metrics, including fractional anisotropy, mean diffusivity, tract density, tract volume and number of tracts, can be used to represent structural connectivity in MDD (Klooster et al., 2020). However, edge weight (which includes a combination of number of tracts, tract length and region of interest (ROI) size) may be a more suitable metric for the measurement of structural connectivity in relation to functional connectivity (Huang & Ding, 2016). While previous studies have reported significant alterations in structural connectivity in MDD, these findings have shown considerable variability, and depending on the network and/or tracts examined, have pointed to increases (de Kwaasteniet et al., 2013), decreases (Davis et al., 2019; Wu et al., 2020), or both increases and decreases (Wu et al., 2011) in connectivity. Supplementing the structural connectivity data with functional indices may provide some clarity regarding the brain changes that are having the greatest impact in depression.

Importantly, the above studies do point to a concordance between the functional brain regions that are dysregulated in MDD and their associated internodal structural connectivity (Greicius, Supekar, Menon, & Dougherty, 2009). However, the indices of functional connectivity do not map directly, one to one, with the white matter connectivity alterations identified by DTI (Greicius et al., 2009). Indeed, it has been suggested that neither functional nor structural imaging modalities are reliable enough alone to reflect the highly interconnected nature of the brain (Kambeitz et al., 2017; Park & Friston, 2013). Functional connectivity may arise from indirect white‐matter pathways (Nth‐order structural connections) or undetectable white matter connections (from DTI imaging techniques) (Koch, Norris, & Hund‐Georgiadis, 2002). Moreover, with DTI, branching or crossing fibers can make it difficult to resolve long‐range interhemispheric connections and therefore imposes limits on the mapping of structural connectivity (Behrens, Berg, Jbabdi, Rushworth, & Woolrich, 2007; Peled, Friman, Jolesz, & Westin, 2006; Wiegell, Larsson, & Wedeen, 2000). As a consequence, recently there has been considerable interest in the application of multimodal approaches that jointly examine the structural and functional integrity of parallel distributed neural circuits implicated in psychopathology (Reddi, 2017).

Some researchers have argued that multimodal fusion techniques may provide a better representation for whole brain connectivity and a better diagnostic classification between groups (Pineda‐Pardo et al., 2014). Instead of inferring structural connectivity from functional connectivity, and vice versa, their interconnected structural‐functional relationship is quantitatively measured.

There are a number of different approaches to combining structural and functional datasets including: (1) joint analysis (data integration) which extracts common features from separate data sources to perform statistical analysis such as correlations (Honey et al., 2009; van den Heuvel et al., 2013); (2) asymmetric data fusion which uses one dataset to constrain and analyze another (i.e., deriving structural connectivity from functional data) (Taylor & Saad, 2013); (3) symmetric data fusion (data fusion based on higher order statistics) which performs separate analyses for functional and structural data but combines them statistically (Bowman, Zhang, Derado, & Chen, 2012; Zhu et al., 2014) and (4) machine learning algorithms which utilize computational models that are automated to improve through iterative optimization (Dyrba et al., 2012; Rosa et al., 2015). Calhoun and Sui (2016) argued that among all multimodal approaches, asymmetric and symmetric data fusion, respectively, provide the most information (Calhoun & Sui, 2016). Data fusion uses statistical methods to combine the effects of different metrics (retrieved from separate complementary modalities) in a single measure. It is a more realistic representation of the real biology of brain networks, instead of studying brain networks from one angle alone (from a single modality).

To date only two studies have been published that combine functional and structural neuroimaging data in a symmetric data fusion approach for the study of MDD, both of which used a joint‐ICA fusion approach. Choi and associates conducted a preliminary study with a small sample size of four MDD participants and nine healthy comparison participants (HC) employed a joint ICA approach for combining functional and structural connectivity (Choi et al., 2008). These researchers reported changes in fractional anisotropy (FA) white matter values and changes in the strength of functional connectivity in MDD patients compared with HC (Choi et al., 2008). Furthermore, there were detectable differences in both the functional and structural connectivity in the “subgenual anterior cingulate cortex (sACC) and perigenual ACC, anterior midcingulate cortex, caudate, thalamus, medial frontal cortex, amygdala, hippocampus, insula, and lateral temporal lobe” (Choi et al., 2008). Ramezani et al. (2015), also reported a joint analysis approach with a small sample size of 25 participants (11 MDD, 14 HC). Their results indicated no detectable differences between MDD and healthy control participants when examining either fMRI or DTI in isolation, but when employing the joint‐ICA fusion approach, detectable differences in hippocampal volume loss were identified (Ramezani et al., 2015). This illustrates the added value of utilizing a combined fMRI and DTI approach for the study of MDD connectivity in a symmetric data fusion approach. Another joint‐model, which has not yet been applied to MDD data in the literature, is the “anatomically‐weighted functional connectivity” (awFC) (Bowman et al., 2012). This approach combines structural and functional connectivity in a mathematically dense approach.

Software packages for the analysis of functional and structural connectivity can substantially speed up processing time and reduce the likelihood of human error (Man et al., 2015). Presently available toolboxes designed to combine functional and structural data include: Graph Analysis Toolbox (GAT) (Hosseini, Hoeft, & Kesler, 2012), Brain Connectivity Analysis Toolbox (BrainCAT) (Marques, Soares, Alves, & Sousa, 2013), Multimodal Imaging Brain Connectivity Analysis (MIBCA) (Ribeiro, Lacerda, & Ferreira, 2015) and Brain Connectivity Toolbox (BCT) (Rubinov & Sporns, 2010). Another toolbox is the “Functional and Tractographic Analysis Toolbox” (FATCAT) (Taylor & Saad, 2013) which extracts functional connectivity (Pearson correlation) and corresponding tractography metrics (i.e., FA, tract count) between functionally‐derived ROI‐pairs. However, it does not combine these two modalities in a fusion approach. A processing toolbox consists of a set of software tools and a recommended (or definitive) pipeline. Some of these toolboxes are overly simplified (fully automated with no control over parameters), while others are highly specialized (i.e., only used for task‐based or only rsfMRI), computationally demanding (i.e., nonlinear fitting of DTI) or have complex workflows (Cusack et al., 2015).

Here, for the first time, we combine a connectivity toolbox (FATCAT) (Taylor & Saad, 2013) and a data fusion method (awFC) (Bowman et al., 2012), into a novel single pipeline known as “FATCAT‐awFC.” This yields a single powerful hybrid pipeline that combines functional and structural connectivity information into a single index, known as anatomically weighted functional connectivity (awFC). While FATCAT uses functionally derived ROIs to extract DTI parameters, awFC fuses both datasets together in a complex approach. The FATCAT‐awFC pipeline preserves the complexity of the relationship between structural and functional connectivity and provides maximal information, while allowing for simple implementation. To the best of our knowledge, this article is the first to design, explore, and compare a unique multimodal fusion approach with unimodal approaches in a large sample of patients with MDD. Using the FATCAT‐awFC pipeline, we expect to find differences between MDD patients and HC in commonly identified RSNs, including the DMN, FPN, DAN, VAN, and LIM (Yeo et al., 2011). We hypothesized that performing a joint functional‐structural connectivity analysis using the FATCAT‐awFC approach may allow us to better discriminate connectivity changes between MDD and HC groups compared to analyzing these changes using a single modality.

MATERIALS AND METHODS

This study was executed as part of the Canadian Biomarker Integration Network in Depression (CAN‐BIND‐1) (Kennedy et al., 2019; Lam et al., 2016; MacQueen et al., 2019).

Participants

The CAN‐BIND‐1 study enlisted a total of 267 participants (164 MDD and 103 Healthy Comparison participants), which were available. Of these, 17 participants were excluded due to high levels of motion in rsfMRI data (Jenkinson, Bannister, Brady, & Smith, 2002), 4 participants were excluded following a visual inspection of fMRI data for artifacts paired with a Jaccard similarity index, and 7 participants were excluded due to missing DTI data. This left a total of 239 participants, 143 MDD and 96 HC (excluded participants = 21 MDD: 7 HC). Data was collected in unmedicated MDD patients prior to the initiation of the selective serotonin reuptake inhibitor (SSRI) escitalopram. Participants were recruited from six sites across Canada: Calgary (Hotchkiss Brain Institute), Hamilton (St. Joseph's Healthcare Hamilton), and Kingston (Providence Care Mental Health Services) (Kennedy et al., 2019; Lam et al., 2016; MacQueen et al., 2019), Toronto (2 sites: University Health Network, and Centre for Addiction and Mental Health), Vancouver (Djavad Mowafaghian Centre for Brain Health). Research ethics approval for the study was obtained from the local ethics boards at each site. Study group demographic information can be found in Table 1. The demographic data was analyzed using gtsummary packages (Sjoberg, 2021) in R software (R Core Team, 2018). The Bonferroni method was applied to correct for of multiple comparisons where appropriate.

TABLE 1

Demographic and clinical characteristics of the study group

Characteristic	Healthy control participants, N = 96a	Patients with MDD, N = 143a	Group comparison p‐valueb
Sex			.9
Female n (%)	62 (64.6%)	94 (65.7%)
Male n (%)	34 (35.4%)	49 (34.3%)
Age, years mean (SD)	32 (10)	33 (12)	.8
Education, years, mean (SD)	18.5 (2)	16.9 (2)	<.001 d
MADRS mean (SD)	1 (2)	29 (6)	<.001 d
Age of onset of MDD, years, mean (SD)	NA	19 (8)
Number of MDE's mean (SD)	NA	4 (3)
Duration of current MDE n (%)
≤1 year	NA	77 (53.8%)
1–2 years	NA	14 (9.8%)
>2 years	NA	42 (29.4%)
Unknown/unreported	NA	10 (7.0%)
Antidepressants n (%)
Drug naïve	0	73 (51.0%)
Past history of antidepressants	0	70 (49.0%)
Comorbiditiesc n (%)
Social anxiety disorder	0	31 (21.7%)
Generalized anxiety disorder	0	32 (22.4%)
Panic disorder	0	23 (16.1%)
Agoraphobia	0	14 (9.8%)
Posttraumatic stress disorder	0	10 (7.0%)
Bulimia nervosa	0	3 (2.1%)
Alcohol abuse (past 12 months)	0	2 (1.4%)
Non‐alcohol substance abuse (past 12 months)	0	2 (1.4%)

Abbreviations: MDE, major depressive episode; MADRS, Montgomery Åsberg Depression Rating Scale; NA, not applicable; n, the number of participants.

n (%); mean (SD).

Pearson's Chi‐squared test; Wilcoxon rank sum test; Fisher's exact test.

The Mini‐International Neuropsychiatric Interview was used to diagnose the DSM‐IV‐TR disorders (Diagnostic and Statistical Manual).

Significant after Bonferroni correction.

Inclusion and exclusion criteria

Inclusion criteria for HC included: 18–60 years of age with no history of psychiatric disorder or unstable medical diagnosis, and able to understand instructions in English. The inclusion criteria for MDD were: 18–60 years of age, currently experiencing a major depressive episode with a duration of three or more months as defined in the Diagnostic and Statistical Manual IV‐TR (American Psychiatric Association, 2000) and as identified by the Mini International Neuropsychiatric Interview (MINI; Sheehan et al., 1998), with a Montgomery‐Åsberg Depression Rating Scale (MADRS) (Montgomery & Åsberg, 1979) score equal to or greater than 24 and sufficient fluency in English to complete study procedures. In addition, MDD participants were required to have been free of psychotropic medications for at least 5 half‐lives prior to baseline testing, and able to comprehend instructions in English.

The exclusion criteria for MDD patients excluded patients with any axis I (aside from MDD) diagnosis which would be considered a primary disorder that could interfere with the study, or bipolar I/II disorder. Additional exclusion criteria included high suicide risk or heightened risk of hypomanic switch, and previous failure to respond to more than four pharmacological interventions. Participants were also excluded if they previously failed to respond to aripiprazole or escitalopram treatments, and/or received psychological treatment within the past 3 months from baseline and planned to continue psychological treatment.

Exclusion criteria common to both groups (HC and MDD) involved individuals with: a history of substance abuse within the past six months, neurological disorders, head trauma, pregnant or breastfeeding, and/or have any other contraindications to MRI. Every participant in the study provided informed written consent and was compensated for their participation. MDD patient comorbidities are listed in Table 1.

Data acquisition

Cognitive Testing: A computerized cognitive test battery, the CNS‐Vital Signs (CNS‐VS) was used to assess participants’ level of cognitive functioning (Gualtieri & Johnson, 2006). Five cognitive subscales of the CNS‐VS were administered: memory, cognitive flexibility, complex attention, processing speed and neurocognitive index (a summary score that consists of the mean of five cognitive variables: complex attention, memory, psychomotor speed, reaction time, and cognitive flexibility) (Iverson et al., 2009).

Images were acquired using receiver head coils on six 3T MR scanners: (One Signa HDxt from GE Healthcare, USA; Three Discovery MR750 from GE Healthcare, USA; One Intera from Phillips, Netherlands; One Trio Tim from Siemens, Germany) [see (MacQueen et al., 2019) for more details].

Functional images were acquired using a whole‐brain T2*‐sensitive blood oxygen level dependent (BOLD) echo planar imaging series with the following parameters: repetition time (TR)/echo time (TE)/flip angle = 2000 ms/30 ms/75°, voxel size = 4 × 4 × 4 mm3, field of view (FOV) = 256 mm for all sites, matrix size = 64 × 64 and, slices = 34–40 for full brain coverage. During rsfMRI acquisition, participants were required to lie still, and keep their gaze on a fixation cross for a scanning time of 10 minutes, with 300 volumes recorded in total.

Anatomical reference scans were acquired across sites following a similar acquisition protocol, although Siemens scanners reported different repetition times from their MPRAGE sequence. The parameters were visually optimized to produce similar image contrast levels across sites. The 3D T1‐weighted scans were acquired using a whole‐brain magnetization‐prepared gradient echo sequence with the following parameters: TR/TE/flip angle: 6.4–7.5 ms/2.7–3.5 ms/8–15° (Exception: Siemens Scanners TR = 1760, 1840 ms), inversion time: 450–950 ms, voxel size: 1 × 1 × 1 mm3, matrix dimensions 240 × 240 and 256 × 256, slice thickness: 1 mm, number of slices: 155–192. A vitamin E pill was taped on the right side of the participants' heads as a fiducial marker. Acquisition time for anatomical data ranged from 3:30 to 9:53 minutes [see (Lam et al., 2016) and (MacQueen et al., 2019) for more details].

The DTI acquisition used single‐shot spin‐echo echo‐planar imaging (EPI) with the following parameters: TR/TE/flip angle: 8000–9000 ms/94 ms/90°, voxel size: 2.5 × 2.5 × 2.5 mm3, FOV = 240 × 240 mm, matrix dimensions 96 × 96, 155–192 slices at 2.5 mm thickness, axial slices = 52–58. A single b‐value (b = 1000 s/mm2) was applied to 30–31 noncollinear gradient directions. Image space reconstruction was completed with an acceleration factor of 2 at individual sites including: GE ASSET, Philips SENSE, or the GRAPPA k‐space method. Acquisition time for DTI data was approximately 5 min [see (Lam et al., 2016) and (MacQueen et al., 2019) for more details].

Data preprocessing

To begin, dicom images were converted to nifti using MRIcron (Rorden & Brett, 2000). An optimized preprocessing pipeline, OPPNI (Churchill, Raamana, Spring, & Strother, 2017; Churchill, Spring, Afshin‐Pour, Dong, & Strother, 2015), was used to perform the following resting‐state preprocessing steps. Principle component analysis (PCA) was used to identify the centroid of the data and measure the Euclidean distance of each volume to the centroid of all volumes. The volume with the least amount of head displacement was chosen based on the smallest Euclidean distance to the centroid. This was considered the reference volume from which the mean distance for all other volumes was calculated. Data was motion corrected using AFNI's 3dVolreg (Cox, 1996), by matching each volume displacement to the reference volume. Basic Censoring (CENSOR—from the OPPNI pipeline) was applied, to remove outlier volumes and replace them with interpolated values from neighboring volumes (Churchill & Strother, 2013). Slice timing offsets were corrected (TIMECOR—from the OPPNI pipeline) for the interleaved acquisition by using AFNI's 3dTshift (Cox, 1996) using Fourier interpolation. This was followed by AFNIs 3dBlurtoFWHM (Cox, 1996), which was used to spatially smooth fMRI images with a 6 mm full width at half maximum smoothing kernel along three directions (x,y,z). A binary mask was created for each functional run using AFNIs 3dAutomask in which non‐brain voxels were excluded and only voxels corresponding to brain areas remained. Afterwards, a neuronal tissue mask (to exclude non‐neuronal tissues such as ventricles and sinuses) was obtained using the PHYCAA+ (first part) algorithm (Churchill & Strother, 2013). In the next few steps, nuisance regressors were calculated using linear regression. A general linear model using second order Legendre polynomial was applied to the functional data to regress low frequency fluctuations, in the range of 0.01–0.1 Hz. Next, the motion parameters (derived from 3dvolreg) were used as motion parameter estimates for PCA, and were regressed from the data. PCA was performed on the six motion parameter estimates (pitch, yaw, rol l, x, y, z), whereby the PC with the largest variance (accounting for 85% variance) was taken to be the first PC (of the six PCs) and regressed out. PCA was able to correct for the maximum variance caused by head motion while simultaneously reducing multicollinearity and preserving power. A global signal removal step was performed that regressed out the first PC (highly correlated with global signal effects) of the fMRI data (Carbonell, Bellec, & Shmuel, 2011). Physiological (i.e., cardiac and/or respiratory) noise was corrected using the PHYLUS, PHYCAA+ (second part) algorithm. A low‐pass filter was then applied (LOWPASS) to remove BOLD frequencies above 0.10 Hz. FSL's FLIRT tool was then applied: first, functional data was aligned to the structural image in native space, second, functional data was transformed to register the structural image to the Montreal Neurological Institute (MNI) template (4 mm resolution). The first 5 functional volumes were discarded to avoid relaxation effects at scan start. The remaining 295 consecutive volumes were used for data analysis.

Motion artifacts (i.e., physiological motion causing ghosting), inhomogeneity (signal intensity changes and image distortions), digital imaging artifacts (i.e., phase wrap‐around artifacts) and hardware related artifacts (radio frequency inferences and spike noise) are confounding factors that affect connectivity (Maknojia, Churchill, Schweizer, & Graham, 2019; McRobbie, Moore, Graves, & Prince, 2017). However, head motion is the most problematic confounding factor that can significantly impact resting state functional connectivity (rsFC), as each voxel relies on the spatial correspondence over a time course. Sub‐millimeter motion may distort functional connectivity estimates (Maknojia et al., 2019). Motion‐related artifacts are also strongly correlated with framewise displacement (FD) measures (Dosenbach et al., 2017). Censoring the data was achieved based on the Jenkinson mean framewise displacement criteria (FDjenk). Volumes were marked as motion contaminated if FDjenk > 0.20 mm. If 125 volumes of data (~5 min or more) were retained, participants were not excluded, otherwise the participant was removed from the sample for not having enough data for the stable estimation of rsFC (Lanka & Deshpande, 2019). Thus, seventeen participants were removed due to gross motion. The Jenkinson volume‐based FD formula was calculated as follows (Jenkinson, 1999; Jenkinson et al., 2002):where R is the radius of the head (R = 80 mm), c represents the coordinates for the center of the volume, and A and b are defined as:

In addition, the correspondence between the functional data transformed to MNI space and the MNI 152 template was calculated using the Jaccard similarity index to evaluate the accuracy of registration.

When acquiring DTI data, rapid switching of applied diffusion gradients can introduce eddy currents, which warp the DTI image in the phase encoding direction (Hecke, Emsell, & Sunaert, 2016). Each participants diffusion‐weighted volumes was aligned to the b = 0 images using an affine transformation (eddy_correct) (Hecke et al., 2016; Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012) to minimize distortion by eddy currents, reduce head motion effects, and improve the signal to noise ratio (SNR). The diffusion tensor model was fit with the weighted least‐squares technique to minimize the influence of outlier volumes. In addition, the brain tissues types (gray matter, white matter and cerebrospinal fluid) were extracted using the FSL (Jenkinson et al., 2012) Brain Extraction Tool (BET) to improve co‐registration. Afterwards, the diffusion data was registered to the standard space FA atlas (1 × 1 × 1 mm resolution; average of 58 FA images) FMRIB58 (Webster, 2012).

Modifications in the awFC approach

Building on the awFC approach proposed by Bowman et al. (2012), our current study utilized independent component analysis (ICA) and FATCAT to extract networks of regions of interest, rather than performing cluster analysis as outlined in Bowman et al. (2012). The singular value decomposition (SVD) clustering process implemented by Bowman et al. (2012) is computationally expensive for a matrix of size n × m and becomes increasingly more complex between each region pair as the number of ROIs increases (Vasudevan & Ramakrishna, 2019). ICA, on the other hand, reveals distinct spatial maps, across healthy and clinical study populations (Juneja, Rana, & Agrawal, 2016; Vergun et al., 2016). ICA is a powerful methodology, and is straightforward to apply with FATCAT's recommended pipeline involving FSL's MELODIC (Griffanti, 2019; Nascimento et al., 2017).

Generated resting‐state networks (RSNs)

In our study, we used an ICA (data‐driven) approach to identify RSNs, which were then thresholded to generate ROIs. It has been suggested that data‐driven approaches are more accurate and more sensitive at detecting the greatest effects between groups (Ma, Wang, Chen, & Xiong, 2007; van de Ven, Formisano, Prvulovic, Roeder, & Linden, 2004). Group ICA (gICA) was used to derive standard RSNs from 239 participants (143 MDD, 96 HC). The rsfMRI data of MDD and HC groups was concatenated in time for each session across participants into a single 4D dataset and decomposed into 20 independent component (IC) maps using the MELODIC gICA in FSL (Smith et al., 2004). Twenty ICs is the typical dimensionality in rsfMRI studies (Taylor & Saad, 2013). Matching is performed based on spatial correlation; to match ICs to the Yeo 7‐network template (Yeo et al., 2011). This was performed using FATCAT's 3dMatch tool (Taylor & Saad, 2013), along with visual inspection. Binarized maps were created for the selected ICs that best matched the standard functional RSNs template. 3dMatch identified and extracted a total of five ICs that matched the five of the standard functional networks, which included the DMN, FPN, LIM, VAN, and DAN. However, through visual inspection, we were able to identify that gICA split the LIM into two distinct components. Consequently, to better match the LIM, we combined the ROIs of both of these components into one network by applying the fslmaths function. These combined ROIs were then defined as the LIM (see Figure 1). The remaining 14 ICA components were not included in the study because they contained non‐gray matter regions, motion artifacts, edge alignment artifacts and other networks that were not of interest.

FIGURE 1

The resting state networks and corresponding regions of interest (ROIs) derived through group independent component analyses of RS fMRI data. CANBIND‐1 resting‐state fMRI data was used to extract ROIs. Five resting‐state networks were identified and extracted from the components (DMN, default mode network; DAN, dorsal attention network; FPN, fronto‐parietal network; Limbic, limbic network; VAN, ventral attention network). Z‐score maps were thresholded and binarized using FATCATs 3dROIMaker to generate network masks (DMN, Z = 5.5; FPN, Z = 9; Limbic, Z =6; DAN, Z = 5.5; VAN, Z = 11). The colored regions depicted represent different ROIs within each network

ROI selection

The binarized spatial maps (derived from ICs) that were identified as RSNs, were stacked in a 4D image file. The 4D stacked image file was then separated out so that each network had its own 4D file (using fslsplit). This step was necessary to set an appropriate threshold for each individual network. Functional connectivity should not be set to the same unified threshold across all networks since this would inaccurately define network ROIs, as each network may vary significantly in noise level (Wang, Adeli, Wang, Shi, & Suk, 2016). Therefore, a different threshold was applied for each component: DMN (Z > 6), FPN (Z > 9), DAN (Z > 5.5), VAN (Z > 11), LIM (Z > 6), with a minimum ROI size restriction of 30 voxels (see Table 2). The levels of thresholding were selected to qualitatively and visually capture the networks observed in Yeo et al. (2011) and in line with commonly identified RSNs in the literature (Sala‐Llonch, Bartrés‐Faz, & Junqué, 2015). The FATCAT tool 3dNetCorr (Taylor & Saad, 2013) was used to generate a functional connectivity matrix using Pearson's correlation for each participant's five RSNs.

TABLE 2

Z‐threshold values for 3dROIMaker and number of ROIs per network

Network	Threshold, Z	Number of ROIs
DMN	6	5
FPN	9	7
VAN	11	5
DAN	5.5	4
LIM	6	3

Abbreviations: DAN, dorsal attention network; DMN, default mode network; FPN, front parietal network; LIM, limbic network; ROIs, region of interest; VAN, ventral attention network.

DTI image processing

Raw DTI dicom images from the scanner were converted into a single 4‐D nifti file using dcm2nii. Tensors were estimated from diffusion data using AFNI's 3dDWItoDT (Taylor & Saad, 2013) using nonlinear fits and a scheme file containing both the b‐value and b‐vectors. The following indices were estimated from the diffusion tensor: Eigenvalues (L1, L2, L3), eigenvectors (V1, V2, V3), FA, mean diffusivity (MD), axial diffusivity (AD) and radial diffusivity (average of two radial eigenvalues; RD), all of which were done in the participant's native space. All parameter estimates have some noise and errors included in their values. Thus, an advantage of probabilistic tractography is its ability to incorporate confidence intervals and uncertainty parameters into the calculation. Uncertainty of the diffusion tensor parameter was calculated using Monte‐Carlo simulation with nonparametric resampling (i.e., bootstrap and direct variants). The variance of the FA and the primary eigenvector (e1) was estimated with FATCATs 3dDWUncert (Taylor & Saad, 2013) using 300 jackknife‐resampling iterations. Together, the DTI parameters and uncertainty measures with target ROIs were used to perform probabilistic tractography.

The 3dROIMaker step (outlined above) (Taylor & Saad, 2013), also returned inflated ROIs for use with the DTI data. The inflated ROIs were necessary to allow regions to maintain contact with the mean FA tract skeleton (defined as FA > 0.2) (Nugent et al., 2019). The ROIs were transformed to each individuals' DTI space, from MNI standard space, using nonlinear transformations (Bowman et al., 2012; Yeatman, Dougherty, Myall, Wandell, & Feldman, 2012). 3dTrackID (Taylor & Saad, 2013) was applied to perform probabilistic tracking between each region pair with the following settings: tract length > 20 mm; turning angle < 60°; Nseed = 5 tract seeds per voxel; Nmc = 1000 Monte Carlo iterations; and a fractional threshold (ftr) = 0.05 (to calculate the number of tracts per voxel, included in the final white matter (WM) ROI: ftr × Nseed × Nmc = 250 tracts/voxel). An FA threshold of 0.2 was set to reduce partial volume effects after warping (Yeatman et al., 2012). The 3dTrackID (Taylor & Saad, 2013) step returned DTI metrics including: white matter volume (physical volume and number of voxels), counts of tracts, FA, MD, L1, RD and AD. Number of tracts and tract length are used later along the FATCAT‐awFC pipeline.

Generating structural connectivity matrix

A Poisson‐regression based adjustment was applied (to reduce the likelihood of false positives due to distance bias): , where is the distance between each region pair, is the unbiased number of tracts (Bowman et al., 2012). We estimated and adjusted for the bias that exists between the number of tracts and physical distance with the effect to more accurately represent structural connectivity strength. To account for indirect structural connectivity, we relied on the awFC approach to calculate all possible second‐order connections (indirect connections) with the following equation: π  = max[π , max(π π )], where π is the probabilities of structural connectivity, i is the starting ROI, j is target ROI, and m is the third connection. This equation calculated the structural connectivity probabilities of direct connections and indirect connections, taking the higher connectivity value to be the neural pathway. For instance, if the connectivity is such that the structural connectivity is higher for indirect connections versus direct connections, we took the indirect connection to be the pathway used to connect the functional regions (see Supplemental for more information).

Functional and structural connectivity combined into one unit (awFC)

Once structural connection probabilities, distance bias, direct/indirect structural pathways were calculated and factored into structural connectivity using the awFC approach, the structural connectivity was ready to be integrated into the functional connectivity (Figure 2). In this study, a multiplicative combination technique was used to derive the fused model, whereby the dissimilarity between region pairs (for connectivity) were multiplied together to generate the fused dissimilarity matrix (Liu, Li, Xu, & Natarajan, 2018). We first computed functional dissimilarity (1 − functional connectivity) and the structural dissimilarity (1 − structural connectivity) between each region pair. To transform the fused dissimilarity metric back to a correlational value (awFC), we performed a simple subtraction: . The dissimilarity metric (which emphasized correlations and differences) was transformed back to a connectivity metric; the awFC (Bowman et al., 2012) (see Supplemental for more information).

FIGURE 2

The FATCAT‐awFC analysis pipeline. awFC, anatomically weighted functional connectivity; FATCAT, functional and tractographic connectivity analysis toolbox

Correction for multiple comparisons

To test the hypothesis that the awFC underlying the RSNs varies between groups, we performed a Mann–Whitney test between each region pair within each RSN and within‐group contrasts between MDD and HC groups. Unless otherwise noted, all reported p‐values for the statistical tests of functional connectivity, structural connectivity, and awFC were corrected for multiple comparisons using the false discovery rate (FDR) criterion proposed by Benjamini and Hochberg (Waite & Campbell, 2006). The significance level was set to p (FDR corrected) < .05. Effect sizes were generated using Cohen's d in the statistics package R (R Core Team, 2018).

Statistical analysis

To evaluate the association between awFC within RSNs and cognitive changes in the MDD group compared to HC, a post‐hoc test was performed on regions with significant associations. First, multicollinearity was assessed using Pearson correlation pair plots, among the cognitive variables (neurocognitive index, memory domain, complex attention, cognitive flexibility, processing speed). Pearson correlation pairwise comparisons were produced using the function ggpairs from the GGally package (Schloerke et al., 2018) in R. Since multicollinearity exists among variables (see Figure S1a), PCA can be applied to reduce information redundancy and preserve important information (Kassambara, 2017; Refaat, 2010). PCA was performed with the R package (R Core Team, 2018) using the princomp function, in which a set of orthogonal PCs were produced corresponding to a linear combination of the original variables (Hair, Black, Babin, & Anderson, 2009). PCs were retained based on two criteria: if they had an eigenvalues >1.0 (Kaiser, 1960) and visually from the “first elbow” of the scree plot. A scree plot was created using the fviz_eig function from the factoextra package (Kassambara, 2017) in R. The PCs that met this criterion were taken to be independent variables in a PC regression model with awFC as the dependent variable. PC regression was performed using lme function in R's nlme package (R Core Team, 2018). PC regression (PCR) was used to evaluate any potential PC effects. If the PC showed a significant effect, it was evaluated further to interpret the results in terms of the original cognitive variables. This was done in order to interpret the data in a more meaningful manner. Only factor loadings greater than 0.40 were considered. Multiple regression analysis was performed to investigate whether cognitive variables interacted with the connectivity pattern within RSNs. They were conducted using the lme function in R's lme package R (R Core Team, 2018), whereby the awFC was taken as the dependent variable and each loadings > 0.40 as the independent variable. The regression was evaluated using the participant within site as a random effect. Age and sex were added as covariates in each statistical model. All the index scores were standard scores, which were mean centered. Multiple comparisons were corrected for using FDR.

RESULTS

Regions of interest

See Table 2 for ROI thresholds selected and resultant number of ROIs. Table 3 presents information on ROIs MNI coordinates of peak voxel and size of ROIs per network.

TABLE 3

Regions of interest (ROIs) defined within the five resting state networks

		MNI coordinates			Volume (# of voxels)
ROI #	Anatomical location	x	y	z
Default mode network
1	Cerebellum/lateral occipital cortex	34	−86	−40	726
2	Posterior cingulate cortex (PCC)	−2	−54	24	559
3	Medial prefrontal lobe	−6	42	12	1200
4	Middle temporal gyrus	−58	−22	−20	168
5	Left inferior parietal lobe	−50	−70	28	278
Frontoparietal network
6	Middle and inferior temporal gyrus	62	−42	−16	182
7	Right parietal lobe and lateral occipital cortex	50	−58	44	703
8	Right frontal lobe	−42	46	−4	43
9	Left frontal lobe	46	22	36	1200
10	Left parietal lobe	−50	−54	44	145
11	Cerebellum	−46	−70	−36	394
12	Frontal pole	46	14	20	127
Limbic network
13	Left cerebral cortex and temporal lobe	−30	−6	−36	400
14	Right cerebral cortex/right temporal lobe	30	−10	−36	400
15	Frontal pole, frontal medial pole	30	42	−12	850
Ventral attention network
16	Left temporal lobe	−38	−18	−8	1169
17	Cingulate gyrus	14	−34	40	962
18	Left DLPFC (frontal lobe)	−30	34	24	215
19	Right DLPFC (frontal lobe)	46	42	0	241
20	Right temporal lobe	46	−10	−16	1293
Dorsal attention network
21	Right superior parietal lobule	50	−30	40	1339
22	Right lateral occipital cortex, inferior temporal gyrus	54	−62	−12	100
23	Left lateral occipital cortex, inferior temporal gyrus	−50	−66	−8	112
24	Left superior parietal lobule	−50	−30	36	1454

Note: All ROIs were derived from FATCAT's “3dROIMaker” command. MNI coordinates and volume of each individual ROI were identified. The functionally defined ROIs covered a number of anatomical structures that were reported.

Abbreviations: ROI, region of interest; RSNs, resting state networks.

Functional connectivity group differences

Compared to HC, the MDD group showed less functional connectivity in the DMN between the posterior cingulate cortex (ROI 1) and cerebellum/occipital regions (ROI 2), (W = 8406.5, p adj = .0332). In addition, reduced functional connectivity in MDD compared to HC within the VAN was found between the left temporal lobe (ROI 16) and the right dorsolateral prefrontal cortex (DLPFC) (ROI 19) (W = 8365, p adj = .0211). We did not find significant functional connectivity differences between MDD and HC groups within the LIM, FPN, or the DAN (see Table 4).

TABLE 4

A Wilcoxon test was used to identify significant connectivity (structural, functional and anatomically weighted functional connectivity) differences between MDD and HC participants between regions of interest (ROIs).

Start ROI	End ROI	SC p‐value (FDR corrected)	FC p‐value (FDR corrected)	AwFC p‐value (FDR corrected)
Default mode network
1	2	.009**	.033*	.032*
1	5	.71	.19	.20
1	3	.71	.24	.24
2	4	.43	.23	.23
3	5	.73	.20	.20
5	2	.013*	.30	.30
Frontoparietal network
6	7	<.001***	.82	.58
8	9	.040*	.47	.47
Limbic network
13	14	<.001***	.76	.76
Ventral attention network
16	17	<.001***	.72	.71
17	18	<.001***	.87	.87
18	16	<.001***	.17	.16
19	17	<.001***	.73	.73
19	16	.64	.021*	.042*
20	19	<.001***	.072	.036*
Dorsal attention network
21	22	<.001**	.98	.90
23	24	<.001***	.76	.75
24	21	<.001***	.87	.86

Abbreviations: HC, healthy comparison; MDD, major depressive disorder; SC, structural connectivity; FC, functional connectivity; awFC, anatomically weighted functional connectivity.

p‐value (FDR corrected) < .05.

p‐value (FDR corrected) < .01.

p‐value (FDR corrected) < .001.

Structural connectivity group differences

Between‐group comparisons of structural connectivity between HC and MDD patients, are presented in Table 4. There were significant differences in all 5 RSNs for MDD compared with HC groups (see Table 4). These were characterized by lower connectivity values for MDD compared to HC for all networks. See Table 4, for a comparative summary of SC, FC, and awFC.

Anatomically weighted functional connectivity group differences

Exploratory analysis of the five RSNs awFC connectivity revealed reduced correlation differences between ROI pairs for MDD groups compared with HC in three regions: one in the DMN and two region pairs in the VAN (see Table 5). Figure 3a,b illustrates the group ROIs that demonstrated significant connectivity differences between MDD and HC groups. Lower awFC connectivity was found in the DMN between the PCC (ROI 1) and cerebellum/occipital regions (ROI 2), (W = 7917, p adj = .0322) for the MDD group compared with the HC (see Table 5). In addition, MDD patients demonstrated lower awFC in the VAN, between the left temporal lobe (ROI 16) and the right DLPFC (ROI 19) (W = 8274, p adj = .0421) compared with HC (see Table 5). Reduced connectivity in MDD compared with HC within the VAN was observed between the right temporal lobe (ROI 20) and the right DLPFC (ROI 19) (W = 8366, p adj = .0361). No other significant differences were found within the remaining RSNs. A summary of the significant regions and p‐values is presented in Table 5, and a summary of the mean and standard deviations of connectivity values are presented in Table 5. In addition, Figure 4 displays boxplots of awFC values between ROI‐pairs within RSNs.

TABLE 5

A Wilcoxon test was performed whereby significant anatomically weighted functional connectivity differences between MDD and HC groups, and their associated structural and functional connectivity (significance set to FDR‐corrected p < .05) are reported. Mean, standard error and effect size are also reported for each group.

						MDD		HC
Start ROI	End ROI	SC padj	FC padj	awFC padj	Effect size	Mean	SE	Mean	SE
Default mode network
1	2	.00905**	0.0331*	0.0320*	0.34	−0.011	0.017	0.061	0.022
Ventral attention network
16	19	0.64	0.021*	0.0421*	0.32	0.36	0.024	0.46	0.027
19	20	<.001***	0.0722	0.0361*	0.36	0.22	0.023	0.33	0.028

Abbreviations: ROI, region of interest; SC, structural connectivity; FC, functional connectivity; awFC, anatomically weighted functional connectivity; padj, FDR corrected p‐values; MDD, major depressive disorder, HC, healthy comparison; SE, standard error.

p‐value (FDR corrected) <.05.

p‐value (FDR corrected) <.01.

p‐value (FDR corrected) <.001.

FIGURE 3

ROIs within RSNs with significant brain connectivity group differences (a) VAN, orange = right DLPFC, blue = left temporal lobe, green = right temporal lobe. (b) DMN, red ROI = cerebellum/lateral occipital cortex, yellow ROI = posterior cingulate cortex (PCC). ROIs, regions of interest; RSNs, resting state networks; A, anterior; P, posterior; R, right; L, left; S, superior; I, inferior

FIGURE 4

Boxplots demonstrated lower anatomically weighted functional connectivity between ROI‐pairs for the major depressive disorder (MDD) and healthy comparison (HC) participants. Boxplots also quantified the strength of connectivity for MDD and HC groups (a) AwFC between the left temporal lobe and the right DLPFC within the VAN (b) AwFC between the right temporal lobe and the right DLPFC within the VAN (c) AwFC between the occiptial lobe/cerebellum and the PCC within the DMN. AwFC, anatomically weighted functional connectivity; DLPFC, dorsolateral prefrontal cortex; VAN, ventral attention network; DMN, default mode network; Occ., occipital; Cerr., cerebellum; PCC, posterior cingulate cortex. Asterisks identify significant between‐group differences following FDR correction (p < .05)

Post hoc analysis

Performing PCA on the cognitive variables resulted in five PCs (each PC is a linear combination of the original variables). Applying the Kaiser‐Guttman rule (Guttman, 1954; Kaiser, 1961), of extracting only PCs with an eigenvalue >1, revealed the first PC had an eigenvalue >1 (). In addition, from the scree plot (Figure S1b), we selected one PC at the marked “elbow,” whereby 62.5% of the variance in the data was explained. Therefore, based on these two criteria, only the first PC was retained. PCR was then performed on the first PC, in which a significant effect was detected. This warranted further analysis to determine the cognitive variable associated with the awFC. Variables (loadings > 0.4) were investigated further. The neurocognitive index (NCI) (loading = 0.40), cognitive flexibility (loading = 0.45), processing speed (loading = 0.49), and complex attention (loading = 0.52) met this criterion. Therefore, multiple linear regressions were performed on each variable individually. Therefore, multiple mixed‐effects linear regressions were performed, separately on each variable.

Applying PCR to the DMN (uncorrected p = .0475) showed significant associations between the first PC and the awFC. Multiple linear regressions revealed that awFC between the PCC and the cerebellum/occipital lobe was significantly associated with changes in MADRS (uncorrected p = .02) and complex attention (uncorrected p = .047). Second, a PCR for the VAN revealed significant associations between the first PC and awFC (uncorrected p = .02). Therefore, multiple linear regressions were performed for variables with loadings > 0.4. The VAN, between the right temporal lobe and the right DLPFC revealed, awFC was significantly associated with complex attention (uncorrected p = .028). However, none of the cognitive associations survived correction for multiple comparisons.

DISCUSSION

We developed a novel pipeline for combining functional connectivity (derived from fMRI) and structural connectivity (derived from DTI) and used it to study awFC connectivity changes in MDD patients. We analyzed a total of 24 ROIs (from five resting‐state functional networks), three of which revealed statistically significant differences in awFC between MDD and HC groups using the FATCAT‐awFC approach. For each region pair we also conducted standard functional and structural connectivity analyses to compare against our novel combined functional‐structural analysis approach (FATCAT‐awFC).

As predicted, the multivariate connectivity analysis was capable of revealing group differences not identified by the univariate analysis within RSNs. We found reduced awFC connectivity within the DMN between the PCC and cerebellum/lateral occipital cortex in the MDD group compared with the HC group. Aberrant connectivity between this ROI‐pair was also found in both the traditional structural and traditional functional connectivity approaches, which supported the changes we found in the awFC measure output by the FATCAT‐awFC pipeline. In MDD, similar findings in resting state functional connectivity (rsFC) have been reported by Liu et al. (2012) who found that rsFC between the PCC (associated with self‐referential processing) and cerebellum was reduced in MDD groups compared to HC. Negative self‐referential processing, is a common feature of MDD and has been associated with MDD severity (Lou, Lei, Mei, Leppänen, & Li, 2019). We also observed decreased awFC between the right DLPFC and the right temporal lobe (encompassing the right temporo‐parietal junction). In this case, when we examined the functional and structural data independently, we did not find differences in functional connectivity, but there were detectable changes in structural connectivity. Past literature (Hwang et al., 2015; Penner et al., 2018) has identified that reduced functional connectivity between the right temporo‐parietal junction and right DLPFC in MDD, is associated with depression severity. Lower connectivity between these regions was associated with more severe depression symptomology (Hwang et al., 2015). We also observed reduced awFC between the right DLPFC and the left temporal lobe (includes the left temporo‐parietal junction). It is interesting to note that we observed significant differences in functional connectivity between this region pair, but structural connectivity had no detectable changes with MDD compared to HC. With our FATCAT‐awFC approach, it appears that significant group differences in the underlying SC do not drive the awFC connectivity differences as much as FC differences do. Our finding were supported by (Hwang et al., 2015; Penner et al., 2018), who also detected decreased functional connectivity between the right DLPFC and the left temporal lobe in MDD patients. According to (Samson, Apperly, Chiavarino, & Humphreys, 2004) the reduced connectivity between the frontal lobe and the left temporo‐parietal junction may be associated with difficulty assessing thoughts that another person possesses. Inaccurately interpreting mental states, often results in reduced social interactions and may contribute to further social isolation in MDD (Weightman, Air, & Baune, 2014). Connectivity changes (between these three region pairs) in MDD reported in our study are well documented in MDD literature, confirming that our method is effective at detecting some of the neural changes associated with MDD.

In contrast to our predictions, we did not observe group differences in awFC network connections of the LIM. These findings may be a result of our use of relatively large ROIs, which encompassed a number of brain regions with varying functions. In our analyses, the LIM that was extracted by 3dMatch consisted of three large ROIs, where the first two were homotopic ROIs, whereby each ROI included the parahippocampal gyrus, temporal fusiform gyrus, inferior temporal gyrus; the third ROI crossed the midline and consisted of the bilateral amygdala, nucleus accumbens, caudate, paracingulate gyrus, frontal medial gyrus, and putamen. Consequently, it may be that because these ROIs were so large, they lacked the specificity necessary to identify localized group differences. In addition, in contrast to expected patterns, we did not observe group differences in awFC network connections in either the FPN, or DAN.

In our study, we conducted a post hoc test that consisted of a PCA and PCR, to study the association between connectivity and cognitive data. We found a trend in the association between complex attention and the awFC within the DMN (PCC to cerebellum/occipital regions). We also found a trend level association between complex attention and the awFC between the right temporal lobe to the right DLPFC. The standard functional Yeo template, classifies the anterior temporo‐parietal junction as a component of the VAN (Yeo et al., 2011). The VAN is involved in regulating emotional salient events (Korgaonkar, Goldstein‐Piekarski, Fornito, & Williams, 2020) and as such. our observation of reduced awFC between this ROI‐pair may point to connectivty changes that underlie the mood dysregulation associated with MDD (American Psychiatric Association, 2013).

Our FATCAT‐awFC approach identified region pairs that were observable in one modality but not the other. The FATCAT‐awFC model identified group differences in connectivity, some of which were only captured using SC, and others identified through FC. Using SC alone resulted in many region pairs having detectable differences in MDD compared to HC. However, this may be a result of false positive findings in SC (Bowman et al., 2012), whereas FC identified just two regions with group differences after correction for multiple comparisons. Furthermore, our FATCAT‐awFC approach revealed a connectivity change in MDD that was undetected using the conventional FC approach alone. Finally, for most comparisons the FATCAT‐awFC approach resulted in p‐values that were lower than those in a single modality, suggesting that FATCAT‐awFC was sensitive to connectivity changes distinguishing groups.

Contributions of the FATCAT‐awFC

The FATCAT‐awFC pipeline was designed to be a more practical solution for combining FC and SC together. With our FATCAT‐awFC pipeline, we hope that researchers will benefit from a faster and more intuitive approach to combining SC and FC, as opposed to using a computationally intense method (awFC method) or a toolbox that does not provide maximal information (FATCAT) alone. The combination of FATCAT with awFC provides for a unique hybrid pipeline that combines the advantages of an intuitive, rapid and efficient toolbox with a computationally intense data‐fusion approach that provides an abundance of information. By combining FC and SC data we are able to better represent whole brain connectivity as opposed to studying it utilizing data from one modality alone. The awFC metric is able to measure the combined effect of SC and FC and may provide us with a more accurate connectivity value as it relates to different neuropsychiatric disorders.

There were two main limitations in this study. The first is that one set of group ROIs (derived from gICA) were generated from the CAN‐BIND participants. This was done to have a consistent number of ROIs across participants; however applying the same group networks across all participant‐level networks implies spatial similarity among all participants (Sohn et al., 2015). However, gICA does not account for inter‐participant variability in functional boundaries, and does not construct participant‐specific spatial maps of networks. The second limitation was the ambiguity of selecting a cut‐off threshold for ROI creation. There is no standard for selecting ROIs, although many studies have selected a threshold of Z > 2.3 (Sohn et al., 2015). ROI selection is a threshold‐dependent process and can have very different effects on the outcomes and conclusions of a study (Sohn et al., 2015). Our ROIs were fairly large and encompassed a number of functionally different regions. Tong and associates (Tong et al., 2016) found that larger ROIs are often accompanied by greater variance within connectivity data, in comparison to smaller ROIs, resulting in a smaller effect size. Consequently, the larger ROI used in this study may have reduced our capacity to detect group differences.

Hebbian theory is summarized as: “neurons that fire together, wire together” (Keysers & Perrett, 2004), which suggests that regions that are functionally connected (temporally synchronous) are in principal structurally connected through fiber tracts. The complimentary nature of FC and SC suggests their combination would provide a more complete picture of connectivity as a whole. Due to the complex nature of data fusion in neuroimaging, many researchers have not utilized this approach. By integrating a toolbox (FATCAT) with a computationally intense technique for combining SC and FC (awFC), we were able to implement a relatively straightforward approach to the combination of functional and structural data. We also showed that FATCAT‐awFC was capable of identifying ROI‐pairs, which would have been missed when only applying unimodal analyses. Hence, the joint approach of FATCAT‐awFC allows for a more detailed understanding of the interconnected nature between structural and functional connectivity and how it relates to depression.

CONFLICT OF INTEREST

Dr. Milev has received consulting and speaking honoraria from Allergan, Janssen, KYE, Lundbeck, Otsuka, Pfizer and Sunovion, and research grants from CAN‐BIND, CIHR, Janssen, Lallemand, Lundbeck, Nubiyota, OBI, OMHF and Pfizer. Dr. Frey has received a research grant from Pfizer. Dr. Strother receives funding from the OBI and CIHR (MOP137097) for neuroimaging analysis in CAN‐BIND and he is the Chief Scientific Officer of ADMdx, Inc., a neuroimaging consulting company. Dr. MacQueen has had consultant payments or honoraria from: Allergen, Pfizer, Lundbeck, Janssen, Johnson & Johnson. Dr. Kennedy has received research funding or honoraria from the following sources: Abbott, Alkermes, Allergan, BMS, Brain Canada, Canadian Institutes for Health Research (CIHR), Janssen, Lundbeck, Lundbeck Institute, Ontario Brain Institute, Ontario Research Fund (ORF), Otsuka, Pfizer, Servier, Sunovion and Xian‐Janssen. Dr. Kennedy holds stock in Field Trip Health. Dr. Lam has received honoraria for ad hoc speaking or advising/consulting, or received research funds, from: Allergan, Asia‐Pacific Economic Cooperation, BC Leading Edge Foundation, Canadian Institutes of Health Research, Canadian Network for Mood and Anxiety Treatments, Canadian Psychiatric Association, Hansoh, Healthy Minds Canada, Janssen, Lundbeck, Lundbeck Institute, MITACS, Movember Foundation, Ontario Brain Institute, Otsuka, Pfizer, St. Jude Medical, University Health Network Foundation, and VGH‐UBCH Foundation. All other authors report no biomedical financial interests or potential conflicts of interest.

Supplementary Figure S1 (a) A pair plot with all the variables (memory domain, complex attention, cognitive flexibility, processing speed and neurocognitive index) showing the Pearson's correlation coefficients demonstrates multicollinearity exists among variables. (b) Scree Plot shows the eigenvalues/variance for each of the principal components. The number of principal components that were retained in our study was determined at the first “elbow” of the plot. Note: Mem_dom, memory domain; comp_attn, complex attention; cogn_flex, cognitive flexibility; proc_spd, processing speed; nci, neurocognitive index; Corr, correlation.

Click here for additional data file.