## Abstract

Adaptive behaviors are built on the arbitrary linkage of sensory inputs to actions and goals. Although the sensorimotor and associative frontostriatal circuits are known to mediate arbitrary visuomotor mappings, the underlying corticocortico dynamics remain elusive. Here, we take a novel approach exploiting gamma-band neural activity to study the human cortical networks and corticocortical functional connectivity mediating arbitrary visuomotor mapping. Single-trial gamma-power time courses were estimated for all Brodmann areas by combing magnetoencephalographic and MRI data with spectral analysis and beam-forming techniques. Linear correlation and Granger causality analyses were performed to investigate functional connectivity between cortical regions. The performance of visuomotor associations was characterized by an increase in gamma-power and functional connectivity over the sensorimotor and frontoparietal network, in addition to medial prefrontal areas. The superior parietal area played a driving role in the network, exerting Granger causality on the dorsal premotor area. Premotor areas acted as relay from parietal to medial prefrontal cortices, which played a receiving role in the network. Link community analysis further revealed that visuomotor mappings reflect the coordination of multiple subnetworks with strong overlap over motor and frontoparietal areas. We put forward an associative account of the underlying cognitive processes and corticocortical functional connectivity. Overall, our approach and results provide novel perspectives toward a better understanding of how distributed brain activity coordinates adaptive behaviors.

**SIGNIFICANCE STATEMENT** In everyday life, most of our behaviors are based on the arbitrary linkage of sensory information to actions and goals, such as stopping at a red traffic light. Despite their automaticity, such behaviors rely on the activity of a large brain network and elusive interareal functional connectivity. We take a novel approach exploiting noninvasive recordings of human brain activity to reveal the cortical networks and corticocortical functional connectivity mediating visuomotor mappings. Parietal areas were found to play a driving role in the network, whereas premotor areas acted as relays from parietal to medial prefrontal cortices, which played a receiving role. Overall, our approach and results provide novel perspectives toward a better understanding of how distributed brain activity coordinates adaptive behaviors.

- corticocortical coupling
- functional connectivity
- gamma-band neural activity
- Granger causality
- magnetoencephaplography
- visuomotor behaviors

## Introduction

Most cognitive functions arise from the dynamic coordination of neural activity distributed over large-scale brain networks (Bressler and Menon, 2010). Adaptive behaviors and the ability to flexibly choose appropriate actions depending on visual cues and internal goals are no exceptions. A key form of visuomotor guidance is the ability to learn arbitrary and causal relations linking visual inputs to actions and outcomes, named arbitrary visuomotor learning (Wise and Murray, 2000). Performance of arbitrary visuomotor mappings relies on the neural activity of the associative and sensorimotor frontostriatal circuits (Wise et al., 1996; Murray et al., 2000; Passingham et al., 2000; Wise and Murray, 2000; Hadj-Bouziane et al., 2003; Petrides, 2005). However, the interplay between these brain regions is still unclear.

A potential key player for uncovering functional cortical networks is the gamma-band (30–150 Hz) neural activity. Animal studies have observed narrow-band low-gamma oscillations in the 30–80 Hz range during sensory stimulation (Ray and Maunsell, 2010) and cognitive processes, such as attention (Fries et al., 2001; Bichot et al., 2005; Buschman and Miller, 2007; Gregoriou et al., 2009; Bosman et al., 2012) and working memory (Pesaran et al., 2002). Human studies have consistently found correlates of cognitive processes in the high-gamma range (from 60 to 150 Hz) in invasive (Brovelli et al., 2005; Crone et al., 2006; Jerbi et al., 2009; Lachaux et al., 2012; Cheyne and Ferrari, 2013; Ko et al., 2013), noninvasive (Vidal et al., 2006; Ball et al., 2008; Darvas et al., 2010), and multimodal neurophysiological studies (Dalal et al., 2009). Correlations with BOLD responses in animals (Logothetis et al., 2001; Niessing et al., 2005; Goense and Logothetis, 2008) and humans (Lachaux et al., 2007; Nir et al., 2007; Scheeringa et al., 2011; Hermes et al., 2012; Ojemann et al., 2013) further suggest gamma-band activity as a proxy for local processing.

At the large-scale level, binding local activations into coordinated spatiotemporal network activity results from a complex interplay between neuronal dynamics and anatomical connectivity. To quantify the degree of coordination, functional connectivity (FC) measures include various forms of statistical dependences between neural signals, such as linear correlation and Granger causality (Brovelli et al., 2004; Ding et al., 2006; Bressler and Seth, 2011). Patterns of FC can then be analyzed through network-based measures using graph theory (Rubinov and Sporns, 2010). Here, we take a novel approach based on gamma-band neural activity and FC analysis to reveal the human cortical networks and FC mediating arbitrary visuomotor mapping. Single-trial gamma-power time courses were estimated for all Brodmann areas (BAs) by combining magnetoencephalographic (MEG) and structural MRI data with spectral analysis and beam-forming techniques. Performance of visuomotor associations was characterized by an increase in high-gamma activity (HGA) and FC over the sensorimotor and frontoparietal network together with medial prefrontal areas. Graph measures and link community analysis (Ahn et al., 2010) further support that arbitrary visuomotor mappings are mediated by overlapping cortical subnetworks dominated by a core circuit, exchanging Granger interdependences, centered over the dorsal frontoparietal network, with a privileged role for the premotor and medial prefrontal areas.

## Materials and Methods

#### Experimental procedure and data acquisition

##### Experimental conditions and behavioral tasks.

Eleven healthy participants accepted to take part in our study (all were right handed; the average age was ∼23 years; 4 were females and 7 males), gave written informed consent according to established institutional guidelines and local ethics committee, and received monetary compensation (€50). We used an associative visuomotor mapping task, where the relation between visual stimulus and motor response is arbitrary and deterministic (Wise and Murray, 2000). The domain of visuomotor control distinguishes two forms of guidance: standard and nonstandard (or arbitrary) mapping. In standard mapping, such as reaching toward an object, the spatial properties of the visual cue (i.e., target location) provide relevant information for visuomotor transformation. In arbitrary mapping, the visual cue and the action are not spatially congruent; and other properties of the object, such as its form or color, are transformed to plan actions (Petrides, 1987; Wise et al., 1996; Murray et al., 2000; Passingham et al., 2000; Wise and Murray, 2000). To investigate arbitrary visuomotor mapping, we asked participants to perform a finger movement associated to a digit number: digit 1 instructed the execution of the thumb, digit 2 for the index finger, digit 3 for the middle finger, and so on (Fig. 1). Maximal reaction time was 1 s. After a fixed delay of 1 s following the disappearance of the digit number, an outcome image was presented for 1 s and informed the subject whether the response was correct, incorrect, or too late (if the reaction time exceeded 1 s). Incorrect and late trials were excluded from the analysis because they were either absent or very rare (i.e., maximum 2 late trials per session). The next trial started after a variable delay ranging from 2 to 3 s (randomly drawn from a uniform distribution) with the presentation of another visual stimulus (Fig. 1*b*). Each participant performed two sessions of 60 trials each (total of 120 trials). Each session included three digits randomly presented in blocks of three trials.

##### Anatomical, functional, and behavioral data acquisition.

