Intervention Models in Functional Connectivity Identification Applied to fMRI

Recent advances in neuroimaging techniques have provided precise spatial localization of brain activation applied in several neuroscience subareas. The development of functional magnetic resonance imaging (fMRI), based on the BOLD signal, is one of the most popular techniques related to the detection of neuronal activation. However, understanding the interactions between several neuronal modules is also an important task, providing a better comprehension about brain dynamics. Nevertheless, most connectivity studies in fMRI are based on a simple correlation analysis, which is only an association measure and does not provide the direction of information flow between brain areas. Other proposed methods like structural equation modeling (SEM) seem to be attractive alternatives. However, this approach assumes prior information about the causality direction and stationarity conditions, which may not be satisfied in fMRI experiments. Generally, the fMRI experiments are related to an activation task; hence, the stimulus conditions should also be included in the model. In this paper, we suggest an intervention analysis, which includes stimulus condition, allowing a nonstationary modeling. Furthermore, an illustrative application to real fMRI dataset from a simple motor task is presented.


INTRODUCTION
Functional magnetic resonance imaging (fMRI) based on blood oxygenation level-dependent (BOLD) signal has become one of the most prominent and powerful tools in cognitive neuroscience [1]. Most fMRI studies found in the literature focus on the detection of neuronal activation and brain mapping via statistical analysis. However, understanding cortical dynamics is a crucial step toward inferring cortical functioning.
Several evidences [2][3][4] suggest that modeling the interactions between different brain structures is paramount to understand the mechanisms guiding specific cognitive behaviour. However, the determination of parameters involved in cortical dynamics is still an open question. A number of techniques are being used to detect patterns of interaction between cortical areas, most using an ad hoc concept. So far, most connectivity studies have investigated temporal correlation as a measure of connectivity [3], even though it is not enough to identify the direction of information flow. In fact, Pearson correlation coefficient in time series analysis is just a measure of linear association. The connectivity mapping via correlation analysis is obtained firstly by selecting a seed voxel, and then Pearson correlation is calculated against all the other brain voxels. In most cases, the selection of the seed voxel is derived from the activation maps. Hence, as the activation detection is based on the similarity between the observed BOLD signal in a voxel and an expected haemodynamic curve, the correlation connectivity analysis is close to an activation mapping considering the seed voxel BOLD signal as the expected curve. Finally, we conclude that the correlation connectivity mapping is not sufficient to provide additional information in relation to the activation analysis based on general linear model (GLM).
Other statistical methods, such as the structural equation modeling (SEM), are more attractive to overcome this 2 International Journal of Biomedical Imaging shortcoming. Büchel and Friston [2] modeled the occipitoparietofrontal network involved in attention tasks using structural equation modeling. Zhuang et al. [5] applied SEM to a bimanual motor coordination experiment. Rowe et al. [6] modeled the prefrontal cortex in a color selection task. An improvement of SEM applied to fMRI analysis is the dynamic causal model (DCM), proposed by Friston et al. [7]. However, these two modeling approaches require a complete prespecification of the connectivity structure. Additionally, as DCM is estimated via Bayesian algorithms, it also requires the prior densities of the parameters of interest. In fact, these models measure the instantaneous connectivity, but require the direction of information flow. Therefore, these approaches are not enough to provide a complete identification of the connectivity pattern. Furthermore, the autocorrelation of BOLD signal is another obstacle for the application of these models and, in most cases, is simply ignored.
Granger causality [8] is a very prominent concept to describe information flow and connectivity. Baccalá et al. [9] and Baccalá and Sameshima [10] introduced a frequencydomain connectivity identification method for EEG using partial directed coherence and this causality concept. Goebel et al. [11] and Roebroeck et al. [12] introduced the concept of Granger causality in fMRI via vector autoregressive modeling (VAR). They have also shown the applicability of this approach to BOLD signals using simulations and illustrating it with real data derived from visuomotor studies. Compared to other approaches described previously, the main advantage of Granger causality identification via VAR models is the fact that prior specifications about connectivity structure are not necessary. Also, if one has prior partial knowledge about this structure, it can be naturally included in the model as a restriction in the parameters to be estimated.
The stationarity condition is one of the main obstacles to VAR modeling application in fMRI analysis. This assumption requires the connectivity structure to be the same during all acquisition times. Although acceptable in resting state or one-condition experiments, it may not be valid in paradigms with more than one condition. In that matter, as the connectivity structure may change according to the stimulus, comparisons between these structures are also an interesting point.
In this paper, we propose the use of a generalization of VAR models via intervention analysis (structural break models), which allows a natural modeling of connectivity and also statistical comparisons of the connectivity structures in experiments with different stimulus conditions.

