Transcutaneous innervation zone imaging from high-density surface electromyography recordings

Objective. It is of great value to accurately localize innervation zones (IZs) to better diagnose and treat neuromuscular diseases, but it is challenging to do so noninvasively from surface electromyography (sEMG) recordings because of the blurring/distorting effect of the low conductive fat tissues. This study aimed to develop an innovative transcutaneous IZ imaging (TIZI) technique to precisely and efficiently localize the IZ distribution directly over the muscle surface in vivo from high-density sEMG recordings (HD-sEMG). Approach. The TIZI technique was implemented by incorporating HD-sEMG recording, signal decomposition, finite element analysis and inverse calculation. The performance of TIZI was evaluated on the flexor digitorum superficialis (FDS) muscle with simulated sEMG signal and experimental signal recorded from both healthy (n  =  3) and stroke participants (n  =  4). The accuracy of imaging was validated by both of the Pearson correlation coefficient (PCC) and localization error (LE) between the TIZI results and the ‘true’ IZ distribution. Main results. In the simulation study, results have shown PCCs of 99.85%  ±  0.11%, 99.79%  ±  0.08%, 99.63%  ±  0.22% and 99.31%  ±  0.54% at the depth of 10, 15, 20 and 25 mm and SNR of 25 dB. PCCs of 98.74%  ±  1.78% and 97.82%  ±  1.20% were respectively obtained for experimental signals acquired from the healthy and spastic FDS muscles. The TIZI provided smaller LEs of 1.4  ±  0.92 mm and 2.02  ±  1.3 mm, compared to LEs of 7.42  ±  2.29 mm and 7.8  ±  1.77 mm from skin observations in healthy and spastic FDS, respectively. Significance. Results have demonstrated the high performance of the proposed TIZI technique by transcutaneous imaging of the IZ distribution of the skeletal muscles. The performance improvement can be attributed to the elimination of the blurring/distorting effect caused by the low conductive fat and high conductive skin tissues. TIZI may provide an advanced neurological tool for the clinical treatment of neuromuscular diseases, such as guiding botulinum neurotoxin injections in spasticity management.

observed, such as the formation of complex endplates through denervation and reinnervation, small angulated fibers and fiber-type grouping [15,16]. These factors working together or alone can impose significant challenges on the IZ localization for the guidance of the BoNT injection.
With significant advances in surface electromyogram (sEMG) signal recording and processing techniques, high-density sEMG (HD-sEMG) array has been employed to non-invasively identify the IZ surface location by recognizing the reversal of signal polarity with a differential setup [17,18]. However, the volume conductor imposes a blurring/distorting effect on the surface potential distribution due to poorly conducting subcutaneous fat and well conducting skin [19,20]. This is analogous to the blurring effect existing in electroencephalogram (EEG) caused by the low conductive skull and relatively high conductive scalp [21]. This distortion effect may lead to an erroneous estimation of surface IZ localization. A three-dimensional innervation zone imaging (3DIZI) technique [22] has been previously proposed to tackle this issue by imaging the IZ in the 3D space of muscles. Its effectiveness has been fully validated based on synthetic and experimental sEMG signals from both healthy and spastic muscles [23][24][25][26]. However, the implementation of the 3DIZI requires a full 3D computational model of the body part that contains the target muscles, bones and all surrounding tissues, thereby incurring substanti al computational effort and consequently restraining its clinical applicability. Clearly, there is an unmet need for technologies that can characterize the IZ distribution in muscles with higher efficiency and maintained accuracy.
In order to bridge this gap, a transcutaneous IZ imaging (TIZI) technique was proposed in this study to image the IZ distribution transcutaneously over the muscle surface. The TIZI requires a simplified model which consists of only two components, i.e. (1) the skin and (2) fat tissues overlaying the target muscles, purposely excluding all the muscles, tendons, compact bone, cancellous bone and other tissues, and as such greatly reduces the model complexity and computational effort. The performance of TIZI was validated with both computer simulations, as well as experimental data from healthy and spastic flexor digitorum superficialis (FDS) muscles. FDS muscle was chosen because it is one of the most frequently used muscles responsible for finger flexion and the clenching of fists [27] and one of the most frequently BoNT injected upper-limb muscles for the treatment of focal spasticity in the post-stroke survivors [28].

