Larger and Denser: An Optimal Design for Surface Grids of EMG Electrodes to Identify Greater and More Representative Samples of Motor Units

Abstract The spinal motor neurons are the only neural cells whose individual activity can be noninvasively identified. This is usually done using grids of surface electromyographic (EMG) electrodes and source separation algorithms; an approach called EMG decomposition. In this study, we combined computational and experimental analyses to assess how the design parameters of grids of electrodes influence the number and the properties of the identified motor units. We first computed the percentage of motor units that could be theoretically discriminated within a pool of 200 simulated motor units when decomposing EMG signals recorded with grids of various sizes and interelectrode distances (IEDs). Increasing the density, the number of electrodes, and the size of the grids, increased the number of motor units that our decomposition algorithm could theoretically discriminate, i.e., up to 83.5% of the simulated pool (range across conditions: 30.5–83.5%). We then identified motor units from experimental EMG signals recorded in six participants with grids of various sizes (range: 2–36 cm2) and IED (range: 4–16 mm). The configuration with the largest number of electrodes and the shortest IED maximized the number of identified motor units (56 ± 14; range: 39–79) and the percentage of early recruited motor units within these samples (29 ± 14%). Finally, the number of identified motor units further increased with a prototyped grid of 256 electrodes and an IED of 2 mm. Taken together, our results showed that larger and denser surface grids of electrodes allow to identify a more representative pool of motor units than currently reported in experimental studies.


Significance Statement
The application of source separation methods to multichannel electromyographic (EMG) signals recorded with grids of electrodes enables users to accurately identify the activity of individual motor units.However, the design parameters of these grids have never been discussed.They are usually arbitrarily fixed, often based on commercial availability.Here, we showed that using larger and denser grids of electrodes than conventionally proposed drastically increases the number of identified motor units.The samples of identified units are more balanced between early-recruited and late-recruited motor units.Thus, these grids provide a more representative sampling of the active motor unit population.Gathering large datasets of motor units using large and dense grids will impact the study of motor control, neuromuscular modeling, and human-machine interfacing.

Introduction
Decoding the neural control of natural behaviors relies on the identification of the discharge activity of individual neural cells.Classically, arrays of electrodes are implanted close to the cells to record their electrical activity.The application of algorithms that separate the simultaneous and overlapping activity of these cells has enabled researchers to study neural processes in multiple areas of the brain (Stringer et al., 2019), such as in the motor or the sensorimotor areas (Churchland and Shenoy, 2007;Gallego et al., 2020).At the periphery of the nervous system, it is also possible to record the activity of individual motor neurons innervating muscle fibers (Duchateau and Enoka, 2011;Heckman and Enoka, 2012;Farina et al., 2016).The motor unit, i.e., a motor neuron and the muscle fibers it innervates, acts as an amplifier of the neural activity, as one action potential propagating along a motor neuron's axon generates an action potential in each of the innervated muscle fibers.The discharge activity of motor units can be identified by decomposing surface electromyographic (EMG) signals into trains of motor unit action potentials (MUAPs) using, e.g., blind-source separation algorithms (Holobar and Farina, 2014;Farina and Holobar, 2016).The multiple observations for source separation are obtained by recording EMG signals with grids of electrodes.This approach usually allows for the reliable analysis of 5-40 concurrently active motor units (Del Vecchio et al., 2017, 2020;Hug et al., 2021a).
While the design of intracortical (Jun et al., 2017;Steinmetz et al., 2018) and intramuscular (Muceli et al., 2015(Muceli et al., , 2022) ) electrodes arrays has scaled up over the years to record larger samples of neural cells, the configuration of grids of surface EMG electrodes has not systematically evolved.Most researchers currently use grids with 64 electrodes arranged in 13 Â 5 or 8 Â 8 montages, the interelectrode distance (IED) between adjacent electrodes (e.g., 4, 8, or 10 mm) being dictated by the size of the muscle to cover.Yet, optimizing these parameters, i.e., grid size and IED, may influence the performance of EMG decomposition.Currently, there are no recommendations on optimal designs for grids of electrodes.
Source separation algorithms are based on the necessary condition that identifiable motor units have a unique representation of their action potentials across the multichannel EMG signals (Farina et al., 2008;Holobar and Farina, 2014;Farina and Holobar, 2016).This implies that the three-dimensional waveform of a MUAP (one time dimension and two spatial dimensions) is unique within the pool of active motor units detected by the grid of electrodes.In practice, the identified motor units are those that innervate larger numbers of muscle fibers, as their action potentials tend to have the largest energy.Conversely, low-threshold motor units usually remain hidden since their energy is close to the baseline noise.Increasing the density of electrodes would increase the spatial sampling of EMG signals (Farina and Holobar, 2016), which in turn should improve the discrimination of MUAPs, allowing the identification of a higher number of motor units.Additionally, increasing the density of electrodes may reveal the hidden low-threshold motor units by sampling their action potentials across a higher number of electrodes, leading to a better compensation of the additive noise in the mixture model of the EMG signal (Farina and Holobar, 2016).
In this study, we combined computational and laboratory experiments to identify the optimal design parameters of grids of surface electrodes with the aim to maximize the number of identified motor units.We first simulated a pool of 200 motor units and the associated EMG signals recorded from grids of electrodes of various sizes and densities.These simulations showed that the greater the size and the density of the grid, the higher the percentage of theoretically identifiable motor units and the relative ratio of theoretically identifiable deep units.We confirmed these theoretical results with experimental signals recorded with a grid of 256 electrodes with a 4-mm IED that was downsampled in the space domain to obtain six grid configurations (surface range: 2-36 cm 2 and IED range: 4-16 mm).Finally, we prototyped a new grid of 256 electrodes with a 2-mm IED and demonstrated that the number of identified motor units further increased with 2-mm IED.The entire dataset (raw and processed data) and codes are available at https://doi.org/10.6084/m9.figshare.22149287.

