Asymmetric directed functional connectivity within the frontoparietal motor network during motor imagery and execution

Both imagery and execution of motor control consist of interactions within a neuronal network, including frontal motor-related and posterior parietal regions. To reveal neural representation in the frontoparietal motor network, two approaches have been proposed thus far: one is decoding of actions/modes related to motor control from the spatial pattern of brain activity; and the other is estimating directed functional connectivity (a directed association between two brain regions within motor areas). However, directed connectivity among multiple regions of the frontoparietal motor network during motor imagery (MI) or motor execution (ME) has not been investigated. Here, we attempted to characterize the directed functional connectivity representing the MI and ME conditions. We developed a delayed sequential movement and imagery task to evoke brain activity associated with ME and MI, which can be recorded by functional magnetic resonance imaging. We applied a causal discovery approach, a linear non-Gaussian acyclic causal model, to identify directed functional connectivity among the frontoparietal motor-related brain regions for each condition. We demonstrated higher directed functional connectivity from the contralateral dorsal premotor cortex (dPMC) to the primary motor cortex (M1) in ME than in MI. We further identified significant direct effects of the dPMC and ventral premotor cortex (vPMC) to the parietal regions. In particular, connectivity from the dPMC to the superior parietal lobule (SPL) in the same hemisphere showed significant positive effects across all conditions, while interlateral connectivities from the vPMC to the SPL showed significantly negative effects across all conditions. Finally, we found positive effects from A1 to M1, that is, the audio-motor pathway, in the same hemisphere. These results indicate that the sources of motor command originating in the d/vPMC influenced the M1 and parietal regions for achieving ME and MI. Additionally, sequential sounds may functionally facilitate temporal motor processes.

To evaluate the functional relationship between multiple brain regions associated with MI and ME, it is important to investigate the directed functional connectivity of the frontoparietal network, including both hemispheres ( Hanakawa, 2016 ). Previously, several studies have identified the directed functional connectivity of ME and MI, such as the relationship between SMA and M1 using dynamic causal modeling (DCM) ( Bajaj et al., 2015a ;Kasess et al., 2008 ;Penny et al., 2004 ) and a frontoparietal network using Granger causality mapping (GCM) ( Bajaj et al., 2015b ;Deshpande et al., 2008 ;Gao et al., 2011 ). DCM is an approach to estimate directed functional connectivity based on a neurophysiological model. However, selecting the best model (i.e., connectivity pattern) from all possible models is often infeasible in DCM, especially because of the many brain regions and related factors of interest . GCM is a data-driven method and can be more efficient than DCM. A variant, termed "conditional GCM " (CGCM: Ding et al., 2006 ) was applied to analyze the fMRI data of the ME/MI conditions ( Gao et al., 2011 ;Wang et al., 2019 ); researchers found asymmetry in the directed functional connectivity in right-handed subjects. However, using GCM/CGCM for many brain regions is not very straightforward as the Granger causality measure typically needs to be measured separately for every single connection.
In the frontoparietal network, asymmetry between brain regions that are contralateral and ipsilateral to a manipulating hand is a characteristic aspect of motor-related processing. In general, M1 in both hemispheres is activated during MI and ME, but the region contralateral to a controlled hand shows relatively higher activity than the ipsilateral region because of inhibition from the contralateral M1 to the one ipsilateral to the hand in adults ( Hanakawa et al., 2005 ). Ipsilateral M1, however, also encodes single finger movements during ME, even though the mean activation is lower, as revealed by decoding analysis ( Diedrichsen et al., 2013 ). The asymmetric properties of motor control can appear not only in the mean activation or decoding accuracies, but also in directed functional connectivity. Using conditional GCM, bidirectional edges of directed functional connectivities were extracted in the left-hand task among the L-SPL, L-IPL, and L-dPMC; therefore, the selected edges for the right-handed and left-handed tasks were not the same ( Gao et al., 2011 ). However, the asymmetric network structures of the various motor-associated regions in the frontoparietal network and related networks, in particular in both hemispheres, have been poorly understood.
Here, we investigated the asymmetric network structure related to motor-associated regions represented by the directed functional connectivity during MI and ME. To investigate the network structure, we developed a delayed sequential movement and imagery (dSMI) task by extending the one proposed by Hanakawa et al. (2008) to acquire enough fMRI data samples under ME and MI. We applied a novel causal discovery approach, DirectLiNGAM ( Shimizu et al., 2011 ) to the fMRI time courses in the motor-related volumes of interest (VOIs) in order to estimate their directed functional connectivity on each condi-tion. DirectLiNGAM is an algorithm that can be used to extract causal links based on a linear non-Gaussian acyclic causal model (LiNGAM: Shimizu et al., 2006 ). It is particularly suitable for analyzing large-scale connectivity structures with many VOIs, as compared to DCM/GCM. Finally, we statistically evaluated the directed functional connectivity within frontal regions, between frontoparietal regions, and between audio-motor regions, and illustrated the symmetric/asymmetric structure of the network.

