Simulation of Motor Unit Action Potential Recordings From Intramuscular Multichannel Scanning Electrodes

Multi-channel intramuscular EMG (iEMG) provides information on motor neuron behavior, muscle fiber (MF) innervation geometry and, recently, has been proposed as a means to establish a human-machine interface. Objective: to provide a reliable benchmark for computational methods applied to such recordings, we propose a simulation model for iEMG signals acquired by intramuscular multi-channel electrodes. Methods: we propose several modifications to the existing motor unit action potentials (MUAPs) simulation methods, such as farthest point sampling (FPS) for the distribution of motor unit territory centers in the muscle cross-section, accurate fiber-neuron assignment algorithm, modeling of motor neuron action potential propagation delay, and a model of multi-channel scanning electrode. Results: we provide representative applications of this model to the estimation of motor unit territories and the iEMG decomposition evaluation. Also, we extend it to a full multi-channel iEMG simulator using classic linear EMG modeling. Conclusions: altogether, the proposed models provide accurate MUAPs across the entire motor unit territories and for various electrode configurations. Significance: they can be used for the development and evaluation of mathematical methods for multi-channel iEMG processing and analysis.


Main
Acronyms  I TRAMUSCULAR EMG (iEMG) modelling supports the interpretation of the iEMG signal generation in human muscles. It allows to vary the parameters of both the motor neuron (MN) pool and the muscle fibers (MFs) in order to test and validate iEMG-based computational methods, such as motor unit (MU) territory estimation and iEMG decomposition.
Different applications require simulation models of different complexities. For example, iEMG decomposition algorithms are often tested using signals simulated with phenomenological EMG models [1]. These approaches involve the convolution of experimental or simulated MN spike trains with experimental motor unit action potentials (MUAPs) and provide known spike trains and adjustable levels of additive noise. However, they do not model the neuromuscular jitter, morphological variability of MUAPs, noise generated by distant MUs, and the electrode geometry.
Compared to the phenomenological approach, biophysical modeling of iEMG includes the calculation of each single fiber action potential (SFAP) as a function of the MF's morphology and of the electrode's position. These approaches provide an infinitely wide dictionary of MUAPs and take into consideration the terminal arborization geometry, neuromuscular jitter, as well as shape and position of the electrode. This flexibility is This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ important in such modelling applications as the analysis of MUAP morphology [2], simulation of neuropathies and myopathies [3], [4], and assessment of accuracy of EMG decomposition algorithms [5], [6].
Existing iEMG simulation methods [7]- [10] focus on singlechannel intramuscular recordings. However, in clinical and basic physiological applications, there is a growing interest in multi-channel electrode systems [11]- [14]. In this paper, we present several methods that allow the simulation of multichannel recordings, which, to the best of our knowledge, have not yet been addressed in other studies. Previously proposed models that are the closest to our approach are the multi-channel surface EMG simulation methods [15], which, however, mainly address the description of the volume conductor properties.
Multi-channel iEMG modeling requires a number of specific characteristics that cannot be obtained by multiple application of an existing single-channel approach. In particular, multi-channel recordings may involve simultaneous observations from different parts of the muscle, thus the geometry of MU territories should be consistent in the whole muscle cross-section. For example, the models described in [7] and [8] provide the simulation of a single MUAP, not considering the mutual arrangement of MU territories. This problem was later addressed in [9] by using a geometric model of MU territory locations in the muscle cross-section. However, this previous approach assigns an exact number of MFs to each MN, which forces the innervation of MFs at large distances from the imposed territory. While this issue minimally influences the single-channel, localized intramuscular recordings, it has a substantial effect on the multi-channel ones.
Also, multi-channel recordings are influenced by the endplate geometry to a larger extent than the single-channel ones. The unimodal Gaussian and uniform distributions used to generate the NMJ locations in previous approaches [9], [16] produce an excessive scattering of single fiber action potentials (SFAPs) propagation delays. This makes the simulated MUAP waveforms more complex and with greater changes across the channels than in experimental observations [2], [17]. It was suggested in [2] that the axial distribution of the neuromuscular junctions (NMJs) should be instead composed of smaller sub-distributions corresponding to the separate branches of the axon.
In this paper, we propose several novel techniques that altogether provide a simulation of MUAPs recorded by intramuscular multi-channel electrodes. These approaches can also be used in combination with previous ones [7], [9], [15], [18] to pursue different simulation goals. More specifically, the new elements we propose include: 1) a new way to generate even distributions of MF coordinates and MU territory centers; 2) a controllable and accurate method for MF-MN assignment; 3) an improved model of MN axon terminal arborization; 4) multi-channel intramuscular electrode modeling; 5) extension to shifting or scanning intramuscular electrode modeling; Although the proposed model focuses on the simulation of MUAPs, it can be easily extended to a full iEMG generation. This requires an appropriate MN pool model for the generation of spike trains, such as that presented in [16].
The remainder of this article is organized as follows. MFs and territory centers distributions in the muscle cross-section are presented in Section II-A. Next, the fiber-neuron (MF-MN) assignment procedure is described in Section II-A3. Section II-B4 explains the modeling of terminal arborizations and of MN action potential (MNAP) propagation delays. Finally, the simulation of MUAPs recorded by multi-channel and scanning electrodes is described in Section II-C6. Examples that demonstrate the performance of the proposed sub-models are presented both in sections "Methods" (II) and "Results" (III).

