Beamformer-based spatiotemporal imaging of linearly-related source components using electromagnetic neural signals

NeuroImage 114 (2015) 1–17

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Beamformer-based spatiotemporal imaging of linearly-related source components using electromagnetic neural signals Hui-Ling Chan a,b, Li-Fen Chen c,d, I-Tzu Chen a, Yong-Sheng Chen a,e,⁎ a

Department of Computer Science, National Chiao Tung University, Hsinchu, Taiwan Department of Neurology and Neurosurgery, Montreal Neurological Institute, McGill University, Montreal, Canada Institute of Brain Science, National Yang-Ming University, Taipei, Taiwan d Integrated Brain Research Unit, Department of Medical Research, Taipei Veterans General Hospital, Taipei, Taiwan e Institute of Biomedical Engineering, National Chiao Tung University, Hsinchu, Taiwan b c

a r t i c l e

i n f o

Article history: Received 5 January 2013 Accepted 14 March 2015 Available online 21 March 2015 Keywords: Beamformer Functional connectivity Temporal correlation Magnetoencephalography

a b s t r a c t Functional connectivity calculated using multiple channels of electromagnetic brain signals is often over- or underestimated due to volume conduction or field spread. Considering connectivity measures, coherence is suitable for the detection of rhythmic synchronization, whereas temporal correlation is appropriate for transient synchronization. This paper presents a beamformer-based imaging method, called spatiotemporal imaging of linearly-related source component (SILSC), which is capable of estimating connectivity at the cortical level by extracting the source component with the maximum temporal correlation between the activity of each targeted region and a reference signal. The spatiotemporal correlation dynamics can be obtained by applying SILSC at every brain region and with various time latencies. The results of six simulation studies demonstrated that SILSC is sensitive to detect the source activity correlated to the specified reference signal and is accurate and robust to noise in terms of source localization. In a facial expression imitation experiment, the correlation dynamics estimated by SILSC revealed the regions with mirror properties and the regions involved in motor control network when performing the imitation and execution tasks, respectively, with the left inferior frontal gyrus specified as the reference region. © 2015 Elsevier Inc. All rights reserved.

Introduction Elucidating the mechanisms and pathways involved in neural communication is increasingly important for understanding how information is processed within the brain (Averbeck and Lee, 2004) as well as for mapping a comprehensive functional connectome (Biswal et al., 2010; Smith, 2012). Previous findings have suggested that temporal correlation of neural activity may be an indication of communication and information flow between cortical neurons or neural assemblies (Salinas and Sejnowski, 2001; Singer and Gray, 1995; von der Malsburg, 1999). Neurons are capable of synchronizing their firings on a time scale of milliseconds to fulfill sensory-motor, perceptual, and cognitive functions (Azouz and Gray, 2000; Engel and Singer, 2001; Engel et al., 1999; König and Engel, 1995; Singer and Gray, 1995). The connectivity of neural networks presents rapid variations. Hence, it is beneficial to estimate connectivity by using magnetoencephalography (MEG) or electroencephalography (EEG) recordings with a high degree of temporal resolution (Schoffelen and Gross, 2009). This study ⁎ Corresponding author at: Department of Computer Science, National Chiao Tung University, 1001 University Road, Hsinchu, Taiwan. Fax: +886 3 5721490. E-mail addresses: [email protected] (H.-L. Chan), [email protected] (L.-F. Chen), [email protected] (I.-T. Chen), [email protected] (Y.-S. Chen).

http://dx.doi.org/10.1016/j.neuroimage.2015.03.038 1053-8119/© 2015 Elsevier Inc. All rights reserved.

proposes a temporal correlation-based source connectivity estimation method, which can be used to quantify the interdependency between a reference signal and source activity calculated from electromagnetic signals. Numerous techniques have been developed to estimate connections among brain regions using MEG or EEG signals. Among them, correlation coefficient (Brazier and Casby, 1952) and coherence (Nunez et al., 1997) are the most commonly used measures to evaluate linear associations in the time and frequency domains, respectively. In measuring the degree of synchronization, phase synchronization methods can be used to estimate the relationship of oscillation phases between two signals (Hindriks et al., 2011; Lachaux et al., 1999; Mormann et al., 2000; Tass et al., 1998). Moreover, generalized synchronization methods measure the level of synchronization and provide the direction of information flow (Arnhold et al., 1999; Rulkov et al., 1995). In theory, these measures of functional connectivity are all related to correlation (Marrelec et al., 2005). Quiroga et al. (2002) showed similar results of connectivity estimation using the above-mentioned measures. In simulation experiments, Silfverhuth et al. (2012) demonstrated that the correlation coefficient technique is sensitive to direct causal connections. Ansari-Asl et al. (2006) and Wendling et al. (2009) further demonstrated that, with regard to coupling model parameters, both correlation coefficient and coherence techniques achieved sensitivity equal or

2

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

superior to other methods based on phase synchronization and general synchronization. Finally, in an EEG study by Guevara and Corsi-Cabrera (1996), the results obtained by coherence were highly similar to those obtained by correlation coefficient. Coherence methods quantify the synchronization between two signals by comparing their spectral power dynamics, which makes these techniques suitable for detecting rhythmic synchronization (Nunez et al., 1997). Compared to coherence methods, correlation methods are better suited to detect the synchronization caused by time-locked responses (Kujala et al., 2007) or transient responses, such as P300 waves. Moreover, correlation methods achieve greater sensitivity than coherence methods when information related to time latency is provided beforehand (Ansari-Asl et al., 2005). To estimate latency, Gevins and Bressler (1988) introduced a measure called eventrelated covariance (ERC), which is the multi-lag cross-covariance between the time series of two channels. Latency is then determined according to the time lag with the maximum covariance. This measure has been applied to predict performance accuracy (Gevins et al., 1987), to study the function of gamma-band brain waves (Gevins et al., 1995; Menon et al., 1996), and to investigate propagation in both visuomotor tasks (Gevins et al., 1989a,b) and working memory (Gevins and Cutillo, 1993). ERC was only applied to estimate connectivity or interdependency between a pair of electrodes. However, at the sensor level, connectivity may be over- or underestimated due to the effects of volume conduction on EEG or field spread on MEG (Hillebrand et al., 2012; Nunez and Srinivasan, 2006; Winter et al., 2007). These effects often result in artificial synchrony by a single cortical source contributing to multiple channels simultaneously (Palva and Palva, 2012). Schoffelen and Gross (2009) reported an example of overestimated connectivity using simulation data generated from uncorrelated sources. In a simulation study by Grasman et al. (2004), an estimated connectivity map displayed problematic distribution when two synchronized sources existed within the brain. In fact, data of a single MEG or EEG channel is a mixture of numerous sources with varying degrees of correlation among them. As a result, interpreting complex data related to MEG or EEG sensor-level connectivity poses significant challenges. Connectivity estimation methods calculated at the cortical level could reduce the artificial synchrony at the sensor level caused by volume conduction and field spread (Palva and Palva, 2012). These methods first estimate cortical source activity and then calculate source connectivity by determining the level of interdependency among the estimated sources. For example, BESA source coherence can be used to calculate coherence between dipole sources estimated by using the dipole fitting method (Hoechstetter et al., 2004). Brookes et al. (2011a) estimated the functional connectivity by calculating the envelope correlation or coherence between two brain source signals obtained by beamfoming spatial filters. In the study conducted by Hipp et al. (2012), the artificial synchrony between two sources was reduced by orthogonalizing signals prior to the calculation of coherence. Beamforming methods can estimate brain source activity with better spatial resolution than other source estimation methods (Dalal et al., 2008; Darvas et al., 2004; Sekihara et al., 2005). For beamforming methods, the neural activity is modeled by an equivalent current dipole and is estimated by a spatial filter for each position in the brain source space. Vector beamformer, such as linearly constrained minimum variance (LCMV) beamformer (Van Veen et al., 1997), calculates the components of dipole source activity in three orthogonal directions. For scalar beamformer, the dipole orientation has to be determined to obtain the beamforming spatial filter either by an exhaustive search, as in the synthetic aperture magnetometry (SAM) method (Robinson and Vrba, 1998), or by an analytical solution, as in the maximum contrast beamformer (MCB) method (Chen et al., 2006). In the experiments conducted by Vrba and Robinson (2000) and Chen et al. (2006), scalar

beamforming methods can provide a better sensitivity and spatial resolution than LCMV beamformer. Estimation of source activity for specific brain regions may include the contribution of multiple neural populations (Gross and Ioannides, 1999). In other words, each neural region contains multiple sources with various orientations and temporal profiles. Therefore, activity measured in individual brain regions may actually comprise multiple components and uncorrelated components may result in underestimated source connectivity. To decrease the influence of uncorrelated components in connectivity estimation, Gross et al. (2001) proposed a vector beamforming method, dynamic imaging of coherent sources (DICS), to calculate the source component with dominant coherence. DICS has been applied to investigate the functional network during reading (Kujala et al., 2007) and to estimate the connection density (Kujala et al., 2008). This paper proposes a beamformer-based imaging method, called spatiotemporal imaging of linearly-related source component (SILSC), which is capable of estimating whole-brain functional connectivity with low artificial synchrony. For each targeted position, SILSC determines a spatial filter, which extracts the source component with the maximum temporal correlation to a given reference signal. The orientation of the dipole at the targeted position is accurately calculated in a closed-form manner. By calculating the correlation value between the reference signal and each cortical region within the brain, SILSC produces a correlation map for further identifying regions significantly correlated to the reference. Following the concept of ERC in determining time latency, SILSC uses a sliding time window to estimate the propagation latency between different regions of the brain. Experiments of simulation data and an MEG experiment involving the imitation of facial expressions demonstrated the feasibility of the proposed method. Methods and materials Beamformer-based correlation imaging The correlation coefficient, Rθ, between the reference signal a(t) and the source activity sθ(t) is defined as follows: Efðsθ ðt Þ−Efsθ ðt ÞgÞðaðt Þ−Efaðt ÞgÞg Rθ ¼
1
1 ; E ðsθ ðt Þ−Efsθ ðt ÞgÞ2 2 E ðaðt Þ−Efaðt ÞgÞ2 2

ð1Þ

where E{⋅} denotes the expectation value and the parameters θ = {r, q} represent the dipole located in position r ∈ ℝ3 with orientation q ∈ ℝ3. Source activity sθ(t) can be estimated up to a scale factor λθ by applying a spatial filter wθ ∈ ℝN on the MEG measurements m(t) ∈ ℝN: T

sθ ðt Þ ¼ λθ wθ mðt Þ:

ð2Þ

In this paper, “T” indicates the transpose of a matrix or vector and N is the number of MEG channels. By substituting Eq. (2) into Eq. (1) and canceling the scale factor λθ, the correlation Rθ can be calculated as follows: n n o o T T E wθ mðt Þ−E wθ mðt Þ ðaðt Þ−Efaðt ÞgÞ Rθ ¼ n
2 o12
1 : E ðaðt Þ−Efaðt ÞgÞ2 2 E wTθ mðt Þ−E wTθ mðt Þ

ð3Þ

DICS DICS is a vector beamforming method which can estimate the coherence or correlation between the reference signal and the source activity in each brain region (Gross et al., 2001, 2002). In vector beamforming, three spatial filters are computed for three

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

3

orthogonal directions [qx, qy, qz] at each location r. The solution of i  h spatial filters Wr ¼ wθx ; wθy ; wθz ¼ wr;q x ; wr;qy ; wr;q z ∈ ℝN3 can

minimum norm and minimum variance criteria, the beamformer wθ can be derived as below (Chen et al., 2006):

be obtained by applying the unit-gain constraint and the minimum variance criterion, as shown below (Van Veen et al., 1997; Gross et al., 2001):

^θ ¼ w

  ^ ¼ ðC þ αIÞ−1 L L T ðC þ αIÞ−1 L −1 ¼ A B−1 : W r s r s r r r r

ð4Þ

In Eq. (4), the covariance matrix Cs ∈ ℝN × N of MEG recordings ms(t) ∈ ℝN is calculated by: n o T Cs ¼ E ðms ðt Þ−Efms ðt ÞgÞðms ðt Þ−Efms ðt ÞgÞ ;