Participants
This study was approved by the Ethical Committee of Advanced Telecommunications Research Institute International and followed the Declaration of Helsinki. A total of 27 adults (three women, mean age = 25.4 ± 5.9 years) participated in this study. Participants received cash remuneration for their involvement in this experiment. All participants had normal or corrected-to-normal vision, no reported history of drugs that affected the central nervous system, and no neurological disease. We assessed the right-handedness of all participants using the Edinburgh Handedness Inventory ( Oldfield, 1971 ). All participants provided written informed consent before the experiment. The data supporting the findings of this study are available from the corresponding author, Takeshi Ogawa, upon reasonable request.

Delayed sequential movement and imagery task
We conducted a dSMI task, which was modified from the original version ( Hanakawa et al., 2008 ) because we needed to collect sufficient trials within a limited fMRI experiment time to optimize the directed functional connectivity with a causal discovery approach. We wrote an in-house program coded for the dSMI task using Psychophysics Toolbox extensions ( Brainard, 1997 ;Kleiner et al., 2007 ;Pelli, 1997 ) in MATLAB (MathWorks Inc., Natick, USA). Briefly, each run incorporated both performance modes (ME, MI, and PL: passive listening) in a pseudorandom manner. Three different classes of visual stimuli (Move: ME; Image: MI; Listen: PL) were alternately presented ( Fig. 1 ). As an instruction stimulus (IS), an Arabic numeral (1-3) was visually presented to indicate the initial finger for the tapping sequence (1 = the pointing finger, 2 = the middle finger, 3 = the ring finger). The initial finger was semi-randomly selected in each trial. After the IS, a cue stimulus (CS) was presented for 2 s. Each participant was instructed to sequentially execute/image tapping following short beep sounds (Move: ME or Image: MI) or to passively listen to the beep sounds (Listen: PL) in one of the modes (ME, MI, and PL). The participants were asked to press the buttons of the response device following the beep sound (duration = 0.2 s; tone = 500 Hz; frequency = 0.667 Hz) sequences in the instruction displaying the Move stimulus (Move trials) during the induction period (IP), and to press the button corresponding to the next finger. In response to the Image stimulus, the participants were asked to image the finger tapping during the IP (7 s) according to the instructed finger as performed by him/her (first-person perspective) without actual finger tapping (Image trials). Then, the participants were required to press the response button that corresponds to the next finger. In the case of Listen stimulus, the participants were asked to respond not only without the movement of their fingers but also without imagining their finger movements during the IP and answer periods. During the IP in the MI and PL conditions, we confirmed the absence of finger movement by the absence of a button response. The Move, Image, and Listen stimuli were pseudo-randomly assigned to the cue stimuli; therefore, the participants could not predict the performance mode in the subsequent trial. The interstimulus interval (ISI) pseudo-randomly lasted for 8, 10, or 12 s. We did not analyze the PL condition in further analyses. Finally, we collected 30 trials in four conditions (right ME: R-ME; right MI: R-MI; left ME: L-ME; left MI: L-MI).  1. Design of the delayed sequential motor imagery and execution (dSMI) task. A participant was instructed to remember the instruction stimuli (IS) that were presented for 2 s (1: pointing finger; 2: middle finger; 3: ring finger) and specified a segment of the cue stimuli (CS), which showed the following conditions: Move, Image, and Listen. Based on the CS, the participants sequentially tapped their fingers (Move: motor execution, ME), or imagined the tapping (Image: motor imagery, MI), or did nothing (Listen: passive listening, PL) during the induction period (IP).