Computational study
A pool of 200 motor units was simulated to test whether increasing the density and the size of surface grids of electrodes would impact the number of theoretically identifiable motor units.The simulations were based on an anatomic model entailing a cylindrical muscle volume with parallel fibers (Farina et al., 2008;Konstantin et al., 2020), in which subcutaneous and skin layers separate the muscle from the surface electrodes.Specifically, we set the radius of the muscle to 25.4 mm and the thicknesses of the subcutaneous and skin layers to 5 and 1 mm, respectively.The centers of the motor units were distributed within the cross-section of the muscle using a farthest point sampling technique.The farthest point sampling filled the cross-section by iteratively adding centers points that were maximally distant from all the previously generated motor unit centers, resulting in a random and even distribution of the motor unit territories within the muscle.The number of fibers innervated by each motor neuron followed an exponential distribution, ranging from 15 to 1500.The fibers of the same motor unit were positioned around the center of the motor unit within a radius of 0.2-9.8mm, and a density of 20 fibers/mm 2 .Because motor unit territories were intermingled, the density of fibers in the muscle reached 200 fibers/mm 2 .The MUAPs were detected by circular surface electrodes with a diameter of 1 mm.The simulated grids were centered over the muscle in the transverse direction, with a size ranging from 14.4 to 36 cm 2 and an IED ranging from 2 to 36 mm.

Participants
Six healthy participants (all males; age: 26 6 4 years; height: 174 6 7 cm; body weight: 66 6 15 kg) volunteered to participate in the first experimental session of the study.They had no history of lower limb injury or pain during the months preceding the experiments.One of these individuals (age: 26 years; height: 168 cm; bodyweight: 51 kg) participated in a second experimental session to test the prototyped grid with an IED of 2 mm.The Ethics Committee at Imperial College London reviewed and approved all procedures and protocols (no.18IC4685).All participants provided their written informed consent before the beginning of the experiment.

Experimental tasks
The two experimental sessions consisted of a series of isometric ankle dorsiflexions performed at 30% and 50% of the maximal voluntary torque (MVC) during which we recorded high-density electromyographic (HD-EMG) signals over the tibialis anterior muscle (TA).The participant sat on a massage table with the hips flexed at 30°, 0°b eing the hip neutral position, and their knees fully extended.We fixed the foot of the dominant leg (right in all participants) onto the pedal of a commercial dynamometer (OT Bioelettronica) positioned at 30°in the plantarflexion direction, 0°being the foot perpendicular to the shank.The thigh was fixed to the massage table with an inextensible 3-cm-wide Velcro strap.The foot was fixed to the pedal with inextensible straps positioned around the proximal phalanx, metatarsal and cuneiform.Force signals were recorded with a load cell (CCT Transducer s.a.s.) connected in-series to the pedal using the same acquisition system as for the HD-EMG recordings (EMG-Quattrocento; OT Bioelettronica).The dynamometer was positioned accordingly to the participant's lower limb length and secured to the massage table to avoid any motion during the contractions.
All experiments began with a warm-up, consisting of brief and sustained ankle dorsiflexion performed at 50% to 80% of the participant's subjective MVC.During the warmup, all participants learnt to produce isometric ankle dorsiflexion without co-contracting the other muscles crossing the hip and knee joints.At the same time, we iteratively adjusted the tightening and the position of the straps to maximize the comfort of the participant.Then, each participant performed two 3-to 5-s MVC with 120 s of rest in between.The peak force value was calculated using a 250-ms moving average window, and then used to set the target level during the submaximal contractions.After 120 s of rest, each participant performed two trapezoidal contractions at 30% and 50% MVC with 120 s of rest in between, consisting of linear ramps up and down performed at 5%/s and a plateau maintained for 20 and 15 s at 30% and 50% MVC, respectively.The order of the contractions was randomized.One participant (S2) did not perform the contractions at 50% MVC.

High-density electromyography
In the first experimental session, four adhesive grids of 64 electrodes (13 Â 5 with a missing electrode in a corner; gold coated; 1-mm diameter; 4-mm IED; OT Bioelettronica) were placed over the belly of the TA.The grids were carefully positioned side-to-side with a 4-mm distance between the electrodes at the edges of adjacent grids (Fig. 1A).The 256 electrodes were centered to the muscle belly and laid within the muscle perimeter identified through palpation.The skin was shaved, abrased and cleansed with 70% ethyl alcohol.Electrode-to-skin contact was maintained with a bi-adhesive perforated foam filled with conductive paste.The grids were wrapped with tape and elastic bands to secure the contact with the skin.The four preamplifiers were connected in-series with stackable cables to a wet reference band placed above the medial malleolus of the same leg.HD-EMG signals were recorded in monopolar derivation with a sampling frequency of 2048 Hz, amplified (150Â), bandpass filtered (10-500 Hz), and digitized using a 400 channels acquisition system with a 16-bit resolution (EMG-Quattrocento; OT Bioelettronica).
In the second experimental session, one ultradense prototyped grid of 256 electrodes (26 Â 10 with a missing electrode in each corner; gold coated; 1-mm diameter; 9 cm 2 area; 2-mm IED; custom-manufactured for this study by OT Bioelettronica; Fig. 1H) was placed over the belly of the TA and the HD-EMG signals were recorded using the same procedure as previously described.