ð5Þ

where ms(t) comprises the segment of MEG measurements during a time h i window Ts, I ∈ ℝN × N is the identity matrix, Lr ¼ lθx ; lθy ; lθz ∈ℝN3 is

ðCs þ αIÞ−1 Lr q Aq ¼ r : ðLr qÞT ðCs þ αIÞ−1 Lr q qT Br q

ð11Þ

In this work, the value of regularization parameter α was set to be 0.003 and 0.03 times the maximum eigenvalue of the covariance matrix Cs for all simulation experiments other than Simulation 5 and the facial expression imitation experiment, respectively. Accuracy of the dipole orientation q strongly influences the estimation of the spatial filter wθ (Hillebrand and Barnes, 2003). Rather than performing an exhaustive search, SILSC provides an analytical solution for dipole orientation q by maximizing the square of correlation between the source signal sθ(t) and the reference signal a(t). After replacing the measurements m(t) by mc(t) in Eq. (2) and then substituting Eq. (2) into Eq. (1), the square of correlation R2θ can be calculated as follows:

the gain matrix for the dipole located at position r, and the two matrices Ar = (Cs + α I)−1L r ∈ ℝN × 3 and Br = L Tr Ar ∈ ℝ3 × 3 depend only on dipole position r. Regularization parameter α is used to compromise between the minimum norm and the minimum variance criteria. Then the covariance csa,r between the reference signal a(t) and the source in the dominant direction of the dipole at location r can be calculated by:

After substituting Eq. (11) into Eq. (12), the dipole orientation q can be calculated by maximizing the square of correlation R 2θ as follows:

n o ^ Tc ; csa;r ¼ σ 1 W r ca

^ ¼ arg max R2θ q

ð6Þ

where σ1{·} indicates the largest singular value of the expression in the braces. In Eq. (6), the vector cca ∈ ℝN contains the covariances between the reference signal a(t) and the MEG recordings mc(t) ∈ ℝN of all channels: cca ¼ Efðmc ðt Þ−Efmc ðt ÞgÞðaðt Þ−Efaðt ÞgÞg;

ð7Þ

where mc(t) comprises the segment of MEG measurements during a time window Tc. The variance cs,r of the source in the dominant direction of the dipole at location r can be calculated by: n o ^ TC W ^ : cs;r ¼ σ 1 W r r c

ð8Þ

In this equation, the covariance matrix Cc ∈ ℝN × N of MEG recordings mc(t) ∈ ℝN is defined by: n o T Cc ¼ E ðmc ðt Þ−Efmc ðt ÞgÞðmc ðt Þ−Efmc ðt ÞgÞ :

ð9Þ

After substituting Eqs. (6) and (8) into Eq. (3), the correlation Rθ between the reference signal and the source activity located at r can be obtained by n o ^ Tc σ1 W r ca Rθ ¼ n n oo1
1 : 2 ^ TC W ^ E ðaðt Þ−Efaðt ÞgÞ2 2 σ1 W r r c

ð10Þ

Note that the covariance matrix can be replaced by the cross spectral density matrix for the estimation of the coherence, instead of the correlation value, and the dipole orientation can be decomposed into two orthogonal components on the tangential plane (Gross et al., 2001). SILSC DICS can compute the correlation value by using the covariance between the reference signal and the source activity in two or three orthogonal directions of a dipole moment. In contrast, the proposed SILSC method is an inverse algorithm which can provide an analytical solution to determine the dipole source orientation and can estimate the degree of correlation between the reference signal and the source activity in each brain region. By applying the unit-gain constraint as well as the

T

2

Rθ ¼ 

wTθ Cc wθ

T

w c c w  θ
ca ca θ : E ðaðt Þ−Efaðt ÞgÞ2

q

¼ arg max q

qT Pr q ; qT Q r q

ð12Þ

ð13Þ

where Pr ∈ ℝ3× 3 equals ATr ccacTcaAr and Q r ∈ ℝ3× 3 equals ATr CcAr. In the derivation of Eq. (13), the term E{(a(t) − E{a(t)})2} is a scalar and can thus be eliminated. The solution to Eq. (13) is the eigenvector corresponding to the maximum eigenvalue of the matrix Q −1 r Pr (Chong and Zak, 2004). Because Q r and Q −1 r Pr are both 3 × 3 matrices, the matrix inverse of Q r and the eigenproblem of Q −1 r Pr can be solved in a closed-form manner. In the proposed method, the activity estimated by the beamformer ŵθ is the source component with the maximum correlation to the reference signal. This source component is a part of overall source activity at targeted location r. The time window Ts for the segment ms(t) of MEG recordings should contain the event-related activity, whereas the time window Tc is the duration to be determined whether the segment mc(t) contains source activity correlated to the reference signal. Therefore, the length of Tc must be the same as that of the reference signal for calculating the correlation. The time windows Ts can have a longer period but it should cover the duration of Tc. Reference signals To probe the targeted neural network within the brain, a reference signal related to the network should be specified or provided beforehand. The reference signal can be the source activity of a specified reference region or an independent component decomposed from MEG measurements. The source activity can be estimated from MEG measurements through inverse methods, such as minimum norm estimates (Hämäläinen and Ilmoniemi, 1994) and conventional beamforming methods (Van Veen et al., 1997; Chen et al., 2006). In this study, MCB was applied on MEG measurements to estimate the temporal profiles of the sources at reference regions (Chen et al., 2006). The reference region can be manually specified according to physiological a priori information or selected as the area with the strongest source activity estimated from the MEG data (Gross et al., 2001). Independent components are decomposed from the measurements m(t) by independent component analysis (ICA) under the assumption that the sources contributing to the signal m(t) are statistically independent (Bell and Sejnowski, 1995; Comon, 1994). In other words, m(t) = Xu(t), where X ∈ ℝN × J represents the mixing matrix and each element ui(t) in

4

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

u(t) = [u1(t), u2(t), …, uJ(t)]T ∈ ℝJ represents one of J independent components. When one of the independent components is specified as the reference signal, SILSC calculates the correlation coefficient between the brain activity in each targeted region and the selected independent component. Correlation dynamics Compared to the reference signal, the correlated source activity has similar temporal profile but may have a time shift because of the propagation delay in the neural network. For each instance of latency, we can apply SILSC to calculate a correlation map of the whole brain. Because correlated activity at different regions may have different latencies, SILSC establishes spatiotemporal correlation dynamics by sliding the time window Tc over the MEG measurements m(t) and produces a series of correlation maps, as its flowchart shown in Fig. 1. The spatiotemporal correlation maps can then be used for further interpretation or connectivity analysis. In this study, the time window Ts was set to be the same with Tc and it followed the sliding process of Tc. For each time interval, the dipole orientation and spatial filter were recomputed for the correlation estimation at each brain region. Reduction of spurious connectivity SILSC reduces the effect of field spread by reconstructing the sources linearly related to the reference signal. However, the source leakage of

spatial filters may cause spurious connectivity. That is, a single dipole source may pass its source activity, with different magnitudes, to the beamforming spatial filters estimated at nearby positions around the dipole source. Spurious connectivity may appear with high correlation values among neighboring regions resulted from the common source activity. In this work, we slightly modified the method proposed by Brookes et al. (2011a) and applied it to detect the brain regions with spurious connectivity. For each targeted position, SILSC generates a set of simulated MEG data originated from two uncorrelated dipole sources, one located at the reference region and the other at the targeted position. The orientations and the powers of these two dipole sources are set to be the same as those calculated from the original MEG measurements using the method described in the SILSC section. For the targeted position and reference region, moreover, the spatial filters estimated from the original MEG measurements are applied to the simulated data set to calculate the source waveforms at these two locations. Then the correlation value between these two waveforms is calculated. Because the underlying two dipole sources are uncorrelated, the correlation value of simulation data is purely originated from the source leakage of the spatial filters and is then used as the reference of thresholding to detect spurious connectivity. For each position in the brain, the correlation value estimated for a time window of MEG measurements is considered only if it is larger than the corresponding correlation threshold calculated from the simulation data. The correlation values passing the

Selection of a reference signal Select a position as the reference region

Apply ICA to decompose MEG signals into independent components Select one of the independent components

or

Reference region

Calculate the source signal at the reference region Reference signal

Specify a segment of the component as the reference signal

Specify a segment of the source activity as the reference signal time Time window 1

SILSC For each time window For each brain region Time window 2

Estimate the dipole orientation by Eq. (13) Calculate the beamformer by Eq. (11) Calculate the correlation between the filtered source activity and the reference signal by Eq. (12)

Time window 3

Correlation dynamics

time MEG recordings

Fig. 1. Flowchart for mapping the spatiotemporal correlation dynamics between the reference signal and the source activity estimated from electromagnetic signals. Step 1 involves the selection of a reference signal from the estimated source activity in the reference region or from independent components obtained from ICA. A sliding time window on MEG measurements enables SILSC to calculate the correlation dynamics between the reference signals and the brain activity in each neural region.

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

thresholding are then examined using one-sample t-test and p-values are adjusted with Bonferroni correction (p b 0.001). In this work, the above-mentioned procedures were applied in the experiments of Simulation 6 and facial expression imitation study. Contrast-weighted correlation Because the correlation coefficient Rθ is a normalized measure, SILSC may obtain high correlation values for those regions having source activity with low amplitudes but similar profiles as the reference signal. As in Chen et al. (2006), SILSC calculates the contrast value Fθ for each position in the brain, which is the ratio of the source magnitude during the active state over that during the control state: T

Fθ ¼

wθ Cc wθ wTθ Cn wθ

!1 2

;

ð14Þ

where Cn is the covariance matrix computed from the MEG recordings during the control state. Then the average of the contrast values across the whole brain and all the time windows is calculated as a reference of thresholding. Only the regions having contrast values larger than this threshold are considered for further processing. To reveal task-related sources having strong amplitudes, the estimated correlation value Rθ is multiplied by the contrast value Fθ and we obtain the contrast-weighted correlation at the targeted position: e ¼ jR j  F : R θ θ θ

ð15Þ

Parcellation To investigate the correlation dynamics of each brain region, the parcellation based on the automated anatomical labeling (AAL) atlas was performed in the facial expression imitation experiment. The individual AAL atlas was obtained by spatially normalizing the CH2 template (Holmes et al., 1998) to the magnetic resonance imaging (MRI) data of the subject and then applying the estimated transformation matrix to the standard AAL atlas. This procedure was performed by using the statistical parametric mapping software package (Ashburner, 2000). For each of the 90 cerebral regions in the individual AAL atlas, we calculated the correlation value Rθ using Eq. (12) and the contrast value Fθ using Eq. (14) at each voxel position in the region. A voxel was removed from the parcellation process if its correlation value was lower than the threshold for spurious connectivity detection, as described in the Reduction of spurious connectivity section, or if its contrast value was lower than the contrast threshold, as described in the Contrast-weighted correlation section. The regional contrastweighted correlation RAAL was then computed as below: s 1X e ; R V v¼1 θv

V

RAAL ¼

ð16Þ

where V is the total number of voxels in the targeted AAL region and Vs is the number of voxels in this region passing both the correlation and contrast thresholding. Experiments To evaluate the performance of SILSC, we applied it to analyze MEG data in six simulation studies and an experiment involving the imitation of facial expressions. Materials A whole-head MEG system (Neuromag Vectorview 306, Neuromag Ltd., Helsinki, Finland) was used to record the magnetic fields produced

5

by electrophysiological activity in the brain. In the six simulation studies, the forward model used to generate the MEG recordings was also based on this MEG system. The data sampling rate was set at 1 kHz for both the simulation studies and the facial expression imitation experiment. MRI data of the participating subjects were acquired on a 1.5 T GE scanner at the Taipei Veterans General Hospital (3D-FSPGR, TR = 8.67 ms, TE = 1.86 ms, matrix size = 256 × 256 × 128, voxel size = 1.02 × 1.02 × 1.5 mm3). In this study, the preprocessing procedure for MEG data comprised five steps: artifact rejection to remove epochs with magnitudes exceeding 6000 fT/cm in any gradiometer channel, ocular artifact rejection to remove epochs with magnitudes exceeding 0.2 mV in any electrooculography channel, signal space projection, bandpass filtering (2–30 Hz), and baseline correction to subtract the average over a prestimulus interval from −200 to 0 ms. Simulations Given K sources with known dipole parameters θk = {rk, qk} and temporal waveforms sθk ðt Þ, k = 1, …, K, simulated MEG data m(t) was generated using the following forward model: mðt Þ ¼

