A Method for Identification of Mechanical Response of Motor Units in Skeletal Muscle Voluntary Contractions Using Ultrafast Ultrasound Imaging—Simulations and Experimental Tests

The central nervous system coordinates movement through forces generated by motor units (MUs) in skeletal muscles. To analyze MUs function is essential in sports, rehabilitation medicine applications, and neuromuscular diagnostics. The MUs and their function are studied using electromyography. Typically, these methods study only a small muscle volume (1 mm3) or only a superficial (<1 cm) volume of the muscle. Here we introduce a method to identify so-called mechanical units, i.e., the mechanical response of electrically active MUs, in the whole muscle ( $4\times 4$ cm, cross-sectional) under voluntary contractions by ultrafast ultrasound imaging and spatiotemporal decomposition. We evaluate the performance of the method by simulation of active MUs’ mechanical response under weak contractions. We further test the experimental feasibility on eight healthy subjects. We show the existence of mechanical units that contribute to the tissue dynamics in the biceps brachii at low force levels and that these units are similar to MUs described by electromyography with respect to the number of units, territory sizes, and firing rates. This study introduces a new potential neuromuscular functional imaging method, which could be used to study a variety of questions on muscle physiology that previously were difficult or not possible to address.


I. INTRODUCTION
The central nervous system controls force production of the skeletal muscles in a quantal manner by successive recruitment of motor units (MUs) [1]. The MU consists of a motor neuron and a few tens to some hundred [2] innervated muscle fibers within a delimited territory (Fig. 1). Upon neural activation, the muscle fibers are electrically depolarized and give rise to a contraction [3] that results in force production through shortening and thickening of the fibers (mechanical twitch), which is a subtle (micrometer amplitude) and tran-The associate editor coordinating the review of this manuscript and approving it for publication was Giulia Matrone . sient phenomenon (millisecond scale) (Fig. 1). To study the MUs' function is central in the understanding of the skeletal muscle physiology in sports medicine applications and the diagnosis of neuromuscular diseases [4], [5]. The analysis is based on recording and extracting activation characteristics of multiple individual MUs and allows access to the neural input from the spinal cord and the muscle fibers' condition.
Today, electromyography (EMG) is the gold standard to study MUs, and there are multiple techniques available. Invasive needle EMG typically samples electrical potentials at the tip of a needle with a sampling volume of about one mm 3 [6] but provides no quantitative depth or spatial information. Scanning EMG is a spatial extension of the needle EMG The motor unit comprise a bundle of muscle fibers located within a territory and they are controlled by a motor neuron coordinating the activation by a firing pattern. Upon activation, the fibers of a MU contract causing a mechanical shortening and thickening of the fibers within that territory (twitch). In the cross-sectional plane (perpendicular to fiber orientation) this contraction will correspond to a radial expansion of the territory.
technique, which gives quantitative information about the MU territory size along the direction of the needle (1-D) [7] and only rough information about the depth. High-density surface EMG uses electrode arrays or matrices on the skin's surface to record the activity and decompose signals into individual MUs [8], [9] within measurement depths up to 1 cm but with poor depth selectivity [10]. Thus, the existing bioelectrical techniques have limitations in their field of view. Ultrasound imaging is a non-invasive technique allowing mechanical information from a large field of view [11]. In skeletal muscle applications, ultrasound has mainly been used for visualization of structural features and detection of contraction onset [12]- [14]. Recently, non-invasive high frame rate (>2000 images per sec) high-resolution ultrasound imaging has been introduced [15] and successfully applied in several medical applications, e.g., in imaging of the brain [16]- [18] and imaging of MUs during externally controlled electro-stimulations [19], [20]. To stimulate single MUs electrically is difficult since the axons' excitability of a given nerve may be similar [21], [22]. Also, information about neural coordination and recruitment cannot be obtained, which is crucial in the study of voluntary contractions. However, there are no methods available that can identify individual MUs under voluntary contractions in ultrasound image sequences. The challenge is the highly complex physiology involved in muscle activation, where tens to hundreds of MUs with overlapping MU territories are active simultaneously with complex firing patterns coordinated by the nervous system.
The rationale for using ultrasound imaging to study MUs is that the mechanical contraction of the fibers of an MU is linked to the electrical depolarization via the excitation-contraction mechanism [23]. In theory, that information on the neural input to the MUs (firing pattern) and their territory size is accessible in both bioelectrical and mechanical domains. In contrast, the bioelectrical and mechanical responses of a contracting MU will differ.
The purpose of this study was to develop a method that identifies the mechanical response of individual active MUs (which we define as mechanical units) in the muscle under voluntary contraction. There were two specific aims. The first aim was to evaluate the method's theoretical performance to identify MUs' neural firing patterns and territories based on simulated data. The second aim was to assess the feasibility of the method on experimental data of voluntary contractions by comparing characteristics of the identified mechanical units with those of MUs described by EMG. We present a spatiotemporal decomposition framework based on ultrafast ultrasound and tissue velocity imaging sequences.