Anatomical MRI images were acquired for each participant using a 3-T whole-body imager equipped with a circular polarized head coil. High-resolution structural T1-weighted anatomical image (inversion-recovery sequence, 1 × 0.75 × 1.22 mm) parallel to the anterior commissure-posterior commissure plane, covering the whole brain, were acquired. MEG recordings were performed using a 248 magnetometers system (4D Neuroimaging, magnes 3600). Visual stimuli were projected using a video projection, and motor responses were acquired using a LUMItouch optical response keypad with five keys. Presentation software was used for stimulus delivery and experimental control during MEG acquisition. Reaction times were computed as the time difference between stimulus onset and motor response. Sampling rate was 2034.5 Hz. Location of the participant's head with respect to the MEG sensors was recorded both at the beginning and end of each session to potentially exclude sessions and/or participants with large head movements. However, none of the participants moved >3 mm during all sessions. Thus, all participants were considered for further analysis.

#### Single-trial HGA at BAs

##### Preprocessing and spectral analysis of MEG signals.

MEG signals were down-sampled to 1 kHz, low-pass filtered to 250 Hz, and segmented into epochs aligned on finger movement (i.e., button press). Epoch segmentation was also performed on stimulus onset and the data from −0.5 and −0.1 s before stimulus presentation was taken as baseline activity for the calculation of the single-trial HGA. Artifact rejection was performed semiautomatically. For each movement-aligned epoch and channel, the MEG signal variance and *z*-value were computed over time and taken as relevant metrics for the identification of artifact epochs. All trials with a variance >1.5 × 10^{−24} across channels were excluded from further analyses. Additional metrics, such as the *z*-score, absolute *z*-score, and range between the minimum and maximum values, were also inspected to detect artifact. Two MEG sensors were excluded from the analysis for all subjects.

Spectral density estimation was performed using a multitaper method based on discrete prolate spheroidal (slepian) sequences (Percival and Walden, 1993; Mitra and Pesaran, 1999). To define the frequency bandwidth of interest for the estimate of single-trial gamma-band activity, we performed time-frequency analyses of the MEG time series at the sensor level. MEG fields were transformed to a planar gradient configuration (i.e., by computing the gradient tangential to the scalp). Time-frequency power estimates were then computed for all sensors over the gamma band from 30 to 150 Hz (in steps of 5 Hz) using 4 orthogonal tapers 0.25 s in duration and 20 Hz of frequency resolution, each stepped every 0.01 s. Time-frequency power estimates were log-transformed and then normalized (i.e., *z*-transformed) with respect to the mean and SD computed over a baseline period from −0.5 to −0.1 s before stimulus onset. Time-frequency maps were averaged across virtual gradiometers and participants. Results showed increases in power with respect to baseline in the high-gamma range from ∼60 to 120 Hz (Fig. 2). Based on the spectral analysis at the sensors level, subsequent analyses focused on the high-gamma band from 60 to 120 Hz. MEG time series were multiplied by *k* orthogonal tapers (*k* = 8) (0.15 s in duration and 60 Hz of frequency resolution, each stepped every 0.005 s), centered at 90 Hz, and Fourier-transformed. Complex-valued estimates of spectral measures *X _{sensor}^{n}* (

*t*,

*k*), including cross-spectral density matrices, were computed at the sensor level for each trial

*n*, time

*t*, and taper

*k*.

##### Source analysis and calculation of HGA.

Source analysis requires a physical forward model or leadfield, which describes the electromagnetic relation between sources and MEG sensors. The leadfield combines the geometric relation of sources (dipoles) and sensors with a model of the conductive medium (i.e., the headmodel). For each participant, we generated a headmodel using a single-shell model constructed from the segmentation of the cortical tissue obtained from individual MRI scans (Nolte, 2003). Leadfields were not normalized. As source location, a 3D grid with regular spacing between the dipole locations of 10 mm was generated for each participant. Individual MRI scans were then warped to the template MRI in MNI space, and the normalization parameters were applied to the dipole grid. Such procedure assures that individual subjects' grid points are located in equivalent brain areas across all subjects according to MNI space. The headmodel, source locations, and the information about MEG sensor position were combined to derive single-participant leadfields.

We used adaptive linear spatial filtering (Veen et al., 1997) to estimate the power at the source level. We used the Dynamical Imaging of Coherent Sources method, a beam-forming algorithm for the tomographic mapping in the frequency domain (Gross et al., 2001), which is well suited for the study of neural oscillatory responses based on single-trial source estimates of band-limited MEG signals (for review, see Hansen et al., 2010). At each source location, Dynamical Imaging of Coherent Sources uses a spatial filter that passes activity from this location with unit gain while maximally suppressing any other activity. The spatial filters were computed on all trials for each time point and session and then applied to single-trial MEG data. Dynamical Imaging of Coherent Sources allows the estimate of complex-value spectral measures at the source level, *X _{source}^{n}* (

*t*,

*k*) =

*A*(

*t*)

*X*(

_{sensor}^{n}*t*,

*k*), where

*A*(

*t*) is the spatial filter that transforms the data from the sensor to source level (for a detailed description of a similar approach, see Hipp et al., 2011). The single-trial high-gamma power at each source location was estimated by multiplying the complex spectral estimates with their complex conjugate, and averaged over tapers

*k*,

*P*(

_{source}^{n}*t*) =

*X*(

_{source}^{n}*t*,

*k*)

*X*(

_{source}^{n}*t*,

*k*)

_{k}

^{*}. Single-trial power estimates aligned on movement and stimulus onset were log-transformed to make the data approximate Gaussian and low-pass filtered at 50 Hz to reduce noise. Single-trial mean power and SD in a time window from −0.5 and −0.1 s before stimulus onset was computed for each source and trial, and used to

*z*-transform single-trial movement-locked power time courses. Similarly, single-trial stimulus-locked power time courses were log-transformed and

*z*-scored with respect to baseline period, so to produce HGAs for the prestimulus period from −1.6 to −0.1 s with respect to stimulation for subsequent FC analysis. Finally, the anatomical position of each source was labeled according to BA using the binary representation of the Talairach–Tournoux atlas (Talairach and Tournoux, 1988) digitized for the Talairach Daemon (Lancaster et al., 2000). Single-trial HGA at each BA was defined as the mean

*z*-transformed power values averaged across all sources within the same BA and cerebral hemisphere (a total of 76 BAs covering both hemispheres). The preprocessing steps, artifact rejection, spectral analyses, and source analysis were performed using FieldTrip toolbox (Oostenveld et al., 2011).

#### Single-trial FC measures between BAs

FC measures characterize statistical dependencies between neural signals, where these dependencies can be undirected (e.g., linear correlation) or directed (e.g., Granger causality). FC analysis provides descriptive statistical measures and makes no assumptions about the underlying structural connectivity. Here, we adopted a comprehensive approach exploiting both undirected (linear correlation and total interdependence) and directed (Granger causality) FC measures for the analysis of statistical dependences between HGA neural signals. FC measures provide relevant information about the coupling and the time-lagged relations only if the temporal relations between signals are preserved during the preprocessing steps. Averaging across trials can destroy induced temporal relations. This is why FC measures were computed on a single-trial basis.

##### Linear correlation.

A conventional undirected connectivity measure is linear correlation. We computed the linear correlation coefficient between the HGAs of pairs of BAs at the single-trial level over a time window ranging from −1 to 0.5 s around movement onset, thus including 300 time samples (i.e., HGA were computed every 0.005 s). Single-trial correlation coefficients were also computed among pairs of HGA that were computed during the prestimulus period, from −1.6 to −0.1 s before stimulus onset. The statistical analyses then searched for significant differences between the movement-related correlations coefficients and those computed in the prestimulus interval.

