Granger causal time-dependent source connectivity in the somatosensory network.

Exploration of transient Granger causal interactions in neural sources of electrophysiological activities provides deeper insights into brain information processing mechanisms. However, the underlying neural patterns are confounded by time-dependent dynamics, non-stationarity and observational noise contamination. Here we investigate transient Granger causal interactions using source time-series of somatosensory evoked magnetoencephalographic (MEG) elicited by air puff stimulation of right index finger and recorded using 306-channel MEG from 21 healthy subjects. A new time-varying connectivity approach, combining renormalised partial directed coherence with state space modelling, is employed to estimate fast changing information flow among the sources. Source analysis confirmed that somatosensory evoked MEG was mainly generated from the contralateral primary somatosensory cortex (SI) and bilateral secondary somatosensory cortices (SII). Transient Granger causality shows a serial processing of somatosensory information, 1) from contralateral SI to contralateral SII, 2) from contralateral SI to ipsilateral SII, 3) from contralateral SII to contralateral SI, and 4) from contralateral SII to ipsilateral SII. These results are consistent with established anatomical connectivity between somatosensory regions and previous source modeling results, thereby providing empirical validation of the time-varying connectivity analysis. We argue that the suggested approach provides novel information regarding transient cortical dynamic connectivity, which previous approaches could not assess.

The human brain can be seen as a multi-layered complex network, stimulated by external inputs, to process affective and cognitive information, thus instructing the execution of appropriate responses1234. The exploration of stimulus-activated neural sources and their connectivity provides deep insights into the brain mechanisms of information processing within the complex brain network56. The human somatosensory system is studied by means of functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG) or electroencephalography (EEG) with tactile sensory stimulation. Unilateral sensory stimulation reliably elicits contralateral somatotopic activation of the primary (SI) and secondary (SII) somatosensory cortex789101112. Some studies reported reliable activation of contralateral SI and bilateral SII in response to unilateral sensory stimuli131415161718192021, and a few studies observed occasional ipsilateral SI activation2223242526. Despite these numerous reports few studies have investigated the quantitative relationship of the cortical activations, let alone the time-dependent Granger causal connectivity. Thus appropriate investigations remain to date rare and sometimes inconclusive. Consequently the information flow among cortices involved in the somatosensory information processing is not well understood, especially the relationship between the contralateral somatosensory cortices (both SI and SII) and ipsilateral SII927.

In the neurosciences, functional connectivity of neural networks is usually estimated by classical methods, such as correlation and coherence, based on time- or frequency domain analysis in sensor space. These functional connective methods do not provide a direction of information flow. When dealing with brain functional data sets, ultimately the algorithms should be capable of dealing with truly multivariate data accounting for direct and indirect influences. The algorithms also need to cope with various noise sources, address both linear and nonlinear systems, and provide the strengths of the directed interactions for each sampling point. Furthermore, they need to provide a rigorous statistical framework. Until now a rigorous mathematical framework comprising all of these features was still lacking. Some methods could estimate time-dependent connectivity patterns without concern for observational noise and non-stationarity2829. Other techniques that enable investigation of directional relationships on very short time scales from multiple realisations of short and transient noisy time series have been proposed3031. An overarching framework for data-based modelling, combining state space modelling323334 for autoregressive models with time-dependent autoregressive coefficients, enabling the estimation of time-resolved renormalised partial directed coherence in the frequency domain, fulfills all the above requirements. The approach quantifies the direction and strength of time-dependent network connectivity without relying on any prior assumption about the nonlinear, non-stationary and stochastic brain signals3536. The performance of this novel approach has also been validated in a simulation study35.