Preprocessing of fMRI data
We preprocessed fMRI data according to previous studies ( Ogawa et al., 2018( Ogawa et al., , 2021. The data were processed using SPM12 software (Wellcome Trust centre for Neuroimaging) as follows: we (i) discarded the first four volumes allowing for T1 equilibration; (ii) corrected the remaining data for slice timing; (iii) realigned the mean image of that sequence to compensate for head motion; (vi) co-registered the structural image to the mean of the functional images and segmented into three tissue classes in the MNI space; (v) normalized the functional images using associated parameters and resampled in a 2 × 2 × 2 mm grid; and vi) spatially smoothed using an isotropic Gaussian kernel of 8 mm full-width at half maximum.

Volumes of interest
We followed a previous study ( Zabicki et al., 2017 ) which used sixteen brain regions and added the bilateral dorsolateral prefrontal cortex (DLPFC) as control regions. Therefore, we defined 20 VOIs ( Table 1 )   were reconstructed from the Destrieux Atlas ( Destrieux et al., 2010 ) in the Wake Forest University Pick Atlas ( Maldjian et al., 2003 ). These masks were resliced to match the preprocessed EPIs with the same spatial resolution. We obtained fMRI time courses within each VOI for each participant. For removal of the spurious variance along with their temporal derivatives, we applied a linear regression including six head motion parameters averaging signals in gray matter, white matter, and cerebrospinal fluid. Additionally, to reduce spurious changes caused by head motion, the data were checked by a method that reduces motion-related artifacts ( Power et al., 2012 ). Finally, we extracted time courses of fMRI signals through a band-pass filter, ranging between 0.008 − 0.1 Hz. Eleven participants were excluded from the fMRI analysis: eight participants did not perform all sessions of the experiment; three participants had excessive head/body movements during the fMRI scans.

DirectLiNGAM
A linear non-Gaussian acyclic causal model (LiNGAM: Shimizu et al., 2006 ) is a causal discovery model that utilizes higher-order distributional statistics to determine directed network connections. The assumptions of LiNGAM are as follows: 1) each node (variable) is linearly related to the other nodes, 2) its residual or external input has non-Gaussian distributions, and 3) the corresponding directed dependency graph is acyclic, i.e., the nodes have a natural ordering. Letting x i denote the i th observation variable ( = 1 , ⋯ , ) , its linear association to the preceding/parent variable can be represented as where is the connection weight, quantifying the direct causal effect, and ( ) is the causal order for each variable. The external influences have non-Gaussian distributions, with non-zero variances, and are independent of each other. Rearranging Eq. (1) into a matrix form, we have where x and e are the vectors of observed variables and external influences , respectively; B = ( ) is an n × n connectivity matrix that can be permuted to a strictly lower triangular matrix (i.e. with all zeros in the main diagonals as well as the upper triangular part). Then, the identifiability of matrix B and the causal ordering/strength among the variables can be rigorously proven ( Shimizu et al., 2006 ) based on the theory of independent component analysis ( Hyvärinen et al., 2001 ). We specifically used the DirectLiNGAM algorithm to estimate the strength of the directed functional connectivity (i.e., nonzero entries of B) among brain regions for each of the MI/ME conditions (see Appendix A). DirectLiNGAM does not require algorithmic parameters and

Fig. 2.
Preprocessing of the input data matrix for DirectLiNGAM. (A) An example of a data matrix in a single VOI for a single subject. The dimension of the data matrix is the number of voxels × 720 samples (6 samples × 30 trials × 4 conditions). (B) Time courses of the principal component analysis (PCA) scores after normalizing the matrix by the mean and variance of the matrix for each VOI to decompose loadings and scores. Red solid, brown solid, red dotted, and brown dotted boxes correspond to the R-ME, R-MI, L-ME, and L-MI conditions, respectively. (C) Input data matrix, x, for the DirectLiNGAM algorithm. x was collected from only the first principal component (PC1) with sign correction using the temporal correlation between the PC1 and a mean time course of the same VOI, concatenated across all subjects for each condition (6 samples × 30 trials × 16 subjects = 2880 samples, i.e., the size of x is 20 VOIs × 2880 samples). is promised to converge to a solution within less repetitive steps if the data strictly follows the assumptions and consists of sufficiently large samples, e.g., > 2000 data points as suggested by previous studies Xu et al., 2014 ). During the optimization process, the algorithm also prunes redundant edges and selects significant edges in an acyclic manner. This is important to concisely visualize a directed acyclic graph (DAG) with multiple regions and to interpret differences across conditions. We used an open-source MATLAB implementation of DirectLiNGAM ( Shimizu et al., 2011 ).
We first prepared a matrix of time courses of the fMRI signal by extracting epochs consisting of 6 samples from the onset of the IP for each VOI (ME and MI conditions contained 30 trials; 6 samples × 30 trials × 4 conditions = 720 samples; matrix size: number of voxels × 720 samples; Fig.2 A). After normalizing the matrix by the mean and variance of the matrix, we conducted a principal component analysis for each VOI to decompose the loadings and scores ( Fig. 2 B). We corrected the sign of the first principal component (PC1) so that it has a positive correlation with the mean time course within the VOI. We collected the corrected PC1 across all VOIs and concatenated the matrix across all subjects for each condition (20 VOIs × 2880 samples [6 samples × 30 trials × 16 subjects], Fig. 2 C). We applied DirectLiNGAM for each condition, and estimated the B matrix, the nonzero entries of which identify the direct causal effects, i.e., directed functional connectivity, between two nodes, as well as the acyclic causal ordering of the nodes.

Statistics of network strength
To evaluate the stability of directed functional connectivity across the subjects, we applied the bootstrapping method to estimate the lower and upper boundary of the strength of the connections with 100 repetitions following a method described by Xu et al. (2014) . We assumed that the fMRI data, x, can be linearly represented by Eq. (2) . To check the validity of the DirectLiNGAM approach, we confirmed the non-Gaussianity of external influences, e = ( I − B )x, of the estimated model using the Kolmogorov-Smirnov test ( p < 0.05, "kstest.m ", in the Statistics toolbox, MATLAB). We then estimated 100 models for each condition (R-ME, R-MI, L-ME, and L-MI) and applied the bootstrapped t -test for each condition ( t -test, p < 0.001 with Bonferroni correction). Additionally, to evaluate the statistical difference of the direct effects between the two conditions, we applied the two-sample t -test ( p < 0.001 corrected using the Bonferroni method).

