Time-Frequency-Domain Copula-Based Granger Causality and Application to Corticomuscular Coupling in Stroke

The corticomuscular coupling (CMC) characterization between the motor cortex and muscles during motion control is a valid biomarker of motor system function after stroke, which can improve clinical decision-making. However, traditional CMC analysis is mainly based on the coherence method that can’t determine the coupling direction, whereas Granger Causality (GC) is limited in identifying linear cause–effect relationship. In this paper, a time-frequency domain copula-based GC (copula-GC) method is proposed to assess CMC characteristic. The 32-channel electroencephalogram (EEG) signals over brain scalp and electromyography (EMG) signals from upper limb were recorded during controlling and maintaining steady-state force output for five stroke patients and five healthy controls. Then, the time-frequency copula-GC analysis was applied to evaluate the CMC strength in both directions. Experimental results show that the CMC strength of descending direction is greater than that of ascending direction in the time domain for healthy controls. With the increase of grip strength, the bi-directional CMC strength has an increasing trend. Meanwhile, the bi-directional CMC strength of right hand is larger than that of left hand. In addition, the bi-directional CMC strength of stroke patients is lower than that of healthy controls. In the frequency domain, the strongest CMC is observed at the beta frequency band. Additionally, the CMC strength of descending direction is slightly larger than that of ascending direction in healthy controls, while the CMC strength of descending direction is lower than that of ascending direction in stroke patients. We suggest that the proposed time-frequency domain analysis approach based on copula-GC can effectively detect complex functional coupling between cortical oscillations and muscle activities, and provide a potential quantitative analysis measure for motion control and rehabilitation evaluation.


Introduction
important results. Hu et al. 34 have proposed an e®ective model-free, copula-based GC (copula-GC) measure that can be used to reveal nonlinear and high-order moment causality, and obtained better performance than GC in the cause-e®ect assessment of neural time series, such as local¯eld potential (LFP) and neural spike trains (NST). 16,31,34 Geng et al. 30 have proposed a conditional GC method based on copula by extending the bivariate copula-GC to multivariate time series analysis. It achieved better performance than linear GC, VGC and KGC in motor imagery EEG data analysis. Dauwels et al. 29 and Gao et al. 32 used copula Gaussian graphical models to characterize the interdependence of multi-channel EEG signal oscillation activities and infer the connectivity in various regions of the brain. Ince et al. 33 proposed the statistical framework of Gaussian copula mutual information (GCMI), which can be applied to various types of neural signals and promote comparative study of neural information in order to better understand the information processing function of brain network. From the above research work, the copula method has a solid theoretical foundation and strong application potential, and provides a new research approach for analyzing nonlinear correlation between neural signal variables.
In this study, a time-frequency domain copula-GC method was proposed to investigate the CMC strength between the cerebral cortex and the contralateral muscles for both stroke patients and healthy controls. The copula-GC was applied to EEG and EMG data collected during controlling and maintaining steady-state force output in both time-domain and frequency-domain. Such studies can provide some new characteristics of the CMC strength and°ow after stroke as well as strengthen the understanding of mechanisms underlying motor dysfunction.

Subjects
Five stroke patients (PA) who had motor de¯cits of the upper limb ( Table 1; mean age AE SD: 47.8 AE 2.3 years; range, 45-51 years; all male; all right handed) and¯ve healthy controls (HC; mean age AE SD, 25.8 AE 1.3 years; range, 24-27 years; all male; all right handed) without any history of neurological or psychiatric disease were recruited in this study. None of the subjects were involved in regular strength or endurance training. S1-S5 are healthy volunteers, S6-S8 are patients with mild stroke, S9 and S10 are patients with severe stroke. This study was performed with the oversight of the Institutional Review Board of Guangdong Provincial Work Injury Rehabilitation Hospital. All subjects gave informed consent prior to their participation.

Experiment paradigm
During the experimental, the subjects sat in an electrically shielded and dimly lit room. In order to control and maintain the steady-state force output, a spring grip meter (EH101, Lynx Mall, China) was used in grip tasks at 5 kg, 10 kg and 20 kg force levels. The experimental paradigm is shown in Fig. 1. The task was performed with each hand separately and the hand order was counter-balanced between the subjects. Each steady-state force output trial started with 25 s rest and then the subjects were instructed to perform speci¯c steady-state grip output for 5 s, and subsequently they relaxed for 25 s. Each steady-state grip task included 10 trials with each hand. After each speci¯c steady-state grip task was completed, the subjects rested for 20 min, and then switched to next steady-state grip output task in order to avoid muscle fatigue. As the stroke subjects were weak in exercise performance, the gripping tasks were only performed at 5 kg and 10 kg force levels with each hand. Control subjects successfully accomplished all the steady-state grip tasks.