##### Covariance-based Granger causality measures.

Granger causality (Granger, 1980; Brovelli et al., 2004; Ding et al., 2006; Bressler and Seth, 2011; Seth et al., 2015) is a directed FC measure. Granger (1963, 1980) proposed a statistical criterion to infer causality from process *X* to process *Y* based on the extra knowledge that can be obtained about the future of *Y* given the past of *X*, in a given context *Z*. This criterion can be formalized as the condition of independence as follows:
where the superindex *i* refers to the whole past in a time window of length *T* of a process up to and including sample *i*, *V* refers to the whole system {*X*,*Y*, *Z*}, and *X ^{i}* indicates the exclusion of the past of

*X*. That is,

*X*is Granger noncausal to

*Y*given

*Z*if the above equality holds. This condition of independence constitutes a strong form of Granger noncausality (Granger, 1980), which can generally be tested within an information theoretical framework by calculating the Kullback–Leibler divergence (Cover and Thomas, 2006) of the two probabilities in Equation 1. The Kullback–Leibler divergence is a conditional mutual information and, for a bivariate context comprising only

*X*and

*Y*, it is defined as follows: where the conditional entropies are as follows: and The conditional entropies in Equations 3 and 4 quantify the uncertainty of

*Y*

_{i+1}given the knowledge about the pasts

*Y*only or given the past of both

^{i}*X*and

^{i}*Y*, respectively. Accordingly,

^{i}*F*

_{X→Y}(Eq. 2) compares the uncertainty in

*Y*

_{i+1}when using knowledge of only its own past or both

*X*and

^{i}*Y*. Such conditional mutual information has been studied in different fields, such as communication theory and econometrics, and it has been formulated under different assumptions (e.g., stationarity vs nonstationarity) and terminology (for details, see Chicharro and Ledberg, 2012). In neuroscience, it is best known as the Transfer entropy (Schreiber, 2000). Although Transfer entropy in its general form (Eq. 2) does not require a model of the interaction and it is powerful in capturing any (linear and nonlinear) conditional dependence between

^{i}*Y*

_{i+1}and

*X*(Vicente et al., 2011; Wibral et al., 2014), its calculation is not trivial when high-dimensional probability distributions need to be estimated over short time windows as often encountered in neurophysiological data (for a review of the estimation of information theoretic measures to calculate Transfer entropy, see Hlaváčková-Schindler et al., 2007). To estimate FC measures on short data segments, such as single trials, we implemented an approximation of Transfer entropy that exclusively considers the first term of its expansion (also referred to as the Gram-Charlier expansion as used in Amari et al., 1995), which quantifies the contribution of second-order statistics (i.e., covariance matrix). In other words, our approach is equivalent to implementing a weaker criterion of Granger noncausality, called noncausality in mean (Granger, 1980). This criterion compares the errors of the optimal linear predictors of

^{i}*Y*

_{i+1}, minimizing mean squared error, which corresponds to the conditional variance of the future given the past, so that Equation 1 is reduced to the following: where σ

^{2}(·|·) denotes the conditional variance. The statistical measure used to test Equation 5 is commonly referred to as the Granger causality. In particular, the classical parametric implementation (Granger, 1963) uses autoregressive models to estimate prediction errors. Given that prediction errors can also be estimated from the covariance matrix (Lütkepohl, 2005; Chicharro, 2011), the conditional variances in Equation 5 determine the first term in the expansion of the entropies of Equation 2. This is because, as previously mentioned, the first term of the expansion quantifies the contribution of second-order statistics to the entropy. Put it into practice, given an interval

*T*defining the portion of data

*X*and

^{i}*Y*used to compute covariance measures and a

^{i}*lag*determining the number of time points of the pasts

*X*and

^{i}*Y*used, submatrices of the covariance matrix of the variables comprising the past and the future, {

^{i}*X*

_{i+1}

*,X*

_{i},X_{i-1},…

*X*

_{i-lag+1},

*Y*

_{i+1}

*,Y*

_{i,}

*Y*

_{i-1},…

*Y*

_{i-lag+1}}, are used to estimate the entropies (Cover and Thomas, 2006). For example, where Σ denotes a covariance matrix and |·| denotes the determinant. Within this approach, stationarity is assumed for the calculation of covariance matrices, whereas Gaussianity and linearity are not, and only determine how close is the in mean criterion (Eq. 5) to the strong criterion formalized in Equation 1. That is, for Gaussian variables, Equation 6 constitutes an exact expression for the entropy because higher-order moments are determined by the second-order moments, and Transfer entropy can be considered equivalent to Granger causality (Barnett et al., 2009). Whenever such equivalence does not hold (e.g., for non-Gaussian and nonlinear systems), the covariance-based approximation of the Transfer entropy is always equivalent to the parametric estimation of Granger causality. In addition, it has the advantage of not requiring fitting autoregressive models as an intermediate step. Given this equivalence, we will refer to our estimator as covariance-based Granger causality.

The relation between information-theoretic measures of dependence and Granger causality measures is not restricted to the case of Transfer entropy and linear Granger causality. Measures of Granger causality appear as terms in the decomposition of total interdependence between two time series. It has been formulated in information theoretic (Marko, 1973; Rissanen and Wax, 1987) and linear system (Geweke, 1982, 1984; Chicharro, 2011) frameworks, and applied to neural signals in the spectral domain (Ding et al., 2006). In the time domain and for finite time series, this decomposition is expressed as follows:
and it represents the total interdependence between *X* and *Y* (Geweke, 1982; Chicharro and Ledberg, 2012). *F _{X,Y}* quantifies the dynamic increase of the total interdependence between two time-series at a given point in time, in contrast to the static interdependence quantified by linear correlation, and for an infinite lag converges to the mutual information rate. Such total interdependence is the sum of three Granger causality measures: two directed measures plus the “instantaneous” Granger causality term, accounting for unconsidered common influences to the processes. The instantaneous Granger causality is defined as follows:
and can be expressed in terms of conditional entropies,
The total interdependence is expressed as follows:
The four measures in the decomposition of Equation 7 can be estimated with our covariance-based approach, implying the assumptions and approximations discussed above. Covariance-based Granger causality measures (

*F*

_{X→Y},

*F*

_{Y→X},

*F*

_{X·Y}and the total interdependence

*F*) were computed between pairs of HGAs on a single-trial basis over a time window

_{X,Y}*T*ranging from −1 to 0.5 s around movement onset and for the baseline period, from −1.6 to −0.1 s before stimulus onset, similarly to linear correlation analysis, giving a total of 300 samples in each window. The

*lag*used for the calculation of the covariance matrices was 30 time points (i.e., 15% of

*T*). The MATLAB (The MathWorks) code developed for the calculation of the covariance-based Granger causality measures can be downloaded from the website of the corresponding author.

##### Relation between linear correlation and covariance-based Granger causality measures.

