1 Introduction

The acquisition of crawling skills is regarded as one of the most crucial developmental milestones in human motor skill development (Chen et al., 2017; Xiong et al., 2021). Various physical abilities can be promoted by crawling movement in infancy, including eye-hand coordination, balance, and spatial concepts (McEwan et al., 1991). Therefore, the evaluation of the crawling function has great application value in the fields of disease diagnosis and rehabilitation treatment (Gao et al., 2018; Xiong et al., 2018a,b, 2021). In particular, as a typical quadruped movement (Cole et al., 2019), there are various kinds of inter-limb coordination modes (ILCMs) during crawling. Figuring out the ILCMs during crawling may help clinicians better evaluate patients’ motor dysfunction and then develop more precise rehabilitation treatment plans.

Early studies on ILCMs during crawling were mainly based on the subjective observation and judgment of observers (Hildebrand, 1967; Freedland and Bertenthal, 1994). In this century, various kinetics-based sensing technologies, i.e., motion capture technology (Patrick et al., 2009, 2012; MacLellan et al., 2012, 2013, 2017; Righetti et al., 2015; Gao et al., 2018; Xiong et al., 2018a,b), inertial measurement unit (IMU) (Vitali et al., 2019), and 3-axis accelerometers (Ma et al., 2017), have been successfully introduced into motion analysis. Patrick et al. (2009) defined ipsilateral phase lag (IPL) value to quantify ILCMs during human crawling. In their definition, when a represented the whole crawling cycle and b represented the phase difference between the moment when the left palm and the left leg contacted the ground, IPL can be calculated as (b/a) * 100%. IPL values closed to 50% indicated trot gait, where diagonal limbs moved in coordination; IPL values closed to 0 or 100% indicated pace gait, where ipsilateral limbs moved in coordination; IPL values closed to 25% or 75% indicated no-limb-pairing gait, where all limbs moved at regular intervals (Patrick et al., 2009). Using a motion capture system, Patrick et al. carried out research on the ILCMs during human crawling and stated that infants would like to adopt trot gait (Patrick et al., 2009). Using 3-axis accelerometers and pressure sensors, Ma et al. investigated hands-knees crawling in adult humans and found that, at low speeds, most adults crawled using trot or no-limb-pairing gait, while they tended to use trot or pace gait as their crawling speed increased (Ma et al., 2017).

In summary, current research on ILCMs was in the preliminary stage, and the research results had certain limitations. First, crawling is a full-body movement related to the coordinated contraction of a series of limbs and trunk muscles. It is universally acknowledged that motion capture systems have high light requirements in the experimental environment, which is usually a fixed experimental site and is easily affected by occlusion, resulting in data loss. Inertial sensors such as IMU and ACC will produce cumulative errors, affecting the accuracy of analysis. Specifically, neither subjective observation nor kinetics-based sensing technologies can reflect the motor control characteristics of the central nervous system (CNS) from the perspective of muscle contraction. Second, most of the existing studies divided the ILCMs into trot gait, pace gait, and no-limb pairing gait according to the IPL value (Patrick et al., 2009, 2012; MacLellan et al., 2012, 2013; Righetti et al., 2015; Chen et al., 2017; Ma et al., 2017; Vitali et al., 2019). However, the IPL was difficult to distinguish modes with subtle differences. For instance, during hands-knees crawling, the trot gait (diagonal limbs moving together) had the same IPL value (50%) as sequential gait (left palm → right knee → left knee → right palm). Third, although relevant studies have found that human beings have different choices for the ILCMs at different crawling speeds (Chen et al., 2017; Ma et al., 2017), there was a lack of exploration of the possible reasons for making the choice.