Grid configurations
During the first experimental session, we recorded EMG signals from the TA with a total of 256 electrodes covering an area of 36 cm 2 over the muscle (10 Â 3.6 cm, 4-mm IED; Fig. 1A).To investigate the effect of electrode density, we down-sampled the grid of 256 electrodes by successively discarding rows and columns of electrodes and artificially generating three new grids covering the same area with IEDs of 8, 12, and 16 mm, that involved 256, 64, 35, and 20 electrodes, respectively (Fig. 1B-D).It is noteworthy that the 8-and 16-mm grids covered a surface of 32 cm 2 because they included an odd number of rows and columns.To investigate the effect of the size of the grid, we discarded the peripherical electrodes to generate grids of 63, 34, and 19 electrodes with a 4-mm IED, covering areas of 7.7, 3.8, and 2 cm 2 over the muscle (Fig. 1E-G).We chose these sizes to match the number of electrodes used in the density analysis, thus comparing grids with similar number of electrodes, but different densities and sizes (in Fig. 1B vs E and C vs F).
During the second experimental session, we recorded EMG signals from the TA with an ultradense grid of 256 electrodes covering an area of 9 cm 2 over the muscle (5 Â 1.8 cm, 2-mm IED; Fig. 1H).Using the same procedure as above, we generated two artificial grids of 64 and 32 electrodes with an IED of 4 and 8 mm, respectively.

HD-EMG decomposition
We decomposed the signals recorded in all the conditions using the same algorithm, parameters, and procedure.First, the monopolar EMG signals were bandpass filtered between 20 and 500 Hz with a second-order Butterworth filter.The channels with low signal-to-noise ratio (PNR) or artifacts were discarded after visual inspection.The HD-EMG signals were then decomposed into individual motor unit pulse trains using convolutive blindsource separation, as previously described (Negro et al., 2016).In short, the EMG signals were first extended by adding delayed versions of each channel.We kept the same extension factor for all the conditions to reach 1000 extended channels, as previously suggested (Negro et al., 2016).The extended signals were spatially whitened to make them uncorrelated and of equal power.Thereafter, a fixed-point algorithm was applied to identify the sources embedded in the EMG signals, i.e., the motor unit pulse trains, or series of d functions centered at the motor unit discharge times.In this algorithm, the contrast function g(10) ¼ log(cosh( 10)) was iteratively applied to the EMG signals to skew the distribution of the values of the motor unit pulse trains toward 0, and thus maximize the level of sparsity of the motor unit pulse train.The high level of sparsity matches the physiological properties of motor units, with a relatively small number of discharges per second (,50 discharge times/s during submaximal isometric contractions).The convergence was reached once the level of sparsity did not substantially vary (with a tolerance fixed at 10 À4 ) when compared with the previous iteration (Negro et al., 2016).At this stage, the motor unit pulse train contained high peaks (i.e., the d functions from the identified motor unit) and lower values because of the activities of other motor units and noise.High peaks were separated from lower values using peak detection and Kmean classification with two classes.The peaks from the class with the highest centroid were considered as the discharge times of the identified motor unit.A second algorithm refined the estimation of the discharge times by iteratively recalculating the motor unit filter and repeating the steps with peak detection and K-mean classification until the coefficient of variation of the interspike intervals was minimized.This decomposition procedure has been previously validated using experimental and simulated signals (Negro et al., 2016).After the automatic identification of the motor units, duplicates were automatically removed.For this purpose, the pulse trains identified from pairs of motor units were first aligned using a cross-correlation function to account for a potential delay because of the propagation time of action potentials along the fibers.Then, two discharge times were considered as common when they occurred within a time interval of 0.5 ms, and two or more motor units were considered as duplicates when they had at least 30% of their identified discharge times in common (Holobar et al., 2010).In principle, the limited level of synchronization between individual motor units results in a few simultaneous discharges between pairs of motor units.A threshold of 30% is therefore highly conservative to ensure the removal of all motor units with a level of synchronization well above physiological values.It is worth noting that most of the motor units identified as duplicates after the automatic decomposition had almost 100% of their discharge times in common.In that case, the motor unit with the lowest coefficient of variation of the interspike intervals was retained for the analyses.At the end of these automatic steps, all the motor unit pulse trains, i.e., the output of the decomposition resulting from the projection of EMG signals onto individual motor unit filters, were visually inspected, and manual editing was performed to correct the false identification of artifacts or the missed discharge times (Del Vecchio et al., 2020;Hug et al., 2021b;Avrillon et al., 2023).The update of the motor unit filters with the corrected discharge times and the recalculation of the motor unit pulse trains always improved the distance between the discharge times and the noise, quantified with the pulse-to-noise ratio (PNR; Holobar et al., 2014).Note that this manual step is highly reliable across operators, as previously demonstrated by Hug et al. (2021b).Duplicates were checked a second time after manual editing, with very rare cases of removal as most of the duplicates were automatically identified after the automatic decomposition.Only the motor unit pulse trains which exhibited a PNR .28 dB after manual editing were retained for further analysis.
We further tested whether decomposing subsets of electrodes within a highly populated grid of 256 electrodes increased the number of identified motor units.Indeed, the lower ratio of large motor units sampled by each independent subset of 64 electrodes could allow the algorithm to converge to smaller motor units that contribute to the signal.For a similar number of iterations, it is likely that these motor units would have otherwise contributed to the noise component of the mixture model of the EMG signal (Farina and Holobar, 2016).Thus, we decomposed the grids of 256 electrodes (4and 2-mm IED; Fig. 1A,H) as four separated grids of 64 electrodes before removing the motor units duplicated between grids.

