# UNDERSTANDING MOTOR UNIT AND MUSCLE ALTERATIONS FOR NEUROLOGIC REHABILITATION

EDITED BY : Ping Zhou, Xu Zhang, William Zev Rymer and Yingchun Zhang PUBLISHED IN : Frontiers in Neurology

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-304-3 DOI 10.3389/978-2-88963-304-3

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# UNDERSTANDING MOTOR UNIT AND MUSCLE ALTERATIONS FOR NEUROLOGIC REHABILITATION

Topic Editors:

Ping Zhou, University of Texas Health Science Center at Houston, United States Xu Zhang, University of Science and Technology of China, China William Zev Rymer, Shirley Ryan AbilityLab, United States Yingchun Zhang, University of Houston, United States

Neurological injuries (such as hemispheric stroke and spinal cord injury, etc.) can result in muscle weakness and spasticity. The affected muscles often experience progressive changes in their intrinsic mechanical properties, giving rise to muscle contracture and associated alterations in muscle internal structural changes. The mechanisms behind the changes are multifactorial, including disuse, autonomic changes, peripheral neuropathy, a loss of central motor neuron trophic influences, and/or their combinations. Regardless of the origins of weakness and other changes, it is very important to understand or quantify complex neuromuscular changes after a neurological injury.

In this Research Topic, we will focus on examination of skeletal muscle changes after neurological injuries because it is the primary organ involved in the generation of force for movement and is also the main effector organ of impairment after neurological injuries, resulting in disability. Of particular interest, the Research Topic will focus on examination of changes in different motor unit components because motor unit is the final common pathway for neuromuscular control and it provides a basic structure-function framework for the examination of neural and muscular disorders. Understanding the changes in different motor unit components (such as motor unit number, size, territory, control properties, etc.) can help identify specific mechanisms causing functional and anatomical changes in the affected muscles, and thereby guiding development of effective treatments for individuals with different neurological injuries.

Citation: Zhou, P., Zhang, X., Rymer, W. Z., Zhang, Y., eds. (2019). Understanding Motor Unit and Muscle Alterations for Neurologic Rehabilitation. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-304-3

# Table of Contents


Zhixian Gao, Lin Chen, Qiliang Xiong, Nong Xiao, Wei Jiang, Yuan Liu, Xiaoying Wu and Wensheng Hou

*33 Inter-Limb Muscle Synergies and Kinematic Analysis of Hands-and-Knees Crawling in Typically Developing Infants and Infants With Developmental Delay*

Qi L. Xiong, Xiao Y. Wu, Jun Yao, Theresa Sukal-Moulton, Nong Xiao, Lin Chen, Xiao L. Zheng, Yuan Liu and Wen S. Hou

*42 Increased Corticomuscular Coherence and Brain Activation Immediately After Short-Term Neuromuscular Electrical Stimulation*

Rui Xu, Yaoyao Wang, Kun Wang, Shufeng Zhang, Chuan He and Dong Ming

*52 Position as Well as Velocity Dependence of Spasticity—Four-Dimensional Characterizations of Catch Angle*

Yi-Ning Wu, Hyung-Soon Park, Jia-Jin Chen, Yupeng Ren, Elliot J. Roth and Li-Qun Zhang

*62 Motor Unit-Driven Identification of Pathological Tremor in Electroencephalograms*

Aleš Holobar, Juan A. Gallego, Jernej Kranjec, Eduardo Rocon, Juan P. Romero, Julián Benito-León, José L. Pons and Vojko Glaser

*77 How Physical Activities Affect Mental Fatigue Based on EEG Energy, Connectivity, and Complexity*

Rui Xu, Chuncui Zhang, Feng He, Xin Zhao, Hongzhi Qi, Peng Zhou, Lixin Zhang and Dong Ming

*90 Muscle Fatigue Post-stroke Elicited From Kilohertz-Frequency Subthreshold Nerve Stimulation*

Yang Zheng, Henry Shin and Xiaogang Hu

*98 Variation of Finger Activation Patterns Post-stroke Through Non-invasive Nerve Stimulation*

Henry Shin, Yang Zheng and Xiaogang Hu

*Ying Chen1,2†, Huijing Hu3†, Chenming Ma1,2, Yinwei Zhan4 , Na Chen2 , Le Li2 \* and Rong Song1 \**

*1Key Laboratory of Sensing Technology and Biomedical Instrument of Guang Dong Province, School of Engineering, Sun Yat-sen University, Guangzhou, China, 2Department of Rehabilitation Medicine, Guangdong Engineering Technology Research Center for Rehabilitation Medicine and Clinical Translation, The First Affiliated Hospital, Sun Yat-sen University, Guangzhou, China, 3Guangdong Work Injury Rehabilitation Center, Guangzhou, China, 4School of Computers, Guangdong University of Technology, Guangzhou, China*

#### *Edited by:*

*Xu Zhang, University of Science and Technology of China, China*

#### *Reviewed by:*

*Bo Yao, Chinese Academy of Medical Sciences and Peking Union Medical College, China Wens Hou, Chongqing University, China*

#### *\*Correspondence:*

*Le Li lile5@mail.sysu.edu.cn; Rong Song songrong@mail.sysu.edu.cn*

*† These authors have contributed equally to this work.*

#### *Specialty section:*

*This article was submitted to Stroke, a section of the journal Frontiers in Neurology*

> *Received: 16 December 2017 Accepted: 22 February 2018 Published: 12 March 2018*

#### *Citation:*

*Chen Y, Hu H, Ma C, Zhan Y, Chen N, Li L and Song R (2018) Stroke-Related Changes in the Complexity of Muscle Activation during Obstacle Crossing Using Fuzzy Approximate Entropy Analysis. Front. Neurol. 9:131. doi: 10.3389/fneur.2018.00131*