1) Distribution of MU Territory Centers in the Muscle
Cross-Section: As previously described, [9], [18], [19], we approximate the muscle by a cylinder where all MFs are parallel to the z-axis and where the xy plane constitutes the muscle crosssection. In order to assign MFs to MNs and simulate MUAPs, the territories of the MUs in the muscle cross-section should be first defined. A MU territory is a circular region in the muscle cross-section determined by two parameters: the xy coordinates of its center and its area. Anatomically, it is the area that includes all the MFs innervated by the MN. In this section, we describe a new method to generate the territory centers, while the areas will be addressed later in Section II-A4.
In the existing models [16], [18], the coordinates of MN territory centers are assumed to evenly fill the muscle crosssection. For this purpose, the territory centers are drawn from a uniform distribution over the muscle cross-section. However, the uniform distribution provides very uneven arrangement of the points (see Figure 2 for illustration). In turn, this leads to uneven and physiologically incorrect densities of the innervated MFs, as shown in [20]. Another simulation model [9] generates centers as randomly altered nodes of a rectangular grid, which partly resolves this problem, but may be inconvenient to implement due to the fact that the generated centers may end up outside of the admissible region.
We propose an alternative approach that consists of using the farthest point sampling (FPS) [21], [22]. FPS is a family of algorithms that fill an enclosed 2D domain by iteratively adding points that are maximally distant from the previously added ones. This property of the FPS algorithms allows to maximally disperse the MU territory centers in the muscle cross-section, thus achieving an even scattering. At the same time, the generated points are quasi-random (e.g., their arrangement solely depends on the form of the region and on the position of the initial point) and do not follow any regular pattern, such as the rectangular grid. A similar approach was proposed in [20], where the Mitchell's Best Candidate method was applied [23].
FPS algorithms provide an additional feature that may be of interest in iEMG simulation. They place each new territory center as far as possible from the previously generated ones. The subsequent assignment of the territory centers to the MUs can be performed either randomly or in a specific order. If it is performed in the order of decreasing MN sizes, the centers of the MUs with similar sizes will be as distant from each other as possible, thus achieving an even distribution of MUs of different sizes in the cross-section. An illustration is provided in Figure 1, where the circles around the territory centers indicate the relative sizes of corresponding MUs. In this example, the muscle crosssection is evenly filled with territory centers of MUs of different sizes. This property is relevant when an even distribution of size-dependent parameters in the muscle cross section is of interest. Possible examples are the distributions of MFs of types I and II [24] (types I and II are predominant in, respectively, smaller and larger MUs [25]), as well as MF diameters [26], [27] (MFs of type I tend to have smaller diameters than those of type II [28]). Even distribution of MUs of different sizes will contribute to an even distribution of such parameters in the muscle cross section. Alternatively, MNs can be assigned to the territory centers in a randomized order if the aforementioned structure is not required.