Analyses Computational study
We first estimated the theoretical percentage of identifiable motor units for each of the simulated conditions.To do so, the simulated MUAPs detected over the entire set of electrodes were compared with each other.The comparisons were done pairwise by first aligning the MUAPs in time using the cross-correlation function, and then computing the normalized mean square difference between the aligned action potentials.Pairs of action potentials with a mean square difference below 5% were considered not discriminable.The 5% criterion was based on the variability of motor unit action potential shapes observed experimentally for individual motor units (Farina et al., 2008).After computing all pair-wise comparisons, we then computed the percentage of action potentials that could be discriminated from all others, i.e., the theoretical percentage of identifiable motor units.This metric is independent from the algorithm used for decomposition and establishes a theoretical upper bound in the number of motor units that can be identified by any decomposition algorithm.For each theoretically identifiable motor unit, we also computed the distance between the center of the territory of the corresponding muscle fibers and the skin surface.

Laboratory study, number of identified motor units
We reported the absolute number of motor units (PNR .28 dB) identified with all the grid configurations.For each participant, the number of identified motor units was then normalized to the maximal number of motor units found across all conditions, yielding normalized numbers of identified motor units N expressed in percentage.For each condition, we calculated the mean and standard deviation of the N values across participants.To investigate the effects of density and size of the grid, we fitted logarithmic trendlines to the relationships between the averaged N values and IED or grid size.We also fitted a logarithmic trendline to the average N values and their corresponding number of electrodes, in which case the conditions involving the same number of electrodes, but different grid size and density, were given a weight of 0.5 in the minimization function.We reported the r 2 and p-value for each regression trendline.To maintain consistency with the computational study investigating the number of theoretically identifiable motor units across grid designs, the trendlines were fitted on the results obtained when the complete grids of 256 electrodes were decomposed as independent subsets of 64 electrodes, which systematically returned the highest number of identified motor units.The trendlines fitted on the results obtained with the decomposition of the 256 electrodes as a whole are reported in Extended Data Figure 4-1.

Laboratory study, properties of identified motor units
To investigate the effects of electrode density and grid size on the properties of the identified motor unit, we used a typical frequency distribution of the motor unit force recruitment thresholds in the human TA (Caillet et al., 2022b), where F th j ð Þ is the force recruitment threshold of the j th motor unit in the normalized motor unit pool ranked in ascending order of F th : The identified motor units were then classified according to this relationship and their measured force recruitment threshold, the first half of the active pool being "early recruited," and the second half "late recruited" (Henneman and Mendell, 1981;Caillet et al., 2022a).For each condition, we reported the percentage of identified motor units that were "early recruited."We did not report this metric when five or fewer motor units were identified in one condition for three or more participants.

Laboratory study, correlation between observations
We assessed how the density of electrodes impacted the information redundancy in EMG signals recorded by adjacent electrodes.To this end, MUAP shapes were identified over the 256 electrodes with the spike-triggered averaging technique.To do so, the discharge times were used as a trigger to segment and average the HD-EMG signals over a window of 50 ms.For each motor unit, we identified the electrode with the highest action potential peak-to-peak amplitude and calculated the average correlation coefficient r between this action potential and those recorded by the four adjacent electrodes with an IED of 4, 8, 12, and 16 mm.We also repeated this correlation analysis for the ultradense grid of 256 electrodes using an IED of 2, 4, and 8 mm.

Results
All the datasets (raw and processed data) and codes used to process the data are available at https://doi.org/10.6084/m9.figshare.22149287.

Computational study
We simulated the discharge activity of 200 motor units recorded by 84 configurations of grids of electrodes (surface range: 14.4-36 cm 2 , IED range: 2-36 mm; Fig. 2).The number of theoretically identifiable motor units increased with the size of the grid, from 46.7 6 7.7% of the motor units theoretically identifiable with a grid of 14.4 cm 2 to 77.8 6 5.5% of the motor units theoretically identifiable with a grid of 36 cm 2 .The number of theoretically identifiable motor units also increased with shorter interelectrode distances.For example, with a grid of 36 cm 2 , the number of theoretically identifiable motor units increased from 63.5% to 83.5% of the motor units with an IED of 36 and 2 mm, respectively (Fig. 2B).Increasing the surface size and the density of the grid of electrodes revealed deeper motor units.The averaged distance of theoretically identifiable motor units from the skin increased with the size of the grid (14.3 6 0.1 vs 16.5 6 0.2 mm with grids of 14.4 and 36 mm 2 , respectively; Fig. 2C), but not with the IED of the grid (15.6 6 1.1 vs 15.5 6 0.9 mm with an IED of 36 and 2 mm, respectively; Fig. 2D).