Participants
Three healthy subjects (one female, mean age 35.7 ± 4.9 years) with no reported history of neuromuscular diseases, and four stroke subjects (two females, mean age 64.3 ± 8.3 years) with upper limb spasticity were enrolled in this study. The average duration of the chronic stroke is 8.5 ± 3.5 years (three ischemic and one hemorrhagic). All subjects were well informed of the risk of the study and gave informed consent. The experiment was approved by the University of Houston and University of Texas Health Science Center-Houston institutional review board. The FDS muscles of the dominant and paretic side were investigated for the healthy and stroke subjects, respectively.

Experiment configuration
The subjects were seated comfortably on a heightadjustable experiment chair with the tested arm held in a neutral position. After the skin preparation, two high-density (8 by 8 sensors) sEMG grids (TMSi, Enschede, The Netherlands) were placed adjacently to form an 8 by 16 (muscle fiber direction) electrode array, as presented by figure 1(c). The longitudinal midline of the combined electrode grids was placed over the estimated direction of FDS muscle fibers, following the recommendation of Norizan et al [27]. The diameter of each individual electrode is 4.5 mm and the interelectrode distance is 8.5 mm. The reference electrode was placed on the lateral epicondyle of the humerus and a ground electrode was attached to the wrist of the contralateral arm using a fully soaked Velcro wrist band (TMSi, Enschede, the Netherlands). The subjects were instructed to perform isometric finger flexions at the proximal interphalangeal joints (PIP) against an adjustable metal plate with torque sensor (TRS-500, Transducer Techniques, Temecula, CA, USA). Three 5 s trials of maximum voluntary contraction (MVC) separated by 2 min of rest were performed. The maximal force value of three attempts was selected as MVC of the FDS muscle to guide submaximal contractions. Then three trials of isometric finger flexion at 20% MVC were performed. Each trial lasted 10 s. The 128-channel monopolar sEMG data was recorded via a Refa-136 amplifier (TMSi, Enschede, the Netherlands) at a sampling rate of 2048 Hz with a 24-bit resolution.

Signal processing
Baseline HD-sEMG signals were examined by evaluating the root mean squared value and the corrupted channels were identified and excluded from further analysis. Raw HD-sEMG data were notch filtered at 60 Hz and bandpass filtered at 10-500 Hz with second order Butterworth filters. A K-means clustering modified convolution kernel compensation (KmCKC) algorithm, recently developed in our lab, was employed to decompose the HD-sEMG signals into constituent pulse trains [29]. In brief, a time instant n 0 with significant muscle activity was selected based on the global activity index as defined in [30]. Then ten firing instants of the motor units (MUs) that were simultaneously active at n 0 were selected.
Using a K-means clustering method, these firing instants were clustered into two groups. The group with the largest number of firing instants was used to construct the initial innervation pulse train (IPT) by using a linear minimum mean square error method. Finally, a modified multi-step iterative convolution kernel compensation (CKC) method [30] was adopted to update the estimated IPT until desirable decomposition accuracy was achieved. The motor unit action potential (MUAP) was constructed by using the spike-triggered averaging method [31]. Since the MUAP generates from the neuromuscular junction (indicated as IZ) and propagates in two opposite directions towards the fiber endings [32], the surface position of the IZ of a particular MU can be localized from the bipolar map of MUAP by checking the reverse phase of the propagating signals. The 128-channel sEMG measurement corresponding to the time instant that the electrode closest to the identified IZ achieves its maximal absolute values was used to solve the inverse problem using the proposed TIZI.