Linear correlation and covariance-based Granger causality measures share common properties. As previously mentioned, linear correlation and total interdependence are symmetric measures quantifying static and dynamic dependencies, respectively. Although these measures are not related by a decomposition analogous to Equation 7, there is a strong relationship between the existence of both types of dependencies. A lack of total interdependence implies a lack of linear correlation; and, if we assume that the future of *X* and *Y* causally depends on their own past, respectively, the opposite relation is also true. This occurs because linear correlation ρ is related to the covariance-based approximation of the mutual information *I(X*_{i+1};*Y*_{i+1}*)*, which is equal to the logarithm of 1 − ρ^{2}, and because conditioning on the past cannot create new dependencies (Chicharro and Panzeri, 2014). It is also clear from Equation 7 that the directed and instantaneous Granger measures are smaller than the total interdependence. Thus, null total interdependence implies the absence of Granger causality measures because they constitute non-negative contributions to the total interdependence. In other words, Granger causality is present if, and only if, both linear correlation and total Granger interdependences are not zero.

#### Numerical simulations

Covariance-based Granger causality measures are not as routinely used as linear correlation analysis in functional neuroimaging. In particular, the relation between the total Granger interdependence with conventional linear correlation analysis and its sensitivity to varying levels of stationarity and noise in the signal is unknown. We therefore performed numerical simulations to compare the accuracy in inferring network FC (i.e., the presence of links between nodes, irrespective of directionality) using three comparable undirected connectivity measures: (1) linear correlation coefficient, (2) partial correlation coefficient, and (3) total Granger interdependence *F _{X,Y}*.

##### Neural mass model (NMM) networks.

Neural signals were simulated using convolution-based neural mass models of local field potentials. Convolution-based NMM consider cortical mesocolumns, rather than a single cell's electrophysiological properties, and they are particularly suited for modeling neurophysiological responses, such as those recorded using the EEG and MEG (David et al., 2006; Moran et al., 2013). The model includes three cell subpopulations: spiny stellate cells in granular layer IV, pyramidal cells, and inhibitory interneurons in extragranular layers (II and III, V and VI). For each node, 13 measured variables include currents and membrane potentials of these three cell subpopulations. A connection between columns A and B corresponds to a link between the pyramidal cell of column A and the spiny stellate cells of column B. Exogenous input can also be included and modeled as an excitation of spiny stellate cells. The generated neural population signals from each node of the network integrate the membrane potentials from three cell subpopulations. We generated NMM signals with two networks depicted in Figure 3*a*. The first network was composed of 5 nodes with unidirectional links drawn in green. The second network was composed of all 10 nodes with unidirectional links shown in red and green (Fig. 3*a*). The goal was to simulate two functional networks associated to baseline (the first network) and event-related (the second) activities, respectively. Using such network functional connectivities, we generated neural signals that were subsequently cut into 50 epochs (trials), each of 150 time points. The number of time points was set to simulate the temporal intervals considered for subsequent HGA analyses. Indeed, the temporal window of interest for FC measures was set 1.5 s (from −1 s to 0.5 s with respect to motor response). HGAs were then expected to contain both ongoing and response-locked modulations, leading to potential latency shifts between areas. To specifically examine the influence of nonstationarities in the network across nodes, such as event-related modulations, we “injected” the networks with an exogenous input via node 1. The input was a Gaussian waveform centered at the 75th sample at each trial (the middle of the trial) and variance equal to 10 samples. The amplitude of the Gaussian input was varied from 0 to 0.04 in steps of 0.005. This injected Gaussian made the data nonstationary. To simulate different levels of stochastic background activity (i.e., noise level), we added Gaussian white noise and varied the amplitude from 0.01 to 0.1 in steps of 0.005. We generated datasets for each network type by varying parametrically the amplitude of the exogenous stimulus and noise level.

##### Statistical inference of FC using linear correlation and Granger causality measures.

For a given dataset with a given input and noise amplitude, we computed three measures of statistical interdependence on a single-trial basis. Given the scope of the numerical simulations, we considered undirected connectivity measures only: (1) linear correlation coefficient, (2) partial correlation coefficient, and (3) total Granger interdependence *F _{X,Y}* (see previous sections). These connectivity measures were computed from NMM signals using different window lengths (

*T*ranged from 20 to 110 in steps of 2 samples). For the total Granger interdependence, the

*lag*was varied from 2 to 0.2 ×

*T*. The maximal

*lag*value was set to 0.2 ×

*T*to avoid small sample size problems in the estimation of covariance matrices, required for the calculation of conditional entropies. Total Granger interdependence was log-transformed to approximate Gaussian distribution.

The goal here was to recover from the simulated neural data the undirected FC (the presence of a link between two nodes, rather than the direction of influence) of the graph depicted in red in Figure 3*a*. This functional network corresponds to the difference between the second and first network. To do so, we performed a paired two-sample *t* test between the connectivity measures computed on the first and second network. More formally, we tested the null hypothesis H_{0}: *X _{ij}*(1) =

*X*(2), where

_{ij}*X*(

_{ij}*w*) is a given connectivity measure between node

*i*and

*j*of network

*w*. The absolute value of the log-transformed

*p*value associated to the

*t*test was tabulated in a sample connectivity matrix of

*n*-by-

*n*, where

*n*is the number of nodes. The sample connectivity matrix was then compared with the known connectivity matrix corresponding to the red network in Figure 3

*a*using a receiver operating characteristic (ROC) analysis and by computing the area under the ROC curve (AUC). The AUC was calculated for all three connectivity measures (linear correlation, partial correlations, and total Granger interdependence) with varying window lengths

*T*(and

*lag*for the total Granger interdependence) on all datasets generated with varying input amplitude and noise levels. To analyze the differences between free parameters and their associated variances, we performed a three-way ANOVA on the AUC computed on the statistical analyses of linear correlation coefficients. The three factors were stimulus amplitude, noise amplitude, and window length

*T*, respectively. All factors showed statistically significant modulations (Table 1). However, the sum of squares of each source (Table 1) shows that the window length is the factor with the largest variance. This suggests that input amplitude and noise level have negligible effect on AUC compared with the effect of window length. We therefore averaged the AUC values across different input amplitude and noise levels and plotted the mean value with associated SEs as a function of window length (Fig. 3

*b*, blue). The graph shows that the mean AUC rapidly increases as window length increases, and it reaches a plateau at ∼70 samples. In general, the AUC is acceptable (>0.95) for

*T*>40 samples. The red curve displays the mean AUC values for the partial correlation coefficient. Here the curve rises more slowly with window length, and it never reaches stable value >95%. These results suggest that linear correlation analysis is more reliable for estimating FC and less sensitive to window length than partial correlation. Finally, we computed the mean AUC averaged over input amplitude and noise levels for the total Granger interdependence analysis. Figure 3

*c*shows the mean AUC values for different values of window length

*T*and

*lag.*Similarly to the previous results, the AUC values increase with longer window lengths. A stable AUC (>0.95) is achieved for window lengths >40 and lags >4 (i.e., >10% of

*T*). These results suggest that the total Granger interdependence can be used reliably for identifying changes in network FC, and it is comparable to linear correlation analysis. However, in addition to measuring the symmetrical dependences, the total Granger interdependence can be decomposed into the directed and instantaneous terms, which can provide directional information of network FC. Overall, the results of the simulations suggest that our method is a good candidate for characterizing cortical FC and may represent an appropriate tool for the analysis of neural data and HGAs.

#### Statistical analysis