In this paper, we apply the time-varying Granger causal connectivity approach to MEG data from healthy humans. We use reliable unilateral sensory paradigms to explore source connectivity patterns within the human somatosensory system in response to tactile stimulation by identifying time-dependent effective connectivity from high-density MEG recordings. We hypothesize that somatosensory activation, as measured by MEG in humans, should qualitatively and quantitatively reflect the effective connectivity of the somatosensory network. The source locations and single-trial waveforms are estimated by Minimum Norm Estimation (MNE) to mitigate the field spread in MEG recordings, which complicates the interpretation of cortical connectivity patterns at the sensor level. The time-dependent effective connectivity patterns among brain sources is calculated by the new approach combining renormalized partial directed coherence with state space modeling. We also employ a rigorous statistical evaluation procedure to guarantee the applicability to measured signals. This study demonstrates the effectiveness of the time-varying Granger causal connectivity to explore the brain effective connectivity in source space in MEG/EEG fields, which provides a new way to investigate the network interplay between stimulus-evoked cortical activations.

Results

The spatial distribution of the grand averages for the healthy subjects is illustrated in Fig. 1(a). We found that the evoked response to right index finger stimulation is most prominent in left central channels consistent with activation from postcentral gyrus of the anterior parietal lobe, where primary and secondary somatosensory cortex (SI and SII) are located, and the averages over the other brain areas look like the non-phase locking activity in the background. Figure 1(b) shows the grand average somatosensory MEG waveforms in all channels for the subjects, with the scalp distribution at the latencies of several typical MEG peaks (71 ms and 201 ms). Significant time points around the peak of somatosensory were selected by visual inspection of the waveform. The grand average MEG waveform consisted of a clear positive and negative peak at 71 ms followed by a negative peak at 201 ms. The scalp distribution at 71 ms extended towards temporal regions with negative maximum at the contralateral temporal region and the central region. The scalp topography at 201 ms showed a clear maximum at the ipsilateral central region. According to the results of the source distribution among all the subjects (Table 1), we selected three sources here. The source locations are SI-l (-35 mm, -29 mm, 47 mm), SII-r (41 mm, 1 mm, -8 mm), and SII-l (-41 mm, 1 mm, -8 mm) at the group level. The fitted sources explained the scalp distribution of MEG with RV of 5.5% and 10.8 ± 4.2% at the group level and single-subject level, respectively.

Figure 1

(a) Spatial distribution of MEG averaged responses. At each sensor location, the traces illustrate signals recorded by two orthogonal gradiometers. (b) Group averaged waveforms from each channel are plotted in different colors and superimposed and scalp topographies at 71 ms and 201 ms peak latencies are displayed.

Table 1

Source distribution of right somatosensory MEG.

Source location	The number of subjects
Left primary somatosensory cortex	21
Right secondary somatosensory cortex	17
Left secondary somatosensory cortex	15
Left secondary auditory cortex	5
Right primary somatosensory cortex	5
Anterior cingulate cortex	2
Left supramarginal gyrus	2
Right premotor cortex	1
Right supramarginal gyrus	1
Colliculus	1
Thalamus	1
Left intraparietal sulcus	1
Medial prefrontal cortex	1
Right secondary auditory cortex	1
Left primary auditory cortex	1
Posterior cingulate gyrus	1

Figure 2 shows the locations and grand average time series of the representative sources (SI-l, SII-r and SII-l) of the subjects. The results indicate that 1 the sources explaining the peak at 71 ms are mainly located in contralateral SI (SI-l) and contralateral SII (SII-l), 2 the sources explaining the peak at 201 ms are mainly located in ipsilateral SII (SII-r). Post hoc comparisons by Fisher’s least significant difference (LSD) test revealed that 1 the peak latencies of SI-l are significantly shorter than the peak latencies of SII-r (P < 0.001), 2 the peak latencies of SII-l are significantly shorter than the peak latencies of SII-r (P < 0.001).

Figure 2

Representative source locations and the corresponding grand average source time courses of somatosensory MEG responses from all participants are presented.