2) Distribution of Muscle Fibers in the Cross-Section of the Muscle:
The distribution of MFs in the muscle crosssection (xy-plane) should be uniform and, preferably, should take into consideration the diameters of MFs. In previous works [16], [18], the MFs' locations were drawn from the uniform distribution within the corresponding MU territory. As it was noted above, this method, combined with previous approaches to the territory centers generation, does not guarantee the global uniformity of MF density in the entire muscle crosssection. Alternatively, in [9] MFs were positioned prior to the generation of the territory centers in the nodes of a regular rectangular grid.
We propose to generate the MFs locations using the aforementioned FPS algorithms family. In our simulation model, we generate them independently from the territory centers, and simulate the MF-MN assignment in following steps, similarly to [9]. Figure 2 provides a comparison of the locations generated by this method with those drawn from the constant distribution. In our simulation we use constant density of 400 MFs per mm 2 . However, it is worth mentioning that the FPS also allows a variable local density of MFs. for each MN n do 5: calculate P a n , P g n , P d n ; 6: w n ← P a n · P g n · P d n ; 7: end for 8: assign MF f to a random MN n with weight w n ; 9: end while : We establish a randomized procedure in which the generated MFs are assigned to MNs (i.e., innervated by this MN) according to the following parameters: the expected number of MFs innervated by the MN, MF's proximity to the territory center and presence of adjacent MFs already assigned to that MN. For each MF-MN pair, the probability of assignment is represented by a score that combines the influence of each of these factors:

3) Fiber-Neuron Assignment
where r P f (n) is the score that characterizes the probability of the f -th MF to be assigned to the n-th MN, given the size of the MN, positions of the MF and MN's territory center, and neighboring MFs innervation; r P a n denotes the a priori probability of assignment, i.e., solely given the size of the MN; r P g n (x f , y f ) denotes the probability of MF with coordinates (x f , y f ) to be innervated by n-th MN, given its territory center location; r P d n (n, f, n c ) is an indicator function returning 0 if any of n c closest neighbors of f -th MF is already assigned to n-th unit, and 1 otherwise. The pseudocode of the MFs assignment procedure is presented in Algorithm 1. Let us consider each multiplier in (1) in a detailed way.

4) A Priori Probability of Assignment:
A priori probability of an MF being innervated by n-th MN is proportional to the total number of MFs that this MN should innervate. We suppose that this number is itself proportional to the size of the MN. Thus, we calculate the a priori probability P a n of MF innervation in the following way: where s n is the size of the n-th MN which can be modeled using an exponential distribution for recruitment thresholds, as proposed in [16].

5) Distribution of MF Coordinates Around a Territory
Center : We assume that the innervation territories of MNs are circular; this assumption is common to most EMG simulation models and is supported by experimental data [29]. We model the P g n as a symmetrical two-dimensional Gaussian distribution: (3) where x and y are the coordinates of the MF; mean μ g i is coincident with the territory center of n-th MN; S n is an out-of-border coefficient (see explanation below in this section); standard deviation σ g n = a n /πC is proportional to the MN's innervation area a n with scattering coefficient C.
We assume that the innervation areas of MNs are proportional to their sizes with a scaling factor A/k: a n = s n /s N · A/k, where A is the area of the muscle cross-section. The value of k sets up the area of the largest MN as a fraction of the muscle cross-section area: k = A/a N . The value of k varies across muscles. In our simulation, we have chosen a value of k = 4. The value of the scattering coefficient C regulates the tightness of the Gaussian distribution of MFs around the territory center. We calculate it assuming that a n is the area of 0.99 confidence circle for the corresponding distribution, giving us C = inv-χ 2 (0.99, 2) = 9.21.
Innervation areas of some MNs may partly lie outside of the muscle border. Coefficient S n takes this fact into consideration and normalizes the corresponding distribution (3) in order to ensure that the number of innervated MFs will still be correct. S n is calculated as a double integral of the original Gaussian distribution in (3) above the domain corresponding to the muscle cross-section. Thus, the probability P g n sums to 1 while integrating over the muscle region. Considering the classification proposed in [26], this approach can be assigned to uniformaugmented territory placement.

6) Adjacency of Muscle Fibers Innervated by the Same
Motor Neuron: Due to the phenomenon of self-avoidance in the arborizations of MNs axons [30], MFs of the same motor unit rarely lie next to each other. This fact is reflected in equation (1) by using an additional factor P d n (n, f, n c ), which equals to zero if at least one of the n c closest MFs is already innervated by n-th MN.
The value of n c also regulates the scattering of MNs' MFs across the muscle cross-section. We suggest n c = 5 for a regular   Figure 3. This distribution follows very closely the one imposed by the model F n = F · s n / N n =1 s n . The resulting innervation areas also lie close to their modelimposed values. They can be calculated both as areas of convex hulls or areas of 0.99 confidence ellipsoids. An example of resulting territories' forms is provided in Figure 4.