Statistical inference of single-trial HGAs and connectivity measures was performed using a linear mixed-effect (LME) model approach. LME models are particularly suited for the analysis of data collected from multiple subjects, where it is important to take into account the variability across participants. The relation between mixed-effect analyses and the two-stage “summary statistics” procedure (Friston et al., 2005) and their combination with permutation tests (Mériaux et al., 2006) has been analyzed in the literature. Briefly, LME models generalize single-subject findings to a population. To do so, they formalize the relation between a response variable and independent variables using both fixed and random effects. Fixed effects model the response variable in terms of explanatory variables as nonrandom quantities. For example, experimental conditions related to population mean may be considered as fixed effects. Random effects are associated with individual experimental units drawn at random from a population, which may correspond to different participants in the study. In other words, whereas fixed effects are constant, random effects are drawn from a prior known distribution. An LME model is generally expressed in matrix formulation as follows:
where *y* is the *n*-by-1 response vector and *n* is the number of observations. *X* is an *n*-by-*p* fixed-effects design matrix and β is the fixed-effects vector of *p*-by-1, where *p* is the number of fixed effects. *Z* is an *n*-by-*q* random-effects design matrix and *b* is a *q*-by-1 random-effects vector, where *q* is the number of random effects; *e* is the *n*-by-1 observation error. The random-effects vector, *b*, and the error vector, *e*, were assumed to be drawn from independent normal distributions. Parameter estimation was performed using the maximum likelihood method.

To test for significant modulations in single-trial HGA and connectivity measures around finger movement with respect to the baseline period, we used a random-intercept and random-slope LME model, which is described by the following:
where *y*(*t*) = [*y _{bl}*(1),

*y*(2), …,

_{bl}*y*(

_{bl}*np*),

*y*(1,

_{mv}*t*),

*y*(2,

_{mv}*t*), …,

*y*(

_{mv}*np*,

*t*)].

*y*(

_{bl}*j*) is a vector containing the baseline neural activity for all trials and sessions (i.e., data from both sessions were concatenated) for subject

*j*= 1,2, …,

*np*where

*np*is the number of participants, at time instant

*t. t*does not refer to trials but time within each trial.

*y*(

_{mv}*j*,

*t*) is a vector including the neural data across all trials and two sessions for subject

*j*at time

*t*with respect to movement onset. The design matrices contain two columns. The first column is a vector of ones to model the intercept; thus, it was eliminated from Equation 12. The second column contains negative ones for baseline trials and ones for event-related trials, therefore modeling the change with respect to baseline, or slope, and it is referred as

*x*and

_{j}*z*in Equation 12. Thus, the first and third terms in the right-hand side of Equation 12 model the intercepts, which correspond to the mean values between baseline and movement-related activity. The second and fourth terms model the slopes, which are the differences between baseline and movement-related activity. The β

_{j}_{1}(

*t*) values are fixed across subjects, whereas the

*b*

_{1j}(

*t*) values model the random variations across subjects. In other words, the parameter β

_{1}(

*t*) models the change in neural activity (i.e., HGA power or FC measures) with respect to baseline at each time point

*t*at the group level; the parameter

*b*

_{1j}(

*t*) models the change in neural activity with respect to baseline for each participant

*j*and therefore explains the across-subjects variability. The across-subjects variability was considered of no interest for the scope of the current analyses. We thus analyzed fixed-effects representative of the entire population. Given the structure of the fixed-effect design matrix, significant differences in movement-related neural activity with respect to baseline can thus be inferred by testing whether β

_{1}coefficients are significantly greater than zero. More formally, the significance of movement-related modulations was inferred using a

*t*test by testing the null hypothesis H

_{0}: β

_{1}≤ 0. Statistical inference was performed for each time point

*t*and each BA for the analysis of HGAs. For the analysis of connectivity measures, statistical inference was performed among pairs of BAs (i.e., no temporal information). To account for the multiple-comparisons problem at the single time-point level, we controlled the false discovery rate (FDR) (Benjamini and Yosef, 1995). The threshold for significance at each time point level was set to

*q*< 0.001. To further assess the validity of our results, we quantified the minimum number of consecutive significant time points required to reject a null hypothesis of absence of a cluster given a chance probability

*p*

_{0}= 0.5 (two possible outcomes, significant or nonsignificant) and kept only those clusters whose duration exceeded a significance level of 0.001. Details of the calculation are given by Smith et al. (2004). The fit of the LME model was performed using the

*fitlme.m*function in Statistical Toolbox of MATLAB (The MathWorks).

#### Characterization of functional cortical networks

##### Hierarchical approach for the characterization of network's nodes and FC.

We defined a functional cortical network as a set of BAs whose HGA and connectivity measures display significant task-related modulations. To assess statistical significance, we used an LME approach (see Statistical analysis). We used random-intercept and random-slope in LME models, where the “intercept” models the average between baseline and movement-related activity, and the “slope” models the difference between baseline and movement-related activity. Significant changes in movement-related HGA and connectivity measures were identified by testing whether the coefficient of the fixed-effect β_{1} was significantly larger than zero. To account for multiple-comparison problems, significance level was set to *q* < 0.001, corrected for the FDR.

Network characterization was performed at four levels of analysis using a hierarchical approach. First, we identified the group of BAs displaying a significant increase in movement-related HGA with respect to mean baseline HGA (averaged from −0.5 to −0.1 s before stimulus onset). Second, among these BAs, we determined which pairs of BAs displayed a significant increase in linear correlation during the movement-related period with respect to baseline. Third, among those significant pairs, we characterized the pairs of BAs exhibiting a significant increase in total Granger interdependence *F _{X}*

_{,}

*. Fourth, among those latter pairs of BAs, we searched for significant increases in directed and instantaneous Granger causality. Such a hierarchical approach has several advantages in terms of exhaustiveness and interpretability of the results. It considers a functional cortical network as composed of links between BAs that show significant effects both in static (i.e., linear correlation) and dynamic dependencies (total interdependence and Granger causality). Such an approach may be especially suited for the interpretation of directed Granger causality measures. In general, even though Granger causality can be used to determine which causal structures are compatible with observed dependencies (e.g., Chicharro and Panzeri, 2014), it does not provide a measure of “true” causality (Pearl, 2009). Therefore, Granger causality is often combined with additional measures that can provide complementary information about the properties of the network under study. If spectral representations of Granger causality are used, the most relevant measures are the magnitude squared coherence among nodes and time-frequency power maps. These can reveal the presence of oscillatory activity (i.e., a peak in the average power spectrum) and the spatial extent (i.e., a peak in the average coherence spectrum). This approach has been used, for example, in Brovelli et al. (2004) where the combined application of power, coherence, and Granger causality measures provided a valuable tool for in-depth investigation of the nature of functional coupling of distributed neuronal assemblies. This approach excludes discriminating potential cases in which significant Granger causality is present among nodes without a clear oscillatory activity and/or coherence. Similarly, in the time domain, as in the current study, we interpret directional Granger causality measures only among nodes previously identified as part of task-dependent functional large-scale network of interest, where BAs display an increase in HGA power, linear correlation, and total Granger interdependence. Furthermore, we checked the robustness of our results and ran all-to-all connectivity analysis. Results showed that the core of the identified functional network is robust to the approach used (results not shown).*

_{Y}##### Detection of functional subnetworks (link communities).