Modeling
A template full forearm model, consisting of FDS, flexor carpi ulnaris, flexor carpi radialis, palmaris longus, flexor digitorum profundus, flexor pollicis longus, extensor muscle, skin, fat and bone, was constructed from the general magnetic resonance image data set of the second healthy subject by following the same procedure as our previous work [22]. To achieve a subject-specific TIZI calculation, the model was adjusted specifically for each subject based on the arm size and ultrasound images acquired with a portable B-mode ultrasound scanner (M Turbo, SonoSite, Bothell, WA), by taking into consideration the thickness of fat tissue. The adjusted model was then meshed into finite element (FE) models with 336 256 tetrahedral elements and 62 052 nodes, using Abaqus 6.12 (SIMULIA, Providence, RI), as shown in figures 1(a) and (b). For TIZI calculation, only the skin and fat geometry were used which contains 190 756 elements and 38 109 nodes; and the full model was established for the purpose of validation. Isotropic conductivity values of 1.00, 0.05 and 0.02 S m −1 were assigned to skin, fat and bone, respectively, and anisotropic conductivity values of 0.1 and 0.5 S m −1 (anisotropy ratio of 5) were assigned to muscle in transverse and muscle fiber direction [19,33,34]. A total of 9390 current dipoles were evenly distributed throughout the 3D space of the FDS, flexor carpi ulnaris and palmaris longus muscles with a spatial resolution of 2 × 2 × 2 mm 3 .

Inverse mapping
The electric potential distribution ∅ inside of the skin-fat volume conductor Ω can be described by the following Laplace's equation [35].
where σ is conductivity. By using the FE method, differential equation (1) can be expressed as: The FE stiffness matrix K was calculated based on our previous study [36]. The nodes in the FE model were partitioned into three parts: the 128 electrode nodes on the skin surface denoted by subscript s, the nodes on the outer surface of muscle denoted by subscript m (the inner surface of fat is equivalent to the outer surface of muscle) and the rest nodes inside the skin-fat volume conductor denoted by subscript r. By breaking down the K matrix accordingly, equation (2) can be written in the partitioned form as:  The linear relationship A · ∅ m = ∅ s can be obtained by eliminating φ r from equation (3), where is the transformation matrix between the ∅ s and ∅ m [35]. Since the number of surface measurements is much smaller than the number of nodes on the muscle surface, the equation A · ∅ m = ∅ s represents an underdetermined problem. The weighed minimumnorm regularization was employed to solve the ∅ m based on the measured ∅ s by minimizing the following equation [37]: where the first term of the above objective function is the residual error and the second term is the regularization term representing the power of the solution. The ∅ m can be estimated as: where W is a diagonal weighted matrix and its ith element, W i , can be calculated as the norm of the ith column of transformation matrix A. The regularization parameter λ was determined as the point with the maximum curvature on the L-shaped curve describing the relationship between the two items in equation (4). For more details about the L-curve method please refer to [38].