Single-trial data of somatosensory source waveforms (SI-l, SII-l and SII-r), including 120 to 180 trials for each source of each subject, obtained from source analysis, were used in the rPDC estimation to assess the relationship among sources. The rPDC order, p, selected using Akaike’s Information Criterion (AIC), ranged from 30 to 60; the order was selected for each analysis optimally using AIC. We also projected the data in the reference interval from -100 ms to 0 ms to the representative sources. Time-varying effective connectivity patterns among somatosensory sources for the subjects, represented as time-frequency regions that have significantly increased rPDC values compared with the reference interval (P < 0.01), are summarized into the following temporally distinct groups within the post-stimulus interval around 0-320 ms (Fig. 3). Significant increases in effective connectivity are observed from SI-l to both SII-l (20-230 ms, 5-25 Hz) and SII-r (20-230 ms, 5-25 Hz), from SII-l to SI-l (90-200 ms, 5-17 Hz) and from SII-l to SII-r (220-320 ms, 5-19 Hz). The time-resolved directed network structure can be inferred from Fig. 3 and Table 2 (see also Fig. 4). The cortical information was mainly transmitted from SI-l to SII-l and to SII-r (20-230 ms, the latter one about 5 ms later than the former). Later during the post stimulus interval, the information flow was observed from SII-l to SI-l (90-200 ms) and then to SII-r (220-320 ms).

Figure 3

Time-frequency representations of rPDC as a measure of time-dependent Granger causal influences within the neural network between SI-l, SII-l and SII-r averaged across all the subjects. The regions circled by blue lines had significantly larger rPDC values than those in the reference interval from -100 ms to 0 ms (P < 0.01).

Table 2

The mean and standard deviation (STD) of rPDC between SI-l, SII-l and SII-r of right somatosensory MEG responses of identified time-frequency regions in Fig. 3.

Source	Mean		STD
from SI-l to SII-l	0.17		0.01
from SII-l to SI-l	0.18		0.01
from SI-l to SII-r	0.16		0.01
from SII-r to SI-l	0.00		0.00
from SII-l to SII-r	0.17		0.02
from SII-r to SII-l	0.00		0.00

Figure 4

Graph obtained from rPDC analysis applied to source time series located at SI-l, SII-l and SII-r for somatosensory evoked MEG.

Discussion

In this study, we applied a novel source connectivity analysis method based on causality inference to assess the time-varying effective connectivity among stimulus-elicited neural sources for healthy participants. The time-varying effective connectivity revealed that the cortical information flow of somatosensory input at early latencies was processed through a complex pattern of both feedforward and feedback interactions between contralateral SI and bilateral SII: 1) from contralateral SI to contralateral SII, 2) from contralateral SI to ipsilateral SII, 3) from contralateral SII to contralateral SI, and 4) from contralateral SII to ipsilateral SII.

Detection of causal relationships between neural processes is commonly hampered by time-dependent dynamics and observational noise contamination. Time-dependent dynamics is of particular interest in EEG/MEG studies because electrophysiological information is believed to be transferred in part via short bursts of oscillatory brain activity373839. Effective connectivity analysis based on the concept of Granger causality40 allows for the evaluation of the direction and strength of causality between neuronal activations, e.g., PDC41. However, Granger causality cannot reveal the dynamic effective connectivity for transient and markedly nonstationary neurophysiological processes42. In order to address these challenges, we employed a combination of renormalized partial directed coherence and non-linear state space modeling allowing for the detection of time-dependent and eventually continuously changing causal influences in multivariate processes without relying on prior assumptions about the nonlinear, non-stationary and stochastic brain signals36. Schelter et al. demonstrated its performance in a simulation study and proved that it combines characteristics that others lacked so far35. In the present study, we used this novel approach to explore the causal influence empirically by testing the performance of the method with somatosensory evoked MEG responses of healthy subjects.