K X

lθk sθk ðt Þ þ nðt Þ;

ð17Þ

k¼1

where l θk ¼ L rk qk is the lead field vector calculated by multiplying the gain matrix Lrk by the dipole orientation qk and n(t) comprises background activity and sensor noise. To simulate the background activity of the brain, 3000 dipole sources were randomly deployed within a sphere of radius 70 mm centered at [0, 0, 40 mm] in the head coordinate system. The temporal activity of the additive background sources was drawn from a zero-mean Gaussian distribution with standard deviation of 5 nAm. The gain matrix for each dipole position was estimated using the spherical head model presented by Mosher et al. (1999) with the forward calculation method derived by Sarvas (1987). In Simulation 1, the two procedures for reference signal selection were examined. In this simulation, four dipole sources {s1, s2, s3, s4} with amplitude-modulated sinusoidal waveforms were placed in three positions {r1, r2, r3} within the brain, as shown in Fig. 2(a). The peak amplitudes of temporal profiles for the four dipoles were set to 80 nAm. Note that the temporal profile of s2 is the same as that of s3 with a 350-ms delay. The orientations of the four dipoles were randomly selected on the tangential plane and the MEG data obtained from the four dipole sources was computed using Eq. (17). For the MEG data in this simulation, we applied MCB to calculate the F-statistic maps showing significance levels of brain activity for all brain regions (Chen et al., 2006). The dipole orientation was estimated through the maximization of the F-statistic value, representing the contrast of source powers between the active and control states. In this simulation, the MEG data during 200–400 ms were used to calculate the source power of active state and the empty room recordings were used to calculate the source power of control state. The peak position with the highest F-statistic value was selected as the reference region. According to the temporal profile of reference region, the source activity power was large from 200 to 400 ms. Therefore, the source activity during this time period was specified as the reference signal used in SILSC to calculate the correlation dynamics. We also applied FastICA (Hyvärinen, 1999; Hyvärinen and Oja, 2000) to the simulated MEG data and calculated the correlation maps for each of the selected components using SILSC. In Simulation 2, the performance of SILSC was evaluated by 700 sets of MEG data (each with ten trials) with respect to the signal-to-noise ratio (SNR) of data ranging from −10 to −1 decibels (dB). In each data set, one dipole source sθg was randomly located on a spherical surface with a radius of 70 mm centered at [0, 0, 40 mm] in the head coordinate system. The orientation of each dipole source was randomly aligned on

6

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

(a) Simulation setup

(c) Correlation maps referenced with source activity at r2 Correlation maps

r1 0.99

… During 200-400 ms

s1 0.8

0.99

During 550-750 ms 0.8

Estimated source activity at the reference region r2

amplitude (10 nAm)

r2 Reference signal s2

6 4 2 0 -2 -4 -6 -8

0 100 200 300 400 500 600 700 800

time (ms)

r3 s3

Generate MEG data by Eq. (17)

s4

amplitude (100 fT/cm)

(b) MEG data

…

4

Estimate source activity and choose the reference region with peak activity

2 0 -2 -4 -6 0 100 200 300 400 500 600 700 800

time (ms)

Decompose MEG data into independent components using ICA

(d) Correlation maps referenced with independent components Correlation maps

Correlation maps

Independent components u1 u2

0.99

0.75

u1 as the reference signal

0.99

SILSC

u3 0.85

u2 as the reference signal

u4

…

0.99

0.8

u3 as the reference signal

0.99

u114 100 200 300 400 500 600 700 800

time (ms)

0.8

u4 as the reference signal

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

the tangential plane. During the time interval of 172–428 ms, the source activity was a mixture of two sinusoidal waveforms: sθg ðt Þ ¼ a½ð1−βÞsin 2π f 1 t þ β sin 2π f 2 t ;

ð18Þ

where the parameter set {a, β, f1, f2} contains the amplitude a (nAm), the constitution ratio β, and two frequencies f1 and f2 (Hz). In Simulation 1, these parameters were set to be {100, 0.3, 17, 5}. The SNR for the dipole source sθg was calculated in dB as follows (Goldenholz et al., 2009): SNR ¼ 10 log10

Q X N l 2θg ;p Ps X ; Q N q¼1 p¼1 P n ðp; qÞ

ð19Þ

where P s represents the mean power of the activity for source sθg , Q represents the number of trials, lθg ;p represents the p-th element of the lead field vector lθg , and P n ðp; qÞ represents the mean power of noise term in the p-th channel during the q-th trial. The reference signal for the correlation calculation had the same waveform as the ground truth source during the period of 172–428 ms. To estimate the time latency between the reference signal and the source activity, the center of a 256-ms sliding window was moved from 296 to 304 ms in 1-ms steps to segment MEG recordings mc(t) for the correlation calculation. The time latency (TE) was determined by the time window having the maximum correlation value. In Simulation 3, the performance of SILSC was examined by comparing the results estimated from the data sets generated in Simulation 2 with those estimated from the corresponding surrogate data. For each MEG data set in Simulation 2, 10 sets of Fourier-transformed surrogate data were generated by randomizing signal phases (Theiler et al., 1992). The reference signal is the same as that used in Simulation 2. The time window for correlation estimation was fixed at the interval of 172–428 ms. The spatial filter, dipole orientation, correlation value, contrast value, and contrast-weighted correlation of each voxel position were computed for each set of surrogate data by using SILSC. In Simulation 4, the performance of SILSC with respect to the underlying correlation and SNR was evaluated by the 7,000 sets of MEG surrogate data used in Simulation 3. In Simulation 4, the correlation errors and localization errors were evaluated from the data sets having the underlying correlation values higher than 0.312, which is corresponding to p = 0.01 with Bonferroni correction when degree of freedom (DOF) is the length of reference signal minus two and the number of voxel positions is 57,280. Four metrics were defined to assess the localization accuracy, correlation estimation, and spatial extent of contrast-weighted correlation dynamics obtained in Simulations 2 and 4. The localization error (LE) was calculated as the distance between the ground truth dipole location and the position with the maximum contrast-weighted correlation. The correlation error (CEg) and the orientation error (OE) at the ground truth dipole position were calculated as the differences between the ground truth values and the estimated ones at the ground truth dipole position. To consider the effect of mislocalization, the correlation error at the peak position of contrast-weighted correlation dynamics (CEp) was calculated as the difference between the ground truth correlation value and the estimated one at the position and time corresponding to the peak of contrast-weighted correlation dynamics. In Simulation 5, the sensitivity and spatial resolution were evaluated for SILSC and DICS methods with respect to the regularization parameter

7

by using the MEG data sets generated with different underlying correlation values. Waveforms 1, 2, 3, and 4 were generated using Eq. (18) with parameters {100, 0.3, 17, 5}, {100, 0.6, 17, 5}, {100, 0.8, 17, 5}, and {100, 1, 17, 5}, respectively. For each of Waveforms 1, 2, 3, and 4, one hundred orientations were randomly selected on the tangential plane to generate 100 sets of MEG data for the dipole source at the voxel with MRI coordinate [101, 106, 94]. Waveform 1 was used as the reference signal, which is correlated to Waveforms 2, 3, and 4 at the level of 0.79, 0.56, and 0.34, respectively. The regularization parameter values were set to be 0.0003, 0.003, and 0.03 times the maximum eigenvalue of covariance matrix Cs, represented by 0.0003×, 0.003×, and 0.03×, respectively. SILSC and DICS were individually applied to the 400 data sets for calculating the correlation values between the reference signal and the source activity at every voxel position on the xy-plane at z = 94. Then LE, CEp, CEg, full widths at half maximum (FWHM), and the number of side clusters were calculated from the correlation maps obtained by applying both methods on the simulation data. For DICS, FWHM was calculated as the number of voxel positions having the estimated correlation values larger than the half of the peak value, whereas the contrast-weighted correlation values, instead of the correlation values, were examined in the calculation of FWHM for SILSC. These voxel positions were then clustered on the basis of spatial adjacency. Except for the cluster containing the voxel having peak correlation value, other clusters were regarded as side clusters. Note that the peak location used in the calculation of LE and CEp was determined from the contrastweighted correlation map and correlation map in SILSC and DICS, respectively. In Simulation 6, the correlation dynamics were calculated for the simulation data (ten trials) comprising two spatially distant sources with or without temporal overlap. The temporal profiles of both sources used in this simulation were the same as the source waveform used in Simulation 2, but the late source was delayed by 128 ms or 256 ms and thus had 50% temporal overlap or no temporal overlap with the early one, respectively. The early and late sources were placed at the voxels with MRI coordinate [101, 106, 94] and [121, 106, 94], respectively. To estimate the correlation at different time latency, the center of a 256-ms time window Tc was moved from 100 to 600 ms in 5-ms steps to segment MEG recordings mc(t). The correlation dynamics were calculated for the voxel positions along the x-axis at y = 106 and z = 94. The correlation dynamics obtained by SILSC were then processed by the spurious connectivity reduction procedure, as described in the Reduction of spurious connectivity section, and the voxels with Bonferroni-corrected p b 0.001 were retained with DOF equal to the size of time window Tc minus two. The voxels having contrast values Fθ lower than the contrast threshold were further removed. Finally, the contrast-weighted correlation values were calculated for the remaining voxels by Eq. (15). Facial expression imitation experiment The feasibility of SILSC was further examined by applying it to analyze the data acquired in a facial expression imitation study, in which the aim was to investigate the role of mirror neuron systems in the processing of facial emotions. Grayscale facial images selected from the Taiwanese Facial Expression Image Database (Chen and Yen, 2007) were used as visual stimuli. The experimental paradigms of imitation and execution conditions are shown in Fig. 3. Each facial stimulus was presented for 2 s followed by a 1.5-s resting block in which a fixation cross was presented at the center of the screen. Two male volunteers without neurological or psychiatric disorders participated in this

Fig. 2. Simulation setup and experimental results of Simulation 1. (a) Four dipole sources, s1, s2, s3, and s4, with four temporal waveforms located at three positions were used to generate a set of MEG data; (b) MEG measurements were generated using Eq. (17) and then averaged. The first dipole, s1, was placed at r1 (blue dot) with its temporal waveform shown in blue. The second dipole, s2, was placed at r2 (green dot) with its temporal waveform shown in green, which has the same profile as the temporal waveform of s3 with a 350-ms delay. The dipole sources, s3 and s4, were placed at r3 (yellow dot) with their temporal waveforms shown in dark green and red, respectively. (c) The correlation maps during 200–400 (left) and 550–750 ms (right) with reference signal specified as the source activity during 550–750 ms in the reference region r2, which is the position with the peak F-statistic value calculated using MCB. (d) Four correlation maps calculated using SILSC with reference signal specified as one of the four independent components, u1, u2, u3, and u4, decomposed from the MEG data.

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

Imitation

8

+ neutral

smile

Execution

2s

+ smile

stimulus

+ n eutral

smile

response

1.5 s

+ smile

+

+ smile

+ neutral

+ neutral

stimulus smile

response

Fig. 3. Paradigm of the facial expression imitation experiment. Under the imitation condition, the visual stimuli consisted of facial images presenting either neutral or happy expressions. Under the execution condition, the visual stimuli consisted of facial images of neutral expressions with or without a red cross overlaid on the nose.

Results Simulation 1: reference signal selection In Simulation 1, a reference signal was selected from either independent components or source activity estimated by MCB. Position r2 had the highest F-statistic value and was specified as the reference region. The source signal at r2 during the period of 550–750 ms was selected

