Electrical cross-sectional imaging of human motor units in vivo

OBJECTIVES
In many neuromuscular diseases, weakness results from a disruption in muscle fibres' arrangement within a motor unit. Limitations in current techniques mean that the spatial distribution of fibres in human motor units remains unknown.


METHODS
A flexible multi-channel electrode was developed and bonded to a clinical electromyography (EMG) needle. Muscle fibre action potentials were localised using a novel deconvolution method. This was tested using simulated data, and in recordings collected from the tibialis anterior muscle of healthy subjects.


RESULTS
Simulated data indicated good localisation reliability across all sections of the electrode except the end sections. A corrected fibre density was estimated up to 1.4 fibres/mm2. Across five recordings from three individuals, between 4 and 14 motor units were detected. Between 1 and 20 muscle fibres were localised per motor unit within the electrode detection area, with up to 220 muscle fibres localised per recording, with overlapping motor unit territories.


CONCLUSIONS
We provide the first direct evidence that human motor units spatially overlap, as well as data related to the spatial arrangement of muscle fibres within a motor unit.


SIGNIFICANCE
As well as providing insights into normal human motor physiology, this technology could lead to faster and more accurate diagnosis in patients with neuromuscular diseases.


Introduction
The spatial distribution of muscle fibres within a motor unit is of fundamental interest in both healthy and diseased muscle. To date, the only means of accurately demonstrating the spatial distribution of muscle fibres has been the glycogen depletion technique (Edstrom and Kugelberg, 1968). These complex experiments can only be performed in animals since they require the histological analysis of explanted muscle cross-sections exposed to highfrequency electrical stimulation for several minutes. These experiments demonstrate that individual muscle fibres within a motor unit are widely separated (Bodine et al., 1988;Bodine-Fowler et al., 1990), and that a single motor unit can cover up to 50% of the cross-sectional area of the muscle (Buchthal et al., 1959). From this, it has been inferred that the fibres in one motor unit overlap with several other units (Fig. 1a). However, since only a single motor unit can be studied per muscle, this has never been shown directly. Furthermore, since only one muscle can be analysed per animal, most studies using the glycogen depletion technique describe 5-10 motor units at best. Consequently, the statistical distribution of muscle fibres within a normal motor unit remains open to debate with random (Edstrom and Kugelberg, 1968), regular (Kugelberg, 1973) and clustered (Buchthal et al., 1957) models having been proposed. Further neurophysiological investigations of the motor unit include determining the territory of a motor unit with a linear multielectrode (Buchthal et al., 1959), the electrical activity of the whole motor unit via macro EMG (Stålberg and Fawcett, 1982), and the statistical distribution of muscle fibres through the territory of the motor unit via scanning EMG (Stålberg and Dioszeghy, 1991).
Changes in the spatial distribution of muscle fibres within motor units are among the earliest changes in many neuromuscular diseases. For example, in Amyotrophic Lateral Sclerosis (ALS), progressive loss of motor nerve axons removes the nerve supply to those muscle fibres in the dependent motor units (Nijssen et al., 2017). Compensatory re-innervation by the surviving axons results in larger motor units containing clusters of closely-packed muscle fibres (Schwartz et al., 1976). Although these changes can be shown using the glycogen depletion technique in animal models, this is only possible in human subjects in a limited capacity via fibre density which measures the number of single fibre action potentials within the limited sampling territory of a single fibre electrode (Sanders et al., 2019) and does not indicate relationship of fibres to one another. Recent developments such as motor unit MRI (Birkbeck et al., 2020) have permitted the size, shape and dimensions of individual human motor units to be determined, but the relatively low spatial resolution (approximately 1 mm 2 voxel size) makes it impossible to detect the position of individual muscle fibres using this method.
Gross changes in motor unit structure can be inferred from intramuscular electromyography (EMG). Since only a single recording surface at the tip of the needle is used (Whittaker, 2012), the result-ing motor unit potential (MUP) represents the temporally and spatially summated activity of every muscle fibre within approximately 2.5 mm of this surface (Dumitru et al., 1997). Not surprisingly, attempting to infer changes in muscle fibre spatial distribution from this composite 1-dimensional signal requires expert interpretation, and even in the best hands, the technique has only moderate diagnostic accuracy and specificity (Narayanaswami et al., 2016). The small sampling area also requires the needle to be moved to several positions within the muscle, increasing the time taken and the discomfort for the patient.
Early attempts to increase the sampling area with linear multielectrode arrays with up to 14 recording surfaces (Buchthal et al., 1957) have not entered routine use due to the large crosssectional area (up to 4x a conventional EMG needle) and the high manufacturing costs, requiring them to be re-sterilised. Furthermore, limitations in recording apparatus meant that only 1 or 2 channels of data could be recorded simultaneously. Hence, although individual fibres' position along the axis of the electrode could be estimated, the spatial distribution of these fibres within the motor unit cross-section could not be determined (E Stålberg et al., 1976).
Substantial improvements in different technical fields over the last 20 years have now made it feasible to address these limitations. Electrodes have been developed with a high density of recording contacts over their shaft (Angotzi et al., 2019;Farina et al., 2008;Jun et al., 2017;Musk, 2019;Rousche et al., 2001;Sohal et al., 2016;Vetter et al., 2004); these have mainly been used in the central Muscle fibres (red cylinders) run approximately parallel to the long axis of the muscle and are grouped into motor units (pink circles). Axis definitions reference for the geometry referred to in all figures and equations. X axis indicates along the direction of the needle, and Y axis is perpendicular to this, radial to the muscle fibre direction. Z axis is perpendicular to the needle and runs along the direction of the muscle fibres. Right; spatial distribution of fibres in a single motor unit revealed by glycogen depletion. Non-depleted fibres stain purple. Depleted fibres (pink areas *) are widely spaced and interdigitate with fibres from several other motor units. (B) Upper panel: 32-channel flexible MicroEMG array comprising two 10 mm thick parylene layers encasing tungsten:titanium conductor prior to bonding. Lower panel: micrograph of MicroEMG array bonded to a conventional EMG needle. A 400 mm wide zig-zag array of 32 recording electrodes (purple circles) each of 25 mm diameter extends 12,800 mm along the needle. (C) Sterile MicroEMG probe interfaced with INTAN 32-channel amplifier via a custom-made PCB. (D) Upper panel: system in use to obtain motor unit recordings from the tibialis anterior muscle of a healthy volunteer. Lower panel: conventional EMG (upper) records muscle fibre action potentials within $1 mm of the needle tip. MicroEMG probe records from up to 32 electrodes spanning a 12.4 Â $2 mm corridor along the needle shaft.
nervous system for recording neural activity, but similar techniques have also been applied to muscle (Muceli et al., 2015). Integrated circuits have been developed incorporating multiple amplifiers and associated digitisation circuitry, allowing multiple channel recording with miniaturised equipment. Finally, computer storage and processing capacity has greatly expanded, making it feasible to deal with high data volumes. Here, we have exploited these advances to develop multi-electrode electromyography arrays with up to 32 contacts capable of simultaneously recording single fibre potentials from multiple motor units. Using this novel system, we can localise up to 220 muscle fibres with $100 mm accuracy in human muscle, in effect producing an electrical cross-sectional image of motor unit structure. Such images provide, for the first time, direct evidence for the overlapping arrangement of motor units, and hold great promise for clinical application.