B. Model of the Axon Branching
A neuromuscular junction (NMJ) is a biological interface between a single MF and its innervating motor neuron axon. In order to innervate all its MFs, an axon splits into smaller branches, forming a complex and uneven tree structure with neuromuscular junctions at its leaves [30].
The motor neuron action potential (MNAPs) originates in the soma of the MN and propagates along the MN axon branches until reaching each of the innervated MFs. The lengths of the paths to each MF vary due to the scattering of the neuromuscular junctions in the muscle. This causes the scattering of MNAP propagation delays, which affects the morphology of the MUAPs. As concluded in [29], temporal dispersion of the MUAPs is to a larger extent due to the spatial dispersion of the NMJs than to the differences in conduction velocities of the MFs.
In this section, we will show how our simulation model calculates the coordinates of neuromuscular junctions and MNAP propagation delays. 1) Structure of the Axon Branching Model: An MN axon can be represented as a root of a tree structure that splits into several branches of smaller radii. This process is then repeated several times within each branch until each MF is reached.
In our model, we suppose that the split is done only twice (see Figure 5). Thus, each MF and its NMJ are assigned not only to a motor neuron but to a specific branch of its axon. We establish such a model in order to constrain the complexity of MUAPs while providing physiologically correct distributions of NMJs along the muscles.
2) Fiber-Branch Assignment : We assume that the number of branches is proportional to the MU's size and that its value defines the number of phases in its action potential. The following expression provides the numbers of branches/phases that correspond to the experimental action potentials for small motor units (1-2 phases) as well as for the largest ones (4-6 phases): where s n is size of n-th motor unit and · stands for rounding to the nearest integer.
For each motor unit, in order to assign each MF to a specific branch, we first define the number of branches B n using equation (4) and then run the k-means clustering algorithm over the MU's MFs coordinates in the cross-sectional plane, looking for B n clusters. In this case, k-means seeks for B n groups of closely-located MFs of the MU. Then the MFs of each group are assigned to a single axon branch.

3) Coordinates of Neuromuscular Junctions:
In the muscle cross-section (xy-plane), the NMJs coincide with their MFs. Along the z-axis, the NMJs are scattered around a certain center point, within a range that varies depending on the muscle [2], [9]. Previous approaches to the modelling of this distribution used the uniform or Gaussian densities [9], which produced correct results in single-channel simulations. However, in the case of multi-channel or scanning recordings, these approaches generate MUAPs that vary too fast between the neighboring channels or subsequent electrode positions, compared to the experimental results [2], [17]. To resolve this problem, we will implement the innervation structure of suggested in [2]. We suppose that the axial positions of the NMJs in a MU should be distributed within narrow sub-bands with different means, each associated to a specific axon branch. Our assumption is illustrated in Figure 6: each axon branch innervates a certain cluster of MFs within the muscle cross-section (left) while axial locations of their NMJs are distributed in different sub-bands (right). We model these sub-bands as Gaussian clusters with mean values scattered across the z-axis of the muscle, and standard deviations much smaller than the entire range of NMJ distribution.
We calculate mean values μ b n and standard deviations σ b n of intra-cluster densities using the following model: where parameters a μ , b μ , a σ , b σ in combination with MF and axon conduction velocities define the dispersion of the MNAP propagation delays, and, thus, the duration of MUAPs. In order to obtain an initial estimate of these parameters, we impose the largest and the smallest MUs to have MUAPs with durations of the main spike [31], [32] of 2.5 ms and 7.5 ms respectively. Considering the mean conduction velocities of their MFs to be 2.5 · 10 3 mm/s and 5 · 10 3 mm/s [7], we can approximately calculate the necessary span of their neuromuscular junctions, giving correspondingly l m in = 6.25 mm and l m ax = 37.5 mm. We also assume that the standard deviation of the cluster centers σ n is larger than the intra-cluster deviation σ b n since a MUAP usually contains several distinct phases. In simulation, we have found that it is convenient to set d σ = σ n /σ b n to 4. Finally, we note that the span of neuromuscular junctions along the z-axis for n-th motor unit can be roughly calculated as 3(σ n + σ b n ). These considerations give us the following system of equations: the solution of which for d σ = 4 gives a σ = 0.4, b σ = 2.1, a μ = 1.7, b μ = 8.3 (all in millimeters). We should consider these values as the upper estimates, since the MUAP lengths are also influenced by the MF conduction velocities and the MNAP propagation delays (see Section II-B4 for details) which were not yet taken into account. According to our observations and for MNAP propagation delays listed in Section II-B4, a set a σ = 0.25, b σ = 1, a μ = 1, b μ = 2.5 produces MUAPs with physiologically correct forms and duration.