A. DEFINITIONS
Here we define a mechanical unit as the spatio-temporal mechanical response of an individual active MU comprising a territory and a twitch. A twitch is here the mechanical thickening of MU' fibers after one depolarization, i.e., the EMG analog to the MU action potential. A twitch train is the fusion of repeated twitches caused by the repeated depolarizations of an MU firing pattern, i.e., the EMG analog to MU action potential train.

B. THE PROPOSED METHOD
The identification of mechanical units is based on 2-D crosssectional velocity image sequence data of skeletal muscle tissue, capturing the thickening of fibers during the contraction phase (Fig. 1). We assume no mechanical connectivity between the fibers of different MUs, and thus the force is transmitted along the fiber direction only. First, the dimension of the data is reduced before separating components with different spatiotemporal characteristics. Next, a unit selection step is applied where we separate the sources that are probable mechanical units from noise. Finally, each unit's spatial territory and firing pattern were estimated. See Fig. 2 for an overview of the method.

1) DIMENSION REDUCTION
The 2-D tissue velocity image sequence can be viewed as a 3-D data matrix that can be transformed by vectorization of the images to a 2-D matrix Y , having dimension n s × n t , where n s = n x × n y , and n x and n y denote the number of pixels for the width and height, and n t is the number of images (time samples). In this work, n x = n y = 128 pixels (corresponding to 40 mm, i.e., 0.3 mm/pixel) and n t = 4000 (i.e., two seconds data at 2000 Hz image rate). To reduce the computational complexity of the decomposition, we reduced the dimension of Y through Singular Value Decomposition (SVD) [24]. We constructed a new Y based on the first k = 100 components of the SVD. This choice is arbitrary, but k needs to be much larger than the expected number of active units in a contraction and will be the upper limit of how many 50300 VOLUME 8, 2020 units that can be identified by our method. k = 100 retained (empirically) ∼99% of the total variance.

2) SPATIOTEMPORAL DECOMPOSITION
The separation is based on spatiotemporal independent component analysis (stICA) [25]. ICA has been widely used in EEG data analysis [26], [27] and in particular, for decomposing electrical sources from MUs in EMG [28]- [30], which is a similar problem to the present one with the recorded signals being essentially linear combinations of the sources. Here, we assume that each component in Y is a linear combination of the unit sources [31]. We use the approach described in [25] to formulate and solve the stICA problem: Y = S T where S is the vectorized spatial matrix with dimension n s × k, T is the temporal matrix with dimension n t × k, and is a diagonal scaling matrix required to ensure that S and T have amplitudes appropriate to their respective cumulative distribution functions. To estimate S and T we use a kurtosisbased cost-function [25], mainly motivated by an MU's mechanical contraction should mainly take place in a small localized region and nowhere else in the image causing the spatial maps to have a sparse spatial distribution of velocities.
The algorithm includes a ratio parameter (α) defining emphasis on spatial sparsity vs. temporal sparsity, which was set to 0.8 based on the empirical evaluation (data not shown, see Discussion).

3) UNIT SELECTION
The output from the stICA gives us k = 100 sparse spatial maps and corresponding twitch trains because we selected that number in the dimension reduction step. To discriminate the probable mechanical units from noise we apply a similar technique as in [28]. First, we extract two features from each of the k estimated twitch trains; kurtosis and skewness of the power spectral density (PSD). A twitch train of a MU should have the majority of the energy within the range 5 to 15 Hz, due to the neural firing rates [1] and the convolution of the firings and mechanical waveform [19]. As a consequence of this convolution, the PSD should have a high kurtosis and skewness. Next, we discriminate the mechanical units from noise using K -means clustering [24], [32] (K = 2) and selecting the class with highest centroid values ( Fig. 2d and 5e). Throughout the text, we define the discriminated subsets as twitch trains T and spatial maps S. VOLUME 8, 2020

