A Novel Muscle Innervation Zone Estimation Method Using Monopolar High Density Surface Electromyography

This study presents a novel method to estimate a muscle’s innervation zone (IZ) location from monopolar high density surface electromyography (EMG) signals. Based on the fact that 2nd principal component coefficients derived from principal component analysis (PCA) are linearly related with the time delay of different channels, the channels located near the IZ should have the shortest time delays. Accordingly, we applied a novel method to estimate a muscle’s IZ based on PCA. The performance of the developed method was evaluated by both simulation and experimental approaches. The method based on 2nd principal component of monopolar high density surface EMG achieved a comparable performance to cross-correlation analysis of bipolar signals when noise was simulated to be independently distributed across all channels. However, in simulated conditions of specific channel contamination, the PCA based method achieved superior performance than the cross-correlation method. Experimental high density surface EMG was recorded from the biceps brachii of 9 healthy subjects during maximum voluntary contractions. The PCA and cross-correlation based methods achieved high agreement, with a difference in IZ location of 0.47 ± 0.4 IED (inter-electrode distance = 8 mm). The results indicate that analysis of 2nd principal component coefficients provides a useful approach for IZ estimation using monopolar high density surface EMG.

where terminal branches of a motor neuron contact muscle fibers. Topographically, the IZ often lies within in narrow band near the mid-point of the muscle fibers. When a nerve impulse activates the muscle fibers, motor unit action potentials (MUAPs) propagate in opposite directions from the IZ toward the muscle tendons. Knowledge of IZ localization is of great importance for both basic science and clinical investigations of muscle morphology in health and disease [1], [2], [3], [4]. For example, IZ localization can be used to guide botulinum toxin injection for spasticity treatment in clinical practice [5], [6], [7], [8]. The characteristics and positioning of recording electrodes in relation to the IZ can affect the spectral and temporal features of surface electromyography (EMG) signals [9], [10], [11], [12], [13], [14].
Surface EMG signals recorded by a linear array or a matrix of channels have been extensively used to estimate IZ location [15], [16], [17], [18], [19]. In addition to classical IZ detection methods based on EMG amplitude [11], [14], [15], frequency [14], [15], and cross-correlation [15], [17] analysis, various algorithms have been developed for IZ detection from electrode arrays, such as the optical flow technique [20], robust linear analysis [21], Radon transform [22], and graph-cut segmentation [23]. Almost all existing IZ detection methods require processing of bipolar or single differential surface EMG signals derived from electrode arrays. Rodriguez-Falces developed an approach for IZ localization by processing monopolar signals based on the observation that in a monopolar montage the maximum MUAP with steepest rising phase is located near the IZ [24]. However, a method for localization of the IZ based on monopolar interference surface EMG signals has not been reported.
Principal component analysis (PCA) is a widely used method in EMG signal processing [25]. For example, PCA has been applied to improve muscle force estimation using monopolar surface EMG signals [26], [27], or to identify regional activation from monopolar surface EMG envelopes [28]. Of particular note, Laguna et al. [29] proved that the 2 nd principal component coefficients of PCA are proportional to the time delay of different channels. As MUAPs propagate in opposite directions from the IZ toward the muscle tendons, the signals from channels closest to the IZ are expected to have the shortest time delays. Based on this relationship, the current study aimed to develop a novel method to estimate the IZ location from monopolar surface EMG signals by analyzing the 2 nd principal component coefficients derived from PCA.
The performance of the proposed method was assessed with both simulation and experimental approaches. High density surface EMG signals in a monopolar electrode configuration were processed for IZ detection by analyzing the 2 nd principal component of PCA. For both simulated and experimental surface EMG signals, the performance of the 2 nd principal component analysis was compared with cross-correlation analysis [15], [17] and root-mean-squared (RMS) amplitude analysis [11], [14], [15], which were applied to bipolar signals. Both the simulation and experimental studies indicate that 2 nd principal component analysis provides an effective and robust way to estimate the IZ from monopolar high density surface EMG signals.