In view of the shortcomings of existing research, this study attempted to carry out research on the pattern recognition of ILCMs during crawling by means of electromyography (EMG) signals, which carry abundant muscle activation information and neuromuscular control information of the CNS. Compared to motion capture systems and inertial sensors, EMG signals do not require a specific experimental site and are not affected by cumulative errors. Therefore, it has been widely used for pattern recognition of human motion intention, such as ankle joint movements (Al-Quraishi et al., 2017), lower limb jump locomotion phases (Lu et al., 2021), hand gestures (Côté-Allard et al., 2019), and muscle forces (Mokri et al., 2022). The following benefits can be achieved by introducing EMG-based motion intention recognition technology into the inter-limb coordination pattern recognition. First, given the successful application of EMG in fine finger motion recognition (Côté-Allard et al., 2019), we believe that it can achieve better accuracy in coarse crawling motion recognition; Second, unlike roughly dividing the ILCMs into trot, pace, and no-limb-pairing gait, using EMG signals recorded from related muscles was expected to achieve more accurate classification; Third, the pattern recognition scheme for ILCMs based on EMG signals can provide a novel crawling motion analysis technology, which was helpful for clinicians to evaluate patients’ motor dysfunction status from the perspective of muscle function and also helped researchers understand the neuromuscular control mechanism during human crawling.

The fundamental purposes of this study are as follows: (1) first introduce EMG-based motion intention recognition technology into crawling motion classification; (2) explore the reasons for the choice of different ILCMs at various crawling speeds based on the classification results and kinetic parameters of self-selected ILCM under the pattern recognition scheme. The innovations and primary contributions can be summarized as follows: (1) Unlike relevant research studies, which roughly divided the crawling modes into three gaits, this study targeted eight ILCMs defined by the sequence in which limbs contact the ground; (2) The feasibility of providing clinicians with an EMG-based ILCM pattern recognition scheme has been verified by pattern recognition experiments in participant-specific way, multi-participant way, and participant-independent way respectively, with three classifiers including BiLSTM, SVM, and KNN; (3) By classifying the participants’ self-selected ILCM at different speeds using the best classification model, the possible reasons why the choice of the ILCM changed with crawling speed were discussed.

2 Materials and methodology

The research flowchart of this study is illustrated in Figure 1, and each step is described in detail below.

Figure 1. Flowchart of the research route.

2.1 The experimental design and crawling data acquisition

Eight ILCMs were defined, as shown in Table 1. Each crawling cycle was initiated by the left palm touchdown and defined as the time interval between two consecutive left palm touchdowns (Patrick et al., 2009; MacLellan et al., 2012; Chen et al., 2017).

Table 1. The eight designed inter-limb coordination modes.

The crawling data acquisition experiment involved the participation of 10 healthy adults, comprising three female subjects and seven male subjects, with an average age of 23.90 ± 0.88 years. None of the participants had a previous diagnosis of neuromuscular disorders. This study was approved by the Ethics Review Committee of Anhui Medical University (No. PJ 2014-08-04). All participants were informed about the study details and provided written informed consent.

A laboratory-made multi-channel system with 30 EMG sensors and 1 pressure sensor was used to collect crawling data. As shown in Figure 2A, EMG electrodes contained bipolar separating silver wires 23 mm in length and 20 mm in width, with a 10-mm interval between them (Li et al., 2023). As shown in Figure 2B, EMG signals were recorded from 15 muscles on each side of the body and 30 muscles from the whole body in total, which were highly related to crawling, including anterior deltoid (AD), adductor longus (AL), biceps brachii (BB), biceps femoris (BF), brachioradialis (BR), extensor carpi radialis (ECR), flexor carpi radialis (FCR), latissimus dorsi (LD), rectus femoris (RF), sartorius (SA), semitendinosus (SE), triceps brachii (TB), trapezius (TR), vastus lateralis (VL), and vastus medialis (VM). The placement of EMG sensors was based on the guidelines of the SENIAM protocol (Hermens et al., 2000). Before placing the sensors, the target muscles were shaved and cleaned with alcohol swabs. To detect crawling cycles, a pressure sensor was affixed with kinesiology tape to the flexor pollicis brevis muscle of the left palm. The company WAAAX manufactured the RP-C18.3-ST thin-film piezoresistive pressure sensor, which has a diameter of 18.3 mm. The sampling rate of EMG electrodes and pressure sensor was set to 1,000 Hz.

Figure 2. (A) Illustration of the homemade data acquisition system; (B) placement of EMG sensors and pressure sensors.