4) FIRING PATTERN ESTIMATION
We standardized T to have zero mean and unit variance and apply a band-pass Butterworth filter (1 to 50 Hz) of order 6 before thresholding at 0.35, empirical choice based on performance on simulated data (not shown). Then, the firing pattern of the units in T was estimated by the time instants of local maxima of the twitch train ( Fig. 2f and 5g).

5) SPATIAL TERRITORY ESTIMATION
The spatial maps in S have continuous values, and binary maps were derived by thresholding the spatial map according to full width at half maximum (FWHM) procedure, corresponding to a threshold equal to 50% of the maximum pixel-value (Fig. 2e,f and Fig. 5f,g). To produce a connected spatial territory, we computed the convex hull of the threshold image. Thereafter, we calculate the corresponding diameter of the spatial territory as √ ((4 × Area/π)).

C. PERFORMANCE EVALUATION
The performance of the proposed method to identify the mechanical units was evaluated in simulated conditions using three descriptors: number of units, similarity between estimated and simulated firing pattern, and similarity between the estimated and simulated territories. Due to the ambiguity of permutations in ICA, we estimated the permutation matrix P of the spatial maps S and twitch trains T to be able to compare estimated and simulated units. The twitch train T may be autocorrelated by some lag because of wave propagation from the centroid pixel for each unit. To be able to compare the estimated and simulated twitch trains, we constructed a crosscorrelated adjusted twitch train T and minimize || T P − T || with respect to P and use the Hungarian Method [33] to estimate P. The number of units' metrics we used was divided into two different metrics. The first metric was the ratio between the number of estimated units and the simulated number of units. The second metric was the ratio between the number of estimated units (given that the estimated are the simulated one) and the simulated number of units. Firing pattern estimation performance was evaluated using the rate of agreement (RoA) between the estimated and the simulated firing patterns [29]. The RoA is calculated as where c j is the number of firings of the j th firing pattern that was identified (tolerance in firing timing ±30 ms), A j the number of false identified firings, and B j the number of unidentified firings. Territory estimation performance was evaluated using sensitivity and specificity of the estimated territory pixels compared the simulated component's pixels. To evaluate the method's performance, we simulated 100 realizations with no synchronization with 10 to 80 (increasing by 10) simultaneously active units at two noise levels (Gaussian white noise of 20 and 40 dB), which resulted in a total of 1600 simulations. We also evaluated the performance for different territory overlaps and different firing synchronizations.

1) SIMULATIONS
A detailed description of the simulations is presented in the Supplementary Information and here we give a short description. The simulation model is based on an EMG simulation model [34], replacing the electric action potential responses with mechanical spatiotemporal twitch responses. We modelled the mechanical response (muscle thickening) of a twitch during MU contraction in the plane perpendicular to fiber direction (cross-sectional). The MU territory is modelled as a circle and the single mechanical twitch response of a MU is modeled using in vivo empirical electro-stimulated MU tissue velocity waveform [19]. The model is based on the assumption of no mechanical connectivity between the fibers of different MUs and thus the force is transmitted along the fiber direction only. The parameter settings were set to mimic a biceps brachii muscle at weak isometric contraction levels. The firing pattern parameters of the MUs was set to a firing rate (FR) randomly between 8 and 13 Hz with an inter-pulse-interval variation modelled as a N (0, 0.2/FR) distribution [1], [35]. Synchronization was obtained by randomly shifting selected firings as proposed in [36] where the variability between synchronized firings was on average zero ms with two ms in standard deviation [36], [37]. The level of MU synchronization was simulated as the percentage of MU firings synchronized with firings of other MUs [36]. The sizes of MU territories were set to be randomly distributed between 2.5 and 10 mm in diameter and their centroids positioned randomly distributed in the simulated region. The simulated tissue velocity image sequences were 100 × 100 × 3000 pixels, corresponding to 40 mm × 40 mm × 3 s with spatial resolution of 0.4 mm/pixel, and 1 kHz frame rate (Fig. 2a).