General linear model
A first-level subject-wise analysis was performed using the general linear model (GLM) in SPM12. We calculated T-statistics for each subject using boxcar regressors, convoluted with a canonical hemodynamic response function without time and dispersion derivatives. We obtained the preprocessed individual EPIs through a high-pass filter, a 128 s cutoff period, to remove the effect of low-frequency drift. Moreover, we regressed out the six parameters of head motion obtained from the realignment procedure. The fMRI data of each session were modeled with two regressors of interest (ME and MI) corresponding to the task trials included in the session. The regressor included all brain activity during the IP (for 7 s) in the analysis. To localize the activated brain regions during the MI and ME conditions, we calculated contrast images includ- We applied the contrast images as the input to grouplevel random-effect analyses with one-sample t -tests. The threshold for statistical significance was set at < 0.001 (uncorrected p-value) and a cluster-based family-wise-error correction of p-value < 0.05.

Directed functional connectivity of fMRI data
Applying the DirectLiNGAM algorithm, we estimated the weight coefficient matrices, B, from data x for each condition (R-ME, R-MI, L-ME, L-MI). By shuffling the subjects' order with 100 repetitions for each condition, we illustrated colormaps of the mean weights over repetitions, such that only significant ones were shown ( Fig. 3 A). According to the colormaps, we drew DAGs to visualize directed effects between two regions ( Fig. 3 B). To indicate graphs of the network, we obtained two indices, such as in-degree (number of edges from other nodes to the node) and out-degree (number of edges from the node), to quantify importance in the network in Fig.3 CD. First, brain regions in the bilateral d/vPMC showed high out-degrees and low in-degrees. This indicates that the bilateral d/vPMC strongly influences other regions, but direct effects from other regions are small (out-degree > in-degree, Fig. 3 C and D). According to the colormaps and DAG, we determined the most dominant connection from the dPMC to M1 contralateral to the hand (L-ME and L-MI: R-dPMC → R-M1; R-ME and R-MI: L-dPMC → L-M1; Fig.3 AB, Supplementary Table 1). We identified positive direct effects from the dPMC to M1 ipsilateral to the hand (R-ME and R-MI: R-dPMC → R-M1; L-ME and L-MI: L-dPMC → L-M1). As a bilateral interaction, we found positive direct effects (all conditions: R-vPMC → L-vPMC, R-SMA →L-SMA; L-ME, R-ME, L-MI: L-dPMC → R-dPMC), and negative direct effects (all conditions: R-vPMC →L-SMA). Additionally, we found strong negative direct effects from the L-dPMC to R-M1 in R-ME, but also weak negative effects in the other three conditions (all conditions: L-dPMC → R-M1). Most positive direct effects in the same hemisphere showed symmetric structure; however, the interhemispheric directed functional connectivity with negative effects showed inhibitions in the target region and an asymmetric structure.
Second, we also found that brain regions in the bilateral A1 regions showed high out-degrees and low in-degrees ( Fig. 3 C and D). In the colormap of the B matrix and DAG, we found a significant positive effect, e.g., in the audio-motor network (L-ME and L-MI: R-A1 → R-M1; R-ME and R-MI: L-A1 → L-M1), and a weak negative connection to the parietal regions. We added A1 into our analysis as a control region based on a previous study ( Zabicki et al., 2017 ), which was not expected motor processes associated with MI and ME. However, our results suggest that A1 may influence not only M1, but also parietal regions similar to the d/vPMC during the dSMI task, because we used a sequential pitch to smoothly lead finger tapping for the participants.
Third, we found that the in-and out-degrees in the parietal regions and the SMA were not obviously different compared to those in the frontal, auditory, and control regions. Connections in the frontoparietal network, in particular from the d/vPMC to SPL/IPS/IPL were identified. We found direct effects from the frontal regions (d/vPMC) to parietal regions (SPL, IPL, IPL), for instance, L-ME (L-dPMC → L-SPL, L-IPS; L-vPMC → L-IPS; R-vPMC → R-IPS), L-MI (L-dPMC → L-IPS; R-vPMC → L-IPS, R-IPS, R-IPL), R-ME (R-vPMC → R-IPL; L-vPMC, L-dPMC → L-IPS), and R-MI (L-dPMC → L-SPL, L-IPS; R-vPMC → R-IPL, R-IPS). Meanwhile, connectivity from the parietal regions to the frontal regions was not selected by DirectLiNGAM. The direct effects from frontal to parietal regions showed an asymmetric structure. In the parietal regions, the bilateral IPS seemed a hub influencing neighboring regions such as the IPL and SPL across all conditions. Thus, to achieve MI and ME, direct effects may show interactions within the frontoparietal network and intraparietal network.
Regarding the usage of bilateral FMG and DLPFC as control regions, the profiles of their in-/out-degrees were the opposite to those of the d/vPMC and A1, such that the in-degrees of these regions were higher than the out-degrees. Additionally, directed connections from the FMG and DLPFC to the frontal motor regions, parietal regions, and A1 were not selected by DirectLiNGAM ( Fig.3 A) across all conditions. This result indicates that brain activity in the FMG and DLPFC receives broad input from many regions in the dSMI task, but in a possibility different process from the motor control.