Throughout the data acquisition experiment, participants crawled on a sponge pad measuring 11.2 m in length and 0.8 m in width. First, each participant crawled in self-selected ILCM at their self-understanding slow, medium, and fast speeds. From the pressure signal on the left palm of the participants, the low, medium, and fast speeds in the self-selected ILCM were calculated to be approximately 2 s/cycle, 1.5 s/cycle, and 1 s/cycle, respectively. Then, each participant was asked to crawl in the eight ILCMs (M1 to M8) at three different speeds, namely slow speed (approximately 3.5 s/cycle), medium speed (approximately 2.33 s/cycle), and fast speed (approximately 1.75 s/cycle). To help participants complete crawling movements at specific speeds and ILCMs, a series of audio files prompting the landing order of each limb was generated for the eight ILCMs at three crawling speeds. Before data collection, the participants learned how to crawl under the alert of audio files until they became proficient. For any combination of speed (fast, medium, and slow) and ILCM (M1 to M8, self-selected), as a crawling trial, pressure signals and EMG signals for at least 15 consecutive crawling cycles were collected. To minimize muscle fatigue, participants were required to rest for about 10 min between each group of 6 crawling trials. Table 2 displays the total number of crawling cycles for each participant across all eight ILCMs (M1–M8) and their self-selected crawling mode.

Table 2. Number of crawling cycles.

It is worth noting that this study proposed the following countermeasures for the artifacts that may be introduced by data collection systems in different scenarios: (1) To avoid power frequency noise, lithium battery was used to power the data acquisition system; (2) The experimental operator checked the interface connection between each electrode and the data acquisition system to ensure that all sensor connections were normal before the experiment begins; (3) The sensor connecting lines located on the same limb were fixed together with adhesive tape to reduce the disturbance; (4) During crawling data collection, the experimental operator observed the real-time signal on the laptop and independently stored the data of each crawling mode. Once the artifact was discovered, the experiment was suspended, and the current experimental data were discarded. Only after the cause of the artifact was found and the fault was eliminated can the experiment be restarted.

2.2 Data preprocessing

Figure 3A illustrates the pressure signal obtained during a consecutive crawling movement. When the left palm touched the land, the left palm entered into the stance phase and the pressure signal reached the maximum value; when the left palm left the land, the left palm entered into the swing phase and the pressure signal returned to zero. Figure 3B shows that the first derivative of the pressure signal during left-palm stance and swing phase was zero. When transitioning from the swing phase to the stance phase, the first derivative became greater than zero; when transitioning from the stance phase to the swing phase, it became less than zero. Therefore, this study utilized the aforementioned characteristics of the pressure signal’s first derivative to segment crawling cycles. In Figure 3, the red asterisks indicated crawling cycle starting points and the purple asterisks indicated swing phase starting points. The crawling cycle was segmented by two adjacent red asterisks.

Figure 3. (A) Illustration of the pressure signal collected from left palm; (B) the first derivative of the pressure signal. Red asterisks indicate cycle starting points, as the beginning point of the stance phase, and purple asterisks indicate the beginning point of the swing phase; (C) illustration of raw EMG signal from brachioradialis; (D) the envelope of EMG signal before cycle normalization and amplitude normalization; (E) the envelope of EMG signal of six muscles after cycle normalization and amplitude normalization.

Various EMG features can be extracted for EMG-based motion intention pattern recognition, including (1) time-domain features, such as mean absolute value, slope sign changes, and waveform length (Hudgins et al., 1993); (2) frequency-domain features based on fast Fourier transform, i.e., intermediate frequency (MDF), average frequency (MNF), autoregressive coefficient (AR), etc. (Abdelouahad et al., 2018); and (3) time-frequency-domain features based on Wigner-Ville transform (WVT), wavelet transform, etc. (Li et al., 2020). Given that the variation in muscle activation intensity can indicate the differences between ILCMs, this study opted to utilize an EMG signal envelope as a feature for pattern recognition. In particular, the EMG signals of each crawling cycle were high-pass filtered, demeaned, rectified, and low-pass filtered to extract the envelope (Li et al., 2023). Then, the envelope amplitude was normalized to unit variance (Teruya et al., 2021), and each crawling cycle length was normalized to 1,000 points. As an example, Figures 3C, D present the original EMG signal fragment from the brachioradialis and its envelope, respectively, before cycle normalization and amplitude normalization. Figure 3E demonstrates the EMG signal envelopes of six muscles after undergoing cycle normalization and amplitude normalization.