D. EXPERIMENTAL TESTS
To test the feasibility of our method on experimental data, we compared characteristics of identified mechanical units with well-known characteristics of MUs from EMG studies, i.e., number of active units, firing rates, and territory sizes.
In particular, during weak voluntary contractions, only a few MUs should be activated [1], and the muscle produces higher force by successive recruitment of more MUs [38] at weak force levels [39], [40]. To test if we could observe this with our method, we compared mechanical units identified in experimental measurements at 1%, 2.5%, and 5% isometric maximal voluntary contraction (MVC).

1) SUBJECTS AND SETUP
Eight healthy subjects (27 to 45 years old, five men and three women) performed isometric elbow flexion with the right arm using a dynamometer (Kin-Com, 500H, Chattanooga Group, Inc., Tennessee, USA). The hip, back, and shoulders were strapped to the chair to ensure that the same body position was maintained during the exercise. The right elbow was positioned on elbow support on the arm support of the dynamometer. The arm support was adjusted so that the elbow joint and the hand were strapped to a plate in a supinated rotation [41]. To determine maximal force, the subjects carried out two MVCs, each lasting three seconds, with at least one min recovery period between each contraction. The highest force from the two trials was selected and used to calculate the force targets. The subjects received visual feedback of the force generation on a monitor, including target force (the straight solid line at respective force level) and generated force (horizontal dotted markers). The subject was instructed to generate a steady force as close as possible to the target without changing the position of the hand, elbow, or shoulder. The subjects were given two minutes to practice on 1 and 5% levels. Subsequently, the subjects performed an isometric elbow flexion at 1, 2.5, and 5% of MVC for two seconds. The subjects gave informed consent, and the project conformed to the Declaration of Helsinki and was approved by the regional ethical review board in Umeå (dnr 2012-300-31M).

2) ULTRASOUND ACQUISITION AND TISSUE VELOCITY ESTIMATION
A research ultrasound system (SonixTouch, Ultrasonix Medical Corporation, Richmond, CA) was used together with a 9 MHz L9-4 linear transducer and a 128 channel DAQ module, to acquire two seconds of data from a plane wave transmit and parallel receive sequence at an image rate of 2000 images/s [15]. The radiofrequency (RF) data was sampled at 40 MHz and reconstructed using sum-and-delay beamforming method [19]. Tissue velocity image sequences were estimated from the reconstructed RF images by using the two-dimensional autocorrelation approach [42] with a 1 mm maximal displacement in the depth direction for all VOLUME 8, 2020 channels between subsequent images and a sliding window of 10 ms. Only the axial component was retained [19]. The reconstructed RF-line signals were band-pass filtered (2-15 MHz with order 6) and a high-pass filter (5 Hz) was applied along the frame time of the velocity data, prior to spatial 2-D median filtering (1 × 1 mm kernel).

E. DATA PROCESSING
All data processing was carried out using MATLAB (2018b, Mathworks, Nattick, MA, USA), apart from modelling wave equations in COMSOL Multiphysics R (COMSOL Multiphysics R v. 5.2., COMSOL AB, Stockholm, Sweden) using the coefficient form PDE interface.