A critical step in the analysis of networks is the detection of communities or groups of nodes potentially corresponding to functional subunits or modules. Communities, however, may overlap (i.e., a given node may belong to more than one group). Large-scale cortical networks may be no exception. For example, a single BA in associative cortical areas, such as BA7 in the superior parietal cortex, may contribute to processing visual inputs in an occipitoparietal network while participating in visuomotor transformations within a frontoparietal network. Link communities (i.e., groups of links rather than nodes) provide an appropriate framework capturing the relationships between overlapping communities while revealing hierarchical organization (Ahn et al., 2010). We adopted such framework for the analysis of adjacency matrices based on linear correlation and Granger causality networks (levels 2, 3, and 4) to identify functional subnetworks. Detection of link communities was based on the calculation of similarity measures between links from undirected and weighted networks. Single-linkage hierarchical clustering was used for the creation of link dendrograms. Link communities were found by cutting the dendrogram at the maximum partition density *D*, which measures the quality of a link partition. The partition density has a single global maximum along the dendrogram in almost all cases because its value is the average density at the top of the dendrogram (a single giant community with every link and node) and it is very small at the bottom of the dendrogram (most communities consist of a single link). Link community analysis was performed using a Python implementation downloaded from https://github.com/bagrow/linkcomm.

## Results

### Behavioral results

We asked participants to perform an arbitrary visuomotor mapping task (Wise and Murray, 2000), where the relation between visual stimulus and motor response is deterministic and highly acquainted (Fig. 1). Each participant performed two sessions of 60 trials each (total of 120 trials). Each session included three digits randomly presented in blocks of three trials. Before MEG acquisition, participants performed 30 trials to familiarize with the task and experimental setup. The average reaction time during MEG recordings was 0.504 ± 0.004 s (mean ± SEM). To quantify potential time-dependent or learning-related modulations of reaction times (RTs) during the execution of the task, we compared the set of RT among participants for the first 50 trials of each session. One-way ANOVA of RTs across trials did not reveal any significant group effect (*F*_{(49,500)} = 0.9, *p* = 0.6638). However, a significant effect was found for session 2 (*F*_{(49,500)} = 1.61, *p* = 0.0071). Pairwise comparison across trials showed that only the first trial was statistically different (α = 0.01, Bonferroni correction, *multcompare.m* function in MATLAB). Therefore, we removed the first correct trial of the second session, in addition to all incorrect and late trials, from subsequent MEG analyses.

### Single-trial high-gamma MEG activity

Time-frequency spectral analysis at the MEG sensor level over the gamma band from 30 to 150 Hz revealed increases in power with respect to baseline in the high-gamma range from ∼60 to 120 Hz (Fig. 2). Tomographic mapping was therefore restricted to the high-gamma band by combining multitaper spectral analysis (Percival and Walden, 1993) and frequency-domain beam-forming algorithm (Gross et al., 2001). Anatomical information about sources was obtained from single-subject MRI and Talairach atlas (Lancaster et al., 2000). Functional and anatomical data were merged to compute single-trial HGA at each BA, defined as the *z*-transformed log-power values with respect to baseline activity averaged across sources within the same BA and cerebral hemisphere, giving a total of 76 BAs covering both hemispheres. Figure 4 shows exemplar single-trial HGA aligned on finger movement for BA 4, comprising the primary motor cortex (left panel), and BA 7, including the precuneus and the superior parietal lobule (right panel), both for the left hemisphere. In this exemplar subject, HGA modulations >3 SDs (*z*-score) are visible on a single-trial basis.

### Visuomotor-related BAs

A functional cortical network is here defined as a set of BAs whose HGA and connectivity measures display significant task-related modulations. We first identified 35 BAs as the network nodes, which displayed significant increases in visuomotor-related HGA with respect to the mean baseline HGA (averaged from −0.5 to −0.1 s before stimulus onset). To assess statistical significance, we used an LME approach, a model particularly suited for the analysis of data collected from multiple subjects. Figure 5 is a statistical map displaying the time course of *t*-values for each BA, grouped in lobes for both hemispheres. The performance of arbitrary visuomotor mappings was associated with a significant increase in HGA over a distributed cortical network covering most of the parietal and frontal areas (Table 2). The largest increase in HGA was observed over the left parietal lobe, both the dorsal (BA 5L and 7L) and lateral areas (BA 39L and 40L), and over sensorimotor (BA 1-2-3 and 4) and premotor regions (BA 6). In prefrontal cortex, the strongest activation was present over the left dorsolateral and dorsomedial cortices (BA 9L) and, bilaterally, over the ventral and dorsal regions of the anterior cingulate area (BA 24 and 32, respectively). Occipital areas showed a significant increase, especially in BA 19, but to a lesser extent due to the alignment of HGA to movement rather than stimulus onset (results not shown). Exemplar time courses of activation of *t* and *p* values at sensorimotor and premotor areas are depicted in Figure 6*a*, *b*. The earliest significant activations arise ∼0.45 s before movement onset (mean reaction time was 0.504 s and SD was 0.118 s, maximal RT was 1 s). If we took the occurrence of the peak of HGA as an index of maximal involvement of a given BA (Table 2, column 4), the earliest activation could be seen between −0.43 and −0.3 s at the border between visual (BA 19L and R, BA 17 and 18R) and temporal areas (BA 37L). Subsequently, parietal areas and frontal areas activated. However, a clear sequential activation is missing. Indeed, the analysis of HGA modulations is not informative about the temporal dependences between BAs. To address this issue, we computed undirected and directed FC measures on a single-trial basis between visuomotor-related BAs.

### FC between visuomotor BAs

We restricted our analysis to BAs displaying a significant increase in HGA with respect to baseline. We first identified the set of pairs of BAs displaying significant increase in linear correlation with respect to baseline (*q* < 0.001, FDR-corrected). Figure 7*a* shows the adjacency or connectivity matrix associated with linear correlation analysis between BAs within the visuomotor-related cortical network (Fig. 5). A significant increase in linear correlation is depicted in the connectivity graph shown in Figure 7*c*. Most of the links were among BAs of the parietal, sensorimotor, premotor, and prefrontal cortices. Node strength was computed as the sum of weights of links connected to each BA, where the weight was the absolute value of the log_{10}-transformed *p* values. Node strength is depicted in color (Fig. 7*c*) and reveals the importance of a given BA in mediating arbitrary visuomotor behavior. BAs with the highest node strength were the superior (BAs 5 and 7) and lateral (BA 40) parietal cortices, the left sensorimotor (BA 1-2-3 and 4), bilateral premotor areas (BA 6), and the bilateral regions of the dorsal anterior cingulate area (BA 24). Several BAs displaying a significant increase in HGA appeared disconnected from the rest of the network (e.g., BA 9L).

Within the set of linearly correlated pairs of BAs, we investigated whether the total Granger interdependence increased with respect to baseline in addition to linear correlation. The total Granger interdependence characterizes the dynamic dependencies, and it quantifies the increase in mutual information between the activity of two BAs at a given point in time with respect to their past, and it is the sum of directed and instantaneous Granger causality measures (see Materials and Methods). Figure 7*b* shows the connectivity matrix for the total Granger interdependence, and Figure 7*d* depicts the links and node strengths. Overall, changes in total Granger interdependence (i.e., dynamic) are more specific than the changes in linear correlation (i.e., static interdependence). Only 15 links between BAs exhibit a significant increase, and they are concentrated over the dorsal frontoparietal network, including the superior parietal cortices, motor (left) and premotor (bilateral) areas, and the dorsal anterior cingulate areas (BA 24). A significant increase is also observed between the left premotor (BA 8L) and anterior cingulate areas (BA 32L). According to node strength, BA 7, 6, and 24 played pivotal roles.

