 Methodology
 Open Access
 Published:
Identifying bidirectional total and nonlinear information flow in functional corticomuscular coupling during a dorsiflexion task: a pilot study
Journal of NeuroEngineering and Rehabilitation volume 18, Article number: 74 (2021)
Abstract
Background
The key challenge to constructing functional corticomuscular coupling (FCMC) is to accurately identify the direction and strength of the information flow between scalp electroencephalography (EEG) and surface electromyography (SEMG). Traditional TE and TDMI methods have difficulty in identifying the information interaction for short time series as they tend to rely on long and stable data, so we propose a timedelayed maximal information coefficient (TDMIC) method. With this method, we aim to investigate the directional specificity of bidirectional total and nonlinear information flow on FCMC, and to explore the neural mechanisms underlying motor dysfunction in stroke patients.
Methods
We introduced a timedelayed parameter in the maximal information coefficient to capture the direction of information interaction between two time series. We employed the linear and nonlinear system model based on short data to verify the validity of our algorithm. We then used the TDMIC method to study the characteristics of total and nonlinear information flow in FCMC during a dorsiflexion task for healthy controls and stroke patients.
Results
The simulation results showed that the TDMIC method can better detect the direction of information interaction compared with TE and TDMI methods. For healthy controls, the beta band (14–30 Hz) had higher information flow in FCMC than the gamma band (31–45 Hz). Furthermore, the betaband total and nonlinear information flow in the descending direction (EEG to EMG) was significantly higher than that in the ascending direction (EMG to EEG), whereas in the gamma band the ascending direction had significantly higher information flow than the descending direction. Additionally, we found that the strong bidirectional information flow mainly acted on Cz, C3, CP3, P3 and CPz. Compared to controls, both the betaand gammaband bidirectional total and nonlinear information flows of the stroke group were significantly weaker. There is no significant difference in the direction of beta and gammaband information flow in stroke group.
Conclusions
The proposed method could effectively identify the information interaction between short time series. According to our experiment, the beta band mainly passes downward motor control information while the gamma band features upward sensory feedback information delivery. Our observation demonstrate that the center and contralateral sensorimotor cortex play a major role in lower limb motor control. The study further demonstrates that brain damage caused by stroke disrupts the bidirectional information interaction between cortex and effector muscles in the sensorimotor system, leading to motor dysfunction.
Background
In the process of human voluntary movement, the motor cortex of the brain sends out instructions to control muscle actions through the motor nerve pathway, and the sensory information of the muscle is fed back to the cortex through the sensory nerve pathway to ensure the accurate execution of the action [1,2,3]. This information interaction can be quantified by the coupling relationship between the EEG signal and the SEMG signal of the effector muscle with the development of noninvasive, hightime resolution scalp EEG acquisition technology. Therefore, the functional corticomuscular coupling (FCMC) has become an important way to reveal the controlfeedback mechanism of the nervous system and to evaluate the motor function and rehabilitation effect on patients with neurological diseases such as stroke [3,4,5,6].
The key challenge to constructing the interrelationship between complex neurophysiological signals is to accurately capture the information flow between the signals. More specifically, it includes two important indicators: direction and strength. Studies were conducted on FCMC from these two aspects. The coherence method is one of the main methods to quantify the functional coupling between the cerebral motor cortex and the effector muscle [7, 8]. However, previous studies confirmed that the information interaction between the motor cortex and the effector muscles was directional [3, 9]. The lack of ability to identify the direction of information interaction limits the application of the coherence method in the analysis of FCMC. Granger causality (GC) and its extension methods used as directional methods to measure the causal relationship between time series have been applied in the analysis of FCMC [2, 3, 10]. The GC method is based on a linear autoregressive model, and its statistical nature is a prediction of stationary timeseries data [11]. However, neurophysiological signals have proved to be nonlinear [12,13,14].Therefore, the effectiveness of GC in analyzing the relationship between nonlinear neurophysiological signals is also questioned [4, 15]. Modelfree methods have been used in recent years to analyze the information interaction between neurophysiological signals so as to address the challenge of nonlinearity. The commonly methods are mutual information (MI) and transfer entropy (TE).
The MI method evaluates the interaction relationship between two random variables X and Y by measuring the shared information between them [16]. MI can detect the linear and nonlinear statistical correlations between two signals, and therefore is widely used in the field of neuroscience [17,18,19]. However, MI is a symmetrical measurement method that cannot determine the direction of information flow. To solve this problem, Vastano et al. proposed timedelay mutual information (TDMI) to detect information transmission in spatiotemporal systems [20]. TDMI was then introduced into the analysis of information transmission between neurophysiological signals [6, 21, 22]. Nonetheless, it is difficult for both MI and TDMI to accurately estimate the probability density function (PDF) and joint probability density function (JPDF) in the calculation process for short and complex time series [23]. On the other hand, different estimation methods also directly affect the accurate establishment of the relationship between signals.
The TE method is also a modelfree method based on information entropy, with the ability to detect linear and nonlinear coupling [24]. Benefited by its asymmetric and transition probability calculation characteristics, TE was considered to be an effective method for detecting causality between neurophysiological signals in recent years [25, 26]. Unfortunately, TE cannot accurately detect the coupling in practical applications when the time series is not long enough [4, 21].
Reshef et al. proposed the maximal information coefficient (MIC) method in 2011 [27]. The generality attribute of MIC meets the requirements of measuring different functional relationships; its equitability attribute ensures that different functional relationships obtain similar measured values at the same noise level. In particular, Reshelf et al. were the first to propose a formula to calculate the nonlinear components of the relationship between two variables, that is, MICρ^{2}, where ρ represents the Pearson correlation coefficient. On the contrary, neither MI nor TE can identify pure nonlinear coupling because the results include both linear and nonlinear coupling. Due to the aforementioned advantages, MIC was widely used in the field of neuroscience [28,29,30]. In our previous study, MIC was first applied to the analysis of linear and nonlinear coupling components in FCMC [31]. However, limited by the symmetry of mutual information, MIC is also symmetric, that is, MIC(X, Y) = MIC(Y, X), so MIC fails to identify the direction of information interaction between signals. To our knowledge, none of the above MICbased studies analyzed directional specificity in information interaction.
To overcome this limitation, a timedelayed maximal information coefficient (TDMIC) method was proposed in this study by introducing a timedelay parameter to capture the information transmission delay between two short time series. The algorithm was first tested with simulated data to verify the effectiveness of this method. Linear and nonlinear systems with short data lengths were constructed to compare the performance of TDMIC, TDMI, and TE (kernel estimator) in identifying the direction of information flow. As an application of experimental data, the TDMIC method was applied to explore the directional specificity of total and nonlinear information flow of healthy controls and stroke patients in FCMC in a specific frequency band. This study provided a new perspective for exploring the characteristics of FCMC.
Methods
Timedelayed maximal information coefficient
For the finite data set D of ordered pairs, the data points {x; y} were distributed in a twodimensional space, and the data space was divided into xbyy grids. In this case, the MI of the two variables was expressed as:
where p(x, y) is the JPDF of time series X and Y, and p(x) and p(y) are the marginal PDF of X and Y, respectively. PDF and JPDF were obtained by calculating the probability of data points in D that fell into each grid.
When the number of grids xbyy was fixed, different grid division methods were used. The maximum value of MI among all grid division methods was determined:
To facilitate comparison across grids with different dimensions, \(I^{*} (D,{\text{x}},y)\) was normalized by log min{x; y}, and then the characteristic matrix M of a set of data D was defined as follows:
After all elements in the matrix M were normalized, the score range obtained was between 0 and 1.
For the data set D of ordered pairs with sample size n, the MIC was defined as the maximum value of the characteristic matrix obtained by all grid partitioning:
where the grid size xbyy was limited with \(B(n)^{{(B(n) = n^{0.6} )}}\) to reduce the calculation efforts. The range of the MIC value was [0, 1]; the higher the score, the stronger the correlation between the two variables.
In addition, RESHEIF et al. defined a natural measure of nonlinearity as follows [27]:
where ρ denotes the Pearson productmoment correlation coefficient. When the NL value was greater than 0, it indicated a nonlinear relationship.
MIC had the characteristic of symmetry, that is, for variables X and Y
Therefore, MIC could not identify the direction of information flow. In this study, a timelag parameter was introduced, and the ability to detect information transmission between two signals was obtained by calculating MIC with different time lags (τ), which was named TDMIC.
where, I_{G} (X, Y, τ) is the MI of the time delay in the case of xbyy (G). When the information of X at time t was decomposed at Y at time t + τ, the JPDF between Y and X had an obvious peak at time t + τ. Naturally, I_{G} (X, Y, τ) was larger than I_{G} (X, Y). Therefore, the sign of the time lag where MIC(X, Y, τ) reached its peak was used to infer the direction of information flow between X and Y.
For the experimental application, to estimate the total flow of information between two physiological time series (EEG and EMG), the cumulative information flow within a certain delay D was estimated using the following equation [6, 26].
In this study, the maximum delay D was set to 40 data points, and the step size k was set to 1 for calculation.
To compare the performance of the algorithm, TDMI and TE methods are also implemented. TDMI method was refered to the code published by Li et al. [21]. TE method with kernel estimator was refered to the code published by Lizier et al. [32]. We also calculated the nonlinear component of TDMIC (NTDMIC) and its cumulative value (C_{NTDMIC}). All algorithms described in this paper were implemented by MATLAB.
Verification of the TDMIC algorithm
In this study, directed linear and nonlinear systems were constructed separately to verify the ability of the proposed algorithm to identify the direction of information flow. Furthermore, the Henon map was used to verify the ability of the algorithm to detect the coupling strength. Considering the randomization of the initial values of X and Y, each model was randomly generated 10 times. Subsequently, the algorithm was applied to the study of FCMC while maintaining ankle dorsiflexion. The data length of both simulation and experimental data was set to 1000 to verify the performance of the proposed algorithm in identifying the information flow between short time series.
Numerical simulation data
Unidirectional dynamical system
First, a unidirectional linear dynamic system was constructed using the following model [33]. The calculations showed that a linear information flow existed from time series Y to X.
Then a unidirectional nonlinear dynamic system was constructed as follows based on the aforementioned unidirectional linear model. A major nonlinear information flow was observed from Y to X.
Bidirectional dynamical system
Second, a bidirectional linear dynamic system was constructed using the following model [21]. That is, a bidirectional linear information flow existed between the time series X and Y generated by the system.
A new bidirectional nonlinear dynamic system was constructed as follows based on the aforementioned bidirectional linear model.
For all models introduced earlier, u_{t} and v_{t} represented two independent and identically distributed (i.i.d) standard Gaussian random variables.
HENON map
Henon map was used to verify the ability of TDMIC to detect the direction and strength of information flow between time series. Two time series (X and Y) with unidirectional coupling relationships were generated using the Henon maps. X and Y were the driving system and the response system, respectively, that is, information flow from X to Y:
where E is the coupling parameter with an interval of [0, 1], and the coupling strength between two time series could be changed by adjusting the value of E.
Experimental data
Ten subjects (mean age, 59.2 ± 7.0 years; range, 50–68 years; 8 male) with chronic stroke (more than 3 months after onset of stroke) and ten healthy controls (mean age, 58.7 ± 7.2 years; range, 46–67 years; 8 male) without any history of neurological disease were recruited. Patient demographics are shown in Table 1. All the subjects were able to complete the experiment as required.
The preparation before the experiment was similar to our previous study [31]. The difference is that we changed the experimental paradigm from the autonomic dynamic dorsiflexion task to a steadystate dorsiflexion task to obtain stable data. There was a crossshaped mark in the center of the computer screen to attract the attention of the participants. After 2 s, a right arrow appeared, prompting the participants to dorsiflexion of the right ankle and maintain this state for 50 s. Figure 1 shows the experimental setup. The participants then rested for 60 s to avoid muscle fatigue. Each participant repeated the aforementioned task five times.
During the study, EEG and EMG signals were simultaneously acquired with an EEG amplifier system (Neuroscan, Australia). Using the international 10–20 system, 26 electrodes were used to record the EEG data (i.e., FP1, FP2, Fz, F3, F4, F7, F8, FC3, FCz, FC4, C3, Cz, C4, CP3, CPz, CP4, P3, Pz, P4, T7, T8, P7, P8, O1, Oz, and O2). The EMG signal from the tibialis anterior (TA) of the right leg was recorded with bipolar electrodes. EEG and EMG data were sampled at 1024 Hz. The electrode wires were fixed with a tape to reduce motion artifacts caused by shaking. The data were maintained for 2–49 s for subsequent analysis to obtain the data under the steady state. Finally, 5 48 slong epochs free of artifacts in each participant were obtained. Data were further cut into 1000 data point segments with no overlapping. 50 Hz power frequency interference was removed, and Bandpass filtering (2–100 Hz) was performed on EEG. Then, the independent component analysis (ICA) algorithm was used to remove artifacts, such as electrooculogram (EOG) and EMG. For EMG, a notch filter was used to remove the 50 Hz power frequency interference, and a 2stage IIR bandpass filter (5–100 Hz) was performed to remove lowfrequency noise.
TDMIC and NTDMIC (the nonlinear component of TDMIC) values were calculated from the beta (14–30 Hz) and gamma (31–45 Hz) bands in the power temporal power map of EEG and EMG obtained by Morlet wavelet transform. Cumulative values of TDMIC and NTDMIC (i.e.C_{TDMIC} and C_{NTDMIC,} represents the total information flow and the nonlinear information flow) were calculated for 40 data point delay. The estimated time delay was close to 40 ms, which was in the range of latencies (20–40 ms) reported between the cortex and muscles [3, 9, 34, 35].
Statistical significance
In this study, the permutation test was used for significance testing. The two original time series were randomly shuffled to generate surrogate data. As for simulated data, the significance level alpha was set to 0.01. For experimental data, the repeatedmeasures analysis of variance (rANOVA, a = 0.05) was performed on TE, TDMI, and TDMIC. Greenhouse–Geisser correction was used to correct the degree of freedom. Bonferroni correction was used for multiple comparisons. All statistical analyses were conducted in SPSS/PC, version 22.0 (SPSS Inc., IL, USA).
Results
Results for numerical models
Figure 2 indicates TDMI and TDMIC values as a function of time lag from two time series generated from the unidirectional models. As shown in Fig. 2, whether it was a linear system or a nonlinear system, both the TDMI and TDMIC curves reached a significant large peak at the positive time lag (linear: τ = 1, nonlinear: τ = 3). The peak values were significantly greater than the significance threshold (a = 0.01). This finding indicated that the direction of information flow recognized by TDMI and TDMIC for unidirectional linear and nonlinear systems was Y to X, which was consistent with the information flow direction of unidirectional models. The TE analysis on the linear system showed that the TE value from X to Y was 3.84 × 10^{–2} and that from Y to X was 2.399 × 10^{–1}. The significance threshold of TE obtained by the permutation test was 9.53 × 10^{–2}. This finding indicated that the TE method recognized unidirectional information flow consistent with the model, that is, from Y to X. On the contrary, for the nonlinear system model, the TE value from X to Y was 3.48 × 10^{–2} and that from Y to X was 1.104 × 10^{–1}. Both the TE values were above the significance threshold. Hence, the TE method was able to identify the causal relationship between X and Y that the model suggested. The results of TE are summarized in Table 2.
Figure 3 indicates TDMI and TDMIC values as a function of time lag from two time series generated from the bidirectional models. Two obvious peaks were observed in TDMIC curves, as shown in Fig. 3(b), which were located at positive and negative lags (linear: τ = ± 1, nonlinear: τ = ± 1). Both peaks were significantly greater than the significance threshold level (a = 0.01). This finding indicated that the bidirectional information flow between X and Y was detected by TDMIC, which was consistent with the models. On the contrary, as shown in Fig. 3a, no obvious peaks above the significance threshold were observed in the TDMI curves. This observation showed that TDMI failed to identify the direction of the information flow of bidirectional systems under this data length (1000). As shown in Table 2, whether it was a bidirectional linear or nonlinear model, both the TE values were below the significance threshold, indicating that the TE method did not recognize a significant information flow between X and Y. This was also inconsistent with the models.
As shown in Fig. 4a, a peak value greater than the significance threshold was observed at the negative lag (τ = –1). This indicated that the direction of information flow recognized by TDMIC for the Henon map was X to Y, which was consistent with the model. In contrast, the information flow direction of the unidirectional Henon map was misinterpreted as bidirectional by TE, as shown in Table 2. Figure 4b shows the ability of TDMIC to detect the coupling strength following the change in the coupling parameter E. The maximum value of TDMIC values also increased monotonically with E. Additionally, a local maximum was observed in Fig. 4b around E = 0.2.
Results for experimental data
Figure 5a and b presents the grand averaged TDMIC and NTDMIC curves as a function of delay in the beta (14–30 Hz) and gamma (31–45 Hz) bands. The EEG signal collected at the Cz position was selected for analysis. Overall, whether it was the TDMIC or the NTDMIC, the betaband information flow from EEG to EMG was stronger than that from EMG to EEG. Interestingly, the ascending information flow (EMG to EEG) in the gamma band was higher than the descending information flow (EEG to EMG). Figure 5c presents the grand averaged TDMI curve as a function of delay in the beta and gamma bands. Compared with TDMIC and NTDMIC, the TDMI curve did not clearly distinguish the ascending and descending information flows. The interval of significance thresholds for TDMIC, NTDMIC and TDMI in the beta band were [0.1546, 0.1554], [0.152, 0.1528] and [0.0030, 0.0035]. In the gamma band, the interval of significance thresholds for TDMIC, NTDMIC and TDMI were [0.1460, 0.1464], [0.1430, 0.1439] and [0.0032, 0.0035], respectively. Both of the TDMIC and NTDMIC values were above the significance thresholds. However, TDMI was below the significance threshold. Additionally, the TE method did not detect significant bidirectional information flow between experimental data. The TE results of the experimental data are listed in Table 3.
Figure 6 shows the grand averaged normalized topographies of C_{TDMIC} and C_{NTDMIC} for controls. The averaged C_{TDMIC} topography of EMG → EEG was similar to that of EEG → EMG with the peak value at similar electrodes: Cz, C3, CP3, P3, Pz and CPz. The difference was that the peak area of the EMG → EEG topographic map was more scattered. In addition, the peak distribution of C_{NTDMIC} in the two directions (EMG → EEG and EEG → EMG) was similar to that of C_{TDMIC}, which was mainly distributed at Cz, C3, CP3, P3, Pz and CPz.
The cumulative values of TDMIC and NTDMIC in the beta and gamma bands for both directions were calculated to further quantify the differences between the controls and the stroke patients. Then, three way rANOVA was performed for each method, with subject (two levels: stroke and healthy control) as a betweensubject factor, with frequency band (two levels: beta and gamma) and direction (two levels: descending and ascending) as withinsubject factors. Figure 7a shows the results of statistical analysis for controls. The results showed that the C_{TDMIC} and C_{NTDMIC} in the beta band were significantly higher in the descending direction than in the ascending direction (C_{TDMIC}: F(1, 18) = 5.07, p = 0.037; C_{NTDMIC}: F(1, 18) = 7.86, p = 0.012). On the contrary, the ascending C_{TDMIC} and C_{NTDMIC} in the gamma band were significantly higher than the descending values (C_{TDMIC}: F(1, 9) = 8.56, p = 0.009; C_{NTDMIC}: F(1, 9) = 11.18, p = 0.004). Figure 7 (b) shows the results of statistical analysis for stroke groups. Different from the control groups, the C_{TDMIC} and C_{NTDMIC} results showed no significant difference in beta or gamma bands in both directions (beta band, C_{TDMIC}: F(1, 18) = 4.20, p = 0.055; beta band, C_{NTDMIC}: F(1, 18) = 0.54, p = 0.473; gamma band, C_{TDMIC}: F(1, 18) = 2.19, p = 0.156; gamma band, C_{NTDMIC}: F(1, 18) = 0.76, p = 0.396).
Furthermore, compared to the controls, the C_{TDMIC} and C_{NTDMIC} results of the stroke groups were significantly weaker both in the beta and gamma bands in the descending direction (i.e., EEG to EMG), as shown in Fig. 8a (C_{TDMIC}, beta band: F(1,18) = 33.0, p = 0.000, Bonferroni; C_{TDMIC}, gamma band: F(1,18) = 27.84, p = 0.000, Bonferroni; C_{NTDMIC}, beta band: F(1,18) = 11.65, p = 0.003, Bonferroni; C_{NTDMIC}, gamma band: F(1,18) = 4.63, p = 0.045, Bonferroni). Similarly, compared to the controls, the results of the stroke groups were significantly weaker both in the beta and gamma bands in the ascending direction (i.e., EMG to EEG), as shown in Fig. 8b (C_{TDMIC}, beta band: F(1,18) = 17.46, p = 0.001, Bonferroni; C_{TDMIC}, gamma band: F(1,18) = 65.68, p = 0.000, Bonferroni; C_{NTDMIC}, beta band: F(1,18) = 4.81, p = 0.042, Bonferroni; C_{NTDMIC}, gamma band: F(1,18) = 10.82, p = 0.004, Bonferroni).
Discussion
This study proposed the TDMIC algorithm to solve the problem of inability to identify causal interactions in MIC applications. The simulation results showed that TDMIC could accurately identify the information flow direction of all models with short data lengths and detect the coupling strength of nonlinear systems. On the contrary, with the same short data length, the performance of TE or TDMI was not as good as that of TDMIC in identifying the direction of information flow. The application of experimental data showed significant bidirectional total and nonlinear information flows in FCMC in the beta and gamma bands. Further analysis showed that the strength of total and nonlinear information flow in the descending direction were significantly higher than that in the ascending direction in the beta band, while an opposite phenomenon was observed in the gamma band. Additionally, strong total and nonlinear information flow mainly acted on the center and contralateral sensorimotor cortex. Further controlled experiments showed that the total and nonlinear information flows in both beta and gamma bands were significantly weaker in stroke group than in healthy control group. This study extended the application of MIC and suggested a new idea for the study of nonlinear coupling components in FCMC.
Compared with TE and TDMI, the TDMIC method could more effectively identify the direction of information flow between short time series, which might be related to the derivation of these algorithms. TE was proposed to explore whether the historical information of the driver could improve the prediction of the state of the recipient [24]. The value of TE (Y to X) between time series X and Y was expressed by the following formula
where \(x_{n}^{k} = \{ x_{n  1} ,x_{n  2} , \ldots x_{n  k} \}\) and \(y_{n}^{l} = \{ y_{n  1} ,y_{n  2} , \ldots y_{n  l} \}\) are k and ldimensional delay vectors, which represent the history of X and Y. The formula showed that TE involved the calculation of highdimensional PDF. This meant that the calculation of TE required long and stable data to accurately construct a highdimensional PDF [21]. In addition, TE was equivalent to GC under Gaussian conditions [36]. Both GC and TE might detect false causality due to incomplete observation of the state of the drive system [11, 37]. In the case of short data length (1000) in this study, the performance of the TE method was not satisfactory, especially for bidirectional linear and nonlinear systems.
As an asymmetric extension of the MI method based on information theory, TDMI also involved the calculation of the PDF. The accuracy of PDF calculations directly affected the validity of TDMI results. Unlike TE, the dimension of the TDMI PDF was only 2, which avoided the problem of highdimensional PDF construction in the TE method. Nevertheless, TDMI still needed long stationary time series data to accurately calculate PDF [6]. Roulston et al. used the standard error formula to prove that MI had obvious errors in the case of short data [38]. In this study, TDMI failed to detect the information flow direction on the bidirectional linear and nonlinear models, thus limiting the application of TDMI to nonstationary EEG signals. The brain proved to be a nonlinear dynamic system [13, 14]. It was difficult to obtain long stationary EEG data in motor task experiments. For instance, in this study, the duration of ankle dorsiflexion was about 1 s. The EEG signal was nonstationary and dynamic during the whole action. Therefore, the application of the TDMI method in shortterm exercise task experiments needed to be carefully evaluated.
Unlike TE and TDMI, TDMIC was an asymmetric extension based on the MIC algorithm. The MIC algorithm ensured that different types of functional relationships were accurately captured by finding the grid division method with the largest MI value [27]. This was different from TDMI in terms of relying on a single PDF calculation method to calculate MI. Especially for complex time series, a single discrete method was not always suitable for different types of functional relations. MIC solved this problem well using the calculation principle. At the same short data length (1000), the performance of TDMIC in identifying the direction of information flow between time series pairs generated by four different models was significantly better than that of TE and TDMI.
The ability to accurately capture the coupling strength between time series is important for evaluating the effectiveness of a new algorithm. The Henon map results showed that the maximum value of TDMIC increased monotonically as coupling strength increase, which was consistent with the trend of the MIC curve. The local maximum we observed in the result was related to the characteristics of the Henon map, namely, it can be interpreted as the minima of the largest subLyapunov exponent [39]. This was also consistent with the previous studies that used Henon map to verify new algorithms [40,41,42]. The value of TDMIC was always greater than the value of MIC. The Henon map had a typical nonlinear unidirectional information flow(X to Y). Therefore, according to the principle of the algorithm, the value of MIC at the time of negative lag was naturally greater than the value at the time lag τ = 0. These results indicated that the TDMIC algorithm could accurately identify the coupling strength between nonlinear dynamic systems.
The direction and strength of the information flow in FCMC needed to be accurately identified to evaluate the motor function and reveal the motion controlfeedback mechanism. Beta and gammaband FCMCs were demonstrated to be associated with movement tasks [8, 43]. The Cz electrode position was considered to be related to leg movement, and therefore EEG signals recorded from the Cz channel were selected for analysis in this study. Significant betaband total information flow was observed in both descending and ascending directions. This observation was consistent with previous findings indicating not only descending motor output information but also ascending somatosensory feedback information [3, 15, 44]. As the cortex and the periphery constituted a closedloop sensorimotor system, the interaction between EEG and EMG was inevitably affected by the bidirectional information flow. Further statistical significance showed that the betaband total information flow in the descending direction was significantly higher than that in the ascending direction, which was consistent with previous findings on the steadystate force output task for the upper limbs [2, 6, 45]. This might be associated with the experimental paradigm of steadystate force output. During the steadystate force output, the task for the upper or lower limbs needed stronger motor control signals than sensory feedback integration. Betaband oscillations affected the transmission of descending control instructions, which were used for force stability and output [1, 46]. Also, a significant bidirectional gammaband total information flow was observed. The difference was that the gammaband information flow in the ascending direction was stronger than that in the opposite direction. This result indicated that the transmission of sensory feedback was the main information flow in the gamma band. The gammaband coupling was confirmed to be related to the generation of dynamic force and the integration of information such as attention, vision, and proprioception [43, 46]. The stronger somatosensory feedback flow in the gamma band might provide evidence for these conclusions.
Additionally, significant bidirectional nonlinear information flow was observed in beta and gamma bands, which might be accounted for by the mechanism of neural signal production. Motor output and somatosensory feedback were mainly produced by nonlinear neuronal interaction in the cortex [47]. Therefore, the bidirectional information flow in FCMC naturally had obvious nonlinear characteristics. Similar to the total information flow, the direction specificity of the nonlinear information flow might also be caused by the experimental paradigm and the different functions of beta and gammaband oscillations. A comprehensive assessment of nonlinear interactions in the sensorimotor system was demonstrated to have clinical significance [48].Our future studies will explore the clinical significance of nonlinear information flow in FCMC.
The coherence between the contralateral sensorimotor cortex and effector muscle (TA) in lower limb tasks has been confirmed by several previous studies [8, 10]. Our study were partially consistent with these previous findings. As shown in the betaband topographic maps, the strong ascending total information mainly flowed to C3, Cz, CP3, P3 and CPz, from where the descending information was output. These electrode positions are generally thought to be associated with the central and contralateral sensorimotor cortex. What differentiates our study from the previous ones is that we observed strong bidirectional total information flows that mainly act on this region. This finding indicated that the central and contralateral positions of the sensorimotor cortex played a major role in motor control and sensory feedback in lower limb motor tasks. A nearinfrared study on gait also confirmed that the medial primary sensorimotor cortices were activated during foot movements [49]. Additionally, the peaks of the topographic map of the total information flow were scattered more in the ascending direction (EMG to EEG) than in the opposite direction (EEG to EMG). This might have to do with the physiological structural difference between the motor control pathway and the sensory feedback pathway. The descending motor output was mainly completed through the corticospinal tract, with direct information transmission. However, the ascending sensory feedback pathway involved the cerebellum, brainstem, and thalamus, with a more complicated information transmission process. The inconsistency in information transmission, which resulted in the positive and negative directions, also showed nonlinear FCMC.
Furthermore, for healthy controls, after separating the nonlinear information flow form FCMC, the bidirectional nonlinear information flow also mainly acted at C3, CP3, P3 and CPz. This was similar to the results of some recent studies on hand tasks [6, 35]. Jin et al. used TDMI to observe a significant nonlinear information flow from the contralateral sensorimotor cortex to the effector muscle during a wrist extension task [6]. However, they did not present a further discussion on the information flow from the effector muscles to the sensorimotor cortex. Recently, Yang et al. used the MSPC method and found the peak of the ascending and descending nonlinear coherences at the CCP3 and C1 electrodes, respectively, during constant contraction of the right upper limb [35]. Our findings indicated that the nonlinear information flow of the contralateral sensorimotor cortex had the dominant role in motor control and sensory feedback regardless of upper or lower limb tasks.
Compared to the healthy controls, we did not observe significant directional differences in the strength of information flow both in the beta and gamma bands. This may be due to structural damage to the patient's brain, which affects normal information interaction [50]. Another possible explanation is that the individual differences among patients caused by factors such as different disease levels, different brain damage locations, and different stroke onset times. Furthermore, compared with the controls, the bidirectional total and nonlinear information flows in both the beta and gamma bands of the stroke group were significantly reduced. This result was consistent with the previous studies [5, 51,52,53]. This weakening of FCMC may be caused by cortical damage or muscle changes resulting from stroke [5, 54]. On the one hand, neural activities through the pyramidal tract were significantly reduced after brain injury, leading to the disassociation of presynaptic and postsynaptic activities, thereby weakening the corticalspinal connection [52]. On the other hand, it has been demonstrated that neuromuscular disorders lead to an increased MU discharge variability and a decreased firing rate after stroke [55]. In particular, Mima et al. previously demonstrated that weak coupling is primarily caused by impaired information flow from the brain to the muscles [56]. The results of weaker descending information flow that we observed in the stroke group was the same. Meanwhile, weaker information flow in the ascending direction was also observed in the stroke group. As mentioned earlier, information flow in the ascending direction plays an important role in somatosensory tasks. Weakening of the ascending information flow may have caused the proprioception disorder of the stroke patient who, in fact, usually suffers proprioceptive dysfunction [57].This study may have demonstrate that the cerebral lesion caused by stroke damages the bidirectional information interaction between the cortex and the effector muscles in the sensorymotor system, and this damage leads to obstacles in limb movement control and proprioceptive feedback.
Conclusions
This study proposed the TDMIC algorithm to address the challenge of accurate identification of information flow in FCMC. Simulation and experimental results showed the effectiveness of the proposed method. This study extended the related research of FCMC on information flow and further explored the frequency specificity and directional specificity of bidirectional nonlinear information flow. The weakening of bidirectional information flow may reflect the underlying mechanism of limb sensorimotor dysfunction after stroke. The proposed method might provide a deeper understanding of the controlfeedback mechanism in motor control and serve a useful tool for the clinical evaluation of motor function. Further studies will recruit more stroke patients for a longterm analysis, focusing on evaluating the effects of different rehabilitation strategies on rehabilitation outcomes.
Availability of data and materials
The datasets of the experiments in the current study are available from the first author on request.
Abbreviations
 FCMC:

Functional corticomuscular coupling
 EEG:

Electroencephalography
 SEMG:

Surface electromyography
 TDMIC:

Timedelayed maximal information coefficient
 TE:

Transfer entropy
 MI:

Mutual information
 GC:

Granger causality
 TDMI:

Timedelay mutual information
 PDF:

Probability density function
 JPDF:

Joint probability density function
 MIC:

Maximal information coefficient
 NTDMIC:

Nonlinear component of TDMIC
 C_{TDMIC} and C_{NTDMIC} :

Cumulative value of TDMIC and NTDMIC
 rANOVA:

Repeatedmeasures analysis of variance
 TA:

Tibialis anterior
 ICA:

Independent component analysis
References
 1.
Baker SN. Oscillatory interactions between sensorimotor cortex and the periphery. Curr Opin Neurobiol. 2007;17(6):649–55.
 2.
Mima T, Matsuoka T, Hallett M. Information flow from the sensorimotor cortex to muscle in humans. Clin Neurophysiol. 2001;112(1):122–6.
 3.
Witham CL, Riddle CN, Baker MR, Baker SN. Contributions of descending and ascending pathways to corticomuscular coherence in humans. J Physiol. 2011;589(15):3789–800.
 4.
Chen X, Zhang Y, Cheng S, Xie P. Transfer spectral entropy and application to functional corticomuscular coupling. IEEE Trans Neural Syst Rehabil Eng. 2019;27(5):1092–102.
 5.
