# 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<EMAIL_ADDRESS>
# 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<EMAIL_ADDRESS>Rong Song<EMAIL_ADDRESS>
*† 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.<PHONE_NUMBER>34,<PHONE_NUMBER>08).

### 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<EMAIL_ADDRESS>Wufan Chen<EMAIL_ADDRESS>Xiaoyun Wang<EMAIL_ADDRESS>
#### 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.  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<EMAIL_ADDRESS>Wensheng Hou<EMAIL_ADDRESS>
#### 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<EMAIL_ADDRESS>Wen S. Hou<EMAIL_ADDRESS>
#### 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),