A Dynamic Entropy Approach Reveals Reduced Functional Network Connectivity Trajectory Complexity in Schizophrenia

1. Introduction

Schizophrenia, or split mind, ranks among the most studied disorders in modern psychiatry and neuroscience. There are manifold reasons for this, but perhaps the two most prominent are the extent of its resultant social problems and the sheer variety of symptoms which it can produce. Symptoms of schizophrenia broadly fall into three categories, positive, negative, and cognitive, and can range from delusion, paranoia, and hallucination to apathy, anhedonia, and social withdrawal. This wide range of symptoms makes schizophrenia a difficult disorder to treat, as treatments for one group of symptoms may have limited or no effect on another. Indeed, it has been suggested that schizophrenia should not be considered a single disorder at all, but rather a family of disorders aggregated by historical precedent [1].

For the past quarter-century, the most prevalent hypotheses on schizophrenia’s causes have focused on the interactions between brain regions. Perhaps the most prominent is the dysconnectivity hypothesis [2], which may be summarized as the proposition that aberrant connectivity, rather than focal abnormalities, is the primary cause of schizophrenia. In this sense, it may be considered a forerunner to the concept of distributed cognition—quite an early forerunner, as it was first proposed at the beginning of the 20th century [3–5]. The past half-century of research into schizophrenia has amassed considerable evidence in its favor [6], with experimental results suggesting that schizophrenia affects connectivity and plasticity from the level of individual synapses [7–11] to interregional communication [12–14] and white matter tracts. Most of these studies, however, base their analyses on anatomically rather than functionally derived atlases and do not examine time courses through a dynamic systems-inspired latent or state space. They may, therefore, miss some of the crucial insights which examining such systems through these lenses can provide.

At the “global”, or whole-brain, scale, most studies have examined static functional connectivity (FC) or structural connectivity (SC). Such studies have reported reduced global connectivity in schizophrenia patients [15,16], particularly between auditory, sensorimotor, and visual networks [17]. In addition, network-based analyses of functional connectivity in schizophrenia indicate a general reduction in organization and efficiency in the structural and functional connectivity of schizophrenia patients [18]. Beyond this, however, findings have been inconsistent. Some reports suggest decreased communication between frontal and temporal areas of the brain [12]; others have found increased connectivity within the default-mode network [19]; still others imply decreased connectivity within and between the default-mode network and cortical regions [20]. Attempts to model the static functional connectivity of schizophrenia patients have provided equally confused results, with both reduced effective connectivity and increased structural connectivity suggested [21]. Overall, then, analysis of structural and static functional connectivity has been unable to conclusively identify the large-scale changes which underly the symptoms of schizophrenia.

One potential reason for the inconclusiveness of static analyses is that they, by nature, neglect the dynamism of real brain activity. For decades, electroencephalography (EEG) has demonstrated that functional microstates constantly and fluidly change [22–24], a fact which was captured in fMRI in 2010 [25,26]. Although fMRI captures microstate alterations at a timescale of seconds rather than the milliseconds of EEG, its superior spatial resolution compared to EEG allowed researchers to precisely identify recurring connectivity states [15,27–30]. Schizophrenia patients have displayed altered dynamics of these states. For instance, a 2014 study suggested that patients have a higher probability of entering states with attenuated cortical-subcortical connections and increased intra-sensory connectivity than controls, and that these same patients also display elevated low-frequency power in thalamo-sensory communication [17]. This same study suggested that the thalamo-sensory hyperconnectivity reported in its static FC analysis may result from these alterations. However, aside from state transition probabilities and average dwell times, this study did not report any metrics designed to capture state dynamics. This lack of metrics designed to capture the dynamics of functional connectivity has proven a problem, not just for schizophrenia research, but for the broader field of functional neuroimaging.