Finally, among those pairs, we searched for significant (*q* < 0.001, FDR-corrected) increases in directed and instantaneous Granger causality (Fig. 7*e*). A significant increase in Granger causality was found from superior parietal (BA 7) to premotor areas (BA 6). In turn, the left premotor areas (BA 6) were found to drive the activity in the anterior cingulate areas bilaterally (BA 24). As highlighted by the node strength, the left premotor cortex (BA 6L) plays a key role in the Granger causality network supporting visuomotor mapping, by acting as a relay node between the parietal and prefrontal areas. No links displayed a significant increase in instantaneous Granger causality.

As a control analysis, we checked whether the changes in FC across areas reflected more than just a change in HGA at the group level. To do so, we analyzed the relationship between changes with respect to baseline in HGA and linear correlation across participants. Given that such across-subjects variability is modeled by the slope-coefficient in the random terms (i.e., *b*_{1}* _{j}*(

*t*) in Eq. 12), we plotted the relationship between the

*t*-values associated with

*b*

_{1}

*(*

_{j}*t*) arising from the linear correlation analysis and the mean

*t*-values averaged over time (from −1 s to 0.5 s around movement onset) for the HGA modulations of the corresponding BAs (Fig. 8

*a*). Each participant is depicted in a different color, and each dot is a given a pair of BAs displaying significant increase in linear correlation. The plot shows that participants displaying high changes in power also show high modulations in linear correlation, and vice-versa. Indeed, the Pearson's correlation between changes in HGA and the linear correlation was highly significant across participants (

*r*= 0.4495,

*p*= 9.33e-59). We then tested whether such relationship was present also at the group level. To do so, we plotted the same relationship for fixed effects (Fig. 8

*b*). In this case, the Pearson's correlation between changes in HGA and FC was

*r*= 0.1483, and it was not significant (

*p*= 0.112). These results therefore suggest that a systematic relationship between HGA modulations and FC measures is present across participants and needs to be taken into account in the statistical analyses. Given that we analyzed fixed-effects only, our results are therefore not contaminated by such confound.

### Link communities of visuomotor networks

A critical step in the analysis of networks is the detection of communities, or groups of nodes potentially corresponding to functional subunits or modules (Girvan and Newman, 2002). Communities, however, may overlap (i.e., a given node may belong to more than one group). Large-scale cortical networks may be no exception. Link communities (i.e., groups of links rather that nodes) provide an appropriate framework capturing the relationships between overlapping communities while revealing hierarchical organization (Ahn et al., 2010). We adopted such framework for the analysis of adjacency matrices based on linear correlation and Granger causality networks to identify functional subnetworks. Figure 9*a* shows the partition density, which measures the quality of a link partition, as a function of dendrogram level for both linear correlation and total Granger interdependence. Link communities were found by cutting the dendrogram at the maximum partition density *D* (Ahn et al., 2010). The number of link communities for the linear correlation and total Granger interdependence was 17 and 5, respectively. For each link community, we computed its strength, defined as the sum of node strengths within each link community with respect to the total node strength of the entire network. Figure 9*b* shows the link community strength and shows that the first link communities account for 67% and 70% of the total strength, whereas the second link communities explain 6% and 12% for linear correlation and total Granger interdependence, respectively. The first link community for the correlation graph is shown in Figure 9*c* and confirms that the dorsal frontoparietal network plays a key role in visuomotor mapping. Node strength is largest for parietal, sensorimotor, and premotor areas. In addition, the parietal (BA 7 and 40L) together with motor cortex (BA 4L) participate in 5 other link communities, suggesting that these BAs participate in different cortical communities (Fig. 10). The first community for the total Granger interdependence is shown in Figure 9*d*, highlighting strong dependencies between the constituents of the dorsal frontoparietal network. The most relevant BAs for such module are the premotor areas BA 6, and the parietal and the left motor areas participate. Overall, community detection analysis suggests that arbitrary visuomotor mapping is mediated by overlapping functional subnetworks of cortical areas, with a core network exerting Granger causal influences centered over the dorsal frontoparietal network, with a privileged role of the premotor and anterior cingulate cortices.

## Discussion

### Cortical network nodes involved in arbitrary visuomotor mapping