2.3 EMG-based crawling pattern recognition scheme

To verify the feasibility of providing clinicians with an EMG-based scheme for accurate recognition of ILCM, the pattern recognition experiments on the eight defined ILCMs were carried out in the participant-specific way, multi-participant way, and participant-independent way, respectively, at four crawling speeds, including low, medium, fast, and mixed speed (EMG data mixed from low, medium, and fast speeds), with three classifiers, namely BiLSTM network, SVM, and KNN.

2.3.1 Three classifiers

BiLSTM, SVM, and KNN were adopted to complete the pattern recognition task for crawling motion with different ILCMs based on the following considerations: (1) SVM model has kernel trick characteristic, KNN has non-parametric nature (Samuel et al., 2019), and these two classifiers have been widely adopted within the realm of EMG-based pattern recognition; (2) both KNN and SVM classifiers were characterized by their ease of implementation and efficient training (Samuel et al., 2019); and (3) LSTM model, which was good at memorizing the timing correlation, has been applied successfully in EMG-based gesture recognition (Chen et al., 2020). As a variation of LSTM, the BiLSTM network had better performance than the regular LSTM network. KNN and SVM were implemented by using Python’s sklearn toolkit. This section mainly introduced the implementation of the BiLSTM network.

Ff=δWfh·ht−1+Wfx·xt+bf    (1) Fu=δWuh·ht−1+Wux·xt+bu    (2) Fo=δWoh·ht−1+Wox·xt+bo    (3) C¯t=tanhWch·ht−1+Wcx·xt+bc    (4) Ct=Ff·Ct−1+Fu·C¯t    (5)

To better understand how the BiLSTM network works, we first figured out its most important part, namely the LSTM unit (Choi et al., 2019; Marentakis et al., 2021). As shown in Figure 4A, an LSTM unit consisted of several main parts, including the input information xt at step time t, memory cell state Ct, temporary memory cell state C¯t, hidden state ht, and three gates (forget gate Ff, update gate Fu, and output gate Fo). The previous memory cell Ct−1, previous hidden state ht−1, and current input xt decided the output hidden state ht and memory cell Ct together. The temporary memory cell state C¯t was decided by a tanh layer based on the previous hidden state ht−1 and current input xt. The functions of the three gates can be summarized as follows: the forget gate Ff decided the information to be thrown away by a sigmoid layer, the update gate Fu decided that the information should be stored in the next LSTM unit and the output gate Fo decided what should be ultimately output. The specific calculation process of the LSTM unit can be seen from formula (1)–(6), where σ (·) was the sigmoid function, W represented weight matrices, and b denoted biases.

The BiLSTM network structure adopted in this study is presented in Figure 4B. The working principle of the BiLSTM network was that the input data sequenceX=X1,X2,…Xt was not only fed to a LSTM network named the forward layer but also was simultaneously fed to another LSTM network named the backward layer by reversing the order of the input sequence of data. In other words, the first sequence X1 was the input of the first unit of the forward layer and also the input of the last unit of the backward layer. The output of the forward layer and the output of the backward layer jointly decided the final output Y1. After the BiLSTM layer, four dense layers were adopted. Python’s Keras toolbox was adopted to implement this network. The adaptive moment estimation (Adam) was used to avoid overfitting. The hyperparameters were listed as follows: learning rate = 0.001, batch size = one-tenth of the training samples, beta_1 = 0.9, beta_2 = 0.999, epsilon = 1e-08, LSTM hidden units = 128, and training epochs = 60.

2.3.2 Three pattern recognition ways

In the participant-specific way, the test data and training data were from the same participant. For each participant, a three-fold cross-validation was performed, where two-thirds of the data were used for training and the remaining one-third was used for testing.

In the multi-participant way, the data from all 10 participants were combined, and three-fold cross-validation was adopted, where two-thirds of the data were used for training and the remaining one-third was used for testing.

In the participant-independent way, the leave-one-out strategy was adopted. Classifiers were trained using crawling data obtained from nine participants, and the well-trained classifiers were applied to test the remaining crawling data.

2.3.3 Performance evaluation and statistical analysis