Differences in direct effects due to modes and laterality
To identify the direct effects between the two conditions, we evaluated the differences between two samples (A: L-ME vs. L-MI, B: R-ME vs. R-MI, C: L-ME vs. R-ME, and D: L-ME vs. L-MI in Fig.4 ; we show only significant edges thresholded at p < 0.001, a two-sample t -test, corrected by the Bonferroni method). Each element showed a difference in the mean direct effect between two conditions if the difference was statistically significant. First, we found higher effects from the dPMC to M1 contralateral to the hand (green dotted boxes in Fig. 4 C and D). Second, we found the frontoparietal connectivity: L-dPMC → R-IPS and R-vPMC → R-IPS in L-ME, and R-dPMC → R-IPS and L-vPMC → R-SPL in R-ME (orange dotted boxes).
Obvious differences in directed functional connectivity are represented in the edges dPMC → M1 and A1 → M1, contralateral to the hand ( Fig. 4 A and B). The direct effects of dPMC → M1 were more strongly affected in ME than in MI. This edge is essential to achieve motor control and related imagery. Second, the strength of the edges (A1 → M1) was also higher in ME than in MI. One possibility is that the dSMI task required finger tapping followed by the beep sound which may enhance the audiomotor process in the ME condition.
Comparing the laterality of ME ( Fig. 4 C) and MI ( Fig. 4 D), the edges dPMC → M1 contralateral to the hand clearly showed differences in laterality. In comparison between ME conditions, we identified edges L-M1 → R-IPL, R-dPMC → R-IPS larger in R-ME, and edges R-IPS → R-IPL in L-ME ( Fig.4 C). In the case of MI, we identified edges L-M1 → L-SPL, R-M1 → L-FMG, and R-IPL → R-DLPFC to be larger in R-MI, and edges L-IPL → L-DLPFC and R-A1 → R-M1 to be larger in L-MI ( Fig. 4 D).

Fig. 3.
Direct acyclic graphs and matrices of directed functional connectivity estimated by DirectLiNGAM. The colormaps (A) and direct acyclic graphs (B) show the mean strengths of directed functional connectivity, B, between two brain regions. Significant edges were connected between two VOIs (bootstrapping with 100 repetitions, threshold of p < 0.001 corrected by the Bonferroni method). VOIs were divided into three parts, frontal motor regions (solid black lines and dots in B), parietal regions (solid pink lines and dots) and control regions (solid green lines and dots) in (A) and (B). The width of the edge in (B) shows the absolute value of the strength of the directed functional connectivity and an arrow of the edge is the direction. Positive weights are colored in red and negative weights are colored in blue. For visualization purposes, a colormap was set between − 1 to 1. (C) and (D): out-/in-degrees were counted with respect to significant functional directed connectivity for each VOI. Bars are colored in blue for R-ME, orange for L-ME, yellow for R-ME, and purple for L-MI. VOIs: volumes of interest; R-ME: right motor execution; L-ME: left motor execution; L-MI: left motor imagery.