Fang Y, Daly JJ, Sun JY, Hvorat K, Fredrickson E, Pundik S, Sahgal V, Yue GH. Functional corticomuscular connection during reaching is weakened following stroke. Clin Neurophysiol. 2009;120(5):994–1002.
 6.
Jin SH, Lin P, Hallett M. Linear and nonlinear information flow based on timedelayed mutual information method and its application to corticomuscular interaction. Clin Neurophysiol. 2010;121(3):392–401.
 7.
Conway BA, Halliday DM, Farmer SF, Shahani U, Maas P, Weir AI, Rosenberg JR. Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man. J Physiol. 1995;489(3):917–24.
 8.
Gwin JT, Ferris DP. Beta and gammarange human lower limb corticomuscular coherence. Front Hum Neurosci. 2012. https://0doiorg.brum.beds.ac.uk/10.3389/fnhum.2012.00258.
 9.
Xu Y, McClelland VM, Cvetković Z, Mills KR. Corticomuscular coherence with time lag with application to delay estimation. IEEE Trans Biomed Eng. 2016;64(3):588–600.
 10.
Artoni F, Fanciullacci C, Bertolucci F, Panarese A, Makeig S, Micera S, Chisari C. Unidirectional brain to muscle connectivity reveals motor cortex control of leg muscles during stereotyped walking. Neuroimage. 2017;159:403–16.
 11.