Laboratory study, grids of 256 electrodes with an IED of 4 mm Number of identified motor units
The motor unit pulse trains automatically identified across all conditions, intensities, and participants were visually inspected and carefully edited when a missing discharge time or a falsely identified artifact were observed.On average, 9 6 4% and 22 6 9% of the motor Figure 2. Results from the 200 simulated motor units (MUs) with 84 configurations of grids of electrodes.A, Each solid line represents a motor unit territory, the scatters being the muscle fibers.Blues lines are the theoretically identifiable motor units with a grid of 21.6 cm 2 and an interelectrode distance (IED) of 18 mm, while the orange lines are the motor units revealed with a grid of 21.6 cm 2 and an IED of 2 mm.Gray lines represent the nonidentifiable motor units.The percentage of theoretically identifiable motor units (B) and their distance from the skin (C) are reported for the 84 configurations.
units automatically identified at 30% and 50% MVC, respectively, were removed after visual inspection and manual editing.Furthermore, when the four grids of 64 electrodes were separately decomposed, 30 6 5% and 24 6 6% of the automatically identified motor units were removed because they were identified in more than one grid (only one pulse train was retained in case of duplicates).The highest number of identified motor units was systematically reached with the separate decomposition of the four grids of 64 electrodes with an IED of 4 mm, with 566 14 motor units (PNR ¼ 34.2 6 1.1) and 45 6 10 motor units (PNR ¼ 34.0 6 0.9) at 30% and 50% MVC, respectively (Fig. 3).At least 82% of the motor units identified in one condition were also identified in the conditions involving a higher number of electrodes.Similarly, 91-100% of the motor units identified in one condition were also identified with the 256-electrode configuration (4-mm IED, 36-cm 2 size; Fig. 1A) with the four grids decomposed separately.
When considering the effect of the size of the grid (IED fixed at 4 mm; Fig. 1A,E-G), we found the lowest number N of motor units with a grid of 2 cm 2 , with 4 6 2 motor units and 4 6 2 motor units at 30% and 50% MVC, respectively (Fig. 5A,C).Additional motor units were then gradually identified with larger grid sizes.The highest number of motor units was observed with a grid of 36 cm 2 , with 56 6 14 and 45 6 10 motor units at 30% and 50% MVC, respectively, with the 4 Â 64-electrode decomposition procedure (Fig. 4A,C).With the 256-Figure 3. Discharge times of the maximum number of motor units identified in one participant (S1) at 30% (A) and 50% maximal voluntary contraction (MVC) (B), with 79 and 58 identified motor units, respectively.The motor units were identified with separated decompositions of the four grids of 64 electrodes (4-mm interelectrode distance).C, Discharge times of the 30 first recruited motor units during the ascending ramp of force (black curve) at 30% MVC (black box in A).
electrode decomposition procedure, 43 6 11 and 25 6 6 motor units were identified at 30% and 50% MVC, respectively (Fig. 4A,C).Finally, we found an increasing logarithmic relationship between the normalized number of motor units N, averaged for each participant, and the size of the grid, with r 2 ¼ 0.99 (p ¼ 3.0•10 À4 ) and r 2 ¼ 0.98 (p ¼ 0.001) at 30% and 50% MVC, respectively (Fig. 5B,  D).It is noteworthy that the parameters of the fits were very similar at 30% and 50% MVC in both analyses.
As both the density and the size of the grid determine the number of electrodes, we finally fitted the relationship between the normalized number of motor units N and the number of electrodes.As observed previously, more motor units were identified with a larger number of electrodes, following a logarithmic tendency with r 2 ¼ 0.98 (p ¼ 0.018) and r 2 ¼ 0.95 (p ¼ 0.016) at 30% and 50% MVC, respectively (Fig. 6).A plateau should theoretically be reached with grids of 1024 and 4096 electrodes (36-cm 2 grids with 2-and 1-mm IED, respectively), with a prediction of 50% and 90% more motor units.
For a fixed number of electrodes, it is noteworthy that the size and the density, although linked, may have different impact on the number of identified motor units (Fig. 6, black crosses).For example, 1.25 times more motor units were obtained with the 64-electrode condition (32 cm 2 , 8mm IED; Fig. 1B) than with the 63-electrode condition (7.7 cm 2 , 4-mm IED; Fig. 1E) for the group of participants at 30% MVC.

