article: Two main approaches in exploring causal relationships in biological systems using time-series data are the application of Dynamic Causal model (DCM) and Granger Causal model (GCM). These have been extensively applied to brain imaging data and are also readily applicable to a wide range of temporal changes involving genes, proteins or metabolic pathways. However, these two approaches have always been considered to be radically different from each other and therefore used independently. Here we present a novel approach which is an extension of Granger Causal model and also shares the features of the bilinear approximation of Dynamic Causal model. We have first tested the efficacy of the extended GCM by applying it extensively in toy models in both time and frequency domains and then applied it to local field potential recording data collected from in vivo multi-electrode array experiments. We demonstrate face discrimination learning-induced changes in inter- and intra-hemispheric connectivity and in the hemispheric predominance of theta and gamma frequency oscillations in sheep inferotemporal cortex. The results provide the first evidence for connectivity changes between and within left and right inferotemporal cortexes as a result of face recognition learning.
  In order to exploit the full potential of high-throughput data in biology we have to be able to convert them into the most appropriate framework for contributing to knowledge about how the biological system generating them is functioning. This process is best conceptualized as first building a nodal network derived from empirically derived knowledge of the biological structures and molecules involved (nodes) and then secondly to use computational-based steps to discover the nature, dynamics and directionality of connections (directed edges) between them. Causality analysis based upon experimental data has become one of the most powerful and valuable tools in discovering connections between different elements in complex biological systems [1]–[6]. However, despite some encouraging successes in various areas in systems or computational biology its development and application have been impeded by a number of issues about the meaning of causality. For example, in clinical science, the current emphasis on how to apply causality approaches mainly resides in resolving the problem of how clearly to define causality itself [7]. A typical problem cited is the so called “Simpson paradox” in which the successes of groups seem reversed when the groups are combined. This demonstrates the ambiguity that can result in determining causal relationships based only on frequency data. However, this issue disappears if we incorporate time into the definition of causality as Granger has done in the field of Economics [8]. Nevertheless there is still no accepted unified way to tackle this issue. Taking altered gene expression data using microarray analysis as an instance, there are three approaches one can use to deal with the time-series data obtained: the simple dynamical system approach, the dynamical Bayesian network approach and the Granger causality approach which is a generalization of the dynamical system one. In [9], we have discussed in detail the pros and cons of applying the latter two approaches and shown potential advantages in using Granger when sufficient repeated measurements are available. With brain activity data from functional magnetic resonance imaging (fMRI) experiments, two prominent techniques have been introduced to address temporal dependencies and directed causal influences: Dynamic Causal (DCM) and Granger Causal (GCM) models. These two models have always been considered to differ radically from each other [10], [11]. DCM establishes state variables in the observed data and is believed to be a causal model in a true sense. On the other hand, GCM is a phenomenological model which just tests statistical dependencies among the observations to determine how the data may be caused [10]–[12]. The importance of the two approaches in interpreting fMRI data is demonstrated in [10] and in 2008 there were around 450 papers published devoted to both approaches and excluding those relating to other types of biological data. The key question we want to address in the current paper is whether we can develop an extended and biophysical constraint approach to share the features of the various approaches mentioned above, and in particular of the two causal models: DCM and GCM? The significance of such an approach is obvious and we would expect that its application could represent a powerful new tool in systems and computational biology, particularly in association with increasingly powerful genomic, proteomic and metabolic methodologies allowing time-series measurements of large numbers of putatively interacting molecules. In this paper, we will show that GCM can be extended to a more biophysical constraint model by incorporating some features of the bilinear approximation of DCM. By setting up a conventional VAR model with additional deterministic inputs and observation variables, we can create a more general model: Extended Granger Causal Model (EGCM) which offers a new way to establish connectivity. The EGCM is first tested in two toy models. With both state and observation variables, the interactions between nodes are successfully recovered using an extended Kalman filter approach and partial Granger to establish causality in both time (DCM) and frequency (GCM) domains respectively. The GCM approach itself is not tailored particularly well for biological experiments where we are often faced with the case of the data being recorded with and without a stimulus present. The time gap between two adjacent stimuli is very short and we would expect the network structure to remain unchanged during the whole experiment although the form and the intensity of the input may be unknown. This scenario is also the case for the gene network data considered in [1], [2], where the authors have treated the two situations separately, although there should be a common and true structure for both. We have therefore also used EGCM in toy models to establish its efficacy in revealing the true network structure when there is an intermittent input to affect state variables. To exemplify the direct application of EGCM to establishing causality in a specific biological system, we have applied it to local field potential (LFP) data recorded in the sheep inferotemporal cortex (IT) of both left and right hemispheres before and after they learn a visual face discrimination task [13]. There is electrophysiological, molecular neuroanatomical and behavioral evidence for asymmetrical processing of faces in the sheep brain similar to humans [14]–[16] although cells in both the left and right IT respond selectively to faces [15]. Learning also alters both local and population based encoding in sheep IT as well as theta-nested gamma frequency oscillations in both hemispheres and there is greater synchronization of theta across electrodes in the right IT than there is in the left IT [13]. There is considerable interest in establishing functional differences between the ways the left and right brain hemispheres interact and process information [17]–[19]. It has recently been hypothesized that the left hemisphere specializes in controlling routine and tends to focus on local aspects of the stimulus while the right hemisphere specializes in responding to unexpected stimuli and tends to deal with the global environment [18], [19]. Establishing altered causal connections and frequency dependency within and between the two hemispheres during face recognition learning will help test this hypothesis.
  All animal experiments were performed in strict accordance with the UK 1986 Animals Scientific Procedures Act (including approval by the Babraham Institute Animal Welfare and Ethics Committee) and during them the animals were housed inside in individual pens and able to see and communicate with each other. Food and water were available ad libitum. Post-surgery all animals received both post-operative analgesia treatment to minimize discomfort and antibiotic treatment to prevent any possibility of infection. The traditional and widely used Granger Causal Model takes the form [20]: (1) where, are coefficient matrices, is the noise, and the model has a vector autoregressive representation with an order up to. In spite of its successful application, GCM requires the direct observation of the state variables and does not include designed experimental effects in the model which form some of its limitations. Here we extend GCM to a more reasonable and biophysical constraint model by incorporating additional deterministic inputs and observation variables, closely following equations which are the features in the Dynamical Causal Model and its bilinear approximation form [10]. The extended Granger Causal Model takes the form: (2) where and are deterministic inputs, are the observation variables which are the function of the state variables, are the coefficients that allow the inputs to modulate the coupling of the state variables, and are intrinsic and observation noise and are mutually independent. Now, if we can recover the state variables from the noise observation variables, all the problems can be considered in the framework of the traditional Granger causality. It' s clear that the formal difference between EGCM and GCM is that we have included observation variables and deterministic inputs which are assumed to be known and will affect the connection of the state variables as well as the state variables directly. However, EGCM also has a strong connection with Dynamical Causal Model. We refer the readers to Discussion section for a detailed discussion on the importance of this extension, in particular the relationship between Eq. (2) and the Volterra type series expansion [10]. For the extended Granger Causal model Eq. (2), we now introduce an algorithm to estimate the state variables as well as all its parameters which will give us the first inspiration of the connection of the state variables. Let, , . Then, the VAR () model can be reduced to a VAR (1) model which takes the form: where is the parameter vector to be estimated. and are both zero-mean uncorrelated Gaussian noise with covariance matrix and respectively. In order to apply the model to real data, we have to estimate both the states and parameters of the model from input variables and noise observations. A widely used method for this dual estimation is extended Kalman filter (EKF) [21], . Here we recursively approximate the nonlinear system by a linear model and use the traditional Kalman filter for the linearized model. Let thenwhere is uncorrelated Gaussian noise with covariance matrix. Define where, . Then, the EKF algorithm for dual estimation consists of two steps: prediction and updating. After recovering the state variables using the EGCM algorithm above, we can define the causality with the idea proposed by Granger. The only difference is that, in our EGCM model, two deterministic inputs and are added to the normal autoregressive representation. Here, we provide the formulation of EGCM causality in both time domain and frequency domains. For simplicity of notation, here we only formulate EGCM for two time series and. To generalize them to more general case of multi time series, we refer the readers to [23], [24]. Assume that and in our EGCM have the following representation: (3) A joint representation in our EGCM that includes the past information of both processes and can be written as: (4) where is the maximum number of lagged observations included in the model. , are prediction errors with variance and are uncorrelated over time. Then, according to the causality definition of Granger, if the prediction of one process can be improved by incorporating the past information of the second process, then the second process causes the first process. So, in the extended model here, we define that if the variance of prediction error for the process is reduced by the inclusion of the past information of the process, then, a causal relation from to exists. This can be quantified as (5) If, there is no causal influence from to and if, there is. Similarly, we can define the causal influence from to as (6) Our EGCM also allows a frequency domain decomposition to detect the intrinsic causal influence which provides valuable information. We define the lag operator to be and assume here that the input is a constant, i. e. to avoid the appearance of nonlinearity. Then, the joint representation of both processes and in equation (4) can be expressed as: (7) Rewrite equation (7) in terms of lag operator, we have: (8) where, , , . Since what we really care about is the causal relationship caused by the intrinsic connection of the state variables rather than the outside driving force, i. e. the input, after fitting the model (7) and getting the covariance matrix of the prediction error, we just go on with the following model: (9) which means that after fitting the EGCM with input to eliminate outside influence, we just focus on the intrinsic causal influence in the frequency domain. After normalizing equation (9) using the transformation proposed by Geweke [25], [26] to eliminate the cross term in the spectra, we assume that we have the normalized equation in the form: (10) Fourier transforming both sides of equation (10) leads to (11) Recasting equation (11) into the transfer function format we obtain (12) After proper ensemble averaging we have the spectral matrix (13) where * denotes the complex conjugate and matrix transpose and is the covariance matrix of the prediction errors in equation (11). Hence, we can define the causal influence from to at frequency as (14) Similarly, we can define the causal influence from to at frequency as well. Note that although here we just provide the definition of pairwise Granger causality for EGCM, it' s obvious that similar methods can be easily applied to the definition of conditional, partial or complex Granger causality in both time and frequency domains [9], [23], [24], [27]–[30]. Since the explicit meaning of the parameters in the EGCM (i. e. the intrinsic coupling among state variables, the strength of the inputs to modulate the coupling and the influence of the inputs on the state variables directly), we can also get an idea of the connection of the state variables and how the inputs affect them from the fitted model before we translate it into a single number.
  Local field potential data were obtained from 64-channel multielectrode arrays implanted in the right and left inferior temporal cortices of three sheep (one sheep only had electrodes in the right hemisphere) as previously described [13]. Recordings were made while the animals were presented with pairs of faces which they were required to discriminate between using an operant response in order to obtain a food reward. In between face pair presentations the animals were presented with a visual fixation stimulus (a white spot on a black screen). Recordings were made during sessions of 20–40 trials where they were either still learning the discrimination or had successfully achieved the learning criterion (correct choice of rewarded face). In both the left and right IT the main oscillatory frequencies recorded are in the theta () and gamma () ranges and these two frequencies are coupled (theta phase and gamma amplitude) [13]. We have previously shown that learning increases theta amplitude, the ratio of theta to gamma, theta/gamma coherence and the tightness of theta phase [13]. With these experimental data, we can directly use our EGCM to detect the global network for all electrodes in both brain hemispheres. However, due to the size of the network, there are at least a few thousand free parameters to fit. To avoid this issue, we adopt another approach here. For each session we randomly select 3 time series from each region respectively and apply our model to detect the network structure for the six electrodes. This procedure is repeated for 100 times for each session (see Fig. 3 for such an example). The visual stimulus to the IT (including feedforward and feedback signals) is impossible to know in the experiment. However, as we have shown above, we are able to make the assumption that the effect of the stimulus can be regarded as a constant input to each electrode. The inclusion of the stimulus signal will certainly make the model more reasonable. A further problem here is that if we intend to reconstruct the connections for each six electrodes (left and right) before and after the stimulus respectively, we could end up with two different structures for the time series (not shown). This is certainly not the case since the duration of the stimuli is quite short (1–3 seconds) and the connections will not change in such a short time. To recover a reasonable structure of the connection in these areas in the brain, we therefore assume here that the connections in each trial don' t change and the time series before and after the stimulus are generated from a unified structure. With the application of our EGCM approach, we can include the intermittent stimulus and obtain a comparatively reliable structure. Fig. 3 (top-panel) shows such an example where three electrodes in the left and three in the right are randomly selected and that inter- and intra-hemisphere interactions are detected. Fig. 3 (bottom-panel) is the corresponding frequency decomposition of the top panel. In Fig. 4 we show the mean connections within and between the left and right IT calculated using EGCM and as a function of learning. The results clearly demonstrate an asymmetry between the hemispheres. The top-panel is an illustration of the bottom panel which summarizes the results of all experimental data for the two sheep. The most noticeable change is a decrease in the number of connections from the left to the right and an increase in connections within the right but not in the left IT. Indeed there was a strong negative correlation between the number of left to right connections and the number within the right IT for both animals (Sheep B, (); Sheep C, () ). These changes occurred as soon as the learning criterion was successfully achieved in successive blocks of trials (i. e. in as little as 5–10 minutes in the case sheep B where learning was successfully achieved within a specific recording session) and were maintained after 1 month or more post-learning. They were not simply the results of stimulus repetition because in Sheep B where recordings were made in repeated sessions of up to 120 trials but where the learning criterion was not achieved for one of the face pairs there were no connectivity changes observed. One of the advantages of the extended approach is that we have a frequency domain decomposition. Brain rhythms, not surprisingly, have also been intensively investigated in the literature [35]. Here we concentrate on the two main frequency bands: theta band () and gamma band () present in our IT recording data and which have been extensively linked to mechanisms of learning and memory [13], [35]. Using the frequency decomposition of our extended model discussed in Methods section, we looked at the following two quantities: Fig. 5A shows the mean and maximum ratio integrating the data from all the three sheep in the experiment and at different stages of learning. From this it can be seen that both the mean and maximum ratios in the right hemisphere IT are about double those in the left hemisphere. This clearly indicates that for the right hemisphere, the theta band interaction is more dominant, i. e. , the right hemisphere deals more with signals of lower frequency. In order to provide a deeper insight into this frequency story, we compute the ratios at different stages of learning. Fig. 5B shows the mean and maximum ratio in two sheep before learning, after learning and a month after learning. An additional set of data is one week after learning for sheep C only. The most noticeable change is the reduction in the interactions in the theta band (low frequencies) in the right IT which occurs after learning and is maintained subsequently. Combining Fig. 4C (right) with Fig. 5C, we see that learning in general changes the connections in the right hemisphere (increasing), however, the increasing interactions are mainly due to the enhancement of the interaction at the high frequency domain.
summary: The right temporal cortex has previously been shown to play a greater role in the discrimination of faces in both sheep and humans. In the frequency domain, analysis of the relative causal contributions of low (theta 4–8Hz) and high (gamma 30–70Hz) frequency oscillations reveals that prior to learning, theta activity is more predominant in right than in left hemisphere processing, and that learning reduces this so that high frequency oscillations gain more control. We have been able to demonstrate that the frequency of connections increases in the right hemisphere and decreases between the left and right hemispheres after learning. The results are obtained based upon a way to combine aspects of both the Granger and Dynamic Causal Models, which can be used to establish significant causal relations in both time and frequency domains and applied to local field potential recordings from multiple (64 channel) electrodes implanted in the inferotemporal cortex of both sides of the brain in sheep in order to establish changes in causal connections within and between the two hemispheres as a result of learning to discriminate visually between pairs of faces. It is anticipated that this new approach to the measurement of causality will not only help reveal how the two brain hemispheres interact, but will also be applicable to many different types of biological data where variations in both frequency and temporal domains can be measured.
keywords: mathematics
  computational biology/systems biology
  mathematics/statistics
title: A Novel Extended Granger Causal Model Approach Demonstrates Brain Hemispheric Differences during Face Recognition Learning