GRANGER CAUSALITY AND CONNECTIVITY
In fMRI analysis, the definition of connectivity can be divided into two concepts: functional and effective connectivity. The first is defined as "correlations between spatially remote neurophysiological events" [4]. In contrast, effective connectivity is related to the "influence of one neural system over another" [4]. Note that effective connectivity implies in functional one, but the reciprocal may not be valid. The activity correlations or synchronisms may be observed due to ex-ternal factors, not only due to synaptic interactions between the areas involved. An illustrative example involving differences between these two concepts of connectivity could be a paradigm with simultaneous stimulation of visual and auditory cortex. The neuronal activity in these areas will be correlated (functional), but not related to neural interactions (effective). However, the simplicity of functional connectivity concept makes it very useful, mainly in cases where the neural activity is measured indirectly, as in fMRI time series.
Granger causality is a very useful concept for the description of brain areas connectivity and direction of information flow identification [10,13]. In the context of time series, Granger [8] defined causality in terms of predictability. This concept was originated in econometrics, focusing the understanding of relationships between financial time series such as prices, indexes, interest rates, and so forth. The basis of this concept is that effect cannot precede cause, suggesting that causality can be detected toward past and future relationships. A signal x t is said to Granger-cause a signal y t if the past values of x t help the prediction of present values of y t . In other words, if the variance of the prediction error of y t , considering all the information until the time t, is less than the one obtained excluding the information of past values of x t , then x t is said to Granger-cause y t .
Goebel et al. [11] introduced Granger causality identification between BOLD time series using VAR models, showing the applicability of this approach in simulated and real fMRI data. Note that Granger causality aims to identify interactions and relationships between signals via precedence and prediction. However, it is more related to functional connectivity than effective one, because it cannot distinguish real influences from prediction power.
Let a k-dimensional multivariate time series composed by k signals measured on time t. In order to measure the prediction improvement of Y t using a collection of p past values of the series where v is the intercept vector (related to the process average), u t is an error vector of random variables with zero mean and covariance matrix Σ given by v and A i are coefficient matrices given by João Ricardo Sato et al. 3 As the error term u t has zero mean, the predicted values of Y t conditional to the past values are given by The VAR model allows an easy way to identify Granger causality. If the coefficient a jli for some i is nonzero, we say that signal y lt Granger-causes the signal y jt . In other words, the past values of the signal y lt help the prediction of the present and future values of the signal y jt . It is important to mention that this kind of relationship is not reciprocal, for example, y lt may Granger-cause the signal y jt , but not necessarily y jt causes y lt , indicating the direction of information flow.
Consider a functional magnetic resonance dataset. Select k voxels in the volume, obtaining a BOLD k-dimensional signal. Using the concept of Granger causality and the VAR modeling, it is possible to verify if the BOLD signal of certain brain areas Granger-causes another areas' BOLD signal, by testing the significance of the estimates of matrix A t . Therefore, we are able to test the functional connectivity and direction of the information flow.
However, VAR modeling is only suitable in cases of stationary time series with coefficients and error covariance matrix invariant on time. Hence, considering fMRI studies, a weakness of this approach is the assumption that both activation and connectivity functions are constant in the whole scanning interval. Let an epoch functional magnetic resonance image experiment with two conditions A and B. It is reasonable to expect functional connectivity under A, but it may not be the same under B. Therefore, we propose the use of intervention VAR (structural breaks model) focusing on the identification of changes in functional connectivity. The intervention VAR model is defined by where the coefficient matrices are given by u t is an error vector of random variables with zero mean and covariance matrix Σ (C) defined by and C indicates the block condition (A or B). For simplicity, assume that for j = 1, 2, . . ., k, l = 1, 2, . . ., k, i = 1, 2, . . ., and δ (B) j , ∂ (B) jli , and ψ (B) jl are the increments on the coefficients during B condition. If at least one of the coefficients δ (B) j , ∂ (B) jli , and ψ (B) jl is nonzero, it implies the existence of structural changes. In other words, we have different coefficient matrices for each condition. Thus, we have a VAR structure for each block, but all the parameters are globally estimated, allowing a statistical test to the connectivity changes. The intervention VAR model can be estimated using an interactive generalized least-square estimator. Consider the following vector and matrices: where vec is an operator that concatenates all the columns of a matrix in a column vector, Δ is a vector of zeros and ones indicating the stimulus condition at time t (the tth element of Δ is zero/one if the acquisition at time t occurs during A/B condition), 1 is a column vector of ones, and ⊗ R is the row-Kronecker product, which is defined as the Kronecker product applied separately for each row, that is, The error covariance matrix is given by Γ. The generalized least-square estimator of the coefficients of the intervention VAR model is given by [14,15] However, as the covariance matrix Γ is unknown, we propose the use of an interactive two-stage least-square estimator. The residuals are estimated on the first step and the covariance matrix Γ on the second one, as an extension of the Cochrane-Orcutt procedure. The variances/covariances in Γ where σ 2 is the estimated residual variance, m is a vector, and C is a contrast matrix corresponding to the following test: Under the null hypothesis, W has an asymptotic chisquare distribution with rank (C) degrees of freedom [14]. Basically, this procedure performs simultaneous tests of the equality between m and linear combinations (C) of parameters in β. Note that the connectivity parameters can be easily tested considering the W statistic. Further, in case of non-Gaussian errors distribution, the martingales central limit theorem [16] implies that the classic asymptotic properties of the generalized least-square estimator are valid.
To finish this section, it is important to highlight some points about the application of the intervention VAR models to connectivity analysis in fMRI. Firstly, although Granger causality identification via VAR models is closely related to interactions, it cannot make a distinction between real influences or predictive power. Secondly, it is important to mention that this approach depends on sampling frequency. As Granger causality identification is based on information contained in past values, low sampling rates result in aliasing and data aggregation. In fact, sampling frequency represents a challenge for connectivity modeling in fMRI, as short acquisition time implies in low signal-to-noise ratio. On the other hand, neural interactions occur in a frequency much higher than fMRI acquisitions (TR). Thus, the connectivity structure identified via VAR models is only related to very low frequencies information, and fast interactions are not detected. Additionally, we would like to emphasize the point that although intervention model does not require the assumption of global stationarity, it depends on this assumption to be valid during each paradigm condition. The validity of this assumption is reasonable in block design paradigms, but it may not be true in event-related ones. The application of the proposed approach to fMRI data with event-related designs could result in an imprecise estimation, such as an average of the information flow intensity across the activation timepoints.