The pattern recognition accuracy was determined by calculating the ratio of correctly identified samples to the total number of samples. To investigate the effects of the independent variables (pattern recognition way, classifier, and crawling speed) on the recognition accuracy, a one-way ANOVA and univariate ANOVA were conducted using IBM SPSS Statistics 26. The significance level was set at 0.05.

2.4 Analysis scheme for the possible reasons for the choice of self-selected crawling mode

As shown in Figure 1, the analysis of the possible reasons why participants’ choice of ILCM varied with crawling speed was based on the classification results of self-selected ILCM and the statistical results of the stance duration, swing duration, and duty factor of the stance phase.

The classification model, performing the best in the pattern recognition experiment on the eight defined ILCMs, was used to classify the participants’ self-selected crawling mode at different speeds. The stance phase durations and swing phase durations for the left palm were calculated according to the pressure signal, as shown in Figure 3A. The duty factor of the stance phase was calculated according to the following formula (7):

Dutyfactor=StancephaseStancephase+Swingphase∗100%    (7) 3 Results 3.1 The pattern recognition results for the eight defined ILCMs

To demonstrate the feasibility of the EMG envelope feature in distinguishing different ILCMs, the t-SNE dimensionality reduction algorithm (Van der Maaten and Hinton, 2008) was carried out on all the 2,736 crawling cycles of M1 ~ M8, and Figure 5 shows the visualization results. It can be observed that the eight defined ILCMs can be clearly distinguished from each other.

Figure 5. The t-SNE dimensionality reduction of EMG envelope samples of the eight inter-limb coordination modes.

Figure 6 shows the classification accuracies averaged across all 10 participants in the participant-specific way. BiLSTM obtained average accuracies of 99.23, 99.09, 99.66, and 98.71% at low, medium, fast, and mixed speeds, respectively. KNN achieved average accuracies of 99.23, 99.69, 99.29, and 99.65% at low, medium, fast, and mixed speeds, respectively. SVM obtained average accuracies of 99.55, 99.69, 99.36, and 98.75% at low, medium, fast, and mixed speeds, respectively.

Figure 6. Classification accuracies were averaged across all 10 participants in the participant-specific way.

Figure 7 shows the classification accuracies in the multi-participant way, and BiLSTM obtained the average accuracies of 98.78, 98.91, 98.56, and 99.08% at low, medium, fast, and mixed speeds, respectively. KNN obtained the average accuracies of 99.44, 99.67, 99.78, and 99.78% at low, medium, fast, and mixed speeds, respectively. SVM obtained average accuracies of 99.44, 99.56, 99.67, and 99.60% at low, medium, fast, and mixed speeds, respectively.

Figure 7. Classification accuracies in the multi-participant way.

Figure 8 shows the classification accuracies in the participant-independent way, and BiLSTM obtained the average accuracies of 95.76, 95.51, 88.56, and 95.42% at low, medium, fast, and mixed speeds, respectively. KNN obtained the average accuracies of 96.21, 97.54, 94.43, and 96.89% at low, medium, fast, and mixed speeds, respectively. SVM obtained average accuracies of 98.31, 98.89, 95.60, and 98.31% at low, medium, fast, and mixed speeds, respectively.

Figure 8. Classification accuracies were averaged across all 10 participants in the participant-independent way.

Considering that the SVM classifier achieved the best performance in the participant-independent way, taking the 302 crawling samples from participant P1 as test data, the confusion matrix for the SVM classifier at mixed speed was given in Figure 9. Only one sample of M1 and one sample of M2 were misclassified as M5, and one sample of M7 was misclassified as M6. The recognition accuracy of some ILCMs, such as M3, M4, M5, M6, and M8, was 100%.

Figure 9. Confusion matrix for SVM classifier at mixed speed in the participant-independent way.

Table 3 presents the results of the statistical analysis examining the effects of pattern recognition way, classifier, and crawling speed on recognition accuracy. Based on the findings, the following conclusions can be drawn: (1) Pattern recognition way had a considerable effect on the recognition accuracy. More precisely, the recognition accuracy in the participant-independent way was notably inferior to that in the participant-specific way (p = 0.000**) and multi-participant way (p = 0.000**). Nonetheless, no significant distinction was found (p = 0.938) between the participant-specific way and the multi-participant way; (2) Crawling speed had a significant impact on recognition accuracy. However, only fast speed had a significant difference with other speeds (low speed, p = 0.015*; medium speed, p = 0.003*; mixed speed, p = 0.008*); (3) Classifier employed had a marked effect on recognition accuracy. BiLSTM obtained significantly lower recognition accuracy than KNN (p = 0.021*) and SVM (p = 0.001*). Nonetheless, no significant difference (p = 0.236) was found between the recognition accuracy achieved using KNN and SVM.