A. Framework
The muscle IZ estimation framework based on the 2 nd principal component is shown in Figure 1, where S = W T X. PCA projects standardized (zero mean and unit variance) EMG signals X (M-by-N matrix, N samples and M channels) from an original space of M channels to a new space of M components which are uncorrelated over the dataset. Each W column (coefficients matrix, M-by-M matrix) contains coefficients of all M channels for a corresponding principal component, and the columns are in descending order of component variance. Analysis of the 2 nd principal component coefficients can provide useful information for muscle IZ location. The rationale is provided below.

B. Rationale
PCA is a useful linear transformation technique that has been used in many EMG studies [25]. Briefly, for standardized electrode array surface EMG signals X, PCA is to perform the eigen decomposition on the correlation matrix R x , which is a M×M matrix with each element representing the covariance between two channels. Then the so-called principal components are the eigenvectors of the correlation matrix, which are sorted based on their corresponding eigenvalues. The first principal component is the eigenvector which has the largest eigenvalue. Here we used a simplified time misalignment data model to indicate that the 2 nd principal component coefficients are related with the time delay of each channel.
In an ideal condition, the M observed standardized signals x i (t) is modeled by equation (1), where s(t) is an unknown, deterministic signal, θ i is the time delay of the i th channel, v i (t) is independent, zero-mean, Gaussian white noise with variance For small θ i , the observed signal can be approximated by From [29], the correlation matrix R x , can be approximated by where n is the sample point, E s and E s are the energy of s and s , respectively. Further, the study [29] points out that the coefficients of the 2 nd principal component of x, denoted by ψ, is approximately proportional to θ , namely: where β is a proportionality factor. A standardized monopolar EMG signal can be approximately described by model (1). The time delays between different channels of EMG signal are very small. Therefore, it can be considered that the time delays θ satisfies equation (4). The channels located near the IZ should have minimum time delay. Based on the approximate proportional relationship between θ and ψ, we can determine the time delays θ and thus estimate IZ location once ψ is available through a routine operation of PCA. Note that although the sign of β is not derived in the model, we can use the derivative of the coefficients of each column to estimate the IZ. Both simulation and experimental approaches were used to test the validity of this strategy for IZ detection.

C. Performance Validation With Surface EMG Simulation
In the simulation study, the Fuglevand motor neuron pool model [30] and a surface EMG model [31] were followed. A motor neuron pool including 120 motor units was simulated. The distribution of the recruitment threshold (RTE) of each motor unit was modelled as an exponential function. Therefore, most of the motor units would be recruited at relatively lower excitation levels, and few had relatively high thresholds. Once the excitatory drive exceeded the recruitment threshold, the motor unit started to discharge at a minimum firing rate of 8 Hz. The last motor unit was recruited at 40% maximum excitation. The firing rate of each motor unit increased linearly until the peak firing rate was reached. In this study, all motor units followed the "onion skin" firing strategy, meaning that the peak firing rates of the later recruited large motor units were lower than the early recruited small ones. The peak firing rates of the largest and smallest motor units was 25 Hz and 35 Hz, respectively. The inter-spike interval of motor unit firing was modeled as a random process with a Gaussian probability distribution function.
The muscle was simulated in a cylindrical shape, with a radius of 8 mm. The thickness of fat plus skin layers was 2.5 mm. All 120 motor units innervated a total of 70,000 muscle fibers. The number of muscle fibers per motor unit was simulated to have a 100-fold range and followed an exponential distribution. A tripole model [31] was used to simulate generation of a MUAP. Briefly, two current tripoles originated at the neuromuscular junction, propagated in opposite directions, and terminated at the fiber-tendon endings. The monopolar signal detected by the electrode was the summation of the contribution from each of the tripoles, i.e., a MUAP was simulated as the sum of its constituent muscle fiber action potentials. In the simulation of a maximum voluntary contraction (MVC), surface EMG signals at a specific channel was generated as a sparse combination of MUAP trains from all 120 active motor units. A sampling rate of 2 kHz per channel and a duration of 10 seconds were used in the simulation.
The surface EMG signals were simulated to be recorded by a 40-channel surface array (5 columns by 8 rows, with a horizontal and vertical inter-channel distance of 5 mm). The channel columns were aligned parallel to the muscle fiber direction ( Fig. 2A). To verify the PCA results, 9 simulations with the IZ located at different row positions were investigated (Table I). For simplicity, the IZ location of each column of the array was simulated at the same row.
The effect of noise on the results was investigated under two situations. First, different signal to noise ratios (SNRs) of 5 dB, 10dB, 15 dB and 20 dB were simulated in each channel with additive zero-mean Gaussian noise (spatially independent). For each SNR level, we ran the simulation 20 times (i.e., 20 repetitions × 5 columns × 9 locations = 900 trials under each SNR level), to test the IZ estimation performance. Second, we simulated more or less contamination of a specific channel among a column of electrodes. As Figure 3 shows, the IZ was located at row 4, and the signals of row 3 was simulated to be more severely contaminated at 0 dB compared with other channels at 20dB. This situation is not unusual in high density EMG recordings [32]. The noise, if not processed, will affect associated bipolar channels, as shown in Fig. 3B.
D. Testing With Experimental Surface EMG 1) Participants and Consent: The experimental surface EMG signals were acquired from the biceps brachii muscle of 9 healthy subjects (28.9 ± 4.8 years, mean ± standard deviation). They were well informed of the experiment procedures and possible risks of the experiments, and gave written informed consent approved by the local ethics committee (ethical approval number: GWIRC-AF/SC-07/2016.20).