Characteristics of identified motor units
We found an increasing logarithmic relationship between the percentage of early recruited motor units for each participant and the density of the grid, with r 2 ¼ 0.91 (p ¼ 2.8 Â 10 À3 ) at 30% MVC (Fig. 7F).Contrary to the density, the size of the grid did not impact the percentage of early recruited motor units, with the percentage ranging from 20% to 29% across all sizes, and the logarithmic trendline returning a negligible slope and a low r 2 ¼ 0.28 (Fig. 7C,G).Such differences were also not observed at 50% MVC, where the percentage of early recruited motor units remained below 10% for all conditions.
To support the above observations made at 30% MVC, grids with the same number of electrodes, but different densities and sizes, were directly compared.62% of the motor units identified with the grids of 64 electrodes (32 cm 2 , IED 8 mm) and 63 electrodes (7.7 cm 2 , IED 4 mm) were identified in both conditions at 30% MVC.28 6 9% Figure 4. Effect of the electrode density on the number of identified motor units N at 30% (A, B) and 50% maximal voluntary contraction (MVC) (C, D).The boxplots in the left column report the absolute number N of identified motor units per participant (gray dots) and the median (orange line), quartiles, and 95% range across participants.In the right column, the normalized number of motor units N logarithmically decreases with interelectrode distance d (4, 8, 12, and 16 mm in abscissa) as . The standard deviation of N across subjects is displayed with vertical bars.Moreover, the quality of the motor unit pulse trains (i.e., decomposition accuracy, estimated by the PNR) increased when increasing the density of electrodes (Extended Data Fig. 4-3 for more details).Two decomposition procedures were considered for the 256-electrode condition; the grid of 256 black electrodes indicates that the 256 signals were simultaneously decomposed and the grid of 256 electrodes of four different colors indicates that four subsets of 64 electrodes were decomposed.To maintain consistency with the computational study, the trendlines were fitted with the 4*64 condition, which returned the higher number of identified motor units (Extended Data Fig. 4-1 for the other fitting condition).It is worth noting that computationally increasing the density of electrodes by resampling the EMG signals with a spatial interpolation did not reveal any previously hidden motor units (Extended Data Fig. 4-2).
Figure 5.Effect of the size of the grid on the number of identified motor units N at 30% (A, B) and 50% maximal voluntary contraction (MVC) (C, D).The boxplots in the left column report the absolute number N of identified motor units per participant (gray dots) and the median (orange line), quartiles, and 95% range across participants.In the right column, the normalized number of motor units N logarithmically decreases with the size of the grid s (2, 3.8, 7.7, and 36 cm 2 in abscissa) as The standard deviation of N across subjects is displayed with vertical bars.Moreover, the quality of the identified motor unit pulse trains (i.e., decomposition accuracy, estimated by the pulse-to-noise ratio) increased when increasing the size of the grid (Extended Data Fig. 4-3 for more details).Two decomposition procedures were considered for the 256electrode condition; the grid of 256 black electrodes indicates that the 256 signals were simultaneously decomposed and the grid of 256 electrodes of four different colors indicates that four subsets of 64 electrodes were decomposed.To maintain consistency with the computational study, the trendlines were fitted with the 4*64 condition, which returned the higher number of identified motor units (Extended Data Fig. 4-1 for the other fitting condition).N of identified motor units at 30% (A) and 50% maximal voluntary contraction (MVC) (B).The discrete results per participant are displayed with gray data points.The average values N per condition are displayed with black crosses.Weighted logarithmic trendlines were fitted to the data and returned (A) , and (B) . Two decomposition procedures were considered for the 256-electrode condition; the grid of 256 black electrodes indicates that the 256 signals were simultaneously decomposed and the grid of 256 electrodes of four different colors indicates that four subsets of 64 electrodes were decomposed.To maintain consistency with the computational study, the trendlines were fitted with the 4*64 condition, which returned the higher number of identified motor units (Extended Data Fig. 4-1 for the other fitting condition).
Figure 7. A, Typical frequency distribution of motor unit force recruitment thresholds in a human TA.The black dashed lines denote the theoretical portions of the population of motor units recruited at 30% and 50% maximal voluntary contraction (MVC).Effect of the grid density (B, D, F) and grid size (C, E, G) on the percentage of early recruited motor units identified at 30% (B, C, F, G) and 50% MVC (D, E).The boxplots report the results per participant (gray dots) and the median (orange line), quartiles, and 95% range across participants.(F) At 30% MVC, the percentage of early recruited identified motor units logarithmically decreases with interelectrode distance d (4, 8, 12, and 16 mm in abscissa) as 44:6 À 13:1log d ð Þ r 2 ¼ ð 0:91; p ¼ 2:8 Á 10 À3 Þ. (G) At 30% MVC, the percentage of early recruited identified motor units does not vary with the size of the grid s (2, 3.8, 7.7, and 36 cm 2 in abscissa), the logarithmic trendline fitting (20:511:2log s ð Þ) returning a negligible slope and low r 2 ¼ 0:28 p ¼ 8 Á 10 À4 À Á .The standard deviation across subjects is displayed with vertical bars.Two decomposition procedures were considered for the 256-electrode condition; the grid of 256 black electrodes indicates that the 256 signals were simultaneously decomposed and the grid of 256 electrodes of four different colors indicates that four subsets of 64 electrodes were decomposed.To maintain consistency with the computational study, the trendlines were fitted with the 4*64 condition, which returned the higher number of identified motor units (Extended Data Fig. 4-1 for the other fitting condition).We did not report the results when five or fewer motor units were identified in one condition for three or more participants.
of the motor units specific to the 8-mm IED grid were early recruited, while 44 6 11% of the motor units specific to the 4-mm IED condition were early recruited.Similar results were obtained with the grids of 35 (36 cm 2 , 12-mm IED) and 34 electrodes (3.6 cm 2 , 4-mm IED), where a higher number of early recruited motor units were specifically identified with denser rather than larger grids.

Correlation between MUAPs from adjacent electrodes
Figure 8 reports the effect of the density of electrodes on the level of correlation r between the profiles of action potentials recorded by adjacent electrodes.The lowest average correlation coefficient r was observed with an IED of 16 mm (r ¼ 0.87 6 0.03 and r ¼ 0.88 6 0.04 at 30% and 50% MVC, respectively).The level of correlation increased with a shorter IED, with r ¼ 0.96 6 0.04 and r ¼ 0.95 6 0.05 between the profiles of action potentials recorded by adjacent electrodes with a 4-mm IED at 30% and 50% MVC, respectively (Fig. 8B,C).
Laboratory study with an ultradense prototyped grid of 256 electrodes with 2-mm IED A total of 31 and 26 motor units (PNR .28 dB) were identified for one participant with the ultradense grid of 256 electrodes (2-mm IED, 9 cm 2 ; Fig. 1H) at 30% and 50% MVC, respectively (Fig. 9B).Note that the signals from four independent subsets of 64 electrodes were decomposed separately.For that participant, more motor units were identified with the ultradense grid of 256 electrodes than with the grid of 64 electrodes covering the same area (Fig. 5A,C).Indeed, 31 and 26 motor units were, respectively, identified at 30% and 50% MVC with the grid of 256 electrodes (Fig. 9C), while 25 (24 6 5 for the group) and 19 (18 6 4 for the group) motor units were identified with the grid of 64 electrodes (Fig. 5A,C).Moreover, fewer motor units were identified when the electrode density of the ultradense grid was decreased (Fig. 9C), with 22 and 13 motor units identified with a 4and 8-mm IED at 30% MVC, respectively, and 21 and 9 motor units identified with a 4-and 8-mm IED at 50% MVC, respectively.At 30% MVC, the rate of increase of N between 4-and 2-mm IED followed the prediction computed in Figure 4B and illustrated by the dash line in Figure 9C.At 50% MVC, the rate of increase of N (Fig. 9C, dotted line) was lower than the prediction.As previously observed, the correlation between adjacent MUAPs increased from r ¼ 0.92 with an 8-mm IED to r ¼ 0.98 with a 2-mm IED at 30% MVC, and from r ¼ 0.85 with an 8mm IED to r ¼ 0.93 with a 2-mm IED at 50% MVC (Fig. 9A).All the motor units identified with the 8-and 4-mm IED were also identified with the 4-and 2-mm IED grids, respectively.Finally, more motor units with an early recruitment were identified when increasing the density from 8-to 4-mm IED (Fig. 9B, blue vs black trains), and from 4-to 2-mm IED (Fig. 9B, red trains).