The HGA of multiple BAs significantly increased during the execution of familiar stimulus–response associations covering a cortical network primarily including frontal and parietal regions (Fig. 5). In the frontal lobe, the motor (BA 4) and premotor areas (BAs 6 and 8) were found to activate, confirming the role of the sensorimotor circuits in the performance of visuomotor mappings necessary for the instantiation of habitual behaviors (Yin and Knowlton, 2006; Graybiel, 2008; Balleine and O'Doherty, 2010). In the parietal lobe, the somatosensory BA 1-2-3 and parietal areas BA 5, 7, 39, and 40 activated. The involvement of such frontoparietal network, including the superior parietal lobe and the dorsal premotor areas, may support visuomotor computations transforming visual information into motor plans (Wise et al., 1996; Wise and Murray, 2000; Corbetta and Shulman, 2002; Culham and Valyear, 2006). More anteriorly, the left dorsal and dorsomedial prefrontal cortex (BA 9L) and the ventral and dorsal regions of the anterior cingulate area (BA 24 and BA 32, respectively) were also significantly active. Our results suggest the involvement of dorsomedial prefrontal territories of BAs 24 and 32, such as the rostral cingulate zone within the medial frontal areas, which has been described as a crucial node of the human motor system and probably corresponding to the cingulate motor area described in nonhuman primates (Picard and Strick, 1996; Amiez and Petrides, 2014). In such case, the increases in HGA observed in BAs 24 and 32 may correspond to the activation of visuomotor-related neural populations of the rostral cingulate zone. In rodents, however, neural interventions of the infralimbic prefrontal cortex, which may correspond to the ventromedial portions of BAs 24 and 32 (Wallis, 2012), have been found to suppress or block habits (Coutureau and Killcross, 2003; Killcross and Coutureau, 2003; Hitchcott et al., 2007; Smith et al., 2012; Barker et al., 2013; Smith and Graybiel, 2013). Therefore, our results may also be interpreted as confirming the involvement of ventromedial prefrontal regions in familiar stimulus–response behaviors, required for the consolidation of habits. Overall, our results confirm that portions of the medial prefrontal cortex are required for arbitrary visuomotor mappings (Murray et al., 2000).

### Cortical network FC of arbitrary visuomotor mapping

To characterize the FC within the visuomotor-related network, we computed undirected and directed measures between active BAs. The results suggest that visuomotor associations are mediated by an increase in FC between areas of the frontoparietal network and the medial prefrontal cortex. The superior (BAs 5 and 7) and lateral (BA 40) parietal cortices, the left sensorimotor (BA 1-2-3 and 4), bilateral premotor areas (BA 6), and the bilateral regions of the dorsal anterior cingulate area (BA 24) are central, as evidenced by node strengths in correlation and Granger causality graphs. In analogy to the concept of cell assemblies, these BAs (in particular, the parietal areas BA 7 and 40L, and motor cortex BA 4L) participated in multiple functional subnetworks (i.e., link communities), suggesting an overlapping and hierarchical organization of visuomotor-related brain networks. The bilateral superior parietal areas (BA 7) played a driving role in the network, exerting Granger causality on the dorsal premotor areas (BA 6). Such directed functional coupling may be supported by direct corticocortical projections from the superior parietal lobe to premotor areas (Wise et al., 1997). Premotor areas acted as relay nodes passing information from parietal to medial prefrontal cortices (BA 24). BA 24 played a receiving role in the network.

### An associative account of arbitrary visuomotor mapping, cortical networks, and FC

As most forms of visuomotor guidance, arbitrary visuomotor mapping relies on internal representations linking visual cues to motor responses and expected outcomes. Contemporary associative learning theory postulates that acquisition is mediated by goal-directed actions selected according to expected outcomes as well as current goals and motivational state (Rescorla, 1991; Dickinson and Balleine, 1994; Staddon and Cerutti, 2003). Consolidation is considered as the gradual formation of stimulus-driven habitual responses (Dickinson, 1985; Dickinson and Balleine, 1993). Accordingly, performance of arbitrary visuomotor mappings may be viewed as a form of acquainted instrumental behavior. Several studies have opposed arbitrary mapping to instrumental conditioning, to emphasize the conditional nature of associations (Petrides, 1987, 2005; Amiez et al., 2012). However, this apparent contradiction may be resolved if we consider that an arbitrary visuomotor learning task is a form of concurrent and multiple schedule in which two or more reinforcement programs, each signaled by their own discriminative stimuli (i.e., the visual cues), are randomly presented one at a time to the participant who can choose between different operant behaviors (i.e., finger movements) (Bouton, 2006). Indeed, several studies in rodents have shown that conditional T-maze tasks in which turns are instructed by discriminative acoustic tones generate habitual behaviors with overtraining (e.g., Yin and Knowlton, 2006; Graybiel, 2008).

We therefore suggest an associative account of the cognitive processes underlying arbitrary visuomotor mappings, assuming that their performance and consolidation form the basis of highly acquainted instrumental behaviors, such as habits, and represent a special form of instrumental behavior. The selection of a motor response (R) under habitual control is thought to be “triggered” by antecedent stimuli (S), rather than their outcome (Dickinson, 1985). Outcomes, however, can control actions through a form of S–R association in which the stimulus properties of the outcome (O^{S}) can select an associated response (R) (i.e., through an O^{S}–R association) (Balleine and O'Doherty, 2010). The associative-cybernetic model formalizes these notions as the coordination of a habit memory module encoding S–R associations with an associative memory module, representing the belief that a given response causes a particular outcome (Dickinson and Balleine, 1994; Dickinson, 2012). We suggest that the activation of the frontoparietal system and Granger causality from the posterior parietal to premotor and motor areas are associated with S–R visuomotor transformations (i.e., the habit memory module). The posterior parietal cortex and premotor areas may then recruit subcortical areas, such as the sensorimotor territories of the dorsal striatum (Yin and Knowlton, 2006; Graybiel, 2008; Balleine and O'Doherty, 2010; Smith et al., 2012; Smith and Graybiel, 2013), such as the putamen, and finally reach sensorimotor cortical area for implementation of motor commands. In turn, FC from dorsal premotor to medial prefrontal cortex may encode O^{S}–R associations (i.e., the associative memory module) in conjunction with associative and limbic portion of the striatum, such as the anterior caudate nucleus and the ventral striatum, which are known to encode the value of outcomes. Accordingly, even during simple stimulus–response transformations, the dorsal-premotor and medial prefrontal network would monitor correct performance. This interpretation is supported by recent findings suggesting that habits require areas of the associative and limbic frontostriatal loops, in addition to sensorimotor frontostriatal circuits (Coutureau and Killcross, 2003; Killcross and Coutureau, 2003; Smith et al., 2012; Smith and Graybiel, 2013). Having access to O^{S}–R associations, the medial prefrontal network would become central for the detection of motor errors (i.e., unwanted S–R associations), such as in performance monitoring (Botvinick et al., 2004; Ridderinkhof et al., 2004; Rushworth et al., 2007; Bonini et al., 2014) and cognitive control (Miller, 2000; Koechlin et al., 2003; Badre and D'Esposito, 2009). The interaction between the habit and associative memory modules may generate the intention to move (Dickinson, 2012) and may reside in the FC between the posterior parietal cortex and dorsal premotor areas, also mediating the cognitive processes related to motor awareness (Desmurget and Sirigu, 2009). Such system view of the cognitive processes mediating arbitrary mapping arising from cortical FC is depicted in Figure 11.

### FC analysis of HGA

We developed a novel approach to map visuomotor-related BAs and their functional coupling. Covariance-based Granger causality measures provide an intuitive measure linking directed and undirected FC measures and can be used on short neural signals (e.g., single trials) and large networks (e.g., hundreds of nodes). A visuomotor-related cortical network was then defined as a set of BAs whose HGA and FC measures display significant modulations. The use of HGA power modulations as a sole correlate of local neural processing has limitations. A potential limitation in comparing event-related and baseline-related FC measures rests on the fact that significant effects could arise from differences in univariate signal-to-noise ratio, rather than changes in time-lagged interactions among signals. Such potential confounds have been also termed “weak asymmetry,” as opposed to “strong asymmetries” referring to FC modulations arising from changes in time-lagged interactions (Haufe et al., 2013). Even though the simulation and neural results suggest that connectivity changes cannot be explained completely in terms of power changes, in general, we cannot completely discard an influence of the signal-to-noise ratio. Another limitation is the loss of information about phase relations between MEG sensors and, consequently, brain sources. This precludes testing current hypotheses concerning candidate physiological mechanisms mediating the creation of selective communication links between brain regions through phase synchrony (e.g., the “communication-through-coherence” hypothesis) (Fries, 2005). Furthermore, complementary coding schemes between brain regions based on cross-frequency coupling (Lisman and Jensen, 2013; Jensen et al., 2014) cannot be evaluated either. Further work is needed to link visuomotor-related HGA and phase-synchrony modulations in lower bands, such as the α (∼10 Hz) and β (∼20 Hz) range, which are known to play a role in visuomotor behaviors (e.g., Hummel and Gerloff, 2005; Jerbi et al., 2007; Rilk et al., 2011). Finally, given recent evidence of FC patterns undergoing temporal dynamics (Bassett et al., 2011; Zalesky et al., 2014), our study provides novel perspectives toward the study of FC dynamics in task-related context on time-scales relevant for behavioral and cognitive processes.

## Footnotes

This work was supported by Projets exploratoires pluridisciplinaires from the CNRS/Inserm/Inria “Bio-math-Info” entitled “Computational and neurophysiological bases of goal-directed and habit learning.” D.C. was supported by the Autonomous Province of Trento, Call “Grandi Progetti 2012” project “Characterizing and improving brain mechanisms of attention—ATTEND”. V.J. was supported by FP7 BrainScaleS and the James S. McDonell Foundation. We thank Faical Isbaine and Sophie Chen for helping in the MEG experiments; and Driss Boussaoud and Pascal Huguet for useful discussion.

The authors declare no competing financial interests.

- Correspondence should be addressed to Dr. Andrea Brovelli, Institut de Neurosciences de la Timone, UMR 7289 CNRS, Aix Marseille University, Campus de Santé Timone, 27 Bd Jean Moulin, 13385 Marseille, France. andrea.brovelli{at}univ-amu.fr