APPLICATION TO REAL DATA
In order to illustrate the usefulness of the proposed approach, the intervention VAR modeling was applied to real fMRI data derived from a motor task study.
Six normal right-handed subjects performed a simple right hand fingertapping task, in an AB periodic block design experiment. The functional magnetic resonance images were acquired in a GE 1.5 T Signa LX MR system equipped with a 23 mT/m gradient, (TE: 40 milliseconds, TR: 2000 milliseconds, FA: 90 • , FOV: 240 mm, 64 × 64 matrix; 15 slices, thickness: 7.0 mm, gap: 0.7 mm) oriented in the AC-PC plane in a single run. There were one hundred volumes collected during five cycles of rest-task performance. Each cycle had the duration of 40 seconds corresponding to 20 volumes, 10 volumes acquired during rest, and 10 volumes during the activation task. The subjects were in the dark room with noise-reducing headphones customized for functional MR. Instructions to begin and finish movements were given via auditory stimuli.
The images were preprocessed considering motion correction and spatial smoothing. The activation brain mapping was obtained using the XBAM software [17]. Spatial normalization transformation to the stereotatic space of Talairach and Tournoux [18] was performed using SPM2 [4]. Since the main interest was the study of motor function, the analyses were limited to the superior slices, reducing the number of multiple comparisons and consequently increasing the sensitivity of the statistical tests.
The group activation maps (cluster P value < .01) are presented in Figure 1 (top). Note that there are significant activations found on the left primary and contralateral motor cortex (BA 4), supplementary and premotor areas (BA 6), and primary sensitive areas (BA 3, 2, 1), which are classically involved in motor control.
Taking into account a significant activation in the supplementary motor area (SMA), the intervention VAR analysis was performed in a bivariate fashion considering a seed in the activation local maxima of SMA against all other voxels in the whole brain. As changes in BOLD signal are not instantaneous, we considered a delay of 2 seconds in the condition specification. Statistical tests for differences in the information flow intensity from SMA to the voxel of interest, between rest and task, were performed in each individual separately. The individual W statistics were mapped to Gaussian quantiles, and the group analysis was performed in SPM2 (on-tailed t test, which is similar to a random effects analysis [4]). The connectivity changes map (voxel P value < .005) is presented in Figure 1 (bottom). The areas with significant changes in the information flow intensity from SMA were the premotor cortex, presupplementary motor area, and primary sensitive cortex. The intervention VAR model was also applied in a trivariate analysis, considering the selected seed in SMA, and voxels (connectivity changes minimum P value) in the pre-motor cortex (PM) and pre-supplementary motor area (PSMA). In Figure 2 (top), diagrams (arrow P value < .05) describing the connectivity structure during rest and fingertap and also differences in information flow intensity are presented. The BOLD signals derived from selected ROI of one subject are presented in Figure 2 (bottom).
When contrasting fingertapping with rest we found significant changes in PM cortex, primary sensitive areas, and PSMA. All these areas are classically involved in movement control [19]. The current understanding of motor control in the literature suggests that PSMA provides the main input to SMA, which is possibly responsible for providing internal representation of movement sequences, and is involved in learning process of new movements. SMA is believed to send information to the PM and primary motor cortices, and is also associated with complex calculation to achieve maximum performance, based on feedback information from sensitive areas.
During the rest, we observed increased connectivity between the SMA and PM with pre-SMA. Some studies demonstrated that PSMA plays an important role in cognitive motor control, which involves sensory discrimination and movement decision making (go/nongo) or motor selection for the action after stimuli [20,21]. In order to start the sequence of fingertapping, a motor decision must be made, based on the instruction previously given to the participant. All areas (pre-SMA, PM, and SMA) are involved in the initiation of movement and modification in their connectivity pattern might be explained by the attention to stimulus presentation and monitoring during the task.

CONCLUSION
Most connectivity studies results rely on analyzing secondorder correlations that do not give additional information  Figure 2: Group connectivity structure during resting, fingertapping, and information flow changes between these two conditions (top). The numbers in the arrows describe the information flow correspondent P value. An illustrative chart of ROI's BOLD signals of one subject is also presented. about neural interactions. Other advanced methods like structural equation modeling (SEM) or dynamic causal models (DCM) could be an attractive alternative. However, they heavily depend on a prior knowledge about involved neural circuitry.
Granger causality concept, in its most general form, is a flexible definition of relationship and temporal order. It could be tested by a simple VAR modeling without any connectivity prespecification. In this paper, we introduced a new approach for connectivity modeling based on Granger causality and intervention VAR models, which enabled us to compare differences in connectivity structures as a statistical hypothesis test. An initial application of intervention VAR models to fMRI data produced biologically plausible results, but further experiments are necessary to reveal its potential as a new tool to investigate neural systems.