4) Delay of MNAP Propagation:
Once all the MFs of a motor unit are assigned to their branches and the z-coordinates of NMJs are generated, we can calculate the delays of MNAP propagation towards each junction. We divide the lengths of each segment of the axon (see Figure 5) by their propagation velocities, thus, the delay for f -th MF assigned to b-th branch of n-th motor unit is: where r x j f are the coordinates of the neuromuscular junction of f -th MF in k-th branch of n-th motor unit. r x b k are the coordinates of the terminal arborization root at the k-th branch of the n-th motor neuron's axon, calculated as the mean of the NMJ locations in 3D space: n are coordinates of the first branching point of n-th motor unit, which is calculated as the mean of arborization roots: r v t is the MNAP propagation velocity in a terminal arborization of the motor neuron axon; We assume that the branches' propagation velocities are much smaller than that of the axon due to their smaller diameter and the absence of myelination in case of terminal arborization. The values that we used in our model are: v b = 10 m/s, v t = 1 m/s (for comparison, typical propagation velocity of a myelinated MN axon is about 50 m/s). To the best of our knowledge, there is yet no experimental data on v t and v b in the literature.

1) Muscle Fiber Action Potential Modeling:
In order to simulate single fiber action potentials (SFAPs) we define a potential induced in an observation point by an elementary current source located at a narrow slice of the MF at coordinate z [7], [33]: r σ r and σ z are the radial and axial conductivities of the muscle tissue (0.063 S/m and 0.33 S/m respectively [7]); During MF contraction, the transmembrane current sources are continuously distributed along the MF rather than being concentrated in certain points. Let us denote this distribution, at time t, as I(z, t). The potential φ p (t) at the observation point p, generated by this distribution can be calculated as the convolution [7]: (9) In order to model the transmembrane current distribution I(z, t), we use the approach proposed in [15] (see their equation (17)), which was also recently implemented in another EMG simulation model [18]. The intracellular action potential model from [7] was used, as suggested in [34]. The MF diameters and conduction velocities were calculated using models presented in [9] and [7].
2) Motor Unit Action Potential Modeling: Action potential of a MU is modeled as a linear sum of its MFs' SFAPs [9]. Taking in consideration the delays and neuromuscular jitter in axon branches (see also Figure 7 for illustration): where Φ np (t) is the MUAP of n-th motor unit, observed in point p; F n is the number of MFs innervated by this motor unit; φ f p (t) is the SFAP of the f -th MF of n-th motor unit, observed in point p; d f is the delay between the MNAP discharge and its arrival to the NMJ of the MF (see expression (7)); ζ is a neuromuscular delay caused by the jitter, randomly drawn for each MF and for each new realization of Φ ip . Delays caused by the neuromuscular jitter are of the order of μ s [35]. To be able to simulate their influence on the SFAPs, while using reasonable sampling frequencies, we use sub-sample waveform shifting presented in [36].

3) EMG in a Single Observation Point:
EMG is modeled as a linear sum of contributions from all MUs [37], while the contributions are the spike trains convoluted with the MUAPs. We formulate the expression for simulated EMG, acquired in observation point p, in a similar way: where y p (t) is the simulated EMG signal in observation point p, U n is a vector of spikes' time instants for the n-th motor neuron; and k is the index of a spike in U n .