In neuroscience studies, the event-related experiment with different stimulation paradigms is a fundamental method to activate corresponding brain areas connected by underlying neuronal networks. The activated network is characterized by the evoked data in epochs, which are mixed with field spread. The field spread often complicates the interpretation of cortical connectivity patterns at the sensor level5. Carrying out connectivity analysis in source space is the best choice to disentangle the above problem, providing direct estimation of the interaction between neuronal sources and attenuating the effect of field spread544. As expected, we observed that somatosensory sources were mainly located at contralateral SI and bilateral SII after right index finger somatosensory stimulation when applying CSST and MNE to somatosensory MEG. We also found that peaks in MEG were mainly generated from these sources as shown in Fig. 1 and Fig. 2. Here, since no significant changes were observed in the post-stimulus interval after 320 ms, only source time courses in the post-stimulus from 0 to 320 were analysed. The validity of our approach is substantiated by three independent bits of evidence. First, the fitted dipoles explained the scalp distribution of MEG with RV of 5.5% and 10.8 ± 4.2% at the group level and single-subject level, respectively. Second, as shown in Fig. 2 the sources explaining the peak at 71 ms are mainly located in SI-l and SII-l, and the sources explaining the peak at 201 ms are mainly located in SII-r. Third, the results of source modelling are consistent with a large number of previous studies2945.

Tactile events are accompanied by transient changes in mu (8-15 Hz) and beta (15-30 Hz) band oscillatory activity4647. The time–varying effective connectivity shows patterns in exactly those frequency ranges. Note that there is no restriction to a frequency or frequency band imposed by the analysis technique. We therefore considered the frequency band from 5 Hz to 30 Hz for the connectivity analysis; both, previous findings about the tactile events and our analysis suggest this as the optimal frequency band.

In MEG studies of the somatosensory cortex, the tactile pulses elicit evoked responses within 50–80 ms originating from contralateral SI48495051. A study recently reported that contralateral SI was the earliest source activated by somatosensory stimuli with peaks at 63 ± 14 ms6. In our study, the initial response to the tactile stimulus was observed at around 70 ms, which is consistent with those previous findings. We also found that contralateral SI was the earliest activated source in response to somatosensory stimuli and that the connectivity from contralateral SI represented the earliest information flow (onset at about 20 ms to contralateral SII), which is in agreement with many previous studies6752535455. Unilateral sensory stimuli are transmitted to contralateral SI from primary afferent fibers of dorsal root ganglion or trigeminal sensory neurons, also supporting our results56. Contralateral SII was activated after receiving somatosensory information from contralateral SI, which is consistent with a large number of studies691327. Our results demonstrate that ipsilateral SII was activated after receiving somatosensory information from contralateral SI (onset at 20 ms) and then from contralateral SII (onset at 220 ms)69. There is considerable anatomical and electrophysiological evidence for interhemispheric connections between the somatosensory cortices via the corpus callosum. In particular the majority of SII neurons display bilateral receptive fields565758. The activation of ipsilateral SII in our study is in accord with previous studies959. Our results indicate that functional information flow follows these anatomical pathways from contralateral SI to bilateral SII, from contralateral SII to SI, and finally from contralateral SII to ipsilateral SII in the somatosensory information processing for healthy subjects.

In conclusion, a combination of renormalized partial directed coherence and state space modelling has been demonstrated to provide a reasonable estimation of time-varying effective connectivity among neural sources on the order of milliseconds from multi-channel MEG recordings. We showed that neural sources that responded to somatosensory stimuli displayed a complex pattern of connectivity that included both feedforward and feedback information flow consistent with current views of even simple sensory responses invoking a network of activity as opposed to sequential hierarchical activation of cortical regions. Our results confirmed that this analysis approach provides somatosensory network connectivity patterns consistent with previously reported anatomical and neurophysiological measures. Our empirical validation of this novel method provides feasibility for applying this method to clinical populations to identify biomarkers of disrupted network connectivity that represent the transitions from normal state to disease state, which is crucial for understanding the mechanism of numerous diseases.

Materials and methods

Human somatosensory MEG study