III. RESULTS
A. SIMULATIONS Fig. 2 demonstrates an example of a simulated tissue velocity sequence and how it is affected by the different modules of the method's processing pipeline. The tissue velocity sequence presents a complicated velocity field (Fig. 2a). After the SVD processing, activity from populations of units can be seen (Fig. 2b). The decomposition yield both components of likely mechanical units (Fig. 2c, top and bottom) as well as noise components (Fig. 2c, middle). The unit selection automatically discriminates the mechanical units from noise (Fig. 2d,e), and finally, their firing pattern and spatial territories (Fig. 2f) are extracted. The territory and firing pattern estimation performance evaluation is illustrated in Fig. 3a, and Fig. 3b respectively. The method found 75-95% of the simulated number of units when simulating up to 20 units, but beyond that, performance decreased (Fig. 4a). At low force levels, the method overestimated the number of units. Note: 20 units corresponds roughly to a contraction level of 5% MVC, assuming 200-300 units in the biceps muscle [43], [44], and full recruitment at around 50% MVC [45]. The estimated spatial territories had a sensitivity of 50-80% of the simulated pixel territories when simulating up to 50 units (Fig. 4b), and specificity was about 100%. The estimated firing patterns had an RoA of 90% to the simulated firing pattern when ten units were simulated (Fig. 4c), and performance decreased when including more units to 50% at 30 units. The higher noise level had a negative impact on spatial territory estimation but not on the other metrics.
We found that MU synchronization (0-20%) did not have any influence on the method's performance. In general, the MU territory overlap (0-100%) did not influence the performance, but there was a slight variation in sensitivity (see Supplementary Information).
In summary, the performance evaluation shows that our method has the potential to identify (overlapping and synchronized) individual units and their territories and firing patterns under isometric contractions at weak force levels up to 5% MVC. Fig. 5 shows an example of how the experimental tissue velocity sequence at 2.5% MVC is affected by the different modules of the method's processing pipeline. The acquired ultrasound data is from a cross-sectional image plane of the biceps brachii ( Fig. 1 and Fig. 5a), and the B-mode images showed no visible deformation or motion during the isometric contractions (Fig. 5b). Still, the tissue velocity showed a complicated velocity field in the whole muscle (Fig. 5b). After the SVD processing, reduction in complexity can be seen (Fig. 5c), and the decomposition yield both components of likely mechanical units (Fig. 5d, top and bottom) as well as likely noisy components (Fig. 5d, middle). The unit selection discriminated the probable mechanical units (Fig. 5e,f). Finally, their firing pattern and spatial territories were extracted, as well as the twitch response of a unit estimated by spike-triggered averaging (Fig. 5g). The effect of the different force levels 1-5% MVC is exemplified in Fig. 6, and the results from all eight subjects and force levels are summarized in Fig. 7. All force levels showed subtle tissue velocity with a complicated pattern in the whole muscle without any visible movement in the B-mode sequences (Supplementary Information). We found that the number of units increased with force level (Fig. 6), where the number of units at 1%, 2.5%, and 5% MVC averaged 7, 9, and 12, respectively (Fig. 7a). We found that the territory diameter (5-6 mm) and contraction duration of the twitch responses (40-50 ms) were similar for all force levels (Fig. 7c-d). We also found that the firing rate (11)(12) was similar for all force levels (Fig. 7b). The units were found at depths from a few mm down to 40 mm without any obvious relation between force level and depth within the muscle (Fig. 6). The results were consistent in all eight subjects (Supplementary Fig. 1 to 8).

B. EXPERIMENTAL TESTS
In summary, the results show that the identified mechanical units increased in number with increasing force while retaining similar territory sizes, firing rates, amplitude, and duration of the twitch response without any relation to the depth within the muscle.

IV. DISCUSSION
The smallest unit of voluntary contraction in skeletal muscle is the MU, and its function has been predominately studied using EMG methods. In this study, we introduce a complementary method to study so-called mechanical units, i.e., the mechanical response of an electrically activated MUs, during voluntary contraction using non-invasive ultrafast ultrasound imaging and spatiotemporal decomposition. We show that the method can identify simulated MUs' spatial territories and firing patterns under low force isometric contractions. Also, we tested the feasibility of the method in experimental measurements from the biceps brachii muscle. We found that the characteristics of the identified mechanical units were similar to those of MUs described by EMG [1], [7].