This study investigated the complexity of the electromyography (EMG) of lower limb muscles when performing obstacle crossing tasks at different heights in poststroke subjects versus healthy controls. Five poststroke subjects and eight healthy controls were recruited to perform different obstacle crossing tasks at various heights (randomly set at 10, 20, and 30% of the leg's length). EMG signals were recorded from bilateral biceps femoris (BF), rectus femoris (RF), medial gastrocnemius, and tibialis anterior during obstacle crossing task. The fuzzy approximate entropy (fApEn) approach was used to analyze the complexity of the EMG signals. The fApEn values were significantly smaller in the RF of the trailing limb during the swing phase in poststroke subjects than healthy controls (*p* < 0.05), which may be an indication of smaller number and less frequent firing rates of the motor units. However, during the swing phase, there were non-significant increases in the fApEn values of BF and RF in the trailing limb of the stroke group compared with those of healthy controls, resulting in a coping strategy when facing challenging tasks. The fApEn values that increased with height were found in the BF of the leading limb during the stance phase and in the RF of the trailing limb during the swing phase (*p* < 0.05). The reason for this may have been a larger muscle activation associated with the increase in obstacle height. This study demonstrated a suitable and non-invasive method to evaluate muscle function after a stroke.

Keywords: fuzzy approximate entropy, obstacle crossing, stroke, gait, electromyography

### INTRODUCTION

Stroke, a leading cause of disability, often leads to functional limitations in the activity of daily living (ADL). Stroke survivors have a high risk of falling during all poststroke stages (1). Mackintosh et al. found that 36% of community-dwelling elderly people with chronic poststroke symptoms reported falling in the past year, which is significantly more than 24% of the healthy controls (2). Rehabilitation intervention offers beneficial effects on motor recovery after a stroke (3) and can reduce the risk of falling (4). A better understanding of motor function impairment in stroke survivors will help design

**4**

effective recovery strategies during rehabilitation to reduce the incidence of falling.

Obstacle crossing is a complex walking ADL and requires sufficient foot obstacle clearance for the swinging limb, stability of the stance limb (5), and coordination of the whole body to prevent the loss of balance. Half of all stroke subjects either lose their balance or make casual foot contact with the obstacle during crossing (6), which indicates that obstacle crossing threatens the safety of patients after a stroke. Many studies based on the kinematic indices have been conducted to analyze stroke patients' gaits. Kerrigan et al. examined the joint angles of stroke patients during level walking to quantitatively define the most commonly used strategy termed circumduction (7). Lu et al. investigated the motor performance in high-functioning poststroke patients during obstacle crossing and found that stroke survivors appeared to adopt a specific symmetric kinematic strategy with an increased pelvic posterior tilt and swing hip abduction (8). Said et al. quantified the modifications of kinematic characteristics in stroke survivors during obstacle crossing and found that stroke survivors had reduced toe-obstacle clearance and closer horizontal distance after clearance with increased crossing time compared to healthy controls (9).

Previous studies based on the kinematic analysis identified a significant number of stroke-related features for obstacle crossing. Further information about muscle function requires electromyography (EMG) signals, which can be recorded from the muscle surface (10). EMG analysis based on time and frequency domains was widely used in previous studies. Zhai et al. proposed a self-recalibrating classifier of hand movements based on the convolutional neural network using short latency dimension-reduced sEMG as an input (11). Chen and Yang successfully reconstructed drawings of trace reconstructions using a novel three-step hybrid model based on the root mean square (RMS) of seven-channel EMG signals and a gene expression program (12). Kisielsajewicz et al. found that the coherence between synergist muscles in the affected upper limb of stroke patients was lower than that of healthy subjects during reaching tasks (13). Our previous results showed greater muscle activation levels, increased muscle co-contraction, and lower mean power frequencies in persons after a stroke compared to controls during obstacle crossing. These findings indicated that abnormal muscle activation patterns might contribute to difficulties in maintaining balance during obstacle crossing (14). Since the generation of EMG signals is non-linear (15), simple linear modeled features, such as the RMS, integrated EMG, and mean power frequency, reported recently are limited in characterizing muscle dynamics (16). Some non-linear methods have been introduced to analyze the EMG signals, including fractal dimension, average maximum finite-time Lyapunov exponents, and recurrence quantification analysis (15, 17, 18). However, these non-linear dynamic methods usually require very large data sets to achieve reliable results. This may lead to spurious results when applied to small data sets from experiments (19). To solve this problem, entropy-based methods, such as approximate entropy (ApEn), sample entropy (SampEn), and fuzzy approximate entropy (fApEn), have been introduced to analyze EMG signals (19–21). For example, Zhang et al. used SampEn to detect the onset of muscle activity and found that it was more robust than the RMS method (21, 22). It was further used to examine the EMG-torque relation in the complexity domain. This demonstrated that complexity analysis is a novel tool to examine neuromuscular changes after stroke (23).

Entropy was first introduced by Shannon and later termed information entropy (24). Kolmogorov then developed K-S entropy based on the information entropy, which was applicable for examining the complexity of systems (25). However, K-S entropy is not useful for the analysis of measured signals because these signals are noise, and K-S entropy is unable to analyze noisy signal (26). Pincus subsequently introduced ApEn, which is applicable to noisy and small data sets (26). Although ApEn has many advantages compared with linear analysis methods, it is biased. To solve this problem, SampEn was then developed based on ApEn (27). SampEn is less dependent on the size of data sets and shows better relative consistency, but SampEn(*m, r, N*) is not defined in the case of small *N* and *r* (27). Chen et al. later developed fApEn as another complexity analysis method. It combines Zadeh's fuzzy sets with entropy-based methods (19). Due to its excellent robustness and consistency (28), fApEn can analyze muscle function in patients with neuromuscular disorders. Ao et al. found that the fApEn values in the elbow muscles were lower compared to healthy controls (29). Sun et al. found that fApEn values increased with force-generating capacity in stroke survivors during robot-aided rehabilitation training sessions (30). To date, there are limited studies on the dynamics of muscle function during complex tasks, such as obstacle crossing following stroke. This is critical to daily living.

In this study, fApEn was used to analyze the EMG signals recorded from eight muscles of the lower limb of poststroke subjects and compared those with healthy subjects when performing obstacle crossings tasks at different heights. This study aimed to investigate the alterations in the complexity of the EMG signals between the two groups and between different heights during the task. It also aimed to identify dynamic muscle function changes after stroke. Our hypothesis was that the complexity of the generated EMG signals would decrease due to muscle damage after stroke. The complexity would increase along with the obstacle height due to the underlying mechanisms of muscle activation. The correlation between the fApEn values of the EMG signals and the clinical scales could provide further details regarding muscle function after stroke.

### MATERIALS AND METHODS

#### Participants

Five poststroke subjects with at least 3 months onset prior to data collection and were capable of stepping across a 30% leg length height obstacle were recruited. In addition, eight healthy subjects of similar heights and gender participated in the experiment as controls. The Fugl-Meyer Assessment (FMA) and Berg Balance Scale for lower extremities were used to evaluate the motor function of the poststroke subjects. The clinical scales assessments were conducted by an experienced physiotherapist. The basic information for the poststroke patients is shown in **Table 1**. This study was approved by the Ethics Committee of the First Affiliated Hospital of Sun Yat-sen University. This study was conducted in accordance to the Declaration of Helsinki. All subjects provided written informed consent prior to enrollment.

#### Apparatus

The kinematic data were recorded by a 6-camera 3D motion analysis system (Vicon Motion Systems, Oxford, UK). Two force plates (AMTI, Watertown, MA, USA) situated in the middle of the path were used to record the force signals. The heightadjustable obstacle was placed between them. A diagram of the two force paths and obstacles is presented in **Figure 1A**. EMG data were recorded from the rectus femoris (RF), biceps femoris (BF), tibialis anterior (TA), and medial gastrocnemius of both sides for all subjects using preamplified wireless transmission modules. **Figure 1B** shows the flow diagram of the procedure


*To avoid indirectly identifiable patient data, the genders of stroke group were presented as four males and one female; and the ages were presented as a range.*

*BBS, Berg Balance Scale; FMA, Fugl-Meyer assessment scale of the motor function in paretic low-extremity; L, left; R, right.*

for data collection. The apparatuses used in this study were the same as our previous experiment. A detailed description was provided in our previous study (14). The sample frequency for Vicon cameras was 100 Hz and 1 kHz for the force plates and EMG modules.

#### Procedure

The subjects' heights, body weights, and leg lengths were first measured and recorded before kinematic data collection. The distance from the anterior superior iliac spine to the lateral malleolus was measured as leg length. This was then used to calibrate the height of the obstacle of each individual. Thirtyfive 15-mm light-reflective markers and silver-silver chloride (Ag-AgCl) electrodes were attached to corresponding positions on each of the subjects. The target skin area was shaved and cleaned with alcohol to obtain better signals before the attachment of the electrodes (14).

The gait trials began after the preparation. The subjects were asked to walk along a volunteered walkway (8 m) at a volunteered speed with bare feet with an obstacle placed at a midway distance. Details of the trials were described previously (14). **Figure 1C** presents the gait cycle during obstacle crossing. Trials in which the subjects touched the obstacles or asked for assistance were ignored, and three successful trials for each height were recorded.

All subjects completed the maximum voluntary contraction tasks and three different height obstacle crossing tasks. No incident of fall was observed during all trials. The trials where help was received from therapist or the obstacle was touched were discarded. Discomfort or feelings of fatigue were not reported by any subjects during the tasks.

#### Chen et al. Stroke-Related Changes during Obstacle Crossing

#### Data Processing

A 20-Hz low-pass fourth-order Butterworth filter was employed to filter the kinematic and kinetic data. When the toe marker was 2 mm off the ground, this was regarded as the toe-off time. The heel strike time could be recognized according to the change of force signal received by the force platforms. The gait cycle was then divided into two phases or a single lower limb: swing phase and stance phase. The raw EMG signals were collected at a frequency of 1 kHz and then were filtered through a fourth-order Butterworth filter with a frequency band from 10 to 350 Hz. As the frequency of mains of power supply was 50 Hz, a digital notch filter was used to subtract the disturbance of strong electromagnetic fields of 50 Hz that were present in the experiment conducted area.

The fApEn of an *N* sample series is computed as follows:

First, for a given *m*, we formed *m*-dimensional vector sequences:

$$\begin{aligned} X\_i''' &= \left\{ \mu\left(i\right), \dots, \mu\left(i+m-1\right) \right\} - 1 / m \sum\_{j=0}^{m-1} \mu\left(i+j\right), \\ \text{where } \mu\_0\left(i\right) &= \frac{1}{m} \sum\_{j=0}^{m-1} \mu\left(i+j\right). \end{aligned}$$

For every *Xi <sup>m</sup>* , the distance *dij <sup>m</sup>* between the two vectors *Xi m* and *Xj <sup>m</sup>* is:

$$d\_{\vec{v}}^{m} = \max\_{k \in \{0, m-1\}} \left| \mu \left( i + k \right) - \mu\_{\boldsymbol{0}} \left( j \right) - \mu \left( j + k \right) - \frac{1}{m} \sum\_{j=0}^{m-1} \mu \left( i + j \right) \right| \quad (1)$$

The definition of similarity degree between two vectors *Xi m* and *Xj <sup>m</sup>* is as follows:

$$D\_{ij} = \exp\left(-\left(\frac{d\_{ij}^{\
u n}}{r}\right)^{n}\right),\tag{2}$$

where *r* and *n* in Eq. 2 determine the width and gradient of the boundary.

The function φ*<sup>m</sup>* averages the similarity and is defined as follows:

$$\left[\phi^{\mathfrak{m}}\right](n,r) = \frac{1}{N-m+1} \sum\_{i=1}^{N-m+1} \ln\left(\frac{1}{N-m+1} \sum\_{i=1}^{N-m+1} D\_{\vec{\eta}}^{\mathfrak{m}}\right). \tag{3}$$

The fApEn is then calculated as follows:

$$\text{fApEn}\left(m, n, r, N\right) = \text{Inq}^m\left(n, r\right) - \text{�wp}^{m+1}\left(n, r\right). \tag{4}$$

Here, *m* = 2 and *r* = 0.15 \* SD(signal) were set according to the previous study (30).

#### Signal Processing and Statistical Analysis

The fApEn values for all the muscles were averaged over three replicates for each subject during each height obstacle crossing. The SD values were also calculated. A two-way (group: control and poststroke × obstacle height: 10, 20, and 30% of leg length) repeated measure of variance (ANOVA) was performed on the fApEn values. A Bonferroni *post hoc* test was used to analyze the fApEn values. Kolmogorov–Smirnov test was applied to the variables. Pearson's correlation coefficient was used to examine the relationship between the clinical scales and fApEn values when the variables were normally distributed. Spearman's correlation coefficient was used when the variables were non-normally distributed. The significance level was set at 0.05. All the data were analyzed in SPSS 19.0 statistical software (SPSS Inc., USA).

#### RESULTS

**Figure 2** shows that the fApEn values of the four lower limb muscles of the trailing limb during the stance phase. **Figure 3** presents the fApEn values of these four muscles during the swing phase. As presented, during this gait cycle, the fApEn values of poststroke subjects were lower than those of healthy controls. In addition, significantly lower fApEn values were found in the RF of poststroke subjects during the swing phase when compared with healthy controls (*p* < 0.05). As shown in **Figures 2** and **3**, most fApEn values of the four lower limb muscles of the trailing limb increased with the height of the obstacle. Furthermore, a significant increase was observed in the BF during the swing phase when the height of the obstacle increased from 10 to 30% of leg length (*p* < 0.05).

**Figures 4** and **5** present the results of the four lower limb muscles of the leading limb during the swing phase and stance phase. As shown in **Figure 4**, during the swing phase, the fApEn values of the BF and RF in poststroke subjects were higher than in the healthy controls, and this result was similar for the TA during the swing phase when the obstacle height was 20 and 30% of leg length. Meanwhile, during the swing phase, fApEn values for all four muscles were lower in poststroke subjects compared with healthy controls. However, these differences between groups were non-significant (*p* > 0.05). Similar to the results of the trailing limb, the increase in the fApEn of the muscles with the obstacle height was also found in the leading limb. In addition, as presented in **Figure 5**, the fApEn value of the BF during the stance phase was statistically significantly greater when the obstacle height was 30% of leg length compared with 10 and 20% (*p* < 0.05). The correlations between fApEn and the two clinical scales were not statistically significant.

#### DISCUSSION

We recorded EMG signals and calculated the fApEn values of four lower limb muscles of poststroke subjects during different phases of obstacle crossing at different heights. Complexity change in muscle activations were then compared between poststroke subjects and healthy controls when they conducted this challenging task.

#### fApEn Values Change after a Stroke

The decreased fApEn values of EMG signals in poststroke subjects could be explained by that muscles were damaged,

Figure 3 | The details of fuzzy approximate entropy (fApEn) values of each height for trailing limb during swing phases. (A) The fApEn values of rectus femoris (RF); (B) the fApEn values of biceps femoris (BF). (C) The fApEn values of tibialis anterior (TA); (D) the fApEn values of medial gastrocnemius (MG). \*Significant effect between groups. The bar (-) indicates significant effect between heights.

and they became disused after stroke. This might lead to their degenerated internal structure (31). The motor unit properties changed because of the reduced corticofugal output from the paretic hemisphere (32). In addition, the number of functioning motor units and the firing rate decreased with reduced discharge variability after a stroke (33, 34). These changes directly

Figure 4 | The details of fuzzy approximate entropy (fApEn) values of each height for leading limb during swing phases. (A) The fApEn values of rectus femoris (RF); (B) the fApEn values of biceps femoris (BF). (C) The fApEn values of tibialis anterior (TA); (D) the fApEn values of medial gastrocnemius (MG).

Figure 5 | The details of fuzzy approximate entropy (fApEn) values of each height for leading limb during stance phases. (A) The fApEn values of rectus femoris (RF); (B) the fApEn values of biceps femoris (BF). (C) The fApEn values of tibialis anterior (TA); (D) the fApEn values of medial gastrocnemius (MG). The bar (-) indicates significant effect between heights.

affect the electrical activity and might lead to reduced muscle force and disability in subtle responses to the perturbations during functional tasks. Here, the reduction in fApEn values was reflected in the decreased complexity of the EMG signals of poststroke subjects, which might be related to the alternations in the properties of motor units. Our findings were consistent with those reported by Ao et al., who found lower fApEn values in elbow muscles of poststroke subjects compared to healthy controls during trajectory-tracking tasks. This could be attributed to a reduced number and firing rate of active motor units (29). Similarly, the decreased complexity of EMG signals was also reported in other patients with neuromuscular disorders, such as Parkinson's disease (35) and cerebral palsy (36), which could be explained by the disease-induced muscle fiber degeneration (37).

During swing phase, the fApEn values for the BF, RF, and TA (20 and 30% leg length) in the leading limb were higher in poststroke subjects than in healthy controls. This might be as a result of abnormal gait during obstacle crossing after a stroke. The swing phase of the leading limb during obstacle crossing caused the subjects to elevate their lower limb to secure sufficient toeobstacle clearance, and this could be a challenge for poststroke subjects, leading to abnormal muscle activation patterns. Indeed, our previous study found that to avoid falling, toe-obstacle clearance of stroke survivors was greater than in healthy controls (38). The increased activation of thigh muscles in the BF and RF was found in stroke survivors (14). This might also contribute to the increased complexity of the EMG signals (36).

### fApEn Values Change with Task Difficulties

When sustaining different levels of maximal voluntary contraction force in the upper limb, the complexity of EMG signals had been demonstrated to increase with increasing muscle contraction forces (23, 37, 39). In line with these findings, our results demonstrated that increasing obstacle heights demanded an increase in muscle contraction forces that in turn led to the recruitment of more motor units and increased firing rates in active motor units (14). Our results suggested that the complexity of the EMG signals increased with greater task demands, and this could also be applied to poststroke subjects. Therefore, safely crossing higher height obstacles requires increased muscle contraction forces and more activated motor units, leading to higher entropy values for the EMG signals (37).

However, there are still discrepancies about the complexity changes with task difficulties. To investigate the effect of task demands on motor entropy, Hong and Newell found that the entropy values of muscle forces decreased as the task demands increased (40). They explained the decreased entropy with increased task demands, but reduced environmental information, revealing a compensatory interaction between tasks and the environment on the force dynamics. Moreover, Murillo et al. found that fuzzy entropy of postural sway in healthy young adults decreased from the stable condition to the mid-level instability condition. This increased again at the highest instability condition at the anterior–posterior axis, which reflects the adaptations of postural control system to the platform instability (41). Therefore, the compensatory and adaptive nature of the motor control system to the task complexity warrants further investigation, especially in stroke survivors. Entropy analysis could be used to evaluate the effects of rehabilitation interventions targeting the motor recovery to restore complex motor tasks in persons after a stroke.

#### Limitations

There were several limitations in this study. First, considering the insufficient strength of the paretic leg during the stance phase, we did not instruct the poststroke subjects to first step over the obstacle with their unaffected side due to safety issue. Thus, we could not compare the paretic side with the unaffected side during the same task. In the future, we should introduce stroke subjects to first step over the obstacle with both affected and unaffected limbs. Second, moderate to high functional level of persons after stroke were recruited in this study. A large-scale study of different types of stroke subjects should be recruited in future study to investigate the influences of group and obstacle height, which may help explore the mechanisms and guide rehabilitation after stroke.

### CONCLUSION

In this study, the stoke-related changes in complexity of lower muscles during obstacle crossing were investigated using fApEn. Results show that the complexity of RF in trailing limb during stance phase decreased in stroke group, which might be associated with the reduced number and firing rate of MU. However, during the swing phase, there were non-significant increases in the fApEn values of BF and RF in the trailing limb of the stroke group, resulting in a coping strategy when facing challenging tasks. During the gait, the complexity of muscle activation increases with obstacle height. That might be because higher obstacles demand greater muscle forces, which causes more motor units to be recruited and triggers higher firing rates of motor units. These findings based on the fApEn values of the EMG signals indicate that the complexity analysis using fApEn could be a suitable and non-invasive method to evaluate muscle function changes after stroke.

### ETHICS STATEMENT

The study was approved by the Ethics Committee of the First Affiliated Hospital of Sun Yat-sen University. The study was conducted in accordance to the Declaration of Helsinki. All subjects provided written informed consent prior to enrollment.

## AUTHOR CONTRIBUTIONS

YC, CM, and NC conceived and designed the study and performed the experiments. YC and HH wrote the paper. RS and YZ contributed to experiments. RS and LL reviewed and edited the manuscript. All authors had read and approved the manuscript.

### ACKNOWLEDGMENTS

The authors would like to thank all the participants of this study.

### FUNDING

This work was supported in part by the National Natural Science Foundation of China (No. 31771016, 81702227), the Guangdong Science and Technology Department (No. 2014B090901056, 2015B020214003, 2017B020210011), and the Project of Science and Technology Program of Guangzhou (No. 201604016034, 201604020108).

### REFERENCES


array. *IEEE J Biomed Health Inform* (2016) 21(6):1562–72. doi:10.1109/ JBHI.2016.2626399


**Conflict of Interest Statement:** All financial, commercial, or other relationships that might be perceived by the academic community as representing a potential conflict of interest must be disclosed. If no such relationship exists, authors will be asked to confirm the following statement.

*Copyright © 2018 Chen, Hu, Ma, Zhan, Chen, Li and Song. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Ultrasonography Monitoring of Trauma-Induced Heterotopic Ossification: Guidance for Rehabilitation Procedures

Qing Wang1,3 \*, Peizhen Zhang<sup>1</sup> , Pengdong Li <sup>2</sup> , Xiangfen Song1,3, Huijing Hu<sup>2</sup> , Xuan Li <sup>4</sup> , Wufan Chen1,3 \* and Xiaoyun Wang<sup>2</sup> \*

*<sup>1</sup> Department of Biomedical Engineering, Southern Medical University, Guangzhou, China, <sup>2</sup> Guangdong Work Injury Rehabilitation Center, Guangzhou, China, <sup>3</sup> Guangdong Provincial Key Laboratory of Medical Image Processing, Southern Medical University, Guangzhou, China, <sup>4</sup> Nanfang Hospital, Southern Medical University, Guangzhou, China*

#### Edited by:

*Xu Zhang, University of Science and Technology of China, China*

#### Reviewed by:

*Ke-Vin Chang, National Taiwan University Hospital, Taiwan Wens Hou, Chongqing University, China*

#### \*Correspondence:

*Qing Wang wq8740@smu.edu.cn Wufan Chen chenwf@smu.edu.cn Xiaoyun Wang xiaoyunwang77@hotmail.com*

#### Specialty section:

*This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology*

Received: *30 June 2018* Accepted: *24 August 2018* Published: *13 September 2018*

#### Citation:

*Wang Q, Zhang P, Li P, Song X, Hu H, Li X, Chen W and Wang X (2018) Ultrasonography Monitoring of Trauma-Induced Heterotopic Ossification: Guidance for Rehabilitation Procedures. Front. Neurol. 9:771. doi: 10.3389/fneur.2018.00771* Traumatic injury is one of varying causes of heterotopic ossification (HO). After HO occurrence, rehabilitation training need alterations to avoid the aggravation of HO. Therefore, monitoring of HO development plays an important role in the rehabilitation procedure. The aims of this study are to evaluate the post-traumatic HO occurring at various joints, to describe the features of HO development in ultrasound images, and to provide a guidance for the orthopedist to make individualized rehabilitation therapy. Eight subjects with the post-traumatic HO were recruited in this study. The joints on the injured side was examined by plain radiography. The joints on the injured side and the corresponding sites on the uninjured sides were scanned by ultrsonography. The HO tissues were segmented automatically using a semi-supervised segmentation algorithm. Then the HO tissues were evaluated in comparison with the corresponding region of the uninjured side. During the development stage of immature HO, ultrasonography was sensitive to observe the involved soft tissue and the calcification of HO. The characteristics of HO tissues in ultrasound image included the hyperechoic mass occasionally accompanied with acoustic shadow and the irregular muscular architecture. It was found that the mean grayscale value of HO was significantly higher (*p* < 0.001) than that of the uninjured side at the middle and late stages. During the development period of HO, the HO grayscale value gradually increased and the mean grayscale of value of mature HO was significantly higher (*p* < 0.05) than that of immature HO. According to the information of HO provided by ultrasound, the orthopedist properly adjusted the rehabilitation treatment. The results demonstrated that the visualization of HO using ultrasonography revealed the development of HO in the muscle tissues around the injured joints and thus provide a guidance for the orthopedist to make individualized rehabilitation therapy. Ultrasound could be a useful imaging modality for quantitative evaluation of HO during the rehabilitation of traumatic injury.

Keywords: heterotopic ossification, ultrasonography, trauma, rehabilitation, diagnosis

### INTRODUCTION

Heterotopic ossification (HO) has been clinically described as lamellar bone formation in the periarticular soft tissues, where osseous tissue should not exist. This aberrant bone formation is commonly associated with orthopedic interventions, trauma, stroke, traumatic brain injury, spinal cord injury, and neurological disorder. Therefore, it is summarized under three ways of acquiring HO: traumatic, non-traumatic (very rare), or neurogenic. However, the definitive pathogenesis of HO is still quite unclear. Early HO causes restriction of joint movement, local pain and swelling. The clinical symptoms of early superficial HO may include erythema and localized warmth (1, 2). At the late stage of HO, the mature osseous tissues lead to severe limitation of the range of motion of the joint, pain in the affected joints, and even nerve or vessel compressed by HO (2).

In the clinical laboratory, biochemical markers such as serum alkaline phosphatase (AP) and bone alkaline phosphatase (BAP) are usually examined to reveal the alterations in bone metabolism. The changes in biochemical markers reflect the bone formation or bone resorption, but could not provide visualization of the change of bone (3). Citak et al. proved that laboratory findings of elevated AP and BAP might not be reliable for early HO detection because of their low diagnostic specificity (4). Therefore, the laboratory examinations are not suggested to be routine examination.

Nowadays, several imaging modalities have been applied to evaluate HO. In department of orthopedics in a hospital, plain radiography is usually used to locate and visualize the heterotopic bone tissues. Plain radiograph is considered as the gold standard of the clinical diagnosis of HO because it is inexpensive and convenient to detect HO. Computed tomography (CT) is sometimes used to provide 3-D information of HO and clearly observe the location and volume of HO. However, the cost of CT exam is high. Patients are exposed to a higher dose of radiation by plain radiography and CT, which only provide the diagnosis of HO in the late stage that the osseous tissues develop into the matured bone. Therefore, these two imaging modalities are unable to detect ossification in the early inflammatory stage. The HO progression needs to be further evaluated (2). Magnetic resonance imaging (MRI) is useful for imaging soft tissue in the early stage of HO. Ledermann et al. proposed a HO grading scale based on MRI characteristics to evaluate the development of HO: grade 1 = fluid attenuation without calcifications, grade 2 = calcification of soft tissue without evidence of bone formation, grade 3 = immature bone formation, and grade 4 = mature bone with cortical formation (5). However, MRI has low resolution and is expensive and relatively insensitive to the bone tissues.

The advantages of ultrasonography are that ultrasound avoids ionizing radiation, and is widely available, repeatable and inexpensive for bedside monitoring (6). It is proved that ultrasound is a sensitive imaging method in evaluation of soft tissue lesions and calcifications (7). Some previous studies utilized ultrasonography to evaluate the muscle tissues by measurements of grayscale value and thickness of the muscles (8, 9). Furthermore, the relationship between ultrasound measurements and the diseases were explored (8, 10). A recent study revealed the negative correlation between the echo intensity of the rectus femoris and muscle strength (11). Other previous studies demonstrated that ultrasound could be effective for detecting immature HO. Bedside ultrasound were applied to diagnose immature HO caused by brain or spinal cord injury and the results suggested that ultrasonography could be a useful first-line imaging modality in the diagnosis of early HO (12–14). Recent studies stated that ultrasonography could distinguish matured HO from the surrounding soft tissues with a high specificity (15–17). Trauma leading to artificial joint replacement and internal fixation of bone fracture without neurological injury, may also cause HO. The rehabilitation outcome of the injured joints may be greatly affected by HO. Therefore, ultrasonography could be applied in orthopedic rehabilitation for providing a guidance for the orthopedist to make individualized rehabilitation therapy. However, quantitative evaluation of HO and serial follow-up ultrasonography for depicting HO progression were rarely reported.

In this study, ultrasonography was applied to observe the calcification in soft tissues in the participants with the posttraumatic HO caused by trauma such as bone fracture and contusion in the rehabilitation center. The aims of this study are to evaluate the post-traumatic HO occurring at various joints, to describe the features of HO development in ultrasound images, and to provide a guidance for the orthopedist to make individualized rehabilitation therapy. The ultrasonographic results were compared with those of plain radiography, which is the gold standard of the clinical diagnosis of HO in orthopedics.

#### MATERIALS AND METHODS

#### Participants

The participants were selected from Guangdong Work Injury Rehabilitation Center, China. The criteria of participants include the following: (I) trauma such as fractures, dislocations and contusion occurring at the joints, (II) ossification in soft tissues surrounding the injured joint but not attaching to the cortical bone visualized by radiograph, and (III) motion restriction and pain of the affected joints. Eight participants (6 males and 2 females, Age: 23–48) with symptoms of HO were recruited. The exclusion criteria include the following: (I) inadequate acquisition of ultrasound images, and (II) osteochondroma. The basic information for the participants is shown in **Table 1**. Four participants had post-traumatic HO at the elbow joint (50%), 3 at knee joint (37.5%), and 1 at shoulder joint (12.5%).

The clinical assessments of the range of joint motion were conducted by an experienced physiotherapist. The study was approved by the Ethics Committee of the Guangdong Work Injury Rehabilitation Center, Guangzhou, China. All subjects provided written informed consent prior to enrollment.

#### Apparatus

In this study, four participants No. 1, 2, 3, and 4 were scanned by a portable ultrasound system (Mindray M5, Mindray, Shenzhen, China) with a linear transducer (Mindray 7L4s, a range of central


TABLE 1 | Background data of the participants.

frequency from 5 to 10 MHz). The other four participants No. 5, 6, 7, and 8 were scanned by a wireless ultrasound probe with a fixed central frequency of 10 MHz (Uprobe-3N, linear transducer, SonoStar, Guangzhou, China). A pilot test was conducted to ensure that the B-mode grayscale images obtained by the two ultrasound systems were of similar image quality.

All participants underwent examination of plain radiography (Aristos VX plus, SIEMENS, Germany). Two participants No. 4 and 5 were also examined by CT to obtain more information of the injured joint to distinguish HO from osteochondroma.

### Procedure of Ultrasound Examination

After examination of plain radiography, all participants were examined by ultrasonography. First, the injured joint was scanned at the site of HO occurrence (**Figure 1**). The ultrasound probe was placed on the scan site with enough coupling gel interposed between the transducer and the skin. The depth of view of the ultrasonographic scan was set at 4 cm for the joints to clearly display the HO lesions and the structure of the muscles and bone. The gain was fixed at a constant intensity for all scans. All the ultrasound images were obtained respectively at the short and long axis of the muscles. Next, the same scanning procedure was performed at the corresponding site on the uninjured side as control.

Because participant No. 8 was with HO on the immature bone formation stage, he underwent the first ultrasound scan 2 days after the first plain radiography and the follow-up ultrasound scan was performed every week. Four follow-up scans were performed to monitor the development of HO. At each time point, five longitudinal and transverse ultrasound images of HO were obtained at the adjacent scan sites with an interval of approximately 0.5 cm, respectively. During this period of ultrasound monitoring, the participant was also underwent rehabilitation therapy. At the end of this study, the participant underwent an additional X-ray examination again.

FIGURE 1 | Longitudinal ultrasound scanning of HO at the posteriorlateral site of the injured elbow joint of participant No. 1.

For the participants with mature HO, the ultrasound scan was performed 2 days after plain radiography without the follow-up ultrasound scan.

### Segmentation and Assessment of HO

Because the edge of the target tissue was not easily distinguished from the background in the ultrasound image, manual segmentation of the HO tissue was time consuming and operator dependent. Therefore, this study applied a semi-supervised segmentation algorithm based on patch representation and continuous min cut (18) to semi-automatically segment the region of interest (ROI) of HO in ultrasound image. Under semisupervision of the clinical expertise, the HO can be accurately and specifically segmented from the surrounding soft tissues according to the texture features of the HO and muscles. **Figure 2A** shows the HO ROI segmented from the background in ultrasound grayscale images of the injured joint. Meanwhile, a same size ROI of the health tissues was selected at the corresponding position on the uninjured side **(Figure 2B)**.

After segmentation of the tissues, the ultrasound characteristics of the images were assessed. The mean grayscale values of the HO tissue was quantitatively evaluated in comparison with the normal muscle tissue.

The segmentation and calculation of mean grayscale value were performed by a self-developed program using Matrix Laboratory (Matlab, version 2016b).

FIGURE 2 | Ultrasound greyscale images of the injured left shoulder joint (A) and the uninjured right shoulder joint (B) of participant No. 2. Red profile in (A) and (B) represents the HO segmented from the surrounding soft tissues and the normal muscle tissue selected on the corresponding position, respectively.

#### Statistical Test

The results of the pilot study showed that there was no significant difference between the grayscale values of the HO tissues in ultrasound images recorded using the two ultrasound systems. Therefore, the effect of ultrasound system on the measured grayscale values was not considered in this study.

The grayscale values were expressed as means and standard deviations (mean ± SD). Due to the small number of the participants, the non-parametric statistical tests were utilized in this study. The Mann–Whitney U-test was performed to test the statistical significance in the grayscale values between the HO and health muscle, and Friedman test was performed to test the statistical significance in the grayscale values measured on different time spots during the HO development period. The significance level was set at 0.05. This statistical analysis was performed by using the Statistical Package of Social Sciences (SPSS, version22, USA).

### RESULTS

The results of plain radiography showed that seven participants were diagnosed with the post-traumatic mature HO and

participant No. 8 was with the immature bone formation. Therefore, for participant No. 8, to track the alterations of the immature bone tissue during the development of HO, the followup ultrasound scans and plain radiography were performed.

### HO in the Immature Bone Formation Stage

The participant No. 8 was in hospital because he suffered swelling and limited range of motion of left knee joint 4 months after surgery for trauma. **Figure 3A** shows discontinuous, faint, and poorly demarcated calcification indicating the immaturity of HO. Approximately 1 month after the first plain radiograph, the second plain radiograph showed that the HO tissues became continuous and distinguishable indicating the more ossified HO tissues, but still not mature enough for HO excision (**Figure 3B**).

In comparison with the plain radiographs, ultrasound images clearly showed the alterations in the anatomic, morphological structure of the involved soft tissues and HO tissues during the rehabilitation procedure (**Figure 3C**). From the transverse ultrasound images, mixed hyperechoic and hypoechoic areas in the swelling affected muscle tissues and loss of the texture pattern of muscular fibers in the early stage of immature bone formation compared with the health muscle tissues on the uninjured joint. During the development of HO, it was found that the echogenicity and homogenesis of the HO tissues increased as the ossification increased, and an acoustic shadowing caused by the HO tissues reduced the smoothness of the profile of the femur. Similarily, the longitudinal ultrasound images clearly visualized the development of the immature HO. During rehabilitation procedure, the discontinuous hyperechoic calcified tissues grew to become continuous lamellar HO tissue (**Figure 3C**). It was noted that the HO tissues in ultrasound image (**Figure 3C**) was consistent with lamellar bone in plain radiograph (**Figures 3A,B**). Ultrasonography is more sensitive to visualize the involved soft tissue and immature bone formation.

### HO in the Mature Bone Formation Stage

The ultrasound images of the injured joints show the irregular muscular architecture and the calcified foci (HO) (**Figure 4**). As **Figure 4C** shown, irregular hyperechoic areas and loss of deep fascia or aponeruroses in the swelling affected muscles were found in three participants No. 1, 5, and 8. Besides HO tissues, loss of the textural structure of the muscular fibers in the swelling affected muscles was found in participant No. 4 (**Figure 4E**). **Figure 4G** shows local muscle evagination due to the HO tissues in the muscles disturbing the muscle structure and causing worse pain in participant No. 6. With less injured muscles in two participates No. 2 and 3, the effect of the mature HO tissues on the texture of muscular fibers sometimes was not obvious (**Figure 4A**). The characteristics of mature HO in ultrasound image might include small calcified foci in size but with high echogenicity and sometimes with acoustic shadow. **Figures 4B,D,F,H** show the normal textural structure and echogenicity of the muscles and bones.

### Analysis of Grayscale Value

**Figure 5A** shows that at the middle stage of the in immature bone formation the grayscale value (86.40 ± 13.42) of the immature HO already increased significantly (p < 0.001) compared with the health muscle tissue (55.1 ± 12.01). At the late stage of HO, it was found that the grayscale value of the mature HO increased greatly (119.09 ± 22.70) and was significantly higher (p < 0.001) than that of to the health soft tissues (70.06 ± 31.09).

FIGURE 4 | Ultrasound images of the mature HO in the different joints. (A) Ultrasound image of the left knee joint of participant No. 3 showing that HO (arrow) with an acoustic shadowing ( ). (C) Ultrasound image of the right knee joint of participant No. 5 showing a hyperechoic area ( ) in the muscles and the HO tissues (arrows). (E) Ultrasound image of the right elbow joint of participant No. 4 showing the hypoechoic HO tissues (arrows) and loss of the textural structure of muscular fibers. (G) Ultrasound image of the right elbow joint of participant No. 6 showing local muscle evagination (dotted line box) by the hyperechoic HO tissues (arrows). (B,D,F,H) Ultrasound images of corresponding position of the uninjured side in (A,C,E,G), respectively.

**Figure 5B** shows that the gradual increase in the grayscale value of the immature HO tissues during the development period. The grayscale value of the mature HO for 5th ultrasound scan (92.18 ± 8.37) was significantly increased in comparison with the 1st, 2nd, 3rd, and 4th scans (p < 0.05).

#### Effect of HO on Rehabilitation Therapy

For the patients with high HO maturity, the strength of the rehabilitation therapy needs to be strengthened appropriately, in order to improve the joint movement function. However, for the patients with low HO maturity, the strength of rehabilitation treatments such as joint loosening and drafting should be soft in order to avoid locally strained muscles and other soft tissues, which might lead to HO worsened.

According to the ultrasound images, the orthopedist obtained the information of HO development of participant No.8. The individualized rehabilitation therapy was performed. The treatment intensity of the traditional rehabilitation training reduced. The participant was anesthetized and underwent with manipulation and arthroscopic surgery. After ∼1-month rehabilitation therapy, the range of extension of the injured knee joint was 5–10, while the range of flexion increased 40◦ reaching up to 95◦ . The muscle strength also slightly increased from grade 4 to grade 4+. Furthermore, with the increase in maturity of HO, no obvious increase of HO size was observed via ultrasonography and plain radiography examinations.

### DISCUSSION

### HO Muscle and its Neuromuscular Function

Although the exact mechanism of HO in traumatic and neurogenic HO is unknown, two common factors precede the formation of HO (19). One factor is trauma or neurological injury. The other is the tissue expression of bone morphogenetic proteins (BMPs), which induce bone formation. In this study, all the participants sustained bone fracture at different joints with HO. We found that the morphologic and textural pattern of HO muscle changed in comparison with normal muscle. Differently, a recent study on spastic muscle induced by spinal cord injury did not find obvious changes in textural pattern of the involved muscle but only shortening of muscle fibers in ultrasound image (20). This may suggest that ultrasonography visualizing the damage in the textural pattern of the involved muscle is able to diagnose HO tissue.

Further, previous study reported that HO formation was related to a series of changes within not only muscles but also nerves and vessels (21). Those alterations may affect the neuromuscular function. After traumatic injury, the sensory nerves detect any damage to the bone, for example, the alterations due to trauma and BMPs. Then, the sensory nerves signal to the central nervous system to start the remodeling program. However, new bone is generated in incorrect places, and HO occurs. Therefore, the regulation of peripheral nerves system is involved in the formation of HO (21, 22).

### Quantitative Assessment of Grayscale Value

Ultrasound images show characteristics, such as size and morphological structure, of the HO tissues. However, these characteristics are individually different in the different joints of different patients. Ultrasound echoes reflect the acoustic impedance, density, stiffness of the tissues. The magnitudes of the echoes are displayed in ultrasound image by grayscale values. During occurrence of HO, primitive mesenchymal cells transform into osteoblasts in the late stage. Previous studies demonstrated that the alterations in the properties of the tissues led to the changes in the ultrasound echoes reflecting the relative diseases (8–11). Therefore, this study chose grayscale value as quantitative parameter, which is relative to the maturity degree of the HO tissues (23). The results show the significant difference in grayscale values not only between the HO tissue and the health muscle tissue but also between immature and mature HO tissues.

FIGURE 6 | The HO tissues at left elbow joint of participant No. 4 visualized by multi-modality imaging. (A) A spiny HO (arrow) in plain radiograph. (B) A slender HO and two small and granular HOs discovered in CT image (arrows). (C) Two small hyperechoic masses and a slender hyperechoic region (arrows) in triceps brachii in ultrasound image.

### Radiographic Classification and Evaluation for HO

Previous studies evaluated and classified HO according to the plain radiograph. Based on radiographic findings, HO in the hip joint is grouped into four classes (Class I to Class IV), called as the Brooker classification (24). The Hastings and Graham classification system (25) was proposed to evaluate HO at the elbow joint into three classes. However, these classification systems were only applied to individual joint and did not obtain good classification results of HO in other anatomical locations.

A recent study proposed an analog scoring method (26) for radiographic classification and evaluation of HO based on normotopic reference bone. Its results showed high correlation (R <sup>2</sup> = 0.89) between the scores of the analog scale and the heterotopic bone volumes measured by micro-CT, i.e., higher score means larger size of HO. In comparison with the analog method, this study evaluated grayscale value of HO in ultrasound image, indicating that higher echogenicity means more mature HO.

#### Limitation of Plain Radiography

In department of orthopedics in a hospital, plain radiography is considered as the gold standard of the clinical diagnosis of HO because of its convenience to detect the calcified tissues. However, plain radiography is a method of projection imaging producing two-dimensional images with x-ray radiation through the body tissues. Therefore, the plain radiograph does not provide any depth information. This limitation leads to misjudgment on HO. In this study, participant No.4 suffered from pain and limited range of motion in the left elbow joint. The radiography revealed a lamellated and slender calcification that intersected the distal humerus (**Figure 6A**). However, 3-D CT image visualized a slender HO and two small and granular HO tissues (**Figure 6B**). Similarily, ultrasonography also displayed two small hyperechoic calcified foci and a lameller hyperechoic region in triceps brachii (**Figure 6C**). It was found that the overlap of the HO tissues and distal humerus or other tissues in the radiograph might cause the misdiagnosis of the HO. Other previous studies slimilarily demonstrated that 3-D CT image help plain radiograph accurately locate the position and number of the HO tissues (27, 28).

### Advantages and Limitations of Ultrasonography in Detection of HO

Previous studies have proven that ultrasound is a useful tool and provides an earlier diagnosis of neurogenic HO than radiography (12–14). This study similarily showed that ultrasonography could depict the changes of related soft tissues and immature HO tissues, that plain radiograph might fail to imagine in orthopedic rehabilitation. The hypoechoic area and loss of the lamellar pattern of muscular fibers in ultrasound image (**Figure 4E**) probably due to the inflammatory edema surrounding HO. Hyperechoic holistic muscle, especially hyperechoic muscular fibers (**Figure 4C**) might be caused by the lack of exercise or amyotrophy.

In addition, ultrasonography could be conveniently performed a bedside monitoring using portable or wireless ultrasound device. There was no significant difference between two devices for imaging and quantitatively evaluation of HO. This study demonstrates that grayscale value analysis method for evaluating HO is independent on ultrasound system. Finally, ultrasonography is suitable for follow-up tracking the development of HO due to its real-time imaging and no radiation.

However, operator dependency limits the use of ultrasonography (29), especially in department of orthopedics. Most orthopedists themselves have no experiences of performing ultrasound scanning and diagnosing HO using ultrasound image. And orthopedists are lack of the cooperation with ultrasound specialist. Another limitation of this study is the number of the proper participants. No participant with early stage of post-traumatic HO was included in this study, only one participant with the immature HO in the middle stage was involved.

### CONCLUSION

In this study, ultrasonography was applied to visualize the development of immature HO and surrounding soft tissues. After the HO extracted from the ultrasound images, the grayscale value was used to quantitatively assess the immature and mature HO tissues in the middle and late stage. The results show the significant difference in grayscale values not only between the HO tissue and the health muscle tissue but also between immature and mature HO tissues. This study suggested that ultrasonography has potentials to be a useful imaging modality for monitoring the development of HO and providing quantitative evaluation on HO. Combination of ultrasonography and plain radiography in diagnosis of HO help orthopedists to make individualized rehabilitation therapy.

### AUTHOR CONTRIBUTIONS

QW and XW conceived and designed the study. PZ and PL performed the experiments. QW and PZ wrote the paper. XS, XL,

#### REFERENCES


and HH contributed to experiments. XW and WC reviewed and edited the manuscript. All authors had read and approved the manuscript.

### FUNDING

This work was supported in part by the Project of Science and Technology Department of Guangdong province (No. 2016A020216017 and No. 2013B021800039), the National Natural Science Foundation of China (No. 81371560), and Natural Science Foundation of Guangdong Province, China (No. 2014A030313329, No. 2015A0303 10527).

### ACKNOWLEDGMENTS

The authors would like to thank all the participants of this study.


and questionable, Bone (2018) 109:147–52. doi: 10.1016/j.bone.2017. 08.011


injury. Bone Joint J. (2015) 97-B:899–904. doi: 10.1302/0301-620X.97B7. 35031

29. Ranganathan K, Hong X, Cholok D, Habbouche J, Priest C, Breuler C, et al. High-frequency spectral ultrasound imaging (SUSI) visualizes early posttraumatic heterotopic ossification (HO) in a mouse model. Bones (2018) 109:49–55. doi: 10.1016/j.bone.2018.01.034

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wang, Zhang, Li, Song, Hu, Li, Chen and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Degraded Synergistic Recruitment of sEMG Oscillations for Cerebral Palsy Infants Crawling

Zhixian Gao<sup>1</sup> , Lin Chen1,2, Qiliang Xiong<sup>1</sup> , Nong Xiao<sup>3</sup> , Wei Jiang<sup>3</sup> , Yuan Liu<sup>3</sup> , Xiaoying Wu1,4 \* and Wensheng Hou1,2,4 \*

*<sup>1</sup> Key Laboratory of Biorheological Science and Technology of Ministry of Education, Chongqing University, Chongqing, China, <sup>2</sup> Collaborative Innovation Center for Brain Science, Chongqing University, Chongqing, China, <sup>3</sup> Department of Rehabilitation Center, Children's Hospital of Chongqing Medical University, Chongqing, China, <sup>4</sup> Chongqing Medical Electronics Engineering Technology Research Center, Chongqing University, Chongqing, China*

#### Edited by:

*Xu Zhang, University of Science and Technology of China, China*

#### Reviewed by:

*Eugene Golanov, Houston Methodist Hospital, United States Rong Song, Sun Yat-sen University, China*

#### \*Correspondence:

*Xiaoying Wu wuxiaoying69@163.com Wensheng Hou w.s.hou@cqu.edu.cn*

#### Specialty section:

*This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology*

Received: *19 May 2018* Accepted: *22 August 2018* Published: *18 September 2018*

#### Citation:

*Gao Z, Chen L, Xiong Q, Xiao N, Jiang W, Liu Y, Wu X and Hou W (2018) Degraded Synergistic Recruitment of sEMG Oscillations for Cerebral Palsy Infants Crawling. Front. Neurol. 9:760. doi: 10.3389/fneur.2018.00760* Background: Synergistic recruitment of muscular activities is a generally accepted mechanism for motor function control, and motor dysfunction, such as cerebral palsy (CP), destroyed the synergistic electromyography activities of muscle group for limb movement. However, very little is known how motor dysfunction of CP affects the organization of the myoelectric frequency components due to the abnormal motor unit recruiting patterns.

Objectives: Exploring whether the myoelectric activity can be represented with synergistic recruitment of surface electromyography (sEMG) frequency components; evaluating the effect of CP motor dysfunction on the synergistic recruitment of sEMG oscillations.

Methods: Twelve CP infants and 17 typically developed (TD) infants are recruited for self-paced crawling on hands and knees. sEMG signals have been recorded from bilateral biceps brachii (BB) and triceps brachii (TB) muscles. Multi-scale oscillations are extracted via multivariate empirical mode decomposition (MEMD), and non-negative matrix factorization (NMF) method is employed to obtain synergistic pattern of these sEMG oscillations. The coefficient curve of sEMG oscillation synergies are adopted to quantify the time-varying recruitment of BB and TB myoelectric activity during infants crawling.

Results: Three patterns of sEMG oscillation synergies with specific frequency ranges are extracted in BB and TB of CP or TD infants. The contribution of low-frequency oscillation synergy of BB in CP group is significantly less than that in TD group (*p* < 0.05) during forward swing phase for slow contraction; however, this low-frequency oscillation synergy keep higher level during the backward swing phase crawling. For the myoelectric activities of TB, there is not enough high-frequency oscillation recruitment of sEMG for the fast contraction in propulsive phase of CP infants crawling.

Conclusion: Our results reveal that, the myoelectric activities of a muscle can be manifested as sEMG oscillation synergies, and motor dysfunction of CP degrade the

synergistic recruitment of sEMG oscillations due to the impaired CNS regulation and destroyed MU/muscle fiber. Our preliminary work suggests that time-varying coefficient curve of sEMG oscillation synergies is a potential index to evaluate the abnormal recruitment of electromyography activities affected by CP disorders.

Keywords: muscle synergy, cerebral palsy, sEMG oscillations, infants crawling, synergistic recruitment

### INTRODUCTION

Cerebral palsy (CP) is a permanent movement disorder caused by brain injury in early childhood with a high prevalence of 2–3 per 1,000 live births (1). The regulation of muscle coordination and motor unit (MU) recruitment are usually affected due to the secondary musculoskeletal morphology alterations and changes in the electrophysiological characteristics of MU followed by cerebral injury (2). Surface electromyography (sEMG) signals regulated by central nervous system (CNS) have been widely used to evaluate MU recruitment patterns and neuromuscular function, and it is widely accepted that motor dysfunction of cerebral palsy destroys the synergistic recruitment of muscular activities (3). Although the impact of CP on synergistic electromyography activities has been observed among the muscle group for limb movement (4), very little is known about how motor dysfunction of CP affects the organization of the myoelectric frequency components due to the abnormal motor unit recruitment patterns.

Various features of sEMG signals have been proposed for assessing the abnormal neuromuscular functions in CP patients. The sEMG characteristics from individual muscle, including parameters in time- and frequency-domain, can be used to measure the abnormal overall sEMG activity of CP. Researchers found the changes of magnitude, frequency content and timing from sEMG signals of individuals with CP. The time- and frequency-domain parameters were also utilized to analyze symmetry, cadence, or smoothness of muscle activity for CP patients (5–7). On the other hand, sEMG characters extracted from multiple muscles, such as muscle synergy and co-activation index, can provide the information about relative muscle activation of muscle groups and insights into motor control. Tang et al. (3) reported that CP recruited fewer muscle synergies during upper limb movements with simplified neuromuscular control strategy. Gross et al. (8) found that coactivation index of the rectus femoris/semitendinosus couple was more sensitive to speed, which could be explained by altered motor control. Furthermore, increasing studies suggested that the motor dysfunction caused by CP affected not only the overall sEMG activity but also the sEMG frequency components or oscillations (9, 10).

Recent studies have demonstrated that the sEMG spectrum profile or frequency components correlated with physiological status of movement. Von et al. (11) reported that high-frequency components of sEMG from tibialis anterior were activated before heel strike during running, while low-frequency components were dominated after heel strike. Liu et al. (12) declared that the highest frequency components of sEMG were more sensitive to muscle fatigue than the raw sEMG signal. Moreover, coherence between two sEMG oscillations was used to evaluate the neural control of movement under different conditions, such as low-alcohol and fatigue (13, 14). In addition, studies also reported that coherence can be affected by voluntary force (15). More recent studies employed sEMG oscillation components for abnormal motor functions analysis, and found that the characteristics of sEMG oscillation components correspond to various neuromuscular damages in CP (10, 16). On other hand, sEMG frequency spectrum is correlated the types of recruited MUs, Wakeling et al. have demonstrated that slower and faster MUs in muscles indeed generate low and high sEMG frequency components, respectively (17, 18). Moreover, sEMG oscillations with specific frequency ranges have been successfully used to evaluate the recruitment patterns of corresponding types of MUs in both animal and human muscles during movement (19, 20).

It is generally accepted that multiple element synergistic organization is a fundamental strategy for motor control, and the elements can be well-organized by CNS to perform limb movement (21, 22). A variety of methods have been employed to extract multiple oscillations from sEMG (12, 23, 24), and multivariate empirical mode decomposition (MEMD) shows better performance in decomposing multi-channel sEMG signals due to its self-adaptability, fully data-driven and generalized multivariate extension (16, 25). These studies have indicated that the resulted multi-scale oscillations are mode-aligned across channels, and analysis of multiple oscillations is at the same scale. To obtain the synergistic pattern of multiple elements, although component-based algorithms, such as principal component analysis (PCA) (26, 27) and independent component analysis (ICA) (28), have been introduced to extract synergies, nonnegative matrix factorization (NMF) is a better choice (29, 30) due to its non-negative constraint for all the matrices (the raw matrix and the obtained matrices). A few works have demonstrated that NMF can be used to extract muscle synergies from multi-channel sEMG signals in upper limb muscles during motion tasks (31), even to obtain activation patterns of muscletendon units and time-vary coefficient curves from high-density sEMG signals during dynamic motion tasks (32).

Inspired by muscle synergy, different types of MUs might also be manifested as specific oscillation synergy patterns. In other words, the synergistic recruitment mechanism should be manifested with coordinated activation of muscle groups and of different types of MUs in individual muscles. Hence, we assumed sEMG frequency components or oscillations might also be recruited with a synergistic pattern for normal motor function, whereas motor dysfunction of CP may affect the organization patterns of oscillation components. To this end, we recorded the sEMG signals from biceps brachii (BB) and triceps brachii (TB) muscles of the upper limbs during infant crawling, extracted sEMG oscillations and analyzed their synergy patterns to explore the abnormal organization of sEMG oscillations in CP.

#### METHOD

#### Participants

Under the approval of children's hospital of Chongqing medical university, 17 infants with typical development (TD, 11.4 ± 1.7 months, 12 males and 5 females) and 12 infants with cerebral palsy (CP, 22.3 ± 5.5 months, 9 males and 3 females) were recruited in this study. Cerebral palsy infants were collected from the department of rehabilitation center in the children's hospital of Chongqing medical university. The TD infants were recruited from local child health clinics. The inclusion criteria for CP infants included: (1) all CP infants were under the age of 3 years old; (2) crawling continuously on their hands and knees during the test; (3) no other diseases that lead to motor function deficits according to the historical records. The TD infants were born at term with normal birth weight and had no neurological impairment. To make sure all subjects can crawl continuously on hands and knees, the Gross Motor Function Measure (GMFM-88) was conducted before experiments. Informed consent forms were obtained from participants' legal guardians before experiments.

#### Experimental Protocol

The infants were encouraged to crawl at their own pace from one end of a mat (360 × 120 cm) to the opposite side. Valid crawling cycle sequence requires infants to crawl continuously on hands and knees more than four cycles. Before the experiment, infants crawled on the mat several minutes to warm up.

As the major activated muscles of the upper limb during crawling, bilateral BB and TB were selected for sEMG recording. After disposable surface electrodes were attached to the muscle belly and bandaged to reduce motion artifacts (**Figure 1**), sEMG signals were collected using a sEMG system (ME6000T8, Mega Electronics Ltd, Finland) with a 15–500 Hz bandwidth and a 1,000 Hz sampling rate.

To specify the crawling phase, 14 markers were attached on the bony landmarks of the bilateral wrist, elbow, shoulder, hip, knee, ankle, the right spine scapula, and the trunk, respectively (**Figure 1**). The kinematic data of infants were recorded at 100 frames/s by a 3D motion capture system (Raptor-E, Motion analysis corporation, USA) with six high-speed digital cameras. Kinematic data and sEMG data were stored in a desktop computer and a laptop computer respectively synchronized offline.

### Data Analysis

As shown in **Figure 2**, to study the organization of sEMG oscillations during crawling, multivariate empirical mode decomposition (MEMD) was employed to extract the multi-scale myoelectric oscillations, and non-negative matrix factorization (NMF) was used to analyze the patterns of synergistic organization within sEMG oscillations. Then, the recruitment coefficients within each refined crawling phase were evaluated.

#### Pre-processing

Firstly, the valid data of more than four consecutive crawling cycles were segmented from raw signals, and sEMG for the BB and TB of both sides and kinematic data were simultaneously recorded (see **Figure 3**). Then, the valid sEMG signals were processed with a zero-lag high-pass filter (4th order Butterworth filter, 20 Hz). Then, a 50 Hz notch filters were adopted to remove power frequency noise.

#### Oscillation Extraction With MEMD

EMD is a fully data-driven analysis approach which selfadaptively decomposes a non-linear and non-stationary signal into several intrinsic mode functions (34). However, for multichannel signals, applying EMD to each channel could produce a different number of misaligned IMFs. Rehman et al. (35) proposed the MEMD method to extend EMD to multivariate signal decomposition. They treated n-variable time series as ndimensional vectors and employ the spherical coordinate system to project n-dimensional vectors along different directions in (n-1)-dimensional space, and the mean value of the envelopes of these n-dimensional projection sequences is obtained as the local mean of multiple time series. After MEMD decomposition, the ndimensional raw signal {v (t)} T <sup>t</sup>=<sup>1</sup> = {v<sup>1</sup> (t), v<sup>2</sup> (t), · · · v<sup>n</sup> (t)} can be decomposed into d layers multivariate IMFs hi(t) and residual r(t), which can be described as

$$\mathbf{v}\left(t\right) = \sum\_{i=1}^{d} h\_i\left(t\right) + r\left(t\right) \tag{1}$$

The IMFs obtained after MEMD analysis have the same number for each channel and orderly align scales across channels (35).

Referred literature (16), we firstly concatenated the valid sEMG signal segments from all subject as a 4-channel sEMG

dataset D, which is shown as following formula:

$$D = \begin{bmatrix} CH\_{11} - CH\_{12} - \cdots - CH\_{1j} \\ \vdots \\ CH\_{i1} - CH\_{i2} - \cdots - CH\_{ij} \end{bmatrix} \tag{2}$$

where CHij is the ith channel of the valid sEMG signal segments for the jth subject. Here, we recorded four channel sEMG signals, and totally 29 subjects conducted the test; that is to say, i = 1, 2, 3, 4; j = 1, 2, ... , 29. To reduce the mode mixing in multivariate IMFs (36), 2-channel Gaussian white noise (with the same length as D) were added to the 4-channel sEMG dataset D to composite a 6-channel dataset. Then, the composite dataset was decomposed by MEMD, which yield a total of 23 IMFs with aligned scales across channels, cycles and subjects. Example of raw sEMG signals and their first nine IMFs of BB in a crawling cycle are illustrated in **Figure 4**, of which the right panel show the corresponding power spectra for raw sEMG and IMFs. The frequency ranges of the first nine IMFs are listed in **Table 1**, the frequency band ranges of IMFs decrease orderly, and corresponding IMFs for different subjects are scalealigned.

TABLE 1 | Frequency ranges of first nine IMFs determined by 3dB bandwidth.


#### Kinematics Data Preprocessing and Crawling Phase Segmentation

Kinematics data were processed with a zero-lag low-pass filter (4th order Butterworth filter, 4Hz) to remove high-frequency noise. As shown in **Figures 3A,B**, crawling cycle phases was determined with the z coordinates of the wrist (ZW) and the shoulder joint angle (the joint angle in sagittal plane) (37). **Figure 3C** showed the examples of segmentation from two subjects during one crawling cycle. The squared of time derivative of Z<sup>W</sup> (v 2 ZW , velocity squared) was applied to segment swing and stance, and the crawling cycle begins with swing. A threshold of v 2 ZW was set at 0.5 (m<sup>2</sup> /s2 ) to decide the onset of limb moving and the end moment of a crawling phase (38). To refine the crawling cycle in detail, the swing phase was divided into forward swing phase (FSP) and backward swing phase (BSP) with the maximum of ZW; and stance was divided into braking phase (BRP) and propulsive phase (PRP) with the minimum of time derivative of A<sup>S</sup> (wS).

#### Synergistic Recruitment Analysis of sEMG Oscillations

As the frequency of IMF9 is below 20 Hz and out of sEMG frequency range (20–500 Hz), the first eight IMFs are included for oscillation recruitment analysis. Envelopes of the first eight IMFs of each channel sEMG signal were extracted (Hilbert spectrum and median filtering) and segmented into cycles according to kinematics data. The envelopes were re-sampled into 200 points for each cycle. The synergies of sEMG oscillation were extracted with NMF, which decompose a non-negative matrix V into two non-negative matrixes including a base matrix W and corresponding recruitment coefficient matrix C (39). Here, eight envelopes of IMFs resulted from a channel of sEMG within a crawling cycle composed the non-negative matrix V. The NMF method can be described as following

$$\mathbf{V}^{m \times n} = \mathbf{W}^{m \times s} \mathbf{C}^{s \times n} \tag{3}$$

In this study, the matrix of V m×n represents the envelopes of m IMFs (m = 8), and n is the sample number (n = 200). Each column of Wm×<sup>s</sup> represents an oscillation synergy with m weight factors (1 ≤ s ≤ m), and s is the number of synergies. Each row of C s×n represents corresponding recruitment coefficients, which shows how each synergy is modulated over time. To maintain the modulation of information within oscillation synergies, C was normalized to its maximum values of the crawling cycle.

The number of oscillation synergies was determined by calculating the variation accounted for (VAF) between the original matrix V and the reconstruction matrix V ′ = WC for each s value (from 1 to 8) (40). VAF is calculated as follow

$$\text{VAF} = 1 - \frac{\left(V - V\right)^2}{V^2} \tag{4}$$

Here, the selection criteria were that the mean of VAFs was larger than 95% and the increment of VAF was <1%.

#### Statistical Analysis

The recruitment coefficient values for each oscillation synergy were averaged over 3–4 valid cycles for each subject. For each TD infant or CP infant whose both side limbs were affected, as the symmetric movement for left and right side, the oscillation synergy and recruitment coefficient of sEMG signals recorded from right and left BB and TB muscles have been averaged firstly. In order to compare the recruitment pattern of each oscillation synergy, independent sample t-test was applied on the recruitment coefficients for each phase. The significance level was set to 0.05. All the data were analyzed in SPSS 22.0 statistical software (SPSS Inc., USA). In addition, Pearson's correlation coefficient (r) between any two oscillation synergies with the same order number was calculated to assess the similarities of oscillation synergy structures. Four kinds of r values were calculated between same within-groups muscles (e.g., BB of TD1 vs. BB of TD2), between different within-groups muscles (e.g., BB of TD1 vs. TB of TD2), between same between-groups muscles (e.g., BB of TD1 vs. BB of CP1) and between different betweengroups muscles (e.g., BB of TD1 vs. TB of CP1).

#### RESULTS

All of the recruited infants in this study can crawl continuously on hands and knees more than 4 valid cycles (TD, 4.9 ± 0.9 cycles; CP, 5.9 ± 1.4 cycles). Their hands-and-knees crawling scores were at least 44 (TD, 50.9 ± 3.3; CP, 63.7 ± 8.5) according to GMFM. Although the crawling motor function scores of CP infants equal or even exceeded TD infants (in this study) after rehabilitation training, their motor development was indeed slower than that of TD infants at the same age. CP infants exhibited abnormal movement behavior and sEMG activities.

### Three Stable Oscillation Synergies With Low-, Medium- and High-Frequency Ranges in Muscles

The mean VAF values acquired in the NMF decomposition of the sEMG oscillations are shown in **Figure 5**. For bilateral BB and TB muscles, the number of oscillation synergies (s) was chosen as 3 according to the selection criteria above.

**Figure 6** shows the three oscillation synergies composed of the IMFs with different weights for BB and TB in TD and CP group. Each synergy of sEMG oscillation owned a dominant frequency range. The main contributors of synergy1 were IMF1, IMF2, IMF3, and IMF4, and their frequency ranged from 81 to 500 Hz (see **Table 1**); synergy2 was dominated with IMF4, IMF5 and IMF6, which located the frequencies from 34– 189 Hz. Synergy3 was comprised of oscillations with much lower frequency, IMF6 and IMF7, and the corresponding frequency range was 23–69 Hz. Furthermore, as listed in **Table 2**, the Pearson's correlation coefficient of synergy 1 was more than 0.883, which indicated that the composition of IMF1∼IMF8 exhibited high structural similarity for any paired muscles. The composition of IMF1∼IMF8 in synergy 2 or synergy 3 showed high similarity for the BB and TB of CP and/or TD group as well, and corresponding correlation coefficients were more than 0.816 and 0.861, respectively. Altogether, there were three stable oscillation synergies of sEMG for BB and TB both in CP infants or TD infants crawling, and each synergy can be characterized with sEMG oscillation composition (i.e., IMF) of specific frequencies.

### Dynamic Recruitment of Oscillation Synergies for Crawling Movement

The recruitment coefficient curves of sEMG oscillation synergies are illustrated in **Figure 7**. During a crawling cycle, each of those three oscillation synergies was dynamically recruited with time-varied coefficient curves. It can be observed, either for TD or CP infants, that they recruited the synergy 1 in BB with a low level firstly, then the recruitment reached a peak level in the second half crawling phase. Whereas the recruitment strength of synergy 2 and synergy 3 reached a peak quickly at the beginning of crawling cycle (i.e., FSP) and decreased then. For the TB, both TD and CP infants almost synchronously recruited those three synergies, and their activation level reached a peak in the middle phase of crawling (i.e., BSP and BRP). However, the motor dysfunction of CP affected the recruitment of oscillation synergies for crawling movement both in BB and TB.

#### The Impact of CP on the Recruitment Pattern of Oscillation Synergies

To quantitatively evaluate the impact of CP on the recruitment pattern of oscillation synergies, the amplitude of recruitment coefficient curves was compared in the refined crawling phases of FSP, BSP, BRP, and PRP. As shown in **Figure 8**, CP recruited significant less synergy3 or low-frequency sEMG oscillations in BB for slow contraction in FSP crawling (**Figure 8C**, p = 0.023), and less synergy 1 or high-frequency sEMG oscillations in TB for the fast contraction in PRP crawling

FIGURE 6 | Oscillation synergies extracted from 8 sEMG IMFs of BB and TB during one crawling cycle in two groups. (A) BB of TD group; (B) BB of CP group; (C) TB of TD group; (D) TB of CP group.

(**Figure 8D**, p = 0.007). During the BSP (p = 0.007) and BRP crawling (p = 0.006), CP maintained a high level of synergy3 in BB, which revealed that the sEMG oscillations with low-frequencies cannot be de-recruited effectively (**Figure 8C**). Furthermore, during CP crawling in FSP (p = 0.002) and BRP (p = 0.047), there was a significantly higher activation level for synergy 2 in BB (**Figure 8B**), whereas less synergy 2 was recruited in TB for BSP crawling (**Figure 8E**, p = 0.029).

FIGURE 5 | Mean VAF values corresponding to different number of oscillation synergies.

### DISCUSSION

To deeply explore the effect of neuromuscular damage of CP on the motor regulation, this study aimed to extract and evaluate the oscillation synergy patterns from sEMG signals with MEMD and NMF during infant crawling. The present work showed that sEMG signals contained stable structures of oscillation synergy, and the motor dysfunction of CP affected the recruitment of oscillation synergies for crawling movement both in BB and TB.


### Synergistic Recruitment of sEMG Oscillations for Crawling Movement

Crawling requires multiple skeletal muscles to participate in the regulation of limb joint flexion and extension. Muscles play different roles and their contraction patterns also vary with crawling phases. The results showed that the activation pattern of muscles can be expressed as synergistic recruitment of multiscale sEMG oscillations. As shown in **Figure 6**, the activation of BB or TB muscles during crawling can be represented by the synergistic combination of the oscillation subsets with different frequency ranges, i.e., high-frequency range (synergy1), medium-frequency range (synergy2), and low-frequency range (synergy3). Moreover, as shown in **Figure 7**, these synergistic sEMG oscillations have been dynamically organized for infant crawling movement. Both for BB and TB, the activation level of high-frequency synergy (C1) keeps to a way of progressively developing pattern, whereas the activation level of mediumand low-frequency synergy (C2 and C3) reach a peak and then descend gradually. These results follow the principle that the recruitment of MUs has been typically graded (41) from slow to fast during dynamic contractions (42) for joint movement, and types of recruited MUs can be manifested as sEMG frequency components (17 , 18). In context, the observations in our study suggest that, during the swing and stance of infant crawling, the organization of MUs activation pattern in BB and TB can be represented as synergistic sEMG oscillations, which have been dynamically regulated for different contraction status and crawling phase. Although there is some time delay for BB and TB contraction due to their different role in crawling movement, the coefficient curve for sEMG oscillations with high-, medium-, and low-frequency are alternatively enhanced to a complementary pattern among these synergies, which reveal that different types of MUs have been coordinately recruited for limb movement in crawling.

Crawling is a periodic movement of flexion-extension, in which the upper limbs postures alternate with FSP, BSP, BRP, and PRP successively. Meanwhile, BB and TB exhibit alternating coordinate relaxation and contraction. During swing, BB gradually contracts to achieve forward swing of upper limbs in FSP; slow MUs are activated, and the recruitment strength of synergy 2 and synergy 3 with lower-frequency range reach a peak. However, TB presents some time delay in crawling movement and kept extension/relaxation in FSP, and then begin contraction to swing arm backward with fast developing of high-, medium-, and low-frequency oscillations, and medium- and lowfrequency oscillations reached their peak activations. To perform backward swing, BB extends in BSP as slow MUs are de-recruited, and the activation of high-frequency oscillations exceeded the medium- and low-frequency ones. During stance, both BB and TB are activated to response loading period, and fast MUs are dominantly recruited accordingly to absorb the strike shock and quickly stabilize the kinematic behavior of the joint. As shown in **Figure 7**, the coefficient curve of synergy 1 (C1) maintains a high level, whereas medium- and low-frequency oscillations (C2 and C3) decrease gradually. In other words, more fast MUs for high-frequency oscillations are activated to compensate the

TABLE 2 | Correlation coefficients between oscillation synergies.

Synergy1

> TD

CP

TD

CP

TD

CP

Synergy2

Synergy3


BB

> TD BB

TB CP BB

TB *Data were expressed as mean* ±

*standard error.*

0.896 ± 0.013 0.907 ± 0.014 0.890 ± 0.012

0.931 ± 0.017 0.910 ± 0.013

0.883 ± 0.014

0.968 ± 0.003 0.912 ± 0.015 0.945 ± 0.003 0.921 ± 0.015 0.929 ± 0.006 0.894 ± 0.015 0.932 ± 0.010 0.848 ± 0.038 0.927 ± 0.006 0.892 ± 0.011 0.929 ± 0.010 0.861 ± 0.021

TB

BB

TB

BB

TB 0.897 ± 0.012 0.860 ± 0.022 0.863 ± 0.011

0.961 ± 0.004 0.816 ± 0.053

0.822 ± 0.021

BB

TB

BB

TB 0.893 ± 0.010 0.883 ± 0.061 0.877 ± 0.011

0.938 ± 0.005 0.868 ± 0.018

0.869 ± 0.011

BB

TB

BB; (D) synergy 1 of TB; (E) synergy 2 of TB; (F) synergy 3 of TB. \*0.01 < *p* < 0.05, \*\*0.005 < *p* < 0.01, \*\*\**p* < 0.005.

de-recruited slow MUs in stance phase of crawling, which is a mechanism for CNS timely and moderately regulating the recruitment of different types of MUs in muscles for joint movement (42).

### Abnormal Recruitment Pattern of Oscillation Synergies Under Motor Dysfunction of CP

Motor dysfunction of cerebral palsy originates from brain injury leading to a series of neurological changes, such as reduced input to the CNS (43) and decreased activity of descending inhibitory system (44). Subsequently, motor dysfunction of CP is also manifested as abnormal myoelectric activities and abnormal movement postures due to the affected cerebral nerve in motor area (2). As illustrated in **Figures 6**, **7**, although three types of synergistic sEMG oscillations with different frequency band are presented during CP infants crawling, there are some changes in recruitment coefficient curves, which suggest that brain injury of cerebral palsy affected the MUs recruitment. Generally, CNS employs a synergistic strategy to organize the multi-element motor system within different scales to simplify motor control (45), and the combination of motor elements with certain weights constructs motor synergy. However, motor dysfunction of CP mainly alters the organization of medium- and low-frequency oscillations of sEMG, especially the intensity proportion between medium- and low-frequency oscillations. Compared with TD, CP infant recruits less low-frequency oscillation (synergy3) during FSP crawling, whereas TB activates less medium-frequency oscillation (synergy2) during BSP crawling. It was reported that the impaired CNS of CP is unable to effectively drive lowthreshold MUs (46–48), our results reveal that the inadequate recruitment of low-threshold MUs induced insufficient activation of low- and medium-frequency oscillations when BB and TB of CP infants perform flexion contraction in FSP and BSP crawling, respectively.

In addition to the aforementioned insufficient activation of low- or medium-frequency sEMG oscillations, it also can be observed that motor dysfunction of CP is unable to de-recruit those undesired sEMG oscillations. As shown in **Figures 7**, **8**, for the BB of CP infants, the coefficient curves of synergy3 (C3) and synergy2 (C2) hold on an abnormal level in BRP crawling, and coefficient curve of synergy3 (C3) also keep on an abnormal level in BSP crawling. These results suggest that, during BRP and BSP, CP infants are unable to inhibit their low-frequency and/or medium-frequency sEMG oscillations in the BB muscles. It is accepted that different MU owns different intrinsic contraction properties, and slow MUs are activated for slow contraction and low-intensity activities, while fast MUs are recruited for fast contraction and high-intensity activities (49, 50). Central nervous system selectively recruits desired MUs and de-recruit undesired MUs to produce synergistic contraction pattern, and organizes different types of MUs at a certain activation level can improve the synergistic control for muscular activities and joint movement. For CP infant, motor dysfunction of injured motor nerve or declined activity of descending inhibitory system (31) cause an insufficient organization of different types of MUs; another observation showed that spastic diplegic CP loss the ability to fully recruit and optimally activate available motor units because of their central activation failure (39). So, our work reveal that, the synergistic recruitment of MUs can be characterized with the coordinated activation of different types of MU with appropriate proportion, and CP cannot effectively regulate the activation intensity among different MUs.

### Synergistic Organization of Muscle Functional Units Exhibited by Oscillation Synergy Patterns

The regulation of motor function by the CNS is a synergistic organization of multiple function units, which has been widely confirmed in muscle groups. Variety of methods have been used to evaluate the intensity, spectrum or frequency component of sEMG activity. However, it is still difficult to analyze the organization patterns of different types of MUs. Multiscale oscillation modes in sEMG signals are a combination of MU activation and organization patterns. The results of this study demonstrated that the synergistic organization of different oscillations in sEMG signals can be used to evaluate the organization of MUs. In fact, sEMG frequency components are related to physiological state and have been applied to assess motor function under different conditions (13, 14, 23). However, these studies only provided an overall effect of frequency components in a certain period of time. We decomposed sEMG signals of individual muscles into multiscale oscillations, which is correlated to MU recruitment pattern, and then, utilized NMF to extract oscillation synergies and their recruitment coefficient curves. The time-varying curves represent how CNS regulate oscillation synergies over time, which provided a novel insight to better understand the dynamic process of the CNS organizing different types of MUs.

The results of this study showed that three stable multiscale oscillation synergies might be underlying intrinsic structure of EMG activity during crawling. The coefficient curves can further reveal the abnormal MU activation patterns during crawling in CP. Some researchers have found that motor dysfunction of CP manifested as abnormal entropy of multi-scale oscillations of sEMG activities in CP, and a major peripheral cause of these abnormalities is abnormal MU recruitment (10, 16). This study revealed that motor dysfunction of CP also manifested as abnormal recruitment of oscillation synergies. For CP, the motor dysfunction is caused by original brain injury and resulted abnormal descending motor commands. It can be observed that, in CP group, the oscillation synergies with low-frequency range (synergy3) were insufficiently recruited when BB muscles performed slow contraction during swing, on other hand, they insufficiently de-recruited the oscillation synergies with low-frequency range (synergy2 and 3) when BB extended and relaxed in BSP and BRP. These results suggested the recruitment coefficient curves of oscillation synergies can exhibit the time-varying organization of sEMG oscillations modulated by CNS, and the curves also can reveal abnormal motor control caused by neuromuscular deficits.

#### CONCLUSION

This is the first time to adopt the synergistic mechanism of motor function to characterize the recruitment pattern of sEMG oscillations. The present results reveal that, the myoelectric activities of a muscle can be manifested as sEMG oscillation synergies with different frequency ranges, and motor dysfunction of CP degrade the synergistic recruitment of sEMG oscillations due to the impaired CNS regulation and destroyed MU/muscle fiber. Our preliminary work suggests that, the synergistic organization for motor control can be manifested rather than muscle group, there are oscillation synergies in sEMG signals. Also, the time-varying coefficient curve of sEMG oscillation synergies is a potential index to evaluate the abnormal recruitment of different types of MUs affected by CP disorders. In future work, we would further investigate the sEMG oscillation synergies in other muscles of upper and lower limbs, and the body size of test sample should be enlarged to verify the influence of CP on the synergistic recruitment of sEMG oscillations.

#### ETHICS STATEMENT

This study is under the approval of children's hospital of Chongqing medical university, and all the tests are conducted

#### REFERENCES


in the department of rehabilitation center, children's hospital of Chongqing medical university.

### AUTHOR CONTRIBUTIONS

ZG analyzed the data and signals. ZG, LC, and WH drafted and revised the work. XW, NX, WJ, and WH interpreted the data. ZG, QX, and YL collected the data. XW and WH proposed the research topic. WH organized the data collection, signal processing, and manuscript draft.

#### FUNDING

This work was supported by the National Natural Science Foundation of China (NSFC 31470953, 31771069, 31800824), the Chongqing Science & Technology Program (cstc2016shmszx130060, cstc2015jcyjB0538), and the Fundamental Research Funds for the Central Universities (No. 106112016CDJXZ238826).

#### ACKNOWLEDGMENTS

We much appreciated the infants and their parents for participating in this study. We would like to thank the department of rehabilitation center, children's hospital of Chongqing medical university for recruiting infants and their help in data collection.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Gao, Chen, Xiong, Xiao, Jiang, Liu, Wu and Hou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Inter-Limb Muscle Synergies and Kinematic Analysis of Hands-and-Knees Crawling in Typically Developing Infants and Infants With Developmental Delay

Qi L. Xiong1,2, Xiao Y. Wu1,2,3, Jun Yao<sup>4</sup> \*, Theresa Sukal-Moulton<sup>4</sup> , Nong Xiao<sup>5</sup> , Lin Chen1,2 , Xiao L. Zheng1,2,3, Yuan Liu<sup>5</sup> and Wen S. Hou1,2,3 \*

*<sup>1</sup> Key Laboratory for Biorheological Science and Technology of Ministry of Education, Chongqing University, Chongqing, China, <sup>2</sup> Chongqing Medical Electronic Engineering Technology Research Center, Chongqing University, Chongqing, China, <sup>3</sup> Collaborative Innovation Center for Brain Science, Chongqing University, Chongqing, China, <sup>4</sup> Department of Physical Therapy and Human Movement Sciences, Northwestern University, Chicago, IL, United States, <sup>5</sup> Department of Rehabilitation, Children's Hospital of Chongqing Medical University, Chongqing, China*

#### Edited by:

*Ping Zhou, University of Texas Health Science Center at Houston, United States*

#### Reviewed by:

*Marco Fidel Avila-Rodriguez, Pontificia Universidad Javeriana, Colombia Masahiro Shinya, Hiroshima University, Japan*

\*Correspondence:

*Jun Yao j-yao4@northwestern.edu Wen S. Hou w.s.hou@cqu.edu.cn*

#### Specialty section:

*This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology*

Received: *30 June 2018* Accepted: *26 September 2018* Published: *16 October 2018*

#### Citation:

*Xiong QL, Wu XY, Yao J, Sukal-Moulton T, Xiao N, Chen L, Zheng XL, Liu Y and Hou WS (2018) Inter-Limb Muscle Synergies and Kinematic Analysis of Hands-and-Knees Crawling in Typically Developing Infants and Infants With Developmental Delay. Front. Neurol. 9:869. doi: 10.3389/fneur.2018.00869* Hands-and-knees-crawling is an important motor developmental milestone and a unique window into the development of central nervous system (CNS). Mobility during crawling is regularly used in clinical assessments to identify delays in motor development. However, possible contribution from CNS impairments to motor development delay is still unknown. The aim of this study was to quantify and compare inter-limb muscle synergy and kinematics during crawling among infants at a similar developmental age, however, clinically determined to be typically developing (TD, *N* = 20) infants, infants at risk of developmental delay (ARDD, *N* = 33), or infants with confirmed developmental delay (CDD, *N* = 13). We hypothesized that even though all of the groups are at a similar developmental age, there would be differences in kinematic measures during crawling, and such differences would be associated with CNS impairment as measured by electromyography (EMG) features. Surface EMG of eight arm and leg muscles and the corresponding joint kinematic data were collected while participants crawled on hands and knees at their self-selected velocity. Temporal-spatial parameters and normalized Jerk-Cost (JC) function (i.e., smoothness of movement) were computed from the measured kinematics. The inter-limb muscle synergy and the number of co-activating muscles per synergy were measured using EMGs. We found that the infants with CDD demonstrated higher normalized JC values (less movement smoothness), fewer muscle synergies, and more co-activating muscles per synergy, compared to infants with TD (*p* < 0.05) and ARDD (*p* < 0.05). Furthermore, the normalized JC values were correlated (*p* < 0.05) with the number of co-activation muscles per synergy. Our results suggest a constrained neuromuscular control strategy due to neurological injury in infants with CDD, and such constrain may contribute to the reduced movement smoothness in infant crawling.

Keywords: infant crawling, EMG, muscle synergy, kinematics, motor developmental delay

## INTRODUCTION

Mobility during hands-and-knees crawling is regularly used in clinical assessments to benchmark delay in motor development because it is an early example of skillful gross motor ability. Clinically, a delay can be quantified relative to typical achievement of gross motor milestones (1), but the extent to which such delay is related to CNS impairment cannot be ascertained. Kinetic and kinematic measures during 4-beat crawling can enrich infant assessment and provide a non-invasive window to CNS function.

Previous studies have measured kinetics and kinematics separately in infants crawling. Early studies used film recording to investigate the movement pattern of limbs while crawling from a small sample size of infants (N = 7) (2), while more recent studies used 3D motion capture to examine inter-limb coordination patterns during crawling. The typical pattern of crawling is with diagonal limbs tending to move together and ipsilateral limbs alternating during crawling on hands and knees (3, 4). Quantitative data concerning muscle activities in human infant crawling is sparse. It has been briefly described as triceps brachii activation throughout the stance phase of the arm during crawling, with quadriceps femoris activated during swing phase of the leg (4). Our previous work demonstrated that muscle co-activations of lower extremities during crawling is correlated to their motor skill development (5), but it only quantified the coordination between antagonist muscles of a single limb, which provides little information about the inter-limb coordination across arm and leg muscles during hands-and-knees crawling.

Muscle synergy analysis is a valid tool to explore the coordination across multiple muscles during locomotion (6), and reflects the CNS control for locomotion as a linear combination of several muscle activation patterns in order to complete functional tasks. Crawling is a self-motivated rhythmic locomotion that involves controlled inter-limb muscle coordination for movement. Thus, quantifying inter-limb muscle synergy during crawling has the potential to explore the underlying factors related to movement abnormalities related to the changes/impairments of the CNS.

Muscle synergy extraction based on surface electromyography (sEMGs) and non-negative matric factorization (NMF) algorithm has been used to explore muscle coordination during locomotion in neurotypical populations as well as a number of pathologic conditions, such as stroke (6), spinal cord injury (7), or cerebral palsy (8). For instance, Dominici et al. concluded that two basic muscle synergies are retained through infant development, and are augmented by new synergies during the development of independent walking (9). Steele et al. found that individuals with CP (age range 3.9–70 years) demonstrated fewer synergies during gait compared with unimpaired individuals (8), similar to the constrained muscle control found in adults following stroke (6, 10). Muscle synergy analysis has also been used to quantify the kinetic feature during hand-and-knee crawling. Chen's study (11) extracted two alternative intra-limb muscle synergies during crawling in healthy adults, with one related to the stance phase and the other related to the swing phase (11). However, muscle synergy analysis during infant crawling has not been systematically investigated either in typical development or neurological disorders.

In order to fill the gap, we simultaneously measured kinetic (i.e., EMG) and kinematic features during crawling with typically developing infants and infants with different risks or severities of developmental delay. Because smoothness and well-coordinated movement are typical features of well-developed human motor behavior (12), we expected that kinematic output, such as the smoothness of movement, would be altered in infants with developmental delay. At the same time, we hypothesized that CNS control in infants with developmental delay is also impaired, and would be manifested in the metrics of muscle synergy. Finally, we hypothesized that CNS development/impairment would be associated with the movement smoothness.

## METHODS

### Participants

We recruited 47 atypically developing infants (age range 8–43 months, 14.21 ± 6.91 months; female: N = 20, male: N = 27) from the Department of Rehabilitation Center, Children's Hospital of Chongqing Medical University. Infants visited the hospital to follow up for the risk of developmental delay (13) due to: (1) premature delivery (gestational age<37 weeks); (2) low birth weight (<2,499 grams), regardless of gestational age; or (3) lack of oxygen to the brain during birth. The age of one infant was 43 months, which is far from the distribution of other infants' age and therefore was excluded as an outlier. The remaining 46 infants (age range 8–32 months, 12.78 ± 4.87 months) were included for data analysis in this study. In addition, 20 developmental-age-matched healthy infants (age range 8–15 months, 10.95 ± 2.25 months; female: N = 9, male: N = 11) were recruited from local child health clinics as "typical development (TD)" controls. They were all full-term with normal birth weight, and no diagnosed health conditions per parent report. All infants were studied at the Department of Rehabilitation Center, Children's Hospital of Chongqing Medical University. The experiments were performed with informed, written consent of the parents or guardians of the infants, and the procedures were approved by the ethics committee of Children's Hospital of Chongqing Medical University (approval number: 065/2011). Partial results (i.e., crawling velocity, cadence, stance phase time) from the 20 TD infants have been published before (5).

#### Clinical Assessment

For all of the participants, the Gross Motor Function Measure (GMFM-88) and Gesell Developmental Scale were assessed by specialist physicians. GMFM-88 measures gross motor function, including lying and rolling, crawling and kneeling, sitting, standing, walking activities. Each function is scaled in the range of 0–100 (1). Gesell Developmental Scale is a set of developmental metrics, which assesses the ages and stages of development in young infants (1).

For the infants who were atypically developing, developmental age (see the column 2 in **Table 1**) was assessed by the gross motor development part of Gesell Scale, and compared to


TABLE 1 | Participant demographic information.

\**Determined by Gesell Developmental Scale.*

their biological age to calculate delayed age (in months) for each of the infants (see the column 1 in **Table 1**). Infants with a developmental delay of ≤3 months were classified as at risk of developmental delay and those with a delay larger than 3 months as having confirmed developmental delay. This resulted in 13 infants (age: 20.15 ± 5.85 months, delayed age: 8.65 ± 4.47 months) with confirmed developmental delayed (CDD), and 33 infants (age: 11 ± 2.29 months, delayed age: 0.3 ± 0.63 months) who were at risk of developmental delay (ARDD). Although the biological age of CDD group is larger than that for TD and ARDD groups (F = 41.50, p < 0.01), the developmental age of all the groups are similar (F = 0.072, p = 0.790), demonstrating clinically comparable level of motor skills in all the 3 groups. Demographic information for all participants is summarized in **Table 1**.

#### Protocol

Infants first became acquainted with the laboratory by spending time on a floor crawling mat (size 360 × 120 cm). Next, they were encouraged to crawl from one end to the other in response to toys or mother's calling. After training, infants wore only diapers. A motion capture system (Raptor-E, Motion Analysis Corporation, USA) was used to record kinematic movement of infants at 100 frames/s with six high-speed digital cameras. Fourteen reflective markers were taped over the shoulder (lateral to the acromion), elbow (lateral epicondyle), wrist (ulnar styloid process), hip (posterior superior iliac spine), knee (lateral joint line), ankle (lateral malleolus), and trunk (scapula).

Simultaneously, a surface EMG system (ME6000, Mega Electronics Ltd, Finland, bandwidth of 15–500 Hz) with preamplified EMG sensor units was used to measure sEMG from bilateral arm and leg muscles, including: left and right triceps brachii (LTB, RTB) and biceps brachii (LBB and RBB), quadriceps femoris (LQF and RQF) and hamstring (LHS and RHS) (see **Figure 1**) by differential electrodes. All of the sEMG was sampled at 1 kHz and synchronized with kinematic data recording by a TTL pulse. In addition, movements of participants were videotaped.

A valid trial was defined as straight crawling without stop or deviation, for at least three complete, consecutive strides. In each of the participants, the number of valid trials collected varied from 2 to 16 (on average 6.80 ± 3.60), depending on the cooperation of the infant.

#### Data Analysis

The first and last strides of each valid trial were excluded from the following data analysis.

#### Kinematic Analysis

#### **Temporal-spatial parameters**

Missing raw kinematic data was constructed using cubic spline interpolation. Then they were low-pass filtered (6 Hz) with a zero lag 4th-order Butterworth filter to remove high frequency noise. In the current study, we defined the left wrist as the start of the crawl cycle, similar to the heel strike in gait analysis). The temporal-spatial crawling parameters were accordingly calculated from the 3D trajectories of the left wrist, including velocity, cadence, and stance phase time (normalized to crawling cycle, SPT), using the methods previously reported (5).

#### **Movement smoothness**

Movement smoothness was quantified by evaluating the endpoint jerk-cost (JC) at the left wrist, defined as:

$$\text{JC} = \int\_0^T \left(\frac{d^3s}{dt^3}\right)^2 dt \tag{1}$$

where T is the total duration of a crawling cycle, and s is the position vector of the limb segment. JC measures the change between acceleration and deceleration during movement. A smaller JC value reflects fewer such switches and thus indicates a smoother movement (14).

For each crawling cycle of each subject, the endpoint JC of the wrist marker was calculated in anterior-posterior (AP) direction (JCx), medial-lateral (ML) direction (JCy), and vertical (VT) direction (JCz) using Equation (4). To account for variations in crawling velocity and to standardize results, all JC values were normalized by the total duration of each crawling cycle, T. Per subject, the normalized JC were then averaged across all valid crawling cycles.

#### EMG Analysis

#### **EMG preprocessing**

The sEMG signals were band-pass filtered (10–400 Hz) using a 4th-order, zero-phase Butterworth digital filter and a 50 Hz digital notch filter for reducing the power interference. The filtered sEMG signals were then divided into segments according to the initiation of each crawling cycle. Segments of filtered sEMG signals were subsequently demeaned, rectified, and lowpass filtered with a zero lag 4th-order low-pass (9 Hz) to

extract envelope. The envelope was then normalized to its peak value during each trial, then resampled from 0∼100% of the crawling cycle at the 1% step increase. Finally, the normalized envelopes of per participant and per muscle were averaged cross all valid cycles. The averaged envelopes composed the EMG data matrix (8 × 101) for an individual subject during crawling.

#### **Non-negative matrix factorization**

A non-negative matrix factorization was applied to each EMG data matrix to extract muscle synergies. This method decomposed the measured EMG data matrices (M) into two components, spatial structure (W, termed the muscle synergies) and temporal structure (C, or relative activation of those synergies), as expressed by the following equation:

$$M^{m \times t} = W^{m \times n} C^{n \times t} + \varepsilon$$

In this equation, W is an m × n matrix where m is the number of muscles (in this study m = 8) and n is the number of muscle synergies. C is an n × t matrix where t is number of time points (101 across the normalized crawling cycle in this study). ε represents the error between the measured EMG data (M) and the reconstructed EMG from W and C. Thus, each column of W represents the relative weighting of muscles in each synergy, and each row of C represents the activation level of each synergy over the gait cycle. Nonnegative matrix factorization was repeated within an iterative optimization, which minimized the sum of squared error between the activations calculated by W × C and the measured EMG data (15). A typical decomposition result is shown in **Figure 2**.

#### **Determining the number of muscle synergies**

We made no a priori assumptions regarding the number of synergies (s) that would be needed to adequately reconstruct the original EMG. The goodness of fit of the data reconstruction was quantified by the variance accounted for (VAF, ranging 0– 1), defined as VAF = 1 − ||ε||<sup>2</sup> /||M||<sup>2</sup> (6, 16, 17). This is a similarity metric that is similar to Pearson correlation coefficient (r 2 ). However, VAF is a more stringent criterion than r <sup>2</sup> because it evaluates both shape and magnitude of the measured and reconstructed curves (17).

For each subject, we determined the least number of muscle synergies that satisfied the following 2 criteria: (1) the overall reconstructed EMGs accounted for at least 90% of the variance across all muscles (VAF>90%); and (2) each reconstructed EMGs accounted for >75% VAF of the measurement from the corresponding single muscle. These criteria are considered conservative to ensure goodness of reconstruction (6). An example of raw EMG signals and the corresponding muscle synergy was shown in **Figure 3**.

#### **Quantifying the structure of muscle synergy**

As an indicator of selective control and coordination, the number of co-activating muscles contributing to a single muscle synergy was calculated. Specifically, muscles in a synergy were defined as active if their normalized weight values exceeded 0.3 (18). Therefore, for each muscle synergy, the number of co-activating muscles varied from 8 (i.e., all the recorded muscles co-activated) to 1 (i.e., no co-activations). The number of co-activating muscles per synergy was calculated for each subject.

### Statistical Analysis

Group difference in the number of muscle synergies, the number of co-activating muscles per synergy, and temporal-spatial parameters were compared using one-way ANOVA (factor of group) with post-hoc Bonferroni corrections for multiple comparisons.

In addition, a 2-way repeated measures ANOVA (withinsubjects factor of directions (AP, ML, and VT), and group) was used for the dependent variable of the normalized JC value. A Bonferroni corrected post-hoc test was used if there was a significant effect.

Spearman rank correlation tests were performed for correlating kinetic indices (the number of co-activating muscles per synergy) and kinematics (normalized JC values). Significance level was set at p < 0.05. All analyses were performed using the statistical software package SPSS18.0. The results that showed a significant effect were marked with an asterisk in all figures.

FIGURE 3 | (A) Example of raw sEMG collected from an infant with typical development (left figure) and the corresponding four muscle synergy identified by NMF algorithm (right figure); (B) Example of raw sEMG collected from an infant with confirmed development delay (left figure) and the corresponding one muscle synergy extracted by NMF algorithm (right figure).

## RESULTS

### Comparison of the Temporal-Spatial Parameters and Normalized JC Values

A one-way ANOVA found no significant effect of group for the temporal-spatial parameters of velocity (F = 0.445, p = 0.643), cadence (F = 0.289, p = 0.750), or stance phase time (F = 0.716, p = 0.493).

The 2-way repeated measures ANOVA showed a significant effect of group (F = 7.591, p < 0.01, observed power = 0.936) and direction (F = 34.301, p < 0.01, observed power = 1) on the normalized JC value. No significant interaction between these 2 factors was found (F = 1.24, p = 0.297). Post-hoc test using Bonferroni corrections revealed higher normalized JC values in the CDD group (averaged across AP, ML, and VT directions) compared to TD (p < 0.01) and ARDD (p < 0.01) groups (shown in **Figure 4A**). Further post-hoc testing showed higher normalized JC values (averaged across TD, ARDD, and CDD groups) in the VT direction compared to AP (p < 0.01) and ML (p < 0.01) direction (shown in **Figure 4B**).

### Comparison of the Number of Muscle Synergies in Infant Crawling

Of the 20 infants with TD measured, two synergies were identified in 60% (12 subjects), three synergies in 35% (7 subjects), and four synergies in 5% (1 subject). Of the 33 infants with ARDD measured, 45.5% (15 subjects) demonstrated two synergies, and 54.5% (18 subjects) three synergies. Of the 13 infants with CDD measured, 15.38% (2 subjects) demonstrated only one synergy and 84.62% (11 subjects) showed two synergies. There was a significant effect of group in the number of muscle synergies (F = 7.194, p = 0.002, observed power = 0.923) during crawling. A significantly reduced number of synergies (1.846 ± 0.375) was identified in infants with CDD during crawling

compared with infants with TD (2.450 ± 0.604) and ARDD (2.450 ± 0.503), respectively. No significant differences were identified in the number of synergies observed between infants with TD and ARDD (**Figure 5**).

## Comparison of the Number of Co-activating Muscles Per Synergy

There was a significant effect of group in the number of coactivating muscles per synergy (F = 4.889, p = 0.011, observed power = 0.786). As shown in **Figure 6**, co-activation levels were significantly higher in the CDD group (5.730 ± 1.129) compared to the TD (4.979 ± 0.501, p = 0.03) and ARDD (4.95 ± 0.787, p = 0.012) group, respectively. There was no significant difference between TD and ARDD group (p > 0.05).

## Correlations of Muscle Synergy and Kinematic Indices

There were no significant correlations found between the number of co-activating muscles per synergy and crawling velocity, crawling cadence or normalized stance phase time.

**Figure 7** reports the number of co-activating muscles per synergy plotted vs. normalized JC value. The number of coactivating muscles per synergy was significantly correlated (r = 0.330, p = 0.007) to the normalized JC value for all infants (**Figure 7**).

### DISCUSSIONS

### Reduced Movement Smoothness During Crawling in Infants With Developmental Delay

During crawling on hands and knees, there was less smoothness in the movements of infants with known developmental delay. A prevailing hypothesis is that the central nervous system (CNS) is organized so that motor output strategy minimizes critical parameters of trajectory such as jerk (19) in order to achieve an accurate and smooth movement. This hypothesis has been supported by observations that those with neurological diseases affecting their CNS, such as stroke (20) and cerebral palsy (21),

have less smoothness than control participants. Extending to the current study, our results of decreased smoothness in infants with developmental delay could be the result of a suboptimal motor command, due to the delayed or impaired development of the CNS (22). On the other hand, the main role of the CNS in generating smooth trajectories has been questioned by other authors (23) who suggested that the intrinsic properties of muscle tissue may be sufficient to produce smooth motion. Considering this hypothesis, the results of decreased smoothness could be related to muscle property differences, such those reported after neurological injury (24). Our results indicate that reduced smoothness of movement in CDD group could also emerge as a result of increased muscle co-activations, which is supported by the significant correlations between the number of co-activating muscles and the normalized JC value (**Figure 7**).

With regards the temporal-spatial parameters, the lack of significant difference between groups was likely because the recruited infants from the 3 groups are at a similar developmental age, as indicated by the clinical assessments.

#### Constrained Neuromuscular Control Strategy During Crawling in Infants With Developmental Delay

The number of muscle synergies identified in infants with CDD was lower compared to infants with TD and ARDD (**Figure 5**). The reduced number of muscle synergies were consistent with the results during walking in individuals with cerebral palsy (8), Parkinson's disease (25), and stroke (6), suggesting impaired muscle coordination. Our results suggest that infants with developmental delay, who were at high risk of cerebral palsy, have less degrees-of-freedom when coordinating muscles and show a constrained neuromuscular control. It is hypothesized that muscle synergies may represent a library of motor subtasks, and the CNS can flexibly combine them to produce complex and natural movements (26). Damage to the CNS, such as in cerebral palsy or stroke, disrupts this combination process, resulting in recruiting less subtasks (synergies). This points to the importance of intact descending control to appropriately recruit a full library of muscle synergies.

Our study also showed that the number of co-activating muscles per synergy was higher in infants with CDD compared to infants with TD and ARDD (**Figure 6**), which implies more muscle co-contractions during crawling in infants with developmental delay. Increased co-contraction of muscles was also shown in individuals with spinal cord injury (7), cerebral palsy and stroke during locomotion (27, 28). Previous studies have shown that a reduction in the descending signals resulted in higher co-activation of muscle (29). Therefore, higher muscles co-activations could indicate that developmental delay in the CDD group was the result of a brain injury in the infant's early life, even if it was not apparent from brain imaging.

#### Clinical Implications

Our results demonstrated that muscle synergy indices (such as the number of synergies and co-activating muscles per synergy) and kinematic output (normalized JC value) were significantly different between infants with confirmed developmental delay (CDD) and typical developing infants (TD), whereas these same variables did not show a significant difference between typical developing infants (TD) and infants at risk of development delay (ARDD). This result validates the concordance between metrics derived in this study and the clinical indicators of motor delay (i.e., Gesell Developmental Scale), suggesting that synergy indices and kinematic output variables are linked in a meaningful way with developmental delay on a group level. However, the group level for these more sensitive metrics can be extended to further understanding of individual participants, especially those at risk for developmental delay later in life.

These results of muscle synergy and movement smoothness analysis during crawling imply a different control strategy between infants with different risks or severity of developmental delay. In spite of their risk factors, the ARDD group was not found to be clinically delayed in the less sensitive clinical measures, but the presence of CNS impairment could become apparent in more sensitive metrics such as muscle synergy assessment, which may be useful in future for the development of more CNS specific rehabilitation plans.

#### Limitations and Opportunities for Future Work

This study quantified and compared the inter-limb muscle synergy and kinematics during crawling between typically developing infants and infants with developmental delay in a novel way. There are a few limitations that need to be acknowledged. Quadrupedal locomotion requires the coordinated behavior of many muscles of the arms, legs, and trunk (30). Because of the small size of infant's limb and the difficulty of measuring locomotion in infants, we measured four primary muscles from each of the arms and legs. In future work, the experimental protocol will be improved by measuring more skeletal muscles such as the gluteus maximus, abdominal and back muscles.

The cohorts of this study were matched on developmental age in order to collect data when they all had similar functional control over their body in order to successfully perform handsand-knees crawling. However, due to motor delays there was a significant difference in the chronological age of the infant groups studied where the CDD group was older. If the hypotheses being tested were relative to the chronological or biological age of the CNS, the TD group could have been matched on chronologic age—but it would be anticipated that in that case an older cohort of TD children would have at least similar, if not better, skill in crawling than the current typical group, and thus may result even bigger difference between TD and CDD groups.

We also recognize that there are other potential confounding factors, such as different risk factors and etiologies for development delay (premature delivery, low birth weight, or lack of oxygen during birth), which may indicate different mechanisms for development delay. However, the investigation of these confounding factors is beyond the interests of the current study.

There are the two extreme data points in **Figure 7**, showing 2 CDD participants who demonstrated both very high JC values and co-activation of all or nearly all of the recorded muscles during crawling, which was visually different than the cluster formed by the TD group. One of them was identified by SPSS boxplot (use a step of 1.5 × interquartile range) as an "out" value. If we exclude this data point and the correlation within the rest 65 infants was still significant (r = 0.298, p = 0.016).). These two extreme data points suggests that they always activate all the muscles and lack the movement smoothness. Some of the infants with confirmed development delay (CDD) will likely receive a diagnosis of cerebral palsy, which is characterized with by the presence of spasticity and decreased selective motor control (31, 32). A future study should include follow-up with the CDD and ARDD infants, in order to ascertain if features of the crawling data are predictive of a later diagnosis of cerebral palsy.

Even considering the limitations above, the potential for assessing motor function and understanding of the state of the neuromuscular system during crawling period is an exciting prospect. Assessment of pathological impairment in motor control during walking can be conducted by gait analysis, which is been widely used in clinics and typically provides quantified

#### REFERENCES


metrics of kinematics and muscle activity (33). For those infants without walking ability, movement abnormalities are typically assessed by screening tests or visual analysis of their movement quality. This study demonstrates that utilizing more quantitative metrics can reveal impaired neuromuscular strategies before the onset of walking skills and provide insight for development of rehabilitation of protocols during infants' crawling stage. The long-term goal of this work is to develop a standardized measure, similar to gait analysis, that can assess motor function in infant crawling on an individual basis.

### CONCLUSION

This study demonstrated that infants with developmental delay demonstrated fewer inter-limb muscle synergies, increased number of muscles that co-contracted, and reduced movement smoothness during crawling on hands and knees, compared to typical developing infants. In addition, more co-activations across inter-limb muscles are considered to be attributable to the reduced movement smoothness in infant crawling. These muscle coordination and kinematic output deficits revealed impaired neuromuscular strategies during the infant crawling stage.

### AUTHOR CONTRIBUTIONS

WH, XW, NX, LC, and XZ designed the work. QX and YL collected the data. QX analyzed the data. JY, TS-M, and WH interpreted the data. QX and JY drafted the manuscript. TS-M and WH helped to create the final report.

#### FUNDING

This work was supported by the National Natural Science Foundation of China (31470953, 31771069), the Chongqing Science & Technology Program (cstc2016shmszx130060, cstc2015jcyjB0538), and the graduate research and innovation foundation of Chongqing, China (Grant No. CYB17038).

### ACKNOWLEDGMENTS

We much appreciate the infants and their parents for participating in this study. Thank you to the Department of Rehabilitation Center, Children's Hospital of Chongqing Medical University for recruiting infants and their help in data collection.

quadrupeds. J Neurophysiol. (2009) 101:603–13. doi: 10.1152/jn.91125 .2008


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Xiong, Wu, Yao, Sukal-Moulton, Xiao, Chen, Zheng, Liu and Hou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Increased Corticomuscular Coherence and Brain Activation Immediately After Short-Term Neuromuscular Electrical Stimulation

Rui Xu1,2, Yaoyao Wang1,2, Kun Wang1,2, Shufeng Zhang1,2, Chuan He1,2 and Dong Ming1,2 \*

*<sup>1</sup> Lab of Neural Engineering & Rehabilitation, Department of Biomedical Engineering, College of Precision Instruments and Optoelectronics Engineering, Tianjin University, Tianjin, China, <sup>2</sup> Tianjin International Joint Research Center for Neural Engineering, Academy of Medical Engineering and Translational Medicine, Tianjin University, Tianjin, China*

Neuromuscular Electrical Stimulation (NMES) is commonly used in motor rehabilitation for stroke patients. It has been verified that NMES can improve muscle strength and activate the brain, but the studies on how NMES affects the corticomuscular connection are limited. Some studies found an increased corticomuscular coherence (CMC) after a long-term NMES. However, it is still unknown about CMC during NMES, as relatively pure EMG is very difficult to obtain with the contamination of NMES current pulses. In order to approach the condition during NMES, we designed an experiment with short-term NMES and immediately captured data within 100 s. The repetition of wrist flexion was used to realize static muscle contractions for CMC calculation and dynamic contractions for event-related desynchronization (ERD). The result of 13 healthy participants showed that maximal values (*p* = 0.0020) and areas (*p* = 0.0098) of CMC and beta ERD were significantly increased immediately after NMES. It was concluded that a short-term NMES can still reinforce corticomuscular functional connection and brain activation related to motor task. This study verified the immediate strengthen of corticomuscular changes after NMES, which was expected to be the basis of long-term neural plasticity induced by NMES.

Keywords: neuromuscular electrical stimulation, corticomuscular coherence, event-related desynchronization, functional connection, brain activation

### INTRODUCTION

Neuromuscular Electrical Stimulation (NMES) is a technique that can generate contractions of paralyzed or paretic muscles by applying electrical current on these muscles (1). Confidential evidence has shown that NMES can increase the maximal voluntary contraction and neural activation assessed by the twitch interpolation technique (2). Poststroke rehabilitation with NMES has been found to effectively prevent muscle atrophy, improve muscle strength (1, 3) and coordination (4). More recently, a study published in Nature Communications revealed the efficacy and mechanisms of brain-actuated functional electrical stimulation via the clinical performance and functional connectivity (5).

The influence of NMES on muscles is easy to understand as NMES is directly applied to muscles. However, the effect in muscles is not enough to realize motor rehabilitation, since the brain plays an important role in motor recovery. Phenomena of event-related desynchronization/synchronization

#### Edited by:

*Ping Zhou, University of Texas Health Science Center at Houston, United States*

#### Reviewed by:

*Qining Wang, Peking University, China Xiang Chen, University of Science and Technology of China, China*

> \*Correspondence: *Dong Ming richardming@tju.edu.cn*

#### Specialty section:

*This article was submitted to Stroke, a section of the journal Frontiers in Neurology*

Received: *30 June 2018* Accepted: *01 October 2018* Published: *23 October 2018*

#### Citation:

*Xu R, Wang Y, Wang K, Zhang S, He C and Ming D (2018) Increased Corticomuscular Coherence and Brain Activation Immediately After Short-Term Neuromuscular Electrical Stimulation. Front. Neurol. 9:886. doi: 10.3389/fneur.2018.00886*

**42**

(ERD/ERS) of EEG could be found at the frontal and parietal areas when limb movements are executed or imagined, which shows a power decrease/increase in the alpha (8–13 Hz) and beta (14–30 Hz) bands (6). The ERD pattern was used to indicate brain activation and sensitive to different movement speed on action observation (7). What's interesting is that NMES applied on muscles also affects Electroencephalogram (EEG) oscillatory (8), which verifies that NMES on the muscles can activate related brain area, and this activation pattern represented by ERD is similar to that under active movement. Lo et al. used near infrared spectroscopy (NIRS) to investigate cortical activation of different-intensity electrical stimulations (9). These studies evaluated the efficacy of NMES from the view of brain activation.

Moreover, corticomuscular coherence (CMC) is a method to estimate neural coupling via Magnetoencephalogram (MEG) or EEG and Electromyogram (EMG). CMC has drawn much attention since it was first discovered by Conway et al. (10). For now, we have known that the strength of CMC is adjusted or affected by attention (11), muscle contraction type (12, 13), muscle contraction force (14, 15), muscle fatigue (16, 17), and motor learning (18). As CMC statistically calculated the synchronization between the brain and muscle signals, it reflects functional connection between the motor cortex and muscles (19). Due to this, CMC of stroke patients has obtained some focus since Mima et al. first revealed that there was significant EEG-EMG coherence only in the unaffected side of the brain (20). Except for the amplitude of CMC, the location is still different for stroke patients and the control. Rossiter et al. found that the CMC of stroke patients were located more widely than healthy people (21). It may verify that brain regions in the contralesional hemisphere were involved to help recover motor functions. In 2017, the result of an interesting study demonstrated that although CMC was reduced in the acute phase after stroke, there was no significant change within the following 4 ∼ 6 weeks despite of improved behavioral performance (22). Maybe CMC is not an efficient marker for early recovery of motor function following stroke. The continuous learning of CMC should help us make CMC more sensitive to the rehabilitation of stroke.

CMC calculation provides a new perspective to study the efficacy of NMES. Lai et al. have done interesting and important exploration on the EEG-EMG coherence affected by long-term sensory electrical stimulation (23, 24). They found that the electrical stimulation causing no muscle contraction and pain increased the EEG-EMG coherence. The accurate CMC during NMES is also necessary as it provides direct information on the effect of NMES, and reflects transient neural plasticity. However, it is difficult to obtain pure EMG, as the stimulation current contaminates EMG severely. Therefore, we designed an experiment to capture the immediate effect of a shortterm NMES and analyzed both functional connection and brain activation via CMC and ERD respectively. We hypothesized that CMC could be strengthened immediately after NMES.

### MATERIALS AND METHODS

#### Participants and Experiments

Thirteen healthy right-handed people (5 females and 8 males; mean age: 21.2 ± 1.1 years old) from Tianjin University participated in the study. The participants had no history of neuromuscular disorders. The study was approved by the ethics committee of Tianjin University. All participants signed informed consent in advance.

The experiment consisted of one long voluntary session (300 s) and three stimulated plus short voluntary sessions (100 s + 100 s) shown in **Figure 1A**. There was a rest for 5 to 10 min between two sessions. There were 30 trials in the long voluntary session, and 10 trials in each short voluntary session. Each voluntary trial started with 2-s wrist flexing, followed by 5-s wrist flexion holding and 1-s relaxing, and ended up with 3-s resting (**Figure 1B**).

Before the experiment, the participant was seated in front of a 17" monitor with his right arm on the table (**Figure 2A**). His right hand was relaxed to make a slight fist. During the experiment, the participant followed the instructions on the monitor (generated by Psychtoolbox within Matlab) to complete each trial: he flexed his right wrist when "Flexing" showed up in the monitor, held the flexed wrist for the "Holding" part, and then relaxed and rested according to the cue (**Figure 1B**). There was a time label at the onset of holding part. During the stimulated session, the participant was seated still like in the voluntary session and his/her wrist was relaxed without any voluntary movement. His/Her right flexor carpi radialis (FCR) was electrically stimulated for 100 s, with the stimulation frequency at 30 Hz and the peak current varying from 7 to 13 mA for different participants (mean: 10.9 ± 2.3 mA). The wrist of the participant was kept flexed under this stimulation. The stimulation intensity was determined at each participant's tolerance with an actual wrist flexion before the first stimulated session: the peak current was raised from 5 mA by 1 mA each time until the participant felt uncozy and asked to stop (the maximal current), and the peak current used in the stimulated session was 1 mA less than the maximal current. The maximal current is listed individually in **Table 1**.

### Data Recording and Preprocessing

EEG and surface electromyography (sEMG) data were acquired simultaneously with a Neuroscan SynAmps2 amplifier, hardware-filtered in the frequency range of 0.015–250 Hz and sampled at 1,000 Hz. EEG data was recorded with 64 electrodes located in the positions following the 10/20 system (**Figure 2B**), while sEMG data was recorded by 2 Ag/AgCl electrodes placed on the surface of FCR (2-cm interelectrode distance). The recorded data were referenced to the nose and grounded at the prefrontal lobe. An additional 50-Hz notch filter was used during data acquisition.

Data analyses were performed using Matlab R2017b (MathWorks, MA, USA), with the toolbox EEGLAB (Swartz Center for Computational Neuroscience; http://sccn.ucsd.edu/ eeglab/). The acquired EEG data at C1, C3, and C5 electrodes were re-referenced using the surface Laplacian technique (25) according to (1), (2), and (3).


FIGURE 1 | Experimental protocol. (A) The complete experiment, consisting of one long voluntary sessions (300 s) and three stimulated plus short voluntary sessions (100 s + 100 s). (B) One trial in voluntary sessions, with 1-s wrist flexing, 5-s flexion holding, 1-s wrist relaxing and 3-s wrist resting.

where V<sup>I</sup> (I = C1, FC1, CP1, Cz, C3, FC3, CP3, C5, FC5, CP5, or T7) indicates the EEG data acquired at the electrode I and VC<sup>1</sup> ′ , VC<sup>3</sup> ′ , and VC<sup>5</sup> ′ were the re-referenced EEG data at C1, C3, and C5. Then a 4th-order zero-phase Butterworth filter was used to obtain filtered EEG data (5 ∼ 45 Hz) and sEMG data (20∼250 Hz). The full-wave rectified sEMG were obtained as the absolute value of the data.

There were 30 trials (in the long voluntary session) before and 10 × 3 trials (in the short voluntary sessions) after the stimulated sessions (**Figure 1A**). A 3072-point data part started from the time label was extracted from the "Holding" part of each trial (**Figure 1B**). The re-referenced C1 EEG, C3 EEG, C5 EEG, and sEMG data before or after stimulation consisted of 30 data parts. Each data part was further divided into 6 segments of 512 points. In total, 180 data segments were used to calculate the EEG-EMG coherence.

### CMC

Denoting the fast Fourier transform (FFT) of the ith segment of C3 EEG by Xi(f) and of the ith segment of rectified sEMG by Yi(f), the coherence (Coh) at frequency f was estimated as:

TABLE 1 | The peak current and C3 EEG-EMG coherence for individuals.


*The values in bold indicate increases after NMES.*

$$\text{Co}h\left(f\right) = \frac{\left|\sum\_{i=1}^{N} X\_i^\*\left(f\right) Y\_i\left(f\right)\right|^2}{\sum\_{i=1}^{N} X\_i\left(f\right) X\_i^\*\left(f\right) \sum\_{i=1}^{N} Y\_i\left(f\right) Y\_i^\*\left(f\right)}\tag{4}$$

where i = 1, . . . , N is the number of data segments available for analysis, and <sup>∗</sup> denotes complex conjugate. The use of 512-point segments with a sampling rate of 1,000 Hz provided a 1.95 Hz frequency resolution in the coherence spectra. The C1 and C5 EEG-EMG coherence was calculated in the same way.

The confidence level for the coherence (26) was calculated as:

$$CL\left(\alpha\right) = 1 - (1 - \alpha)^{\frac{1}{N-1}}\tag{5}$$

where N is the number of data segments and α is the desired level of confidence. We considered coherence to be significant above the 95% confidence limit (α = 0.95). As there were 180 segments for coherence calculation, CL was 0.0166 according to (5).

The significant coherence Coh<sup>s</sup> used in this study was calculated as (6). This calculation neglected the small differences of CMC below CL. The maximal value of Coh<sup>s</sup> was also used to indicate the strongest corticomuscular connection before or after NMES.

$$\text{Coh}\_{s}\left(f\right) = \begin{cases} \text{Coh}\left(f\right) - 0.0166 \text{ if } \text{Coh}\left(f\right) > 0.0166\\ 0 & \text{else} \end{cases} \tag{6}$$

The mean curve of significant coherence Mcoh was obtained by (7).

$$M\_{\rm Colh} \left( f \right) = \frac{1}{K} \sum\_{i=1}^{K} Coh\_s \left( f \right) \tag{7}$$

where K is the total number of participants.

#### Area of Significant Coherence

The C1, C3, and C5 EEG-EMG coherence values below CL were set to zero according to (6). Only the significant coherence was used in the area calculation. The area of significant coherence (ACoh) within 5 ∼ 45Hz can be used to estimate the strength of corticomuscular coupling (18), and it was calculated as:

$$A\_{Coh} = \sum\_{f=5Hz}^{45Hz} Coh \,\left(f\right) \tag{8}$$

#### Center of Gravity for Frequency

To detect the frequency shifts of the coherence spectrum, we calculated the Center of Gravity for the frequency (CoG<sup>f</sup> ), that is, the frequency at which coherence is concentrated and balanced. The CoG<sup>f</sup> of C3 EEG-EMG coherence was obtained by (9).

$$CoG\_f = \frac{\sum\_{i=1}^{n} f\_i \cdot Coh\_s \text{ ( $f\_i$ )}}{\sum\_{i=1}^{n} Coh\_s \text{ ( $f\_i$ )}} \tag{9}$$

where i=1,. . . ,n indicates the number of significant bins with its respective frequency f<sup>i</sup> and coherence Coh<sup>s</sup> .

#### Median Frequency of sEMG

In order to exclude the effect of muscle fatigue on CMC, we also calculated the median frequency of sEMG before and after the stimulated sessions. The median frequency is defined as the frequency that divided the spectrum into two equal areas. It has been widely used in the studies related to muscle fatigue (27, 28). The median frequency of sEMG during different sessions were calculated and compared to indicate the fatigue states in this study.

#### ERD

The re-referenced and filtered EEG data at C3 channel was downsampled to 200 Hz for ERD analysis. The event-related spectral perturbation (ERSP) method allowed us to inspect the spectral power changes of EEG in the view of time-frequency domain. Therefore, ERSP was calculated as:

$$ERSP\left(f,t\right) = \frac{1}{n} \sum\_{i=1}^{n} \left(F\_i(f,t)^2\right) \tag{10}$$

where i = 1, . . . , n is the number of trials, and F<sup>i</sup> f , t is the spectral estimation of the ith trial at frequency f and time t (29). Short-time Fourier transform (STFT) was used to perform time-frequency analysis with a Hanning window. The number of windows was set to 200 with the length of 512 points. ERSP was calculated using 30 trials of data and the data length was 10 s for each trial, with 2 s before movement onset and 8 s after (1-s flexing, 5-s holding, 1-s relaxing and 1-s resting). Baselinenormalized ERSP was calculated relative to the baseline period (before movement onset).

Frontiers in Neurology | www.frontiersin.org

In order to investigate the difference of brain oscillation before and after NMES, the ERD at C3 within 1 s after the movement onset was extracted as follows:

$$ERD\left(f\right) = \sum\_{t=0s}^{1s} ERSP\left(f, t\right) \tag{11}$$

The first second after movement onset indicated the wrist flexion period, excluding the holding part.

#### Statistical Analysis

All the features mentioned above before and after NMES, including the maximal values and areas of Coh<sup>s</sup> , CoG<sup>f</sup> , and ERD values, were compared with Wilcoxon signed rank test. The significance was calculated two-tailed. For the CoG<sup>f</sup> , only data with significant coherence both before and after NMES were considered.

### RESULTS

## CMC

We calculated the EEG-EMG coherence of all the subjects. The C3 EEG-EMG coherence before (blue line) and after (red line) NMES of each participant is listed in **Figure 3**. There were some participants who did not present significant EEG-EMG coherence before or after NMES, such as P9, P11, P12, and P13.

The grand average of significant coherence is shown in **Figure 4**. The peak values of mean coherence after NMES were larger than those before NMES in **Figures 4A–C**. It was obvious that NMES had different influence on these three channels of coherence. The maximal values and areas of significant coherence (listed in **Table 1**) were calculated for further statistical analysis. The Wilcoxon signed rank test was used and the result indicated in **Table 2** that the maximal value and area of significant coherence for C3 EEG-EMG coherence after NMES was significantly larger than those before NMES (Max: p = 0.0020; Area: p = 0.0098). Although areas of C1 and C5 EEG-EMG coherence after NMES also increased, there was no significant difference.

There were only 9 participants who showed significant C3 EEG-EMG coherence both before and after NMES. The CoG<sup>f</sup> of these 9 participants, its average and median values were listed in **Figure 5**. The average frequency was increased after NMES (Avg.: from 23.7 to 27.8 Hz), but there was no significant frequency shift after NMES according to the result we obtained.

The median frequencies of sEMG are shown in **Figure 6**. The one-way repeated measure ANOVA was used to compare these median frequencies, and no significant difference indicated that the fatigue state of the muscle remained the same given the sensitivity of the median frequency.

#### ERD

**Figure 7** shows the averaged ERSP of C3 EEG before and after NMES. There were obvious ERD patterns (blue area in the figure) in both mu (8∼13 Hz) and beta (14∼30 Hz) rhythms at the beginning and end of the movement part. The average ERD of certain areas are listed in **Table 3**. It shows that the ERD patterns seems to be weakened in the "holding" part between 1 and 6 s and the strongest ERD patterns occur mainly in mu rhythms. However, these changes are not significant according to the result of a two-way (time: Flexing, Holding, and Relaxing; frequency: mu and beta rhythms) repeated measure ANOVA.

In order to investigate and compare the brain activation before and after NMES, we calculated the ERD values of C3 EEG at different frequency bins according to (5). The Wilcoxon signed rank test was used, and the significant differences between two conditions were shaded by gray blocks in **Figure 8**. The blue and red lines represented normalized ERSP before and after NMES respectively. The ERD patterns were significantly stronger in three sub-beta frequency bands.

### DISCUSSION

This study for the first time compared the coticomuscular coherence before and immediately after short-term motoriallevel NMES. We designed an experiment especially for exploring the effect of NMES on the functional connectivity between the brain and muscles. It was clear that the coticomuscular coherence during active movement after NMES was significantly stronger than that before NMES. The result illustrated that NMES strengthened the interaction between the brain and muscles. We also calculated the ERD patterns before and after NMES. The analysis indicated that the ERD patterns were strengthened after NMES. It seems that NMES has a positive influence on the interaction between the brain and muscles and the activation of the brain.

There are many studies working on the rehabilitation effect of NMES. Most of them focus on the comparison of features after NMES training. For example, Sota et al. compared some gait parameters, such as the time of 10-m walking and range of motion for ankle joint, pre- and postintervention to investigate the characteristics of NMES responders (30). CMC combined with clinical functional test was used to estimate the effect of sensory NMES with motor training (24). Although the effect remained after NMES is very important, the instantaneous body response during NMES is also a key point of NMES studies. However, as the stimulation pulses affected and contaminated EMG severely, the studies on the effect during NMES are limited to comparing the features free of EMG, such as walking speed (31), ERD (8) and steady-state somatosensory evoked potential (SSSEP) (32) of the brain, and muscle thickness (33). These studies analyzed the effect of NMES on the whole body, or the brain and muscles separately, without considering motor control based on the interaction between the brain and muscles.

The mechanism of motor control can be revealed by CMC. The application of NMES was certain to cause extreme contamination of sEMG, so we had to compare the CMC before and after NMES to guarantee the data quality and result reliability. There was no studies on the effect-remaining time of NMES, but we believed that the effect of NMES should be the strongest immediately after NMES except for the effect during

NMES. In this case, the stimulated session and the voluntary session after stimulation was divided into three equal parts individually. We tried to use this paradigm to approach the condition during NMES.

line indicates significant coherence before NMES, while the red line indicates significant coherence after NMES.

Although CMC has been widely studied and used (34, 35), its generation is still under debate. From the perspective of coherence, a significant coherence between two subsystems can be achieved by either one-way information flow, reciprocal communication, or the third rhythm generator affecting both (36). The result of our study verified that there were at least two directions of information flow: one was from the brain to muscles, sending cortical motor command; and the other one was from muscles to the brain, caused by NMES. As SSSEP was observed during NMES (32), the regulation of brain activities by NMES was determined. It was possible that the significantly increased CMC was the residual effect of SSSEP.

An interesting detail found in this study was that NMES increased C1, C3, and C5 EEG-EMG coherence average, but only the change in C3 EEG-EMG coherence was significant. It was deduced that the variation of coherence caused by NMES could be used to locate the cortex area in charge of the executed movement. However, whether this change varied with

TABLE 2 | Maximal value and area of significant coherence.


*The values in bold are with significant difference compared to the values before NMES.* \**p* < *0.05.*

the location of NMES or the contracted muscles was unclear in the present study.

Higher CMC often indicated better communication between the brain and muscles, and higher beta band CMC indicate good motor performance (37). Moreover, beta-band CMC was deemed to be related to motor tasks and performance (12, 14, 15, 38). Our main result based on **Table 2** was that NMES increased

C3 EEG-EMG coherence, which was consistent with the newly published work of Pan et al. (24). They reported an increase of CMC after 4-weeks sensory electrical stimulation. Neural plasticity was believed to contribute to CMC increase of stroke patients, and it was crucial for learning new motor tasks (39). For healthy participants in our study, they did not learn new movement, but learn new muscle contraction patterns (NMES). This should also be regarded as learning a new motor task. According to Hebbian and homeostatic plasticity (40, 41), the CMC change after a short-term NMES reflected a transient plasticity and it could go back to former state without repetition of stimulation.

CMC could be affected by many factors. In this study, muscle contraction type and muscle force were not considered because the tasks in voluntary sessions were the same. The wrist flexed according to the clues on the screen and the data processed were extracted during the "Holding" time. The muscle contraction type in this study was static contraction, and the muscle forces in all the voluntary sessions, which was to keep the wrist flexed, should be the same. No precision requirement guaranteed that there was no difference in the attention (11). The median frequency ofsEMG was analyzed to indicate no significant fatigue states. Therefore, the muscle contraction type, muscle force, attention and muscle fatigue were excluded.

The significant CMC is not a universal phenomenon for every person(18, 24). In our study, there were 4 participants who did not show significant CMC before or after NMES. Their peak currents of NMES (10, 15, 10, and 13mA) were relatively

#### TABLE 3 | Average of ERSP.


larger (the median value of these currents is 10mA). A larger stimulation current meant this participant was less sensitive to the stimulation, and he needed stronger stimulation to generate muscle contraction. Therefore, we deduced that the participant with insignificant CMC was most likely insensitive to NMES. However, this should be further verified with specifically designed experiments.

ERD patterns was often used to indicate the brain activation. In **Figure 7**, the weakening of ERD occurred for the holding part. This phenomenon was also shown in (42). This implied that the maintenance of the current sensorimotor state was related to the ERD rebound. The reasons may be that holding a posture was easier to execute than dynamic motor tasks, and the brain completed the static motor task in a low activation level. In order to obtain an obvious ERD variation, we compared ERD of data within the first second after the movement onset. It was found that NMES could induce a stronger ERD pattern. Vidaurre et al. found a stronger ERD pattern during NMES than motor imagery (MI), and successfully used NMES-induced patterns to decode MI (43). In our study, ERD was also strengthened in beta band after NMES. The beta ERD was linked more closely to the primary motor cortex (44, 45). Therefore, the significant cortical beta rhythm suppression showed brain activation related to motor control.

NMES can increase the excitability of human corticospinal (CS) pathways to muscles, which is usually estimated by the motor evoked potentials (MEPs). Mang et al. compared the MEPs induced by transcranial magnetic stimulation (TMS) before and after an NMES session and found that the MEP amplitude after NMES was significantly larger than that before NMES (46, 47). Whether the increase in cortical excitability is due to changes at the spinal level, cortical reorganization, or both is unclear. Such increases can strengthen CS pathways damaged by injury or disease and result in enduring improvements in function (1, 48). Here, we hypothesized that the higher CMC and stronger ERD were caused by the strengthened CS pathways.

Our study did not consider the effect caused by stimulation intensities. A study of healthy participants found that higher NMES current intensities led to greater sensorimotor network activation, and this may be attributable to increased attentional/pain processing and to increased sensorimotor integration (49). Therefore, a maximal tolerated intensity was used in our study in order to obtain significant changes of CMC before and after NMES. However, it was still unclear that how the stimulation intensity influenced CMC or whether there was a difference in CMC for sensory- and motorial-level NMES. The comparison will help understand the function of the sensorimotor circuit. The experiment was designed to approach the condition during NMES, but how CMC changed during NMES still needed to be verified.

The study was undertaken among healthy participants. Therefore, whether NMES would have the same effect in stroke patients needs to be studied further. However, NMES has been used to strengthen CS pathways in stroke rehabilitation (1, 48). It is hypothesized that the strengthened CS pathways will induce a stronger CMC in stroke patients. To be noted, the EEG channels of stroke patients for CMC calculation should be different from that of the healthy, as the lesion may be located at C3 or C4. In order to explore the effect of NMES on patients' CMC, some stroke patients will be recruited for the future work.

## CONCLUSION

As the estimation of NMES real-time efficacy was limited to brain or bodies separately, we designed an experiment with NMES and repeated voluntary wrist flexion, to explore the instantaneous effect after NMES. The result showed a significant increase of EEG-EMG coherence caused by NMES. Additionally, the significant increment was located in C3 position. The strengthened beta ERD indicated stronger brain activation related to motor function after NMES. Therefore, NMES not only strengthened brain activation, but it also induced a stronger connection between the brain and muscles. This result will help understand NMES-induced corticomuscular connection, and predict the body change during NMES. Based on transient neural plasticity, the immediate change after NMES lays a basis of long-term neural rehabilitation.

### AUTHOR CONTRIBUTIONS

All authors contributed to conception and design of the study and were involved in drafting and critically revising the

### REFERENCES


manuscript. Additionally, RX carried out data processing and feature comparison, interpreted the results and prepared the first draft paper. YW completed EEG and EMG pre-processing. KW helped calculate ERD patterns. YW, SZ, and CH carried out the experiment. DM provided interpretation of the results and worked up the draft paper into the final version. All authors gave final approval for publication.

### FUNDING

This work was supported by National Key Research and Development Program of China under grant (No. 2017YFB1300300), National Natural Science Foundation of China (No.81630051, 91520205, 61603269, 81601565, 81571762, and 31500865), and Tianjin Key Technology R&D Program (No. 15ZCZDSY00930).

#### ACKNOWLEDGMENTS

We would like to thank all the participants and Dr. Bin Gu for his help on the experimental platform.


subacute phase after stroke. Clin Neurophysiol. (2017) 128:2217–26. doi: 10.1016/j.clinph.2017.08.033


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Xu, Wang, Wang, Zhang, He and Ming. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Position as Well as Velocity Dependence of Spasticity—Four-Dimensional Characterizations of Catch Angle

Yi-Ning Wu<sup>1</sup> , Hyung-Soon Park <sup>2</sup> , Jia-Jin Chen<sup>3</sup> , Yupeng Ren<sup>4</sup> , Elliot J. Roth<sup>5</sup> and Li-Qun Zhang6,7 \*

 *Department of Physical Therapy and Kinesiology, University of Massachusetts Lowell, Lowell, MA, United States, Department of Mechanical Engineering, Korean Advanced Institute of Science and Technology, Daejeon, South Korea, Department of Biomedical Engineering, National Cheng Kung University, Tainan, Taiwan, <sup>4</sup> Rehabtek LLC, Glenview, IL, United States, <sup>5</sup> Department of Physical Medicine and Rehabilitation, Northwestern University, Chicago, IL, United States, Department of Physical Therapy and Rehabilitation Science, Department of Orthopaedics, University of Maryland, Baltimore, MD, United States, <sup>7</sup> Department of Bioengineering, University of Maryland, College Park, MD, United States*

#### Edited by:

*Ping Zhou, University of Texas Health Science Center at Houston, United States*

#### Reviewed by:

*Lining Zhang, PLA General Hospital, China Amanda Fang Lin, Rosalind Franklin University of Medicine and Science, United States*

> \*Correspondence: *Li-Qun Zhang l-zhang@som.umaryland.edu*

#### Specialty section:

*This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology*

Received: *28 June 2018* Accepted: *25 September 2018* Published: *26 October 2018*

#### Citation:

*Wu Y-N, Park H-S, Chen J-J, Ren Y, Roth EJ and Zhang L-Q (2018) Position as Well as Velocity Dependence of Spasticity—Four-Dimensional Characterizations of Catch Angle. Front. Neurol. 9:863. doi: 10.3389/fneur.2018.00863* We investigated the muscle alterations related to spasticity in stroke quantitatively using a portable manual spasticity evaluator.

Methods: Quantitative neuro-mechanical evaluations under controlled passive elbow stretches in stroke survivors and healthy controls were performed in a research laboratory of a rehabilitation hospital. Twelve stroke survivors and nine healthy controls participated in the study. Spasticity and catch angle were evaluated at 90◦ /s and 270◦ /s with the velocities controlled through real-time audiovisual feedback. The elbow range of motion (ROM), stiffness, and energy loss were determined at a slow velocity of 30◦ /s. Fourdimensional measures including joint position, torque, velocity and torque change rate were analyzed jointly to determine the catch angle.

Results: The catch angle was dependent on the stretch velocity and occurred significantly later with increasing velocity (*p* < 0.001), indicating position dependence of spasticity. The higher resistance felt by the examiner at the higher velocity was also due to more extreme joint position (joint angle) since the spastic joint was moved significantly further to a stiffer elbow position with the higher velocity. Stroke survivors showed smaller ROM (*p* < 0.001), higher stiffness (*p* < 0.001), and larger energy loss (*p* = 0.005). Compared to the controls, stroke survivors showed increased reflex excitability with higher reflex-mediated torque (*p* < 0.001) and at higher velocities (*p* = 0.02).

Conclusion: Velocity dependence of spasticity is partially due to joint angle position dependence with the joint moved further (to a stiffer position where higher resistance was felt) at a higher velocity. The "4-dimensional characterization" including the joint angle, velocity, torque, and torque change rate provides a systematic tool to characterize catch angle and spasticity quantitatively.

Keywords: quantitative evaluation, spasticity, muscle, stroke, catch angle

## INTRODUCTION

Spasticity commonly occurs to patients with neurological disorders, such as stroke, spinal cord injuries, cerebral palsy, and multiple sclerosis (1–3). Spasticity is commonly defined as "a motor disorder characterized by a velocity-dependent increase in tonic stretch reflexes (muscle tone) with exaggerated tendon jerks, resulting from hyperexcitability of the stretch reflex, as one component of the upper motor neuron syndrome" (4). Various measures have been used to assess muscle alterations associated with spasticity. In the clinical setting, spastic muscle is usually evaluated by grading the resistance to a passive stretch felt by a clinician using the Ashworth scale, modified Ashworth scale (MAS), and the Tardieu scale (5, 6). The felt resistance could be caused by a combination of neural and peripheral origins (i.e., biomechanical factors such as soft tissues or muscle properties). Although clinical measures are convenient to use, they can be subjective, less sensitive and qualitative rather than quantitative to varying degrees. Previous studies raised questions about reliability of the MAS assessment of spasticity (7–10). On the other hand, the Tardieu scale has been suggested as an alternative to the MAS (6). Tardieu scale is conducted using various stretch velocities rather than using only one velocity in MAS while determining the angle where resistance felt (i.e., catch angle). It is argued that the MAS does not differentiate between spasticity and contracture, while the Tardieu scale is not confounded by the presence of contracture (11). However, with either scale, determinations of the catch angle (12, 13) and range of motion (ROM) are influenced by stretch velocities and stretch force and subject to errors in reading the joint angle during assessments.

A quantitative assessment with controlled passive stretches is needed to improve the reliability of the clinical measures. Well-controlled quantitative measures, based on motorized mechanical perturbations and electrophysiological approaches, are mostly used in laboratory settings (5, 14–16), but size and ease-of-use issues limit their applications in clinical settings (17–19). Several portable devices have also been developed and spasticity evaluations were performed by deriving viscous neuro-mechanical properties of the limb from passive movement kinematics and joint reaction torques (19–22). However, those measurements did not translate easily to the common clinical assessments of ROM and catch angle. Reflex threshold measured in joint angle during passive movement has been used effectively to evaluate spasticity by investigating the onset of muscle activation to applied disturbance (23–25).

In spite of the current development of spasticity quantification, in clinical setting thus far, clinicians evaluate spasticity based on how much resistance they feel as well as where they feel the reactive resistance while manipulating the joint quickly. The relation with regard to velocity dependence between stretch-induced muscle activation onset and the resistance (catch) felt by the clinicians in the stroke survivors has not been investigated thoroughly. Furthermore, it is not clear whether catch angle is also joint angle position dependent. In other words, whether joint angle might play a role in the increasing resistance felt by clinicians at a higher velocity and being judged as velocity-dependent spasticity is uncertain. A comprehensive but simple way considering stretch velocities, reflex-mediated muscle torque and joint angle is needed to assist clinicians understand the muscle alteration due to neurological disorders and interventions. Therefore, the purpose of this study was to introduce an innovative and quantitative way to depict the spasticity according to the concept of Tardieu scale and further examined the contribution of joint angle position dependence to the catch felt by the examiner.

### MATERIALS AND METHODS

#### Subjects

Twelve chronic stroke survivors (53.0 ± 8.5 years old, ten males and two females) who had a stroke more than 1 year (9.3 ± 5.6 years) and nine healthy controls (51.4 ± 24.9 years old, nine males) were included in this study. The stroke survivors with elbow flexors spasticity were recruited in the study. The subjects who had shoulder or elbow contractures were excluded from the participation. **Table 1** depicts the characteristics of the stroke survivors. The healthy controls had no history of neurophysiologic or musculoskeletal disorders. The study was approved by the Institutional Review Board of Northwestern University. All subjects gave written informed consent before the experiment.

#### Instrumentation

The manual spasticity evaluator (MSE) used in this study was set up as a portable device to assess spastic muscles. A torque sensor (Transducer Techniques, CA, USA) and a hollow-shaft potentiometer (Vert-X51, Contelec AG, Switzerland) comprising MSE, were used to measure joint torques and joint positions respectively (26, 27). Adjustable braces and supports were used to position the forearm and upper arm properly with respect to the MSE (**Figure 1**). Two mechanical stops were used to restrict the moving range of device that prevents over-stretching of the spastic joint. Biceps and triceps muscles activations were monitored using surface EMG electrodes (Bagnoli-8, Delsys Inc., Boston, USA). The torque, position and EMG signals were


the MSE rotation axis. A torque sensor and a hollow-shaft potentiometer were used to measure the joint torque and joint position, respectively. Surface EMG electrodes were used to record EMG signals of biceps brachii and triceps brachii.

sampled at 1,000 Hz with a 16-bit resolution. A custom data acquisition program was used to provide real-time audio and visual feedback to help an examiner control the peak stretch velocity and the peak stretch torque (terminal torque) (26). When the examiner passively moved the subject's forearm, the target velocity and torque profile with boundary lines (e.g., ±10% of target velocity or target torque), as well as the instantaneous joint velocity and torque were displayed on the monitor. At fast velocities (90◦ /sec and 270◦ /sec), the data acquisition program provided audio feedback for controlling peak velocities during passive stretching. At the slow velocity of 30◦ /sec, audio feedback was used to indicate that the designated torque limit had been reached and to stop the passive stretching. By doing so, the applied stretch force could be consistent while quantifying the spastic muscle each time.

A test-retest reliability investigation of the MSE was performed on the healthy controls by the same rater twice with 1 day apart. Excellent reproducibility was found in measures derived from fast stretching (ICC= 0.88) and from slow stretch (ICC= 0.89 for passive ROM, 0.84 for stiffness and 0.75 for energy loss). Excellent intra-rater reliability was also shown while using MSE in the pediatric population (28, 29).

#### Experimental Procedures

In a quiet room, a therapist assessed the spastic elbows of the stroke survivors using MAS (30). The therapist followed the procedure described previously (31) with the exception of the body position. In the current study, the subjects sat upright comfortably instead of lying down. The therapist stabilized the upper arm by holding it proximal to the elbow and moved the forearm in a quick passive motion (∼1 sec) throughout the available elbow range of motion from the end of flexion to maximal extension.

During the experiment, the subject sat next to the MSE with the elbow flexion axis aligned with the rotation axis of the MSE. The shoulder was positioned at 75◦ of abduction and the forearm and upper arm were secured to the supporting braces. Surface EMG electrodes were placed on the biceps brachii and triceps brachii with the reference electrode placed on the lateral epicondyle.

In the clinical practice, the modified Tardieu R1 is the angle of the catch thought to be due to induced stretch reflex at an as fast as possible velocity. The passive ROM (R2) is graded under a slow velocity, which would not trigger the stretch reflex (12). In the current study, we chose 30◦ /sec as the slow velocity to measure R2 and chose two fast velocities (90 and 270◦ /sec) to detect the catch angle (R1) using MSE. The slow velocity of 30◦ /sec was chosen because it did not induce stretch reflex during manual tests that may confound the R2 measurement. Two high velocities (90 and 270◦ /sec) were chosen because 90◦ (right angle) and its multiples are relatively easier for the rater to perceive during manual tests.

Initially, the torque and position offset were recorded with the subject's elbow in the neutral position, defined as the position where subjects felt the most comfortable, not being stretched or restrained (79.6◦ ± 10.9◦ elbow flexion for stroke survivors, and 75.0◦ ± 6.5◦ elbow flexion for healthy controls with full extension defined as 0◦ elbow flexion). To determine the passive ROM and stiffness of the elbow, we then moved the elbow at 30◦ /sec until reaching a pre-defined torque (3 Nm) or a mechanical stop. Three trials of passive stretch separated by 1 min were performed at each of 90◦ /s and 270◦ /sec to evaluate the velocity-dependence of spasticity and catch angle. With practice, the examiner was able to control the peak velocity of passive stretch to match the target velocity as shown on the display (±10% of target velocity). One stretch cycle was defined from full flexion to full extension then back to full flexion. Since spasticity may be altered by repeated stretching (32), the elbow was not stretched more than three cycles in each trial. We also instructed the subjects relax during the tests. If voluntary muscle contractions to assist the stretch were detected by the examiner and shown as intermittent or continuous EMG activities of triceps brachii muscles, and/or EMG activities of the biceps brachii preceding the stretch initiation, the trial was discarded. The examiner would then let the subject rest before stretching the joint again.

The MAS scores, the biomechanical measures and MSEmeasured R1 and R2 were taken by the same examiner.

#### Data Analysis and Statistics

Biomechanical Measurement: Torque and position signals were filtered with a low-pass cutoff frequency of 50 Hz. The derivative of torque with respect to time, dτ (t)/dt was calculated from the acquired torque signals. All the biomechanical measures captured at 30◦ /sec were determined within the torque limits of 3 Nm, including the passive ROM [also described as the R2 angle (12)], elastic stiffness (K), energy loss, and elbow flexor torque at several joint angles. Stiffness, K, was defined as the slope of the torqueangle relationship of ascending arm at 70◦ of elbow flexion, a common range among the subjects. The energy loss was the area between the ascending and descending limbs of the torque–angle curve during rotation (5). Since different subjects had different

ROMs and limb inertias, the energy loss was normalized by the transverse inertia and the individual's range of motion. The transverse inertia was estimated from the length of the forearm, the perimeter of the elbow, and maximum forearm and wrist circumference (33, 34). Joint angle position dependence of the resistance torque was evaluated at three selected elbow flexion angles (45◦ , 60◦ , and 75◦ ) in the two groups and normalized by the body weight and height (35). The angles were near the end of stretch and around the common range (70◦ ).

EMG Signal Processing: The EMG signals were filtered with a passband of 10–450 Hz. To determine the onset of muscle activations (reflex-mediated responses) for verifying the catch angle, EMG signals were presented in a linear envelope (LE) form: with the filtered EMG data full-wave rectified and then low-pass filtered at 10 Hz. The onset of the reflex EMG was determined when the amplitude of the EMG LE was larger than the mean plus three standard deviations (SD) of the background EMG (36). The background EMG was measured during the quiescent period before the passive stretch. The elbow flexion angle corresponding to the reflex EMG onset was thus determined.

Catch Angle Determination and Characterization of spastic muscles: Catch angle is where a sudden occurrence of increased muscle activations in response to a quick passive stretch, which leads to an abrupt stop or increased resistance (torque) before the joint rotation reaches the end of ROM (37). This behavior can be captured using MSE as shown in **Figure 2**: the elbow flexor torque increased with elbow flexors being stretched (flexion angle decreased in **Figures 2A,B**). When the passive stretch triggered a stretch reflex (EMG firing shown in **Figure 2C**), the elbow flexors contracted strongly causing the abruptly increased torque rate (dτ (t)/dt) as indicated by the second positive peak in **Figure 2D** where the catch angle was determined. **Figure 2E** shows that the examiner responded to the abruptly increased resistance with a decreased velocity to avoid over-stretching the joint (value of velocity changed toward zero). We then developed a systematic way to determine R1. Since dτ (t)/dt was affected considerably by the inertia during the initial acceleration (the first peak of dτ (t)/dt), the local minimum of velocity was selected as the first landmark (the dashed vertical line in **Figure 2E**), which occurred shortly after the catch. Next, the joint angle corresponding to the peak dτ (t)/dt in the 300-ms window preceding this landmark was determined as the catch, R1 (26). The ratio R1/R2 was derived through dividing the angular displacement between R1 and flexion angle by the overall ROM, to represent the portion, which is free from the catch. A four-dimensional display, including the variables of joint angle, velocity, torque (τ ), and dτ (t)/dt, was developed to illustrate the events involved in the catch that provides a more comprehensive quantification of spastic muscles group around the tested joint (**Figure 2F**) at various stretch velocities. The width of the shaded area in **Figure 2F** represents velocity. Note the increased width as the high velocity was maintained. The torque increased as the elbow was moved into extension indicated by the dashed arrow. When a catch occurred, dτ (t)/dt increased abruptly and reached a local maximum (the relatively hot color of the dτ (t)/dt line). The abruptly increased resistance and the examiner's reaction to the catch resulted in a

quick velocity reduction to a local minimum (choke, the narrow shaded area).

(F) was developed to illustrate the aforementioned events.

Statistics: Since the data was not normally distributed, nonparametric statistics (Friedman test) were used for comparisons of catch angles at different stretch velocities with a significance level of p < 0.05. To investigate differences in the biomechanical measures (ROM, stiffness, torque at three joint positions, energy loss) and the velocity-dependent properties of muscle compared between the control and CVA groups, a Mann-Whitney Utest with a significance level of p < 0.05 was conducted. Spearman correlation was used for correlating the MAS with all biomechanical measures as well as the catch angle. Correlation was significant at p < 0.05. The intra-class correlation coefficient was chosen as the test statistic to evaluate the test-retest reliability. The two-way mixed model of intra-class correlation coefficient was used. An intra-class correlation coefficient ≥ 0.75 indicated significant reproducibility (38). SPSS software (SPSS Inc., Chicago, Illinois) was used to perform all statistical analysis.

#### RESULTS

#### 4-D Characterization of Catch Angle

**Figures 3A,B** show the representative result of spasticity and catch angle evaluations at velocity of 90◦ /s and 270◦ /s respectively. **Figure 3C** shows the representative curve in a healthy control, the dτ (t)/dt remained at a constant lower value (below 10 Nm/s) throughout the available range of motion. In addition, for healthy controls the joint resistance was much lower when compared to stroke survivors (p < 0.05).

### Dependence of R1 on Stretch Velocity and Velocity-Dependent Properties of Spastic Muscles

During the passive elbow extension, the catch angle in stoke survivors was significantly smaller (elbow more extended, p < 0.001) with increasing stretch velocity (**Figure 4A**). This indicates that the catch angle occurred later at faster stretch velocities. Furthermore, the R1/R2 was significantly higher at higher stretch velocities (**Figure 4B**, p < 0.001). As expected, catch was not observed in any healthy controls. In the representative case shown in **Figures 3A,B**, catch occurred at 78.2◦ and 58.8◦ elbow flexion at the stretch velocities of 90◦ /s and 270◦ /s, respectively.

As a feature of spasticity, the peak resistance torque during a stretch increased with the increasing stretch velocity in the stroke survivors (p < 0.005, **Figure 5**). Healthy controls also showed increased peak resistance at the velocity of 270◦ /s when compared to the velocity of 90◦ /s (p = 0.005). The slope of the relationship between the peak torque and the stretch velocity was significantly higher in stroke (9.34 nu × 10−<sup>5</sup> ± 4 × 10−<sup>5</sup> Nkg−1deg−<sup>1</sup> s) than that in healthy controls (4.99 × 10−<sup>5</sup> ± 2.95 × 10−<sup>5</sup> Nkg−1deg−<sup>1</sup> s; p = 0.02).

### Biomechanical Measures of the Spastic Muscles

ROM measured at a controlled torque of 3 Nm was significantly reduced in the stroke survivors as compared to that of the healthy controls (74.2◦ ± 21.5◦ vs. 107.6◦ ± 8.7◦ , p < 0.001; **Figure 6A**). During extension with a 3 Nm torque limit, the stroke survivors stopped earlier at larger flexion angles (30.0◦ ± 17.6◦ ) compared to healthy controls (10.2◦ ± 10.8◦ ; p < 0.01). **Figure 7** shows the examples of restricted ROMs for the severely spastic (MAS=3), mildly spastic (MAS=1) and healthy controls that were 36◦ -93◦ , 3 ◦ -104◦ , and −1 ◦ -106◦ elbow flexion, respectively.

a representative stroke survivor at 90◦ /sec (A), at 270◦ /sec (B), and a healthy control at 270◦ /sec (C). The X-axis and Y-axis correspond to the elbow flexion angle and torque, respectively. The resistance torque generated by the stretched flexor muscles is positive. The dashed arrow indicates the direction of movement. Velocity is proportional to the width of the gray shaded area. The color of the line gives dτ (t)/dt, with the hot color (red) corresponding to a high rate of dτ (t)/dt. The two small arrows in (A) indicate the local maximum of dτ (t)/dt where the catch occurred and local minimum of velocity resulting from the examiner's reaction to the abruptly increased torque, respectively.

Stiffness measured at a prescribed elbow flexion angle of 70◦ was significantly larger in stroke survivors when compared to healthy controls (**Figure 6C**; 0.058 ± 0.028 Nm/deg vs. 0.017 ± 0.008 Nm/deg, p < 0.001). The stiffness for the severely spastic, mildly spastic and healthy controls was 0.162 Nm/deg, 0.042

Nm/deg, and 0.014 Nm/deg respectively (**Figure 7**, slope of the lines).

Stroke survivors lost more energy (24.14 ± 8.07 Jdeg−1kg−1m−<sup>2</sup> ) than healthy controls (12.66 ± 5.38 Jdeg−1kg−1m−<sup>2</sup> ; p = 0.005; **Figure 6B**). In addition, stroke survivors showed higher torque compared to healthy controls at 45◦ , 60◦ , and 75◦ of elbow flexion (stroke vs. healthy: 0.0116 ± 0.001 Nkg−<sup>1</sup> vs. 0.0062 ± 0.0008 Nkg−<sup>1</sup> , 0.0088 ± 0.0031 Nkg−<sup>1</sup> vs. 0.0045 ± 0.00084 Nkg−<sup>1</sup> , and 0.0069 ± 0.0021 Nkg−<sup>1</sup> vs. 0.0035 ± 0.001 Nkg−<sup>1</sup> , respectively; **Figure 6D**).

### Correlations Between the MAS, R1 and Biomechanical Measures

**Table 2** shows the correlations among the variables evaluated. The catch angle showed significant correlations with ROM (r = −0.685, p = 0.014) and with energy loss (r = −0.679, p = 0.047). The MAS showed significant correlations with ROM (r = −0.897, p = 0.001), as well as stiffness of the joint (r = 0.828, p = 0.006). Reflex-mediated EMG responses at different stretch velocities did not show any significant correlations with MAS, R1 or biomechanical measures (p > 0.05).

#### DISCUSSION

This study demonstrates 4-D plot, a comprehensive, and systematic way, to investigate the group of spastic muscles around the elbow. Spasticity-related biomechanical characteristics, including joint ROM, joint torque, stiffness and energy loss, at various controlled velocities can also be acquired at the same single setting. During a controlled slow stretch, the MSE assesses biomechanical properties of the joint including the ROM, stiffness and energy loss using a controlled slow stretch and determine the catch angle at controlled fast stretch velocities.

Convenient spasticity quantification has been a challenge. In order to evaluate stretch reflex responses accurately, the way to elicit spasticity should be standardized. Many factors including the pre-activation of muscles, position of the joints involved, and applied stretch torque and velocity, as well as the experience of clinicians may result in different outcomes and interpretations. Clinical measures, such as the MAS or Tardieu scales, have been used to identify the catch angle during passive stretch. However, the angles determined using these scales may not be accurate; they generally occur later than a biomechanicallydetected catch and suffer from poor inter-rater reliability (39). As shown in **Figures 2**, **3**, the examiner reacted to the abruptly increased joint resistance by slowing or stopping the passive stretch. However, the catch occurred up to 300 ms prior to this slowing or stopping point. Therefore, peak dτ (t)/dt, instead of the stopping point, should be used as the indicator of the catch angle. In the current study, the instantaneous velocity change along with dτ (t)/dt were used to determine the catch angle reliably and to minimize potential human error. As one can see in **Figures 2**, **3**, the human's reaction characterized as local minimum of speed (choke) was further away from the catch determined by dτ (t)/dt. The discrepancy between human's reaction and the true catch implies the potential human errors in the subjective clinical measures of "catch angle." As seen in the examples in **Figures 2**, **3**, the differences ranged from 3 to 8◦ . In a clinical setting, the catch angle reading is usually from eyeballing of a goniometer moving with the joint, which may introduce even a larger error in determining the catch angle. It should be noted that another peak of torque change rate, which occurred earlier during the passive stretch, was when an examiner overcame the resistance from the limb inertia and was not related to catch angle.

Velocity-dependent increase in muscle tone is a key attribute to spasticity as shown in the current study as well that the stretch velocity has obvious effects on the normalized peak torque of catch (larger slope in **Figure 5**). However, the velocity dependence of spasticity might be partially due to joint angle position dependence. The delayed catch angle associated with fast velocities in our study showed the joint angle position dependence of the increased resistance. At a faster stretch velocity, the joint was quickly stretched further into an angle position where higher resistance existed. Assuming reflexmediated torque developed 60 ms after the stretch reflex was triggered, the joint would have been moved 10.8◦ further in 60 ms at 270◦ /s as compared to the stretch at 90◦ /s ((270◦ /s-90◦ /s) × 0.06 s = 10.8◦ ). The extra 10.8◦ moved the joint to a stiffer position, which might make the examiner feel higher resistance at

FIGURE 5 | Velocity-dependent property of the muscle in subjects post-stroke (solid circles) and healthy controls (solid triangles). The dotted lines demonstrate the dependence on the velocity. Subjects post-stroke had a stronger velocity dependence compared to healthy controls (*p* = 0.02) and had higher torque at all velocities than healthy controls (*p* < 0.001). Asterisks (\*) indicates significant differences between two velocities.

reduced ROM (A), higher energy loss (B), increased stiffness (C), and increased resistance torque at several angles (D), represented as mean ± SE (standard error of mean). Asterisks (\*) indicate significant difference (*p* < 0.05) between the two groups (D). Mann-Whitney U test with significance level of *P* < 0.05 was used to determine the differences.

a faster velocity. Similar results were reported that greater catch angles were associated with higher applied angular velocities (28, 40). Velocity-dependence of the catch angle further confirms that a standardized method to evaluate spastic muscles is essential since the interpretation of catch angle can be confounded by the stretch velocities.

In the current study, passive properties were measured under real-time feedback control by moving the elbow slowly without



*MAS, modified Ashworth Scale; R1, catch angle; EMG, electromyography onset; EL, energy loss; K, stiffness.* \* *indicates P* < *0.05* \*\**indicates P* < *0.01.*

eliciting a reflex response. The lack of a reflex response was corroborated by silence in the EMG signals of the stretched muscle. Significant changes in the passive properties of the spastic elbow of stroke survivors were observed when compared to healthy controls, including increased stiffness and flexors resistance, decreased ROM, and increased energy loss (**Figure 6**). Similar changes were also found in the spastic ankle of stroke survivors with hemiparesis (5). In general, the changes in stiffness and ROM were consistent with what have been reported previously (14, 41, 42). Since the supporting braces fixed to the MSE might hinder the ROM near the end of flexion, the value of elbow ROM shown in current study was smaller than the observations in previous studies (43–45). In addition, **Figure 7** also shows that the reduced ROM was not only in one end, the stroke survivors may lose the range toward flexion or extension that can be relevant to daily functions.

The correlation between passive stiffness and MAS demonstrates that the MAS is more closely related to the passive stiffness of the joint than to joint spasticity, even though it has been commonly used for assessing spasticity in both clinical and laboratory settings. MAS only includes as single stretch velocity and is scored by the amplitude of joint resistance that potentially make MAS reflect passive stiffness over spasticity. Because the felt joint resistance could be from either spastic responses and/or passive stiffness (46, 47) and without various stretch velocities those could not be distinguished. The consequence of this ambiguous assessment may mislabel patients who have increased passive stiffness alone as having spasticity and being treated with inappropriate interventions. Damiano et al. indicated that evaluating patients at different velocities may help to distinguish passive stiffness from spasticity (46), which we adopted in our developed method for spasticity quantification. Administering the Tardieu scale using the MSE could provide proper spasticity characterization under various controlled velocities. The intra-rater reliability of clinical assessments can be improved using the MSE that contains accurate sensors and

provides real-time audiovisual feedback instead of the examiner's subjective manipulation and scoring.

#### LIMITATIONS OF THE STUDY

The brace of the system we used might prevent the all range of motion toward the end of flexion due to the contact of muscle bulks and brace. When there was no compromise using this system to assess the elbow flexors spasticity, one should interpret the passive flexion ROM with cautions. The sample size of the current study is relatively small. A larger size of sample should be considered in a future study.

#### ETHICS STATEMENT

This study was carried out in accordance with the recommendations of Northwestern University Institutional Review Board with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Northwestern University Institutional Review Board.

#### AUTHOR CONTRIBUTIONS

Y-NW, L-QZ, H-SP, and YR designed and carried out the experiment. Y-NW and L-QZ performed the data analysis.

#### REFERENCES


Y-NW, H-SP, YR, and L-QZ wrote the manuscript with support from J-JC and ER.

#### FUNDING

This work was supported in part by grants from the NSF (award no. IIP-0750515) and NIH (award no. R01-HD044295).

for quantification of spasticity. Ann Biomed Eng. (1999) 27:815–29. doi: 10.1114/1.234


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Wu, Park, Chen, Ren, Roth and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Motor Unit-Driven Identification of Pathological Tremor in Electroencephalograms

Aleš Holobar <sup>1</sup> \*, Juan A. Gallego<sup>2</sup> , Jernej Kranjec<sup>1</sup> , Eduardo Rocon<sup>2</sup> , Juan P. Romero3,4 , Julián Benito-León5,6,7, José L. Pons <sup>8</sup> and Vojko Glaser <sup>1</sup>

<sup>1</sup> Faculty of Electrical Engineering and Computer Science, University of Maribor, Maribor, Slovenia, <sup>2</sup> Neural and Cognitive Engineering Group, Centre for Automation and Robotics, Spanish National Research Council, Arganda del Rey, Spain, <sup>3</sup> Neurorehabilitation and Brain Damage Research Group, Experimental Sciences School, Universidad Francisco de Vitoria, Madrid, Spain, <sup>4</sup> Brain Damage Unit, Hospital Beata María Ana, Madrid, Spain, <sup>5</sup> Department of Neurology, University Hospital 12 de Octubre, Madrid, Spain, <sup>6</sup> Center of Biomedical Network Research on Neurodegenerative Diseases, Madrid, Spain, <sup>7</sup> Department of Medicine, Faculty of Medicine, Complutense University of Madrid, Madrid, Spain, <sup>8</sup> Neural Rehabilitation Group, Cajal Institute, Spanish National Research Council, Madrid, Spain

#### Edited by:

Ping Zhou, University of Texas Health Science Center at Houston, United States

#### Reviewed by:

Kemal Sitki Turker, Koç University, Turkey Amit N. Pujari, University of Hertfordshire, United Kingdom Andrés Úbeda, University of Alicante, Spain

> \*Correspondence: Aleš Holobar ales.holobar@um.si

#### Specialty section:

This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology

Received: 29 June 2018 Accepted: 28 September 2018 Published: 29 October 2018

#### Citation:

Holobar A, Gallego JA, Kranjec J, Rocon E, Romero JP, Benito-León J, Pons JL and Glaser V (2018) Motor Unit-Driven Identification of Pathological Tremor in Electroencephalograms. Front. Neurol. 9:879. doi: 10.3389/fneur.2018.00879 Background: Traditional studies on the neural mechanisms of tremor use coherence analysis to investigate the relationship between cortical and muscle activity, measured by electroencephalograms (EEG) and electromyograms (EMG). This methodology is limited by the need of relatively long signal recordings, and it is sensitive to EEG artifacts. Here, we analytically derive and experimentally validate a new method for automatic extraction of the tremor-related EEG component in pathological tremor patients that aims to overcome these limitations.

Methods: We exploit the coupling between the tremor-related cortical activity and motor unit population firings to build a linear minimum mean square error estimator of the tremor component in EEG. We estimated the motor unit population activity by decomposing surface EMG signals into constituent motor unit spike trains, which we summed up into a cumulative spike train (CST). We used this CST to initialize our tremor-related EEG component estimate, which we optimized using a novel approach proposed here.

Results: Tests on simulated signals demonstrate that our new method is robust to both noise and motor unit firing variability, and that it performs well across a wide range of spectral characteristics of the tremor. Results on 9 essential (ET) and 9 Parkinson's disease (PD) patients show a ∼2-fold increase in amplitude of the coherence between the estimated EEG component and the CST, compared to the classical EEG-EMG coherence analysis.

Conclusions: We have developed a novel method that allows for more precise and robust estimation of the tremor-related EEG component. This method does not require artifact removal, provides reliable results in relatively short datasets, and tracks changes in the tremor-related cortical activity over time.

Keywords: pathological tremor, EEG decomposition, surface EMG decomposition, Parkinsonian tremor, essential tremor

**62**

### INTRODUCTION

The role of cerebral cortex in the generation of pathological tremor has been widely studied in essential as well as in Parkinsonian tremor. Accumulated evidence suggests that tremor-related cortical activity exists in both types of tremor (1– 5). Moreover, because of the significant coupling between the cortical activity and the activity in the affected muscles, motor cortex is thought to contribute to tremor generation (2, 3, 6–15).

To the best of our knowledge, all of the existing studies assessed corticomuscular coupling by computing the coherence between cortical activity, recorded with EEG or magnetoencephalograms (MEG), and an estimate of muscle activity derived from the surface EMG (2, 3, 6–13, 15, 16). Since the coherence function reveals a linear relationship between two signals at a given frequency (17), coherence at the tremor frequency is assumed to indicate tremor-related EEG activity. Although robust to noise at different frequencies, coherence only provides an indirect measure of corticomuscular coupling, and does not enable tracking changes in tremor properties over a short time scale. Furthermore, it requires off-line processing of relatively long EEG and EMG recordings, which need to be cleaned of artifacts beforehand. This limits the comparison of the tremor-related cortical activity across conditions and diseases.

Besides the coherence function, the cortical tremor component could also be potentially identified using blind source separation (BSS) algorithms. For example, Delorme et al. (18) identified, using independent component analysis (ICA) techniques, a number of components in the EEG activity of healthy subjects performing a working memory task. However, no group has demonstrated the feasibility of separating the tremor component from other brain activity. Notably, all the ICA/BSS algorithms proposed so far build on general assumptions about the EEG properties such as independence of the identified components, and would not exploit the specific characteristics of tremor, such as the relationship between sensorimotor cortical activity and muscle activity. As a result, they would suffer from large inter-trial and inter-subject variability in convergence toward the specific (tremor-related) EEG component.

With the exception of Gallego et al. (15), where we used cumulative motor unit spike trains (CST) to characterize the neural drive to the muscle, the authors of all of the studies mentioned above used the rectified surface EMG as an estimator of muscle activity. However, recent studies have shown that the CST of several (e.g., ≥5) motoneurons that innervate a muscle provide a more accurate representation of the synaptic inputs to a motoneuron pool than the EMG envelope (19– 21). In addition, rectification of the surface EMG may or may not enhance the detection of synaptic inputs to the pool depending on the muscle contraction level (22). Indeed, as shown in Farina et al. (22) rectification is preferable over the raw EMG only at low contraction levels. Therefore, methods based on the CST rather than on the traditional surface EMG analysis are likely to identify tremor-related cortical activity more reliably.

In this study, we present and validate a novel method to identify tremor-related cortical activity. This method builds on the assumption that the tremor-related EEG component is stochastically phase-locked to the motor unit firings in a tremulous muscle. This implies close to linear relationship between the cortical activity and the firings of the pool of motor units that form a muscle (15, 21, 23). This relationship has been experimentally demonstrated before by the existence of EEG-EMG coherence at the tremor frequency (2, 3, 6–13, 15, 16, 24).

In our method, motor unit spike trains are identified from non-invasively recorded multichannel surface EMG recordings using the Convolution Kernel Compensation (CKC) technique (25–27). These firings are then used to construct the phaselocked estimator of the tremor-related activity in EEG. By using simulated data, we show that our method tracks reliably the tremor activity for a wide range of physiologically realistic conditions. We then apply this method to recordings from nine essential tremor (ET) and nine Parkinson's disease (PD) patients, and show that the method outperforms the traditional coherence approach when detecting tremor-related cortical activity.

### MULTICHANNEL EMG-DRIVEN IDENTIFICATION OF TREMOR-RELATED ACTIVITY

We first assume that EEG signals are linear mixtures of brain oscillations (rhythms) and noise (**Figure 1**). Then, using this data model, we build a linear minimum mean square error (LMMSE) estimator of the cortical tremor activity that exploits the large synchronization of motor unit firings in pathological tremor.

#### Data Model

Assume the mixing model depicted in **Figure 1**, and denote the M EEG channels by **y** (n) = - y<sup>1</sup> (n), y<sup>2</sup> (n). . . yM(n) T , where the n-th sample of the i-th channel appears in the i-th row of **y** (n). The model inputs, sj(n), represent the brain rhythms in the EEG. For example, in a normal condition, these inputs would reflect the alpha, beta or gamma rhythms. As in other BSS/ICA EEG studies, artifacts, such as blinking artifacts, would also be considered a model input. In the case of pathological tremor, one or more of these inputs reflects tremor activity.

The mixing model in **Figure 1** can be expressed in a matrix form as.

**y** (n) = **As**(n) + ω(n) (1)

Where

$$\overline{\mathbf{s}}\left(n\right) = \left[s\_1\left(n\right), s\_1\left(n-1\right)\dots s\_1\left(n-L+1\right), s\_2\left(n\right)\dots\right]$$

$$s\_2\left(n-L+1\right)\dots s\_l\left(n-L+1\right)\right]^T\tag{2}$$

contains blocks of L samples from J sources and the noise vector ω (n) = - ω<sup>1</sup> (n), ω<sup>2</sup> (n). . . ωM(n) T is modeled as a zero-mean

FIGURE 1 | Proposed mixing model of the EEG. CST is coupled to the tremor EEG component s1(n) via an unknown function f(n). The other sources [s2(n) to sN(n)] represent the non-tremoric brain rhythms in the EEG. The impulse responses amj project the j-th source s<sup>j</sup> (n) to the m-th EEG channel ym(n), whereas ωm(n) denotes the additive noise in the m-th EEG channel.

ergodic Gaussian process, spatially and temporarily independent of the activity in **s**(n). The mixing matrix

$$\mathbf{A} = \begin{bmatrix} \mathbf{a}\_{11} \cdots & \mathbf{a}\_{1f} \\ \vdots & \ddots & \vdots \\ \mathbf{a}\_{M1} & \cdots & \mathbf{a}\_{Mf} \end{bmatrix} \tag{3}$$

contains stationary impulse responses **a**mj = - amj(0) · · · amj(L − 1) that convolve the j-th input source and add the result to the m-th EEG channel.

The model described in Equation (1) is typically underdetermined, with more inputs than measurements. However, low energy inputs can always be modeled as physiological noise and, thus, included in ω (n). We can also extend **y** (n) by adding F-1 delayed repetitions of each EEG channel:

$$\overline{\mathbf{y}}\left(n\right) = \left[\mathbf{y}\_1\left(n\right), \mathbf{y}\_1\left(n-1\right) \dots \mathbf{y}\_1\left(n-F+1\right), \mathbf{y}\_2\left(n\right) \dots \right. \\ \left. \left. \left. \left. \left. \left. \mathbf{y}\_n \right\|\right|\right) \dots \mathbf{y}\_2\left(n-F+1\right) \right] \right. \\ \left. \left. \left. \left. \left. \mathbf{y}\_n\left(n-F+1\right) \right|\right) \right] \right. \\ \tag{4}$$

This increases the number of outputs in model (1) and, more importantly, compensates potentially different time delays of same source in different EEG channels (see EMG-Driven Decomposition of EEG section for details). This is important because the same rhythm may be present at different EEG electrodes at close but different time lags (this assumption is further confirmed by the results in Results section). The extended input vector **s** and mixing matrix **A** now change to:

$$\overline{\mathbf{s}}\left(n\right) = \left[s\_1\left(n\right)\dots s\_1\left(n-L-F+2\right), s\_2\left(n\right)\dots\right]$$

$$s\_2\left(n-L-F+2\right)\dots s\_f\left(n-L-F+2\right)\right]^T \qquad \text{(5)}$$

$$\overline{\mathbf{A}} = \begin{bmatrix} \mathbf{A}\_{11} & \cdots & \mathbf{A}\_{1f} \\ \vdots & \ddots & \vdots \\ \mathbf{A}\_{M1} & \cdots & \mathbf{A}\_{Mf} \end{bmatrix} \tag{6}$$

with

$$\mathbf{A}\_{m\circ} = \begin{bmatrix} \mathbf{a}\_{m\circ} & \cdots & \mathbf{0} \\ \vdots & \ddots & \vdots \\ \mathbf{0} & \cdots & \mathbf{a}\_{m\circ} \end{bmatrix}$$

For the purpose of mathematical derivations in **Appendix**, we will further represent the EEG recordings **y**(n) as analytic signals. Note that this does not alter the assumed mixing model, and can be readily fulfilled by applying Hilbert transform to the EEG signals.

#### EMG-Driven Decomposition of EEG

By temporarily neglecting the impact of noise ω (n), the LMMSE estimator of the input sj(n) is given by (25)

$$\begin{split} \hat{\mathbf{s}}\_{\circ}(n) &= \mathbf{c}\_{s\_{\circ}\overline{\sf s}}^{H} \mathbf{C}\_{\overline{\sf s}}^{-1} \overline{\sf v}\left(n\right) = \mathbf{c}\_{s\_{\circ}\overline{\sf s}}^{H} \overline{\sf A}^{H} \left(\overline{\sf A}\,\mathbf{C}\_{\overline{\sf s}}\,\overline{\sf A}^{H}\right)^{-1} \overline{\sf A}\overline{\sf s}\left(n\right) \\ &= \mathbf{c}\_{s\_{\circ}\overline{\sf s}}^{H} \mathbf{C}\_{\overline{\sf s}}^{-1} \overline{\sf s}\left(n\right), \end{split} \tag{7}$$

where superscript <sup>H</sup> denotes conjugate transpose, **c** <sup>s</sup>j**y**=<sup>E</sup> sj(n)**y**(n) is the cross-correlation vector between sj(n) and **y**(n), and **c** <sup>s</sup>j**s**=<sup>E</sup> sj(n)**s**(n) is the cross-correlation vector between sj(n) and **s**(n). Matrices **C<sup>s</sup>** and **C<sup>y</sup>** denote the correlation matrices of **s**(n) and **y**(n), respectively. This LMMSE estimator is Bayesian optimal in the minimum square error sense, also in the presence of noise, but requires a priori knowledge of the cross-correlation vector **c** sj**y** . In experimental conditions **c** sj**y** is not known and needs to be estimated from the measurements.

To obtain an estimate of **c** sj**y** , we developed a method that is based on the assumption that, in an affected muscle, the motor unit firings are phase-locked to the tremor-related cortical activity. This assumption follows from the observations of significant coherence (linear relationship) between the EEG and the rectified EMG (1–10, 13, 14, 16) or, more directly, between the EEG and the motor unit CST (15). Remarkably, the number of motor units needed to represent population behavior during tremor is not very large, which in practice will favor the implementation of the proposed method. For example, in Negro and Farina (21) and Gallego et al. (15), the coherence between the CST and the EEG was nearly maximal with only five motor units and did not increase significantly when more motor units were added to the CST.

To obtain an estimate of **c** sj**y** , we first define the firing pattern of the r-th tremor-affected motor unit:

$$t\_r(n) = \sum\_{k=0}^{K\_r} \delta \left( n - k \frac{f\_s}{f\_r} - d\_r - \Delta d\_{rk} \right), \qquad r = 1, \ldots, R,\tag{8}$$

where δ(τ ) denotes the unit-sample impulse, d<sup>r</sup> is the common time lag (in samples) of the pulses in tr(n) due to the transmission from motor cortex, 1drk is the intrinsic motor unit firing time variability (in samples), and is frequently defined as the k-th realization of Gaussian random variable 1drk ∼ N(0, σ1d<sup>r</sup> ), f<sup>r</sup> is the motor unit firing frequency, f <sup>s</sup> is the sampling frequency and K<sup>r</sup> is the number of firings in the observed time interval. Note that, in order to simplify analytical derivations, we ignored the doublets in the motor unit firing pattern (8). This has small impact on the results presented herein, as doublets can always be modeled as two different instances of the k-th motor unit firing with different 1drk values. Furthermore, as shown in Results section, when respecting its physiologically induced range, 1drk has almost negligible impact on the proposed tremor estimation.

The cross-correlation between the EEG component sj(n) and **y** (n) can be estimated as (see derivation in the **Appendix**)

$$
\hat{\mathbf{c}}\_{\circ \overline{\mathbf{y}}} \approx \frac{\alpha}{N} \sum\_{r=1}^{R} \sum\_{\eta=1}^{N} t\_r(\eta) \,\overline{\mathbf{y}}(\eta) \tag{9}
$$

where N is the number of signal samples and factor α is defined in the **Appendix**). Due to the amplitude ambiguity of source components extracted by BSS (28), the scalar factor α can be neglected and the EEG-component that is phased-locked to the firings of the J motor units expressed as

$$
\hat{\mathbf{s}}\_{\circ}(n) = \hat{\mathbf{c}}\_{s\_{\circ}\overline{\mathbf{y}}}^H \mathbf{C}\_{\overline{\mathbf{y}}}^{-1} \overline{\mathbf{y}}(n) \tag{10}
$$

By knowing sˆ<sup>j</sup> (n) cross-correlation vector **c**ˆsj**<sup>y</sup>** in (9) may be recalculated as

$$
\hat{\mathfrak{c}}\_{\circ\_{\vec{\mathcal{Y}}}} \approx \mathcal{F}^{-1}\left(\lg\left(\hat{\mathbb{S}}\_{\circ}(f)\right)\right) \overline{\mathbf{y}}(\eta) \tag{11}
$$

where Sˆ j f = F sˆ<sup>j</sup> (n) is Fourier transform of sˆ<sup>j</sup> (n) and g (x) = x∗ |x| denotes the element-wise product of element x with its absolute value. Equations (10, 11) are then iteratively recalculated until the convergence is reached. In our study, iterations were stopped when the second norm of sˆ<sup>j</sup> (n) changed for <0.1%. These iterations emphasize the spectral peaks in the extracted EEG component and suppress the wideband frequency components with low energies.

Note that (9) is only true when none of the other oscillatory inputs in **s**(n) is a higher harmonic of the input sj(n). In the opposite case, (10) would identify both the tremor-related EEG component and its higher harmonics as one joint input (see explanations in the **Appendix**).

The presented method still needs a good approximation of the motor unit firing times, tˆ <sup>r</sup> (n), in order to accurately estimate **c**ˆsj**y**. This can be obtained from high-density surface EMG recordings using the CKC decomposition technique (25, 27, 29–31), which has already been demonstrated to be highly robust to high levels of motor unit synchronization (27).

#### SIMULATIONS AND EXPERIMENTAL RECORDINGS

The presented method was validated on a set of synthetic signals and on experimental recordings from 18 tremor-affected patients.

#### Simulations

First, we tested the proposed method in a simple model that generated EEG-like oscillations as mixtures of sinusoidal sources with time-varying amplitudes. Our goal was to test the method's ability to accurately reconstruct such EEG-like sources from their convolutive mixtures, and to study its sensitivity to its three main parameters: extension factor F in (4), motor unit firing variability 1drk and signal-to-noise ratio (SNR).

The simulated signals comprised 10 (J = 10) mutually orthogonal sinusoids s<sup>j</sup> (n) and their first higher harmonics as input signals:

$$s\_{\hat{j}}(n) = \mathbf{a}(n) \cdot \left(\mathcal{B} \cdot \sin\left(2\pi f\_{\hat{j}} n - \phi\_{\hat{j}}\right) + H\_1 \cdot \sin\left(4\pi f\_{\hat{j}} n - \phi\_{\hat{j}}\right)\right),\tag{12}$$

where a(n) is an amplitude modulator generated by filtering white noise with a second order low-pass Butterworth filter with cut-off frequency of 1 Hz. The amplitude B was set equal to 1, whereas the amplitude of the first harmonic, H1, was varied across simulations and was set to 0, 0.2, 0.4, 0.6, 0.8, and 1. This way we simulated the experimentally observed ratios between the basic tremor frequency and its first harmonic (9, 13, 32). The frequency f<sup>j</sup> of the oscillatory inputs was set to 5+j/2 Hz with j = 1, 2 . . . 10, and the phase φ<sup>j</sup> was randomly selected from the interval [0, 2π]. The sampling frequency was set to 1024 Hz and each simulation lasted 30 s. We assumed that the first oscillatory input, s1(n), represented the tremor-related component we wanted to detect.

Next, we simulated the spike trains of ten motor units, tr(n) with r = 1, 2 . . . 10, by finding the local maxima of the first generated oscillatory source, s1(n). We imposed a corticospinal delay d<sup>r</sup> = 10 ms in Equation (8) to simulate the physiological delays (due the transmission from motor cortex to the output of the motoneuron pool) between the tremor-related EEG component s1(n) and the simulated motor unit firing patterns t<sup>r</sup> (n). Finally, the firing variability 1drk of each individual motor unit t<sup>r</sup> (n) was modeled as Gaussian random variable 1drk ∼ N(0, σ1d<sup>r</sup> ) (33, 34) (see **Appendix**). The standard deviation σ1d<sup>r</sup>

was set to 0, 10, and 20% of the average inter-spike interval in t<sup>r</sup> (n), respectively (35–38). Ten Monte Carlo simulation runs were performed for each value of H<sup>1</sup> and σ1d<sup>r</sup> , resulting in a total of 180 (10 × 6 × 3) simulation runs.

We computed 15 synthetic channels **y**(n) in Equation (1), with the mixing matrix **A** having dimensions 15 × 50. The length of the impulse responses **a**mj(l) in (3) was set to L = 5 samples. To simulate different delays in the representation of the oscillatory inputs s<sup>j</sup> (n) across the 15 channels, each element **a**mj had one randomly selected element set to a non-zero random value from a normal distribution N(0, 1), whereas the remaining four elements were set to zero. In each simulation run, five different realizations of zero-mean random Gaussian noise were added to the measurements, with the SNR ranging from 0 to 20 dB in steps of 5 dB. This resulted in 900 (180 × 5) simulated sets of signals. Note that our simple generative model does not realistically represent actual EEG signals; our goal was to use the simulation results to identify the range of parameter values that was adequate for the experimental data analysis. All the simulations were performed in Matlab version 8.6.0.267246 (R2015b).

#### Experimental Recordings

We recorded data from nine ET patients (four females and five males; age, mean ± SD: 70 ± 6 years, range 61–79 years) and nine PD patients (three females and six males; age, mean ± SD: 64 ± 14 years, range 44–88 years) at Hospital 12 de Octubre, Madrid, Spain. In the ET patients, tremor severity ranged from mild (two patients) to severe (three patients), with a mean score of 36 ± 12 (mean ± SD; range 20–51) according to the Fahn-Tolosa-Marin scale (39). In the PD patients, tremor severity also ranged from mild (five patients) to severe (two patients), with a mean score of 12 ± 6 (mean ± SD; range 5–23) according to the UPDRSIII scale (40). All the participants included in the study gave their written informed consent after full explanation of the procedure. The study, which was conducted in accordance with the principles of the Helsinki declaration of 1975, was approved by the ethical standards committee on human experimentation at the University Hospital "12 de Octubre" (Madrid).

The experimental protocol comprised three repetitions of a 30 s long postural task, during which patients kept their arms outstretched, parallel to the ground for 30 s. During these tasks, we measured EEG, multichannel surface EMG, and wrist kinematics. EEG was recorded with 32 passive Au or active Ag/AgCl electrodes (depending on the session) placed on a cap that fulfilled the extended 10/20 system (g.Tec GmbH, Graz, Austria). Electrodes were placed in the following positions: AFz, F3, F1, Fz, F2, F4, FC5, FC3, FC1, FCz, FC2, FC4, FC6, C5, C3, C1, Cz, C2, C4, C6, CP5, CP3, CP1, CPz, CP2, CP4, CP6, P3, P1, Pz, P2, and P4. The ground was placed on the AFz position, with linked earlobes used as a reference. Before positioning the electrodes the skin was slightly abraded with abrasive paste (Meditec–Every, Parma, Italy) and conductive gel (Meditec– Every, Parma, Italy) was put under the electrodes. The recorded signals were amplified (gUSBAmp, g.Tech GmbH, Graz, Austria), band–pass (0.5–60 Hz) and notch filtered (50 Hz), to remove power line interference, and sampled at 256 Hz with 24 bit resolution.

Right wrist kinematics were recorded with inertial measurement units (IMUs) comprising three-axis accelerometers, magnetometers and gyroscopes (Tech MCS, Technaid S.L., Madrid, Spain). These sensors were fixed with surgical tape over the hand dorsum and the distal third of the forearm (on the dorsal side, close to the wrist), respectively, with one of their axes aligned to that of the wrist. Data were sampled at 100 Hz by a 12-bit A/D converter and low pass filtered (<20 Hz). Wrist kinematics was assessed as the difference between the measured accelerations in the axis parallel to the wrist (41).

Surface EMG signals were recorded from the right wrist flexors and extensors with 13 × 5 electrode grids (LISiN– OT Bioelettronica, Torino, Italy, 8 mm interelectrode distance). The electrode grids were centered over flexor carpi radialis and extensor digitorum communis, respectively. Before the placement of the electrode grid, the skin was lightly abraded using abrasive paste (Meditec–Every, Parma, Italy) and cleansed afterward. Electrical conductivity was ensured by filling each of the electrodes in the grids with conductive gel (Meditec– Every, Parma, Italy). A soaked bracelet placed around one of the wrists was used as reference. The surface EMG signals were amplified as bipolar recordings along the direction of the fibers, band-pass filtered (3 dB bandwidth, 10–750 Hz), and sampled at 2,048 Hz by 12–bit A/D converter (LISiN–OT Bioelettronica, Torino, Italy). We synchronized the EEG, EMG, and IMU recordings using a common clock signal, which was fed into each acquisition systems. The rising edge of the first and last clock signal pulses were identified using a purposely-developed Matlab script. Data were then resampled to 2,048 Hz using a routine that incorporated an anti-aliasing filter.

#### Data Analysis

In the experimental recordings, individual motor unit firing patterns were identified from the multichannel surface EMG using the CKC algorithm (25, 27, 29). The pulse-to-noise ratio metric (PNR) (42) was used to assess the accuracy of firing estimation for each identified motor unit. Only reliably identified motor units (PNR > 30 dB; accuracy of motor unit firing estimation >90%) were used for further analysis (42), whereas all the remaining motor units were discarded.

We estimated the tremor EEG component with Equation (7), using the firings of all the identified (experimental recordings) or simulated (simulations) motor units to estimate the crosscorrelation **c**ˆsj**<sup>y</sup>** between the oscillatory components and the EEG signals, as defined in Equation (10). We tested different extension factors, from F = 1 to F = 15. Due to the amplitude ambiguity (see **Appendix**), the estimate sˆ<sup>1</sup> (n) was further normalized to yield a unit norm. Finally, Equations (10, 11) were iteratively applied until the convergence criterion was reached (the second norm of sˆ<sup>j</sup> (n) changed for <0.1%).

The delay between the motor unit CST and the estimated tremor EEG component sˆ<sup>1</sup> (n) was estimated as the argument of the cross spectrum at the basic tremor frequency f<sup>b</sup> (17, 43).

In the simulations, we further assessed the accuracy of our method by computing the normalized mean square error (NMSE; Equation 13) and the cross-correlation coefficient (CC) between the simulated and estimated tremor inputs, sˆ<sup>1</sup> (n) and s<sup>1</sup> (n), after both signals were aligned in time:

$$\text{NMSE} = \frac{\sum\_{n=1}^{N} \left( s\_1 \left( n \right) - \hat{s}\_1 \left( n \right) \right)^2}{\sum\_{n=1}^{N} \left( s\_1 \left( n \right) \right)^2} \cdot 100 \tag{13}$$

In the experimental data, we compared the coherence between the extracted tremor component and the CST in each muscle to the coherence between the CST and the Laplacian-filtered EEG (15) (without any artifact rejection applied). In all these cases the CST was smoothed by convolving it with a 25 ms long Gaussian window. The 99% confidence limit of the coherence function was obtained as (17):

$$1 - (0.01)^{1/(L-1)} \tag{14}$$

where L is the number of disjoint 1-s segments used in the spectral estimation.

Finally, we computed the relative power H1/(H1+B) of the first tremor harmonic with respect to the basic tremor frequency in the estimated tremor EEG component, sˆ<sup>1</sup> (n).

Due to their non-Gaussian distribution (Lilliefors test, p > 0.05), the non-parametric Kruskal–Wallis test was used to compare the differences between the ET and PD patient groups, whereas the Wilcoxon signed rank test was used for paired comparisons. The significance level was set to p < 0.05 and p < 0.01, respectively (see Results section for details).

FIGURE 2 | Estimation of the simulated tremor component. (A) Representative example showing that the estimated tremor component (blue line) was similar to the simulated source (red line). In this example, simulation parameters were SNR = 10 dB and σ <sup>1</sup>d<sup>r</sup> = 10%. (B,C) impact of different values of the extension factor F on the estimated cross-correlation coefficient (CC), the NMSE between estimated and reference tremor component, and the H1/(H1+B) ratio. Results are averaged over 10 simulation runs. Mean values are depicted as thick blue lines, standard deviations as black dashed lines. In the H1/(H1+B) ratio plots, the simulated reference values are depicted with red dashed lines.

## RESULTS

negligibly small.

#### Simulated Data Analysis

**Figure 2A** shows a representative example of detection of the tremor component in the simulated signals, demonstrating that the proposed method accurately extracts the tremor-related component from the synthetic signals. **Figures 2B,C** show the correlation between the estimated and simulated sources, the NMSE and the estimated H1/(H1+B) ratio as a function of extension factor F, for two different SNR and σ1d<sup>r</sup> values. Note the high correlation and small NMSE between the simulated and the estimated tremor-related component, and the accuracy with which the ratio H1/(H1+B) was detected when extension factor was set to F = 3 or higher. For a SNR of 20 dB, extension factors from F = 2 to F = 5 were optimal, whereas for lower SNRs, extension factors from F = 4 to F = 9 were optimal (**Figure 2C**). In both cases the estimated H1/(H1+B) ratio was largely independent from the extension factor in the interval F = 3 to F = 10. All metrics degraded slightly when the model looked too many samples back in the past (F > 10). Based on these results and on the coherence analysis of experimental data (**Figure 6**), we selected an extension factor F = 8 for further analyses. Note that **Figure 6** indicates that our results would hold across a broad range of values of F, from F = 5 to F = 10.

**Figures 3A,B** summarize the NMSE, CC and H1/(H1+B) ratio as functions of the SNR and motor unit firing time variability σ1d<sup>r</sup> at two different simulated H1/(H1+B) ratios. In both cases, the NMSE decreased and the CC increased as the SNR increased, whereas they did not change significantly with σ1d<sup>j</sup> , which suggests that the simulated intrinsic variability in motor unit firing did not affect the source estimation. The estimated H1/(H1+B) ratio did not change significantly with the SNR or σ1d<sup>r</sup> and was always very close to the simulated H1/(H1+B) ratio. Namely, when averaged over all simulated SNRs, σ1d<sup>r</sup> and H1/(H1+B) ratios, the difference between the simulated and estimated H1/(H1+B) ratio was 0.01 ± 0.05.

The delay between the estimated and simulated sources, sˆ<sup>j</sup> (n) and s<sup>j</sup> (n) was largely independent of the SNR, σ1d<sup>r</sup> , and the simulated H1/(H1+B) ratios, averaging 0.4 ± 1.4 ms across all combinations of parameters. When a 10 ms corticospinal delay was imposed between the motor unit firings and their cortical drive s<sup>j</sup> (n), the estimated delay averaged 11.0 ± 1.6 ms. This implies a 1.0 ± 1.6 ms difference with the simulated 10 ms delay. Despite this estimate being quite accurate, we want to note that the current simulations do not generate signals as complex as the recorded EEGs. Nor did the simulations incorporate the delays due to propagation of the motor unit action potentials along the muscle fibers or due to EMG decomposition with CKC (25–27). All these factors contribute to the unknown global delay between 0 and ∼15 ms and cannot be easily estimated in experimental conditions. Thus, we do not expect the experimental estimates of corticospinal delay to be as accurate as those obtained based on the model.

extensors (red squares) and flexors (blue squares) along with the corresponding smoothed CST (red and blue solid lines; the blue trace is inverted for representation purposes). Each square denotes a single motor unit firing. Firings of different motor units are depicted with different vertical offset; (C) power spectrum of the smoothed CST (red solid line for wrist extensors and blue solid line for wrist flexors) compared to that of the wrist acceleration data (black solid line).

TABLE 1 | Summary of the properties of the motor units, identified by multichannel EMG decomposition.


The table shows the number of identified motor units with PNR > 30 dB, their number of firings, and the PNR, for right wrist extensors (EXT R) and flexors (FLE R) of 9 the ET and 9 PD patients. Results are reported as mean ± SD (range).

#### Experimental Recordings

**Figure 4** shows an example of EMG decomposition in a representative ET patient, along with the smoothed CST. **Table 1** summarizes the number of motor units that were detected from the surface EMG with PNR > 30 dB and then used for the identification of the tremor EEG component. On average, 7.7 ± 5.2 and 8.6 ± 6.3 motor units per contraction were identified for the ET and PD patients, respectively. The average number of firings per motor unit was 160 ± 100 for the ET and 218 ± 130 for the PD patients (**Table 1**). In each contraction, all the accurately identified motor units per muscle were used to estimate the tremor EEG component (see **Appendix**). Since we recorded EMGs from the wrist extensors and flexors of the right arm, two estimates of tremor EEG components were extracted per each task repetition.

**Figure 5** shows a representative example of the estimated tremor EEG component. **Figure 5A** depicts the estimated tremor EEG component and how it relates to the smoothed CST, both in the time (left plots) and the frequency domain (right plot). The estimated EEG component exhibits clear tremorrelated activity with peaks both at the basic frequency and the first harmonic of that observed in the pooled motor unit firings and the wrist kinematics. **Figure 5B** shows time-frequency domain contour plots of the extracted EEG component, the smoothed CST and the wrist kinematics, as reference. In the presented case, the tremor-related EEG activity preceded the

pooled motor unit activity and the observed wrist movements, which manifested simultaneously, but this was not always the case.

We investigated how the number of EEG samples in the time domain (extension factor F) influenced the accuracy of the estimation of the tremor-related cortical activity. To this end, we computed the coherence between the tremor EEG component and the smoothed CST for increasing values of F (from F = 1 to F = 15). As shown in **Figure 6A**, the coherence first increased, but it saturated around extension factor F = 8, in agreement with the simulations with lower SNR ratios. Note, however, that the increase in coherence as more EEG samples were included was not significant after F = 5. **Figure 6B** demonstrates that the proposed method significantly outperforms the classical coherence between Laplacian-filtered EEG and spatially averaged rectified EMG, low pass filtered at 15 Hz.

**Figure 7** shows an example of how the proposed method performs compared to the traditional approach of computing the coherence between an estimate of muscle activity (in this case the smoothed CST) and the spatially filtered EEG signals. We calculated the coherence function for the entire recordings (1 s windows with 50% overlapping), and found no significant coherence values in any EEG channel. In contrast, the coherence between the CST and the EEG component extracted using the proposed method (extension factor: F = 8) was significant. These findings generalized to the all the tested ET and PD patients (**Figure 6B**). Our proposed method thus outperforms classic coherence approaches.

In 52% of cases studied (28 of 54), the tremor-related EEG component preceded the CST by 11.0 ± 6.4 ms, whereas in the remaining 48% of cases, the CST preceded the extracted EEG component by 11.0 ± 5.9 ms. All the delays were clustered on the interval between −30 ms and + 30 ms, and we observed no significant difference between the delays in PD and ET patients (P > 0.05, Kruskal–Wallis test). These latency values are in agreement with previous studies (7–10, 13), notwithstanding the limitation listed in the Discussion section.

Finally, a significant difference was observed in the H1/(B+H1) ratio of the extracted tremor EEG between PD and ET patients (**Figure 8**), also in agreement with previous studies (9, 13).

### DISCUSSION

In this study, we derived and validated a new method for the extraction of the tremor-related EEG activity in the case of pathological tremor. The method builds on the physiological coupling between the tremor-related cortical activity and the neural drive to the muscle (the output of the motoneurons that innervate a muscle). In particular, our method combines the motor unit spike trains identified in the decomposition of high-density surface EMG recordings to build an estimator of the tremor-related EEG component. We applied it to EEG recordings to demonstrate its feasibility, but it could also be used for analyzing magnetoencephalographic (MEG) data.

The proposed method was tested on simulated data and on recordings from 9 PD and 9 ET patients. In the simulations, our method detected the simulated tremor component with great accuracy, as indicated by the low NMSE and high crosscorrelation values. The small difference between the simulated and estimated H1/(B+H1) ratio (**Figures 2**, **3**, global average error of 0.006 ± 0.053 for simulated H1/B, ranging from 0 to 1) further demonstrates the fidelity of the estimated tremor component. Our method also yielded very accurate estimates of the delay between the motor unit population activity and the simulated EEG (the average error was 1.0 ± 1.6 ms for a simulated delay of 10 ms).

In the experimental data, the extracted tremor EEG component exhibited clear similarities with the recorded kinematics and motor unit population activity, both in the time and frequency domains. The ground truth about the estimated EEG tremor component is not available in experimental conditions. However, we believe our method performed well because the estimated EEG component exhibited significantly larger coherence with the identified population of motor units than the spatially filtered EEG signals, which is the standard approach (1–10, 13, 16). This observation indicates that the proposed method is likely to help studying the neural mechanisms of tremor. Indeed, our method always identified tremor-related activity in the EEG, while in many of the investigated cases (34 of 54) we did not find significant coherence between the spatially filtered EEG signals and the identified population of motor units. Note that this observation is in agreement with reports that several tens of second long recordings are needed to obtain robust results in standard

recordings. The top right plot shows the spectrum of the smoothed CST form right wrist extensor of an ET patient performing the postural task (EXT). The top left plot shows the coherence between the CST and the tremor EEG component as estimated by the proposed method, while the remaining 12 plots show the coherence between the Laplacian filtered EEG signals and the CST. The label on top of each plot indicates the central EEG electrode. Red dashed lines represent the 99% confidence limit. Note the increase in coherence yielded by our method compared to the traditional coherence approach. Coherence calculated by traditional approach was not significant in any of the electrodes.

coherence analysis (2, 6–10, 13), whereas our datasets were only 30 s long.

Movement artifacts are an important potential confound when studying corticospinal coupling using coherence techniques. We performed two complementary analyses to discard the presence of movement artifacts. First, we tested whether the tremor components were present across many spatially filtered EEG channels, as it would be the case if they resulted from movement artifacts. We calculated the coherence between the CST and each spatially filtered EEG channel. As reported above, in 34 out of 54 cases we did not find significant coherence between the spatially filtered EEG signals and the identified population of motor units. In the remaining 20 cases, significant coherence at the tremor frequency or at its higher harmonics was observed on one or two EEG channels only. As a second control, we examined whether the EEG-CST delays depended on the basic tremor frequency. Finding a significant association between these two parameters would indicate a potential mechanical coupling. Our results ruled out this possibility: the EEG-CST delays lied within the −30 to 30 ms interval. These values did not overlap with the range of delays potentially indicating an artifact (from ∼45 to 100 ms; interval defined by the maximum and minimum tremor frequencies, 11 and 5 Hz, respectively). Therefore, our control analyses indicate that the identified tremor component is unlikely to originate from movement artifacts. Note that these extensive tests are necessary every time the presented methodology is used as, similar to the classical coherence analysis, the movement artifacts could completely mask any tremor-related activity in cortex.

The only parameter that needs to be chosen in our method is the extension factor F in Equation (4). In simulated conditions,

the optimal extension factor was dependent on the SNR (of the input sources in the simulated EEG signals) with the optimal values between F = 2 and F = 5 for SNR of 20 dB, whereas larger values, between F = 4 and F = 9, proved to be optimal in the case of lower SNRs (**Figure 2**). We decided to choose extension factor F = 8 for subsequent analyses. This choice was confirmed during the analysis of the experimental signals, because the coherence between the estimated tremor-related EEG component and the smoothed CST reached the plateau region at F = 8, whereas the overall increase in coherence was not significant after F = 5 (**Figure 6A**). The observation that in all the cases studied F = 1 yielded significantly worse results (lower coherence), indicates that the extension of the convolutive model (1) helps in coping with either the existence of different delays in the representation of brain rhythms across different EEG channels, or with the convolutive nature of EEG mixtures. Since the computational complexity of the proposed method increases with the square of the extension factor F, it is to our advantage that the preferred value of F is relatively small.

Regarding the neurophysiological results of this study, we found that the relative power of the first tremor harmonic compared to the basic tremor frequency is greater in PD than ET patients, regardless of the investigated muscle (**Figure 8**). This is in agreement with other studies using EEG-EMG coherence (9, 13). The observed delay between the estimated tremor EEG component and the pooled motor unit firings also agrees with previously reported values. Several studies in ET and PD patients reported a bidirectional interaction between the primary sensorimotor area of cortex and the affected muscles, with an efferent and afferent delay between 10 and 30 ms (7–10, 13). In our dataset, the EEG activity preceded the motor unit firings in half of the cases, and in the other half followed it. This is likely due to the fact that the primary motor and sensory cortices are next to each other and the limited spatial resolution of the EEG makes their activities hard to disentangle. We want to emphasize that these results were obtained using significantly shorter datasets (30 s vs. the typically ≥60 s long signals employed in other studies), and avoiding the need of manually discarding epochs with artifacts.

Results in **Figure 6B** suggest that the strength of our method derives from the direct use of the CST in the identification of the tremor-related cortical activity. One of the reasons for this is that the CST provide a more accurate representation of the common synaptic input to the muscles than rectified EMG as it eliminates the influence of frequency components introduced by the motor unit action potentials (22, 44). In the case of tremor, the CST has most of its power at the frequency of the tremor and its harmonics (15, 45). Thus, our approach averages out artifacts and other non-physiological factors (42).

Our corticospinal latency results are consistent with previous studies (7–10, 13). However, they must be interpreted with caution. The convolutive mixing models used to represent the EMG and EEG recordings, which are critical for accurate source separation (25, 27, 42), may introduce a temporal uncertainty to the reconstructed spike trains and tremor-related EEG components. We estimate this uncertainty to be about ±5 ms for each reconstructed source. Moreover, the propagation of the motor unit action potentials along the muscle fibers from the innervation zone to the uptake electrodes may introduce additional few ms delay. This could potentially further decrease the accuracy of EEG-CST delay estimation. In the current study, we used arrays of several tens of surface electrodes, whereas many previous studies were based on bipolar EMG recordings. The propagation of the motor unit action potentials may differ substantially across these two setups, and may also be muscle specific. Thus, the delays estimated in our study cannot easily be compared to the ones in other studies.

The availability of our method to automatically assess the accuracy with which each motor unit spike train is identified is also of critical importance because this accuracy is then reflected in the extracted tremor-related EEG component (see **Appendix**). Our group demonstrated in Holobar et al. (42) that motor units with PNR > 30 dB exhibit accuracy > 90% in identification of their firing patterns and in this study, we carefully utilized this knowledge to increase the accuracy of EEG component identification. In the future, we will investigate the minimal number of EEG channels required for accurate detection of tremor-related EEG activity, since it is likely that not all the EEG channels included in this study contribute significantly to tremor identification.

The proposed method is also computationally efficient. The most time consuming step is its first stage (surface EMG decomposition), which typically requires a few minutes of processing time on regular PC for 30 s long measurements. EEG decomposition in Equations (10, 11) is performed quickly.

The method does require multichannel EMG recordings from a muscle, increasing the experimental costs. However, multichannel EMG acquisition demonstrated significant progress in the recent years and became an important source of information in neurophysiology, neurology, sport sciences, prosthetics and ergonomics, to name just a few major scientific fields. Thus, it is likely that the price of multichannel acquisition systems will decrease in the near future.

We limited our study to the EEG decomposition of pathological tremor. The latter is a specific neurological disorder that is characterized by clear spectral peaks in acquired EEG, EMG, and inertial data. It is currently unclear to what extent the presented methodology is applicable to investigations of other types of pathological tremor (e.g., dystonic or cerebellar tremor) or to other disorders, such as multiple sclerosis, stroke and traumatic brain injuries and overactive thyroid, especially as tremor frequently accompanies these disorders. All these questions need to be systematically addressed in separate studies.

In conclusion, we have presented a novel method for estimating tremor-related cortical activity. This method uses pooled motor unit firings to directly extract the tremor component from cortical recordings. Based on the presented results, we believe that our method constitutes a significant step forward in the current state-of-the-art as: (a) it is the first method that directly extracts the tremor component from EEG recordings; (b) it successfully tracks time changes in the tremorrelated cortical activity and has a potential for online tremor detection.

### AUTHOR CONTRIBUTIONS

AH, JG, ER, JR, JB-L, JP, and VG participated in conceptualization and data acquisition. AH, JK, and VG participated in formal analysis and methodology. AH, ER, JB-L and JP participated in project administration and resources. All authors contributed to writing the original draft and revised the manuscript. All authors read and approved the final version of the manuscript and agreed for all aspects of the work.

### FUNDING

This study was supported by the Commission of the European Union, within Framework 7, under Grant Agreement number ICT-2011.5.1-287739 NeuroTREMOR: A novel concept for support to diagnosis and remote management of tremor and by Slovenian research agency (contracts J2-7357 and P2-0041). JG received funding from the Commission of the European Union (FP7-PEOPLE-2013-IOF-627384). ER received funding from the Ministry of Economy and Competitiveness (grant ESSENTIAL, DPI2015-72638-EXP). JR was supported by the Ministry of Economy and Competitiveness (grant DPI2015-68664-C4-3- R (MINECO/FEDER), NeuroMOD Therapy development and evaluation of motor and cognitive impact for Parkinson's disease rehabilitation). JP received funding from the Commission of the European Union (EXTEND project, 779982).

#### REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Holobar, Gallego, Kranjec, Rocon, Romero, Benito-León, Pons and Glaser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

#### APPENDIX

#### Estimation of csj(n),s in the Case of Phase-locked Motor Unit Firings

Assume that tremor-related EEG component in its analytic form can be expressed as a complex exponential function:

$$s\_{\vec{f}}(n) = e^{i\left(2\pi \frac{f\vec{f}}{f\_s}n + \varphi\_{\vec{f}}\right)}\tag{A1}$$

Where i is imaginary unit, f<sup>j</sup> and φ<sup>j</sup> are the frequency and the phase of s<sup>j</sup> , respectively, and f <sup>s</sup> is the sampling frequency. Then the cross-correlation **c**s<sup>j</sup> ,sλ (d) can be approximated by the following sample mean:

$$\begin{split} \mathfrak{c}\_{\boldsymbol{\varsigma},\boldsymbol{\varsigma},\boldsymbol{\varsigma}}(d) &= \frac{1}{N} \sum\_{n=0}^{N-1} s\_{\boldsymbol{f}}(n) \, s\_{\boldsymbol{\lambda}}^{H} \, \Big( n + d \Big) \\ &= \frac{1}{N} \sum\_{n=0}^{N-1} e^{i \left(2\pi \frac{f\_{\boldsymbol{f}}}{f\_{\boldsymbol{s}}} n + \varphi\_{\boldsymbol{\jmath}}\right)} e^{-i \left(2\pi \frac{f\_{\boldsymbol{\lambda}}}{f\_{\boldsymbol{s}}} (n+d) + \varphi\_{\boldsymbol{\lambda}}\right)} \\ &= \frac{1}{N} e^{i \left(\varphi\_{\boldsymbol{j}} - \varphi\_{\boldsymbol{\lambda}} - 2\pi \frac{f\_{\boldsymbol{\lambda}}}{f\_{\boldsymbol{s}}} d\right)} \sum\_{n=0}^{N-1} e^{i 2\pi n \frac{f\_{\boldsymbol{j}} - f\_{\boldsymbol{\lambda}}}{f\_{\boldsymbol{s}}}} \end{split} \tag{A2}$$

As well-known from Fourier analysis, the sum over N in (A2) converges to zero when N → ∞ and fj6= f <sup>λ</sup> . When fj= f <sup>λ</sup> (A2) yields

$$\mathbf{c}\_{s\_j s\_j}(d) = e^{-i2\pi \frac{f\_j}{f\_s}d} \tag{A3}$$

Define now the r-th motor unit spike train as

$$t\_r(n) = \sum\_{k=0}^{K\_r - 1} \delta \left( n - k \frac{f\_s}{f\_r} - d\_r - \Delta d\_{rk} \right) \tag{A4}$$

where d<sup>r</sup> is the common time lag (phase) of the pulses in tr(n), 1drk is the k-th realization of Gaussian random variable 1drk ∼ N(0, σ1d<sup>r</sup> ), f<sup>r</sup> is frequency of motor unit firings, and K<sup>r</sup> is the number of firings in the observed time interval. Then

$$\begin{split} \left| \mathfrak{c}\_{t,s\_{\mathbb{J}}} \left( d \right) \right| &= \frac{1}{N} \sum\_{k=0}^{K\_{\tau}-1} e^{-i \left( 2\pi \frac{f\_{j}}{f\_{s}} \left( k \frac{f\_{s}}{f\_{r}} + d\_{r} + \Delta d\_{rk} + d \right) + \varphi\_{\mathbb{j}} \right)} \\ &= \frac{1}{N} e^{-i \left( 2\pi \frac{f\_{j}}{f\_{s}} \left( d\_{r} + d \right) + \varphi\_{\mathbb{j}} \right)} \sum\_{k=0}^{K\_{\tau}-1} e^{-i 2\pi k \frac{f\_{j}}{f\_{r}}} \frac{f\_{j}}{c} \cdot e^{-i 2\pi \frac{f\_{j}}{f\_{s}} \Delta d\_{rk}} \\ \end{split} \tag{A5}$$

When EEG tremor is coupled to motor unit firings, we have fj

$$\{f\_{\vec{f}} = af\_r \text{ for } a \in \mathbb{Z}. \text{ Thus, } e^{-i2\pi k \frac{f\_r}{f\_r}} = 1, \text{ } \forall k, \text{ and (A4) simplifies to}\}$$

$$\mathbf{c}\_{t\_r,\mathbf{\tilde{y}}}\left(d\right) = \left(\frac{1}{N}e^{-i\left(2\pi\frac{f\_j}{f\_s}d\_r + \varphi\_j\right)}\right)\left(\sum\_{k=0}^{K\_r-1} e^{-i2\pi\frac{f\_j}{f\_s}\Delta d\_{rk}}\right)e^{-i2\pi\frac{f\_j}{f\_s}d} \text{ (A6)}$$

Since 1drk ∼ N(0, σ1d<sup>r</sup> ) and σ1d<sup>r</sup> is relatively small with respect to f<sup>s</sup> , the imaginary component of the sum in (A6) converges to fj

zero, so that PKr−<sup>1</sup> k=0 e −i2π fs 1drk = β≤ K<sup>r</sup> when K<sup>r</sup> → ∞, where β is a real number. When EEG tremor is not coupled to motor unit firings, i.e., f<sup>j</sup> 6= af<sup>r</sup> for a ǫ Z, (A5) converges to zero when K<sup>r</sup> → ∞. Therefore, for sufficiently large K<sup>r</sup> (A5) and (A6) yield:

$$\mathfrak{c}\_{t,\mathbf{s}}\left(d\right) \approx \alpha \mathfrak{c}\_{\mathfrak{s}\_{j,\mathbf{s}}}\left(d\right) \tag{A7}$$

with α = β N e −i 2<sup>π</sup> fj fs dr+ϕ<sup>j</sup> . In practice, The higher the number

of motor unit firings K<sup>r</sup> , the better the estimate (A6). We can increase K<sup>r</sup> by increasing the length of the signal's time support, but this comes at the cost of long signal acquisitions. In the case of pathological tremor, motor unit firing patterns are highly synchronized and active motor units share approximately the same tremor-related firing rate f<sup>r</sup> = fT, ∀r and initial delays d<sup>r</sup> = dT, ∀r (see **Figure 4**). In such case, the CST of R motor units can be modeled by

$$\text{CST}(n) = \sum\_{r=1}^{R} t\_r \ (n) = \sum\_{r=1}^{R} \sum\_{k=0}^{K\_r - 1} \delta \left( n - k \frac{f\_s}{f\_T} - d\_T - \Delta d\_{rk} \right) \text{(A8)}$$

This increases the number of motor unit firings in (A5) to

$$\mathbf{c}\_{\text{CST},\mathbf{s}\_{\text{j}}}(d) = \left(\frac{1}{T}e^{-i\left(2\pi\frac{f\_{j}}{f\_{s}}d\_{\text{r}}+\varphi\_{\text{j}}\right)}\right)\left(\sum\_{r=1}^{R}\sum\_{k=0}^{K\_{r}-1}e^{-i2\pi\frac{f\_{j}}{f\_{s}}\Delta d\_{rk}}\right) \tag{49}$$

$$e^{-i2\pi\frac{f\_{\text{j}}}{f\_{\text{s}}}d} \tag{40}$$

and, thus, improves the estimate (A6).

In noise-free conditions (A6), (A8) and (7) fully justify the approximation (9). In the presence of noise, however, further analytical derivations of the noise residual are required to verify whether the proposed estimator (10) is truly an LMMSE estimator and, thus, Bayesian optimal. Nevertheless, the results on both synthetic and experimental signals confirm the effectiveness of the proposed tremor estimation.

# How Physical Activities Affect Mental Fatigue Based on EEG Energy, Connectivity, and Complexity

Rui Xu1,2, Chuncui Zhang1,2, Feng He1,2 \*, Xin Zhao1,2, Hongzhi Qi 1,2, Peng Zhou1,2 , Lixin Zhang1,2 and Dong Ming1,2 \*

*<sup>1</sup> Lab of Neural Engineering & Rehabilitation, Department of Biomedical Engineering, College of Precision Instruments and Optoelectronics Engineering, Tianjin University, Tianjin, China, <sup>2</sup> Tianjin International Joint Research Center for Neural Engineering, Academy of Medical Engineering and Translational Medicine, Tianjin University, Tianjin, China*

#### Edited by:

*Xu Zhang, University of Science and Technology of China, China*

#### Reviewed by:

*Zehong Cao, University of Technology Sydney, Australia Chin-Teng Lin, National Chiao Tung University, Taiwan Bo Yao, Chinese Academy of Medical Sciences and Peking Union Medical College, China*

\*Correspondence:

*Dong Ming richardming@tju.edu.cn Feng He heaven@tju.edu.cn*

#### Specialty section:

*This article was submitted to Neurotrauma, a section of the journal Frontiers in Neurology*

Received: *30 June 2018* Accepted: *09 October 2018* Published: *31 October 2018*

#### Citation:

*Xu R, Zhang C, He F, Zhao X, Qi H, Zhou P, Zhang L and Ming D (2018) How Physical Activities Affect Mental Fatigue Based on EEG Energy, Connectivity, and Complexity. Front. Neurol. 9:915. doi: 10.3389/fneur.2018.00915* Many studies have verified that there is an interaction between physical activities and mental fatigue. However, few studies are focused on the effect of physical activities on mental fatigue. This study was to analyze the states of mental fatigue based on electroencephalography (EEG) and investigate how physical activities affect mental fatigue. Fourteen healthy participants participated in an experiment including a 2-back mental task (the control) and the same mental task with cycling simultaneously (physical-mental task). Each experiment consisted of three 20 min fatigue-inducing sessions repeatedly (mental fatigue for mental tasks or mental fatigue plus physical activities for physical-mental tasks). During the evaluation sessions (before and after the fatigue-inducing sessions), the states of the participants were assessed by EEG parameters. Wavelet Packet Energy (WPE), Spectral Coherence Value (SCV), and Lempel-Ziv Complexity (LZC) were used to indicate mental fatigue from the perspectives of activation, functional connectivity, and complexity of the brain. The indices are the beta band energy *E*β, the energy ratio *E*α/β, inter-hemispheric SCV of beta band *SCV*<sup>β</sup> and LZC. The statistical analysis shows that mental fatigue was detected by *E*β, *E*α/β, *SCV*β, and LZC in physical-mental task. The slopes of the linear fit on these indices verified that the mental fatigue increased more fast during physical-mental task. It is concluded form the result that physical activities can enhance the mental fatigue and speed up the fatigue process based on brain activation, functional connection, and complexity. This result differs from the traditional opinion that physical activities have no influence on mental fatigue, and finds that physical activities can increase mental fatigue. This finding helps fatigue management through exercise instruction.

Keywords: coherence, complexity, EEG, mental fatigue, physical activities

### INTRODUCTION

Mental fatigue is defined as a "state of reduced mental alertness that impairs performance" (1, 2). It is believed to exist in the nervous system and affect the mental activities of people, such as the motivation and attention (3). It often occurs when working on cognitively demanding tasks for a prolonged period of time (3, 4) and causes difficulties for people in maintaining task performance at an adequate level (5). Furthermore, in some cases, mental fatigue can lead to vital consequences for

**77**

people, such as drivers, pilots and surgeons. A naturalistic driving study found that mental fatigue causes 12% of car crashes and 10% of near-misses (6, 7). Sleepiness is usually tied to mental fatigue. Another study of French drivers showed that 46% of the drivers experiencing near-misses reported related sleepiness and 5.2% of the car accidents were caused by sleepiness (8). Fatigue can affect the drivers' performance regarding EEG and functional near-infrared spectroscopy (fNIRS) (9, 10). Mental fatigue is also a contributing factor for some medical conditions, such as cardiovascular diseases (11), hypothyroidism (12), and fibromyalgia (13). Therefore, it is of great importance to find out the mechanism (3, 14) and the assessment (15, 16) of mental fatigue. However, the level of mental fatigue is difficult to identify. Usually, mental fatigue is detected by significant change of fatigue indices.

Physical activities usually come with physical fatigue. Physical fatigue (or muscle fatigue) is another kind of fatigue which is caused by physical activities and defined as the inability to maintain a required force level after prolonged use of muscle (17). It is a complex, multifactorial phenomenon influenced by the characteristics of the task being performed (18). Physical fatigue is believed to be developed gradually soon after the onset of the sustained physical activity (19). The common protocol to quantify muscle fatigue is to interrupt the fatiguing session with short maximal contraction to estimate the decline in the maximal fore capacity. Additionally, EMG may be used as an indicator of muscle fatigue.

The studies on mental fatigue have been limited by the cognitive tasks (3, 4, 14, 20), such as N-back tasks, serial-7 subtraction arithmetic tasks, Wisconsin Card Sorting Test, and forward digit span. These tasks are all with high-intensity mental activities, which can effectively induce mental fatigue. Additionally, there are almost no physical activities except for the necessary responses of the cognitive tasks. It seems that an isolation exists between mental fatigue and physical activities. A recent study used simulated driving task to induce mental fatigue, and classified mental states based on generalized partial directed coherence of EEG (21). There are only few studies on the interaction between mental fatigue and physical activities. It is verified that mental fatigue decreased exercise tolerance through higher perception of effort (22). Tanaka et al. investigated the effect of mental fatigue on physical activities using alpha-band event-related synchronization (ERD) of magnetoencephalography (MEG) and found that the mental fatigue suppresses activities in the right anterior cingulate cortex during physical fatigue (23). Similar results were also found by Mehta et al. that mental fatigue impaired motor performance and muscle capacity when exploring the prefrontal cortex activation with fNIRS (17). Additionally, Simth et al. analyzed the time-motion data of mentally-fatigued athletes, and observed a reduction in low-intensity activity velocity, which was deduced to be mediated by an increased perception of effort rather than cardiovascular or metabolic mechanisms (24). Therefore, it is commonly accepted that mental fatigue impairs physical performance, probably by increasing the effort perception. However, how the physical activities affect the mental fatigue is still unclear. Finding out the effect of physical activities on mental fatigue will help deeply understand mental fatigue, and on the other hand, it will benefit the prevention of mental fatigue.

EEG is a promising method to estimate mental fatigue (25). EEG energy has been proposed to be a valid and reliable indicator of mental fatigue (26). The EEG energy ratio was proposed in order to consider the variation of EEG energy in more than one frequency band and it can be used as a fatigue indicator during driving (27, 28). The energy ratio of α/β was believed to be a more reliable fatigue detection index than the energy index, since it showed a clear fatigue-increasing process as the ratio between the slow and fast wave activities increased (28, 29). More recently, α/β was also used to detect the fatigue variability between watching 2D and 3D TV (30). Another important parameter is the spectral coherence value (SCV) of EEG. The SCV is often used to investigate the functional connectivity of two signals. It has been regarded as an effective indicator of mental fatigue recently. The inter-hemispheric SCV of beta band decreases after

participant during the experiment.

a long-term cognitive task (31), which may indicate the weakened cooperation of these brain regions during mental fatigue (32). Mental fatigue is also reported to decrease brain complexity (33, 34). The decrease of brain complexity can be interpreted as the decrease of brain's capability to continue a task. Brain complexity can be assessed by many parameters or methods, such as the inherent fuzzy entropy method (35–37) and Lempel Ziv Complexity (LZC). LZC was first proposed by Lempel and Ziv to assess the system's complexity (38). It has been widely used to identify the complexity in EEG (39). It was noted that LZC was more sensitive than the conventional spectral parameters of EEG to reflect the mental activity (40), which distinguished the schizophrenia, the depression and the healthy controls. More recently, this method was also used to recognize the poststroke depression (41). Therefore, LZC can be used to detect the brain complexity and should be able to reflect the mental states.

Cycling is a very common exercise in daily life. Cyclingbased movement can provide a safe and effective way for walking training and lower limb coordination training (42). Usually cycling studies related to fatigue are limited to peripheral fatigue, estimated by Electrocardiogram (ECG) and Electromyogram (EMG) parameters (43, 44). In the present study, the mental fatigue was generated by the n-back tasks. The n-back task is a working memory task, which activate higher-order cognition, such as language, reasoning and problem-solving (45). The reliability and validity of n-back task to cause mental fatigue have been confirmed (14, 23, 46). Physical activities (cycling) were added with the 2-back task to generate mental fatigue during activities.

Most studies focused on mental fatigue only, neglecting the effect of physical activities. This interaction will help us control fatigue and diseases with fatigue syndromes. The purpose of this study is to investigate how physical activities change mental fatigue based on EEG and to compare the sensitivity of different mental fatigue indicators. As physical fatigue increases, the nervous system will make more efforts to maintain the motor task. Therefore, we expect that physical activities will strengthen mental fatigue by occupying more neural resources. This paper is organized as follows: In section Materials and Methods, the experimental design is presented, and data processing methods on mental fatigue detection are introduced. section Results shows the results of the calculation and analysis. Then the obtained results are discussed in section Discussion. Conclusions are drawn in section Conclusions.

### MATERIALS AND METHODS

#### Participants

Fourteen healthy participants (9 females and 5 males; mean age: 22.4 ± 1.6) without any chronic fatigue syndrome or motor


dysfunction, participated in the experiment. The experiment protocol was approved by the Ethics Committee of Tianjin University and the participants have signed a consent form before experiment.

### Experimental Design

The experiment was undertaken in the EEG lab of Tianjin University. It consisted of a mental task and a physical mental task as in **Figure 1A**. There are 2 min rest sessions before and after the fatigue-generated sessions for data acquisition. During the mental task, the participant was seated in a chair in front of the computer. The mental task included three 2-back sessions based on 26 upper-case letters. The letters occurred on the screen randomly, which was generated by Psychtoolbox within Matlab. During this task, each letter was presented for 0.5 s at the center of the screen every 3 s as shown in **Figure 1B**. Participants had to judge whether the presenting letter was the same as the one that had appeared two presentations before. If it was the same (target stimulus), they were to press the left button with their right index finger; if it was different (non-target stimulus), they were to press the right button with their right middle finger. The participants were instructed to perform the task trials as quickly and as correctly as possible during the show of this letter. The participants were trained with the 2-back session before the experiment in order to fully and correctly understand

the experiment. The 2-back session lasted for 20 min and was repeated for 3 times. The ratio of the numbers of target stimulus and non-target stimulus was 1:2.

A cycling part was added to the mental task described before to form the physical-mental task. A spinning bike (Yowza B300, Yowza Fitness, US) with a power measure was used for cycling. The maximum cycling power of every participant was measured before experiment: the participant first cycled for 2 min at 30 W and the power was increased at 10 W/min until the participant failed to maintain cycling (when the spinning speed is <60 rounds/min). This power was the maximum cycling power for the participant. The experimental power (shown in **Table 1**) was chosen as 40% of the maximum power, which is feasible for the participants to maintain for the whole course of the experiment. The participant cycled at the experimental power and completed the 2-back session simultaneously. The physical-mental task lasted for 20 min and was repeated 3 times.

The participants refrained from intense mental and physical activities, consumed a normal diet and beverages (excluding caffeinated beverages), and maintained normal sleeping hours on the day before the experiment. The mental task and physical-mental task were separated by 3 days to exclude the cross interference and the order of the two tasks were randomly arranged.

### Data Acquisition

The resting EEG data were recorded four times during the mental or physical mental task when the participant (with eyes open) was seated without any obvious mental or physical activities as in **Figure 1C**: one (t0) before the task and the other three (t1, t2, and t3) after each fatigue-inducing session (shown in **Figure 1A**). In this way, the mental fatigue caused by mental or physical-mental tasks remained and there was no EMG disturbance in t0, t1, t2, and t3.

EEG data were collected with a Neuroscan SynAmps2 amplifier (sampling rate: 1,000 Hz). The electrodes were placed on the scalp according to the extension of the international 10–20 electrode positioning system (47) with the reference at right mastoid. Eye movements and blinks were monitored by recording the horizontal and vertical Electrooculogram (EOG) with two bipolar pairs of electrodes. The EEG data in F3, F4, FZ, C3, C4, CZ, P3, P4, PZ, T3, T4, T5, T6, O1, O2, and OZ

were analyzed in this study. These channels were selected as the representing channels from the frontal, central, parietal, temporal and occipital individually. The channels of F3, F4, FZ, C3, C4, CZ, P3, P4, and PZ are determined based on a coherence analysis during mental fatigue (31), and T3, T4, T5, T6, O1, O2, and OZ are supplemented for temporal and occipital areas.

### Data Processing

The purpose of the pre-processing was to obtain clear EEG data and to increase the computing speed of feature extraction. The EEG data were re-referenced to bilateral mastoids, downsampled to 500 Hz and filtered at 0.5∼45 Hz with a 4th-order Butterworth zero-phase digital filter. The EOG interference was removed using Independent Component Analysis (EEGLAB) (48). As the EOG artifacts are in larger amplitude than pure EEG and separated by several seconds randomly, the ICA component with these artifacts can be picked out and deleted using EEGLAB. In order to obtain the steady EEG data during t0, t1, t2, and t3, the 1 min data starting from 30 s after the onset of resting EEG recording were extracted and analyzed.

The energy, interhemispheric SCV and complexity features were then calculated based on the processed EEG data.

#### Wavelet Packet Energy

Wavelet Packet Decomposition (WPD) is generalized from the Wavelet Decomposition (WD). The advantage of WPD is that


*The values in bold are significantly different.*

both the detail and approximation coefficients are decomposed (49), that is, precise frequency information is obtained for high frequencies. Daubechies ("db10") was used as the mother wavelet in this study. The 60 s resting EEG were decomposed by an 8-level WPD (28 ≈ fs/2, fs = 500 Hz).

After WPD, the wavelet packet coefficients (WPC) of 256 frequency bins were obtained. The WPE across frequencies were calculated as:

$$WPE\left(f\_i\right) = \|WPC\_i\|\_2, i = \|0, 1, \ldots, 255\tag{1}$$

where || ||<sup>2</sup> is 2-norm computation, i is the node number of the 8th level, f<sup>i</sup> is the corresponding frequency of the ith node.

$$f\_i = \frac{i+1}{256} \cdot \frac{f\_s}{2} \tag{2}$$

The power of EEG data is subdivided into four frequency bands: delta (0.5∼4 Hz), theta (4∼8 Hz), alpha (8∼16 Hz), and beta (16∼32 Hz) bands. The EEG energies and energy ratios of different frequency bands are important in indicating mental fatigue (28, 29). The energy ratio is more reliable than the band energy, and it considers the energy variation of more than one band. As the energy of delta and theta rhythms is not very sensitive to mental fatigue (50), the relative energy of resting EEG in beta band E<sup>β</sup> and the energy ratio of alpha and beta bands Eα/β were calculated as in Equations 3, 4.

$$E\_{\beta} = \frac{\sum\_{j=16}^{32} WPE\left(f\_j\right)}{\sum\_{j=1}^{46} WPE\left(f\_j\right)}\tag{3}$$

$$E\_{\partial/\beta} = \frac{\sum\_{j=8}^{16} WPE\left(f\_j\right)}{\sum\_{j=16}^{32} WPE\left(f\_j\right)}\tag{4}$$

#### Spectral Coherence Value

The SCV of EEG can be used to estimate the relationship between two channels of EEG at any given frequency. If x and y are EEG data of two different channels, the SCV of x and y is estimated as

$$\text{SCV}\_{\text{x,y}}\left(f\right) = \frac{\left|\text{S}\_{\text{xy}}\left(f\right)\right|^2}{\text{S}\_{\text{xx}}\left(f\right) \cdot \text{S}\_{\text{yy}}\left(f\right)}\tag{5}$$

where Sxx and Syy are the power spectral densities of x and y and Sxy the cross spectrum of x and y. The 60 s data was segmented into 59 segments by a Hamming window of 2 s with an overlap of 1 s. The cross- and auto-spectrum were obtained by the average spectrum of these 59 segments.

In order to estimate the interhemispheric functional connectivity, the SCV of EEG signals in F3-F4, C3-C4, P3-P4, T3-T4, T5-T6, and O1-O2 electrode pairs were calculated. The beta band SCV was obtained by Equation 6.

$$SCV\_{\beta} = \sum\_{f=16\text{Hz}}^{32\text{Hz}} SCV\left(f\right) \tag{6}$$

#### Lempel-Ziv Complexity

The detailed method is described in the work of Nagarajan et al. (51). The complexity c(N) was normalized to N/log2N resulting in

$$
\lambda = c \left( N \right) / \frac{N}{\log\_2 N} \tag{7}
$$

The 60 s data was segmented into 20 segments. For each segment, the number of points N = 1,500.

For each segment, the data was binarised by comparison of each data point x(i) (i = 1, 2, . . . , N) with its median value M<sup>d</sup> in the following way:

$$s(i) = \begin{cases} 0, \ x(i) \le M\_d \\ 1, \ x(i) > M\_d \end{cases} \tag{8}$$

The normalized complexity λ k j for segment k channel j was obtained according to Equation 7. The LZC value λ<sup>j</sup> for each channel j was the average of λ k j of 20 segments.

#### Data Analysis

The energy parameters of E<sup>β</sup> and Eα/β, interhemispheric beta band SCV and LZC were compared to indicate mental fatigue. The signed rank test was used to detect the significant difference between the states before and after mental or physical-mental task in each channel or channel pair.

In order to estimate the varying trend of the parameters, the linear fit of the parameters of channels or electrode pairs with

significant difference was performed to fit a linear polynomial curve. These parameters are the sum of the average values in the channels with significance, i.e., E<sup>β</sup> in C3, P3, PZ, T3, T5, and OZ, Eα/β in C3, P3, and T4, SCV<sup>β</sup> in P3-P4, and LZC in

TABLE 3 | Eα/β for mental and physical mental tasks.


*The values in bold are significantly different.*

F8, P3, T3, and T5. As the variation of each channel is also important to provide more accurate information, these indices for each channel with significance are also calculated over time of the tasks.

### RESULTS

#### Data Processing

As the EOG artifacts are very strong in FP1 and FP2 (they are the nearest electrode positions to the eyes), the EEG in FP1 of one participant is shown in **Figure 2** to illustrate the efficacy of the filtering in our study. The red rectangle indicates where the EOG artifact occurs. **Figure 2B** shows that the common filtering can only remove the noise in certain frequencies, and it can not wipe off EOG artifacts. However, EOG artifacts can be removed successfully by ICA as shown in **Figure 2C**.

#### EEG Energy

The average of the beta band energy for all the channels of 14 participants in mental and physical-mental tasks is shown in **Figure 3**, **Table 2**. For the mental task, there is no significant difference in E<sup>β</sup> **Figure 3A**, but in **Figure 3B**, E<sup>β</sup> decreased significantly after physical mental task in Central (C3, p = 0.017), Parietal (P3, p = 0.025, and PZ, p = 0.030), Temporal (T3, p = 0.020, and T5, p = 0.007), and Occipital (OZ, p = 0.042) areas. **Figures 4**, **Table 3** show the energy ratio of alpha and beta bands Eα/β for all the channels. Significant increase of Eα/β only occurs in Central (C3, p = 0.025), Parietal (P3, p = 0.049), and Temporal (T4, p = 0.035) areas after physical mental task (**Figure 4B**).

In view of the average values, the E<sup>β</sup> in all the channels has an obvious decrease after physical mental task (**Figure 3B**), while only E<sup>β</sup> in 8 channels decreases slightly in mental tasks (**Figure 3A**). Eα/β has a similar variability trend that Eα/β increases obviously in most channels after physical mental task (**Figure 4B**).

#### Interhemispheric SCV

In order to estimate the interhemispheric connectivity of the brain, we computed the SCV<sup>β</sup> between the two symmetric channels. The SCV<sup>β</sup> values remains similar before and after mental task (**Figure 5A**, **Table 4**). However, there is a significant decrease in Parietal (P3-P4, p = 0.042) area after physical mental task (**Figure 5B**). The interhemispheric connectivity decreased significantly after physical mental task, but no significance was found after mental task.

### Lemple-Ziv Complexity

The mean LZC of all the participants for all the channels are shown in **Figure 6**, **Table 5**. There is no significant decrease between LZC before and after the mental task, while after physical-mental task there exist significant decreases in Frontal (F8, p = 0.049), Parietal (P3, p = 0.042), and Temporal (T3, p = 0.019, and T5, p = 0.035) areas. To be noted that, the LZC average values before and after mental task are quite similar in **Figure 6A**, but the values decrease after physical mental task in all the channels in **Figure 6B**. This indicates that the brain complexity was significantly influenced by the physical mental task.

#### Linear Fit

During the experiment, the resting EEG were collected for t0, t1, t2, and t3 of the tasks. The features of these periods are shown

TABLE 4 | *SCV*β for mental and physical mental tasks.


*The values in bold are significantly different.*

and fitted in **Figure 7**. The slopes of these features indicate the variability rate of these features when the task is proceeding. It is shown that the slopes of E<sup>β</sup> (**Figure 7A**), SCV<sup>β</sup> (**Figure 7C**), and LZC (**Figure 7D**) in physical mental task are smaller than those in mental tasks, while the slope of Eα/β (**Figure 7B**) in physical mental task is larger than that in mental task. This shows that the variability of these features caused by physical mental task is larger than that caused by mental task.

The indices in each channel with significance are shown in **Figure 8** (A, C, E, and G for mental tasks, and B, D, F, and H for physical mental tasks). It has been indicated that the significant change between values in t0 and t3 only occurs in physical mental tasks. This result shows that all the channels varies in a similar trend, and there is a larger drop or increase in t1 than t3.

### DISCUSSION

This study investigated the effect of physical activities on mental fatigue through specifically designed experiments and different EEG parameters. These EEG features, including beta



*The values in bold are significantly different.*

band energy, energy ratio, SCV, and LZC, estimate the fatigue states from the perspectives of brain activation, interhemispheric connectivity and brain complexity. From our results, the mental fatigue causes a significant index change in physical mental task, and physical activities speed up the fatigue process. This result reveals the effect of physical activities on mental fatigue, which differs from the traditional opinion that physical activities have no influence on mental fatigue, and help instruct exercise for people with fatigue.

#### Indices of Mental Fatigue

The E<sup>β</sup> and Eα/β calculated in this study has been previously used to measure mental fatigue (28, 30). The significant increase of these indices was found with mental fatigue in this study, which confirmed the findings by Eoh et al. (28) and Chen et al. (30). Jap et al. had paid more attention on energy ratios as they combined the power of different frequency bands (29) and provided a measure for greater magnitude of changes in brain activity throughout driving (52). The significant decreases were shown in physical-mental task. It was convinced that the mental fatigue occurred in both tasks, but these parameters were not so sensitive that there was no significant difference in mental task.

Zhang et al. estimated the cortical functional connectivity in Frontal, Central, and Parietal regions during mental fatigue and found that the SCV<sup>β</sup> decreased in Parietal regions (31). Similar results were also found in the present study (**Figure 5**). Beta band EEG is the main EEG wave reflecting excitatory state of the brain cortex, which is associated with increased alertness and excitement. Therefore, the decreased SCV<sup>β</sup> of beta band are related to mental fatigue. The significant decreases occurred in P3-P4 electrode pair after physical-mental task. It was deduced that the mental fatigue after physical mental task caused a significant decrease in SCVβ.

Brain complexity was validated to decrease as the mental fatigue level increased (33). The decrease of brain complexity may be explained by neurons' functional isolation with greater autonomy of brain components (53). In this study, LZC decreased significantly in Frontal, Parietal, and Temporal regions in the physical-mental task. However, there was no significant changes of LZC in mental task. It seemed that LZC was very sensitive to the brain complexity variation caused by physical activities and this variation was distributed almost in the whole scalp.

SCV is commonly used to characterize the synchronization and functional coupling of two signals. A study has provided evidence that it is an effective and reliable way to quantify brain response to mental fatigue (31). Additionally, the complexity is another perspective to estimate mental state. LZC has excellent performance in analysis of mental disorders (41), oppositely it should be sensitive to mental fatigue when the brain is less activated.

### Interaction of Mental and Physical Activities

It is well known that mental fatigue has an impact on cognitive performances (3, 5) and physical performances (22, 23), and it even causes muscle fatigue (54). However, to the authors' knowledge, how physical activities affect mental fatigue has not been thoroughly investigated. The research that EEG-EMG coherence can predict muscle fatigue (55) demonstrates that muscle fatigue affects the brain and muscle activities at the same time. Mashiko et al. found that playing a rugby football match can cause mental fatigue (56), which is consistent to the result of the present study. However, the method they used for the evaluation of mental fatigue was the Profile of Mood State (POMS) test. An interesting study investigated the brain activities during exercise in different temperatures (57). The main finding is the hyperthermia-associated mental fatigue assessed by the shift in EEG power distribution. However, there was no pure physical task in the study as the control.

This present study firstly investigated the mental fatigue in mental and physical-mental tasks with Eβ, Eα/β, SCVβ, and LZC, and estimated the variation of the parameters during the task using linear fit. The result indicated that the physical activities (cycling) are able to produce mental fatigue, causing significant differences in Eβ, Eα/β, SCVβ, and LZC. The result of the linear fit in **Figure 7** demonstrates that both mental and physical-mental tasks can increase Eα/β and decrease Eβ, SCV<sup>β</sup> and LZC, but the variation is more rapid in physical-mental task. Therefore, the physical activities speed up the fatigue process. It is deduced that the control of the movement is a kind of mental activity that can cause mental fatigue. As with physical fatigue increased, the participant should make more effort to complete the task. More attention was naturally paid on cycling. The task with attention

is often used to generate mental fatigue, and attention is highly related to mental fatigue (58). Therefore, the significant changes of these indices may be related to subjective effort and attention for physical activities.

An interesting phenomenon is shown in **Figure 8**, where an obvious decrease or increase occurs in t1, and the variation is retarded in the following t2 and t3. It is deduced that as the task is proceeding, the participant becomes familiar with this task, and finds an easier way to complete the task with fewer neurons joining in. Therefore, the mental fatigue increases slowly in t2 and t3.

### Limitations and Future Work

Although some valuable findings are obtained, there are still several limitations in this study. The study only compared the parameters in mental and physical-mental tasks. In order to determine the effect of physical activities only, a physical task without 2-back task should be analyzed in future studies. EEG connectivity was also estimated by isolated effective coherence (59, 60). The connectivity with direction may be a new attempt to explore the mental fatigue states.

## CONCLUSIONS

This study has investigated the energy parameter E<sup>β</sup> and Eα/β, connectivity parameter (SCVβ) and complexity parameter (LZC) to indicated fatigue in mental and physical-mental tasks. Different from the traditional view that mental fatigue is caused by mental tasks, the present study works on mental fatigue affected by physical activities. According to the statistical results, the participants are more fatigued after physical-mental task than after mental task. The linear fit results also show that physical activities speed up the fatigue process. Thus, the physical activities enhance the mental fatigue significantly. The result of this study provides a new perspective of the cause of mental fatigue. Further, this may help explain why the mental fatigue can impair physical performances: the mental fatigue leads to a decrease in the ability of the

motor control. Therefore, this study helps understand the mechanism of the interaction between mental fatigue and physical activities.

AUTHOR CONTRIBUTIONS

All authors contributed to conception and design of the study and were involved in drafting and critically revising the manuscript. Additionally, RX, CZ, and FH carried out the experimental work and data processing. RX prepared the first draft paper. XZ, HQ, PZ, and LZ provided interpretation of the experimental results. DM and RX worked up the draft paper into the final version. All authors gave final approval for publication.

### FUNDING

This work was supported by National Key Research and Development Program of China under grant (No. 2017YFB1300300), National Natural Science Foundation of China (No.81630051, 91520205, 61603269, 81601565, 81571762, and 31500865), and Tianjin Key Technology R&D Program (No. 15ZCZDSY00930).

#### REFERENCES


#### ACKNOWLEDGMENTS

The authors would like to thank all the participants who participated in the experiment.


headache prevention, and depression treatment. IEEE Access. (2017) 5:10612– 21. doi: 10.1109/ACCESS.2017.2675884


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Xu, Zhang, He, Zhao, Qi, Zhou, Zhang and Ming. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Muscle Fatigue Post-stroke Elicited From Kilohertz-Frequency Subthreshold Nerve Stimulation

Yang Zheng, Henry Shin and Xiaogang Hu\*

*Joint Department of Biomedical Engineering, University of North Carolina at Chapel Hill and North Carolina State University, Raleigh, NC, United States*

Purpose: Rapid muscle fatigue limits clinical applications of functional electrical stimulation (FES) for individuals with motor impairments. This study aimed to characterize the sustainability of muscle force elicited with a novel transcutaneous nerve stimulation technique.

#### Edited by:

*Yingchun Zhang, University of Houston, United States*

#### Reviewed by:

*Cliff Klein, Guangdong Province Work Injury Rehabilitation Hospital, China Zhen Yuan, University of Macau, China*

> \*Correspondence: *Xiaogang Hu xiaogang@unc.edu*

#### Specialty section:

*This article was submitted to Stroke, a section of the journal Frontiers in Neurology*

Received: *11 September 2018* Accepted: *21 November 2018* Published: *04 December 2018*

#### Citation:

*Zheng Y, Shin H and Hu X (2018) Muscle Fatigue Post-stroke Elicited From Kilohertz-Frequency Subthreshold Nerve Stimulation. Front. Neurol. 9:1061. doi: 10.3389/fneur.2018.01061* Method: A hemiplegic chronic stroke survivor was recruited in this case study. Clustered subthreshold pulses of 60-µs with kilohertz (12.5 kHz) carrier frequency (highfrequency mode, HF) were delivered transcutaneously to the proximal segment of the median/ulnar nerve bundles to evaluate the finger flexor muscle fatigue on both sides of the stroke survivor. Conventional nerve stimulation technique with 600-µs pulses at 30 Hz (low-frequency mode, LF) served as the control condition. Fatigue was evoked by intermittently delivering 3-s stimulation trains with 1-s resting. For fair comparison, initial contraction forces (approximately 30% of the maximal voluntary contraction) were matched between the HF and LF modes. Muscle fatigue was evaluated through elicited finger flexion forces (amplitude and fluctuation) and muscle activation patterns quantified by high-density electromyography (EMG).

Result: Compared with those from the LF stimuli, the elicited forces declined more slowly, and the force plateau was higher under the HF stimulation for both the affected and contralateral sides, resulting in a more sustainable force output at higher levels. Meanwhile, the force fluctuation under the HF stimulation increased more slowly, and, thus, was smaller after successive stimulation trains compared with the LF stimuli, indicating a less synchronized activation of muscle fibers. The efficiency of the muscle activation, measured as the force-EMG ratio, was also higher in the HF stimulation mode.

Conclusion: Our results indicated that the HF nerve stimulation technique can reduce muscle fatigue in stroke survivors by maintaining a higher efficiency of muscle activations compared with the LF stimulation. The technique can help improve the performance of neurorehabilitation methods based on electrical stimulation, and facilitate the utility of FES in clinical populations.

Keywords: transcutaneous electrical nerve stimulation, electromyography, muscle fatigue, stroke, kilohertz-frequency

### INTRODUCTION

Functional electrical stimulation (FES) can electrically activate muscle fibers with electrodes placed at the muscle belly, and can be used to promote functional improvement in paralyzed limbs due to neurological disorders (1, 2). Even though FES has shown some success in assisting individuals with neurological disorders, the clinical impact is hampered by a critical limitation which involves rapid muscle fatigue onset during repeated stimulations (3, 4). The increased muscle fatigability can arise from factors including the violation of the size principle of motor unit (MU) recruitment (5, 6) and the highly synchronized activation of muscle fibers. Furthermore, muscles from the paralyzed limbs are more fatigable compared with the muscles of the contralateral side or intact individuals due to central (7) or peripheral fatigue (8–10), which can further limit clinical applications of FES.

Previous studies have explored various techniques to reduce muscle fatigue. For example, prolonged force production has been observed by alternately stimulating muscle synergists (11) or using randomly modulated stimulation parameters to activate different muscle fibers separately (12–14). Similarly, multi-electrode stimulation techniques have been used to asynchronously activate different muscle portions (15–18). More recently, stimulation at the proximal segment of nerve bundles has been investigated (19, 20), in order to recruit a range of MUs, especially the more fatigue-resistant MUs. The nerve stimulation technique can result in less fatigable contractions (21) and more dexterous movements (22). In addition, the nerve bundle stimulation technique can activate afferent fibers with a low current amplitude and a high frequency (23), resulting in a more physiological recruitment order of MUs. However, the activations of different MUs are still highly synchronized, and the H-reflex activity is still secondary compared with the M-wave when large muscle forces are required.

In a previous study (24), clustered narrow pulses with a carrier frequency of kilohertz (high frequency, HF) targeting the proximal nerve bundles can lead to asynchronized activation of muscle fibers and reduced force fluctuations, compared with wide pulses delivered at a low frequency (LF). Furthermore, the HF stimulation technique can reduce muscle fatigue in healthy subjects (25). However, it is not known to what degree the HF stimulation technique can improve the force sustainability in paralyzed limbs of individuals with neurological disorders. Accordingly, a hemiplegic stroke survivor was recruited in the current case study. The elicited finger flexion forces and the EMG activity from the finger flexor muscles were compared between the HF and LF stimulation on both the affected (paretic) and contralateral sides.

### METHODS

#### Experiment and Subject Subject Background

A 54-year-old chronic stroke survivor with severe motor impairment in the arm and hand was recruited. The subject suffered from an ischemic stroke in the left corona radiate 2 years prior to the testing. The Chedoke-McMaster Stroke Assessment of the hand component was 2, and the arm component was 3. No cognition impairment was reported. Two experimental sessions, with one on each side, were performed with an interval of 1 month. This study was carried out in accordance with the recommendations of the Institutional Review Board (IRB) of the University of North Carolina at Chapel Hill with written informed consent from the subject. The subject gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the local IRB. Additional written informed consent was obtained from the subject for the publication of this case report.

#### Experiment Setup

The subject sat in an arm chair at a comfortable height with his back supported during the experiment. The hand was restrained, and the forearm was supported (**Figure 1A**). The individual fingers were secured to four miniature load cells (SM-200N, Interface), respectively for force measurements with a sample frequency of 1 kHz.

An 8 × 16 channel high density EMG electrode pad (OT Bioelettronica) was placed at the finger flexor muscles to obtain the EMG activities (**Figure 1A**). The electrode diameter was 3 mm with an inter-electrode distance of 10 mm. The center of the electrode grid was aligned with the midline between the olecranon process and the styloid process. Monopolar EMG signals were amplified using EMG-USB2+ (OT Bioelettronica). The sample rate, gain and bandwidth of the amplifier were set at 5,120, 200, and 10–900 Hz, respectively.

Sixteen (2 × 8 array) gel-based skin-surface stimulation electrodes (10 mm in diameter) were placed beneath the short head of the biceps brachii along the ulnar/median nerves. A custom-made MATLAB user interface was used to control a stimulator (STG4008, Multichannel Systems) and a switch matrix (Agilent Technologies) to deliver the stimulation trains with different parameters to any pair of the 16 electrodes.

#### Stimulation Paradigm

The sample stimulation trains for the LF and HF mode are illustrated in **Figure 1B**. In the LF mode, biphasic pulses were delivered at 30 Hz. In the HF mode, narrow pulses with a 60-µs pulse width were clustered in bursts with a 20-µs pulse interval, leading to a carrier frequency of 12.5 kHz. Different clusters were also delivered at 30 Hz. A pilot test demonstrated that a single 60 µs pulse evoked no EMG activity nor force responses.

The maximum voluntary contraction (MVC) of individual fingers was first obtained, when the subject maximally flexed the individual fingers against the load cell. Before the fatigue test, we first searched across all electrode pairs using the LF stimuli with a pulse duration of 600 µs to identify the electrode pair that can elicit the strongest muscle contraction with no pain sensation. Then, the current amplitude was adjusted until 30% of MVC was obtained for at least one finger. In the subsequent experiment, the current amplitude and the electrode pair were the same between the two stimulation modes. Muscle fatigue was evoked by delivering 20 (for the affected side) or 30 (for the contralateral side) 3-s stimulation trains with 1-s intervals (**Figure 1C**). Since muscles of the contralateral side are generally

the HF mode (C). The EMG signal of a representative channel recorded on the contralateral side in the HF mode (D).

less fatigable than that of the affected side, additional stimulation trains were delivered to ensure that fatigue can be induced. Two trials with one stimulation mode in each trial were performed in a random order. The resulting trial order was that the LF mode was first tested on the affected side, and the HF mode was first tested on the contralateral side. The initial contraction forces from the two modes were matched such that 30% MVC was elicited in the same finger. A 10-min rest was provided between trials, which was sufficient for the subjects to recover from muscle fatigue with fully recovered forces under the stimulation protocol (26).

#### Data Analysis Finger Force

To compare the overall force output, the raw forces from all fingers were summated at individual sample points, and the force-time integral was calculated based on the summated force. To further quantify the decline of forces between the two modes, the raw force data from each 3-s stimulation train were first averaged for individual fingers. Then, the average forces from all fingers were summated to represent the contraction force level of individual stimulation trains. Next, the force level were fitted with an exponential function, i.e.y = a + be<sup>τ</sup> <sup>t</sup> , where a represents the

force plateau level, b is the absolute force decay, and τ <0 reflects the rate of force decline. A direct comparison of the τ value cannot accurately reflect the force decline rate because the force plateau under two different modes can be different. Therefore, we estimated the number of stimulation trains (termed as the 50% peak stimulation) for the force to decline below 50% of the initial force level. A larger 50%-peak stimulation represents a slower force decay rate.

To estimate the force variability, the forces were first bandpass filtered (20–40 Hz) to eliminate the low-frequency offset and the high-frequency noise. Then, the standard deviation of the filtered force was calculated to represent the magnitude of the force fluctuation for individual stimulation trains. As in our previous study (25), the summated force fluctuation across all fingers was further normalized by the absolute force level of individual stimulation trains to eliminate the effect of force decline.

#### EMG Activity

Within each 3-s stimulation train, the EMG segments 5 ms prior to and 35 ms after the stimulation onset were averaged to obtain the average EMG for individual channels. Stimulation artifacts were then identified with a threshold-crossing method based on the EMG, and the EMG segment between two consecutive stimulation artifacts were extracted to calculate the EMG area, estimated as the integral of the absolute EMG over time (**Figure 1D**). The start and end point of the EMG integral calculation remained the same across the two stimulation modes. The estimated EMG areas were then averaged across all 128 channels to represent the overall level of the EMG activity. To evaluate the efficiency of muscle force generation, the force-to-EMG area ratio (F–EMG ratio) was calculated.

#### Results

**Figures 2A,B** illustrate the summated force of all four fingers between the two stimulation modes for the affected and contralateral sides. Each impulsive force was elicited by a stimulation train, resulting in 20 and 30 force pulses for the affected and contralateral sides, respectively. With matched initial force levels, the force declined notably as successive LF stimulation trains were delivered. In contrast, the force declined more slowly and stayed at a higher level under the HF stimulation. Compared with the LF mode, the force-time integral under the HF stimulation was 24.38 and 45.63% larger in the affected and contralateral side, respectively, indicating more force outputs under the HF stimulation (**Table 1**).

The changes of the contraction force are illustrated in **Figures 2C,D** for the affected and contralateral sides, respectively. For all conditions, the force decreased exponentially with the stimulation train. The average R-squared value of the exponential fit across all conditions was 0.9432 ± 0.0386. For both the affected and contralateral sides, the force decay was smaller and the force plateau was higher after successive HF stimulation trains, compared with that in the LF stimulation (**Table 1**). Under the LF stimuli, only 9 stimulation trains were needed for the force to decline below 50% of the initial force level for both the affected and contralateral sides. On the contrary,


the 50%-peak stimulation under the HF stimuli was 14 for the affected side and was even larger than 30 for the contralateral side (**Table 1**). **Figures 2G,H** illustrate the relative force decay between the affected and contralateral sides under different stimulation modes. Under the HF stimulation, the force had a higher plateau on the contralateral side compared with the affected side. Under the LF stimulation, on the contrary, the force declined continuously for both sides and the decay rate was similar between the two sides.

The normalized force fluctuation increased with successive stimulation trains (**Figures 2E,F**), which was consistent with our previous study (25). The increase of the normalized force fluctuation was larger under the LF stimulation compared with that under the HF stimulation for both the affected and contralateral sides.

The initial EMG activity distribution and the average EMG of the channel with the highest intensity from the first and the last 3-s stimulation train are shown in **Figure 3**. Similar contour lines between the two modes demonstrated that similar muscles or muscle portions were activated between the two modes for both sides. The amplitude of the EMG activity under the HF stimulation was smaller compared with the LF mode. The EMG activity intensity decreased substantially over time under the LF stimulation (**Figures 3D,I**) while the EMG activity under the HF stimulation did not show an obvious declining trend. The F-EMG ratio decreased in both stimulation modes, indicating a decreased efficiency of the muscle force generation after fatigue. However, the F-EMG ratio under the HF stimulation was consistently larger than that under the LF stimulation over time, indicating a consistently higher efficiency of force generation under the HF stimulation.

### DISCUSSIONS

In this case study, a HF stimulation protocol was tested on a hemiplegic stroke survivor to verify whether it can reduce muscle fatigue compared with the conventional LF stimulation. The results showed that the elicited force declined more slowly and the force plateau was higher under the HF stimulation for both sides, resulting in a more sustainable force output at higher levels. The force fluctuation increased more slowly, and was smaller with successive HF stimulation trains compared with the LF stimuli, indicating a less synchronized activation of the muscle fibers. The efficiency of the muscle activation measured as the F-EMG ratio (27) was higher after successive HF stimulation trains. These results indicate that the HF stimulation has the potential

to reduce muscle fatigue in stroke survivors, compared with the LF stimulation, by maintaining a higher efficiency of muscle activations.

The faster fatigue onset during electrical nerve stimulation can arise from the differences between physiologically activated MUs and electrically excited MUs. First, the recruitment order differs between the two activations (5, 6). FES preferentially excites the larger axons with lower resistance that innervate the more fatigable muscle fibers (28), although a random recruitment order has also been observed (29). Second, the electrical activation of MUs are highly synchronized and the twitch forces are time-locked to the stimulation. This is contrary to voluntary contraction where different MUs discharge asynchronously at a range of different firing rates. Therefore, the highly synchronized and physiologically reversed recruitment of MUs during FES replace the normal physiological recruitment of MUs and their discharge regulation, resulting in a fast decline of the elicited forces.

The mechanisms underlying reduced muscle fatigue in the HF stimulation may be multifactorial. First, the HF stimulation can temporally disperse the activation of different motor units (24). Each narrow pulse can only induce a subthreshold depolarization of the axon membrane and different axons might need different numbers of subthreshold depolarizations to reach the threshold, resulting in an enlarged delay between the onsets of action potentials of different axons. Since previous studies have demonstrated that asynchronized activations of different muscle fibers can delay the onset of muscle fatigue (15–17), the decreased synchronization level of motor units with the HF stimulation can thus be a possible mechanism for the higher fatigue-resistance. Second, the HF mode stimulation could elicit more H-reflex activities compared with the LF stimulation. At low levels of stimulation intensity, the Ia afferent fibers are preferentially activated due to their intrinsic properties and their larger diameters compared with the motor axons. Within each burst of the HF stimulation, the current was delivered intermittently. Since the Ia afferent fibers can be activated more easily than the motor axons, the HF stimulation burst could first activate Ia afferent fibers compared with the single wide pulse in the LF mode. As a result, a less number of action potentials may propagate antidromically along the motor axons, and more H-reflex activities (30) can be observed. H-reflex can lead to more physiological recruitment of the motor units. The initial EMG distribution showed similar patterns between the two stimulation modes, indicating similar muscles or muscle portions were activated. However, the H-reflex activity was still considered as a possible mechanism because the surface EMG grid can only capture the activities of superficial muscles and not the EMG activities of deep muscles through the H-reflex pathway. Lastly, a potentially changed axon excitability can also lead to reduced force generation (31). Compared with the LF mode, the HF mode possibly leads to fewer axon drop-out due to a smaller increment of the axon threshold associated with fatigue.

Even though the HF stimulation can prolong the force output for both the affected and contralateral sides, the amount of increased force outcome on the affected side was smaller compared with the contralateral side (**Table 1** and **Figure 2G**). The possible reason is that the muscle fibers of the paralyzed limbs after neurological disorders are innately more fatigable than muscle fibers on the contralateral side (7–10). An additional piece of support for this lies in the major differences in the decline of F-EMG ratio between the two sides. Although the F-EMG ratio after successive HF stimulations plateaued in the contralateral limb, which is consistent with our previous study (25), the affected limb showed continuous decline of F-EMG ratio. This difference can arise from a relatively higher force level of the affected side, given that the subject was severely impaired. The other possible explanation is the disrupted muscle force transmission to the tendon of the affected muscles following stroke (32). The transmission failure of the muscle

#### REFERENCES


force to the tendon counteracts the influence of the HF stimulation on muscle fiber activation. Therefore, even though the affected muscle still showed a decrease in the rate of force decline with HF stimulation, the improvement of fatigue resistance was smaller compared with that on the contralateral side.

#### CONCLUDING REMARKS

The current study showed that the HF stimulation at the proximal segment of the nerve bundles can reduce muscle fatigue in a stroke survivor compared with the conventional LF stimulation. The outcomes demonstrated the clinical significance of the technique to elicit sustainable muscle contractions at high force levels for individuals with neurological disorders.

#### AUTHOR CONTRIBUTIONS

YZ and XH conceptualized the study, YZ and HS performed the experiment, YZ analyzed the data, YZ and HS drafted the manuscript. YZ, HS, and XH revised the manuscript and approved the final version.

functional electrical stimulation for leg extension. Neurol Res. (2000) 22:703– 4. doi: 10.1080/01616412.2000.11740743


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zheng, Shin and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Variation of Finger Activation Patterns Post-stroke Through Non-invasive Nerve Stimulation

Henry Shin, Yang Zheng and Xiaogang Hu\*

*Joint Department of Biomedical Engineering, University of North Carolina at Chapel Hill and North Carolina State University, Raleigh, NC, United States*

Purpose: A transcutaneous proximal nerve stimulation technique utilizing an electrode grid along the nerve bundles has previously shown flexible activation of multiple fingers. This case study aimed to further demonstrate the ability of this novel stimulation technique to induce various finger grasp patterns in a stroke survivor.

Methods: An individual with chronic hemiplegia and severe hand impairment was recruited. Electrical stimulation was delivered to different pairs of an electrode grid along the ulnar and median nerves to selectively activate different finger flexor muscles, with an automated electrode switching method. The resultant individual isometric flexion forces and forearm flexor high-density electromyography (HDEMG) were acquired to evaluate the finger activation patterns. A medium and low level of overall activation were chosen to gauge the available finger patterns for both the contralateral and paretic hands. All the flexion forces were then clustered to categorize the different types of grasp patterns.

#### Edited by:

*Yingchun Zhang, University of Houston, United States*

#### Reviewed by:

*Hongbo Xie, Queensland University of Technology, Australia Yun Peng, NuVasive, United States*

> \*Correspondence: *Xiaogang Hu xiaogang@unc.edu*

#### Specialty section:

*This article was submitted to Stroke, a section of the journal Frontiers in Neurology*

Received: *17 September 2018* Accepted: *03 December 2018* Published: *13 December 2018*

#### Citation:

*Shin H, Zheng Y and Hu X (2018) Variation of Finger Activation Patterns Post-stroke Through Non-invasive Nerve Stimulation. Front. Neurol. 9:1101. doi: 10.3389/fneur.2018.01101* Results: Both the contralateral and paretic sides demonstrated various force clusters including single and multi-finger activation patterns. The contralateral hand showed finger activation patterns mainly centered on median nerve activation of the index, middle, and ring fingers. The paretic hand exhibited fewer total activation patterns, but still showed activation of all four fingers in some combination.

Conclusion: Our results show that electrical stimulation at multiple positions along the proximal nerve bundles can elicit a select variety of finger activation patterns even in a stroke survivor with minimal hand function. This system could be further implemented for better rehabilitative training to help induce functional grasp patterns or to help regain muscle mass.

Keywords: proximal nerve stimulation, neuromuscular electrical stimulation, stroke, electromyography, finger flexion

## INTRODUCTION

Following a stroke, a majority of individuals have paresis due to a loss of excitatory input and subsequent complications, such as disuse atrophy (1) and altered spinal organization (2–4). This loss of voluntary control of muscle activation often limits activities of daily living. Neuromuscular electrical stimulation (NMES) has been widely utilized both in the clinic and in research settings to help restore atrophied muscle and lost functions (5–7). Electrical stimulation has been particularly

**98**

successful with post-stroke survivors for functional recovery (8–10). Research in NMES also aims to restore functional activation of muscles, such as the restoration of hand grasps (11).

Traditionally, NMES uses large electrode pads, targeting the distal branches of the nerve, known as the motor point stimulation (12). Although stimulation of the motor point is straightforward methodologically, NMES is limited to localized muscle activation, which limits its functional efficacy and also leads to rapid muscle fatigue (13). Advances in NMES techniques to alleviate these issues involve various multielectrode techniques, which can stimulate multiple small regions of the muscle to help distribute the current and potentially activate more muscle fibers (14, 15). Crema et al. has also demonstrated flexible activation of multiple fingers using a multi-electrode array across the forearm and hand (16). Other approaches to NMES involve stimulation of the nerve bundle prior to branching and innervating a muscle, which has shown to allow for a larger area of muscle activation and potentially reduce long-term fatigue effects (17–19).

Recent developments have demonstrated the capabilities of an alternative non-invasive transcutaneous electrical nerve stimulation method targeting the ulnar and median nerves proximal to the elbow to flexibly activate individual and multiple fingers (20, 21). In addition, this technique shows the ability to delay the force decline (22, 23). A stimulation electrode grid placed along the two nerves allows us to activate different muscles or muscle portions to elicit varied desired movements, but manually switching between different electrode pairs is timeconsuming. To shorten this process, an automated electrode pair searching method has been developed and tested on intact control subjects (24). This new method can further categorize the total available sets of finger activation patterns across the entire electrode grid, providing valuable information on electrode selection and the force generation capacity of stroke muscles. However, the efficiency of this method has not been tested on stroke survivors. Therefore, this case study recruited a control subject and a stroke survivor with severe weakness of the right arm, and evaluated the available finger activation patterns of the subjects. Our results showed varied activation of multiple fingers from both subjects. Further development of this stimulation technique can provide valuable alternatives to current rehabilitation for the restoration of hand movements.

### METHODS

#### Case Report

A 54-year-old male who had a left hemisphere ischemic stroke 2 years ago was recruited. The participant had limited voluntary motion in the arm and hand with significant muscle atrophy but had no cognitive impairments. The average ratio of the subject's maximum finger forces between hands was 0.076, and the subject's Chedoke-McMaster Stroke Assessment hand score was 2, both indicating severe impairment. A 35-year-old male participant was also recruited as a neurologically-intact control subject for comparison. This study was carried out in accordance with the recommendations of the Institutional Review Board (IRB) of the University of North Carolina at Chapel Hill with written informed consent from the subject. The subject gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the local IRB. Additional written informed consent was obtained from the subject for the publication of this case report.

### Experimental Setup

To compare the proximal nerve stimulation method, the stroke subject's paretic and contralateral sides were tested on two separate occasions. The control subject was tested on the dominant arm. For each experiment, a 2 × 8 stimulation electrode grid was placed along the medial side of the upper arm below the biceps muscle where the ulnar and median nerves are more superficial (**Figure 1A**). All 16 stimulation electrodes were individually connected to a switch matrix (34904A, Agilent Technologies), which could be programmatically controlled. The switch matrix was then connected to a multi-channel stimulator (STG4008, Multichannel Systems), which could also be digitally controlled to deliver any range of current amplitudes between 0 and 16 mA, with a resolution of 20 µs. A brief cycle of 200 µs pulse width, 4 mA amplitude, and 30 Hz stimulation was delivered to every electrode to identify notably uncomfortable electrode combinations, which were then disabled. Following the stimulation setup, the skin of the anterior forearm was cleaned to reduce skin impedance for recording high density electromyography (HDEMG). An 8 × 16 HDEMG array (OT Bioelettronica)with a 10 mm interelectrode distance, was placed over the flexor compartment of the forearm (**Figure 1B**), and the 128 EMG channels were band-pass filtered at 10–900 Hz, with a gain of 500, and sampled at 5,120 Hz (EMG-USB2+, OT Bioelettronica). Lastly, each of the four fingers was individually secured to a uni-axial force transducer (SM-200N, Interface Inc.). Each finger was secured just above the metacarpophalangeal (MCP) joint (**Figure 1A**). The rest of the wrist and thumb were restricted to minimize force contamination. The force was recorded at 1,000 Hz. All the data recording and stimulation control were unified in a custom MATLAB GUI (Mathworks).

### Automated Stimulation Procedure

Once the setup was completed, the subjects were asked to perform maximum voluntary contractions (MVC) with each of the fingers individually and all 4 together. The stimulation procedure was composed of four steps:

#### Initial Medium-Level Grid Stimulation

An initial current level was chosen which can elicit some noticeable finger force without excessive contraction. For the paretic and contralateral sides of the stroke subject, the current levels were 6.5 and 4.5 mA, respectively, and 5 mA for the control subject. At this initial current level, all the different pair permutations were automatically switched and stimulated to test all the stimulation locations (120 maximum pairs). The bipolar stimulation consisted of trains of matched biphasic 200 µs pulse width and 30 Hz pulses. The stimulation was active for 0.5 s, and at rest for 1 s, while both the force and EMG were simultaneously recorded. Each pair was repeated

3 times before the stimulation was switched to the next pair (**Figure 1D**).

#### Max Force Selection and Activation Range Estimation

Once all available electrode pairs were stimulated, the subject was given a minute of rest while the GUI identified the stimulation electrode pair which resulted in the strongest average force across all the fingers within a single repetition. This electrode pair was then used to estimate the current-force relation across all of the other electrode pairs. Stimulation at 1 mA intervals was conducted to determine a rough estimate of the minimum and medium-high activation levels (relative to the current of the initial search) needed for each tested hand.

#### Current-Force Relation

Randomly chosen levels of current between these ranges (Paretic: 2–8 mA, Contralateral: 2–5 mA, Control: 2–7 mA) were stimulated using the same previous parameters, but with an increased 2-s rest between successive stimulations to reduce possible fatigue. Only the force was recorded with each stimulation, and the peak averaged force was calculated. The resultant Current-Force curve was then normalized to the corresponding averaged MVC.

#### Low-Level Grid Stimulation

The current value which activated around 5% MVC was selected to represent a level where low levels of finger motions were available. As the paretic MVC was already low, a value was chosen which was close to the lower take-off region of the current-force relation. The chosen values for the paretic and contralateral sides were 6 and 4.3 mA, respectively. The low-level selected for the control subject was 4 mA. The entire electrode grid underwent the automated stimulation procedure at these new current levels.

#### Data Analysis

The data were processed to simplify its comparison across electrode pairs. First the 30 ms of HDEMG data after each stimulation pulse were aligned and averaged to form a single compound muscle action potential (CMAP), which was again averaged across the 3 repetitions for each electrode pair. The Area-Under-the-Curve (AUC) of each CMAP was calculated as a measure of overall activity of a single EMG channel. These AUC Values were then placed in a 2D array which corresponded to its physical location on the forearm, and this overall heat map was used to compare the muscle activity. Additionally, the force data were smoothed using a 100-ms window with 1-ms steps and averaged across the 3 repetitions. Examples of the processed HDEMG and force data are shown in **Figures 1B,C**, respectively. Any electrode pairs which did not produce at least 0.1N of force in any single finger were excluded from further analysis.

Hierarchical clustering was utilized to categorize the different grasp patterns based on the force data. Since the force data was retained as a 1,500 × 4 array, a 2D correlation coefficient was calculated between the averaged force data of different electrode pairs. This value was considered as the distance between two electrode pairs, and then the complement of this correlation distance (1–Corr. Coef.), also known as the dissimilarity, was calculated between stimulation locations. Using an inconsistency cutoff of 1.1, the initial hierarchical clusters were then used as a starting point to further refine the force patterns caused by each stimulation location. The silhouette coefficient was used as a measure of cluster validity, and therefore each member of each group was shuffled until the average silhouette coefficient of each group was maximized. For each resultant final cluster group, the average ratio of force between each finger was used to threshold whether or not the finger was active. Each of these clusters were therefore labeled by its finger activation.

Lastly, the EMG Activity AUC Maps of each force cluster were correlated with each other to quantify the similarity between the EMG activation of each force cluster. A high correlation indicates that electrode pairs within a force cluster also produces similar EMG activity, and inversely, a low correlation indicates that the electrode pairs within a cluster may produce similar force, but through different muscle portions.

#### Results

The maximum voluntary forces obtained from the subject showed a large disparity between sides. On the contralateral side, the subject's individual finger forces were 19.9, 21.1, 30.1, and 21.9 N for the Index, Middle, Ring, and Pinky fingers, respectively (Average: 23.2 N). For the paretic side, the finger forces were 0.7, 3.6, 2.4, and 0.3 N (Average: 1.8 N). The finger maximums for the control subject were 26.7, 26.3, 20.6, and 26.9 N (Average: 25.1 N). These values were in a similar range as the stroke subject's contralateral side. Although initially obtained to normalize the elicited forces, due to the very low forces of the paretic side, the absolute forces were reported.

The clusters obtained from the contralateral and paretic hands, as well as the control subject, are shown in **Figure 2**. The labels on the top of each cluster indicates which of the four recorded fingers were active. The Contralateral hand and the control subject showed a variety of single and multifinger activation patterns which were mostly an activation of the Index, Middle, and Ring fingers, but also a few Pinky. Similarly, the Paretic hand also resulted in several clusters of activation patterns, although fewer than the contralateral side. The paretic hand resulted in relatively more clusters with only a single finger being active, but still had a few two and three finger activations. Overall, the contralateral hand and the control subject clusters show that the electrode grid had strong preference to the median nerve (Index-Middle-Ring), whereas the paretic-side electrode grid may have had a more evenly divided placement between the two desired nerve bundles.

The results of the AUC Correlation analysis are shown in **Figure 3**. **Figure 3A** shows two examples of EMG activity with high cluster correlation and low cluster correlation. **Figure 3B** illustrates the individual cluster AUC Correlation for two sides of the stroke subject and the control subject. These results suggest that for each force cluster in **Figure 2**, there is a high variability of EMG correlation. Some force clusters also have high EMG activity correlation, whereas other force clusters may have more varied EMG activity, and therefore lower correlation.

## DISCUSSION

In the current case report, an individual with chronic stroke associated muscle weakness was tested with our proximal nerve stimulation system alongside a neurologically-intact control subject to evaluate the capabilities to elicit specific finger activations in highly paretic muscle. Overall, our results demonstrate that our stimulation system is able to activate various different fingers on both the contralateral and paretic sides of this subject.

As a survey of the available finger activations of the stimulation method, our results highlight a few important aspects of the activated finger patterns. Similar to previous results (24), a majority of the finger activations were those of the Index, Middle, and Ring fingers. These correspond to the median nerve, and therefore it can be concluded that the placement of the electrode stimulation array was preferential to the median nerve, especially for the contralateral hand and in the control subject. The sets of force clusters from these two conditions also demonstrate similar ranges of single and multi-finger activation patterns. As for the paretic hand, the activation of the Pinky finger suggests that more of the ulnar nerve was also accessible. As the paretic biceps muscle was also visibly atrophied when the stimulation electrodes were applied, the underlying nerve bundles may have been more easily accessible. Additionally, the results suggest that different electrode pairs are able to activate different portions of the corresponding nerves. Although the different clusters are a post-hoc attempt at organizing the finger forces generated by each electrode pair, in reality, each elicitable finger pattern lies along a continuum of finger activation levels. Different electrode pairs impose a unique electrical field onto the nerve, and thus activates a unique subset of motor and/or sensory axons. As shown in the contralateral force clusters, many of these subset axons project to muscles spanning multiple fingers, but a small number of the electrodes can partially activate a single finger. The different clusters help to identify which sets of electrode pairs can elicit desirable finger grasp patterns. Additionally, although anatomical landmarks were used during the stimulation grid placement, small variations in the location of the electrodes could inevitably lead to different sets of activation even within the same subject. Therefore, although there may be

FIGURE 2 | All Finger Force Clusters. The hierarchical clustering results of a control subject and the two sides of a stroke subject is visualized here. As shown in Figure 1, each finger is individually plotted from the force profiles of each electrode pair in a cluster. The average of all of the force profiles in a cluster are also shown as the darker solid line within each axes. The labels above each cluster represent which fingers were deemed "active" based on its relative ratio to the other fingers. The letters for each finger are, I-Index, M-Middle, R-Ring, and P-Pinky.

inter-session differences in the exact number and range of force clusters, the general similarity between the control subject and the contralateral side suggest that the stimulation grid is able to activate similar patterns of finger movement.

Although the obtained number of clusters are not necessarily indicative of any physiological correlation with muscle function, it is important to note that the paretic side does have notably fewer number and variety of force clusters (12 on contralateral/control vs. 7 on paretic side). This may be due to several factors that occur in paretic muscle after stroke. The first is due to the significant atrophy and weakness of the right arm and hand seen in the participant. As the overall stimulated force level was still very low, the limited muscle mass may not have been able to generate observable levels of force in some activation patterns. Alternatively, various losses in the motor unit (MU) numbers as well as reinnervation of existing MUs are also common occurrences after stroke (25, 26). This may alter the various subsets of activation available through nerve stimulation. Further studies are needed to confirm these possibilities, as the lower number of clusters may also simply be due to the chance involved with electrode placement. Clearly, additional testing involving a large stroke cohort is necessary.

Along with the overview of the different force clusters, as an estimation of the total available activation patterns, the EMG AUC Map correlation within each cluster also provides further insight into the actual EMG activity from each stimulation. **Figure 3** shows that within each force cluster, there may be a wide distribution of similar and dissimilar EMG activity. A force cluster that has high EMG AUC correlation, may imply a set of electrode pairs which consistently activate the same portion of muscles. Contrastingly, a force cluster that has a low EMG AUC correlation may suggest that the set of electrode pairs in the cluster is able to activate different portions of muscles even with the same resultant physical activity. This feature may have far reaching effects regarding prolonging use with stimulation redundancy. Various groups have shown that, distributed muscle activation using multiple pads is able to lead to reduced fatigability in NMES (27–29). Being able to alternate between different electrode pairs which recruit different muscles with similar functional outputs could reduce muscle fatigue in the stimulation. Further investigation on the relation between the force variability and EMG variability is necessary to better understand the different electrode pair choices and their impacts on force production.

### CONCLUDING REMARKS

The current study demonstrates the variety of finger activation patterns that are accessible through our proximal nerve stimulation method. Both the contralateral and paretic sides of a stroke subject were able to be successfully stimulated to produce a number of multi-finger movements along with a few isolated single fingers. The contralateral hand, in particular, was able to elicit a similar variety of finger activation patterns as seen in the control subject. Further development of the technique can also be combined with radial nerve stimulation to also include hand opening, which is just as important for stroke survivors. The automated electrode searching with the force clustering can help rapidly identify the feasible sets of electrode pairs, which can allow us to develop an auto-calibration method between

#### REFERENCES


sessions/days, which can then be applied to any future uses for rehabilitation.

#### AUTHOR CONTRIBUTIONS

HS and XH conceptualized the study design. HS and YZ performed the experiment. HS analyzed the data and drafted the manuscript. HS, YZ, and XH revised the manuscript and approved the final version of the manuscript.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Shin, Zheng and Hu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.