He Z, Maekawa K. On spurious Granger causality. Econ Lett. 2001;73(3):307–13.
 12.
Stam CJ. Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. Clin Neurophysiol. 2005;116(10):2266–301.
 13.
Zuo XN, Di Martino A, Kelly C, Shehzad ZE, Gee DG, Klein DF, Castellanos FX, Biswal BB, Milham MP. The oscillating brain: complex and reliable. Neuroimage. 2010;49(2):1432–45.
 14.
Ioannides AA, Mitsis GD. Do we need to consider nonlinear information flow in corticomuscular interaction? Clin Neurophysiol. 2010;121(3):272–3.
 15.
Yang Y, Dewald JPA, van der Helm FCT, Schouten AC. Unveiling neural coupling within the sensorimotor system: directionality and nonlinearity. Eur J Neurosci. 2018;48(7):2407–15.
 16.
Shannon CE. A mathematical theory of communication. Bell Syst Tech J. 1948;27(4):623–56.
 17.
Chen CC, Hsieh JC, Wu YZ, Lee PL, Wu YT. Mutualinformationbased approach for neural connectivity during selfpaced finger lifting task. Hum Brain Mapp. 2008;29(3):265–80.
 18.
Madeleine P, Xie Y, Szeto GPY, Samani A. Effects of chronic neckshoulder pain on normalized mutual information analysis of surface electromyography during functional tasks. Clin Neurophysiol. 2016;127(9):3110–7.
 19.
