Analysis of corticomuscular connectivity during walking using vine copula

1.   Introduction

Walking is a complex task that depends on the complex control of the central nervous system ^[1]. In the neural control of human movement, the cerebral cortex directly or indirectly controls the activity of spinal neurons, causing muscle fibers to contract and produce movement. Studies have shown that walking training can change the state of brain activity and the excitability of the lower extremity cortex and the spinal cord ^[2]. It can activate the direct cortical-spinal pathway through the primary motor area, reflected in the release of cortical potentials and activation of several peripheral motor units ^[3]. This potential synchronous activity reveals the connectivity between the cortex and the muscle. This connectivity indicates that the information flow between neuromuscular entities is related to cortical commands sent to the muscles and to muscle contraction feedback. Therefore, the study of these types of connectivity is conducive to understanding how the brain controls muscles, the impact of muscle activity on brain function, and the potential nature of specific physiological conditions.

The control connectivity between the cerebral cortex and muscles is manifested at different levels of connectivity between electroencephalography (EEG) signals, surface electromyography (sEMG) signals, and EEG–sEMG signals ^[4]. EEG signals are the result among the interactions of neurons in the brain. The sEMG signal is an electrical signal generated by the activity of muscle fiber neurons and reflects the functional response of the muscle to the brain. EEG and sEMG signals can directly reflect the motion intention of the human body and, therefore, reveal the connectivity between the cortex and muscles during movement. Thus, the synchronous analysis of EEG and EMG signals can aid in estimating the functional relationship between the cerebral cortex and muscle tissue and can be used to describe the corticomuscular functional state ^[5,6,7].

Traditional connectivity analysis methods are based on coherence analysis and analyze the EEG and sEMG signals from the time, frequency, time–frequency, and causality domains ^{[8,9,10,11,12]}. Boonstra et al. ^[13] confirmed that neural drives regulate the frequency domain characteristics of sEMG based on the spectral analysis of sEMG captured during different walking tasks. Using the consistency analysis of sEMG and EEG during steady walking in healthy people, Jensen et al. ^[14] proved that the activities of the motor cortex and corticospinal tract during walking directly lead to muscle activity. Most of the research on EEG–sEMG connectivity is based on linear algorithms. In addition, Granger causality (GC) is commonly used in the field of neuroscience to calculate the connectivity between EEG and EMG. GC has been applied to calculate the coupling value between EEG and EMG signals, and it was found that there is bidirectional coupling between the cerebral cortex and muscle. However, interactions between signals generated by neural activity are highly nonlinear and nonstationary, and existing analysis methods have some shortcomings in this regard. They can only describe linear connectivity relationships between signals, and it is not easy to quantitatively analyze complex causalities, such as nonlinearity and higher-order effects. Continuous improvements in nonlinear analysis methods have gained interest, resulting in their increased application in analyzing the connectivity between the cortex and muscles ^[15,16].

To overcome the aforementioned problems, a copula-based method is proposed in this study to perform corticomuscular coupling analysis. Copula theory ^[17], originating from the field of mathematics and statistics, has advantages in analyzing the dependent structure describing multiple variables, which can be applied to the correlation analysis of multiple time series in neuroscience ^[18,19]. Dauwels et al. ^[18] used copula theory to establish a Gaussian diagram model using multichannel EEG signals to infer the interactions between different brain regions. Hu et al. ^[19] used it to study the correlation between mixed data consisting of discrete neural spike sequences and continuous local field potential sequences, which revealed the nonlinear and higher-order causality of coupling between neural signal time series. As problems become more complicated, constructing a high-dimensional copula has become challenging ^[20]. The vine copula function ^[21], based on the copula function and the graphic tool, Vine, was introduced by Bedford to overcome the problem of increasing complexity. The joint distribution is decomposed by the conversion formula between the conditional distribution and the joint distribution, which decomposes multidimensional variables into multiple copula functions and their marginal distribution functions in the form of a vine. High-dimensional problems are more intuitive when using this model than when using traditional copulas ^[22]. In addition, the vine copula model can take many forms, including the canonical vine (C-vine), drawable vine (D-vine), and regular vine (R-vine) ^[23,24,25]. Different copula functions provide different approaches to data analysis. Compared with the C- and D-vine structures, the R-vine structure has no fixed or special structural restrictions, which can better and more flexibly describe the dependencies among multiple variables. Thus, we decided to use the R-vine copula to analyze the coupling between EEG and sEMG. The model decomposes multidimensional variables into multiple copula functions and their marginal distribution functions in the form of a vine; thus, the causality between the brain and muscle multivariate variables can be determined.