as the reference signal with which to calculate correlation dynamics. Correlation maps for the time windows during 200–400 and 550–750 ms reveal that the positions with peak correlation values are at the ground truth positions r3 and r2, respectively, as shown in Fig. 2(c). A total of 114 independent components were decomposed from the preprocessed MEG measurements using FastICA (Hyvärinen and Oja, 2000). Four of the components {u1, u2, u3, u4} were selected as the reference signals through visual inspection. The correlation maps for these reference signals are presented in Fig. 2(d). The peaks of the correlation maps obtained with reference signals specified as u1, u2, u3, and u4 are accurately located at the ground truth positions r1, r3, r3, and r2, respectively. For position r3, the source activity calculated by MCB was similar to the combination of temporal waveforms of s3 and s4, as shown by the pink waveform in Fig. 4. When u2 was selected as the reference signal, the activity calculated by the SILSC beamformer had a similar profile as s4, as shown by the gray waveform in Fig. 4. Simulation 2: influence of SNR After applying SILSC to the 700 sets of simulation data generated in Simulation 2, we used the position and time window having maximum correlation between the source signal and the reference signal to evaluate the accuracy of localization and time latency. The estimated dipole orientation and correlation value at the ground truth position were used to calculate OE and CEg shown in Figs. 5(a) and (b), respectively. The linear fitting lines demonstrate that high SNR resulted in low orientation error and correlation error for the SNR ranging from −10 to −1 dB. The estimated correlation value at the peak position of correlation dynamics was used to calculate CEp shown in Fig. 5(c). The values of CEp are larger than those of CEg, that is, the correlation at

amplitude (10 nAm)

experiment. Under the imitation condition, images of neutral and smiling faces were randomly displayed on the screen with equal chances. Participants were asked to imitate the facial expression of the images. Under execution condition, neutral face images with or without a red cross overlaid on the nose were randomly displayed to the subjects with equal chances. Participants were instructed to smile when presented with a face image with a red cross on the nose and otherwise to remain neutral. This experiment was approved by the Institutional Review Board of Taipei Veterans General Hospital, Taipei, Taiwan, and both subjects provided written informed consents. After applying the preprocessing procedure, 54 and 50 trials of the MEG data remained for Subject 1 and 60 and 57 trials remained for Subject 2 under imitation and execution conditions, respectively. In this study, SILSC was applied to analyze the trials when subjects smiled in response to smiling faces in imitation condition and to neutral faces with a red cross in execution condition. The left inferior frontal gyrus (IFG) was selected as the reference region with which to identify strongly correlated regions of the brain and associated time latencies. The coordinates of this reference region in MNI space were [−64, 12, 14], which were obtained from the group analysis results presented in our previous report (Wang et al., 2008). This position is located at the pars opercularis of left IFG, which is also known as Brodmann area 44. The source activity in this region was calculated using MCB, in which the prestimulus interval during −100 to −200 ms was specified as the control state. The reference signals for the imitation and execution conditions were the segments of source activity in the left IFG during 100–200 ms and 180–280 ms, respectively (Chen, 2009). Note that all of the trials were used to be the reference signals in the calculation of correlation values by SILSC. The gain matrices used for MCB and SILSC were calculated in the same way as those used in simulation studies. The co-registration between the coordinate systems of individual MR volume and the MEG device was performed by aligning three landmarks (nasion, left pre-auricular, and right pre-auricular points) located in both coordinate systems. The regional contrast-weighted correlation of each cerebral region in the AAL atlas was calculated by using the parcellation procedure described in the Parcellation section.

6 4 2 0 -2 -4 -6 -8

SILSC MCB

200

300

400 500 time (ms)

600

Fig. 4. Estimated source signal at position r3 in Simulation 1. The activity calculated using the SILSC beamformer referenced against component u2 is shown in gray; activity calculated using MCB is shown in pink.

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

corr = −0.55, p = 0.00

corr = −0.19, p = 0.00

9

corr = −0.31, p = 0.00

corr = 0.06, p = 0.11

corr = 0.02, p = 0.52

0.05

LE (mm)

CEp

2

6

0.3 0.2 0.1

0 −10

−1

0 −10

SNR (dB)

−1

(a)

−1 SNR (dB)

(b)

4 2

0 −2

0 −10

SNR (dB)

2 TE (ms)

0.1

4 CEg

OE (degree)

0.4

0 −10

−1 SNR (dB)

(c)

(d)

−10

−1 SNR (dB)

(e)

Fig. 5. Performance evaluation for SILSC in Simulation 2. SNR is in the range of −10 to −1dB. Panels (a)–(e) illustrate OE, CEg, CEp, LE, and TE of the data from 700 simulations, respectively. The black line in each panel is the linear fitting line. The correlation coefficient between each evaluation measure and the SNR as well as the corresponding p-value is shown on top of each panel.

peak position (CORRp) is more underestimated than that at the groundtruth position (CORRg), as shown in the first row of Table 1. Figs. 5(d) and (e) demonstrate that LE and TE have low dependence on SNR, respectively. Most of the data sets reveal TE smaller than 2 ms and LE lower than 4 mm.

estimated by SILSC with small regularization parameter values had higher correlation errors, particularly when the data sets comprised the sources with high underlying correlation values. On the contrary, high regularization parameter values resulted in large LE and FWHM, particularly when the data sets had low underlying correlation values.

Simulation 3: surrogate test

Comparison between SILSC and DICS Figs. 8(a)–(b) illustrate the cross-sectional spatial extent of the average correlation values estimated by DICS and SILSC with respect to the regularization parameter. The underlying correlation for this simulation data was set to be one. When the regularization parameter value was lower, the correlation value was more underestimated by SILSC and DICS. The spatial extent of average contrast-weighted correlation values estimated by SILSC, as shown in Fig. 8(c), reveals smaller FWHM than that of correlation values estimated by either DICS or SILSC. Fig. 9 illustrates the average of LE, CEp, CEg, FWHM, and the number of side clusters obtained from the correlation maps calculated by using SILSC and DICS. The performance difference between SILSC and DICS was examined by paired t-tests. In most of the parameter settings, SILSC showed better localization accuracy than DICS. When the regularization parameter value was set to be 0.003× and 0.03×, CEp of SILSC was significantly larger than that of DICS. However, the peak locations of correlation maps estimated by DICS were deviated from the ground truth position under these conditions. When the regularization parameter value increased, the FWHM values obtained by DICS and SILSC became larger. Note that the contrast-weighted correlation map estimated by SILSC was significantly more focal than the correlation map estimated by DICS, particularly for higher regularization parameter values. On the other hand, the contrast-weighted correlation maps estimated by SILSC revealed no side clusters, whereas 55%, 28%, and 24% of data sets with low underlying correlation revealed side clusters on the correlation maps estimated by DICS with regularization parameter value 0.0003×, 0.003×, and 0.03×, respectively.

In Simulation 3, the performance of SILSC was tested by comparing the correlation estimation results of the 700 sets of simulation data and the corresponding surrogate data with zero latency. Table 1 shows the effect of phase shuffling on LE, OE, and absolute correlation value. After phase shuffling, average LE and OE showed 21-mm and 1.6-degree increases, respectively. Comparing the evaluation measures of the original data and those averaged from the associated 10 surrogate data sets, both localization and orientation errors using the surrogate data were significantly higher than those using the original data (p b 0.01, paired t-test). Simulation 4: influence of underlying correlation and SNR A further simulation study was performed to examine the influence of the SNR and underlying correlation between the reference signal and source activity. Figs. 6(a)–(d) illustrate OE, CEg, CEp, and LE with respect to the underlying correlation and Figs. 6(e)–(h) illustrate those with respect to the SNR. The correlation coefficient between each evaluation measure and underlying correlation or SNR as well as its significance level are shown at the top of each panel. Underlying correlation is positively correlated to CEg and CEp, as shown in Figs. 6(b)–(c), whereas increasing SNR reduces OE, as shown in Fig. 6(e). Moreover, the underlying correlation has more influences on CEp than on CEg. LE has low dependence on either underlying correlation or SNR, as shown in Figs. 6(d) and (h). When the underlying correlation was higher than 0.6, 97% of the results showed underestimated correlation values (CEp N 0). When the underlying correlation increases, the correlation value at peak position is more underestimated but the distance of mislocalization does not change significantly. Simulation 5 Influence of underlying correlation and regularization parameter The data sets with higher underlying correlation resulted in lower LE and FWHM, but higher CEp, as shown in Fig. 7. The correlation values

Simulation 6: delayed source The first rows of Figs. 10(a) and (b) illustrate the correlation dynamics estimated along the x-axis for the simulation data comprising two spatially distant sources without and with temporal overlap, respectively. The gray dashed lines indicate the ground truth position and time point for the early source and the black dashed lines indicate those for the late source. The second rows show the correlation dynamics after the reduction of spurious connectivity. The voxels close to the reference

Table 1 Results of performance evaluation in Simulations 2 and 3.

Original data Surrogate data

LE⁎ (mm)

OE⁎ (degree)

CORRg⁎

CORRp⁎

Ground truth CORR⁎

2.65 ± 1.05 23.72 ± 18.81

1.60 ± 1.38 3.21 ± 2.54

0.94 ± 0.02 0.29 ± 0.06

0.88 ± 0.05 0.32 ± 0.05

1.00 0.28 ± 0.06

⁎ The values estimated using the original data and the mean values estimated using the associated surrogate data have significant differences (p b 0.01, paired t-test).

10

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

corr = 0.04, p = 0.05

corr = 0.33, p = 0.00

6

corr = 0.42, p = 0.00

0.1

0.2

0.05

0.1

corr = −0.03, p = 0.10

0

0 −0.05 −0.1

0.4 0.6 0.8 underlying correlation

LE (mm)

2

CEp

CEg

OE (degree)

8 4

0 −0.1 −0.2

0.4 0.6 0.8 underlying correlation

(a) corr = −0.44, p = 0.00

0.2

0.05

0.1

0

0.4 0.6 0.8 underlying correlation

corr = −0.05, p = 0.01 0.1

4 2

(b)

6

6

0.4 0.6 0.8 underlying correlation

(c)

(d)

corr = −0.12, p = 0.00

corr = 0.01, p = 0.74

0 −10

0 −0.05

−8

−6 −4 SNR (dB)

−0.1 −10

−2

LE (mm)

2

CEp

CEg

OE (degree)

8 4

0 −0.1

−8

(e)

−6 −4 SNR (dB)

4 2

−0.2 −10

−2

6

−8

(f)

−6 −4 SNR (dB)

0 −10

−2

(g)

−8

−6 −4 SNR (dB)

−2

(h)

Fig. 6. Influence of underlying correlation and SNR in Simulation 4: Each panel displays the results obtained from the 2864 simulation data sets having underlying correlation above 0.3120, which is corresponding to p = 0.01 adjusted with Bonferroni correction. Panels (a)–(d) illustrate OE, CEg, CEp, and LE with respect to the underlying correlation, respectively, whereas, panels (e)–(h) illustrate those with respect to the SNR, respectively. The black line in each panel is the linear fitting line. The top of each panel shows the correlation coefficient between each evaluation measure and the underlying correlation or SNR as well as the corresponding p-value.

The contrast-weighted correlation maps averaged across 30–70, 80–120, 130–170, 180–220, and 230–270 ms are shown in Fig. 11(b). During 130–170 ms, the strong correlation was revealed in the left SM1, left AG, left IPL, right anterior insula, and right amygdala of both subjects under imitation condition. Under execution condition, Subjects 1 and 2 had strong correlation revealed in the bilateral occipital lobes. During 180–220 ms, the occipital lobes of both subjects also revealed strong correlation under execution condition. During 230–270 ms, the left SM1 and right somatosensory cortex of Subject 2 showed significant correlation under execution condition.

region were removed in both cases. The correlation dynamics were further thresholded to retain the voxels with Bonferroni-corrected p-values smaller than 0.001, as shown in the third rows. The voxels with contrast values lower than the contrast threshold were further removed and the contrast-weighted correlation values for the surviving voxels are shown in the bottom rows. In each case, the position and time with peak contrast-weighted correlation is at around the ground truth.

Facial expression imitation experiment Discussion

0.2