2) Experimental
Protocol: Two electrode arrays (ELSCH064NM2, Bioelettronica, Torino, Italy) were used in the study (Fig. 2B and 2C). Each electrode matrix consists of 64 channels with an 8 mm inter electrode distance (IED) arranged in a grid of 5 columns by 13 rows (one column consists of 12 channels and the other four of 13 channels). After skin preparation, one channel array was placed over the lateral side (Matrix 1) and the other over the medial side (Matrix 2) of biceps brachii parallel to the fiber direction and fixed with elastic straps. A ground electrode was placed at the elbow. Subjects sat in a height adjustable chair with the back fully against the backrest and the dominant forearm  (right side in all cases) on a custom-made force measuring device (Fig 2D). A load cell (CZL-3 T, Leitai, Bengbu, China) was used to measure elbow flexion force, which was displayed on a second monitor in front of the subject as real-time feedback. Each subject performed two to three MVCs and the largest value of the trials was used for further analysis. During the elbow flexor MVCs, subjects were strongly encouraged to give their best effort and to minimize extraneous movements. The surface EMG signals were recorded by a signal amplifier in monopolar configuration (EMG-USB2, 12-bit A/D converter, Bioelettronica, Torino, Italy) at a sampling frequency of 2048 Hz and a gain of 1000.

E. Signal Processing
For experimental surface EMG signals, the two most lateral and two most medial columns were excluded from signal processing because they tended to have relatively low signal quality. This left 6 columns (78 channels) for IZ estimation (Fig. 2B). For each signal a length of 2s was processed. Both the simulated and experimental signals were firstly standardized (zero mean and unit variance), and then decomposed by PCA in a monopolar configuration. The spline interpolation was applied to the 2 nd component coefficients along the rows to determine the IZ location for each column. For simulated surface EMG signals, the identified IZ from PCA analysis was compared to the IZ location used in simulation. For both simulated and experimental surface EMG signals, the identified IZ from PCA analysis was compared with that estimated from cross-correlation and RMS amplitude analysis of the bipolar signals constructed from the monopolar configuration [11], [14], [15], [17]. The spline interpolation was also applied to the correlation coefficient to determine IZ location. For the RMS-based method, RMS values of row(i)−row(i+1), and row(i)−row(i+2) were calculated (i=1,2,3…row number−2). The IZ was identified as at row(i+1) if row(i)−row(i+2) had the smallest value. The IZ location was identified at between row(i) and row(i+1) if row(i)−row(i+1) had the minimum value. All analyses were implemented using MATLAB R2020b. Figure 4 shows an example of muscle IZ estimation by analyzing the 2 nd principal component of PCA when the IZ was simulated at row 4 for each column of the simulated electrode array surface EMG signals. When the SNR was 20 dB, two principal components explained more than 90% of the data set variance (Fig. 4A). Figure 4B shows the coefficient distribution of the 2 nd principal component, indicating that the 4 th row had the smallest coefficient across columns, suggesting the location of the IZs. This was confirmed by spline interpolation of the coefficients for each column. Figure 4C shows an example of spline interpolation for column 3, where the minimum point of the interpolated line was closest to row 4, indicating the IZ location.