4) EMG in a Single-Channel Electrode:
Due to the fact that a metallic electrode is a conductor, the electric potential is constant across its volume. Its value can be approximated by an average potential across a surface that coincides with the position of the electrode. Therefore, the SFAP detected by the electrode can be calculated as an integral of the potential φ(t) over the electrode's surface. Due to the linearity of (10) and (11), the same applies to MUAPs and the overall signal detected by the electrode.
In our model, we approximate the recording surface by a number of elements with an observation point associated to the center of each element. Therefore, the electrode potential is equal to the sum of element potentials weighted by their areas. Thus, iEMG signal recorded by an electrode is a linear combination of the signals in the observation points: where e p is the area of the electrode's element associated to the p-th observation point, its sign depends on the polarity of the corresponding amplifier input; y p (t) is the EMG signal calculated at observation point p using equation (11); in the following, E and Y(t) will be referred to as electrode matrix and observation vector. See examples in the appendix. MUAP recorded by an electrode modeled by matrix E can be calculated using (11) and (12) while setting U n = {0}:

5) EMG in a Multi-Channel Electrode:
A multi-channel EMG signal can be calculated using a stack of electrode matrices each corresponding to one of the M channels:

6) EMG in a Shifting
Electrode: Shifts can be modeled as a combination of translations and rotations of the electrode along a specified trajectory in the muscle. This trajectory can be approximated by a number of nodes D linked by successive rigid transformations.
The current position of the electrode on the trajectory curve can be specified by a continuous path parameter 0 ≤ λ ≤ D − 1, where D denotes the overall number of nodes in the trajectory.
The signal, as a function of the current position of the electrode, can be calculated as follows: where I d = I · δ(λ − (d − 1)), δ(·) is the Dirac delta function and I is an identity matrix of size P ; Y d (t) is the observation vector in the trajectory node d.
The signal acquired in positions located between the trajectory nodes can be linearly interpolated, given a sufficiently fine trajectory sampling. One can express the signal acquired in a specific position λ on the trajectory as: whereÎ d (λ) is a weighted identity matrix determined as follows: Kernel (17) is equal to one when λ = d − 1 (i.e. when calculating the signal exactly in trajectory node d) and linearly weights the neighbouring nodes d and d + 1 when d − 1 < λ < d. An example of a MUAP captured by a fine-wire electrode, that moved transversally to the MFs, is shown in Figure 8.
Path parameter λ can be a function of force or time since usually electrode shifts occur due to either muscle deformation during contraction or other factors that can be described as functions of time.

III. SIMULATION RESULTS
In order to demonstrate the capabilities of the proposed model, we provide several application examples. Figure 8 shows MUAPs detected by an array of 16 equidistant point electrodes spaced by 1 mm gap and placed at the angle of 30 degrees to the MFs in the plane that passes through the central axis of the muscle. In total, 15 channels were obtained by the consecutive differentiation of the signals. We have chosen an MN for which the center of the territory lied close to the midpoint of the electrode array. The length of the muscle in this simulation was 50 mm, and the end-plate zone was positioned around the center of the muscle, as stated in formula (5). The array was placed in one of the halves of the muscle and didn't cross the MU's end-plate zone.

A. Multi-Channel MUAP Example
From Figure 8, we notice several relevant results. First, the amplitude of the MUAPs in each channels is inversely proportional to the distance between the channel's electrodes points and the center of the MU's territory. Second, the centers of energy of the MUAPs shift to the right as the distance between the channel and the end-plate increases, due to the simulation of SFAP conduction. Third, the transformations between MUAPs in the neighboring channels are consistent and the MUAPs have physiologically correct shapes and durations (approximately 5 ms).

B. Single MU Territory Scanning Simulation
Using equation (16), we can simulate a "scan" of a MU territory. An example is presented in Figure 9 where a MUAP is simulated at 10 equidistant nodes positioned along a straight trajectory passing through the center of a MU's territory. The overall duration of the generated MUAP is approximately 5 ms. This result is qualitatively in agreement with experimental observations [2], [17].

C. MUs Territory Assessment
Using multi-channel MUAP modeling, it is possible to simulate experimental studies that aim at the assessment of MU innervation territories. As an example, we have simulated a procedure similar to the experimental study described in [38]. For this purpose, we have simulated a 10 mm long array of 11 equally spaced intramuscular electrodes, thus comprising 10 differential channels with consecutive differentiation, inserted into a 10 mm wide muscle at an angle of 90 • with respect to the MFs. The diameters of the MU territories were estimated from the simulated recordings as the length of the array comprising electrodes in which the MUAP's peak was greater than four times the standard deviation of the baseline noise [12], [38] (the SNR was set to 15 dB). Only MUs with recruitment thresholds that are below of 50% MVC were considered for this analysis. In order to establish this limit, the contraction force was calculated using the model proposed in [16].
As shown in Figure 10, the estimated diameters generally correlated with their true values. However, for most of the MUs, the electrode array did not cross the center of the territory, which resulted in an underestimation of the diameters (e.g., see MUs 40, 43, 55 in Figure 10). Moreover, this technique limits the estimated diameters to be multiples of the inter-electrode distance and therefore has low resolution. The proposed simulation model can be used to test and evaluate more complex algorithms for the MU territory assessment.