Results of the general linear model of event-related activity
To confirm brain activations between left and right ME/MI, we performed GLM analysis to compare contrasts. First, we compared contrasts between R-ME and L-ME; the results are summarized in Table 2 (activation maps are illustrated in Supplementary Fig. 1 C, E, G, and I). The brain regions that emerged from a contrast of R-ME > L -ME included the postcentral gyrus (M1), insula including the Rolandic operculum (RO) and Heschl gyrus (A1), putamen, SMA, and thalamus in the left hemisphere contralateral to the moving hand. In the opposite contrast, L-ME > R -ME, we identified clusters in the precentral gyrus, thalamus, and putamen in the right hemisphere contralateral to the moving hand. In the case of MI, we found similar tendencies (R-MI > L -MI: L postcentral gyrus [L-M1], L insula; L-MI > R -MI: R postcentral gyrus [L-M1]). These results were consistent with previous findings that brain activity in the contralateral hemisphere is predominantly higher than that in the ipsilateral hemisphere. Therefore, the clusters were asymmetrically colocalized in comparisons between lateralities ( Hanakawa et al., 2008 ).
Next, we compared contrasts between ME and MI such as modes for each hand ( Table 3 ; activation maps are illustrated in Supplementary  Fig. 1 D, F, H and J). The brain regions that emerged from a contrast of R-ME > R -MI included the left postcentral gyrus (L-M1), which means that L-M1 is more activated at R-ME than at R-MI. Similarly, we identified clusters with a contrast of L-ME > L -MI including those in the right precentral (R-M1) and right thalamus. On the other hand, opposite contrasts (R-MI > R -ME, L-MI > L -ME) showed clusters in the prefrontal (L IFG Tri: left inferior frontal gyrus, triangle part; L MFG: left middle frontal gyrus; R MFG), parietal (L IPL; L SPL), and occipital (L SOG: left superior occipital gyrus; L MOG: left middle occipital gyrus; R MOG; R MFG 2) regions. Brain activity corresponding to ME in M1 was higher than that of MI, but brain activity corresponding MI was lateralized in Fig. 4. Differences of the mean connectivity strengths between two conditions represented in colormaps. From two samples estimated from bootstrapping, we illustrated significantly different edges with a two-sample t -test ( p < 0.001 corrected using the Bonferroni method). In (A) L-ME vs L-MI, (B) R-ME vs R-MI, it shows differences between ME and MI. (A) and(B) colormaps show ME as blue and MI as red. In (C) R-ME vs L-ME, (D) R-MI vs L-MI, it shows differences between R and L.

Discussion
In this study, we demonstrated the asymmetric representation of the direct functional connectivity during ME and MI. By applying Di-rectLiNGAM and GLM to the fMRI data, we characterized functional asymmetry in the d/vPMC and A1, parietal regions, and prefrontal control regions. We found higher direct functional connectivity from the contralateral dPMC to M1 compared to that from the ipsilateral dPMC to M1, as well as activations in M1. Because of the beep sounds guiding the finger tapping, A1 showed high out-degrees that predominantly affected M1 as well as d/vPMC. Regarding the parietal regions, strong connections were observed mainly from the d/vPMC, as well as intraconnections. The frontal control regions were affected by other regions, but not vice versa. However, higher activations during MI were distributed in the left hemisphere than ME in both cases of the right and left hand, such that MI showed functional asymmetry. Here, we would like to discuss the network structure with directed functional connectivity and its brain activations across hands and modes.

Functional roles of d/vPMC evidenced by activation patterns and causal discovery during ME and MI
Our results of the GLM in Tables 2 and 3 showed that the brain activity in motor regions is generally contralateral to the movement/imagery side, and that brain activity in M1 during ME was higher than that during MI. These results were consistent with the findings of previous neuroimaging studies ( Kawashima et al., 1993 ;Hanakawa et al., 2003;Gao et al., 2011 ). From the viewpoint of causal discovery analysis, strong positive direct effects (dPMC → M1) in the same hemisphere contralateral to the tapping hand during MI and ME were selected ( Fig. 3 ). Previously, positive directed functional connectivity has been found from the dPMC/SMA to M1 (CGCM: Gao et al., 2011, Wang et al., 2019DCM: Kasess et al., 2008 ), therefore, the positive directed functional connectivity shows active inference from a source region to a target region. Gao et al. demonstrated that the in-out degrees (difference between in-and out-degrees) of the M1 contralateral to the performing hand were higher than those with the right or left hand, regardless of ME or MI tasks. However, our results demonstrated a relatively higher outdegree than in-degree and were inconsistent with Gao et al., 2011 . Our approach dealt with more nodes, not like CGCM/DCM, and might apply a more specific parcellation than the previous method. Consequently, the dPMC may be a core to generate neural representation associated with ME/MI and broadly affect not just the M1 but also the parietal and frontal regions.
We also found positive effects (dPMC → M1) ipsilateral to the tapping hand (Supplementary Table 1). Our task requires the participants to tap or imagine tapping with one hand; therefore, we did not expect to identify the ipsilateral connections between the dPMC and M1 in this analysis. Regarding comparisons of the right and left hemispheres, there were higher contralateral connections than ipsilateral ones in both ME The MNI coordinate corresponds to the peak voxels within each cluster. Clusters were set at a threshold of p < 0.001 and the cluster level family-wise-error was set at p < 0.05. Note: and MI. In a previous study by Gao et al. (2011) , functional connectivity was estimated by CGCM during ME/MI of the right/left hand similar to our task; however, they did not include ipsilateral M1 into the model. Wang et al. (2019) also estimated the functional connectivity of five regions using CGCM during ME and MI in young and older subjects; they found bilateral connections between the PMC and M1. We successfully identified the ipsilateral connections (dPMC → M1) despite the high dimensional network compared to previous studies. This suggests that neural representations in the ipsilateral pathway may partly help to achieve ME and MI.