EEG and EMG data recording
EEG and EMG data were acquired with BrainAmp (Brain Products, Germany) ampli¯ers,¯ltered in the frequency range of 0.5-150 Hz and the sampling frequency was set to 1000 Hz during the acquisition. The EEG signals were referenced to binaural mastoid and recorded using an EEG cap with 32 scalp positions using the international 10-20 standard system. The EMG signals were recorded from the ulnar wrist°exor (UWF),°exor digitorum super¯cialis (FDS), radial wrist°exor (RWF), brachioradialis muscle (BM), musculus biceps brachii (MBB) and triceps ( Fig. 2(b)). Before the electrodes were attached, the hair was needed to clean, and the skin surface was cleaned with alcohol. EEG signals from C3/C4 channels were located in the primary motor zone, and EMG signals from the FDS and UWF channels were selected to calculate the copula-GC values which are further used for analysis in our study.

Copula
In probability theory, copula is a function that links a univariate marginal distribution to a multivariate joint distribution. With copula, the edge distribution can be separated from its joint density distribution, and only the statistical dependency between variables is concerned. Sklar's theorem 35 is the core of copula's statistical theory, pointing out that any multivariate distribution can be expressed as a copula function evaluated on each marginal distribution.
Formally, let X ¼ ðx 1 ; . . . ; x N Þ be a vector random variable with corresponding probability distribution F de¯ned on R N . The copula associated with F is a distribution function C : ½0; 1 N ! ½0; 1 that satis¯es Assuming that F has N th order partial derivatives, its probability density function can be obtained from the distribution function by di®erentiation: where cðuÞ is the copula density function, and p i ðx i Þ is the individual marginal probability density function. The copula function is independent of the marginal distribution and represents the overall dependent structure between the variables.

Copula-based Granger causality
GC is a statistical measure of directional in°uences between two time series. Although increasingly used in various¯elds to detect causal e®ects, GC has not been well adapted to time-varying volatility time series and reveals nonlinear causality. Hu et al. 34 have developed an e®ective copula-based Granger causality (copula-GC, cGC) method for detecting nonlinear and high-order causal relationships without explicit models. Let fX ¼ x t ; Y ¼ y t g be a sample of a stationary stochastic process, and the null hypothesis of causality X ! Y is described as where h refers to the conditional joint density of ðX; Y Þ, f and g are the conditional marginal densities of Y and X, respectively. According to Sklar's theorem, the conditional joint density h can be further written in the form of a conditional copula density function c: where u ¼ Fðy tþ1 jy n t Þ and v ¼ Gðx m t jy n t Þ. F and G represent the condition marginal distributions of Y and X, respectively. Thus, substituting Eq. (5) in Eq. (4), we obtain GC of X ! Y given by It is evident from Eq. (6) that the copula-GC is represented in terms of conditional copula, and it does not involve explicit models.

Time-frequency analysis of EEG-EMG based on copula-GC
The information on the interaction of EEG-EMG signals in the time-frequency domain can reveal the functional relationship between the brain and muscles. In order to quantitatively analyze the interaction between EEG and EMG, the copula-GC value of EEG-EMG is¯rst calculated in time domain. Next, EEG and EMG data were¯ltered into 49 sub-band signals, and then the copula-GC value is computed in each sub-band to perform functional coupling analysis in frequency domain. The detailed steps of the identi¯cation method are described as follows: (i) Data Preprocessing EEG and EMG data were exported to MATLAB 9.3.0 (R2017b, Mathworks Inc, Natick, MA, USA) for subsequent preprocessing and analysis steps. EEG and EMG data are vulnerable bioelectrical signal. They are susceptible to noise in the environment, including power frequency interference and baseline drift. The premise of EEG processing and analytical research is to improve the purity of the signal. The 50 Hz power frequency interference has been¯ltered out in the acquisition software in the acquisition process and a high-pass¯lter was used to remove baseline drift. In this paper, independent component analysis (ICA) was employed to remove the electrooculogram (EOG) signals and then clear the artifacts in the data by using the threshold denoising approach on each layer after wavelet decomposition.

(ii) Time-Domain Analysis
The copula-GC value of EEG-EMG data can re°ect the number of information exchanges between the cerebral cortex and motoneurons. According to Eq. (6), the calculations of copula-GC are performed on the preprocessed EEG and EMG data, where the lag orders of the variables are determined by regression analysis using the Bayesian information criterion (BIC). The copula-GC value of EEG-EMG data can indicate the information exchange between the cortical and spinal cord activities during the execution of a movement.