A. Simulation Results
Similar analysis was performed with a lower SNR. With SNR at 5dB, the first two principal components explained more than 70% of the dataset variance (compared to more than 90% with SNR at 20 dB) (Fig. 4D). With decreased SNR, the IZ can still be estimated from the distribution of the 2 nd principal component coefficients (Fig. 4E), as confirmed by spline interpolation analysis (Fig. 4F). Figure 5 shows a similar analysis for simulated electrode array surface EMG signals when IZ was located between two adjacent channels (row 4 and row 5 in this example). For different SNR levels, the first two principal components explained most of the variance (Fig. 5A, B, C, SNR = 20dB, Fig. 5D, E, F, SNR = 5dB). The distribution of the 2 nd principal component coefficients is also shown (Fig. 5B, E).  For each column, row 4 and row 5 had similar minimum coefficients, suggesting the IZ was located between these rows. This was confirmed by spline interpolation of the coefficients for each column. Figures 5C and 5F show examples of spline interpolation, where the minimum point of the interpolated line was between row 4 and row 5, indicating the IZ location.
For the first noise condition, the performance of different IZ estimation methods is summarized in Table II.
The overall accuracy of each method decreased with reduced SNR levels. The cross-correlation method using bipolar configuration achieved slightly better performance than the PCA based method using monopolar configuration. Only PCA and cross-correlation based methods achieved more than 80% accuracy when the SNR was reduced to 5dB. The mean difference in IZ location estimated from either PCA or correlation was small compared with the real IZ location for different SNR  levels. In contrast, there was a large discrepancy between the real IZ location and the estimated one when the SNR was low.
For the second noise condition, Fig. 6 shows representative performance of different IZ estimation methods, where the true IZ was located at row 4. The 2 nd principal component analysis accurately located the IZ at row 4, whereas the crosscorrelation method located the IZ to the second bipolar pair (correlation coefficient between the bipolar pair row2−row3 and row3−row4), corresponding to row 3 under monopolar configuration (i.e., 1 IED discrepancy). The RMS amplitude analysis located the IZ between row 5 and row 4. Therefore, the PCA based method was less sensitive to low quality of a specific channel, compared with the RMS and correlation methods for IZ estimation. Figure 7 shows an example muscle IZ estimation from experimental electrode array surface EMG signals. Consistent with simulation analysis, the two principal components explained more than 85% of the variance. Figure 7A shows the distribution of the 2 nd principal component coefficients. For each column the smallest coefficients were closest to row 5 or row 6. The spline interpolation of the coefficients for each column was used to determine IZ location. An example is shown in Figure 7D, where the minimum point of the interpolated line was closest to row 6, indicating the IZ location at row 6. For experimental surface EMG signals, the cross-correlation and RMS based IZ detection was also implemented for performance comparison. Figure 7B and 7C  demonstrate the distribution of correlation coefficients between adjacent bipolar signals, and the RMS amplitude of each bipolar channel. For most columns, the minimum coefficient was found between bipolar pair r5−r6 and r6−r7, suggesting the IZ was located approximately at row 6. The spline interpolation of the correlation coefficients for each column was also performed. An example is shown in Figure 7E, where the minimum point of the interpolated line was closest to bipolar pair r5−r6 and r6−r7 (the 5th bipolar pair). However, in most columns, the minimum RMS amplitude was found associated with row 12. An example is shown in Figure 7F, where the minimum point was at r12 − r11, misidentifying the IZ location between row 12 and row 11. An inconsistency in IZ location obtained from the 2 nd principal component and cross-correlation based methods was also noted in this example. For instance, according to the spline interpolation of the 2 nd principal component coefficients the IZ of column 2 was located around row 6 (Fig. 8A). However, the cross-correlation coefficients indicated that the IZ was at row 8 (Fig. 8B). To understand the inconsistency, we checked the raw monopolar surface EMG signals from column 2 (Fig. 8C). It was found that a segment of signal at row 8 was of poor quality, which caused a mis-location of the IZ using the cross-correlation based method. However, the analysis based on the 2 nd principal component was not sensitive to this artifact.