This study begins to bridge that gap by applying a previously developed framework of measuring the entropy rate of individual fMRI scans [31] to a dataset consisting of medicated schizophrenia patients and demographically matched controls. Rather than attempting to explicitly track interactions between brain regions, this framework defines a basis for a state space within which to plot each subject’s time course. This approach allows direct application of tools and metrics from dynamic systems analysis and information theory to subject time courses and thus reduces the challenge of comparing subject trajectories. The decision to quantify fMRI dynamics in a state space is not novel [32], but the present framework specifically seeks to define a basis, i.e. a minimal linearly independent spanning set, for this space to make most efficient use of the Shannon entropy. The results of this analysis suggest substantially disordered state dynamics in the patient population, with patients displaying simpler trajectories in five out of eight state space dimensions, one of which correlates with cognitive performance. Overall, this represents a positive step towards quantifying the dynamic connectivity alterations which underly psychiatric disease and towards the discovery of a clinically useful basis space in which to diagnose and predict patient treatment.

2. Materials and Methods

A high-level overview of the pipeline is displayed in Figure 1. Detailed explanations of each step follow.

A version of the Kozachenko-Leonenko estimator [63,65,66] is employed to estimate the entropy rate for each subject. Specifically, in each subject, this estimator computes the entropy rate of each source $s_i$ over the course of the scan. The temporal independence which ICA enforces between sources ensures that the entropy rate of source $i$ has minimal effect on the entropy rate of source $j$ $\left(i\ne j\right)$. The joint entropy rate of each subject may thus be estimated by summing the subject-wise entropy rates of all sources.

Data Collection

This study utilizes control and schizophrenia patient data from the Function Biomedical Informatics Research Network (FBIRN) repository [33] preprocessed according to the description given in [34]. To summarize, the statistical parametric mapping package (SPM12) was used to correct for subject head motion, slice timing differences, to warp subject anatomy to the Montreal Neurological Institute (MNI) echo planar imaging (EPI) template space, to resample the collected data to $3\times 3\times 3$ mm³ isotropic voxels, and to smooth the resampled fMRI images with a Gaussian kernel with a full width at half maximum (FWHM) of $6$ mm. Subjects with head motion greater than $3°$ were excluded from study, as were subjects whose full brains could not be normalized due to incomplete imaging data. These criteria led to a final dataset of 151 schizophrenia (SZ) patients and 160 healthy controls (HC).

Estimation of the Functional Networks

Spatial functional networks were estimated using NeuroMark’s adaptive independent component analysis (adaptive-ICA) [34], which extends spatially constrained independent component analysis [35,36] to map known fMRI network templates to novel subject data. This requires balancing two at times competing goals: to maximize the spatial independence of networks in each subject, and to ensure that the network maps in each subject correspond to known group-level templates. Here we use the multi-objective optimized ICA with reference (MOO-ICAR) approach which maximizes two competing objective functions in turn until a solution is achieved. This allows adaptive-ICA to capture subject-unique characteristics while maintaining comparable functional networks across datasets. It should be noted that this method allows the capture of both the internal structure of brain functional connectivity networks and the extent of inter-network connectivity via static and sliding-window functional connectivity estimates.

Estimation of the Functional Networks

Before estimating the functional network connectivity (FNC), Du et al., 2020 chose to remove noise sources from each functional network’s subject-level time series. Removal of noise sources involved four steps: first, the removal of linear, quadratic, and cubic trends in the data; second, multiple regressions of the six realignment parameters and their temporal derivatives to control for in-scanner motion; third, de-spiking to remove outliers; and fourth, band-pass filtration to select for signals in the 0.01-0.15 Hz frequency bands. Once these steps were completed, subject-level static functional network connectivity (sFNC) was computed via Pearson correlation. Other measures of statistical similarity could have been used; for instance, mutual information has been proposed due to its sensitivity to nonlinear interactions [37,38]. However, Pearson correlation’s simplicity, interpretability, and ease of computation means it remains the dominant method for estimating functional connectivity.