Probe design
We designed flexible multi-electrode arrays (which we named 'MicroEMG') comprising 32 recording surfaces sandwiched between two 10 mm-thick layers of parylene-C, a biocompatible polymer with sufficient tensile strength to survive insertion through the skin without breakage ( Fig. 1B upper panel) (Sohal et al., 2014). We used a tungsten:titanium (W:Ti) alloy for the electrode surfaces and interconnecting tracks (W:Ti 80:20 ratio) to provide greater flexibility than tungsten alone. Electrode arrays were fabricated on 3-inch diameter silicon wafers using standard semiconductor processing techniques. The recording electrodes were 20 mm Â 30 mm stadium shapes to match the recording area of clinical single-fibre EMG needles. The recording electrodes were arranged into 2-parallel rows with 400 mm spacing and offset by 200 mm to form a zig-zag pattern ( Fig. 1B lower panel). From the first to the last recording electrode, the total length of the zig-zag array was 12.40 mm. These dimensions were chosen such that the array was long enough to record several motor unit cross sections simultaneously from a single insertion, but with sufficient spatial resolution to discriminate potentials recorded from individual muscle fibres (50-100 mm diameter). Each recording electrode was connected to a bond-pad at the top of the flexible array via W:Ti tracks. The width of each track covered as much of the available area in the array as possible in order to minimise the chance of breakage, thinning progressively from the tip until they were 3 mm wide with 2.8 mm spacing as they entered the bond-pad region. To allow these flexible arrays to penetrate muscle, they were bonded to the shaft of a 480 mm diameter clinical EMG needle (TECA Elite, Natus Neurology, USA) using a 1180-M UVcurable biocompatible adhesive (Dymax MD). The 20 mm thick arrays occupied approximately 1/3 of the needle's circumference (Fig. 1B), creating an approximately circular probe with diameter $500 mm. Bench-top electrical testing showed conduction in 330 out of 352 measured recording sites, a yield of 94% working channels. The average impedance in the functional channels was 502 kO ± 109 kO measured at 1 kHz, and no leakage current was detectable between adjacent channels (with a current measurement resolution of 10 nA).
The bondpad region of the probe was clamped against the contacts of a custom-made printed circuit board housing a Intan RHD2132 32-channel low-noise amplifier chip with integrated analogue-to-digital converter which sampled each channel at 20 kHz with 16-bit precision (Intan Technologies, USA) (Fig. 1C).
The needle array and subsequent signal processing was designed for a position perpendicular to the muscle fibres, with the plane of the electrode grid arranged in a radial orientation (Fig. 1D).