2.   Materials and methods

2.1.   Data acquisition

Six healthy subjects (three males and three females; age: 24–26 years; height: 160–180 cm; weight: 48–70 kg) were recruited for the experiments. All subjects read and signed an informed consent form approved by the institutional review board. Simultaneous 64-channel EEG signals were recorded using a portable wireless EEG amplifier, NeuSen. W64 (Neuracle Inc., China), and 6-channel lower-limb sEMG signals were recorded using a Trigno™ Wireless EMG (Delsys Inc., Natick, MA, USA). The sampling rate of the EEG signal was 256 Hz, and all electrode impedances were found to be below 15 kΩ. The sampling frequency of the sEMG was 1000 Hz. All data records were synchronized using labels. Table 1 lists the sEMG electrodes (tibialis anterior (TA), semitendinosus (SEM), vastus medialis (VM), bilateral), and the EEG electrode locations (Fz, F3, F4, Cz, C3, C4, P3, P4, POz) used for the analysis. Subjects were asked to walk forward and backward on a straight section of level ground according to metronome instructions (60 and 120 bmp) and up and down an 8-step staircase in 2 min time blocks, as shown in Table 2. Figure 1 illustrates the experimental setup. The subjects repeated each action 30 times in each group. Each movement sequence (see Table 2) took approximately 6 s to complete. The experiment was repeated 10 times. After every motion type was completed, the participant rested for 5 min to prevent fatigue. During the experiment, the subjects were instructed to fold their hands across their chest and focus their gaze on one point. They were instructed to relax their bodies, refrain from moving their head, and avoid blinking excessively to minimize artifacts in the EEG signals.

2.2.   Vine copula

Copula Theorem: Sklar's theory ^[26] allows F to be a joint distribution function of n-dimensional random variables$X = \left\{{x}_{1}, {x}_{2}, \dots , {x}_{n}\right\}$, with marginal distribution functions ${F}_{i}\left({x}_{i}\right)(i = \mathrm{1, 2}, \dots , n)$. If all marginal functions are continuous, there is a unique copula function C, such that

If the inverse functions of ${F}_{i}(i = \mathrm{1, 2}, \dots , n)$ are ${F}_{i}^{-1}(i = \mathrm{1, 2}, \dots , n)$, then by setting${u}_{i} = {F}_{i}(i = \mathrm{1, 2}, \dots , n)$, the copula function can be calculated as

Similarly, the joint density function $f$ can be defined as

where $c\left({u}_{1}, {u}_{2}, \dots , {u}_{n}\right) = \frac{\vartheta C\left({u}_{1}, {u}_{2}, \dots , {u}_{n}\right)}{\vartheta {u}_{1}{\vartheta u}_{2}\dots {\vartheta u}_{n}}$, which is an n-dimensional copula density function, and ${f}_{i}\left({x}_{i}\right)$ is a marginal density function.

Tree definition: If N is a set of nodes and E is a set of edges, then a tree can be represented by a graph $T = (N, E)$; the tree is constructed using N and E connections and has no cycles.

R-vine definition: $T = \{{T}_{1}, {T}_{2}, \dots , {T}_{n-1}\}$ is a vine on n elements if:

1) ${T}_{1}$ is a tree with nodes ${N}_{1} = \{1, \dots , n\}$ connected by a set of edges denoted ${E}_{1}$.

2) For $i = 2, \dots , n-1$, ${T}_{i}$ is a tree with a node and ${N}_{i} = {E}_{i-1}$ edge set${\mathrm{ }E}_{i}$.