While the static FNC provides valuable information on the extent of inter-network communication, its poor time resolution makes it unable to capture the dynamics of this communication. The two most notable methods proposed to circumvent this problem are the sliding time window approach [25,27,39] and coherence-based connectivity [40–43]. The present study uses the sliding window approach. As the name suggests, this method slides a window over the time series of each ICN in small steps, thus segmenting the total time series into many short, overlapping time series. The functional network connectivity of each time series window is computed in the same way as for static FNC, and the resulting $N\times N$ connectivity matrices concatenated into an $N\times N\times T$ array ($N$ being the number of functional networks and $T$ the number of time series windows). This study convolved a normal distribution of zero mean and standard deviation of three $\left(\mathcal{N}\left(\mu =0,\sigma =3\right)\right)$ with a rectangle 40 times-to-repetition (TRs) long [34] to generate its selection window.

Number of Temporally Independent Sources

Most clustering or source separation algorithms require users to specify the number of sources for which the algorithm should search. This poses a problem in neuroimaging analysis, as the number of recurrent connectivity states which the human brain expresses has not been fully established. As such, researchers do not know a priori the appropriate model order to input into their separation algorithm. It is common to circumvent this problem by testing several numbers of recurrent states, usually within the range of four to eight [44–48] and determining the optimal number by applying a validity index or comparing the results of runs with different source counts. In this instance, eight sources provide the best balance between model order and goodness of fit.

Temporal Functional Network Connectivity Profiles (tFNCPs)

The goal of this study is to identify recurring, overlapping, temporally invariant functional brain patterns and to capture their dynamics. Both clustering and source separation algorithms may be used to achieve this goal, dependent on how these states are hypothesized to manifest in time-resolved fMRI images. A clustering algorithm, such as k-means, may be appropriate if each time-resolved fMRI “snapshot” primarily consists of a single state $S$. If, on the other hand, each time-resolved snapshot contains a mixture of states $\mathit{S}$, a blind source separation algorithm may be more applicable. As this study presumes that each time-resolved fMRI image consists of a linear mixture of underlying “source” connectivity states $\mathit{S}$, a linear source separation algorithm is required. In addition, an entropy-based analysis is made considerably simpler when the time courses $\mathit{T}\mathit{C}$ are statistically unrelated. As such, we required a linear source separation algorithm which maximizes the statistical independence of each predicted source. Independent component analysis (ICA) is the obvious solution to this problem, as it has been proven to minimize the statistical dependencies between sources in neuroimaging data [49] and to isolate functionally meaningful communities in spike train data [50].

Most ICA algorithms accept two-dimensional inputs, with one dimension representing input variables and the second representing samples. For instance, the time series of $N$ functional networks should produce a data array of $N\times T$, with $T$ being the number of samples taken over the course of recording. A dFNC array, on the other hand, has three dimensions: $N\times N\times T$, which standard ICA algorithms cannot process. We circumvent this by converting the upper triangle of each sample dFNC$\left(t\right)$ array to a vector. Repeating this process across all samples converts the data array from $N\times N\times T$ to $\frac{N\left(N-1\right)}{2}\times T$, which can be decomposed into maximally independent sources by any standard ICA algorithm.

Before maximizing independence, researchers typically whiten the data via principal component analysis (PCA). This serves to reduce the dimensionality of the input data and to improve the estimation of independent components, both of which optimize and stabilize ICA outputs. We apply InfoMax [51–53] for the independent maximization step (ICA). This algorithm was run 150 times and the results consolidated via the ICASSO method [54] to identify stable, repeatable tFNCPs and time courses.

Entropy Analysis

We used a Kozachenko-Leonenko entropy estimator [55–58] to estimate each state’s subject-level rate of entropy production. This produces an $S\times N_{subjects}$ array of entropy rate values, with $S$ being the number of sources and $N_{subjects}$ the number of subjects. The substates’ temporal independence allows us to estimate the joint entropy of all sources by summing the entropy rate across substates [59]:

$$H\left(C_1,\dots ,C_N\right)=\sum _{j=1}^{N}H\left(C_j\right)$$

Note that this calculation requires that source signals be statistically independent. Dependencies between sources must be accounted for when estimating source entropy, which makes evaluating joint entropy quite difficult when dealing with statistically related signals. ICA’s minimization of inter- source relationships was a major motivation for its selection.