Signal processing
Pre-Processing. Signals underwent band-pass filtering between 500 Hz-10KHz. Residual mains noise was removed using an algorithm that removed components synchronous to 50 Hz (Riddle and Baker, 2005). Broken channels were identified visually by very low amplitude signals and excluded from further analysis.
Motor unit extraction. Motor unit potential templates (MUPTs) were extracted from complex multi-channel data using a template matching algorithm, which utilised a multiresolution Teager energy operator employing pseudocorrelation (Sedghamiz, 2018;Sedghamiz et al., 2015). The channel with the highest SNR was selected for MUPT extraction, and the onset time of each potential recorded ( Fig. 2A). Templates were inspected visually, and those which appeared to be noise transients on the basis of morphology or firing frequency were removed from further analysis. For each motor unit, a 10 ms window was selected either side of the onsets of each MUPT, and across all the working electrodes of the array (Fig. 2B) to produce the multi-electrode EMG for each MUPT.

Muscle fibre localisation
In order to reduce noise, ten consecutive motor unit discharges were averaged for analysis, effectively forming a moving average window. The averaged multi-electrode signal was processed with a top-hat transform (Maire et al., 2000) to improve the delineation of transmembrane currents with similar locations and low interspike intervals (ISIs) while avoiding noise transients being incorrectly identified as fibre action potentials. Peaks in the rectified multi-electrode signal amplitude (Fig. 2C) representing the electrodes with minimal radial distance to a fibre were identified using a peak-finding algorithm (Natan, 2013). Peaks were analysed in both time and space. The five EMG channels (an electrode 'cluster' as shown in Fig. 3) surrounding each peak were selected. For each cluster, the 0.5 ms of EMG before and after each peak was used in the next step of fibre localisation.
The fibre localisation methodology was based upon the model used by (Nandedkar and Stålberg, 1983) to simulate single fibre action potentials (SFAP) recordings using a point electrode. This transmembrane current generating function (G z ) used by Nandedkar & Stålberg was derived from a volume conduction model (Rosenfalck, 1969) and was used to estimate the SFAP that would be recorded at a specific location. As the multi-electrode allows for the recording of the SFAPs at multiple sites near the active fibre, the muscle fibre transmembrane current generating function could be estimated. The method used the convolution of the simulated generator function (Gz) and a weighting function (W).
Where: z = axial separation r = radial separation I = strength of current source r r = radial conductivity r z = axial conductivity K = r z =r r We defined the x axis as being parallel to the needle, and y axis as being perpendicular to the needle along the plane of the electrode surfaces. Further review of the geometry involved is described in Fig. 1A. The weighting function defined by Eq. (1) describes the potential that would be recorded by a point electrode located at a radial distance r and axial distance z from a cylindrical current source/sink on the surface of the fibre. If the transmembrane current generating function (G z ) can be considered to be a propagating transmembrane current function moving at a constant velocity V along the z axis, this becomes a function in time t with distance d.
The recorded SFAP V ðr;tÞ at a radial distance r and at some arbitrary location of z, is: where ⁄ is the convolution operator, and the transmembrane current passes an electrode at coordinates r, z at t = 0.
Using a multi-electrode, there are additional electrodes lying on a line of constant z at varying radial distance depending on the electrode separation from the muscle fibre. The SFAP recorded at each of these n sites may be used such that for each electrode, the deconvolution of the weighting function and the recorded SFAP should result in the same transmembrane current function.
where r(n) is the radial distance between the fibre and the nth electrode and is the deconvolution operator. The estimation is achieved by searching the solution surface in Â and y to minimise the variance of the estimates of G z . The deconvolution was implemented by the matrix method; in this case, the calculated weighting function for each electrode was constructed as a Toeplitz (diagonal constant) matrix. The matrix contains a sequentially time-shifted vector of the weighting function (Eq. (1)). This type of matrix has a well-behaved inverse (W À1 Þ such that the estimated G z for each electrode can be calculated by: The radial distance term can be transformed to Cartesian coordinates in x,y. Using a suitable range of x,y coordinates, a solution surface may be defined with an error term: Overview of signal analysis pipeline. After acquisition of multi-channel EMG signals (outer window is a multi-channel motor unit potential (MUP), inner window is a single channel MUP) (A) MUPs were extracted from the channel with the highest SNR (marked with *), MUPs were each averaged to produce a multi-electrode MUP (distance indicates x-axis distance) (B). Peaks in signal amplitude around the needle electrode (illustrated to simplify from 3 dimensions to 2 in C) were chosen as the starting point for muscle fibre localisation. Putative fibre location is indicated by the yellow square. Fibres were localised in space on a moving-average window across multiple firings. Fibre positions were tracked through time (Three motor units shown in D) and error values for localisation accuracy were calculated prior to k-means clustering of fibre positions. The electrodes and the resulting deconvolved electrode potentials are shown in red, blue and green. A search was performed through possible fibre locations near to the recording electrode cluster. At each location, the Euclidean distance between the putative fibre and nearby recording electrodes was calculated. The potentials recorded at these electrodes were then deconvolved with a tissue function for these distances (R1,2,3) to produce n putative generator waveforms. The variance across these generator waveforms was then calculated, and the location of the fibre defined as the position with minimum variance. E ðrðnÞÞ ¼ 1=ðmean varðG ðzþVt;nÞÞ Þ An unconstrained non-linear optimisation algorithm (MATLAB function fminsearch) was used to find the value of the x,y coordinates which minimised the variance of the convolved multielectrode potentials for each group of moving averaged motor unit discharges (Fig. 2). A value of 5 was used for K as suggested by (Nandedkar and Stålberg, 1983). The variance function was constructed to measure the variance at each time point, and then measure the average of the variances across the window.

Fibre data visualisation
To summarise this localisation across time, the average number of fibres found at each motor unit discharge was calculated (Fig. 2D). This average fibre count was then used as k to produce a k-means clustering of all the fibre positions. The median position of each cluster was taken to be the normative position of the fibre in calculating fibre density and nearest neighbour distances. The localisation estimate confidence was calculated using the mean consecutive position change. This would mitigate any 'slow drift' of needle position over the 5-minute recording related to operator technique, compared to large 'jumps' in position related to fibres within poor quality fibre potentials. For each motor unit, the distance between each localised fibre and its nearest neighbour was calculated and displayed as a cumulative probability plot. The fibre density was calculated based upon the fibre locations within the area 1 mm off the needle's axis (defined as the point along the x axis through the midpoint of the two rows of electrodes with y = 0). The current clinical definition of fibre density is the mean number of fibre action potentials within the uptake radius of a single fibre electrode, with position optimised to be as close to a single fibre as possible, across twenty fibres (Stålberg, 1979). Estimates of this uptake radius vary, but range from 236-421 mm (Gath and Stalberg, 1979). We aimed to produce a comparable measurement based on the local fibre arrangements, and so measured the average number of fibres within an arbitrary 400 mm distance of each localised fibre, per motor unit (herein termed the needle fibre density). The motor unit territory was measured as the maximum Euclidean distance between any two fibre pairs. Both measures were again confined to muscle fibres localised to within the area 1 mm off the needle's axis.

Modelling
In silico experiments were performed using simulated data to test the likely performance of our electrodes and our fibrelocalisation algorithm. We modified a previously published model (Stashuk, 1993), which defines interdigitated motor unit territory, recruitment & firing rates, and realistic attenuation of EMG signal through tissue. Motor units were activated with inter-discharge intervals selected at random from a gamma distribution (Stein, 1965) of 9th order at a rate of around 10 Hz. The estimated transmembrane current potential was based upon the product of transmembrane currents, using the linear quadrupole model Stålberg, 1983, Andreassen andRosenfalck, 1981). We assumed that each fibre was cylindrical and orientated perpendicular to our array, in a medium with cylindrical anisotropy. We adjusted this model's parameters to be physiologically faithful in terms of the change of signal through increasing distances, with most of the signal attenuated at distances exceeding 1 mm, as previously described in single fibre experiments (Stalberg and Antoni, 1980). Gaussian and mains noise were then added to the signal to achieve a signal to noise ratio of 5-10 dB.
All results were computed using a consumer laptop (MacBook Pro Intel i5 processor, 2016). The process took between 10 and 25 minutes depending upon the length of the recording and number of motor unit discharges.
Single fibre localisation accuracy. A single fibre was simulated in positions across a 7 mm Â 1.5 mm (10.5 mm 2 ) grid in increments of 50 mm surrounding a simulated 32 electrode EMG needle. Error was measured by the Euclidean distance between the predicted and actual location.
Dual fibre resolution. Two fibres were modelled within a 0.8 mm range of the needle's axis (defined as the midpoint between the two rows of electrodes), at radial separation of between 50 mm and 1.5 mm (50 mm steps) and random orientations. This was repeated with ISIs of 0.5 ms, 1 ms, and 2 ms. For each linear distance, 100 trials were performed, for a total of 9,000 trials. The estimated number of fibres was recorded for each trial.
Multiple fibre count. Between 5 and 20 fibres were placed with positions selected from the uniform distribution within the central 5 mm Â 1.4 mm (7 mm 2 ) search space. Fibre positions were independently selected i.e. there was no constraint on their proximity to one another. This central region was chosen to mitigate the edge effects of recording at the needle's furthest spatial extent. The mean nearest neighbour distance varied from 840 mm to 320 mm. Fifty trials were performed at each fibre count, with recordings of 10-second duration simulated and the estimated fibre count compared to the modelled number of fibres.