3) For $i = 2, \dots , n-1$ and $\left\{a, b\right\}ϵ{E}_{i}$, in which $a = \left\{{a}_{1}, {a}_{2}\right\}$, $b = \{{b}_{1}, {b}_{2}\}$, and $\#(a, b) = 1$# represents the cardinality of the collection.

Therefore, the one n-dimensional R-vine structure can be expressed using $n-1$ trees (T), and tree ${T}_{j}$ has $n+1-j$ nodes and $n-j$ edges, which refer to the copula $n-j$ density functions. The edges that share one node in tree ${T}_{j}$ become nodes of and are connected by other edges in tree ${T}_{j+1}$. Equation (2.4) defines the R-vine density function.

Here, $x = ({x}_{1}, \dots , {x}_{n})$, $e = \left\{a, b\right\}$ is the edge between nodes a and b in the vine copula tree structure, ${D}_{e}$ refers to the set of variables in the edge ${C}_{e, a}, {C}_{e, a}|{D}_{e}$ represents the binary copula function that characterizes the edge $\left\{a, b\right\}$, and ${f}_{i}$ is the inverse function of ${F}_{i}(i = 1, \dots , n)$.

Copula function selection: There are two model selection criteria: the Akaike information criterion (AIC) and the Bayesian information criterion (BIC).

Here, $L\left(\widehat{\theta }\right)$ is the estimated likelihood, m is the number of adjustable parameters in the copula function, and n is the number of data points.

2.3.   Processing steps used to construct the vine copula

The sEMG and EEG signals were first preprocessed. Comprehensive experiments were performed using MATLAB's (Mathworks Natick, MA, USA) script data analysis feature. The EEG signal was filtered using a 0.5–50 Hz band-pass filter, after which the common average method in MATLAB was used for rereferencing. Independent component analysis (ICA) was used to remove artifacts from the EEG signals and improve signal purity and accuracy. The sampling rate of the sEMG signal was reduced to 256 Hz. The sEMG signal was filtered using a 0.5–50 Hz band-pass filter, and interference was eliminated using the ICA algorithm. The sEMG signal was then denoised using a wavelet threshold. The dependencies among channels were then analyzed for each movement.

The vine copula can be understood as a joint distribution function in which all marginal distributions approximately follow U [0, 1]. The most important action taken was that the preprocessed data were normalized and transformed to follow uniform distributions via probability integral transformation and were used as the input variable of the R-vine copula (called copula data), which are also the source data for establishing the vine copula model. The Kendall-τ correlation coefficient measurement method ^[27], widely used to measure dependence based on ranks, was selected to determine clear dependence structures between channels. Thus, the Kendall–τ correlation coefficient between each channel was obtained. Before choosing the appropriate copula function to describe the correlation between EEG and sEMG, the maximum likelihood estimation (MLE) method was used to estimate the parameters of the binary copula function. A well-fitted copula function for each class having the lowest AIC and BIC criteria was selected. When the two criteria selected different copula functions, the final decision was made by the AIC. The best pair-copula function among common binary copula functions, such as Gauss copula, Student-t copula, Clayton copula, Gumbel copula, or Frank copula, was selected to measure the correlation between two channels in the R-vine copula structure. After obtaining the R-vine structure for six movements, the pair-copula species and the estimated parameters in the R-vine copula model were obtained. The overall process applied the maximum spanning tree algorithm (see ^[28]) to construct the dependency structure of different movements until the number of edges in the tree was less than 2.

3.   Experimental results and discussion

3.1.   Model results

The first layer tree (Tree 1) of the different movements was estimated using the maximum spanning tree algorithm, as shown in Figure 3. The determination of the structure of Tree 1 highlights the direct and indirect dependence of each channel. Figure 3 clearly shows the EEG–EEG correlations. The closer the channel is to the tree structure, the higher its contribution to the movement process and the stronger its correlation with other channels. In the dependent structure, C3, C4, P3, and P4 play a significant role in connectivity. It can be seen that a fixed dependent structure of C4, Cz, C3, P3, and P4 exists under the WBS and WBF states and that there is a fixed dependent structure according to the strong and weak arrangement of correlations in other states. Signals from all participants exhibited similar traits. The performance of sEMG–sEMG also differed under different motion states. From the sEMG channels highlighted with the red markings in Figure 3, the correlations within each leg are significantly better than the correlation between legs, such as VM and SEM, indicating that these results reflect strong linkages among muscles. The correlations between EEG and sEMG can also be described in the walking state. For example, VM, SEM, and C4 were tightly connected under WF and WBS, as highlighted in Figure 3. Significant movement changes in Tree 1 suggest that the cortex and muscles are involved in specific tasks but that differences in their specific connectivity exist. Figure 3 highlights the correlations between EEG–EEG, sEMG–sEMG, and EEG–sEMG.