The human somatosensory evoked MEG responses were collected as a part of a separate study. Prior to the study, the written informed consent was provided by each subject. The experimental protocols were approved by the Human Research Review Committee at University of New Mexico Health Sciences Center, Albuquerque, New Mexico. The methods were carried out in accordance with the approved guidelines. The MEG data from twenty-one healthy subjects (12 male, 9 female), aged from 12 to 20 years (15.2 ± 2.6, age ± SD), were collected in a magnetically shielded room (Vacuumschmelze – Ak3B) at the Mind Research Network in Albuquerque using a 306-channel whole-head MEG system (Elekta Neuromag) with a sampling rate of 1000 Hz and a band pass filter between 0.1-330 Hz to avoid drifts and aliasing. Prior to the experiment, fiducial points (left and right preauricular points and nasion) and head shape data were collected and checked by the Polaris system. Participants sat upright during the task and were monitored at all times by an audio and video link between the magnetically shielded room and control room. The 306-sensor MEG system measured the magnetic field distribution around the whole head of the seated subjects. The system was fully equipped with real-time motion correction and artifact rejection software. As long as the subject’s head remained within the MEG helmet, movement could be corrected to an optimal head location using the movement correction algorithm provided with the Neuromag software package. This motion corrected dataset could then be further analysed in our study without concern of subjects’ movement across epochs. The stimuli were presented while the subjects were quietly sitting with the head situated in the MEG helmet. The somatosensory stimuli were generated by allowing a short pulse of compressed air to fill an air bladder that was attached to the subjects’ right index finger. The air pulse was controlled by the Presentation software and a compressed air regulator (located outside the shielded room). The air puff stimulus lasted 50 ms with an ISI of 1.0–1.4 seconds. The air bladder pressure applied to the index finger was monitored and recorded simultaneously with MEG collection for offline analysis and interpretation of the results. We collected 120 to 180 trials for the right index finger of each subject.

The 60 Hz powerline noise was removed60. Raw MEG data were filtered for noise sources such as eye blinks and excessive movement, then corrected for head motion using the Neuromag MaxFilter software61. Heartbeat artifacts were then removed by projecting electrocardiogram (ECG) data from MEG sensor waveforms using the signal-space projection (SSP) method62. The data for each stimulus condition were obtained from 100 ms prior to the onset of the stimulus to 600 ms following stimulus onset. The data were baseline-corrected and subjected to a 50 Hz low-pass filter during signal processing.

Structural MRI data acquisition

Structural MRIs were obtained for use in mapping source locations from all the subjects. Sagittal T1-weighted anatomical MR images were collected using a Siemens TIM Trio 3 Tesla MRI system with a multiecho 3D MPRAGE sequence [TR/TE/TI = 2530/1.64, 3.5, 5.36, 7.22, 9.08/1200 ms, flip ange = 7°, field of view (FOV) = 256 mm*256 mm, matrix = 256*256, 1 mm thick slice, 192 slices, GRAPPA acceleration factor = 2].

Source time courses calculation

Sources were localised for each subject using cortical-start spatiotemporal (CSST) multidipole analysis with integrated Multiple Signal Classification (MUSIC)63. CSST is an objective multidipole, multistart procedure in which initial dipole locations are randomly selected from a predefined cortical volume and a nonlinear simplex search is performed for each of these initial configurations. Initial dipole locations were selected from within a predefined head volume, which was defined by a subsample of points taken from within the cortical volume, as determined by coregistered structural MRI. The error is minimized using a reduced chi-square statistic to obtain a final set of dipole configurations which most fully explain the data63. CSST source localization was calculated using 3,4,5,6, and 7 dipole models, based on the averaged responses occurring between 0 ms and 320 ms after the onset of the stimulus for responses to somatosensory stimuli. A shorter time window after stimuli was used for dipole modeling of somatosensory responses to increase the power in order to detect somatosensory evoked sources. The Nelder-Meade minimization procedure was carried out 1500 to 8000 times, depending on the number of dipoles in the model, to help to ensure that the procedure would reach a global minimum. The dipole model that best explained the data was selected for source time course analysis.