Comparison Tests

We utilized both a difference-of-means permutation test [60] and the Kolmogorov-Smirnov two-sample test to search for differences between groups. Student’s t-test was also employed when the Jarque-Bara test indicated normally distributed data. Multiple-comparison correction consists of the false-discovery rate [61].

Regression Analysis

To correct for possible confounds in the data, a multiple linear regression is employed alongside standard hypothesis tests. This regression separately models the effects of site, age, gender, and diagnosis on subject-level joint and source entropies. Additional regression analyses examine the effects of clinical (PANSS positive and negative) and cognitive (CMINDS composite) scores on subject entropy distributions, again while correcting for site, age, and gender effects. As cognitive scores are highly correlated (Figure 2), it was decided to examine only the most general score, namely the CMINDS composite score [62], to avoid effect cross-contamination.

3. Results

3.1. Static Functional Network Connectivity

As a pre-analysis sanity test, we compared static functional network connectivity (sFNC) in the same fashion as [34] for comparison purposes. Figure 3.a displays the mean sFNCs of the FBIRN dataset and 3.b the joint entropy distributions of controls and patients. Notably, healthy controls display visibly greater mean absolute FNC values than schizophrenic patients across the majority of the sFNC matrices’ $\frac{N\left(N-1\right)}{2}$ correlation coefficients. This is due to a reduction in sFNC variance in patients, visible in bar chart format (Figure 3.c) and in the absolute difference between control and patient means (3.d, 3.e). This visible reduction in patient sFNC values provides initial support for the dysconnectivity hypothesis, as it suggests that schizophrenia is broadly characterized by reduced inter-network communication compared to controls. This, in addition to the block-diagonal structure of the sFNC matrices, suggests that the data have been correctly reconstructed.

3.2. Joint Entropy

We estimate each tFNCP’s subject-level entropy rate with a variation of the Kozachenko-Leonenko estimator [63,64] based on k-nearest-neighbor distances between sample points [65,66]. This estimation algorithm is necessary due to the unavailability of source signals’ complete probability distribution, which is necessary for an exact calculation of signal entropy. Fortunately, the statistical independence of the tFNCPs prevents any dependence between their respective entropy rates. As such, the joint entropy rate of each subject can be estimated simply by summing all tFNCPs’ entropy rates within that subject. The resulting group-level distributions of the joint entropy are displayed in Figure 2, with schizophrenia patients on the left and healthy controls on the right. Healthy controls display elevated joint entropy rate relative to patients, with the Kolmogorov-Smirnov test demonstrating that this elevation is highly significant $\left(p=4.506\times {10}^{-8},D^{*}=0.2188\right)$. A difference-of-means permutation test confirms this finding $\left(p=9.999\times {10}^{-5},G=-0.4095\right)$. Student’s $t$-test was not employed because the joint distribution fails the Jarque-Bara test, thus violating Student’s assumption of normality. Group-level means and standard deviations are listed in Table 1.

3.3. tFNCP Entropy

Having confirmed that subject joint entropy rate differs between populations, an obvious question is whether these differences are evenly distributed amongst tFNCPs or whether they concentrate within a subset of them. To answer this question, we compare the group-level entropy rates of each tFNCP using the Kolmogorov-Smirnov two-sample test, with the false discovery rate [61] correcting for multiple comparison. A Student’s t-test is also applied to entropy rate distributions which pass the Jarque-Bera test for normality. This approach reveals decreased patient entropy rates in five of eight tFNCPs (Figure 4), with no qualitative differences between Kolmogorov-Smirnov two-sample test and Student’s t-test results. Row 4.a of Figure 4 displays the connectivity matrices of each tFNCP, row 4.b the tFNCP connectograms, and row 4.c the group-level entropy distribution boxplots. Group entropy means, variances, and test statistics are listed in Table 2.

3.4. Multiple Linear Regression