The best pair-copula function was selected from the commonly used binary copula functions to determine the best dependent vine structure between two channels. After obtaining the R-vine structure for six movements, the pair-copula type and the estimated parameters in the R-vine copula model were obtained. The Kendall–τ correlation coefficient had positive and negative values, indicating a certain dependence between channels and differences in the degree of dependence. In the R-vine copula, the copula function in the first-level tree structure describes the unconditional correlation between the channels. Additionally, the remaining copula functions in the second- to 14^th- level trees describe the conditional correlations. The correlation coefficient reflects the correlation between channels. This method considers the influence of other variables on the relevance of the investigated variables and improves the accuracy of estimating the true connectivity. Although several copula functions exist, only a few are easy to solve for high-dimensional problems. First, the type and number of pair-copula functions were calculated for each channel in the pair-popula parameter estimation for each action, as shown in Table 3. The most appropriate copula function was selected from a common family, among which the Clayton and Student-t functions were particularly popular. The Student-t copula has "fat tail" correlation characteristics. The Clayton copula is sensitive only to one end and, therefore, describes asymmetric correlations. According to the binary function characteristics in the vine copula structure, the copula data in this experiment have both asymmetry and fat tail correlation. The R-vine copula model chooses the most suitable copula function from several binary copula functions to describe nonlinear correlations among data. The various copula functions provide very different results when considering the correlations among data. However, for this same reason, the copula function is more adaptable to different movements.

The determination of the rattan structure not only affects the fitting results of the entire copula model, but also affects certain special dependent information between studied sequences. Therefore, C-vine and D-vine copula models were also established based on EEG–sEMG data to closely compare the fitting effects of different vine copula models. Tables 4–6 list the goodness of the overall fit for each of the three vine copula models. The larger the LogLik (log likelihood) and the smaller the AIC and BIC values are, the lower the error between the constructed model and the actual conditions, the better the model description capability, and the more successful the model construction. Comparatively, the R-vine copula model has a better fit, proving that it is suitable to use R-vine to construct a copula model for analysis of the movements described previously.

3.2.   Connectivity pattern

As shown in Figures 4–6, we used chord diagrams to intuitively visualize the connectivity between the cortex and the muscles by analyzing the Kendall–$\tau$ correlation coefficient of the EEG and sEMG signals. A chord graph is a visualization method that displays relationships between data and mainly describes complex dependencies. Node data are arranged radially along the circumference, and nodes are linked by weighted (wide) arcs. The nodes in the figure represent 15 data channels, the size associated with a node represents the number of channels that depend on the current channel, and the width of the edges represents the connection strength between channels.

The intercortical connectivity is shown in Figure 4, and the full bands highlight the existence of interconnectivity in the cortex for different movements. It can be seen that a speed change has a certain influence when walking on the ground. Compared with WS, the connectivity of Cz, C4, P4, and POz increased during WF, whereas the F3-related connectivity weakened. C3, C4, F3, and F4 played primary connectivity roles in WBS. However, their connectivity decreased with a change in speed in the WBF experiments. Additionally, interconnectivity in the unilateral cerebral cortex was significant during forward walking (WS, WF), whereas the connectivity between the left and right cerebral cortex increased during backward walking (WBS and WBF). Additionally, the intercortical connectivity during WUS and WDS was closer than that when other movements were considered. The interconnectivity in the cortex changes for different motion states, indicating that the activity of the cortex in the walking state varies according to the specific task.