A. METHODOLOGICAL CONSIDERATIONS AND SIMULATION PERFORMANCE
The input to the method is a tissue velocity sequence measuring the axial component in a cross-sectional image plane to the fiber orientation. During isometric voluntary contraction, the muscle is in a steady state of repeated shortening and thickening of the fibers. Thus, no global motion of the tissue can be assumed, and the mechanical response of the activated MU can be seen as radial vibrations. Therefore, axial or lateral velocity components should contain similar information. Since the axial component allows superior velocity resolution estimation (approximately 20 micrometer displacements at 10 MHz) [15], this was chosen. If the orientation of the fibers VOLUME 8, 2020 of the MU is not perfectly perpendicular to the skin's surface, the velocity component and the estimated twitch response will be a superposition of both thickening and shortening motion components of the fibers. While this could limit the interpretation of detailed mechanical parameters of the twitch such as duration and amplitude, it should not have any critical influence on firing pattern or MU territory estimation, as long as an oscillatory motion of the twitch trains can be identified and the projection is approximately perpendicular to the fibers.
The application of dimension reduction by SVD is a standard procedure in ICA pre-processing to reduce computational complexity and achieve whitening [31]. The data dimension was reduced to k = 100 components, which resulted in approx. 99% of the variation in the data is conserved. This implies that a maximal limit of 100 components can be decomposed, but, in this study, we consider low force contractions where only some tens of MUs should be active [45]. The data decomposition was based on a spatiotemporal approach (stICA) using sparsity as separation cost function (kurtosis). The cost function includes a weight in likelihood function related to temporal and spatial sparsity of the decomposed units. Although the temporal twitch train signals have non-Gaussian distributions, we used emphasis on spatial sparsity (α = 0.8). This is motivated by the fact that the spatial territory of an MU is located within a small localized part of the tissue. Thus, an MU's mechanical contraction should mainly take place in that region and nowhere else in the image, causing the spatial maps to have a sparse spatial distribution of velocities. Moreover, the method had the best performance with this α on the simulated signals out of the tested α ranging from 0.0 to 1.0 with increments of 0.2 (data not shown).
The selection of the most probable MU twitch trains from the noise signals was based on K-means clustering. Two assumptions were applied: First, K = 2 implies that there is always both a cluster with probable signals and noisy signals, even if there would be no active MUs. Second, the features describing the temporal signals of the decomposed units of interest are expected to have a dominant peak in the frequency domain and that it should be skewed with higher density towards the lower frequencies. Since the unit selection is based on temporal characteristics as compared to the spatial emphasis in the separation step, they complement each other and ascertains the identification of mechanical units based on a spatiotemporal framework. The unit selection approach is similar to [28] but could likely be improved, e.g., by using more or optimized features and also to use a non-parametric clustering approach.
The territory and firing patterns were estimated using thresholding. The spatial territory estimation is affected such that a lower threshold includes a larger territory and thus has a direct impact on the sensitivity and specificity. According to the simulation results (Fig. 4), territories were on average underestimated, which might be possible to correct by using a lower threshold (see Spatial territory estimation in the Methods section). The firing pattern threshold at 0.35 was empirically set, and it may influence the number of detected peaks. Simulation results on firing pattern showed a 90% RoA for ten simulated units, which means that 90% agreement of true and estimated firings. The RoA decreased by an increasing number of units, with 75% for 20 units and 50% for 30 units. Both thresholding parameters could have been optimized using the simulated data. However, this work is aimed at providing proof-of-concept and not to optimize performance. A potential solution in the future is to optimize such thresholds with respect to territory and firings measurements using simultaneous scanning EMG.
In general, individual firings of different MUs can be synchronized [46], and this means that about 2-4 out of 20 firings could coincide over two-second contraction with 10 Hz firing rate (10-20% synchronization). Although the overall performance evaluation (Fig. 4) did not take this into account, our results evaluating the effect of synchronization (0-20%) showed that it did not have any substantial impact on the method's performance (Fig. 10, Supplementary Information). We believe it is due to the prioritized weight on spatial sparsity. In addition, our results evaluating the effect of MU territory overlap (0-100%) showed that it did not have any substantial impact on the method's performance (Fig. 10, Supplementary Information). This result may be explained by the simultaneous temporally and spatially constraints in the decomposition procedure.
The performance evaluation results show that up to 30 units can be identified (see first paragraph and Fig. 4). This corresponds to a contraction level of about 5% MVC, if translated to a biceps brachii comprising about 200-300 units [43], [44], and assuming full MU recruitment at about 50% MVC [45]. Thus, the method has, in its present form, the potential to be used on low-level contraction forces.

B. SIMULATION MODEL
The simulation model is based on an EMG simulation model [34], where a mechanical twitch response replaces the electrical response. The rationale of the response replacement is the so-called excitation-contraction coupling, which acts as a direct link between the action potential excitation of the fiber membrane and the subsequent mechanical contraction of the sarcomeres micro-filaments causing the contraction twitch. This replacement is based on two key assumptions. First, the tissue velocity sequence is considered a linear combination of the velocities of the individual contracting MUs. Although there are indications of non-linear coupling at higher force and high firing rates (>20 Hz) based on electro-stimulation and mechanomyography (MMG) studies [19], [47], this work deals with low force contractions with low firing rates and few active MUs. At low force contractions, MMG studies [48], [49] support that the tissue velocity field can be approximated by linear superposition. Second, extracellular matrix, endo-, perimysium that is part of transferring the force to the tendons are not included, which implies that only longitudinal force transmission is assumed. No interaction and cross-talk in motions from adjacent active MUs and lateral force transmission components are assumed. This assumption may have resulted in higher performance on simulated data than a simulation model that includes fascia. In future work on developing methods to identify mechanical units in high force contractions, a state-of-the-art full biomechanical model should be used, e.g., as in [50].

