Ideally, if scanning time was not a constraint and image acquisition was carried out indefinitely, the number of recorded time observations (N) would eventually exceed the number of voxels (p) and a simple linear estimate of the association between one voxel i and the remaining gray matter voxels in the brain could be obtained by ordinary least squares (OLS). Under this theoretical scenario, the multiple regression equation would give the predicted time series for the target voxel (Y^i)

and a simple measure of functional connectivity for this voxel would be provided by the multiple correlation coefficient

which quantifies the degree of similarity between the expected (Y^i) and observed (Y_i) time series of voxel i, and which is given by the square root of the coefficient of determination (i.e., R², a standard ouput from regression analyses).

However, in a real fMRI dataset the number of voxels is much larger than the number of time points (N≪p) and OLS is not feasible. This limitation may be overcome by drastically reducing the number of variables through dimensionality reduction or by considering ROIs from a brain parcellation. Here, though, there is also the alternative of using a regularization approach.

Specifically, regularization through ridge regression allows obtaining Y^i by setting a restriction on the parameter estimates of Equation 1

(where ct stands for constant value). Such restriction may be restated as a constrained least-squares minimization regulated by a Lagrange multiplier (λ ≥ 0)

where y_i,1, …, y_i,N stand for the individual components (time points) of Y_i.

Here, the selection of an adequate value for λ (the regularization parameter) will be important in order to achieve a good balance between bias and variance (i.e., to find a model that avoids overfitting while not being not too constrained). Furthermore, some aspects will have to be considered for the validity of ridge regression in the framework of regularized brain connectivity (RBC). Apart from the usual rescaling of variables (time series) to unit variance, λ values will have to remain fixed through all voxels and individuals, otherwise estimates of Equation 2 will not be comparable. Unfortunately, due to the lack of independence between observations in time series, the habitual cross-validation methods will not be suitable. Since there is no easy alternative to cross-validation in the current framework, we have decided to report results using a wide range of λ values (0.5, 1, 5, 10, 50, 100, and 500). All ridge regressions have been carried out with functions contained in the glmnet R library (Friedman et al., 2010).

To highlight the multivariate nature of the RBC, results obtained have been compared to those provided by simple averages of bivariate correlations. To do so, the weighted Global Brain Connectivity (GBC) method (Cole et al., 2010) has been also applied to the resting state images. Specifically, averages of absolute values of all bivariate correlations involving each specific voxel have been calculated. Results from the RBC have been also compared to connectivity estimates obtained through supervised principal component regressions (Bair et al., 2006; Hastie et al., 2009), a dimensionality reduction approach applied to quantify brain connectivity from fMRI images (Salvador et al., 2017).

A sample of N = 148 patients with a diagnosis of schizophrenia according to DSM-IV criteria (i.e. excluding patients with schizoaffective and other schizophrenia related disorders) were recruited from three hospitals from Germanes Hospitalàries located in the Province of Barcelona, Spain (Hospital Bennito Meni C.A.S.M., Hospital Sant Rafael and Hospital de la Mercè). All patients but two were taking antipsychotic medication (atypical N = 109, typical N = 9, both N = 24, unknown N = 3, equivalents of Chlorpromazine: 508.82 mg (mean), 517.55 mg (SD). A second sample of N = 74 healthy controls were recruited from non-medical hospital staff, their relatives and acquaintances, plus independent sources in the community.

Controls reporting a history of mental illness and/or treatment with psychotropic medication or with a psychotic first-degree relative were not included. All individuals in both samples were right handed, aged 18–65, with no history of brain trauma or neurological disease, and not having shown alcohol/substance abuse in the last 12 months. All subjects gave written informed consent before participation and the study procedures were approved by the Comité de Ética de la Investigación de FIDMAG Hermanas Hospitalarias and adhered to the Declaration of Helsinki.

Participants underwent a single MRI session in a 3.0 Tesla Philips Ingenia machine located in the Hospital Sant Joan de Déu (Barcelona, Spain) in which a resting state functional MRI (fMRI) sequence and a T1 structural image for anatomical reference were acquired. Parameters for the resting fMRI bold sequence were: TR = 2,000 ms, TE = 30 ms, flip angle = 70°, in-plane resolution = 3.5 mm × 3.5 mm, FOV = 238 mm × 245 mm, slice thickness = 3.5 mm, inter-slice gap = 0.75 mm, number of volumes = 256. Slices (32 per volume) were acquired with an interleaved order parallel to the AC-PC plane. The T1 image was acquired using a Fast Field Echo (FFE) with TR = 9.90 ms; TE = 4.60 ms; Flip angle = 8°; voxel size = 1 mm × 1 mm; slice thickness = 1 mm; slice number = 180; FOV = 240 mm). fMRI preprocessing steps included movement correction, spike scrubbing, regression of noise-independent components, non-linear normalization to the Montreal Neurological Institute space, regression of noise from ventricles and white matter, and low-frequency filtering in the 0.1–0.02 Hz interval (Salvador et al., 2017). Specifically, for the regression of noise-independent components, individual independent component analyses were previously run with Melodic, a module included in FSL (Smith et al., 2004) and those components showing clear noise patterns (most frequently edge effects due to movement) were selected. Time series of the selected components were regressed out from the time series of each voxel, and residuals were kept as the denoised time series.