Figure 5 highlights the intermuscular connectivity analyzed in this study. It shows that the intraleg connectivity was higher than the interleg connectivity. No obvious difference was observed for walking forward (WS, WF) at different speeds, except for a weakening of the intraleg connectivity and an enhancement of the interleg connectivity between the VM and SEM muscles. The obvious increase in connectivity observed between the WBS and WBF movements occurred primarily because of the connectivity between the TA and SEM muscles. For the WUS and WDS movements, all muscles exhibited strong ipsilateral muscle connectivity. During any exercise, the correlation with SEM was greater than that for any other muscle. Additionally, the muscle of the left leg was more closely connected than that of the right leg, especially between SEML and other muscles, indicating that the SEM location may be the best electrode position for measuring the correlation between EEG and sEMG in this experiment.

Figure 6 shows the corticomuscular connectivity based on different frequency bands. First, we found that EEG–sEMG connectivity in different EEG frequency bands could be reflected in the gamma band. However, these were more obvious in the theta, alpha, and beta bands. Theta, alpha, and beta exhibited an obvious connectivity during WS and a weakened connectivity during WF, mainly reflected in the alpha rhythm. The strong connectivity between WBS and WBF is mainly reflected in theta, alpha, and gamma. Compared with WBS, we found that the connectivity between the cortical area in the WBF movement and the calf muscle increased in theta but decreased in alpha and that there was no obvious change in the connectivity between beta and gamma. The overall connectivity was lower when walking backward than forward. Clearer corticomuscular connectivity was observed for the WDS movement in each frequency band than for WUS. Among these six movements, WS, WBS, and WDS best revealed the connectivity between the cortex and muscles and the influence of speed change on connectivity. Second, a comparison of the EEG–sEMG connectivity analysis in muscles indicated that the connectivity of the cortex and muscles decreased or increased for different movements. Figure 6 shows that the connectivity strength of the cortex and muscles in the theta and beta rhythms at the VM muscle during WF movement was lower than that during WS, whereas SEM showed the opposite. For all movements, the corticomuscular connectivity on SEML was more prominent than that on SEMR. The SEML muscle is noticeably connected to each cortical area, whereas the connectivity of SEMR varies based on the physical task. Particularly, the cortex and muscles associated with SEMR are weakly connected during WDS, whereas SEMR is closely connected with Fz, F3, and Cz during WUS. In short, owing to the different muscle functions during walking, different forms of coupled oscillations exist between cortical neurons and muscle motor neurons. Different rhythms of the EEG are involved in various functional coupling oscillations in the process of motion control, and different rhythms and muscles lead to varied EEG–sEMG connectivity.

3.3.   Construction of corticomuscular network

We analyzed corticomuscular connectivity defined by 15 nodes (nine EEGs and six sEMGs) based on the links between pairs of nodes. The degree of connectivity for each node was calculated for binary networks to measure the extent to which the cortical nodes influenced the muscle nodes. The correlation of Tree 1 of the vine copula model representing the full bands was applied as the weighted boundary value of the network. A threshold selection strategy ^[29] was used to remove the weaker weighted edges in the adjacency matrix, that is, weaker edges in the adjacency matrix were set to zero. At the same threshold (threshold = 0.04), the resulting adjacency matrix was transformed into a simplified corticomuscular function network, as shown in Figure 7. The red dots indicate the anterior muscles (TA and VM, bilateral), and the yellow dots indicate the posterior muscles (SEM, bilateral). There are significant differences in the connectivity of the network models for different activities. Figure 7 shows that the cortex and muscles are actively involved in movement under specific tasks, as indicated by the connectivity changes. The activity of the corticomuscular network in the walking state changes according to the task. Additionally, the connectivity within one leg was significantly better than that between both legs. As the speed changed, the corticomuscular connectivity during walking became obvious.