Jeong J, Gore JC, Peterson BS. Mutual information analysis of the EEG in patients with Alzheimer’s disease. Clin Neurophysiol. 2001;112(5):827–35.
 20.
Vastano JA, Swinney HL. Information transport in spatiotemporal systems. Phys Rev Lett. 1988;60(18):1773–6.
 21.
Li S, Xiao Y, Zhou D, Cai D. Causal inference in nonlinear systems: Granger causality versus timedelayed mutual information. Phys Rev E. 2018;97(5–1):052216.
 22.
Endo W, Santos FP, Simpson D, Maciel CD, Newland PL. Delayed mutual information infers patterns of synaptic connectivity in a proprioceptive neural network. J Comput Neurosci. 2015;38(2):427–38.
 23.
Chai B, Walther D, Beck D, FeiFei L. Exploring functional connectivities of the human brain using multivariate information analysis. Adv Neural Inf Proces Syst. 2009;2009:270–8.
 24.
Schreiber T. Measuring information transfer. Phys Rev Lett. 2000;85(2):461–4.
 25.
Vicente R, Wibral M, Lindner M, Pipa G. Transfer entropy—a modelfree measure of effective connectivity for the neurosciences. J Comput Neurosci. 2010;30(1):45–67.
 26.
Hinrichs H, Noesselt T, Heinze HJ. Directed information flow—a model free measure to analyze causal interactions in event related EEGMEGexperiments. Hum Brain Mapp. 2008;29(2):193–206.
 27.