Following the selection of the optimal source model, the single-trial waveforms of each dipolar source were calculated within a realistic head model with the minimum norm estimate (MNE) software64. The inverse solution yielded estimates of continuous time series of cortical currents. For each patient, the realistic cortical surface and three layers (inner skull, outer skull and skin) were reconstructed from the anatomical MRI images using the Freesurfer software (Compumedics, Charlotte, NC). The boundary element model (BEM) was then constructed with the reconstructed surfaces. The co-registration of MEG and MRI images was achieved by matching the recorded positions of three fiducial points (nasion, left and right preauricular points) with the locations of these points from the MRI images. The lead field matrix relating MEG sensors to the cortical distributed dipoles was computed with the BEM model using MNE. The dipole model, cortical surface and lead field matrix were then used in the MNE software to extract the single-trial timecourses of sources. MNE was applied before time-frequency decomposition here.

Renormalized partial directed coherence combined with state space modelling

Here we briefly introduce the time-dependent Granger causal connectivity analysis method employed in this study. A time continuous multivariate dynamical process Z(t) can only be observed as a multivariate, time discrete sampled signal35

where g(⋅) denotes the observation function with parameter set v; η(t) is a Gaussian distributed independent measurement noise with a given variance. Assuming a linear observation function, we obtain the following model:

for some appropriately chosen variances and that are optimally determined in the estimation process and where C represents the linear observation matrix. A reasonable assumption is that the parameter matrix A(i) should change more slowly than the (stochastic) dynamics itself. The model is then augmented to the over-arching state space model as following:

The a(i) are the matrix entries of A(i) rearranged into a vector. The causal influences can be represented as directed edges in a network, in which the nodes represent the processes. Thus the matrix A(i) contains the interactions between the components of the original process Z(t). The information about the network structure is also contained in this matrix. Since we do not make any assumption about the origin of Z(t), it can model the sensor space as well as the source space equally well.

In networks, influences with a certain delay are typically relevant. This can be accounted for in state space modelling by including previous time steps,

up to a maximum time lag p. This maximum p can be determined relying on a priori knowledge or based on model selection criteria, such as Akaike’s Information Criterion (AIC) used in our study3536. The higher order process can be rewritten as a first order process by introducing

The matrix assumes the specific structure

In the above form, the trajectory of the Z(i) and the a(i) can be determined purely based on the observations Y(i). The dual Kalman filters rely on the knowledge of the exact dynamics, i.e., a precise knowledge of C, , , and estimated by expectation maximisation algorithm3336. The expectation maximization algorithm is an iterative algorithm that converges in the sense of maximum likelihood to the best estimator of the underlying dynamical process Z(i) and the parameters a(i)65.

The renormalized partial directed coherence (rPDC) is derived as follows66:

This is a frequency domain measure for Granger–causality66 that quantifies the direction and the strengths of network connections. X(ω) = [R(FT(A)), I(FT(A))]′ with R the real and I the imaginary parts, FT the Fourier transform, and N the number of data points analysed. The normalization by (V(ω)) − 1 is given by the inverse of N-times the covariance matrix

of the estimates of X. The covariance matrix of the estimated parameters is determined in the expectation maximization algorithm. The mathematical details were demonstrated in our previous work3536.

Statistical analysis

In order to test whether the rPDC within the stimulation period was significantly different from the reference period, two populations were extracted for comparison: the collection of rPDC for the somatosensory stimulation period λ(n,f) at each time-frequency point (n,f) from 0 to 320 ms with a sampling rate of 1000 Hz and the reference population λ(τ,f|τ∈reference period) from -100 to 0 ms for 21 subjects. The null hypothesis is that there is no difference between the means of the two populations.

Author Contributions

B.S., C.G., J.W. and D.L. formulated the idea of the paper and supervised the research. J.S. and T.Z. provided the MEG data samples, and made useful suggestions for the study. B.C. preprocessed the data. L.G. processed and analyzed the data. L.S., B.S. and B.C. helped the data analysis. B.S. and L.S. constructed the computational method. L.G. and L.S. wrote the manuscript. C.G., B.S. and J.S. revised the manuscript.

Additional Information

How to cite this article: Gao, L. et al. Granger causal time-dependent source connectivity in the somatosensory network. Sci. Rep. 5, 10399; doi: 10.1038/srep10399 (2015).