In addition, we used the GC method to construct a corticomuscular network that was used for comparison in evaluating the performance of the proposed method. The average weighted clustering coefficient and weighted feature path length of multiple sets of data based on GC and the proposed method were calculated to measure the characteristics of the corticomuscular functional network. A t-test was used to evaluate any significant difference between the average weighted clustering coefficient and weighted feature path length for different movements, and the results are shown in Figure 8. While the size and density of the six movement networks in the corticomuscular functional network were the same, the network characteristics of the movements were different. Different actions produced significant differences. In contrast, the proposed method can better reflect the significant differences between different actions than GC. Speed was also found to have no significant influence on the brain muscle network. In the network constructed by the proposed method, it is particularly obvious that the weighted clustering coefficient of WBF is the largest among all actions, indicating that the network corresponding to the WBF movement is more complex. WBS exhibited the most significant difference to the other movements, and there were no significant differences among the other movements. This may be due to the conscious and highly concentrated movement control of the brain in the backward walking mode and the abnormal muscle stretching and contraction state induced by specific walking movements, leading to significant changes in the connectivity between the muscle and the cortex. It can be argued that the experimental actions in this study can effectively describe corticomuscular connectivity under specific walking conditions.

According to the network characteristics, a network can be generally divided into random, regular, small-world, and scale-free networks. The average weighted clustering coefficient and feature path length of the corticomuscular function network (exp), regular network (ordered), and random network (random) were compared, as shown in Figure 9. The weighted clustering coefficient of the corticomuscular function network was much larger than that of the random network, and there was no significant difference compared to the weighted clustering coefficient of the regular network. Additionally, the average weighted characteristic path length of the corticomuscular function network was longer than that of the random network and shorter than that of the regular network. As the small-world network has the characteristics of both regular and random networks, it has a higher clustering coefficient and a shorter characteristic path. Therefore, this comparison highlights that the constructed network is meaningful and has small-world characteristics.

3.4.   Limitations and future considerations

This study had several shortcomings. First, the walking speed used in this study was self-selected based on a metronome, unlike walking at a constant speed on a treadmill. Each healthy subject completed 10 experiments for each movement, and the number of subjects represents a limited sample size. We plan to include a sample set of subjects who have suffered a stroke in future experiments. Future work should also consider a more rigorous assessment of the role signals play in corticomuscular connectivity at different frequencies. Additionally, the method of corticomuscular analysis used in this study can be improved to describe the motion response state in a specific walking motion and to quantitatively describe the two-way information transmission characteristics of the cortex and muscles in the process. Furthermore, connectivity analysis of the cortex and muscles helps to explore the control mechanisms involved in different movements. Additionally, the two-way information transmission characteristics of the cortex and muscles in the process are quantitatively described during walking and can provide a theoretical basis for walking rehabilitation evaluations. However, more research is required to determine the exact interactions that lead to corticomuscular connectivity and to describe the functional role of the motor cortex in periodic bilateral movements such as walking.

4.   Conclusions

Neurophysiological analysis has gained popularity in the field of neurorehabilitation science. Walking training remains a challenge with respect to lower-limb nerve rehabilitation systems. This study sought to determine the connectivity between signals from the cortex and the muscle during walking using simultaneous EEG and sEMG recordings from healthy subjects asked to perform specific tasks (level ground and stair walking). To model the connectivity problem between EEG and sEMG, the R-vine copula was initially used to study and analyze the related data for network construction. The vine copula used in this study accurately depicted the effective EEG–EEG, sEMG–sEMG, and EEG–sEMG connectivity. The results indicated that connectivity of the cortex and muscles was present during lower limb movement in humans. The method proposed in this paper was used to construct a corticomuscular functional network to verify that the constructed network is meaningful and has small-world characteristics. The analysis demonstrated that the vine copula method is feasible and can effectively describe the connectivity among cortical muscles under specific walking conditions. With the continuous exploration of analysis methods, research on the connectivity of bioelectric signals will become more in-depth and will continue to mature. This should help to further clarify the principle of corticomuscular control in human gait movement and improve its suitability for evaluation and application of rehabilitation treatments.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 61971169, 62071161, 62061044 and 60903084), Zhejiang Provincial Key Research and Development Program of China (No.2021C03031), and Zhejiang Provincial Xinmiao Program of China (Nos.2020R407054, 2020R407055).

Conflict of interest

No potential conflict of interest was reported by the authors.