The wide range of ages and seven collection sites of the FBIRN dataset raise the possibility of substantial age and site effects, potentially large enough to affect the outcome of the previously described hypothesis tests. To account for this, we employed a multiple linear regression analysis to separate the effects of site, age, gender, and diagnostic status (“patient” vs. “control”) on subject joint and tFNCP-level entropy rate. We used similar analyses to examine the link between subject-level entropy, clinical scores, and cognitive scores. Clinical scores consist of the PANSS positive and negative scores, while the CMINDS composite score is used as a surrogate for cognitive scores due to high correlation among the subscales (Figure 2).

3.4.1. Diagnostic Effects

Multiple linear regression diagnostic results, after controlling for site, age, and gender, remain consistent with statistical hypothesis tests. Site effects are both substantial and highly significant, but do not change the outcomes of hypothesis tests. Age and gender effects are found to be small and statistically insignificant.

3.4.2. Symptomatic Effects

Multiple linear regression finds no statistically significant relationship between patient PANSS scores and entropy rate in any tFNCP. Effect sizes are small, with only four of eighteen examined relationships displaying effect sizes of order ${10}^{-2}$ or higher.

3.4.3. Cognitive Effects

One tFNCP, number five (5) displays a small but highly significant $\left(\beta =0.0633,p=0.0053\right)$ relationship between subject CMINDS composite score [62] and entropy rate. Interestingly, this tFNCP is characterized by strong anticorrelations between two higher cognitive domain regimes: default-mode and subcortical networks oppose visual networks and sensorimotor networks, with cognitive control networks split between these two regimes. The default mode is the most prominent of these domains, being strongly anticorrelated with both opposing domains. Joint entropy rate is also positively related to the CMINDS composite score, although this relationship does not reach statistical significance $\left(p=0.1203\right)$.

4. Discussion

This article presents a novel method for evaluating the entropy of the dynamics of temporally independent FNCs and the relationship of that entropy to psychiatric disease. Of the eight temporally independent FNC profiles identified, five display significant entropic reductions in schizophrenia patients, with the remaining three showing no such reduction. One of these five diagnostically variable tFNCPs displays a significant positive relationship with cognitive score, with this single tFNCP including strong contributions from default mode networks, cognitive control networks, and other higher cognitive domains.

Initial inspection of the group-level static functional network connectivity broadly supports previous work showing reduced global connectivity in schizophrenia patients [67]. Positive correlations between auditory, sensorimotor, and visual networks appear notably attenuated in patients compared to their healthy counterparts. This reduced connectivity strength and structure may reduce patients’ ability to rapidly reconfigure network relationships [68,69].

Extant research in dynamic functional connectivity consistently suggests an attenuation in the complexity of connectivity dynamics in psychiatric disorders, particularly in schizophrenia. Prior work [32] has shown that patients explore less of their connectivity state space than comparable healthy controls do, and that their trajectories within this state space are both more rigid and more repetitive than controls as well. These findings imply a simplification of dynamics in patients, in line with the reduction in entropy described in this article. Such simplification may underlie the cognitive rigidity and inflexibility known to characterize some psychiatric disorders. Interestingly, the findings of Miller et al., 2016 and the present article contrast with other research which has found that intake of psychedelic substances increase activity diversity, turbulence, and information flow compared to placebo [70–72]. This suggests that hallucinatory experiences in the psychiatric and psychedelic realms possess fundamentally different operational mechanisms. Future research should attempt to identify the mechanisms in each of these populations. The finding of reduced entropy in patients may also align with the finding that psychiatric patients operate in a less metastable regime than their healthy counterparts [73,74]. Similarly, it may imply that patients occupy a reduced functional manifold compared to controls [75]. The concentration of entropy declines in specific tFNCPs may suggest an alteration in the connectivity’s energy landscape as well, as described in [72]. Which, if any, of these hypotheses can best predict the subject trajectories in the tFNCP basis space may prove a valuable line of research.