Human recordings
Studies were approved by the Ethics Committee of the Newcastle University Medical Faculty. Five recordings were made from the tibialis anterior muscles in 3 healthy controls, none of whom had any history of neuromuscular disease.
Signals were amplified (gain 200) and sampled to disc at 20 kSamples/s/channel using Intan software, together with the analogue output of a force transducer. Digital data were transmitted over a low voltage differential signal (LVDS) serial peripheral interface (SPI) bus to an interface card and then to a laptop computer. This stored data on a hard disc for later off-line analysis using customwritten MATLAB code (MathWorks, USA). The system was battery powered to ensure mains isolation, and a skin reference was used.
Subjects reclined on a couch with the force output at the ankle joint recorded using a Loadcell Amp PRO force transducer (Elane Electronics, USA). Probes were connected to the amplifier headstage at the bedside, and the probe was inserted into the muscle approximately perpendicular to the skin, with the plane of the electrode grid arranged in a radial orientation (Fig. 1D). Tibialis anterior was selected since the orientation of its muscle fibres is known, allowing the probe to be oriented such that the fibres ran approximately perpendicular to the plane of the array. Subjects were asked to maintain 10% of the maximum voluntary contraction (MVC) strength under visual feedback for between 1 and 5 minutes. Subjects reported mild discomfort as the probe was inserted, and once in place the probe was held still within the muscle within a single position. Recordings were well-tolerated, and no bruising or bleeding was observed.