2

0.15

x

x

1.5

x

1 0.5 0 0.34

0.1

This paper proposes a novel algorithm, SILSC, to estimate the correlation dynamics between a given reference signal and brain source activity using electromagnetic signals. The results of simulation studies demonstrate that SILSC is capable of achieving high accuracy in localization and correlation estimation, as well as low latency error. We also applied SILSC to estimate the correlation dynamics of brain activity from MEG data obtained in a facial expression imitation experiment. Individual region of the brain contains numerous neurons and results in multiple source configurations (Gross and Ioannides, 1999). In other words, multiple sources located in the same position produce a variety of temporal profiles and orientations. For example, van der

x x x x

0.05 0

0.56 0.79 1.00 underlying correlation

(a)

−0.05 0.34

0.56 0.79 1.00 underlying correlation

(b)

60 FWHM (voxels)

2.5 x

CEp

LE (mm)

For the left inferior parietal lobule (IPL), left angular gyrus (AG), right insula, right amygdala, left primary sensorimotor (SM1), and left occipital lobe, their temporal dynamics of regional contrast-weighted correlation values are illustrated in Fig. 11(a). For the occipital lobe, activity of six AAL regions were considered, including calcarine fissure, cuneus, ligual gyrus, superior occipital gyrus, middle occipital gyrus, and inferior occipital gyrus. In the imitation condition, the left SM1 area of Subject 2 had peak earlier than that in the execution condition. The left occipital lobes of both subjects showed stronger correlation under execution condition, whereas the left IPL, left AG, right insula, and right amygdala showed stronger correlation under imitation condition.

x x

x

40

α=0.0003x α=0.003x α=0.03x

x 20

0 0.34

0.56 0.79 1.00 underlying correlation

(c)

Fig. 7. Influence of underlying correlation and regularization parameter in Simulation 5. The means of LE, CEp, and FWHM were obtained from the correlation maps calculated by performing SILSC with α set to be 0.03×, 0.003×, and 0.0003× on the 100 MEG data sets generated using various underlying correlation values. Gray X marks indicate the significant differences (p b 0.001, one-way ANOVA) among the results obtained with three regularization parameter values.

DICS

SILSC 1 correlation

1 correlation

contrast−weighted correlation

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

0.5

0 −40 −20 0 20 40 displacement along x−axis (mm)

0.5

0 −40 −20 0 20 40 displacement along x−axis (mm)

(a)

11

SILSC α = 0.03x α = 0.003x α = 0.0003x

15 10 5

0 −40 −20 0 20 40 displacement along x−axis (mm)

(b)

(c)

Fig. 8. Spatial extents of correlation values estimated by using DICS and SILSC with respect to the displacements along the x-axis. The distribution was the average of correlation or contrastweighted correlation values obtained from 100 simulation data sets with the underlying correlation value equal to one. The values of regularization parameter α were set to be 0.0003×, 0.003×, and 0.03×.

Meij et al. (1993) applied dipole fitting to analyze epilepsy data. The estimated dipole sources in different time periods were located in the same position; however, they presented significant differences in orientation. In Simulation 1, we simulated this situation by placing two sources {s3, s4} with different orientations and temporal waveforms in the same position, r3. When the independent component u2 was used as the reference signal, the temporal waveform estimated by the SILSC beamformer located at the peak position of correlation map is similar to the waveform of s4. However, the temporal waveform calculated by MCB located at the same position presents a mixture of the temporal waveforms of s3 and s4, as shown in Fig. 4. The beamformer with the orientation estimated by SILSC is capable of extracting the source component of brain activity with the maximum correlation to the reference signal. In Simulations 2 and 4, the influences of SNR on the performance of correlation and orientation estimation as well as localization were evaluated. The simulation study carried out by Gross et al. (2001) has also demonstrated that high SNR resulted in low correlation error. In the simulation studies of Quraan et al. (2011) and Hillebrand and Barnes

LE (mm)

α=0.0003x

CEp

0 0.15 0.1 0.05 0

10

CEg FWHM (voxels)

0.15 0.1 0.05 0 1000 800 600 400 200

5

* 0.34

0.56

0.79

1.00

0.34

+ 0.34

+

0.56

0.56

0.79

0.79

1.00

1.00

+

+

+

0.56

0.79

1.00

0.5 0

0 0.15 0.1 0.05 0

+ 0.34

number of side clusters

α=0.003x

10 5

(2003), large source amplitude can reduce the localization error. In general, deep sources result in low SNR of MEG data. The simulation study carried out by Zhang et al. (2013) demonstrated that localizing deep sources by beamforming methods had large localization errors compared to the localization of superficial sources. In our study, high SNR of MEG data resulted in high accuracy of orientation and correlation estimation, as shown in Figs. 5(a)–(c) and 6(e)–(g). However, the SNR of MEG data did not significantly affect the localization error, as shown in Figs. 5(d) and 6(h). These results indicate that SILSC has higher noisetolerance on localization accuracy than the conventional beamforming methods. However, the noise affects the performance on orientation and correlation estimation for either SILSC or conventional methods. SILSC is sensitive to detect the sources correlated to the reference signal. In Simulation 3, the peak correlation value estimated from the original MEG data is significantly greater than that obtained from the corresponding surrogate data, as shown in Table 1. Moreover, Simulations 4 and 5 demonstrated that the underlying correlation between source activity and the reference signal affected the accuracy of correlation estimation but not significantly affected the accuracy of localization

0.34

0.56 0.79 1.00 underlying correlation

+ 0.34

+ 0.34

0.15 0.1 0.05 0 1000 800 600 400 200

0.34

+ 0.34

0.5

*

10

0

+ 0.56

0.56

+ 0.79

+ 0.79

0.79

1.00

*

0.15 0.1 0.05 0

1.00

1.00

+

+

+

0.56

0.79

1.00

+ 0.34

0

0.15 0.1 0.05 0 1000 800 600 400 200

0.5

0.56 0.79 1.00 underlying correlation

0

α=0.03x

+

5

+ 0.56

+

SILSC DICS

+

0.34

0.56

0.79

1.00

+

+

+

0.34

0.56

0.79

1.00

0.34

0.56

0.79

1.00

+

+

+

+

0.34

0.56

0.79

1.00

+

+ * 0.34

0.56 0.79 1.00 underlying correlation

Fig. 9. Comparison of DICS and SILSC on the influence of underlying correlation and regularization parameter in Simulation 5. The black circle and black cross indicate the mean value of LE, CEp, CEg, FWHM, and the number of side clusters obtained from 100 sets of simulation data by using SILSC and DICS, respectively. Panels from left to right illustrate the results calculated by using SILSC and DICS with the regularization parameter set to be 0.0003×, 0.003×, and 0.03×. Gray cross marks (p b 0.01, paired t-test) and star marks (p b 0.05, paired t-test) indicate the significant differences between the results obtained by using SILSC and DICS.

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

voxel (x−axis)

300

400

500

600

p<0.001 with Bonferonni correction (DOF=254, |R|>0.3090, n=550) 80 100 120 140 160

200

300

400

500

600

normalized contrast−weighted correlation dynamics (F>2.8768, n=351) 80 100 120 140 160

200

300

400 time (ms)

500

600

voxel (x−axis)

200

300 400 500 600 after spurious connectivity reduction (n=2092)

voxel (x−axis)

200

correlation dynamics (n=10201) 80 100 120 140 160

80 100 120 140 160

voxel (x−axis)

voxel (x−axis)

80 100 120 140 160

voxel (x−axis)

80 100 120 140 160

voxel (x−axis)

correlation dynamics (n=10201)

voxel (x−axis)

12

(a)

200

200

300 400 500 600 after spurious connectivity reduction (n=2289)

300

400

500

600

p<0.001 with Bonferonni correction (DOF=254, |R|>0.3191, n=570) 80 100 120 140 160

200

300

400

500

600

normalized contrast−weighted correlation dynamics (F>2.5909, n=355) 80 100 120 140 160

1 0.5 200

300

400 time (ms)

500

600

0

(b)