Discussion
This study systematically investigated how the design parameters of grids of surface EMG electrodes (grid size and electrode density) impact the number and the properties of the motor units identified with EMG decomposition.Using a combination of computational and experimental analyses, we found that larger and denser grids of electrodes than conventionally used reveal a larger sample of identified motor units.As most of the motor units that were not identified with less dense and smaller grids had an early recruitment threshold, we concluded that denser grids allow to identify smaller motor units.This is because of a better spatial sampling of MUAPs over the grid, which in turn improves the discrimination of motor units with a unique set of MUAPs among active motor units.These results clarify the direction for designing new grids of electrodes that could span across the entire surface of the muscle of interest while keeping a high density of electrodes, with IED as low as 2 mm.Identifying large sets of small and large motor units is relevant in many research areas related to motor control, such as the investigation of synergies (Hug et al., 2023), neuromuscular modeling (Caillet et al., 2022c), or human-machine interfacing (Farina et al., 2023).
The number N of identified motor units increased across participants with the density of electrodes (Figs. 4,8C), the size of the grid (Fig. 5), and the number of electrodes (Fig. 6).On average, 30 and 19 motor units were identified with the "conventional" 64-electrode grid (8-mm IED, 32 cm 2 surface area) at 30% and 50% MVC, respectively, which is consistent with several previous studies using similar grid designs (Del Vecchio et al., 2020).By increasing the density of electrodes and size of the grid to reach a total of 256 electrodes separated by a 4-mm IED, we identified on average 56 and 45 motor units at 30% and 50% MVC, respectively.We even reached 79 and 59 Figure 9. Results for the ultradense prototyped grid (2-mm interelectrode distance (IED), 5 Â 1.8 cm, 256 electrodes).A, Description of the ultradense grid, where gray circles represent the electrodes.On average, the correlation between the profiles of MUAPs detected over electrodes separated by an IED of 2 mm (orange), 4 mm (blue), and 8 mm (purple) reached r ¼ 0.98, 0.96, and 0.92 at 30% maximal voluntary contraction (MVC), respectively, and 0.93, 0.88, and 0.85 at 50% MVC, respectively.B, Series of discharge times for motor units identified at 30% (left) and 50% MVC (right).The dark ticks represent the discharge times identified with a grid of electrodes with an 8-mm IED.The discharge times in blue were additionally identified with a grid of electrodes with a 4-mm IED, and the discharge times in red were additionally identified with a grid of electrodes with a 2-mm IED.All the pulse trains identified with one grid were also identified with the denser grids.C, Effect of electrode density on the number of identified motor units at 30% (scatters) and 50% MVC (triangles).The trendlines from the density analysis in Figure 4B,D are also reported (red dotted lines).To maintain consistency with the other results, the grid was decomposed as four independent subsets of 64 electrodes, as explained in Materials and Methods, to identify the higher number of motor units.
motor units for one subject (Fig. 3), which is substantially more than the numbers of motor units usually reported in studies with similar methods, and twice those obtained with grids of 64 electrodes in this study.Our computational and experimental analyses showed that the size of the grid is a key factor contributing to the higher number of identified motor units (Figs.2B, 5).According to our simulations, increasing the size of the grid increases the number of theoretically identifiable motor units, i.e., the number of motor units with unique sets of MUAPs across electrodes (Fig. 2B).These differences between MUAPs result from the anatomic and physiological differences between adjacent motor units, such as the length of their fibers, the spread of the end plates, or their conduction velocity, as well as from the properties of the tissues separating the fibers from each recording electrode (Farina et al., 2014).Larger grids better sample these differences across electrodes, revealing the unique profiles of each motor unit action potentials (Farina et al., 2008).The density of electrodes was also a critical factor to increase the number of identified motor units (Figs. 4, 9C).Dense grids especially allowed to better identify early recruited motor units.Classically, the decomposition algorithms tend to converge toward the large and superficial motor units that contribute to most of the energy of the EMG signals (Farina and Holobar, 2016).Conversely, action potentials of the smallest motor units tend to have lower energy and are masked by the potentials of the larger units.These factors explain the lowest representation of low-threshold motor units in available HD-EMG datasets (Caillet et al., 2023).Increasing the density of electrodes would therefore enable to better sample the action potential profiles of these early recruited motor units across multiple electrodes, enabling their identification.However, we observed that increasing the density did not reveal additional early recruited motor units during contractions at 50% MVC (Fig. 7D).This is potentially because of the higher energy of the MUAPs of the motor units recruited between 30% and 50% MVC.Additionally, we also showed in one subject that synthetically increasing the density of electrodes by resampling EMG signals with spatial interpolation does not have the same effect as with denser grids.In this example, 4 and 19 motor units were identified from the interpolated grid with a 4-and 2-mm IED, respectively, versus 19 and 24 motor units with the experimentally recorded signals.All the motor units identified with the interpolated grid were also identified with the experimentally recorded signals (Extended Data Fig. 4-2).
The number of identified motor units N monotonically increased with the density of electrodes (Fig. 4B,D), the size of the grid (Fig. 5B,D) and the number of electrodes (Fig. 6), following significant logarithmic trendlines.Remarkably, very similar logarithmic tendencies were obtained at both 30% and 50% MVC in all the analyses.Altogether, these trendlines suggested that the normalized number of identified motor units N would grow with an electrode density beyond a 4-mm IED.We experimentally tested this hypothesis by designing a new prototyped grid of 256 electrodes separated by an IED of 2 mm.As predicted, more motor units were identified with a 2-mm than with a 4-mm IED, following at 30% MVC the same rate of increase as predicted by the logarithmic trendlines (Fig. 9C) between 4and 2-mm IED.This increase may plateau with higher electrode densities, as the level of correlation between the profiles of MUAPs detected over adjacent electrodes tended to 1 (Fig. 9A).Therefore, the high level of similarity between signals recorded from adjacent electrodes in ultradense grids (IED , 2 mm) may limit the percentage of identifiable motor units (Farina and Holobar, 2016).According to these results, we consider that optimal designs of surface grids of electrodes for identifying individual motor units would involve a surface that covers the muscle of interest with an IED as low as 2 mm.
Another important factor for the accuracy of the discharge times estimated for each individual motor unit is the quality of the motor unit pulse trains, estimated by the PNR (Holobar et al., 2014) or the silhouette value.In this study, we found that the quality of the identified motor units (i.e., decomposition accuracy) increased when increasing the density of electrodes or the size of the grid, with PNR reaching on average 37-38 dB across participants with the grid of 256 electrodes (Extended Data Fig. 4-3).A greater average PNR implies the need of less manual editing following the automatic decomposition (Hug et al., 2021b).The better estimates of motor unit pulse trains depend on the better signal-to-noise ratio following the inversion of the mixing matrix, since the pulse train of each motor unit is computed by projecting the extended, whitened signals on the separation vector (Holobar and Farina, 2014;Farina and Holobar, 2016;Negro et al., 2016).Likewise, the PNR substantially increased after we computationally increased the number of electrodes by spatially resampling the EMG signals.This practical result is of interest for most of the physiological studies that require a lengthy processing time to visually inspect and manually edit the discharge times estimated from the pulse trains of all the motor units (Hug et al., 2021b).
Finally, we increased both the total number and the percentage of early recruited motor units identified by independently decomposing subsets of 64 electrodes within the grids of 256 electrodes, compared with the simultaneous decomposition of all available observations (Fig. 7B,  C).This was likely because of the lower ratio of large motor units sampled by each subset of electrodes, allowing the algorithm to converge to smaller motor units that contributed to the signal (Fig. 7B,C).Importantly, it should be noted that the simulation results were obtained independently of a specific decomposition algorithm, as previously proposed by Farina et al. (2008).On the other hand, the experimental results are based on a specific algorithm.Interestingly, however, the simulation and laboratory results were fully consistent and in agreement, indicating that the difference in shape of the spatially sampled MUAPs is the main factor influencing EMG decomposition.
In conclusion, by increasing the density and the number of electrodes, and the size of the grids, we increased the number of theoretically identifiable and experimentally identified motor units from the surface EMG signals.The identified motor units had pulse trains with high PNR, limiting the manual processing time.Moreover, we identified a higher percentage of early recruited motor units, which are classically filtered out with the conventional grid designs.In this way, a maximum of 79 motor units (PNR .28 dB; mean: 36 dB), including 40% of early recruited motor units, were identified, which is substantially greater than the samples previously reported with smaller and less dense grids.From these results, we encourage researchers to develop and apply larger and denser EMG grids to cover the muscle of interest with IEDs as small as 2 mm.This approach should increase the sample of motor units that can be experimentally investigated with noninvasive techniques.

