|Home | About | Journals | Submit | Contact Us | Français|
We investigated neural regulation of emotional arousal. We hypothesized that the interactions between the components of the prefrontal-limbic system determine the global trajectories of the individual’s brain activation, with the strengths and modulations of these interactions being potentially key components underlying the differences between healthy individuals and those with schizophrenia. Using affect-valent facial stimuli presented to N = 11 medicated schizophrenia patients and N = 65 healthy controls, we activated neural regions associated with the emotional arousal response during fMRI scanning. Performing first a random-effects analysis of the fMRI data to identify activated regions, we obtained 352 data-point time series for six brain regions: bilateral amygdala, hippocampus and two prefrontal regions (Brodmann Areas 9 and 45). Since standard statistical methods are not designed to capture system features and evolution, we used principal component analyses on two types of pre-processed data: contrasts and group averages. We captured an important characteristic of the evolution of our six-dimensional brain network: all subject trajectories are almost embedded in a two-dimensional plane. Moreover, the direction of the largest principal component was a significant differentiator between the control and patient populations: the left and right amygdala coefficients were substantially higher in the case of patients, and the coefficients of Brodmann Area 9 were, to a lesser extent, higher in controls. These results are evidence that modulations between the regions of interest are the important determinant factors for the system’s dynamical behavior. We place our results within the context of other principal component analyses used in neuroimaging, as well as of our existing theoretical model of prefrontal-limbic dysregulation.
Over the past decade, significant research efforts have isolated and emphasized the role of the prefrontal cortex (PFC), the hippocampus and the amygdala in emotion regulation, fear conditioning and extinction.
The predominant view is that the central amygdala (CE) forms the main excitatory component of the arousal response (Sotres-Bayon et al., 2004); in both human and animal studies, damage to the amygdala prevents the acquisition and expression of fear. When an emotionally salient stimulus is detected, its signal is transmitted via thalamic pathways to the lateral amygdala (LA), then to CE (either directly or via more complex intra-amygdalar connections). CE in turn projects to a set of regions that control specific autonomic, endocrine and behavioral responses. It has been suggested that, in addition, the basal amygdala provides contextual contributions, integrating information from the LA, the hippocampus and the medial prefrontal cortex (mPFC) (Sotres-Bayon et al., 2004).
The primary inhibitory pathways are the medial prefrontal cortex (Izquierdo et al., 2005; Phelps et al., 2004; Rosenkranz et al., 2003; Sotres-Bayon et al., 2006) and the hippocampus (Corcoran, 2005; Sotres-Bayon et al., 2004). Damage to the mPFC is known to generally induce emotional and cognitive changes, which seem to be very finely-tuned and subregion-specific — e.g., the dorsal part of the nucleus accumbens is involved in attention and cognitive control, and the ventral part, in emotional regulation (Bush, 2000). The classical view is that neural activity in the mPFC regulates not only the amygdala-mediated fear responses via direct projections to the LA or the intercalated cells, but also the activity in the hippocampus, via projections to CA1 (Curnu Ammonis, area 1 — see appendix) and its own activity (Mora and Myres, 1977; Ferrer et al., 1993).
In addition to the roles of the amygdala and PFC in emotion regulation, which have been confirmed by a wide variety of studies, the potential contribution of the hippocampus has been increasingly explored. Structural MRI studies found decreased hippocampal volumes in depressed (Caetano et al., 2004) and schizotypal patients (Dickey et al., 2007), and correlated the volume loss with the length of the illness. This is consistent with the hypothesis that hypercortisolemia induced by stress vulnerability could result in hippocampal neurotoxicity in conditions such as bipolar disorder and schizophrenia (Rǎdulescu and Mujica-Parodi, 2008). The damage is believed to be restricted to particular areas or features (Sousa et al., 2000; McEwen, 2001), such as the synaptic plasticity within CA1 (Kemp and Manahan-Vaughan, 2007) or the hippocampus-PFC interactions (see Appendix). Newer imaging (Malaspina et al., 2004; Tamminga et al., 2003) and animal studies (Lodge and Grace 2006, 2007) suggest a complex hippocampal modulation of the dopamine systems, and are focused on isolating the signature of this modulation in conditions of emotional dysregulation.
The dynamics of prefrontal-limbic regulation are therefore important in understanding healthy emotion, but perhaps even more crucially important in establishing the etiology and neurobiology-based diagnosis of mental illness. A mild to moderate dysregulation of emotional arousal may be associated with substance abuse (deArcos et al., 2008), or may underlie pathologies such as personality disorders (Newhill et al., 2004). A more severe dysregulation of the inhibitory feedback loops is believed to lead to illnesses such as bipolar disorder or schizophrenia (Aleman and Kahn, 2005). Schizophrenia has many neurobiological features suggesting an underlying dysregulation of emotional arousal, including limbic (Williams et al., 2004a), endocrine (Ritsner et al., 2004; Tandon et al., 1991), and autonomic (Dawson et al., 1994; Mujica-Parodi et al., 2005) abnormalities. While neurotransmitter hypotheses (Joyce and Milan, 2005; Yamamoto and Hornykiewich, 2004) of schizophrenia are well-established and form the bases for development of newer antipsychotics (Bilder et al., 2002; Kapur and Mamo, 2003), it is possible that the neurotransmitter abnormalities are consequent to hyperarousal (Finlay and Zigmond, 1997; Jackson and Moghaddam, 2004; Moghaddam and Jackson, 2004).
Indeed, one known aspect of schizophrenia dysfunction is the difficulty that patients have in recognizing facial emotions, particularly negative expressions, such as anger and fear (Russell et al., 2007). It is probable that the deficits of schizophrenia patients in facial affect recognition correspond to structural and neurophysiological abnormalities in the brain regions involved in these processes (Holt et al., 2006). However, the nature of this deficit and its neural bases are the subject of ongoing research, and conclusions do not always agree.
While some imaging studies of schizophrenia patients revealed volume reductions (Namiki et al., 2007) and attenuated responses of the amygdala (Johnston et al., 2005, Aleman and Kahn, 2005) and mPFC (Takahashi et al., 2004), other studies found elevated hippocampal and amygdalar activity in patients during the passive viewing of faces (Holt et al., 2006), or during discriminating emotions (Gur et al., 2002; Kosaka et al., 2002). Indeed, the quantitative tests usually employed in conjunction with structural and imaging data are statistical tests of mean activation in different brain regions — and are intrinsically limited in power, since they can’t capture time-evolutions and network features. Consequentially, their results may appear sometimes contradictory, as the mean statistics change with time and with the specific circumstances. Relatively little research has been conducted to reconcile these discrepancies among different studies. Relating them to interactions between brain areas and to the concept of time-evolutions could unify the results (Hempel et al., 2003; Fahim et al., 2005; Williams et al., 2004b; Aleman and Kahn, 2005; Das et al., 2007). Since the human brain can be viewed as a dynamical system evolving under constant external perturbations, we propose that a dynamical, system-based analysis is more appropriate, comprehensive and consistent. We emphasize in particular the concepts of “time” and ”brain network,” in contrast with previous studies, which have largely focused on static measures of one variable at a time (BOLD amplitude, connectivity).
Our approach is faithful to this perspective of network dynamics, and is a natural continuation of our prior work. We have previously used cross-correlation coefficients and other related measures to investigate prefrontal-limbic regulation in healthy controls (Mujica-Parodi et al., 2009) and schizophrenia patients (Rǎdulescu and Mujica-Parodi, 2008). However, cross-correlations describe only pairwise relationships (i.e., relate the evolution of only two network components at a time), while the existing literature implicates more potentially relevant regions with significant dynamics. Our own imaging work has found six such possibly inter-related regions. It became therefore necessary to start approaching higher-dimensional data, in order to simultaneously address the dynamics of the entire network — and not constrain ourselves a priori to relationships between node-pairs.
Recent research has produced considerable evidence on network interactions between certain brain regions. Some of the pathways that constitute the neurophysiological support of these modulations have already been isolated by basic and clinical research (see Appendix). The results we present in this paper support this existing evidence. We show that the ROI activations in two populations (of healthy controls and schizophrenia patients) exhibit very specific patterns, which we interpret as marks of the underlying intrinsic dynamics of the system.
Our first hypothesis is that, although the initial network is six-dimensional, the relevant phenomena may happen in a lower-dimensional space. We expect that this may generally be the case in ROI networks, and that reducing the dimension of such complex systems would make their analysis more tractable. According to this hypothesis, the dimensionality of a network defined in an a priori or exploratory manner can be encompassed in only a few relevant components, which would thus suffice to encode the disparate contributions of each node to the overall behavior. In the present study, we investigate whether this is the case in the six-dimensional ROI network (bilateral amygdala, bilateral hippocampus, and two prefrontal regions) established by our previous work (Mujica-Parodi et al., 2009, Rǎdulescu and Mujica-Parodi, 2008), for both healthy controls and patients.
Our second hypothesis is that the prefrontal-limbic system has a strongly deterministic component which is identifiable mathematically, i.e., that the interactions within the system and with the environment determine its evolution over time. New techniques such as Dynamic Causal Modeling (Friston et al., 2003) make this assumption. A positive conclusion to this hypothesis (i.e., evidence that a dynamical systems analysis of such a network is in this case applicable) would establish a more solid basis for such brain connectivity methods. This would implicitly encourage and validate top-down computation of dynamic invariants, such as Lyapunov exponents, or attractor dimension, which are ideal mathematical characterizations of systemic homeostatic regulation.
The third hypothesis of our paper is that there are discernable control versus patient group differences, the identification and understanding of which are essential in establishing better diagnostic and treatment options. While identifying and explaining some of the patterns seen in our imaging data, we concurrently attempt to find similarities and differences between the control and patient populations. In the light of our dynamical approach, we further speculate that our six-dimensional system may have a lower dimensional attractor, which would be a dynamical signature for the homeostatic regulation of healthy controls and schizophrenia patients, respectively. Therefore, differences in regulation would be reflected by differences in the geometry and shape of the attractor.
We tested 11 patients (8 males, 3 females; age μ = 35, σ = 17.1) and 65 control subjects (28 males, 37 females; age μ = 26 yrs, σ = 8.4. We used a sample of convenience obtained from a prior study on individual variability among N = 65 healthy controls (Mujica-Parodi et al., 2009). This was complemented by acquisition of data using the identical paradigm on a smaller number of patients (N = 11), who were much more challenging to recruit and test. In our statistical analysis, we compensated for the size imbalance by using size-independent, non-parametric comparison tests.
All 76 subjects were between the ages of 18-50, with normal hearing and correctable vision. Subjects were excluded if they had a history of neurological or cardiac illness, substance abuse, presence of metal in the body, were pregnant, or taking medication with known interference to accurate assessment of arousal measures (including benzodiazepines). The screening of subjects prior to participation included a medical history and physical exam, conducted by study physicians. Diagnoses for patients and screening for healthy subjects were confirmed using the Structured Clinical Interview for DSM-IV — SCID-I (Ventura et al., 1998). Patients were recruited from the Stony Brook University Hospital Adult Inpatient, Outpatient, and Day Units; control subjects were recruited from the general community. Ideally, we would use unmedicated patients. However, due to the practical and ethical difficulties involved in taking psychotic patients off of their medications, we used instead a medication-heterogeneous sample to avoid medication-specific neural effects. Symptom severity for these patients was assessed using the Positive and Negative Syndrome Scale (Kay et al., 1989) (N = 9; positive symptoms: μ = 23.3, σ = 7.5; negative symptoms: μ = 26.4, σ = 6.3; general symptoms: μ = 46.8, σ = 16.9; total scores: μ = 96.6, σ = 29.1). This study was approved by the Stony Brook University Institutional Review Board. All subjects provided informed consent. Patients had an additional capacity evaluation signed by their treating psychiatrists.
Subjects were scanned on a 1.5T Philips MRI. Time series were acquired using two blocks of 136 T2*-weighted echoplanar images covering the frontal and limbic areas of the brain, with TR = 2500ms, SENSE factor = 2, TE = 45ms, Flip angle = 90, Matrix= 64×64, 3.9mm × 3.9mm × 4mm voxels, and 30 contiguous oblique coronal slices (Mujica-Parodi et al., 2009).
All statistical analyses employed during image processing were performed using the general linear model in the Statistical Parametric Mapping software (SPM99). The raw functional BOLD images were first realigned to the first volume to remove movement-related artifacts using a sinc interpolation. The motion correction algorithm in the SPM99 software package is capable of correcting for motion within 3mm. Realigned images were then spatially normalized into 3 × 3 × 3 mm3 using an affine transformation with a set of 7 × 8 × 7 basis functions and a customized template that was created using the data for the first 12 subjects; the incomplete brain coverage and oblique nature of our slices required us to use a custom template for normalization. For each subject, the scalp was removed from a low-resolution EPI image, using the Brain Extraction Tool (BET) available in MRIcro software, at a fractional intensity threshold of 0.5. These skull stripped images were then registered and normalized to each other and the average image was smoothed with a Gaussian kernel of 8mm full-width half maximum and registered to the EPI template provided by SPM99 to generate the final template. The realigned and normalized time series were then smoothed with a Gaussian kernel of 8 mm full-width half maximum (Mujica-Parodi et al., 2009).
We used a boxcar design with four conditions (Anger, Fear, Happy and Neutral) convolved with the canonical hemodynamic response function. We presented four stimulus runs during each scan. Each run consisted of four blocks, each block depicting 10 neutral, angry, fearful and happy faces taken from the Pictures of Facial Affect (Eckman and Friesen, 1978). Each face was presented for 2s (20s total per block), alternating with a Rest condition, in which we presented a white fixation cross on a dark background. During all conditions, subjects passively viewed the stimuli. To insure that they were alert in the scanner, subjects were addressed by microphone by the staff between runs. After the scan, all subjects completed an exit interview, confirming their attention to the stimuli and alertness throughout the scan.
The amygdala and hippocampus were defined a priori as limbic regions of interest (ROI). Separate ROI masks for the bilateral amygdala and hippocampus were traced on the standard T1 template provided with the SPM package, using established anatomical landmarks to delineate our ROIs (Hastings et al., 2004). Given the lack of certainty regarding the specific functionally-relevant areas for the prefrontal cortex (Sotres-Bayon et al., 2006), we defined these ROIs post hoc, by performing a random effects analysis over the entire control population. This showed two clusters centered in the ventromedial (Brodmann Area 45) and the dorsolateral (Brodmann Area 9) areas of the prefrontal cortex, as well as bilateral activation in the amygdala and hippocampus. The maximally-activated voxel for each cluster was then selected for the group and used to extract time series for each individual subject. Obtaining the time series in this manner allowed us to exploit the advantages of both the single voxel and cluster methods, which were optimized for signal and reliability, respectively. We first extracted the time-course for the whole experiment for a cluster-defined single voxel in each of the ROI. Signal values over the peristimulus intervals — 11 time-points, 3–33s from the onset of each block — were obtained for each condition block (the canonical hemodynamic response function, convolved with the 20s block design, was 30 seconds, and thus introduced a lag for each condition).
For each individual, we thus created a 352 data-point time series for each of the six ROIs, spanning four identical repeat runs of Rest Anger Rest Neutral Rest Happy Rest Fear (each condition block corresponds to 11 data points, adding up to 88 data-points per each eight-block run). Since we were interested in the time series dynamic variability, rather than in the ROI activation amplitudes, used these expanded 352 data-point time evolutions, which consider all emotional conditions, including rest. We then further processed these inclusive series in ways justified by our preliminary analyses, as described in Section 2.3.
Our analysis was done in a six-dimensional space, in order to simultaneously study the ensemble of all mutual interactions and modulations between regions.
First, we performed a preliminary descriptive statistics, and calculated for each individual the temporal mean activation of each ROI over the entire scan (i.e., over 352 data-points). We had two reasons in mind: (1) to estimate the variability of the mean activation among all subjects; (2) to determine whether we can simply use these mean activations to distinguish between the control and patient groups. For each subject, we obtained a six-dimensional point which represents his/her mean brain activation, irrespective of the condition. We analyzed the group distributions of these mean activations, first using standard descriptive statistics, and then by using a principal component analysis.
First, we performed a rank test on each of the six ROI mean activations, but didn’t find any significant statistical differences between the control and patient groups. This was hardly surprising, given the limitations of mean and maximal activation statistics. As explained in the Introduction section, these standard statistics are not expected to capture the network behavior or the evolution in time of the ROI activations (which, according to our hypothesis, underlie the group differences). Of the existing methods that could capture simultaneously the evolution of all nodes in our network, we chose to use PCA as a computationally straight-forward, yet effective tool for data-reduction.
We performed a PCA on the six-dimensional distribution of means, separately for each one of the two groups (controls and patients). We found that for controls, most variability in performance appeared along PC1 ~ (0.12, 0, 0, 0.99, 0, 0)1 (stdev=45.58) and along PC2 ~ (0, 0.16, 0.98, 0, 0, 0) (stdev=14.23). Within the patients group, PC1 ~ BA9 (stdev= 4.5), while all other principal components were comparatively insignificant (see Figure 1).
The goal of our analysis was to transcend the variability in the subjects’ means (as expressed by this preliminary statistics), and to pin down the actual intrinsic features of the system’s dynamics. The high variability in activation levels among subjects, even within the same group, was addressed in two different ways, each offering a different perspective on our data (Section 3.1). We studied (1) contrasts of activation minus rest and (2) prototypical group trajectories. Using both approaches, we performed principal component analyses in the six-dimensional phase-space of our data. Firstly, we searched for any constraints imposed on the system’s behavior by the intrinsic dynamics, as they might be captured by the contrasts or the prototypical group trajectories (Section 3.2). Secondly, we were interested whether the PCA analysis of contrasts (calculated for individuals, rather than prototypes) offered any potential in differentiating between the control and patient groups (Section 3.3).
We created contrasts by subtracting from each arousal block the preceding block of rest. This generated new time series of 176 data-points (each of the four runs contained four contrasts of 11 data-points per contrast). After performing this normalization, a new analysis revealed very small variations in the means of all variable components (see Figure 2).
For a fixed region of interest, we defined two types of prototypical time series: “temporal templates” and “group prototypes”. Firstly, taking a longitudinal approach: by averaging over a large number of identical runs for the same subject, the average time series can be considered to be representative for that individual (“template”), in the sense that it captures some important features of that individual’s general brain dynamics. Secondly, taking a cross-sectional approach: for a large enough group of subjects, the population-average time series can be considered to be representative (“prototypical”) for the whole group, in the sense that it is expected to capture some important common features of the brain dynamics within that population. The noise due to subject variability is expected to be in large part cancelled out by averaging.
In this context, we note that prototypes are affected by the sample size, as is the case with basic statistics, so the technique does not compensate for small samples. In general, it may be difficult to provide enough repeated identical runs to create temporal templates, or enough subjects in each group to create reliable group prototypes. This is true for our imaging paradigm, which employed four repetitions and a total of N = 76 subjects. For our further analyses we used a combination of the two averaging techniques, which we called “temporal group prototypes”. For example, Figure 3a shows the prototypical control performing under a typical Rest Anger Rest Neutral Rest Happy Rest Fear run (the average was taken over the four identical runs, and spans 88 data-points). The corresponding prototype curves for patients were not as smooth (Figure 3b), which could either mean that the patient data were more chaotic, or could be a sign that the patient population (N = 11) was too small to create a reliable prototype even after some initial temporal averaging. We ruled out the former (an intrinsic higher irregularity in patients) by using the approximate entropy (Pincus, 1991, 1993; Pincus and Singer, 1996) to compare the complexity of the time series in the two populations (Table 1). The computation did not show enough significant differences to justify the group differences in and of themselves.
The variability had to be therefore due at least in large part to the small size of the patient sample. With this in mind, the temporal averaging was made more effective by disregarding the emotional condition, and further averaging over the four Rest Arousal blocks in each run. This way we obtained shorter temporal group prototypes (22 data-points), where the noise due to variability was further eliminated, in patients as much as in controls, rendering the prototypes comparably smooth (Figure 4).
While Arousal minus Rest contrasts have been frequently used in the study of fMRI time series, population-averages are not as common. Imaging research has been focused on studying statistical features of populations, especially in the context of identifying psychiatric conditions based on these statistics. We argue that strong population common features can be advantageously brought out by group prototyping in a clean and simple way, avoiding a complicated search through the thicket of less obvious statistical group properties. Any results obtained on prototypical series could be used to focus on more specific dynamics within each group and between groups, and could be further subjected to more detailed statistical tests. This is the approach we have taken in our study, in order to verify whether our prototypical results remained valid as non-prototypical group features. A few recent studies have used group averaged imaging data in similar ways in order to correlate them best with theoretical modeling work (Claassen et al., 2008).
We considered the six-dimensional 88 data-point temporal group prototypes shown in Figure 3 (by averaging each ROI’s time series over the four repeat run, then averaging over the control and patient populations, respectively). On these prototypes, we performed a principal component analysis. For both groups, the first two principal components were much larger than the other four, so that the trajectories can be viewed as almost two-dimensional. In the case of controls, these two principal directions were: PC1 = (0.43, 0.37, 0.43, 0.33, 0.41, 0.44) (with loading 0.189, accounting for 95% of the total variance) and PC2 = (0.36, 0.27, -0.65, 0.50, 0.19, 0.26) (with loading 0.006, 3% of the total variance). The other loadings were ≤ 0.001, i.e., 2% of the total variance altogether. In the case of patients, the two principal components were: PC1 = (0.51, 0.47, 0.30, 0.28, 0.40, 0.40) (with loading 0.277, accounting for 87% of the total variance) and PC2 = (-0.83, 0.29, 0.30, 0.25, 0.07, 0.23) (with loading 0.017, 5% of the total variance). The other loadings were ≤ 0.01,corresponding to a portion of 8% of the total variance. In Figure 5, we identified the PC1-PC2 sections for the two groups, in order to show the control and patient two-dimensional trajectories in the same phase-plane.
In order to verify that the reduced dimensionality result is validated by our other proposed approach, we considered contrasts of Arousal minus Rest, instead of prototypes. For each subject, we considered the first the 176 data-point contrast time series, and second the 44 data-point time series (averaged over the four runs), each containing four contrasts of 11 data-points per contrast. We performed similar PCA analyses on these contrast series, without additionally employing any group prototyping. Figure 6a and Table 2b show the PCA loadings for the distribution of 176 contrast data points for all controls (green) and all patients (red). With the same color coding, Figure 6b and Table 3 show the PCA loadings separately for each of the four contrasts, in all controls and all patients. As before, the loadings decay rapidly, so that in all cases, only the first and the second can be considered relevant to the subjects’ trajectories, accounting together for between 65% and 78% of the total variance. More precisely: in controls, PC1 and PC2 represented respectively 51% and 16% in the first contrast, 53% and 16% in the second, 51% and 14% in the third and 50% and 15% in the fourth. In patients, PC1 and PC2 represented 60% and 16% in the first contrast, 66% and 12% in the second, 55% and 19% in the third and 55% and 21% in the fourth.
Aside from the dimensionality result common to control and patients, the plots and tables in Section 3.1 also suggest a few group differences. First, the PCA loadings in patients seem to be higher than the ones in controls. Second, the PC1-PC2 planes seem consistently different between the two groups, even when analyzed for identical contexts: the left and right amygdala PC1 coefficients are substantially higher in the case of patients, the prefrontal and hippocampal PC1 coefficients are higher in the healthy controls (much higher for the prefrontal areas, and moderately higher for the bilateral hippocampus). In order to verify these results, we perform below more detailed group statistics.
A Kruskal-Wallis rank test on the loading values for controls and patients revealed that the apparent difference seen in the plots was in fact not significant, so it could not be used as a group differentiator.
However, a similar rank test applied on the PCA coefficients revealed significant group differences. For each individual, we used the long 176 data-point contrast, without further temporal averaging. The most prominent difference was observed in the right amygdala coefficients, which were significantly larger in patients (p = 0.02, H = 5.2, df = 1). Other marginally significant differences appeared for the left amygdala (higher coefficients in patients, p < 0.09, H = 2.76, df = 1) and for the Brodmann Area 9 (higher in controls, p = 0.07, H = 3.22, df = 1). The other coefficients were not significantly different (p > 0.10).
Previous imaging work has already proposed the use of principal component analysis as an exploratory, model-free technique in interpreting the shape of hemodynamic responses or differences between two responses (Friston et al., 1995). We would like to place the unique aspects of our approach within this context, address some limitations encountered by prior studies when using PCA, and explore possible extensions.
When applied to MRI, ordinary PCA has generally run into serious difficulties because of the extremely high number of dimensions in the data relative to the number of observations (Viviani et al., 2004). Because of these difficulties, PCA is often limited to regions of interest identified previously by the experimental model (Friston et al., 1993, Friston, 1997), as we have done ourselves in our analysis. This limits somewhat the exploratory quality of the analyses, but it makes the technique more meaningful and the results easier to interpret.
Such ROI based approaches have been employed to study the dynamics of a network of brain regions, and the results have been interpreted in terms of functional connectivity, defined as the correlation of a neurophysiological index measured in different brain areas. For example, Friston et al. (1993) uses a recursive principal component analysis on a large PET data-set to extract two independent principal components referred to as the “intentional brain system” and the “monotonic time effects”.
More generally, research in the last few years has placed a special emphasis on machine learning algorithms, adapted for more robust and performant pattern recognition in various forms of data. Not surprisingly, these studies have reported that the characteristics of the information extracted from the data are correlated with the capabilities of the classifier. However, complex multivariate techniques may be computationally expensive, especially when they have high predictive and classifying performance. The presence of redundant features in the data — natural or artificially added to compensate lack of information (Huang, 2007) — can seriously increase the computation time and cost, as well as degrade the algorithm’s accuracy. Since it involves a fast and inexpensive algorithm, Principal Component Analysis has been successfully used to reduce high-dimensional spectral data and improve the predictive performance of well known machine learning algorithms such as support vector machines (Cao et al., 2003; Lei and Govindaraju, 2005) neural networks (Khoshgoftaar and Szabo, 1994; Dorn et al., 2003) or others (Deegalla and Bostrom, 2006).
In imaging data in particular, multivariate techniques such as support vector machines and neural networks have been shown to be effective in classifications and pattern analyses of time-series (Wang et al., 2007; Weygandt et al., 2007). This stays true for clinical studies of psychiatric patients, where neural networks have been successfully used on fMRI data (Luo and Puthusserypady, 2007), but also on other clinical measures (Florio et al., 1994; Zou et al., 1996; Politi et al., 2005; Ho et al., 2005). Since PCA can be used for direct, immediate results, but also in conjunction with these other effective techniques, it is very important to establish its working parameters and limitations in the context of imaging data, as well as of other clinical measures.
Our main result is the identification of an important common feature of all subjects, both controls and patients: the evolution of the six-dimensional brain network we study is constrained to remain very close to the two-dimensional plane of its first two principal components PC1 and PC2. Our second result concerns the distinctions between controls and patients, which appear in the positions of the PC1-PC2 planes — in the sense that the left and right amygdala coefficients of the first principal component (which is by far the largest) are significantly higher in the case of patients, and the coefficients of Brodmann Area 9 are somewhat higher (trend) in controls. Our results support, in both controls and patients, brain connectivity models that assume causality of activation between regions. Since it is unlikely that such a robust property of the system’s trajectories can be accidental, we further suggest that it is the reflection of a deterministic trend in the system, as hypothesized in the Introduction.
Our main dimensionality result is very suggestive, particularly in the light of its relevance to inter-regional brain modulations. Firstly, it agrees unexpectedly well with theoretical hypotheses; our own formal model of brain regulation under stress (Rǎdulescu, 2008) proposes a two-dimensional Hopf bifurcation as the underlying essence of behavioral phenomena present in fear extinction. This matching experimental result encourages a whole continuing path of parameter estimation, system reconstruction and additional clinical validation of the model, in order to further investigate the properties of healthy and pathological brain modulations. Secondly, the result suggests that, when studying the evolution of a network of brain regions, the most important features of the trajectories may not appear in the pre-defined variables (the ROIs themselves), but rather in a subset of combinations of such variables (in other words, in a “subnetwork” of the original network). Therefore, when analyzing and interpreting the corresponding time series, it is important to work in the right eigenbasis in order to obtain the most relevant results.
Our second result, that the principal component coefficients differentiate between controls and patients, is quite unique in nature, since very few clinical studies of schizophrenia currently embrace a systems approach to prefrontal-limbic dysregulation. This is consistent with our previous work (Rǎdulescu and Mujica-Parodi, 2008), which suggests that patients show a predominantly excitatory response, with inadequate inhibitory modulation via the PFC. Our result implicitly suggests distinct attractor positions for the patients and controls. Future work should concentrate more on the temporal evolutions, on defining the dynamics and refining the geometric structure of these attractors (see Section 5).
The experimental conclusions agree well with the ones of our theoretical paradigm: although both healthy and patient prefrontal-limbic systems seem to constrain the trajectories to live in a two-dimensional plane, there are other intrinsic differences in the evolution of these trajectories that tease apart the two populations. These differences are not to be found in the mean activation statistics of different ROIs. Even more importantly, these differences are not best represented by the evolution of a particular ROI, but by measures that concern a combination of them (such as, in this case, the first principal component).
Our analysis captured an important dynamical feature of our network: all subject trajectories were almost embedded in a two-dimensional plane. Moreover, the largest principal component of the six-dimensional distribution of data-points differentiated significantly between the control and patient populations. These results are evidence that modulations between brain regions could be determinant factors for the system’s dynamical behavior and evolution.
Aside from the a priori choice of the ROIs and the typical pre-processing of the scans using the canonical hemodynamic response function (see the Section 2.2), our approach in this paper was designed to make as few dynamical assumptions as possible. The results support our theoretical model of neural vulnerability to stress in schizophrenic patients versus controls, which was previously developed independently of our experimental data (Rǎdulescu, 2008).
This encourages future efforts in the direction of bridging closer together the theoretical and data-driven approaches. A simultaneous dynamical system analysis would provide an interactive context for estimation and mutual validation of both these models. The ideal approach should adopt in parallel bottom-up and top-down methods, as the most conservative combination.
One aspect of our bottom-up system identification effort consists of using dynamic causal modeling (Friston et al., 2003), to estimate the strengths of the linear and quadratic modulations between brain regions. These should be sufficient for reconstructing the most important dynamics. Our top-down analysis toolbox includes computing dynamic invariants (Lyapunov exponents, attractor dimension) for both the theoretical and clinical systems, in both healthy controls and patients, and compare the results between models and between populations. These parameters are good measures of the complexity and chaos in the system as a whole (Faure and Korn, 2001), while maintaining the ability to characterize these by only a few numerical values — a key feature for clinical use.
Future work will expand the use of these techniques and include a larger number of patients, preferably unmedicated, to be able to separate out for symptom profiles and severity effects within the illness.
This research was supported by funding from the National Institutes of Health 5-MO1-RR-10710 (Stony Brook University General Clinical Research Center) and the Office of Naval Research N0014-04-1-005 (LRMP). The authors wish to thank John Stewart, R.N., Mark Sedler, M.D., and Andrew Francis, M.D. Ph.D. for their assistance with patient recruitment and screening. The authors declare no competing financial interests.
1All PCA coefficients are expressed as vectors of coefficients for the basis (LA, RA, BA9, BA45, LH, RH), such that (a1, a2, a3, a4, a5, a6) represents the six-dimensional vector a1*LA+a2*RA+a3*BA9+a4*BA45+a5*LH+a6*RH
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.