Reshef DN, Reshef YA, Finucane HK, Grossman SR, Sabeti PC. Detecting novel associations in large data sets. Science. 2011;334(6062):1518–24.
 28.
Zhang Z, Sun S, Yi M, Wu X, Ding Y. MIC as an appropriate method to construct the brain functional network. Biomed Res Int. 2015;2015:1–10.
 29.
Su L, Wang L, Shen H, Feng G, Hu D. Discriminative analysis of nonlinear brain connectivity in schizophrenia: an fMRI Study. Front Hum Neurosci. 2013. https://0doiorg.brum.beds.ac.uk/10.3389/fnhum.2013.00702.
 30.
Bhattacharya J, Pereda E, Ioannou C. Functional associations at global brain level during perception of an auditory illusion by applying maximal information coefficient. Phys A Stat Mech Appl. 2018;491:708–15.
 31.
Liang T, Zhang Q, Liu X, Lou C, Liu X, Wang H. Timefrequency maximal information coefficient method and its application to functional corticomuscular coupling. IEEE Trans Neural Syst Rehabil Eng. 2020;28(11):2515–24.
 32.
Lizier JT. JIDT: an informationtheoretic toolkit for studying the dynamics of complex systems. Front Robot AI. 2014. https://0doiorg.brum.beds.ac.uk/10.3389/frobt.2014.00011.
 33.