Figure 6 .
Figure 6.Effect of the number n of electrodes on the normalized numberN of identified motor units at 30% (A) and 50% maximal voluntary contraction (MVC) (B).The discrete results per participant are displayed with gray data points.The average values N per condition are displayed with black crosses.Weighted logarithmic trendlines were fitted to the data and returned (A)N ¼ À104137log n ð Þ r 2 ¼ 0:98; p ¼ 0:018 À Á, and (B)N ¼ À113138log n ð Þ r 2 ¼ 0:95; p ¼ 0:016 À Á.Two decomposition procedures were considered for the 256-electrode condition; the grid of 256 black electrodes indicates that the 256 signals were simultaneously decomposed and the grid of 256 electrodes of four different colors indicates that four subsets of 64 electrodes were decomposed.To maintain consistency with the computational study, the trendlines were fitted with the 4*64 condition, which returned the higher number of identified motor units (Extended Data Fig.4-1 for the other fitting condition).

Figure 8 .
Figure8.Effect of the electrode density on the correlation r between the profiles of motor unit action potentials (MUAP) detected over adjacent electrodes (A) at 30% (B) and 50% maximal voluntary contraction (MVC) (C).The profile of the MUAP detected over the red electrode was compared with those detected over the four adjacent electrodes separated by a 4 mm (orange), 8 mm (blue), 12 mm (green), and 16 mm (purple) interelectrode distance (A).The boxplots denote the correlation coefficient r per participant (gray dots) and the median (orange line), quartiles, and 95% range across participants.