Results
Single Fibre Localisation. Across all positions surrounding the electrode array, fibres were localised with a mean error of 176 mm ± 28 mm (Fig. 4A). An increased error was seen near the ends of the electrode array, where fewer electrode potentials were within the recording radius of the fibre to allow accurate deconvolution. Fibres within 1 mm radius of 4 different recording electrodes were accurately localised to < 100 mm. However, when a fibre was placed more than 0.8 mm from the needle's axis, the transmembrane current potential was often indistinguishable from noise, and so no fibre was detected or localised.
Dual Fibre resolution. At ISIs of 2, 1 and 0.5 ms, the fibre number estimate was correct in 99.03%, 98.23% and 82.6% of trials respectively over all fibre spatial separations (50 mm-1.5 mm). At higher ISIs of 1 and 2 ms, fibre spatial separation did not affect the error rate, while for the lower ISI of 0.5 ms, this rate was 91.04% for fibres spatially separated by over 0.5 mm. At separations below 0.5 mm, the accuracy rate fell to 59.5% (Fig. 4B).
Multiple fibre localisation. The fibre number was reliably estimated up to 10 simulated fibres (a fibre density of 1.4/mm 2 ). Above this level, the fibre density was underestimated (R 2 = 0.91). For example, at 15 simulated fibres (density 2.1/ mm 2 ), a mean of 13.94 fibres (±0.16) were detected, while at 20 simulated fibres (density 2.8/mm 2 ), a mean of 17.26 fibres (±0.27) were detected (Fig. 4C). This reflected the increasing probability that muscle fibres and their respective transmembrane currents were coincident in tie and space. The algorithm appeared robust, with limited variability between simulations (the average standard deviation across all trials was 1.06 fibres).