Diks C, DeGoede J. A general nonparametric bootstrap test for Granger causality. In: Global analysis of dynamical systems. 2001, pp 391–403.
 34.
Raethjen J, Govindan RB, Binder S, Zeuner KE, Deuschl G, Stolze H. Cortical representation of rhythmic foot movements. Brain Res. 2008;1236:79–84.
 35.
Yang Y, SolisEscalante T, Yao J, van der Helm FC, Dewald JP, Schouten AC. Nonlinear connectivity in the human stretch reflex assessed by crossfrequency phase coupling. Int J Neural Syst. 2016;26(8):1650043.
 36.
Barnett L, Barrett AB, Seth AK. Granger causality and transfer entropy are equivalent for Gaussian variables. Phys Rev Lett. 2009;103(23):238701.
 37.
Smirnov DA. Spurious causalities with transfer entropy. Phys Rev E. 2013;87(4):1–12.
 38.
Roulston MS. Estimating the errors on measured entropy and mutual information. Phys D Nonlinear Phenom. 1999;125(3–4):285–94.
 39.
Schiff SJ, So P, Chang T, Burke RE, Sauer T. Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble. Phys Rev E. 1996;54(6):6708–24.
 40.
Stam CJ, van Dijk BW. Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets. Phys D Stat Mech Appl. 2002;163(3–4):236–51.
 41.