Computer simulation
A series of computer simulations was conducted to validate the performance of the TIZI. A group of dipoles located inside of an ellipsoid with sizes of 5 and 4 mm in parallel and normal to the muscle fiber direction of FDS were activated to simulate the IZ [39,40]. In this study, IZs were set at five locations along the muscle fiber direction and four depths of 10, 15, 20 and 25 mm from the skin surface.
The noise-free synthetic surface measurements were generated by solving the forward problem with the assumed IZ and then different levels of Gaussian white noise (GWN) [41] were added to mimic the experimental, noise-contaminated signals. Phinyomark et al reported that the signal-to-noise ratio (SNR) of the sEMG signal acquired from the forearm and upper arm was in the range of 24.89~32.06 dB [42]; therefore the GWNs were added with a mean of zero and SNRs of 15, 20, 25 and 30 dB. The SNR was defined as 10 log 10 ( [43,44], where σ s and σ n are the respective standard deviations of the synthetic noisefree measurement and the added noise. Fifty sets of independent noise with a mean of zero and standard deviation of σ n were generated and added to the synthetic noise-free measurement.

Simulation
For simulation, the true potential distribution over the muscle surfaces, denoted by P t , was calculated by solving the forward problem based on the defined IZs. To further determine if the TIZI can tackle the blurring/distorting issue, the simulated sEMG were also projected onto the muscle surface, denoted by P p , using the weighted average of the k-nearest neighbor method with weights defined by Shepard's method [18]. The reconstructed muscle surface potential of TIZI (P r ), along with P t and P p , were normalized to their own maximal absolute values. A uniform threshold was applied to these three potential distributions for the subsequent comparison. The performance of TIZI was evaluated by the Pearson correlation coefficient (PCC) between P p and P t (PCC pt ) and that between P r and P t (PCC rt ), as well as by the localization error (LE) defined as the Euclidean distances between the centers of the P p and P t (LE pt ) and that between P r and P t (LE rt ). PCC pt and PCC rt were calculated based on equations (6) and (7) [45,46].
where P p , P t and P r are the mean values of P p , P t and P r , respectively, M denotes the number of muscle surface nodes at which either P p or P t is higher than the preset threshold and N denotes the number of muscle surface nodes at which either P r or P t is higher than the preset threshold. The mean and standard deviation of PCC rt were calculated based on the reconstruction results of 250 (5 axial locations × 50 independent noise per location) noise-contaminated measurements for each simulated IZ depth where five different positions were considered.

Experiment validation
For experimental validation, the potential distribution over the muscle surface was obtained using the wellestablished and validated 3DIZI approach [23][24][25][26]. Briefly, 3DIZI was first employed to reconstruct the strength of the distributed current dipoles from HD-sEMG recordings. The WMN was employed to solve the linear equations L · X = ∅ s describing the relationship between the strength of the distributed source dipoles X and the HD-sEMG measurement ∅ s , where the L represents the lead filed matrix. The 'true' potential distribution over the muscle surface (P t ) was finally calculated by solving the forward problem with the reconstructed dipole intensities. Similarly, the PCC pt , PCC rt , LE pt and LE rt were utilized for performance evaluation.  Figure 2 summaries the mean and standard deviation of the PCC rt s corresponding to five noise levels and four IZ depths. High PCC rt s of 99.85% ± 0.11%, 99.79% ± 0.08%, 99.63% ± 0.22% and 99.31% ± 0.54% were respectively achieved with IZ depths of 10, 15, 20, 25mm and a SNR of 25 dB, which is very close to the SNR of 25.81dB observed in the sEMG of the forearm [42]. With an increasing noise level, the average PCC rt decreased and variation of PCC rt increased, respectively. However, even with a relatively high noise level (a SNR of 15 dB), satisfactory PCC rt s of 99.72% ± 0.22%, 99.57% ± 0.29%, 99.53% ± 0.24% and 99.0% ± 0.71% can be achieved for simulated IZs at the depth of 10, 15, 20, 25mm, respectively. Figure 3 shows a typical computer simulation result of a pre-defined IZ located in the FDS (15 mm deep from the skin surface) with a SNR of 25 dB. With a threshold of −0.75, a PCC rt of 99.73% was achieved, while PCC pt was as low as 87.56%. The significant blurring/distorting effect caused by the fat and skin was manifested in two aspects: (1) a larger blue region observed on P p (under given threshold) compared to that was found on P t ; (2) a larger LE pt of 14.2 mm compared to a 0.2 mm of LE rt .

Experiment results
After the sEMG signal decomposition by rejecting the MUs with pulse-to-noise-ratio lower than 30 dB [47], the coefficient of variance of the inter spike interval greater than 0.3 [48,49], one-sided MUAP propagation and symmetrical distribution of MUAPs without propagation [18], 6.5 ± 2.1 and 5.2 ± 2.0 MUs were yielded from the three healthy and four stroke subjects, respectively. Further, the MUs not belonging to FDS determined by the 3DIZI method were excluded from the subsequent TIZI. As a result, in total 8 (2.8 ± 1.5 per subject) and 7 (1.8 ± 0.5 per subject) MUs were yielded for three healthy and four stroke subjects, and their average depths were 12.21 ± 2.63 and 12.68 ± 3.04 mm, respectively. PCC rt s of 98.74% ± 1.78% and 97.82% ± 1.20% were respectively achieved for the healthy and spastic FDS muscles. For healthy subjects, TIZI reduced the LE from 7.42 ± 2.29 mm (LE pt ) to 1.4 ± 0.92 mm (LE rt ); for stroke subjects, the LE was reduced from 7.8 ± 1.77 mm (LE pt ) to 2.02 ± 1.3 mm (LE rt ). Figure 4 shows the experiment results from one representative healthy subject with LE pt of 9.8 mm and LE rt of 1.7 mm. Figures 5 and 6 represent the TIZI results from two stroke subjects. In one subject (figure 5), the TIZI reduced the LE from 6.8 mm (LE pt ) to 0.5 mm (LE rt ). Compared to figure (a) in figures 4 and 5, the blue part of the HD-sEMG recording presented in figure 6(a) lacks a clear boundary. This blurring sEMG recording map, together with the high similarity of bipolar action potentials in each column in figure 6(b), indicate a deeper MU (20.71 mm). According to the simulation results, under this circumstance, the PCC rt was usually relatively low (99.85% ± 0.11% (10 mm) versus 99.31% ± 0.54% (25 mm) with a SNR of 25 dB). However, a PCC rt of 97.63% was achieved by the TIZI and reduced the LE from 9.8 mm (LE pt ) to 2.6 mm (LE rt ).

Discussion
The commonly observed spasticity after stroke is one of the major sources of disability in post-stroke survivors. The aim of this study was to propose an accurate and effective IZ characterization approach to guide the BoNT injection for the treatment of poststroke spasticity. Results have shown a clear blurring/ distorting effect from the skin surface measurements, and therefore surface observations sometime might be unreliable and should not be used directly for the IZ estimation. The proposed TIZI, on the other hand, was proven able to successfully reduce this disturbance and provide an accurate estimation of muscle IZs. In a representative example presented in figure 3, the TIZI technique improved the PCC from 87.56% (PCC pt ) to 99.73% (PCC rt ) and reduced the distance error from 14.2 (LE pt ) to 0.2 mm (LE rt ). Similarly, the LE of 8 MUs was reduced from 7.42 ± 2.29 mm (LE pt ) to 1.4 ± 0.92 mm (LE rt ) using TIZI in healthy FDS muscles.
Drastic changes in muscle properties may occur after stroke, including muscle atrophy, reduced muscle fascicle length, increased fat tissue, as well as disordered control of motor units and compromised motor unit activation [25]. These functional and anatomical alterations in muscle properties together with the morphological change of neuromuscular junction may increase the difficulty in HD-sEMG decomposition and dramatically affect the accuracy of the TIZI. It is therefore necessary to examine the robustness of TIZI performance in spastic muscles. A PCC rt of 97.82% ± 1.20% and a LE rt of 2.02 ± 1.3 mm were obtained. Though these results are slightly less precise compared to those under healthy conditions, a nonparametric Mann-Whitney test for PCC rt and LE rt of the healthy and spastic FDS muscles demonstrates that no significant difference was observed between the two groups, suggesting that the imaging performance can be reliably obtained in both healthy and spastic muscles. Meanwhile, we also checked the IZ location and spread, yet no further difference was observed. The lack of discernable difference may be due to the small sample size of this study.
Though our previously developed 3DIZI can address this blurring/distorting issue by fully considering the electrical properties of the fat and skin during the process of constructing the FE stiffness matrix K, the substantial load required for this approach dramatically limits its application in practice. The proposed TIZI has shown the comparable accuracy of 3DIZI but with higher efficiency in IZ imaging. This improvement can be largely attributed to the simplified model which only contains the skin and fat tissues. The load required by 3DIZI to construct a full forearm model consisting of several layers of small muscles in close proximity [50] is intense, as it requires the segmentation, construction of 3D closed surface and smoothing for each anatomical part. By comparing the two FE models required by 3DIZI and TIZI presented at the beginning of the Result section, the modeling effort required by TIZI is about one fourth of that required by 3DIZI by estimation. The efficiency improvement also attributed to the reduced load to calculate the matrix required by inverse problem solving. For the 3DIZI approach, the time required by lead field calculation mainly depends on the number of dipoles defined, besides the numbers of elements and nodes in the FEM. The forward problem needs to be solved the number of defined dipoles times. As a result, it took 16.03 ± 0.28 h to calculate the lead field matrix L required by solving the inverse problem of 3DIZI. In contrast, for the TIZI, we only need to calculate the transformation matrix A once, and the calculation time depends entirely on the numbers of elements and notes in the FEM. It took only 0.53 ± 0.03 h to calculate the transformation matrix A required by solving the inverse problem of TIZI. In the current study, the TIZI calculation was performed on the adjusted template model corre sponding to each subject's arm size and ultrasound images. By taking advantage of the dramatically improved modeling efficiency of TIZI, a real subject-specific model constructed directly from each subject's MRI images may markedly improve the imaging accuracy. In addition, the effects of subcutaneous fat thickness on the imaging outcome should be analyzed in the future study. The high efficiency of TIZI is expected to benefit the fast and precise IZ estimation in patients and enable clinical friendly application in cases such as guiding BoNT injection.
The crosstalk effect due to a combination of high frequency, propagating and non-propagating signal components [51,52] might exist in the forearm. As such, EMG signals recorded from the electrodes placed on top of the FDS may be contaminated by multiple myoelectric sources in close proximity, which may consequently lead to an erroneous estimation of the signal origins, i.e. the IZs. However, this problem was tackled in the current study by integrating sEMG signal decomposition and 3DIZI imaging. As the IZ of each motor unit is expected to point to an anatomical location that belongs to only one muscle, we applied the 3DIZI method to each decomposed MU, identified the corresponding IZ, and excluded the MUs which fell outside of the FDS space from the TIZI calculation.
As there are currently no in vivo tools available that can reliably measure the 3D muscle activities, the TIZI results of the current study were compared to those of the 3DIZI, which has been well validated using intramuscular recordings in both healthy and spastic muscles [23][24][25]. Although proven accurate and reliable, it should be noted that the 3DIZI results are only an estimation of IZ location rather than its real 3D location. In the future study, the performance of TIZI should be further validated by comparison with the muscle surface potential distribution generated by dipoles located inside of the IZ detected with the intramuscular wire electrode.
It is noteworthy that despite its proven accuracy, TIZI does not provide an accurate estimation about the MU depth, which can be obtained from 3DIZI. However, in cases of superficial or thin muscles, such as the FDS in the forearm, TIZI is expected to work as accurately as 3DIZI with remarkably reduced modeling and computation burdens.
In this study, the KmCKC algorithm was employed to decompose the HD-sEMG data into constituent pulse trains prior to the TIZI reconstruction; therefore, the final imaging results largely depended on the decomposition performance. Due to the signal superimposition and cancellation at higher contrac-tion levels, the KmCKC algorithm works better with low contraction forces, i.e. 10%-30% MVC [29,30]. We acknowledge that at this contraction level, only a portion of FDS motor units are activated; yet this portion is expected to be a representative sample of whole motor unit pool. Moreover, the KmCKC has been successfully applied to decompose force with a contraction level of 50% MVC [18] and the CKC and the progressive fastICA peel-off algorithm even used to decompose the sEMG signal with a contraction level of 100% MVC [53]. Therefore, it is reasonable to expect the generalization of TIZI to higher contraction levels. Moreover, due to the requirement of voluntary contractions, the TIZI method may not be applicable to the patient cohorts who are not able to exert force or have difficulty in maintaining the isometric contraction. Nonetheless, the M-wave induced by maximally stimulating the musculocutaneous nerve does not require voluntary contraction and decomposition. Therefore, based on M-wave data, imaging the global IZ distribution by using the TIZI would be a great choice for this type of patient cohorts and is under our consideration for future study.

Conclusion
In this study, a novel TIZI technique was developed to image the IZ distribution transcutaneously over the muscle surface from HD-sEMG recording with high accuracy and computation efficiency. This approach was validated with both simulated sEMG signals as well as with HD-sEMG signals recorded from three healthy and four spastic FDS muscles. The high accuracy and efficiency for modeling and computation make it a promising and clinically friendly technique for the diagnosis and management of neuromuscular diseases, such as guiding BoNT injections.

ORCID iDs
Ping Zhou https://orcid.org/0000-0002-4394-2677 Yingchun Zhang https://orcid.org/0000-0002-1927-4103 (c) Projected muscle surface potential distribution P p from the HD-sEMG recordings. (c) 'True' muscle surface potential distribution P t calculated based on 3DIZI results. (d) Reconstructed muscle surface potential distribution P r by using TIZI. Missing circles in (a) indicate the reference channel and the discarded channel due to the bad signal quality. Red dots in (c)-(e) represent the geometric centers of the blue potential circles.