Experimental results
Constant force level recordings. Across five recordings from three individuals, between 4 and 14 motor units were identified per recording. For each recording, the median signal: noise ratio was calculated over all channels using the Spurious Free Dynamic Range (Adimulam and Srinivas, 2016;Pozar, 2001). This median varied from 39.6 to 177 (range 20.3-308.7).
Motor units fired at frequencies between 7.52 and 16.11 Hz. We were able to localise up to 20 muscle fibres per motor unit (range 1-20), up to a total of 220 fibres per recording (Fig. 5C). There was no correlation between the SNR of a recording and the number of fibres localised (r = 0.14, p = 0.384).
The median nearest neighbour distances of muscle fibres within a motor unit varied across all of the recordings between 0.54 and 1.21 mm. There was also no correlation between the SNR of a recording and the median nearest neighbour distances (r = -0.09, p = 0.57). The mean fibre density per motor unit varied between 0.625 and 1.8 fibres/mm 2 , while the mean clinical needle fibre density measure was 1.22 (SD ± 0.29) across all motor units and all recordings. The motor unit territory size varied between 2.29-9.72 mm (median 8.05 mm).
Examining the spatial distribution of fibres within motor units, many motor units were seen to be overlapping-no areas were seen with only a muscle fibre population from a single motor unit. (Fig. 5A). While on some recordings there were 'silent areas' of the electrode array without any muscle fibres, these did not correspond to areas of high impedance nor low SNR.
Fibres were sometimes localised to within the dimensions of the array itself, (e.g. Fig. 5A, second row). Where only a low number of fibres were localised from a motor unit (Recording 1 & 4), the fibres tended to localise with low radial distance to the needle's axis. To estimate the statistical distribution of fibres within the muscle, the fibre positions across all recordings were overlaid onto a single map and converted to an absolute position from the midpoint centre of the needle electrode (coordinates 5,0 on Fig. 5A). The fibres were seen to be uniformly distributed up to 4 cm in the x-axis (parallel to the needle), and 0.75 mm in the y axis (per-pendicular to the needle) (Fig. 6). This hypothesis was formally tested using the Kolmogorov-Smirnov test. In both Â and y directions, the test failed to reject the null hypothesis, with p-value of Â = 0.086, and p-value of y = 0.0714. The fibre spatial characteristics are summarised in Table 1.