Ma XF, Huang XL, Du SD, Liu HX, Ning XB. Symbolic joint entropy reveals the coupling of various brain regions. Phys A. 2018;490:1087–95.
 42.
Vlachos I, Kugiumtzis D. Nonuniform statespace reconstruction and coupling detection. Phys Rev E Stat Nonlinear Soft Matter Phys. 2010;82(1 Pt 2):016207.
 43.
Omlor W, Patino L, HeppReymond MC, Kristeva R. Gammarange corticomuscular coherence during dynamic force output. Neuroimage. 2007;34(3):1191–8.
 44.
Campfens SF, Schouten AC, van Putten MJ, van der Kooij H. Quantifying connectivity via efferent and afferent pathways in motor control using coherence measures and joint position perturbations. Exp Brain Res. 2013;228(2):141–53.
 45.
Xie P, Cheng S, Zhang Y, Liu Z, Liu H, Chen X, Li X. Direct interaction on specific frequency bands in functional corticomuscular coupling. IEEE Trans Biomed Eng. 2020;67(3):762–72.
 46.
Mehrkanoon S, Breakspear M, Boonstra TW. The reorganization of corticomuscular coherence during a transition between sensorimotor states. Neuroimage. 2014;100:692–702.
 47.
Vlaar MP, SolisEscalante T, Vardy AN, van der Helm FCT, Schouten AC. Quantifying nonlinear contributions to cortical responses evoked by continuous wrist manipulation. IEEE Trans Neural Syst Rehabil Eng. 2017;25(5):481–91.
 48.
He F, Sarrigiannis PG, Billings SA, Wei H, Rowe J, Romanowski C, Hoggard N, Hadjivassilliou M, Rao DG, Grunewald R, et al. Nonlinear interactions in the thalamocortical loop in essential tremor: a modelbased frequency domain analysis. Neuroscience. 2016;324:377–89.
 49.
Miyai I, Tanabe HC, Sase I, Eda H, Oda I, Konishi I, Tsunazawa Y, Suzuki T, Yanagida T, Kubota K. Cortical mapping of gait in humans: a nearinfrared spectroscopic topography study. Neuroimage. 2001;14(5):1186–92.
 50.
Chen X, Xie P, Zhang Y, Chen Y, Yang F, Zhang L, Li X. Multiscale information transfer in functional corticomuscular coupling estimation following stroke: a pilot study. Front Neurol. 2018;9:287.
 51.
Chen X, Xie P, Zhang Y, Chen Y, Cheng S, Zhang L. Abnormal functional corticomuscular coupling after stroke. Neuroimage Clin. 2018;19:147–59.
 52.
Meng F, Tong KY, Chan ST, Wong WW, Lui KH, Tang KW, Gao X, Gao S. Cerebral plasticity after subcortical stroke as revealed by corticomuscular coherence. IEEE Trans Neural Syst Rehabil Eng. 2009;17(3):234–43.
 53.
Mima T, Toma K, Koshy B, Hallett M. Coherence between cortical and muscular activities after subcortical stroke. Stroke. 2001;32(11):2597–601.
 54.
Gardiner R. The pathophysiology and clinical implications of neuromuscular changes following cerebrovascular accident. Aust J Physiother. 1996;42(2):139–47.
 55.
Thomas CK, Butler JE, Zijdewind I. Patterns of pathological firing in human motor units. Adv Exp Med Biol. 2002;508:237–44.
 56.
Mima T, Hallett M. Corticomuscular coherence: a review. J Clin Neurophysiol. 1999;16(6):501–11.
 57.
Hughes CM, Tommasino P, Budhota A, Campolo D. Upper extremity proprioception in healthy aging and stroke populations, and the effects of therapist and robotbased rehabilitation therapies on proprioceptive function. Front Hum Neurosci. 2015;9:120.
Acknowledgements
The authors would like to acknowledge the participation of all subjects in this study. The authors would also like to thank Qian Zang for her contribution to the statistical analysis.
Funding
This research was supported in part by the National Key Research and Development Program of China (2017YFB1401200), Hebei province postdoctoral scientific research project (B2019005001), Key Research and Development Program of Hebei Education Department (ZD2020146) and national natural science foundation of china (61673158).
Author information
Affiliations
Contributions
Tie Liang developed, computational model, performed experiment design,data collection, analysis, and writing of the manuscript. Qingyu Zhang contributed to experiment design, data collection, analysis. Xiaoguang Liu and Bin Dong assisted with data collection, analysis. Xiuling Liu and Hongrui Wang contributed to the funding, experiment design, provision of resources and writing. All authors read and approved the final version of the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
The study was in accordance with the Declaration of Helsinki and approved by the ethics review committee of Affiliated Hospital of Hebei University (HDFYLL2020091). All subjects signed informed and written consent before any procedure of the experiments.
Consent for publication
The authors consent this article for publication of Journal of NeuroEnginnering and Rehabilitation.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Liang, T., Zhang, Q., Liu, X. et al. Identifying bidirectional total and nonlinear information flow in functional corticomuscular coupling during a dorsiflexion task: a pilot study. J NeuroEngineering Rehabil 18, 74 (2021). https://0doiorg.brum.beds.ac.uk/10.1186/s1298402100872w
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s1298402100872w
Keywords
 Functional corticomuscular coupling
 Timedelayed maximal information coefficient
 Information flow
 Nonlinear coupling