D. Multi-Channel iEMG Decomposition
Finally, we present the results of the application of a decomposition algorithm to the simulated signals. In order to generate an iEMG signal for decomposition, we have implemented the motor neuron pool model proposed in [16].
A linear array of five electrodes (1-mm interelectrode distance) providing four differential recordings, was simulated. We measured the average power of the simulated signal in all channels at maximal net excitation in order to obtain a reference value for the calculation of the standard deviation σ of the additive noise. MUAPs whose maximal absolute value exceeded 4σ in at least one of the four channels were considered detectable and the corresponding MUs were included to the annotation of the signal. The SNR in all channels was set to 15 dB and a trapezoidal contraction reaching 20% MVC was generated.
The simulated iEMG signal was decomposed by MTL, the multi-channel version of the algorithm proposed in [39], [40]. The decomposition was evaluated using the classification phase sensitivity and positive predictivity [41] averaged across 11 detected MUs, which resulted in, respectively, 0.93 ± 0.04 and 0.97 ± 0.03. Therefore, the simulated MUAPs fit the assumptions that the used decomposition algorithm makes about experimental recordings.

IV. DISCUSSION
The presented approaches to the modeling of electrode and innervation geometries can be applied independently and/or in combination with previous methods. For example, the proposed FPS-based method to generate the locations of MFs and MU territory centers can replace or be replaced by other methods [9], [18]. The proposed fiber-neuron assignment procedure can be used within other simulation models that assume an even distribution of MU territories [9], [20]. The NMJ scattering and axon branching model can also be used independently, replacing common uni-modal Gaussian distribution models in other simulation strategies. Finally, the proposed electrode model can be used in combination with other innervation geometry models.
While the presented simulation model has several parameters to tune, these parameters are designed to reflect the physiology of the motor system and thus can be selected according to known physiological variables. We suggest that the proposed approaches, being simple and physiology-based, provide greater precision and flexibility than previous models, specifically for multichannel and scanning electrodes.

V. CONCLUSION
We have described a new model for MUAP simulation that includes the possibility to simulate multi-channel intramuscular electrodes, their arbitrary positioning, and gradual shifting during acquisition. The model also includes new methods for establishing uniform distributions of MFs and territory centers in the muscle cross-section, for tuning the fiber-neuron assignment, and for controlling the complexity of the MUAPs.
This model can be used to simulate a wide range of experimental studies and computational methods, such as for MU territory estimation, conduction velocity measures, denervation, and reinnervation. A full iEMG simulation has also been obtained by convolution of the modeled MUAPs with spike trains provided by a previous model of the motor neuron pool behaviour. This extension can be used for the assessment of decomposition algorithms in terms of their robustness towards the variations in MUAP characteristics.
In conclusion, we proposed a new model for the simulation of intramuscular multi-channel and scanning EMG recordings that has a wide range of potential applications in the test of computational methods applied to intramuscular EMG signals. The corresponding code can be accessed in the authors' online repository [42].

APPENDIX A EXAMPLES OF THE ELECTRODE MATRIX FOR BASIC SIMULATION CASES
A fine wire electrode can be approximated by a pair of points with equal areas. In the case of bipolar acquisition, the resulting signal is equal to the difference between potentials observed in the two points (see expression (12)): A signal acquired by an array of point electrodes with consecutive differentiation can be represented as follows (see expression (14) Assuming that the electrode's trajectory is approximated by only two nodes, an EMG signal from a fine-wire electrode, before the shift (λ = 0) and after the shift (λ = 1) can be expressed as (see equation (15)): where the lower index of y corresponds, as previously, to the electrode element, while its upper index denotes the trajectory node. Applying expression (16) to the previous example, we can calculate the signal acquired at 1/4-th of the way (λ = 0.25):