Discussion
We present a novel intramuscular multi-electrode array capable of determining the firing dynamics and location of up to 220 muscle fibres in a human muscle in vivo from a single insertion. The arrays are of sufficient length to allow simultaneous crosssectional imaging of multiple motor units from both deep and superficial regions of the muscle.
Until now, the only means of directly determining muscle fibre spatial distribution in single motor units has been the glycogen depletion method. Although variations of these experiments have been performed in humans (Garnett et al., 1979), it remains the case that fibre spatial distribution in single human motor units cannot be estimated using this technique. Conventional electromyography can provide information on human motor unit structure, and by using a variety of needle electrodes with spatial scales ranging from individual muscle fibres, i.e. single fibre EMG (Ekstedt and Stålberg, 1971) to entire motor units, e.g. macro EMG (Stålberg, 1980) or scanning EMG (Dieszeghy and Diószeghy, 2002), it is possible to gain information on the density of fibres within motor units and their spatial extent. However, these techniques are time-consuming and require multiple insertions into the muscle to provide sufficient sampling. More importantly, since only a single recording surface is used, the resulting voltage vs time signal represents the spatially and temporally summated activity of every muscle fibre within the pick-up area of each electrode. It is therefore a priori impossible to deconvolve this complex 1-dimensional Motor Unit Potential into the constituent 2-dimensional pattern of muscle fibre potentials.
Modelling studies: Because the spatial distribution of muscle fibres in human motor units is not known, we had no gold standard against which to compare the results of our human experiments. We therefore modelled our fibre-localisation algorithm's performance using single, paired, or multiple muscle fibres at a variety of biologically relevant spatial separations and inter-spike intervals. Our simulations suggested that fibres would be accurately localised within a corridor extending approximately 0.8 mm on either side of our arrays, at temporal separations above 0.5 ms and radial spatial separations above 0.5 mm. Below these separations, our algorithm could fail to separate fibres, tending to under-estimate the true fibre density. Fibre localisation was least accurate at the tips of the electrode, where muscle fibre potentials were recorded by fewer electrodes.
Human recordings: Using our arrays, a large number of motor units and their constituent muscle fibres could be studied from relatively brief recordings. Accurate localisation of up to 220 fibres could be obtained in vivo from a 60-second recording without the need for needle movement. Critically, by localising fibres from multiple motor units simultaneously, the relationship between fibres in neighbouring units could be determined. We found that fibres from multiple motor units overlapped in any given region of the array. Although an overlapping pattern has long been inferred from glycogen depletion and EMG studies, ours is the first direct evidence of this arrangement in either an animal model or a human. Previous studies of human motor units have estimated between 10-15 motor units within a given the pick-up area of a concentric EMG electrode (Edstrom and Kugelberg, 1968), and our data of 4-14 motor units per recording is consistent with this.
One important limitation in these early experiments relates to the relatively high noise level encountered during recordings. This had a number of impacts. Firstly, the narrower filter range may have distorted the shape of the recorded electrode potentials used for localisation. This is because muscle tissue acts to attenuate higher-frequency components of a signal, and therefore the highpass filter may reduce the contributions of the potential at greater z-distances from the electrode, along the length of the fibre, affecting the weighting function described in Eq. (1). Secondly, several noise transients were isolated via the MUPT extraction process, and these necessitated manual removal. In the future we would seek to improve the noise environment for future iterations.
Depending upon the reliability of the MUP template matching used, it is also possible that MUPTs may have been either split or merged. If split, the MUPT would be expected to produce additional motor units with constituent fibres located at similar positions around the needle electrode. If merging occurs, MUPTs would affect localisation due to multi-electrode potentials from  multiple MUPs being merged. We feel this is unlikely for several reasons. Firstly, the very short duration of single-fibre action potentials means that coincidence of single-fibre potentials from different motor units is rare. Secondly, we have shown in our results that fibres belonging to different motor units frequently occupy differing positions within the sample area, with limited overlap (Fig. 5). Despite this, improvements in the MUP template matching algorithm will also be sought for the future. The fibre localisation results are consistent with our modelling work, localising multiple fibres throughout the centre of the recording area. Beyond 0.8 mm of the electrode array, far fewer muscle fibres were localised, consistent with the attenuation of single-fibre signal through tissue we have seen in our modelling. Fibres found beyond 0.8 mm likely represent muscle fibres with increased fibre diameter, with correspondingly increased amplitude transmembrane currents which were better distinguished from noise.
It is possible to calculate nearest neighbour distance distributions from single fibre positions, as described in glycogen depletion studies of cat soleus and tibialis anterior (Bodine-Fowler et al., 1990). Perhaps not surprisingly our results suggest larger median nearest neighbour distances in humans than in cats (0.54-1.21 mm compared to $0.1 mm). However, compared to whole muscle sampling of other methods (Bodine-Fowler et al., 1990), MicroEMG uses a narrow and relatively short sampling corridor size. The edge effects of this sampling region may mean that the true nearest neighbour to a fibre lies outside the detection boundary, thus overestimating the nearest neighbour distance (Floresroux and Stein, 1996).
The ability to localise fibres in 2-dimensions relative to the recording array allows a true fibre density for each motor unit to be calculated. This has previously been estimated as 2.6 fibres/ mm 2 in motor units of biceps brachii, based on post-mortem studies in neonates where an average of 163 fibres per motor unit were found (Christensen, 1959), and a motor unit territory of 9 mm diameter (Buchthal et al., 1959;Hilton-Brown and Stålberg, 1983). Our data show lower fibre densities of between 0.625-1.8 fibres/mm 2 . This lower range may be due to failure to detect all of the fibre potentials near to the electrode, or the inability of our algorithm to separate temporally and spatially adjacent fibres. Conversely, using data from neonatal muscle in which individual fibres have a lower diameter may over-estimate fibre density in adult muscles. It may be that the true value lies somewhere between these extremes.
Comparing the values of mean needle fibre density, across all recordings, we produced a value of 1.22 ± 0.29 fibres. This compares favourably to the literature, where normative data indicate a fibre density of 1.36 ± 0.11 (Hilton-Brown et al., 1985), although we showed much more variability. The motor unit territory dimensions were also within the range reported for previous multielectrode studies (Buchthal et al., 1959(Buchthal et al., , 1957Schwartz et al., 1976;E. Stålberg et al., 1976), although most of these were performed in the biceps brachii muscle. These motor unit dimensions are highly likely to be underestimates, since the fibres of the outer extent of the motor unit may lie beyond the detection range of the electrode on one or both sides, and the electrode position was not optimised to transect the entire motor unit.
It is reasonable to question the location of fibres that appear to be 'within' the dimensions of the needle. We suspect that in these cases the majority of the muscle fibre does occupy this position on either side of the array, with a short segment being deformed away from this position by the array itself. We therefore consider that these positions are a reasonable representation of the iso-electric angle over the whole length of the muscle fibre.
On detailed review of our fibre localisation algorithm, we found an important limitation. We noted apparent clustering of fibres near to the ends of the electrode array. Our modelling showed that localisation error is increased in these areas, likely due to the reduced number of electrode potentials which can be used for signal deconvolution in the search space.
In fact, fibres that are located beyond the length of the needle produce significant electrode potentials on only the closest electrodes. Normally, our algorithm always uses 5 electrode potentials to localise the fibre in space. Where fewer electrode potentials are available, the optimal search position remains biased towards the middle electrode of the 5 electrode potentials. Referring to Fig. 4A, fibres located between x = 0-0.75 mm will all result in estimated positions surrounding x = 0.75 mm. This has the net effect of clustering fibres with true position x < 0.75 mm at position x = 0.75 mm, and a similar effect will occur at the opposite end of the needle (x > 6 mm).
Our fibre localisation algorithm often failed to distinguish fibres with low separations in both space and inter-spike-interval (Fig. 4C). Hence, it is possible that adjacent fibres within the same motor unit would be interpreted as a single fibre, leading to overestimation of nearest neighbour distance and an underestimate of the fibre density.
In healthy tissue, data from glycogen depletion experiments indicate that the number of adjacencies is either the same as would be expected from random innervation (Brandstater and Lambert, 1973) or fewer adjacencies than random (Willison et al., 1980), with 12% of fibres occurring in groups (Bodine et al., 1988;Brandstater and Lambert, 1973). However, in reinnervated muscle, adjacencies occur more often (Schwartz et al., 1976), which is potentially problematic for our technique which aims to improve diagnosis of neuromuscular diseases such as amyotrophic lateral sclerosis. To mitigate this in the future, we aim to study the transmembrane currents produced by the fibre localisation algorithm. Where fibres are grouped and temporally coincident, we would expect to see significantly larger transmembrane currents owing to the adjacency of multiple fibres activating simultaneously.
There are further limitations with our technique; the images produced represent a $2 mm diameter core through multiple motor units rather than demonstrating an entire motor unit cross-section. Additionally, the modelling of fibre localisation is based upon a tissue convolution function that is homogenous and anisotropic. Heterogeneity in the tissue media may lead to localisation errors, and so this technique may be less accurate in sarcopenic or cachectic muscle where more connective tissue or fatty infiltrate can be present.
Despite these limitations, the ability to rapidly gather information on fibre density and spatial organisation in vivo during voluntary muscle contraction provides a powerful alternative to the glycogen depletion technique for the study of human motor physiology and pathology.
For the purposes of this study, we have optimised our algorithm for correctly identifying muscle fibre positioning surrounding the needle electrode. However, since our technique might enable reconstruction of the muscle fibre transmembrane current potential shape and current amplitude, this raises the exciting possibility that other muscle fibre properties may be estimated.
Our system moves electromyography from the interpretation of a highly abstract 1-dimensional motor unit potential to the creation of a 2-dimensional image, promising a step-change in both diagnostic accuracy and ease of use. Our findings in healthy human muscle also provide the first baseline data against which diseaseinduced changes in muscle fibre spatial distribution can be compared.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.