C. EXPERIMENTAL TESTS
The number of detected units at the force levels 1%, 2.5%, and 5% MVC were on average 6, 9, 12 (Fig. 7). The number of units increased with force level, which is in line with previous research on MU recruitment [38]. Assuming that the biceps have 200-300 MUs [43], [44] and that all units are recruited at 50% MVC [45], we roughly estimate that about 4-6, 10-15, and 20-30 units may be active in the whole muscle at the force levels in this study. The territories of the identified units were about 6 mm in diameter for all force levels and are similar to previously reported findings of MU territory size in ultrasound imaging electro-stimulation studies [51], [52], and scanning-EMG studies [7]. The firing patterns had about 8-13 Hz firing rate and were similar in all force levels. The performance evaluation on simulation results with an RoA of 50-90% indicated that the method has a high agreement in the estimation of firing patterns for a few units. Moreover, our findings are similar to motor neuron firing rates at weak force levels [1]. Although it is known that VOLUME 8, 2020 firing rate, in addition to MU recruitment (number of units), contributes to force production at high force levels [53], [54], there is no such consensus at the low force levels studied here [39], [40]. The estimated duration of the twitches was about 50 ms in all force levels and is similar to the contraction phase in previous reports using electro-stimulated muscle activation at weak contraction levels [19], [20]. In summary, these results were similar to the known characteristics of MUs [19], [20], [55]- [57].
From literature, it is known that MU territories generally overlap, but in our experimental tests, we see very few cases with overlap (see Supplementary Information). This result could be a failure of our method to detect such overlapping MUs in real data even though the simulations showed it had minimal effect on performance. However, considering the very low contraction level in our data, few MUs should be active, and the probability of MU overlap should thus be small when the active units are distributed within the whole muscle as they were in our data.
MUs are expected to be successively recruited with increasing force. To assess whether individual mechanical units could be identified at increasing force levels, we colorcoded and manually labeled individual units that were present at the different force levels (Fig. 11, Supplementary Information). We found that 26% of the active units at the first level were present in all subsequent force levels. In addition, 44% of the active units at 5% MVC were active at 1 or 2.5% MVC ( Table 2, Supplementary Information). Thus, a unit active at first level was not necessarily active at the consecutive levels. This may be caused by a failure of the method to detect all active units. However, given these weak force levels and that each force level was recorded in a separate experiment (with two minutes rest between each), these findings could be explained by a simultaneous agonist-antagonist activation strategy used by the motor system (co-contraction) to facilitate movement accuracy and produce weak forces [58].
A key difference between simulated and experimental data is the presence of fascia, i.e., there is a possibility that it may cause artifacts since vibrations and motions will be coupled through the fascia in the muscle. The fascia transfers the force produced by the fibers of an MU to the tendons via both longitudinal and lateral force transmission. The latter implies potential cross-talk between adjacent MUs that could influence our decomposition [59]. However, we did not see any relation between the position of the detected territories and the regions with strongly reflecting fascia in the B-mode images. The different processing steps of the method could also potentially have generated these mechanical units. However, the ICA does not include any spatial dependency (on nearby pixels) but, still, the units that were extracted were spatially confined in a regionally local territory, which indicates that they are mechanical units generated by an MU. In addition, the estimated twitch responses based on averaging of the firing pattern had high SNR, and this indicates that the units' detected firings were derived from an oscillatory physiological source. Taken together, the results of the experimental tests are promising but do not provide explicit evidence that the method can detect the MUs. A validation with gold standard intramuscular EMG is required. Table 1 summarizes the comparison between the state-ofthe-art EMG techniques and our ultrasound-based method. The EMG techniques are fundamentally different from our method, i.e., EMG samples electrical activity and our method samples mechanical activity. The added value of the electrical activity is the MUs depolarization features, whereas the mechanical activity's added value is MUs contractility features. However, some information may (in theory) be the same, e.g., MU territory and firing pattern. The spatial resolution is much higher in ultrasound than in EMG, and using a high-frequency probe can sample a spatial resolution down to 50 micrometers.