B. Experimental Results
IZ positions in each column were independently obtained for each subject using the 2 nd principal component, crosscorrelation and RMS based analyses. Fig. 9 shows an example of identified IZ for each subject at column 4 using different methods. The PCA and cross-correlation based methods achieved more consistent performance than RMS analysis. In total, 54 columns of experimental signals were processed. The mean difference in estimated IZ location between the PCA and correlation methods was 0.47±0.4 IED (inter-electrode distance = 8 mm).

IV. DISCUSSION
We developed a novel muscle IZ estimation method from monopolar high density surface EMG signals, by analyzing the 2 nd principal component coefficients. PCA is an efficient and automated method that renders a statistical description of complex data by revealing hidden structures. Many studies have applied PCA to process multi-channel or high-density surface EMG signals for examination of muscle activation and intrinsic activation patterns [25], [26], [27], [28]. In principle, PCA can characterize components which account for a large amount of the total data variance, representing a common temporal pattern across multiple surface EMG channels. The rationale of applying the 2 nd principal component for IZ estimation was proven in a simplified time misaligned data model. The model shows that the 2 nd principal component coefficients are linearly related with the time delay of different channels [29]. Rodriguez-Falces reported that the IZ location corresponds to the electrode which recorded the highest slope of the rising phase of the monopolar MUAP [24]. In PCA analysis, channels with smaller 2 nd principal component coefficients could be expected to have steeper rates of rise and shorter latencies.
The linear relationship has been used to estimate the time delay of repetitive biomedical signals [29]; they achieved similar results in a more efficient way compared with the conventional maximum likelihood method. Time delay estimation is important for calculating muscle fiber conduction velocity. Conduction velocity of EMG signals is an important indicator of muscle integrity and has been examined in studies of muscle structure [33], fatigue [34], and pathology [35], [36]. Note that in the time misalignment data model, the signal was assumed to have fixed amplitude and shape across all channels. For the simulation and experimental high density surface EMG signals, although signal amplitude and shape varied considerably over time, the 2 nd principal component coefficients still performed well for IZ estimation.
While most of the previous IZ estimation methods using high density surface EMG required processing of bipolar signals, the method reported here is based on a monopolar montage. For experimental high density surface EMG signals, since the IZ was not known a priori, the cross-correlation and RMS methods using bipolar signals were implemented for a comparison with the 2 nd principal component-based approach. High agreement was achieved between the PCA and correlation-based methods. The mean difference between the two approaches was found to be less than one-half IED. There are several factors that may collectively contribute to the difference between the two methods. The cross-correlation based method requires processing of bipolar surface EMG signals for IZ detection, which is less representative of the muscle compared with monopolar signals used by the 2 nd principal component-based method. However, the monopolar montage is more vulnerable to crosstalk from neighboring muscles, especially during strong contractions. In addition, the effect of noise on estimation of the IZ may differ between the two methods. For example, as demonstrated in this study, if one channel in the bipolar configuration is of poor quality, then differentiation will induce artifacts in the adjacent channels leading to errors in IZ location. In contrast, the 2 nd principal component-based method is less sensitive to this kind of noise. Among the three tested methods, the performance of the RMS method is not as good as the other two since surface EMG amplitude can be affected by various factors such as noise and artifacts. The proximal/distal electrodes near tendons may induce misidentification, as previously discussed [15].
Although promising results were achieved for IZ estimation with the 2 nd principal component-based approach, limitations of the current study need to be acknowledged. For both simulation and experimental analyses, the IZ was simplified as a narrow band perpendicular to muscle fibers. However, IZ may have more complex shapes and may differ between muscles [9]. There are reports that some muscles may have multiple IZs [1], [37], [38]. With multiple IZs, the method applied here might not work well since signals originating from different IZs would overlap causing the IZ location to be misidentified. Indeed, this is also a limitation in most previous IZ estimation methods [39]. In the current study, the IED of the electrode array was 8 mm, corresponding to a spatial resolution of 4 mm for IZ detection. A smaller IED is required in order to increase the resolution.
In summary, we developed a novel high density surface EMG processing method to estimate a muscle's IZ location using monopolar high density surface EMG, by analyzing the coefficient distribution of the 2nd principal component of PCA. The performance of the developed method was validated by both simulation and experimental approaches. A primary feature of the developed method is characterization of IZ location based on monopolar signal processing. It can be implemented automatically, and provide a useful tool for muscle IZ estimation in both research and clinical laboratories.