Prior to the calculation of the RBC, GBC, and NRC maps a common gray matter mask was applied to the normalized fMRI images. Then, fMRI volumes were resampled to a 4 mm × 4 mm × 4 mm voxel in order to reduce computational costs. Finally, values from the individual RBC, GBC, and NRC maps where Fisher transformed before carrying out the group analyses.

A subsample of N = 74 patients matched for gender, age and premorbid IQ to the sample of N = 74 controls was selected from the original sample of patients before carrying out group comparisons in RBC, GBC, and NRC. Premorbid IQ was estimated with the Word Accentuation Test (Del Ser et al., 1997). Positive and Negative Syndrome Scale (PANSS) scores (Kay et al., 1987) were used to calculate values of the three Liddle Syndromes (Positive, Negative, and Disorganization) (Liddle, 1987) for all N = 148 patients, and regressions between RBC, GBC, and NRC and clinical severity as measured by the three Syndromes were performed in this larger sample. In the regression analyses, age, gender, type of antipsychotic medication (atypical and/or typical) and dose of medication (equivalents of Chlorpromazine) were considered as nuisance covariates. In all analyses statistical significance was derived from permutation tests carried out using the randomize function included in FSL (Smith et al., 2004) and applying the Threshold-Free Cluster Enhancement (TFCE) method.

When RBC maps from patients were compared to those of healthy controls only areas of significant reductions in connectivity were found. As shown in Figure 1, this pattern was consistent through all regularization values. Tests performed with λ = 50 produced the most extensive differences (Figure 2) including bilateral connectivity reductions in the supplementary motor area, paracentral lobule, precentral and postcentral gyri, precuneus, dorsal and anterior cingulate, ventromedial and right dorsolateral prefrontal cortex, cuneus, calcarine, lingual, lateral occipital areas and parts of the cerebellum.

All λ values applied led to much more widespread patterns of differences than those provided by both the GBC and NRC methods, which were mainly restricted to connectivity reductions in clusters of moderate size in the supplementary motor area / paracentral lobule and postcentral gyrus (Figure 1). As shown in Figure 2, though, reductions in the NRC were clearly larger than those observed with the GBC. A cluster of increased connectivity in patients was also observed in the right Rolandic operculum (Figure 1).

When RBC maps were correlated with scores of the three Liddle Syndromes, a decreasing pattern of connectivity was observed for the Negative Syndrome (Figure 3A). This negative association involved a cluster located in the left occipital cortex (clusters shown in Figure 3 are based on analyses using λ = 50, which led to most extensive differences in the group comparison). Another cluster of negative association with the Negative Syndrome was observed for the GBC, which was also in the occipital cortex, although with a more medial position. The GBC also showed two clusters of negative association with the Disorganization Syndrome, one of them located in the right fusiform area, and the other in the right posterior insula (Figure 3B). The later was also the site for the only significant cluster observed with the NRC, which was of small size and was also negatively correlated with the Disorganization Syndrome. No significant associations were found between any of the three connectivity measures and the Positive Syndrome.