E. APPLICATIONS
The ability to non-invasively identify simultaneously active MUs during a contraction would allow new research opportunities on neuromuscular physiology, i.e., strategies of the central nervous system on MU recruitment. This has been identified as an important function where, e.g., changes in the activity of different intramuscular regions has been shown to prevent fatigue in healthy subjects [41], and also to be related to pain intensity in patients with fibromyalgia [62]. In clinical applications, the ability to identify MUs in the whole muscle would allow larger accessibility than current diagnostic methods. The assessment of spatial territory and the number of active units could allow opportunities to diagnose reinnervation diseases such as Poliomyelitis, in which some of the MU fibers get innervated by other MUs [63]. Combining mechanical unit twitch information with the corresponding electrical excitation information of EMG methods could open up new ways to diagnose myopathies and triadopathies [64], i.e., a critical illness where the mechanical coupling is defect due to myosinopathy [65]. Altogether, our method could allow the study of a variety of questions that previously were difficult or not possible to address.

V. CONCLUSION
Electrophysiological methods, including needle and surface EMG, are well-established techniques that have been used for over 60 years to help our understanding of muscle physiology and diseases. Our present study presents a new and non-invasive method to identify mechanical units, i.e., the mechanical response of an activated MU. Simulations demonstrate proof-of-concept at low contraction forces, and tests on experimental data show that individual mechanical units can be identified with characteristics similar to MUs described by EMG. Our method complements the EMG methods by adding 2-D spatial position and unit territory size of simultaneously active units in the whole muscle. This added value could allow the study of a variety of questions related to MUs and skeletal muscle tissue that previously were difficult or not possible. ERIK STÅLBERG made his thesis on propagation velocity in single muscle fibers, in 1966. He became a Full Professor in clinical neurophysioloy at Uppsala University, in 1991. He has studied microphysiology in normal and diseased muscle and has developed a number of EMG methods such as SFEMG, Macro EMG, Scanning EMG, and methods for quantitation of conventional EMG; these methods are used to study the motor unit from different aspects.

KAREN-HELENE STÖVERUD
received the master's degree in hydrogeology from Utrecht University, The Netherlands, in 2009, and the Ph.D. degree from the Center for Biomedical Computing, Simula Research Lab, University of Oslo, Oslo, Norway, in 2014. She is currently a Postdoctoral Researcher at the Department of Radiation Sciences, Umeå University, Sweden. Her current research interests include mathematical modeling of cerebrospinal fluid circulation and patient specific computational hemodynamics.
JUN YU received the B.S. degree in mathematics from Hangzhou University, China, in 1982, and the Ph.D. degree in mathematical statistics from Umeå University, Sweden, in 1994. He has also been working as a Professor with the Centre of Biostochastics, Swedish University of Agricultural Sciences, before 2012. Since 2012, he has been the Chaired Professor in mathematical statistics with Umeå University. His current research focus is on statistical learning with sparsity, nonconvex learning, compressive sensing, hierarchical spatiotemporal modelling, spatiotemporal marked point processes, and wavelet theory for sparse signal processing. In terms of data analysis tools, it includes intelligent data sampling using compressive sensing, multimodal biomedical image processing, spatiotemporal tree growth models, statistical classification, and pattern recognition of painted vehicle bodies. He has been working with many research projects with real applications, where ongoing projects include biomedical engineering, cancer therapy assessment, forestry, icing impact on train transportation, and industrial quality control.
CHRISTER GRÖNLUND received the Ph.D. degree from Umea University, Sweden, in 2007. Since 2012, he has been an Associate Professor in biomedical engineering with the Department of Radiation Sciences, Umea University. He is the principal investigator of a research group working on method development for electrophysiological, biooptical, and ultrasound imaging. The medical areas of application include vascular assessment for cardiovascular disease risk prediction, and neuromuscular diseases and disorders for diagnostics and sports medicine. VOLUME 8, 2020