(iii) Frequency-Domain Analysis
Through the preprocessing of Step 1, the EEG and EMG data were denoted as X and Y , respectively. In order to describe the CMC characterization in speci¯c frequency-domain, X and Y , which ranged from 1 Hz to 50 Hz, are divided into sub-band data with a frequency interval of 1 Hz using a FIR¯lter. In each speci¯c sub-band, the copula-GC values are then calculated based on the de¯nition of Eq. (6). The bilateral copula-GC values between EEG and EMG data in each sub-band at a given frequency f are marked as cGC X!Y ðf Þ and cGC Y !X ðf Þ. According to the de¯nition of copula-GC, the greater the copula-GC value is, the more information is transferred in this sub-band. And to quantitatively analyze CMC characterization in a speci¯c frequency band, a parameter called coupling strength (CS) which refers to the coherence and TE analysis 11,18 was employed to quantitatively calculate the CMC strength in both directions. The CS of descending direction is marked as CS X!Y , which represents the information of the cerebral cortex transit to the motor neuron. Similarly, The CS of ascending direction is de¯ned as CS Y !X , which shows the feedback of the motor neuron to the control of cerebral cortex.
where cGC X!Y ðf Þ and cGC Y !X ðf Þ are the bilateral copula-GC value at the frequency f , respectively, and Áf denotes frequency resolution.

Copula-GC values for PA and HC in time-domain
The copula-GC values in descending and ascending directions were calculated in the time-domain for all healthy controls and stroke patients. Figure 3 showed average bidirectional copula-GC values for all motor tasks across all subjects. Healthy controls showed higher copula-GC values in descending direction than those in the opposite direction. Additionally, we analyzed the di®erence among di®erent gripping tasks at 5 kg, 10 kg and 20 kg force levels and found that the copula-GC values of righthanded (skilled) are higher than those of left hand in both directions. It can also be noticed that the interaction strength of stroke patients had a decreased trend compared to healthy controls for all motor tasks. This may be due to the intracranial hemorrhage leading to the inability of the motor neurons in the motor area to control movements stably, resulting in the obstruction of information transmission between EEG and EMG. In order to further describe the di®erence of CMC strength for PA and HC, according to the calculated copula-GC value, the di®erence between the ascending pathway and descending pathway with respect to each steady-state grip task was assessed by two sample t-test using SPSS 24.0 for windows (IBM SPSS Inc., USA). Results indicated that there was a signi¯cant di®erence (p < 0:05) among PA and HC in both directions for each task. For comparative purpose, linear GC was also applied to the same data. Figure 4 illustrates the results by GC using the MVGC toolbox. 37 The results show that the copula-GC method has achieved better performance than GC in the cause-e®ect assessment of neural time series. As shown in Fig. 4, it can be observed that the overall value of copula-GC is greater than that of GC. This may be due to the reason that the GC method can only re°ect the linear causality between time series, while copula-GC can detect the linear and nonlinear causality. In addition, it revealed no signi¯cant di®erence in either direction between HC and PA subjects using GC. This is similar to the results obtained in applications to cortical LFP time series. 34 These results suggest that the copula-GC method is a sensitive measure which can be used to detect signi¯cant causal in°uence between HC and PA subjects.

Coupling strength values for PA and HC at each frequency band
Di®erent EEG rhythms may be involved in di®erent ways during motion, the oscillatory responses of di®erent frequency bands may di®er for di®erent movements. Therefore, to further demonstrate the di®erences of CMC between PA and HC subjects at theta, alpha, beta and gamma frequency bands, we calculated the copula-GC values in each speci¯c frequency band for all motor tasks across all subjects according to the detailed description in the Sec. 2.2.3. According to Eqs. (7) and (8), the calculated CS values in both two directions are shown in Figs. 5 and 6, respectively. Figure 5 showed the mean CS values in two directions for HC (S1-S5), and Fig. 6 demonstrated the results for stroke patients (S6-S10). Table 2 showed the mean and standard deviation of the CS for all motor tasks across all subjects.  The results show that the CS is signi¯cant in both descending and ascending directions at beta and gamma bands for all grip tasks. Meanwhile, the CS in descending direction was slightly higher than that in ascending direction for HC. In particular, the results also demonstrated that PA had weaker coupling at the beta band in descending direction compared to HC, while the opposite result was observed in ascending direction. Therefore, the CS in descending direction was slightly lower than that in ascending direction for HC. In addition, the CS di®erences between the ascending pathway and descending pathway at the beta band with respect to each steady-state grip task was also assessed using two sample t-test. Results indicated that there was a signi¯cant di®erence (p < 0:05) in the CS among PA and HC in both directions for each task.