Table 3. The results of multi-variable statistical analysis for crawling pattern recognition.

3.2 Statistical results of stance duration, swing duration, and duty factor of stance phase

Figure 10 illustrates the statistical results of swing durations, stance durations, and the duty factors of the stance phase of the left palm. The results show that, as crawling speed increased, there was a substantial decrease in the duration of the stance phase (Low: 1.437 ± 0.255 s; Medium: 1.009 ± 0.220 s; High: 0.717 ± 0.135 s), while the duration of swing phase remained unchanged or slightly shortened (Low: 0.601 ± 0.149 s; Medium: 0.531 ± 0.078 s; High: 0.496 ± 0.089 s). On the whole, the duty factor of the stance phase decreased with the crawling speed (Low: 70.40 ± 6.54%; Medium: 64.97 ± 5.32%; High: 58.98 ± 4.81%).

Figure 10. Illustration for swing durations, stance durations, and the duty factors of the stance phase of the left palm in self-selected crawling mode.

3.3 Classification results of self-selected ILCM

To figure out the possible reasons why the participants’ self-selected ILCM changes with crawling speed, the self-selected ILCM was classified using the KNN classifier, which was trained in the multi-participant way at mixed speed, and the results were shown in Figure 11.

Figure 11. Classification result of self-selected crawling mode.

Regardless of the crawling speed, M2 accounted for the largest proportion. At low, medium, and fast speeds, 87, 92, and 68 cycles were classified as M2, respectively. Meanwhile, the proportion of M3 changed with the crawling speed. Concretely, the crawling cycle number of M3 was 9, 0, and 23 at low, medium, and fast speeds, respectively. That is to say, when participants crawled at low speed or fast speed, although the trot gait M2 was the most adopted one, some participants preferred to choose M3 instead of M2, especially at fast speed. In addition, some ILCMs, such as M4, M6, M7, and M8, were adopted by nobody at any speed.

4 Discussion 4.1 The clinical application value of the EMG-based crawling pattern recognition scheme

Humans or quadruped animals have multiple ILCMs. Hildebrand divided the ILCMs of quadruped into two categories: symmetrical gaits and asymmetrical gaits (Hildebrand, 1976). Furthermore, symmetrical gaits were divided into four ILCMs in general: pace, lateral sequence, trot, and diagonal sequence (Hildebrand, 1976). Owaki et al. concluded that there were nine ILCMs in quadruped animals, i.e., lateral-sequence walk, diagonal-sequence walk, trot, pace, pronk, canter, bound, transverse gallop, and rotary gallop (Owaki and Ishiguro, 2017). Bellardita et al. found that wild-type mice had four ILCMs, namely, walk, trot, bound, and gallop (Bellardita and Kiehn, 2015). As for human crawling motion, IPL value has usually been used to classify ILCM into pace, trot, and no-limb-pairing modes in most studies adopting observation method or kinetics-based sensing technologies (Patrick et al., 2009, 2012; MacLellan et al., 2012, 2013; Righetti et al., 2015; Chen et al., 2017; Ma et al., 2017; Vitali et al., 2019). However, when ILCMs were simply considered as these three modes, many details of human crawling motion were ignored.

Unlike relevant research studies, eight ILCMs were defined in this study, and the EMG-based pattern recognition scheme was first introduced into the classification of crawling ILCMs. The experimental results demonstrated that the eight defined ILCMs could be distinguished with relatively high accuracies at different crawling speeds using three classifiers. Additionally, the classification performance of the EMG-based crawling pattern recognition scheme is very robust and less affected by crawling speed. Even in the participant-independent way, the SVM classifier achieves above 98% recognition accuracy at low, medium, and mixed speeds. In contrast to the rough classification based on IPL value, the EMG-based pattern recognition scheme could classify more d