Audio-motor network during ME and MI
Regarding the audio-motor network, we found (i) high out-degrees and low in-degrees, (ii) positive effects (A1 → M1) bilaterally, and (iii) negative effects (A1 → parietal regions), in particular, connections within the same hemisphere. We did not expect to find connections related to A1 when we designed this study because we assumed that A1 was a control region. The decoding accuracies of the modality (ME vs. MI) in A1 obtained by multi-voxel pattern analysis were significant in a study by Zabicki et al. (2017) . This suggests that the voxel pattern in A1 contains information that can distinguish between ME and MI. Lima et al. (2016) reviewed the roles of the SMA and pre-SMA in auditory processing and auditory imagery and concluded that activating sound-related motor representation might facilitate behavioral responses to sounds such as 'when' ventral pathway. On the other hand, Morillon et al. (2015) found poorly established anatomical/functional connections between auditory and premotor/motor areas, either weak or absent in primates. Our results indicate that there are direct effects on the audio-motor pathway, to M1 bilaterally and positively, and to the parietal regions negatively in the same hemisphere. Additionally, we found negative effects of L-A1 →L-SMA in all conditions, which means that A1 did not facilitate the motor responses of the SMA. In our task, the participants were instructed to follow the beep sounds to tap their finger or imagine finger tapping, without the inclusion of spatial information in sounds. Periodic sounds associated with motor control may facilitate the preparation of the next finger movement following the sound sequence in M1 and inhibit spatial audio-motor processes, such as the 'where' dorsal pathway ( Morillon et al., 2015 ), according to negative effects from A1 to the parietal regions observed in L-MI and R-MI ( Fig. 3 A).

Functional asymmetry of the parietal regions and left DLPFC
Regarding networks associated with the parietal region, we found functional asymmetry, frontoparietal connections (negative effects: bilateral-vPMC → bilateral SPL; positive/negative effects: bilateral dPMC → bilateral SPL, L-IPS; Fig. 3 ), and that the left parietal regions were highly activated in MI compared to ME ( Supplementary  Fig. 1). In a previous study using CGCM, bidirectional directed functional connectivity between the dPMC and parietal regions (SPL, IPL) in the same hemisphere was identified from fMRI data during ME and MI ( Gao et al., 2011 ). There are functional similarities and differences in PMC across the dorsal-ventral axis, such as the prediction of sequential finger movement ( Nambu et al., 2015 ), representation of chucking sequential finger movement ( Yokoi and Diedrichsen, 2019 ), encoding of spatial location ( Gallivan et al., 2011 ), and linguistic motor planning of the fingers ( Hanakawa et al., 2005 ). Particularly, a specific role of the vPMC is to encode target locations depending on head/eye-centered coordinates and arm/limb-centered coordinates ( Gallivan et al., 2011 ). Hanakawa et al. (2005) illustrated neural representations in the vPMC to mediate translation from the linguistic stimulus to sequential mo-tor planning. Regarding the dSMI task, our results illustrated that the L-vPMC inhibited a spatial motor process in the R-SPL, and that the bilateral dPMC facilitated the motor process of the parietal regions in the same hemisphere. This suggests that the vPMC is functionally associated with encoding effectors and spatial coordination, and that left parietal regions play a role in the sequential motor processing relative to right parietal regions. Additionally, we found common connections within the parietal regions (L-SPL → L-IPL; R-SPL → L-IPL; L-IPL → L-SPL, R-IPS, L-IPL; R-IPS → R-SPL, R-IPL), and higher brain activation in the left parietal regions during MI than ME ( Tables 2 and 3 ). L-IPS had positive effects on the L-SPL, R-IPS, and L-IPL across all conditions. Interestingly, these results indicate the functional asymmetry of the parietal regions. This suggests that the SPL may be associated with the process of somatosensory modality and visual information and its integration. The parietal regions are essential for encoding the spatial processing of target or effectorspecific coordinates, such that neural representation of gestures/actions with the hand/foot are predicted by classification methods (three righthand actions: aiming, extension-flexion, and squeezing, Pilgramm et al., 2016 ; three right-hand actions and their imagery, Zabicki et al., 2017 ; touching left/right and gazing left/right, Gallivan et al., 2011 ). Our results indicated that the activity in the left parietal regions occurred by MI/ME not only with the right hand but also with the left hand.
Regarding the left DLPFC, we identified clusters satisfying the condition MI > ME ( Table 3 ). In early neuroimaging studies, functional asymmetry was found in the M1 and left prefrontal regions ( Kawashima et al., 1993 ), precuneus, and middle frontal gyrus (Hanakawa et al., 2003). The results of DirectLiNGAM showed that the DLPFC and FMG had high in-degree and low out-degree with weak connections, indicating that they may collect sensory information. In particular, MI induced not only motor-related activity in the frontoparietal network, but also activity in the prefrontal regions which integrates information from the body and the environment and participates in higher-order motor control ( van der Meulen et al., 2014 ). It is likely that the frontal regions do not affect the motor process, but rather play a role in monitoring motor processes during ME and MI.