Discussion
It is generally believed that post-stroke motor de¯cits arise essentially from the impairments in neural network that controls movement. CMC is believed to be a direct phenomenon of interactions between physiological oscillatory activities from the brain and the controlled muscles. The cerebral cortex controls the movement of muscle tissue through the spinal cord and peripheral nerves, enabling the limb to perform certain motor functions, and the movement of the limbs a®ects the activity of the cerebral cortex through the a®erent nerves. The CMC strength in descending direction indicated the amount of information that the cerebral cortex transmits to the control muscle neurons, and that in ascending direction represented the amount of information that muscle neurons feed back to the cerebral cortex. 13 It can be seen that the CMC is bilateral, 14,18 and it has a closed-loop structure.
In time-domain, the results (Fig. 3) indicated that the coupling between in both directions is more signi¯cant as the grip strength increases. One possible explanation is that the functional coupling of the cortex muscles re°ects the relatively steadystate of control of the cerebral motor cortex. The motor cortex consumes less energy to maintain steady-state motion. As the grip strength increases, more neurons need to be synchronized, resulting in the motor cortex requiring more energy to maintain steady-state motion. This is consistent with the existing research results. 38 The CMC of right-handed (skilled) is greater than that of left-handed, and it may be due to the reason that the brain has stronger control over the dominant hand. In addition, the CMC strength in both directions at each steady-state grip task was also decreased for PA compared to HC. This di®erence may be due to the motor control dysfunction caused by the damage of the cerebral motor function area and thus prevent them activate the motoneuron and motor cortex exactly. 39 The exact underlying mechanism behind the appearance of abnormal coordination patterns during post-stroke recovery has not been conclusively determined, but might be associated with a loss in cortical control and an increased usage of undamaged, indirect descending direct corticospinal pathway to the proximal and distal muscles via the brainstem. 24 In the frequency domain, the results from Figs. 5 and 6 and Table 2 showed that the CMC in both directions was mainly re°ected at the beta frequency band, indicating that the control movement of the upper limb was dominated within the whole beta band, which is consistent with the literatures. 40,41 It also indicated that the coupling oscillation at the beta band can re°ect the maintenance function of the motor cortex for stable motion output. 42,43 Similarly, the CMC di®erences between PA and HC mainly occurred at the beta band, and this was consistent with the literatures. 10,24 However, it was inconsistent with the literature 6 in which the CMC di®erences between PA and HC mainly occurred in the gamma band. This can be explained by the fact that di®erent experiment paradigm was used. The CMC differences demonstrated that the CS at the beta band in descending direction was slightly higher than that in ascending direction for control groups, while the opposite result was observed for stroke patients. The reason is that the CS for PA in descending direction at the beta band decreased, and that in ascending direction was increased compared to HC. One explanation is that intracranial haemorrhage may cause the neurons in the cerebral motion area unable to control the stable motion, leading to more information feedback to the cerebral motor cortex to activate more neurons. 44,45 There were some limitations in this initial study, which was that the conclusion by the proposed time-frequency domain copula-GC was obtained from a limited number of subjects with small size of EEG and EMG data. However, we believed the current discoveries can provide some biomarkers and biological insights that might be useful for the rehabilitative evaluation of motor function impairment after stroke. Further studies will recruit more subjects to further validate the proposed method. Meanwhile, although the copula-GC method has the ability to detecting nonlinear and high-order causal relationships, it cannot separate linear and nonlinear causal components. Recently, Schaeck et al. 46 has proposed a robust time-varying generalized partial directed coherence (rTV-gPDC) method, which can independently obtain the linear and nonlinear causalities, and contrastively analyze the linear and nonlinear relationships in time-frequency domain to study the CMC phenomenon and pathological mechanism of the diseases related with movement.

Conclusion
In this study, a time-frequency domain copula-GC method was proposed to assess the CMC strength of both PA and HC during steady-state grip tasks. Results showed that the CMC phenomenon is bi-directional. In time-domain, the CMC in descending direction is stronger than that in ascending direction for healthy controls. Meanwhile, the CMC increases with the increasing grip strength, and the CMC of righthanded (skilled) is larger than that of the left hand in both directions. It can also be noticed that the CMC for PA tends to be lower than that for HC in both directions. As for the frequency domain, the strongest prede¯ned CS was observed at the beta frequency band during the steady-state grip output. In addition, the CS in descending direction for healthy controls was slightly larger than that in ascending direction. Moreover, PA show a lower CS in descending direction compared to healthy controls, but the CS in ascending direction is increased for PA. Therefore, the CS in descending direction is lower than that in ascending direction for PA. From the preliminary results of this study, the proposed time-frequency domain copula-GC method can e®ectively characterize the coupling characteristics of cortical muscles at di®erent frequency bands and in di®erent directions of delivery, and can provide a theoretical basis for understanding motor control processes and the pathological mechanisms of movement disorders.