Of the eight tFNCPs identified in this study, five display significant, positive, linear relationships between entropy rate and the diagnosis of schizophrenia. One of these further displays a small but highly significant positive relationship between general cognitive function (CMINDS score) and entropy rate. Within these five tFNCPs, the most notable feature of the functional domain structure is the stability of the visual-sensorimotor domain block. In four of the five tFNCPs with significant entropy-diagnosis relations, these domains strongly correlate with one another. Further, of these four tFNCPs, three extend the block to include the auditory domain. Only tFNCP eight does not display a strong visual-sensorimotor domain block. Instead, tFNCP eight mostly displays control and visual domain interactions with multiple other domains, mostly visual.

Of the five tFNCPs with altered entropy in schizophrenia patients, each displays a different block to anticorrelate, or oppose, the dominant visuo-sensorimotor-auditory domain alignment. Cerebellar and subcortical networks primarily make up tFNCP four’s opposition, along with most of the default-mode and cognitive control domains. Three default-mode and cognitive control networks align with the visuo-sensorimotor blocks as well. tFNCP five, which appears to influence CMINDS score, recruits the subcortical and default-mode networks in opposition to the sensorimotor-visual-auditory block. tFNCP seven appears to possess three distinct domain-level modules, with cognitive control and default-mode networks correlated with themselves while subcortical networks oppose the sensorimotor-auditory bloc. tFNCP eight, as noted, shows links between cognitive control, default-mode, and visual domain areas.

When observing all eight tFNCPs as well as the diagnostically significant five, two notable trends emerge. First, no tFNCP displays strong intra-domain correlation. High correlation values instead concentrate between domains, as seen in the cerebellar-subcortical block of tFNCP three (3). In fact, only two tFNCPs display high intra-domain correlation: numbers 4 and 5. Even within these tFNCPs, this is the exception rather than the rule, with tFNCP 4 showing high intra-domain correlation only in sensorimotor and visual domains and tFNCP 5 only in the default-mode networks. This consistent lack of strong intra-domain correlation stands as a notable contrast to the static FNC of both groups, the strongest connections of which are all within domains.

The second notable trend across all eight tFNCPs is that cognitive control networks seem to straddle the line between the sensorimotor-visual-auditory block and its opposition. This domain is generally divided between more evident functional modules, but its strongest correlations are consistently towards the sensorimotor-visual block even in cases where most of its networks align against it. Its nonaligned status is most apparent in tFMCP seven, in which the cognitive control and default-mode networks form a third block which largely correlate with neither the sensorimotor nor cerebellar functional modules. The cognitive control domain also displays a distinct tendency to split into two largely separate networks, with relatively weak autocorrelation and strong correlation or anticorrelation with other domains being present in four of five significant tFNCPs.

5. Conclusions

Overall, this article proposes a novel means to quantify entropy in neuroimaging signals. Its basis in well-established methods for identifying functionally relevant spatial maps [34,42,76,77] and maximizing independence of temporal signals [49] lends credence to its findings of reduced dynamic connectivity trajectory complexity in schizophrenia patients. This reduction appears specific to identifiable tFNCPs, one of which significantly relates to cognitive function. It thus represents a potential link between the fields of neuroimaging and psychological analysis and is a substantial step in the ongoing search for biomarkers of psychiatric disease.

Future work with this and related methods should focus on bridging the gap between computational psychiatry, neuroimaging, and dynamical systems models. The ability to directly identify the time courses and complexity of specific tFNCP trajectories means that similar methods may be able to complement theoretical work on cognitive manifolds [75] or confirm extant studies on the transition networks and energy landscapes of the human chronnectome [78–82]. Similar methods may also be able to identify functionally related subgroups within heterogeneous disease categories, as has been done with certain psychological models [1,83]. This, in turn, could lead to developing biomarkers of symptom dimensions directly, which may provide a more reliable and intuitive classification system for psychiatric distress than the extant Diagnostic and Statistical Manual of Mental Disorders (DSM-V) or International Classification of Diseases (ICD). We thus anticipate extensive use of this and related methodologies in future studies as part of the wider quest to catalogue and diagnose psychiatric and neurological dysfunction.