Limitations of DirectLiNGAM and the dSMI task
The DirectLiNGAM algorithm is a suitable approach to investigate causal relationships among multiple nodes ( Shimizu et al., 2011 ). The advantages of this algorithm are as follows: i) a causal structure can be fully and uniquely identified, provided all assumptions hold with the sufficient observation sample size, ii) no requirement of prior knowledge for the network structure, and iii) much less complexity of computation comparing to DCM which requires plenty of conditional independence tests . Conventional methods such as DCM or GCM require large computational costs in treating ≥ 4 time series ( Zhang et al., 2008 ;Wang et al., 2019 ). However, DirectLiNGAM is assumed to be DAG and to prune weakly contributed edges with sparseness; therefore, there is a possibility of missing bidirectional directed functional connectivity. Moreover, the LiNGAM algorithm does not consider la time series like GCM. In this study, comparisons of a cyclic model or a LiNGAM model with hidden variables/confounders/lags were beyond the study's scope. In the future, it is necessary to evaluate how the different bidirectional effects, latent confounders, or time lags are represented in GCM and DirectLiNGAM. We proposed the dSMI task based on a previous study ( Hanakawa et al., 2008 ). In our task, the IP duration was fixed at 7 s; therefore, the beep sound rang 10 times. This means that participants could notice answers when they saw the IS in the MI condition. Before the fMRI experiment, we instructed the participants not to immediately predict an answer when the trial started. Unfortunately, there was no way to confirm that the participants performed correct sequential finger imagery in the MI condition. One possibility is to decode brain signal patterns during MI to control a brain machine interface/brain computer interface (BMI/BCI) system. However, this was beyond the study's scope and might be proved in the field of BMI/BCI in the future.

Conclusions
By applying the causal discovery approach to the fMRI data, we found that DAG represented the frontoparietal motor network during ME and MI. To support the estimation of the directed functional connectivity, we also confirmed spatial brain activity with GLM during ME and MI. We developed a dSMI task based on that described by Hanakawa et al. (2008) to acquire enough fMRI data samples under the ME and MI conditions. We found predominant directed functional connectivity (d/vPMC → M1), which is a core of the motor control. Additionally, we identified audio-motor connections from A1. The parietal regions and DLPFC showed functional asymmetry during ME and MI. The parietal regions may not only receive movement-related information from the d/vPMC but may also communicate with neighboring regions. MI can cause a higher cognitive load than ME. Our approach has the limitation of ignoring bidirectional interactions and the assumption of independence between variables. It is necessary to compare several methods of the representation of causalities of multiple brain regions. Our findings might help to understand large-scale interactions across bilateral frontoparietal networks during ME and MI.

Data and code availability
The data that support the findings of this study are available from the corresponding author, Takeshi Ogawa, upon reasonable request.

Author contributions
Takeshi Ogawa designed this study and performed data acquisition. Takeshi Ogawa and Hideki Shimobayashi analyzed the data. Takeshi Ogawa, Hideki Shimobayashi, Jun-Ichiro Hirayama, and Motoaki Kawanabe discussed data analysis and results. Takeshi Ogawa and Jun-Ichiro Hirayama wrote the original draft.

Declarations of Competing Interest
None.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2021.118794 .

DirectLiNGAM algorithm
We briefly review the DirectLiNGAM algorithm ( Shimizu et al., 2011 ) for self-consistency. The algorithm sequentially finds an exogenous variable (i.e. with no incoming edges from the other ) and then subtract its effects from the other variables; thus, it recursively determines the causal ordering and connection strengths of the LiNGAM model. It relies on the fact that is exogenous if and only if is independent of its residuals ( ) = − ; when is regressed on for all ≠ . In practice, an exogenous variable is therefore sought by minimizing the sum of pairwise independence measures evaluated between a variable and all its residuals ( ) . The following statistic is typically used: where denotes the set of the subscripts of variables , i.e., = { 1 ⋯ } , and ̂ ( , ( ) ) denotes a kernel-based estimate of mutual information ( Bach and Jordan, 2003 ). The entire algorithm is summarized next, which we duplicate from the original paper: Direct LiNGAM algorithm ( Shimizu et al., 2011 ) (1) Given a -dimensional random vector x, a set of its variable subscripts U and a × data matrix of the random vector as X, initialize an ordered list of variables ∶= ∅ and ∶= 1 .
(2) Repeat until − 1 subscripts are appended to : a perform least squares regressions of on for all ∈ ∖ ( ≠ ) and compute the residual vectors r ( j ) and the residual data matrix R ( j ) from the data matrix X for all ∈ ∖ . Find a variable that is most independent of its residuals: (1) Append the remaining variable to the end of .
(2) Construct a strictly lower triangular matrix B by following the order in , an estimate the connection weights by using some conventional covariance-based regression such as least squares and maximum likelihood approaches on the original random vector x and the original data matrix X.