Fig. 10. Correlation dynamics estimated by SILSC in Simulation 6 with the two source waveforms having (a) no temporal overlap and (b) 50% temporal overlap. The gray dashed lines indicate the ground-truth position and time corresponding to the early source whereas the black dashed lines indicate those corresponding to the late source. The correlation dynamics after spurious connectivity reduction are illustrated in the second rows. The third rows show the correlation dynamics corresponding to p b 0.001 adjusted with Bonferroni correction. The fourth rows illustrate the dynamics of contrast-weighted correlation values normalized by its maximum (n = number of voxels × number of time windows, DOF = degree of freedom for Student's t-distribution, R = correlation estimated by SILSC, F = contrast).

and orientation estimation. The localization error was negatively correlated to the underlying correlation only when the regularization parameter value was large, as shown in Fig. 7(a). Otherwise, the localization error did not show significant correlation to the underlying correlation, as shown in Figs. 6(d) and 7(a). Figs. 6(b)-(c) and 7(b) demonstrate the positive correlation between the correlation error and the underlying correlation. This relation conforms with the result presented in Gross et al. (2001) that absolute coherence error increased with the underlying coherence. Figs. 2 and 4 demonstrate that the correlation map and the correlated source component estimated by SILSC were affected by the given reference signal. For the cortico-cortical connectivity estimation, the selection of reference region is essential to the accuracy of correlation map estimated by SILSC. The selection of reference region can be knowledge-driven, data-driven, or both. The knowledge-driven selection is better suited for the analysis of data acquired using widelyused experimental paradigms, such as the motor task (Gross et al., 2001), visual task (Siegel et al., 2008), Stroop task (Astolfi et al., 2007), and resting (Cauda et al., 2010). In data-driven selection methods, commonly the region with statistical significance or maximum amplitude is chosen as the reference region (Gross et al., 2001; Friston et al., 2011). In our facial expression imitation study, we used a hybrid procedure for reference region selection: the IFG was first specified according to the physiological knowledge and then the position with peak contrast in the IFG was selected as the reference region (Wang et al., 2008). If it is difficult to determine the reference region, one can consider calculating all-to-all pairwise correlation matrix by sequentially using every brain sources estimated at sample positions as the reference signal (Hillebrand et al., 2012; Schoffelen and Gross, 2011). According to the matrix, the position having significantly high connection density can

then be selected as the reference or nodal point for further analysis (Kujala et al., 2007, 2008). Improving the accuracy of estimated source activity at the reference region would also benefit the performance of SILSC. Compared to minimum norm estimates and dipole fitting methods, beamforming methods can achieve higher sensitivity and spatial resolution (Dalal et al., 2008; Darvas et al., 2004; Sekihara et al., 2005). In this study, MCB was performed to calculate the source activity of the reference region. However, beamformer has a limitation in its suppression of spatially separated but temporally correlated brain sources (Brookes et al., 2007; Hillebrand et al., 2005). Advanced beamforming methods have been proposed to retain the estimated power of correlated sources and thus can achieve higher accuracy than the conventional beamforming methods in the estimation of the reference signal. Huang et al. (2004) replaced the covariance matrices of electromagnetic signals to higher-order ones. In dual-core beamforming methods, spatial filters were calculated by considering a linear combination of lead fields for two sources (Brookes et al., 2007; Diwakar et al., 2011a). Nulling beamformer added the null constraints in the derivation of spatial filter to suppress the influence of correlated sources (Dalal et al., 2006; Hui et al., 2010; Quraan and Cheyne, 2010). Diwakar et al. (2011b) further proposed an enhanced dual-core beamformer, which can calculate the spatial filter by utilizing both the dual source model and the null constraints. The selection of regularization parameter value is a tradeoff between the spatial resolution of beamforming results and the sensitivity to noise (Brookes et al., 2008; Woolrich et al., 2011). A larger regularization parameter value leads to more restriction on the norm of spatial filter and less constraint on the variance of beamformer output, resulting in more leakage from noise and the activity at non-targeted locations. Because

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

left IPL

regional contrast-weighted correlation

8 6 4 2 0

50

100

150

200

13

left AG

250

300

350

8 6 4 2 0

s1-imi s1-exe s2-imi s2-exe

50

100

right insula

150

200

250

300

350

300

350

300

350

right amygdala

5

15 10 5

0

50

100

150

200

250

300

350

0

50

100

left SM1 8 6 4 2 0

50

100

150

200

150

200

250

left occipital lobe

250

300

350

50 40 30 20 10 0

50

time (ms)

100

150

200

250

time (ms)

(a) 130-170 ms

180-220 ms

230-270 ms

imitation

execution

imitation

80-120 ms

execution

Subject 2

Subject 1

30-70 ms

(b)

0

.15

Fig. 11. Contrast-weighted correlation dynamics obtained in the facial expression imitation experiment. Panel (a) displays the temporal dynamics of the regional contrast-weighted correlation in the six regions: the left IPL, left AG, right insula, right amygdala, left SM1, and left occipital lobe. Panel (b) displays the spatial maps of the contrast-weighted correlation values averaged across 30–70, 80–120, 130–170, 180–220, and 230–270 ms.

the beamformer output contains both the activity of the targeted location and the nearby sources, the more leakage means that the effect of spatial averaging is stronger and in a wider range. The results of Simulation 5 demonstrated that a larger regularization parameter value caused larger spatial extent and localization error but smaller correlation error, as shown in Figs. 7(a)–(c). Based on the results of Simulation 5, the regularization parameter value used in Simulations 1–4 and 6 was set to be

0.003 × to avoid large correlation and localization errors. Chen et al. (2006) suggested that the scalar beamformer with small regularization can have high accuracy for the data with high SNR. Because real data usually contains high level of noise and a large number of background sources, a larger regularization parameter value, 0.03×, was used for analyzing the facial expression imitation data. One way of objective determination of this value is to use Bayesian principle component analysis to

14

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

automatically choose the amount of regularization required (Woolrich et al., 2011). The first and fourth rows of Fig. 9 show that SILSC outperformed DICS in terms of localization accuracy and spatial extent. As shown in the fifth row, moreover, DICS might obtain side clusters, causing the ambiguity of localization. For all of the 400 data generated in Simulation 5, SILSC obtained no side clusters. The peak position and ground truth position were within the single main cluster. Therefore, the locations of the peaks in the correlation maps estimated by using SILSC are more accurate and reliable compared to those estimated by using DICS. Beamforming is a scanning method performed for each grid point within a region of interest. In Simulations 2, 3, and 4, SILSC was applied on the grid points regularly sampled every four voxels along each dimension of MRI data. Because the dipole sources were randomly deployed without restriction to these grid points, average LE depended on the spatial sampling rate. In this study, the largest distance between the dipole position and its closest sampling grid point was about 4 mm. In Simulation 2, the LE was generally lower than 4 mm, as shown in Fig. 5(d), demonstrating the high localization accuracy of SILSC. In our facial expression imitation study, the left IPL, left AG, left SM1, right insula, and right amygdala revealed the higher regional contrastweighted correlation during 100–200 ms under imitation condition than that under execution condition, as shown in Fig. 11(a). The pars opercularis of the IFG, specified as the reference region in this study, and the IPL are the core of mirror neuron system (Rizzolatti and Craighero, 2004; Molenberghs et al., 2012). The AG has also been reported to have mirror properties (Ramachandran and Oberman, 2006), responsible to combine the visual and somatosensory information. The SM1 is a critical component of extended mirror neuron system (Pineda, 2008). The frontoparietal circuit consisting of IPL, AG, and SM1 regions has been implicated in action understanding and action execution (Binkofski and Buccino, 2006). The insula and amygdala are involved in the extended mirror neuron system (van der Gaag et al., 2007) and are associated with the recognition of facial expression and emotion. The connectivity between the IFG and anterior insula has been observed in the fMRI study of emotional facial expressions carried out by Jabbi and Keysers (2008). Lee et al. (2006) also demonstrated the correlation between the magnitude of facial movement and the insula activity during emotion imitation. These results demonstrated the correlation dynamics of IFG-referenced mirror neuron system during the process of facial expression imitation, suggesting that the proposed method is capable of investigating the dynamic changes of corticocortical functional connectivity in emotional processing. For both subjects in our study, the left occipital lobe revealed higher contrast-weighted correlation during 100–220 ms under execution condition than that under imitation condition. The high correlation under execution condition was also shown in the left SM1 after 200 ms, which was later than the peak observed under imitation condition. For the execution condition, we selected the reference signal from the activity of left IFG during 180–280 ms. Therefore, the contrastweighted correlation dynamics revealed the correlated cortical activity progressing from the occipital lobe to the left IFG, and finally to the left SM1. Apart from the activation of occipital lobe for primary visual processing, our results conform with those of Nishitani and Hari (2002) that the left IFG and left M1 were activated sequentially under the spontaneous execution condition. The left IFG is activated for motor planning (Johnson and Grafton, 2003), whereas the SM1 is responsible to the execution of action. Therefore, the regions highly correlated to the left IFG may comprise the motor control circuit under execution condition. In addition, Bortoletto et al. (2011) reported that the motor representations of intended actions can facilitate the encoding of the visual inputs of observed actions. The direct influence of motor planning on the visual processing may explain the connection between the left IFG and occipital lobe observed in this study. Under imitation and execution conditions, the regions correlated to the pars opercularis of the left IFG are with mirror properties and related

to motor control, respectively. This difference may result from the functional segregation of pars opercularis. The dorsal sector of pars opercularis is with mirror properties and responsible to understand the observed action whereas the ventral sector is a premotor area for the movement of hand and face (Carr et al., 2003; Molnar-Szakacs et al., 2005). Therefore, SILSC can be used to investigate the functional segregation of a region by specifying it as the reference region for different experimental conditions. Two kinds of reference signals were used in this work: the source activity at the reference region and the independent component. In our facial expression imitation experiment, the left IFG was selected as the reference region for the calculation of correlation dynamics because it was reported as a crucial region in mirror neuron system and motor control. When there are no regions of interest, ICA is an alternative way to select reference signals. In this case, SILSC can provide the correlation map corresponding to the specified independent component. At the same position, SILSC may obtain different source activity when different independent components are specified to be reference. SILSC can obtain the source component which is the most correlated to the selected independent component and reduce the influence of the other source components, as shown in Fig. 4. In Simulation 1, the reference signals were selected from independent components by visual inspection. When analyzing real data, visual inspection tends to be inconsistent and subjective (Naeem et al., 2006). Numerous component selection methods can be employed to select reference signals, such as: template-based approach (Wessel and Ullsperger, 2011; Yang et al., 2011), clustering analysis method (Milanesi et al., 2009), and atlas-based statistical analysis method (Chan et al., 2012). The source localization algorithm presented by Hild and Nagarajan (2009), called column correlation, uses correlation values between the lead field and the spatial map of independent component to express the confidence level of a neural source existing at the targeted location. When using an independent component as the reference signal, the correlation value estimated by SILSC may represent the confidence level that a source related to this independent component activates at the targeted location. The source distribution of an independent component can also be calculated by conventional MEG/EEG source imaging algorithms, such as dipole fitting (Makeig et al., 2004) and minimum norm estimation (Yang et al., 2011; Chan et al., 2012), or by electromagnetic spatiotemporal independent component analysis (EMSICA) (Tsai et al., 2006), which simultaneously estimates time courses and source distributions of independent components. Compared to these methods, SILSC can further calculate the correlation dynamics corresponding to each independent component, which detects the correlated sources with latencies. In source coherence/correlation imaging, the artificial interactions due to the source leakage are inevitable for all inverse algorithms (Sekihara and Nagarajan, 2008; Sekihara et al., 2011). As conventional beamformers, the source activity estimated by SILSC also contains leakage from other sources at non-targeted locations. To reduce the artificial interactions caused by leakages, this work detected the spurious connectivity by using the correlation estimated from the simulation data comprising two uncorrelated sources, which has been applied to the studies of Brookes et al. (2011a,b). Alternatively, the component of reference signal could be removed from the estimated source activity at the targeted region using a linear projection before the calculation of correlation and coherence (Brookes et al., 2012). The estimation of beamforming spatial filter could use the power of non-zero phaselagged interactions, obtained by eigenvector decomposition of the imagery part of cross-spectral density, to reduce the spurious connectivity caused by volume conduction or magnetic field spread (Drakesmith et al., 2013). Note that these two methods face the risk of losing true zero-lag interaction. SILSC uses a single-step procedure to quantify the linear interdependencies between brain sources and the specified reference signal. However, interactions between brain regions can be linear or nonlinear

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

(David et al., 2004). To estimate nonlinear interactions, numerous measures have been proposed, such as generalized synchronization (Stam and van Dijk, 2002), phase locking value (Lachaux et al., 1999), phase lag index (Stam et al., 2007), and mutual information (Jeong et al., 2001; Liu and Ioannides, 2006). These nonlinear measures for the estimation of source connectivity have to be evaluated in a twostep procedure: source activity of each brain region is first estimated by an inverse algorithm, and then the measures are employed to quantify the interactions between regions. Comparing to the nonlinear measures, the linear measures, including correlation and coherence, might be more sensitive to the changes in amplitude, phase, relative timing, and SNR (Schoffelen and Gross, 2009; Sakkalis, 2011). Moreover, both SILSC and DICS methods can estimate source activity and thus can be the inverse algorithm in the two-step procedure of nonlinear interaction estimation methods. Compared to the conventional inverse algorithms, SILSC or DICS is more preferable because the influence of uncorrelated source components is reduced at the first place. In the experiment carried out by Kujala et al. (2007), the reading-related network was detected by the phase synchrony between the source components estimated by DICS. It is also possible to calculate the nonlinear measures on the correlated source components estimated by SILSC. In conclusion, this paper presents a source imaging algorithm, SILSC, which is capable of estimating the spatiotemporal dynamics of correlated brain sources from electromagnetic signals. According to the results of our simulation studies, the localization accuracy of SILSC is robust to noise and underlying correlation of electromagnetic data. SILSC is sensitive to detect sources that are correlated to the reference signal. In addition, SILSC can facilitate the visualization of the spatiotemporal distribution of these sources. When source activity is used as the reference signal, SILSC can provide the temporal dynamics of cortio-cortical interactions, which is essential to construct the human functional connectome. The feasibility of SILSC was demonstrated by using MEG data acquired in an experiment of facial expression imitation. Acknowledgments This work was supported in part by the Taiwan Ministry of Science and Technology Graduate Student Study Abroad Program under grant NSC-102-2917-I-009-038, the UST-UCSD International Center of Excellence in Advanced Bioengineering sponsored by the Taiwan Ministry of Science and Technology I-RiCE Program under grant MOST-103-2911-I-009-101, the Taipei Veterans General Hospital under grant V102E3-004, and by the Taiwan Ministry of Science and Technology under grants NSC-102-2410-H-010-003-MY2 and MOST103-2221-E-009-131. References Ansari-Asl, K., Bellanger, J.-J., Bartolomei, F., Wendling, F., Senhadji, L., 2005. Timefrequency characterization of interdependencies in nonstationary signals: application to epileptic EEG. IEEE Trans. Biomed. Eng. 52 (7), 1218–1226. Ansari-Asl, K., Senhadji, L., Bellanger, J.-J., Wendling, F., 2006. Quantitative evaluation of linear and nonlinear methods characterizing interdependencies between brain signals. Phys. Rev. E 74 (3). Arnhold, J., Grassberger, P., Lehnertz, K., Elger, C.E., 1999. A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D 134 (4), 419–430. Ashburner, J., 2000. Computational Neuroanatomy. Ph.D. thesis. University College London, London, UK. Astolfi, L., Cincotti, F., Mattia, D., Marciani, M.G., Baccala, L.A., de Vico Fallani, F., Salinari, S., Ursino, M., Zavaglia, M., Ding, L., Edgar, J.C., Miller, G.A., He, B., Babiloni, F., 2007. Comparison of different cortical connectivity estimators for high-resolution EEG recordings. Hum. Brain Mapp. 28 (2), 143–157. Averbeck, B.B., Lee, D., 2004. Coding and transmission of information by neural ensembles. Trends Neurosci. 27 (4), 225–230. Azouz, R., Gray, C.M., 2000. Dynamic spike threshold reveals a mechanism for synaptic coincidence detection in cortical neurons in vivo. Proc. Natl. Acad. Sci. U. S. A. 97 (14), 8110–8115. Bell, A.J., Sejnowski, T.J., 1995. An information maximization approach to blind separation and blind deconvolution. Neural Comput. 7 (6), 1129–1159.

15

Binkofski, F., Buccino, G., 2006. The role of ventral premotor cortex in action execution and action understanding. J. Physiol. Paris 99 (4–6), 396–405. Biswal, B.B., Mennes, M., Zuo, X.-N., Gohel, S., Kelly, C., Smith, S.M., Beckmann, C.F., Adelstein, J.S., Buckner, R.L., Colcombe, S., Dogonowski, A.-M., Ernst, M., Fair, D., Hampson, M., Hoptman, M.J., Hyde, J.S., Kiviniemi, V.J., Kötter, R., Li, S.-J., Lin, C.-P., Lowe, M.J., Mackay, C., Madden, D.J., Madsen, K.H., Margulies, D.S., Mayberg, H.S., McMahon, K., Monk, C.S., Mostofsky, S.H., Nagel, B.J., Pekar, J.J., Peltier, S.J., Petersen, S.E., Riedl, V., Rombouts, S.A.R.B., Rypma, B., Schlaggar, B.L., Schmidt, S., Seidler, R.D., Siegle, G.J., Sorg, C., Teng, G.-J., Veijola, J., Villringer, A., Walter, M., Wang, L., Weng, X.-C., Whitfield-Gabrieli, S., Williamson, P., Windischberger, C., Zang, Y.-F., Zhang, H.-Y., Castellanos, F.X., Milham, M.P., 2010. Toward discovery science of human brain function. Proc. Natl. Acad. Sci. U. S. A. 107 (10), 4734–4739. Bortoletto, M., Mattingley, J.B., Cunnington, R., 2011. Action intentions modulate visual processing during action perception. Neuropsychologia 49 (7), 2097–2104. Brazier, M., Casby, J., 1952. Crosscorrelation and autocorrelation studies of electroencephalographic potentials. Electroencephalogr. Clin. Neurophysiol. 4 (2), 201–211. Brookes, M.J., Stevenson, C.M., Barnes, G.R., Hillebrand, A., Simpson, M.I.G., Francis, S.T., Morris, P.G., 2007. Beamformer reconstruction of correlated sources using a modified source model. Neuroimage 34 (4), 1454–1465. Brookes, M.J., Vrba, J., Robinson, S.E., Stevenson, C.M., Peters, A.M., Barnes, G.R., Hillebrand, A., Morris, P.G., 2008. Optimising experimental design for MEG beamformer imaging. Neuroimage 39 (4), 1788–1802. Brookes, M.J., Hale, J.R., Zumer, J.M., Stevenson, C.M., Francis, S.T., Barnes, G.R., Owen, J.P., Morris, P.G., Nagarajan, S.S., 2011a. Measuring functional connectivity using MEG: methodology and comparison with fcMRI. Neuroimage 56 (3), 1082–1104. Brookes, M.J., Woolrich, M., Luckhoo, H., Price, D., Hale, J.R., Stephenson, M.C., Barnes, G.R., Smith, S.M., Morris, P.G., 2011b. Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc. Natl. Acad. Sci. U. S. A. 108 (40), 16783–16788. Brookes, M.J., Woolrich, M.W., Barnes, G.R., 2012. Measuring functional connectivity in MEG: a multivariate approach insensitive to linear source leakage. Neuroimage 63 (2), 910–920. Carr, L., Iacoboni, M., Dubeau, M.-C., Mazziotta, J.C., Lenzi, G.L., 2003. Neural mechanisms of empathy in humans: a relay from neural systems for imitation to limbic areas. Proc. Natl. Acad. Sci. U. S. A. 100 (9), 5497–5502. Cauda, F., Geminiani, G., D'Agata, F., Sacco, K., Duca, S., Bagshaw, A.P., Cavanna, A.E., 2010. Functional connectivity of the posteromedial cortex. PLoS One 5 (9), e13107. Chan, H.-L., Chen, Y.-S., Chen, L.-F., 2012. Selection of independent components based on cortical mapping of electromagnetic activity. J. Neural Eng. 9 (5), 056006. Chen, I.-T., 2009. Beamformer-based Spatiotemporal Imaging of Correlated Brain Activities. Master's thesis. National Chiao Tung University, Hsinchu, Taiwan. Chen, L.-F., Yen, Y.-S., 2007. Taiwanese Facial Expression Image Database. Brain Mapping Laboratory, Institute of Brain Science, National Yang-Ming University, Taipei, Taiwan. Chen, Y.-S., Cheng, C.-Y., Hsieh, J.-C., Chen, L.-F., 2006. Maximum contrast beamformer for electromagnetic mapping of brain activity. IEEE Trans. Biomed. Eng. 53 (9), 1765–1774. Chong, E.K.P., Zak, S.H., 2004. An Introduction to Optimization. 2nd edition. Wiley. Comon, P., 1994. Independent component analysis, a new concept? Signal Process. 36 (3), 287–314. Dalal, S.S., Sekihara, K., Nagarajan, S.S., 2006. Modified beamformers for coherent source region suppression. IEEE Trans. Biomed. Eng. 53 (7), 1357–1363. Dalal, S.S., Guggisberg, A.G., Edwards, E., Sekihara, K., Findlay, A.M., Canolty, R.T., Berger, M.S., Knight, R.T., Barbaro, N.M., Kirsch, H.E., Nagarajan, S.S., 2008. Five-dimensional neuroimaging: localization of the time frequency dynamics of cortical activity. Neuroimage 40 (4), 1686–1700. Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., 2004. Mapping human brain function with MEG and EEG: methods and validation. NeuroImage 23 (Supplement 1), S289–S299 (0). David, O., Cosmelli, D., Friston, K.J., 2004. Evaluation of different measures of functional connectivity using a neural mass model. Neuroimage 21, 659–673. Diwakar, M., Huang, M.-X., Srinivasan, R., Harrington, D.L., Robb, A., Angeles, A., Muzzatti, L., Pakdaman, R., Song, T., Theilmann, R.J., Lee, R.R., 2011a. Dual-core beamformer for obtaining highly correlated neuronal networks in MEG. Neuroimage 54 (1), 253–263. Diwakar, M., Tal, O., Liu, T.T., Harrington, D.L., Srinivasan, R., Muzzatti, L., Song, T., Theilmann, R.J., Lee, R.R., Huang, M.-X., 2011b. Accurate reconstruction of temporal correlation for neuronal sources using the enhanced dual-core MEG beamformer. Neuroimage 56 (4), 1918–1928. Drakesmith, M., El-Deredy, W., Welbourne, S., 2013. Reconstructing coherent networks from electroencephalography and magnetoencephalography with reduced contamination from volume conduction or magnetic field spread. PLoS One 8 (12), e81553. Engel, A.K., Singer, W., 2001. Temporal binding and the neural correlates of sensory awareness. Trends Cogn. Sci. 5 (1), 16–25. Engel, A.K., Fries, P., König, P., Brecht, M., Singer, W., 1999. Temporal binding, binocular rivalry, and consciousness. Conscious. Cogn. 8 (2), 128–151. Friston, K.J., Li, B., Daunizeau, J., Stephan, K.E., 2011. Network discovery with DCM. Neuroimage 56 (3), 1202–1221. Gevins, A., Bressler, S., 1988. Functional Topography of the Human Brain. Hans Huber Publishers, Bern, pp. 99–116. Gevins, A., Cutillo, B., 1993. Spatiotemporal dynamics of component processes in human working memory. Electroencephalogr. Clin. Neurophysiol. 87 (3), 128–143. Gevins, A.S., Morgan, N.H., Bressler, S.L., Cutillo, B.A., White, R.M., Illes, J., Greer, D.S., Doyle, J.C., Zeitlin, G.M., 1987. Human neuroelectric patterns predict performance accuracy. Science 235 (4788), 580–585. Gevins, A.S., Bressler, S.L., Morgan, N.H., Cutillo, B.A., White, R.M., Greer, D.S., Illes, J., 1989a. Event-related covariances during a bimanual visuomotor task I: methods and analysis of stimulus-locked and response-locked data. Electroencephalogr. Clin. Neurophysiol. 74 (1), 58–75.

16

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17

Gevins, A.S., Cutillo, B.A., Bressler, S.L., Morgan, N.H., White, R.M., Illes, J., Greer, D.S., 1989b. Event-related covariances during a bimanual visuomotor task II: preparation and feedback. Electroencephalogr. Clin. Neurophysiol. 74 (2), 147–160. Gevins, A., Leong, H., Smith, M.E., Le, J., Du, R., 1995. Mapping cognitive brain function with modern high-resolution electroencephalography. Trends Neurosci. 18 (10), 429–436. Goldenholz, D.M., Ahlfors, S.P., Hämäläinen, M.S., Sharon, D., Ishitobi, M., Vaina, L.M., Stufflebeam, S.M., 2009. Mapping the signal-to-noise-ratios of cortical sources in magnetoencephalography and electroencephalography. Hum. Brain Mapp. 30 (4), 1077–1086. Grasman, R.P.P.P., Huizenga, H.M., Waldorp, L.J., Bocker, K.B.E., Molenaar, P.C.M., 2004. Frequency domain simultaneous source and source coherence estimation with an application to MEG. IEEE Trans. Biomed. Eng. 51 (1), 45–55. Gross, J., Ioannides, A.A., 1999. Linear transformations of data space in MEG. Phys. Med. Biol. 44 (8), 2081–2097. Gross, J., Kujala, J., Hämäläinen, M., Timmermann, L., Schnitzler, A., Salmelin, R., 2001. Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc. Natl. Acad. Sci. U. S. A. 98 (2), 694–699. Gross, J., Timmermann, L., Kujala, J., Dirks, M., Schmitz, F., Salmelin, R., Schnitzler, A., 2002. The neural basis of intermittent motor control in humans. Proc. Natl. Acad. Sci. U. S. A. 99 (4), 2299–2302. Guevara, M.A., Corsi-Cabrera, M., 1996. EEG coherence or EEG correlation? Int. J. Psychophysiol. 23 (3), 145–153. Hämäläinen, M.S., Ilmoniemi, R.J., 1994. Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32 (1), 35–42. Hild, K.E., Nagarajan, S.S., 2009. Source localization of EEG/MEG data by correlating columns of ICA and lead field matrices. IEEE Trans. Biomed. Eng. 56 (11), 2619–2626. Hillebrand, A., Barnes, G.R., 2003. The use of anatomical constraints with MEG beamformers. Neuroimage 20 (4), 2302–2313. Hillebrand, A., Singh, K.D., Holliday, I.E., Furlong, P.L., Barnes, G.R., 2005. A new approach to neuroimaging with magnetoencephalography. Hum. Brain Mapp. 25 (2), 199–211. Hillebrand, A., Barnes, G.R., Bosboom, J.L., Berendse, H.W., Stam, C.J., 2012. Frequencydependent functional connectivity within resting-state networks: an atlas-based MEG beamformer solution. Neuroimage 59 (4), 3909–3921. Hindriks, R., Bijma, F., van Dijk, B.W., Stam, C.J., van der Werf, Y.D., van Someren, E.J.W., de Munck, J.C., van der Vaart, A.W., 2011. Data-driven modeling of phase interactions between spontaneous MEG oscillations. Hum. Brain Mapp. 32 (7), 1161–1178. Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M., Engel, A.K., 2012. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15 (6), 884–890. Hoechstetter, K., Bornfleth, H., Weckesser, D., Ille, N., Berg, P., Scherg, M., 2004. BESA source coherence: a new method to study cortical oscillatory coupling. Brain Topogr. 16 (4), 233–238. Holmes, C.J., Hoge, R., Collins, L., Woods, R., Toga, A.W., Evans, A.C., 1998. Enhancement of MR images using registration for signal averaging. J. Comput. Assist. Tomogr. 22 (2), 324–333. Huang, M.X., Shih, J.J., Lee, R.R., Harrington, D.L., Thoma, R.J., Weisend, M.P., Hanlon, F., Paulson, K.M., Li, T., Martin, K., Miller, G.A., Canive, J.M., 2004. Commonalities and differences among vectorized beamformers in electromagnetic source imaging. Brain Topogr. 16 (3), 139–158. Hui, H.B., Pantazis, D., Bressler, S.L., Leahy, R.M., 2010. Identifying true cortical interactions in MEG using the nulling beamformer. Neuroimage 49 (4), 3161–3174. Hyvärinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10 (3), 626–634. Hyvärinen, A., Oja, E., 2000. Independent component analysis: algorithms and applications. Neural Netw. 13 (4), 411–430. Jabbi, M., Keysers, C., 2008. Inferior frontal gyrus activity triggers anterior insula response to emotional facial expressions. Emotion 8 (6), 775–780. Jeong, J., Gore, J.C., Peterson, B.S., 2001. Mutual information analysis of the EEG in patients with Alzheimer's disease. Clin. Neurophysiol. 112 (5), 827–835. Johnson, S.H., Grafton, S.T., 2003. From ‘acting on’ to ‘acting with’: the functional anatomy of object-oriented action schemata. Prog. Brain Res. 142, 127–139. König, P., Engel, A.K., 1995. Correlated firing in sensory-motor systems. Curr. Opin. Neurobiol. 5 (4), 511–519. Kujala, J., Pammer, K., Cornelissen, P., Roebroeck, A., Formisano, E., Salmelin, R., 2007. Phase coupling in a cerebro-cerebellar network at 8–13 Hz during reading. Cereb. Cortex 17 (6), 1476–1485. Kujala, J., Gross, J., Salmelin, R., 2008. Localization of correlated network activity at the cortical level with MEG. Neuroimage 39 (4), 1706–1720. Lachaux, J.P., Rodriguez, E., Martinerie, J., Varela, F.J., 1999. Measuring phase synchrony in brain signals. Hum. Brain Mapp. 8 (4), 194–208. Lee, T.-W., Josephs, O., Dolan, R.J., Critchley, H.D., 2006. Imitating expressions: emotionspecific neural substrates in facial mimicry. Soc. Cogn. Affect. Neurosci. 1 (2), 122–135. Liu, L., Ioannides, A.A., 2006. Spatiotemporal dynamics and connectivity pattern differences between centrally and peripherally presented faces. Neuroimage 31 (4), 1726–1740. Makeig, S., Debener, S., Onton, J., Delorme, A., 2004. Mining event-related brain dynamics. Trends Cogn. Sci. 8 (5), 204–210. Marrelec, G., Daunizeau, J., Pelegrini-Issac, M., Doyon, J., Benali, H., 2005. Conditional correlation as a measure of mediated interactivity in fMRI and MEG/EEG. IEEE Trans. Signal Process. 53 (9), 3503–3516. Menon, V., Freeman, W.J., Cutillo, B.A., Desmond, J.E., Ward, M.F., Bressler, S.L., Laxer, K.D., Barbaro, N., Gevins, A.S., 1996. Spatio-temporal correlations in human gamma band electrocorticograms. Electroencephalogr. Clin. Neurophysiol. 98 (2), 89–102. Milanesi, M., James, C.J., Martini, N., Menicucci, D., Gemignani, A., Ghelarducci, B., Landini, L., 2009. Objective selection of EEG late potentials through residual dependence estimation of independent components. Physiol. Meas. 30 (8), 779–794.

Molenberghs, P., Cunnington, R., Mattingley, J.B., 2012. Brain regions with mirror properties: a meta-analysis of 125 human fMRI studies. Neurosci. Biobehav. Rev. 36 (1), 341–349. Molnar-Szakacs, I., Iacoboni, M., Koski, L., Mazziotta, J.C., 2005. Functional segregation within pars opercularis of the inferior frontal gyrus: evidence from fMRI studies of imitation and action observation. Cereb. Cortex 15 (7), 986–994. Mormann, F., Lehnertz, K., David, P., Elger, C.E., 2000. Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients. Physica D 144 (3–4), 358–369. Mosher, J.C., Leahy, R.M., Lewis, P.S., 1999. EEG and MEG: forward solutions for inverse methods. IEEE Trans. Biomed. Eng. 46 (3), 245–259. Naeem, M., Brunner, C., Leeb, R., Graimann, B., Pfurtscheller, G., 2006. Seperability of fourclass motor imagery data using independent components analysis. J. Neural Eng. 3 (3), 208–216. Nishitani, N., Hari, R., 2002. Viewing lip forms: cortical dynamics. Neuron 36 (6), 1211–1220. Nunez, P., Srinivasan, R., 2006. Electric Fields of the Brain: The Neurophysics of EEG. Oxford University Press. Nunez, P.L., Srinivasan, R., Westdorp, A.F., Wijesinghe, R.S., Tucker, D.M., Silberstein, R.B., Cadusch, P.J., 1997. EEG coherency I: statistics, reference electrode, volume conduction, Laplacians, cortical imaging, and interpretation at multiple scales. Electroencephalogr. Clin. Neurophysiol. 103 (5), 499–515. Palva, S., Palva, J.M., 2012. Discovering oscillatory interaction networks with M/EEG: challenges and breakthroughs. Trends Cogn. Sci. 16 (4), 219–230. Pineda, J.A., 2008. Sensorimotor cortex as a critical component of an ‘extended’ mirror neuron system: does it solve the development, correspondence, and control problems in mirroring? Behav. Brain Funct. 4, 47. Quiroga, R.Q., Kraskov, A., Kreuz, T., Grassberger, P., 2002. Performance of different synchronization measures in real data: a case study on electroencephalographic signals. Phys. Rev. E 65 (4). Quraan, M.A., Cheyne, D., 2010. Reconstruction of correlated brain activity with adaptive spatial filters in MEG. Neuroimage 49 (3), 2387–2400. Quraan, M.A., Moses, S.N., Hung, Y., Mills, T., Taylor, M.J., 2011. Detection and localization of hippocampal activity using beamformers with meg: a detailed investigation using simulations and empirical data. Hum. Brain Mapp. 32 (5), 812–827. Ramachandran, V.S., Oberman, L.M., 2006. Broken mirrors: a theory of autism. Sci. Am. 295 (5), 62–69. Rizzolatti, G., Craighero, L., 2004. The mirror-neuron system. Annu. Rev. Neurosci. 27, 169–192. Robinson, S., Vrba, J., 1998. Functional Neuroimaging by Synthetic Aperture Magnetometry. Tohoku Univ. Press, Sendai, pp. 302–305. Rulkov, N.F., Sushchik, M.M., Tsimring, L.S., Abarbanel, H.D.I., 1995. Generalized synchronization of chaos in directionally coupled chaotic systems. Phys. Rev. E 51 (2), 980–994. Sakkalis, V., 2011. Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG. Comput. Biol. Med. 41 (12), 1110–1117. Salinas, E., Sejnowski, T.J., 2001. Correlated neuronal activity and the flow of neural information. Nat. Rev. Neurosci. 2 (8), 539–550. Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol. 32 (1), 11. Schoffelen, J.-M., Gross, J., 2009. Source connectivity analysis with MEG and EEG. Hum. Brain Mapp. 30 (6), 1857–1865. Schoffelen, J.-M., Gross, J., 2011. Improving the interpretability of all-to-all pairwise source connectivity analysis in MEG with nonhomogeneous smoothing. Hum. Brain Mapp. 32 (3), 426–437. Sekihara, K., Nagarajan, S.S., 2008. Adaptive Spatial Filters for Electromagnetic Brain Imaging. Springer. Sekihara, K., Sahani, M., Nagarajan, S.S., 2005. Localization bias and spatial resolution of adaptive and non-adaptive spatial filters for MEG source reconstruction. Neuroimage 25 (4), 1056–1067. Sekihara, K., Owen, J.P., Trisno, S., Nagarajan, S.S., 2011. Removal of spurious coherence in MEG source-space coherence analysis. IEEE Trans. Biomed. Eng. 58 (11), 3121–3129. Siegel, M., Donner, T.H., Oostenveld, R., Fries, P., Engel, A.K., 2008. Neuronal synchronization along the dorsal visual pathway reflects the focus of spatial attention. Neuron 60 (4), 709–719. Silfverhuth, M.J., Hintsala, H., Kortelainen, J., Seppanen, T., 2012. Experimental comparison of connectivity measures with simulated EEG signals. Med. Biol. Eng. Comput. 50 (7), 683–688. Singer, W., Gray, C.M., 1995. Visual feature integration and the temporal correlation hypothesis. Annu. Rev. Neurosci. 18, 555–586. Smith, S.M., 2012. The future of FMRI connectivity. Neuroimage 62 (2), 1257–1266. Stam, C.J., van Dijk, B.W., 2002. Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets. Physica D 163 (3–4), 236–251. Stam, C.J., Nolte, G., Daffertshofer, A., 2007. Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Hum. Brain Mapp. 28 (11), 1178–1193. Tass, P., Rosenblum, M.G., Weule, J., Kurths, J., Pikovsky, A., Volkmann, J., Schnitzler, A., Freund, H.J., 1998. Detection of n:m phase locking from noisy data: application to magnetoencephalography. Phys. Rev. Lett. 81 (15), 3291–3294. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., Farmer, J.D., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physica D 58 (1–4), 77–94. Tsai, A.C., Liou, M., Jung, T.-P., Onton, J.A., Cheng, P.E., Huang, C.-C., Duann, J.-R., Makeig, S., 2006. Mapping single-trial EEG records on the cortical surface through a spatiotemporal modality. Neuroimage 32 (1), 195–207. van der Gaag, C., Minderaa, R.B., Keysers, C., 2007. Facial expressions: what the mirror neuron system can and cannot tell us. Soc. Neurosci. 2 (3–4), 179–222. van der Meij, W., Wieneke, G.H., van Huffelen, A.C., 1993. Dipole source analysis of rolandic spikes in benign rolandic epilepsy and other clinical syndromes. Brain Topogr. 5 (3), 203–213.

H.-L. Chan et al. / NeuroImage 114 (2015) 1–17 Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867–880. von der Malsburg, C., 1999. The what and why of binding: the modeler's perspective. Neuron 24 (1), 95–104. Vrba, J., Robinson, S.E., 2000. Linearly constrained minimum variance beamformers, synthetic aperture magnetometry, and MUSIC in MEG applications. Conference Record of the Thirty-fourth Asilomar Conference on Signals, Systems and Computers vol. 1, pp. 313–317. Wang, H.-C., Chen, L.-F., Shen, H.-P., Lee, P.-S., 2008. Brain smiles while seeing a smiling face: happy mirror neuron system. Annual Meeting of Taiwan Occupational Therapy Association. Taipei, Taiwan. Wendling, F., Ansari-Asl, K., Bartolomei, F., Senhadji, L., 2009. From EEG signals to brain connectivity: a model-based evaluation of interdependence measures. J. Neurosci. Methods 183 (1), 9–18.

17

Wessel, J.R., Ullsperger, M., 2011. Selection of independent components representing event-related brain potentials: a data-driven approach for greater objectivity. Neuroimage 54 (3), 2105–2115. Winter, W.R., Nunez, P.L., Ding, J., Srinivasan, R., 2007. Comparison of the effect of volume conduction on EEG coherence with the effect of field spread on MEG coherence. Stat. Med. 26 (21), 3946–3957. Woolrich, M., Hunt, L., Groves, A., Barnes, G., 2011. MEG beamforming using Bayesian PCA for adaptive data covariance matrix regularization. Neuroimage 57 (4), 1466–1479. Yang, L., Wilke, C., Brinkmann, B., Worrell, G.A., He, B., 2011. Dynamic imaging of ictal oscillations using non-invasive high-resolution EEG. Neuroimage 56 (4), 1908–1917. Zhang, J., Raij, T., Hämäläinen, M., Yao, D., 2013. MEG source localization using invariance of noise space. PLoS One 8 (3), e58408.