EMG Modeling

1. Introduction

The aim of this chapter is to describe the approaches used for modelling electromyographic (EMG) signals as well as the principles of electrical conduction within the muscle. Sections are organized into a progressive, step-by-step EMG modeling of structures of increasing complexity. First, the basis of the electrical conduction that allows for the propagation of the EMG signals within the muscle is presented. Second, the models used for describing the electrical activity generated by a single fibre described. The third section is devoted to modeling the organization of the motor unit and the generation of motor unit potentials. Based on models of the architectural organization of motor units and their activation and firing mechanisms, the last section focuses on modeling the electrical activity of a complete muscle as recorded at the surface.

A mathematical model of a system describes the relations between a number of physical variables involved in the system. A mathematical model is a set of equations that can be implemented on a computer to study and to simulate the behaviour (response) of the system under specific conditions. EMG models presented in this chapter are structure based or structural, which means that they describe elements of the real biological structure and characterize them in a reductional way in order to represent the system’s elements, behaviours or mechanism that are of importance. In the EMG models outlined here, the input variables or parameters are those that describe the anatomical, physiological, and functional properties of the biological structure under study (single fibre, motor unit, or entire muscle), whereas the output parameters are typically the extracellular generated potentials and/or specific quantitative measurements of these potentials.

Models of EMG activity are useful to address the “forward problem”, that is, how specific mechanism and phenomena influence the generated potentials, as well as the “inverse problem”, that is, how the extracellular potentials provide information about the underlying mechanism and phenomena. Accordingly, a desirable feature of an EMG model is that it allows studying the effect of the model's (input) parameters on the waveform of the potential, providing insight into the relationships between the anatomical and/or physiological properties of the fibre and the shape of the potential.

2. Modeling electrical conduction in skeletal muscle

Striated muscle is composed of a large number of striated muscle cells, also called muscle fibers. These elongated, cylindrical cells are arranged parallel to one another, and each one is surrounded by a plasma membrane called the sarcolemma. Muscle contraction is created via the repeated activation of several groups of muscle fibers, each of which is governed by a single motorneuron through its axon (Lieber, 2010). Figure 1(a) shows a portion of a muscle fiber that is attached, at the neuromuscular junction, to the terminal branch of its axon.

2.1. Depolarization and repolarization of a muscle fiber membrane

A muscle cell (fiber) is activated by electrical impulses coming from the motorneuron, which brings about the fiber's depolarization and the generation of a transmembrane voltage (normally defined as the extracellular electrical potential minus the intracellular one). Under normal conditions, the extracellular potential is practically zero, and so the transmembrane voltage can be considered to be practically the same as the intracellular action potential (IAP) (Burke, 1981; Plonsey and Barr, 2000). Membrane depolarization starts at the neuromuscular junction and extends along the muscle fiber towards both ends of the cell. As a result, two IAPs propagate without attenuation and with constant velocity v along the muscle fiber towards the tendons, where they extinguish (Lieber, 2010).

The plasma membrane of a muscle fiber has the well-studied property of actively maintaining a nearly constant potential difference between the intracellular and extracellular environment. This voltage is normally referred to as the resting potential and has a value of about -80 mV [Fig. 1(b)]. The negative sign indicates that the interior is more negative than the exterior (negative polarization). In Fig. 1(a) the polarization of the fiber membrane is represented by a number of layers of negative signs. When the fiber is at rest, the number of negative layers remains unchanged. After fiber activation, two potential profiles (one at each side of the neuromuscular junction) arise along the fiber membrane. These profiles are normally referred to as waves of excitation, excitation sources or simply IAPs [Fig. 1(a)]. In Fig. 1(a) it can be seen that, along the spatial profile of the IAP, the number of negative-signed layers changes progressively with axial distance (the x-axis in Fig. 1). This reflects the fact that the IAP profile has gradual depolarization and repolarization transitions, as shown in Fig. 1(b).

2.2. The electric volume conduction

Knowledge of muscle fiber excitability and IAP propagation would be clinically useless if the electromyographist had to penetrate muscle fibers to obtain information about membrane processes. Moreover, recording of the IAP in situ from human muscle fibers has not yet proved feasible in EMG studies (Ludin, 1973). All electromyography (EMG) techniques are based on the fact that local electrophysiological processes result in a detectable flow of the transmembrane current at a certain distance from the active sources (i.e., muscle fibers). This flow of current in the tissue (i.e., the volume conduction), allows EMG measurements to be made at a distance from the sources.

The so-called Principle of Volume Conduction can be considered as a three-dimensional version of Ohm’s law, which establishes that an electric current I, flowing between two points connected through a resistance R, generates a potential difference V between these points: V = I R. In the case of living tissue, the electrical impedance is the inverse of the electrical conductivity σ. So, the potential recorded at a point P₀ (x₀, y₀, z₀) within an infinite volume with uniform conductivity σ_i produced by a current I_s injected in the same volume at a point P (x, y, z) can be calculated as

$$V_{P_0}=\frac{1}{4\pi {\sigma }_i}\frac{I{}_s}{r_i}$$E1

where r_i is the shortest distance between the points P₀ and P.

From inspection of (1), two main conclusions can be drawn. First, the potential recorded at a certain point is proportional to the strength of the current source, a feature highly desirable for electrodiagnostic medicine. Second, both r_i and σ_i are in the denominator of the equation (1). Thus, assuming a constant transmembrane current, the potential decreases with increasing radial distance and with increasing conductivity.

Figure 1.

a) Schematic representation of a portion of muscle fiber in which two excitation sources [IAP(x)] are propagating with velocity v from the neuromuscular junction (NMJ) to the fiber ends. The polarization of the fiber membrane is represented by several layers of negative signs. The number of negative-signed layers within the fiber region delimited by the intracellular action potential (IAP) changes gradually with axial distance, but it is constant within the regions where the fiber is at rest. The transmembrane ionic electric current, I_m(z), is also indicated. (b) Spatial profile of the IAP with its depolarization and repolarization phases. L and T_in are the spatial extension and temporal duration of the IAP, respectively.

The principle of volume conduction is valid only as an intuitive approach to an understanding of the generation of an extracellular potential within a muscle. The simplicity of equation (1) hides important aspects that need to be clarified. First, the bioelectrical source cannot be described as a single injected current at a certain point, but it is rather a compound of multiple sources (see Section 3.1.2). Second, muscle fibers are of finite length, which implies that the assumption of an infinitely large volume conductor is never satisfied in practice. This will have important consequences: it will give rise to non-propagating components (see Section 2.2). Third, as muscle fibers can often be considered parallel to each other, conductivity of the muscle tissue in the longitudinal direction (σ_x) is higher than in the transversal (σ_r), i.e., the volume conductor is anisotropic with anisotropy ratio K_an = σ_x/σ_r.

3. Modeling the electrical activity of the single muscle fiber

In the last half of the 20th century EMG studies have directed attention to the calculation of single fiber action potentials (SFAPs) produced by excitable fibers and especially by fibers of finite length. The development of SFAP models was possible only after the principles of volume conductivity had been determined and the modelling of bioelectrical sources was well established.

3.1. Modeling the bioelectrical sources of muscle fibers

3.1.1. The principles of bioelectricity

The principle of the electromotive surface proposed by Helmholtz in 1853 provided the basis for the electrical potential theory of volume conductors. In addition, he introduced the concept of an electrical double layer into the theory of electricity and suggested its use for the solution of certain boundary problems in electrical potential theory. The concept of the double layer (or dipole) source, however, went unused for about 30 years until Wilson et al. (1933) demonstrated its appropriateness for modelling the excitation source of a single circular cylindrical fiber. After Wilson, it was not until 1974 that Plonsey definitively clarified the underlying electrostatic principles of the single and double layer sources and related them to the concepts of monopole and dipole, respectively.

Traditionally, two kinds of sources (generators) have been considered in the literature (Lorente de No, 1947; Plonsey, 1974), namely, the “monopole” and the “dipole”. A monopole is a single (point) source or sink of current within a conducting medium. It is quite rare that problems in bioelectricity involve monopoles, since all bioelectric generators involve at least source and sink combinations (Plonsey, 1974).

3.1.2. A distributed presentation of the excitation source (excitation function)

Because of its propagation along the fiber, an IAP does not merely exist as a function of time: it also spreads out along the fiber as a function of space. The length of the IAP profile along the fiber, L, is defined by the product of the IAP duration, T_in, and the propagation velocity v [Fig. 1(b)]. In fact, the formation of an electrical field around a fiber depends on the spatial extension of the IAP along the fiber (Dimitrova and Dimitrov, 2006; Rodriguez et al., 2011). In addition, the spatial profile of the IAP is smooth [Fig. 1(a)], implying that the electric properties of the fiber membrane affected by the IAP change gradually with axial distance. Accordingly, a correct presentation of the excitation function consists of a sequence of cylinders (fiber portions) of equally-infinitesimal length dx, each cylinder containing a density of sources, as represented in Figs. 2(b) and (c). If these sources are considered dipoles, then each of these cylinders should be represented by a double-layer disk [Fig. 2(b)], whereas if the sources are regarded as monopoles, then cylinders should be modeled as single-layer disks.

Figure 2.

a) Representation of the IAP spatial profile and its first spatial derivative, ∂IAP/∂x. Schematic representations of the IAP as a sequence of double layer disks (b) (each disk comprising a density of dipoles), and as a sequence of lumped (point) dipoles (c) lying along the axis of the fiber.

From the above it follows that the calculation of the extracellular potential generated by a single excited fiber, Φ_e, can be reduced to the sum of the potentials produced by a sequence of double (or single) layer disks distributed along the IAP spatial course (Dimitrova and Dimitrov, 2006). The specific mathematical derivation by which the extracellular potential is expressed in terms of double layer disks is presented below.

3.1.3. Calculation of the extracellular potential on the basis of dipoles

The dipole-based presentation of the source was first introduced by Wilson et al. (1933) and subsequently developed by Plonsey (1974). It is based on the hypothesis that the variation of the membrane electrical potential across an infinitesimal portion of the fiber membrane produces, in the extracellular medium, an electrical field that can be assumed to be equivalent to that produced by a lumped dipole (Wilson et al., 1933). So, the potentials produced by a double layer disk and a point dipole whose moment is proportional to the disk area are almost identical. This provides the basis for representing the portion of the fiber affected by the IAP as a sequence of dipoles distributed equidistantly along the IAP spatial profile (Dimitrova and Dimitrov, 2006; Rodriguez et al., 2011), as shown in Fig. 2(c). The strength of each of the dipoles (or dipole moment) is determined by the spatial derivative of the potential profile along the fiber ∂IAP(x)/∂x. Orientation of the dipoles is determined by the sign of ∂IAP(x)/∂x.

Let us consider a fiber element of infinitesimal length dx of Fig. 2(b) lying within the region occupied by the action potential. A current emerges from this differential fiber element into the extracellular region. If we assume that all the current is concentrated along the fiber axis, then this current can be expressed as dx, where$$\overline{p}$$ is the dipole current per unit length. Since this current emerges essentially from a point into an unbounded space, it behaves like a point source of current (dipole generator) that lies in an extensive conducting medium. Then, the contribution to the extracellular potential generated from this component can be expressed as

$$d{\Phi }_e=\frac{1}{4\pi {\sigma }_e}\cdot \frac{d\left(1/r\right)}{dx}\cdot \overline{p}(x,t)dx$$E2

where σ_e is the conductivity of the extracellular medium and r is the distance from excitation source to the recording point, P₀ [see Fig. 3(e)]. If the element$$\overline{p}$$dx is located at the coordinate (x, y, z) and P₀ is located at (x₀, y₀, z₀) then

$$r={\left[{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2\right]}^{\begin{array}{c}\end{array}1/2}$$E3

Normally, the coordinate origin is placed on the fiber axis, whereupon y = z = 0. The total field produced in the extracellular medium by the propagation of a single dipole along the fiber is found simply by integrating with respect to x (i.e., summing up the contributions to the potential from the propagation of this single dipole current element). The result is

$${\Phi }_e\left(x_0,y_0,z_0,t\right)={\displaystyle \underset{x=-\infty }{\overset{x=\infty }{\int }}}\frac{\overline{p}(x,t)}{4\pi {\sigma }_e\cdot {\left[{\left(x-x_0\right)}^2+y_0^2+z_0^2\right]}^{\begin{array}{c}\end{array}3/2}}\begin{array}{c}\end{array}dx$$E4

At this point, we should remember that the source of excitation is actually a “distributed source”, which means that we do not have just one dipole but a sequence of dipoles distributed equidistantly along the IAP profile. The distribution of the dipole moments (strengths) along the fiber axis is determined by the function ∂IAP(x)/∂x. Thus, to be strictly correct, the actual source of excitation is a linear axial dipole source density

$$\overline{p}=-\pi a^2{\sigma }_i\cdot \frac{\partial IAP}{\partial x}\cdot {\overline{a}}_x$$E5

where σ_i is the intracellular conductivity and a is the fiber radius. If we now substitute (5) into (4), we obtain the dipole-based expression for the extracellular potential

$${\Phi }_e\left(t\right)=\begin{array}{c}\end{array}-\frac{a^2\cdot {\sigma }_i}{4\cdot {\sigma }_e}\cdot {\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\begin{array}{c}\\ \end{array}\frac{\partial IAP\left(x,t\right)}{\partial x}\cdot {\overline{a}}_x\cdot \frac{\partial }{\partial x}\left(\frac{1}{r\left(x\right)}\right)dx\begin{array}{c}\end{array}}$$E6

In (6) it can be seen that the term ∂(1/r)/∂x represents the scalar potential generated by a propagating dipole and the term ∂IAP/∂x corresponds to the distribution of moments of the collection of dipoles that form the excitation function.

Calculation of the extracellular potential on the basis of monopoles can also be derived following the steps outlined above, as shown in Plonsey and Barr (2000).

3.2. Models of the extracellular potential generated by a single muscle fiber

3.2.1. The first SFAP models

Lorente de Nó (1947) was the first to obtain an expression for the potential of the external field of a nerve in a volume conductor as a function of an action potential. However, he did not provide a physical interpretation for the concepts of double and single layer sources. Nevertheless, his contribution to the understanding and modelling of the extracellular fields produced by an excitable fiber in a volume conductor was essential for subsequent researchers.

In 1968, Clark and Plonsey proposed a SFAP approximation based on formal solutions of Poisson’s or Laplace’s equation with the corresponding boundary conditions using a method of separation of variables. The solution was based on the Fourier transform technique and modified Bessel functions. The formal solution, however, gave no opportunity for transparent physical interpretations (Andreassen and Rosenfalck, 1981). In addition, the mathematical expressions describing SFAPs using the volume conductor theory were quite complex, computationally time consuming and, therefore, somewhat unsuitable for simulating motor unit potentials (MUPs). In response to these limitations, simplified models were presented. The dipole and tripole models approximated the transmembrane current by two and three point sources, respectively. These models, however, presented important shortcomings. The dipole model was not able to describe the effects of the excitation origin and extinction correctly (George, 1970; Boyd et al., 1978; Griep et al., 1978). The tripole models failed to correctly describe SFAPs close to fibers (Griep et al., 1978; Andreassen and Rosenfalck, 1981).

3.2.2. Assumptions of the core-conductor model and the line-source model: the first approach to convolutional models

Clark and Plonsey (1966, 1968) and then Andreassen and Rosenfalck (1981), provided thorough compilations of the approaches that relate intracellular and extracellular potentials. Both works stressed the necessity of finding the conditions to simplify existing SFAP models. These conditions are summarized in the assumptions of the core-conductor:

1. Axial symmetry is assumed. That is ∂/∂φ=0 (where φ is the azimuth angle) so that, at most, all field quantities are functions of the cylindrical coordinates only. In fact, transmembrane currents as well as intra- and extra-cellular potentials are considered to be functions only of x (i.e., the core-conductor is linear).

2. For a fiber in an extracellular medium of considerable extent, it is assumed that the resistance of the extracellular medium is practically 0, and so the influence of the medium surrounding the active fiber is neglected.

3. The excitation source is assumed to be distributed along the axis of the fiber. Since in general fiber radius is many times smaller than fiber length, this approximation is normally justified.

On the basis of the volume conductor theory, Plonsey (1974) was the first to show that a SFAP can be expressed as a convolution of an excitation source and a weight function. An important step towards the simplification of this model was made by Andreassen and Rosenfalck (1981) who assumed that the transmembrane current was distributed and concentrated along the axis of the fiber (a line-source model). In addition, the authors theorized that the error caused by such a source simplification would be less than 5% provided the radial distance was 5 fiber radii or greater from the fiber axis. However, both the approach of Clark and Plonsey and that of Andreassen and Rosenfalck were still far from being simple and intuitive as they included an intricate mathematical formulation.

3.2.3. Towards an easy formulation of the SFAP convolutional model

After the key contributions of Clark and Plonsey and then Andreassen and Rosenfalck, great effort was directed towards including in SFAP models the effects of the finite size (Gootzen et al., 1991), inhomogeneneity (Dimitrova and Dimitrov, 2001), and frequency dependency (Albers et al., 1988) of the volume conductor. All of these authors were more interested in improving the accuracy of the SFAP model at the expense of reducing the physical transparency and the time efficiency of their approximations. However, these latter two considerations were precisely two features that scientists in other areas of EMG research and practice were increasingly demanding from SFAP models.

The key step towards an easy and intuitive presentation of the SFAP convolutional model was made by Nandedkar and Stalberg in 1983. The authors substituted the complex weight function used by Andreassen and Rosenfalck (that involved Bessel functions) for the simpler expression of a potential generated by a point current source (i.e., a monopole) that propagates along the fiber axis from the neuromuscular junction towards the tendons.

Although the Nandedkar and Stalberg model was found computationally more efficient and more accurate than the dipole and tripole models, it still lacked a suitable approach for dealing with the excitation onset and extinction. To overcome this problem, the excitation wave should be represented with two stacks of double-layer disks distributed equidistantly along the fiber axis. Following this idea, Dimitrov and Dimitrova (1998) replaced the stacks of distributed current dipoles by dipoles lumped along the axis of the fiber (line source model). The weight function was computed as the potential produced by two lumped current dipoles propagating in opposite directions from the endplate toward the fiber ends.

Let us assume that a unit dipole originates at time zero at the neuromuscular junction and propagates along a fiber with a constant velocity v. The potential produced by this unit dipole at the electrode is IR(t) [Fig. 3(a)]. However, the source of excitation is actually a sequence of n dipoles distributed equidistantly along the IAP profile, whose strengths (moments) are given by the function ∂IAP(x,t)/∂t. The first dipole originates at time zero, and each subsequent dipole at an interval ∆t. Let the amplitude of these dipoles be a₁, a₂,…, a_n [Fig. 3(b)]. The first dipole originates at time zero, propagates towards the tendons, generating a potential a₁ IR(t). The second dipole originates at time ∆t, propagates towards the tendons, generating a potential a₂ IR(t-∆t) [Fig. 3(c)]. Hence the total potential recorded by the electrode is

$$SFAP(t)=a_{1}IR\left(t\right)+a_{2}IR\left(t-\Delta t\right)+\mathrm{...}+a_{n}IR\left(t-\left(n-1\right)\Delta t\right)$$E7

As ∆t tends to zero and n tends to infinity, the total potential produced at the point P₀ by the propagating excitation function represents the convolution of dIAP(x,t)/dt and IR(t):

$$SFAP(t)=\frac{dIAP\left(t\right)}{dt}*IR\left(t\right)$$E8

Figure 3.

Presentation of the potential generated by a single fiber (SFAP) (d) as the output of a linear timeshift invariant system whose impulse response IR is the potential produced by a moving unit dipole (a) and whose input signal is a function of distributed dipoles whose strengths (moments) are determined by dIAP(x,t)/dt (b). In (c) the contribution of each dipole to the total SFAP is shown separately. (e) Schematic representation of a muscle fiber innervated by the axon of a motorneuron.

Thus, an active fiber can be considered as a linear timeshift-invariant system for the generation of an extracellular potential. Specifically, the SFAP is the output of the system with the impulse response IR(t) [Fig. 3(d)] and input signal dIAP(x,t)/dt. Since the activation of a skeletal muscle fiber gives rise to two oppositely-directed propagating IAPs, IR(t) is the sum of potentials generated at the observation point by two dipoles moving in opposite directions from the neuromuscular junction to the fiber ends, where they disappear (McGill et al., 2001; Dimitrov and Dimitrova, 1998). This sum of potentials can be described as

$$SFAP{(t)}_{}{=}_{}K\cdot {{\displaystyle {\int }_0^t\frac{dIAP(\tau )}{d\tau }\cdot \frac{1}{v}\cdot \left[I{R}_1\left(t-\tau \right)+I{R}_2\left(t-\tau \right)\right]\cdot d\tau }}_{}{=}_{}K\cdot \frac{1}{v}\cdot \frac{dIAP(t)}{dt}*IR\left(t\right)$$E9

where * means convolution and IR = IR₁ + IR₂, 0 ≤ t ≤ t_max. Constant K is equal to (a² σ_i)/(4 σ_e). The expressions for IR₁ and IR₂ and can be found in Dimitrov and Dimitrova (1998).

3.3. Models of electrodes for recording single-fibre potentials

The advent of single fiber electromyography (SFEMG) allowed investigators to analyze the shape peculiarities of the SFAP. The routine single fibre (SF) electrode consists of a montage of a 25-μm recording port (leading-off surface) referenced to the cannula (Ekstedt, 1964; Stålberg and Trontelj, 1979). When modelling SF electrodes, the leading-off surface is normally considered as a detection point, whereas the cannula is simulated through a plate electrode (Stalberg and Trontelj, 1992). As established by Ekstedt (1964), the conditions for recording an SFAP in a voluntarily activated muscle using an SF electrode are: (1) that the fibre is close to the electrode and thus gives rise to a potential with a short peak-to-peak interval (rise-time, RT), and (2) that any other fibres in the same motor unit that have coincident action potentials are sufficiently remote from the electrode that their contribution to the recorded signal is small.

4.1. Anatomy and physiology of the motor unit

4.1.1. Motoneuron, motor unit fibers, and motor unit territory

The motor unit is the entity that serves as a functional building block for the production of force and movement, both in reflex and voluntary contractions. Sherrington defined the motor unit as the set comprising a single motoneuron axon and the many muscle fibers to which that axon runs and which are hence innervated by it. The motor unit fibers (MUFs) are the different muscle fibers innervated by a single motoneuron. The numbers of muscle fibers supplied by a single motor neuron, through branching of its axons, is often referred to as the motor unit size or innervation ratio, and we will refer to this number as the motor unit fiber number (MUFN). In addition, we can consider the extent of the muscle cross-section that these fibers occupy, which is the motor unit territory (MUT), and the spatial distribution of these fibers within the motor unit territory area (MUTA). Subsequently, the motor unit can be characterization in terms of the motor unit fiber density (MUFD), measured as the number of MUFs per unit of area of its territory.

The MUFN varies greatly from one muscle to another and from one motor unit to another within the same muscle. The glycogen-depletion technique allows identification of the MUFs of a single motor unit within a muscle cross-section. Studies using this technique have demonstrated the large range of variation in MUFN (Burke and Tsairis, 1973), which extends from just a few fibers per motor unit in some muscles to thousands in others. Studies also show great variation within a single muscle: up to 100-fold from the smallest to the largest motor units.

Figure 4.

A schematic representation of the motor unit cross-section showing the MUFs of two motor units and delineating their MUTs which may overlap and may occupy only a fraction of the muscle cross-section.

By delineating the boundary enclosing the MUFs, we obtain a picture of the MUT. MUTs have been shown to be distributed over a localized region of a muscle's volume, with an approximately elliptical shape in cross-section [4] (see Fig. 4). That is, in terms of the muscle cross-section, only a fraction of its area is occupied by each motor unit (Bodine et al, 1988). This spatial feature was confirmed in humans using multilead-electrode-EMG, and using scanning-EMG (Stålberg and Antoni, 1980), with typical cross-sectional corridors between 5 and 10 mm long. Research based on glycogen-depletion techniques demonstrates that motor unit fibers of different motor units are intermingled, hence motor unit territories overlap within the muscle cross-section.

Defining the MUFD as the ratio of the MUFN to the MUTA, it is found that the ranges for MUFD are independent of the MUFN, while MUTA is strongly correlated with MUFN (Kanda and Hashizume, 1992). In reconstructions of the three-dimensional structure of motor units, MUFD is almost constant within a single motor unit when comparing the values obtained over the different cross-sections of the unit along the longitudinal axis of the muscle.

4.1.2. Spatial distribution of motor unit fibers

Another important feature related to motor unit fibers is their spatial distribution within the motor unit territory. Quantitative statistical studies in glycogen-depleted motor units by means of Monte Carlo simulation analysis (Bodine et al, 1988), suggest that MUFs are not homogeneously distributed: there are alternating regions of high and low MUF density within the MUT. High-density MUF areas are not usually located in the center of the MUT, as had been previously assumed. The various studies indicate that MUFs are arranged in small subclusters separated by holes with relatively few fibers, and researchers have suggested that if any clustering is present, it should be related to the axonal branching pattern established during development. The scanning EMG has confirmed the presence of silent areas within motor units of human muscles (Stålberg and Antoni, 1980), which most likely correspond to the holes observed in glycogen-depletion studies. Hence these results seem to agree with the long-range distribution findings in non-human vertebrates by the glycogen-depletion technique.

4.1.3. Motor end-plate zone, fiber diameter, and initiation of depolarization

Electrophysiologically, the motor end-plate is the region where the motor end-plate potential is generated, hence where the intracellular action potential of the muscle fiber starts. Thus, the relative longitudinal position of the motor end-plates of the set of muscle fibers belonging to a single motor unit is determinant in the synchronization of the single fiber action potentials. In addition, there will be different times for the initiation of the depolarization of the different motor end-plates, depending on the length of the axonal sprout innervating it. Besides, SFAPs will propagate at different conduction velocities through the muscle fiber, mainly dependent on muscle fiber diameter. All these factors (spatial configuration of end-plates, initiation of depolarization, and muscle fiber conduction velocity) affect synchronization of the SFAPs contributing to the MUP, and ultimately this will affect the shape, amplitude, and duration of recorded MUPs.

Neuromuscular junctions tend to reside in the middle part of muscle fibers, as the connection is established while a muscle fiber is still growing in both directions. The three-dimensional reconstruction of the motor end-plate zones leads to a two-dimensional membrane lying in the muscle volume. Measurements show that the width of this membrane, which corresponds to the variability in the longitudinal position of the individual motor end-plates, ranges between 6 and 10 mm in the biceps brachii (Aquilonious et al, 1984).

Muscle fiber conduction velocities can be measured in situ (Stålberg, 1966) and have a normal distribution of values for the whole of a given muscle. It is assumed that conduction velocity is directly proportional to fiber diameter, with histological analyses also showing a normal distribution for the diameters of muscle fibers.

Finally, a delay in the initiation of depolarization is caused by two factors: the axonal propagation delay, which is clearly dependent on the length of the axonal terminal branch and its propagation velocity; and the neuromuscular junction transmission delay, which has an average value, but also some variability (the “jitter”).

4.2. Models for the motor unit cross section

The first computational muscle architecture models found in the literature were restricted to the simulation of a single motor unit. In essence, if individual motor units are modeled, of the triad: MUFN, MUTA, and MUFD, two quantities can be arbitrarily fixed, while the third will be a subsidiary quantity. After considering the dimensions and number of fibers of the motor unit, the models must deal with the placement of the individual MUFs within the MUT, in order to follow a certain spatial distribution.

Figure 5.

Representation of a simulated motor unit including the most relevant anatomical and physiological factor to reproduce the architecture and physiology of the motor unit.

The first motor unit model (Griep et al, 1978) was proposed in order to study the properties of the motor unit potential (MUP) by means of computer simulation. Other researchers, (Miller-Larsson, 1980; Gath and Stålberg, 1982; Hilton-Brown et al, 1985), were concerned with the modeling of fiber density and spatial distribution of muscle fibers of a single motor unit, in order to determine the influence of architectural changes produced by neuromuscular disease on clinical electrophysiological recordings. Other single motor unit models have been developed more recently. Other studies (Nandedkar et al, 1988; Stålberg and Karlsson, 2001) proposed a muscle model to study the correlation between anatomical parameters and MUP signals by means of simulations, allowing the study of the MUP variations under different pathological conditions such as denervation and re-innervation.

Two main approaches are used to model the spatial distribution of MUFs: random location and random selection from a muscle fiber grid. In the first approach, a number of MUFs are placed following a given spatial distribution, usually normal or uniform. In the second approach, a predefined grid of evenly distributed muscle fibers is created, and a number of them are selected to be innervated by the motor unit under simulation.

4.3. Models of the physiological parameters affecting the temporal dispersion of a MUAP

In all the above-mentioned models, muscle fiber conduction velocity, delay of initiation of depolarization, and longitudinal position of the motor end-plate are modeled as random variables from statistical distributions. It is important to note that, as individual motor units are being modeled, simulated distributions refer to the MUFs and not to all the fibers of the muscle bundle.

Muscle fiber diameter is modeled as a normally distributed variable with mean and standard deviation fixed to match corresponding values observed in real counterparts of the simulated muscle. Muscle fiber conduction velocity (MFCV) is modeled as a normally distributed variable, obtained from a linear transformation of the muscle fiber diameter. The most widely used relationship between fiber diameter and conduction velocity (Nandedkar and Stålberg, 1983) is:

v = 3.7 + 0.05 ( d – 55 )

With v being the MFCV in m/s and d the muscle fiber diameter in µm. This equation assumes a linear relationship between both quantities, and a MFCV of 3.7 m/s for a muscle fiber diameter of 55 µm. Both these values are the central values of the corresponding distributions for the human biceps brachii muscle.

The delay in initiation of depolarization is usually modeled as a normal or uniform random variable, although in some cases its influence on simulation outcomes is assumed to be negligible and consequently it is not modeled and set to zero for all the muscle fibers.

Finally, the longitudinal position of the motor end-plates is modeled as a normal or uniform random variable, emulating the narrow width that the motor end-plate zone occupies within the longitudinal section of the muscle. Uniform distributions usually lead to overly complex MUPs, which seems to call for normal distributions. However, a specific model accounting for the motor unit fractions observed in scanning-EMG recordings, modeled the motor end-plates of small groups of motor unit fibers as narrower normal distributions with the mean of the distributions again distributed uniformly (Navallas and Stålberg, 2009). This model allows for a motor end-plate zone that is wider, whilst keeping the dispersion, and hence the MUP complexity, locally low.

4.4. Models of recording electrodes for intramuscular EMG

Point recording models, which are accurate enough for the simulation of single fiber EMG recordings, must be extended in order to simulate other needle electrodes where electrode poles are not small enough to be considered as points. Two main approaches are available. One is the analytical approach, which requires calculation of the integrals in order to simplify the calculations (Nandedkar and Stålberg, 1983; Dimitrov and Dimitrova, 1998). The other is the discrete elements approach, where the poles are modeled as grids of points where the potentials are calculated individually, and lately averaged to give the potential of the pole. In the case of concentric needle EMG, two poles must be modeled: the core, which represents the active recording region; and the cannula, used as the reference potential. The core, which is a small plane elliptical region, can be directly modeled using either approach. The cannula, a cylindrical region, which averages the potential over a much broader region than the core, can be simplified as a one dimensional cable structure coincident with the needle axis. In the case of macro EMG recordings, the active area is the cannula itself, and the electrode may be modeled accordingly.

There are other effects related to the needle insertion procedure that can also be modeled. Whenever the needle electrode is inserted, fibers in the way of insertion can be displaced by the needle shaft. This effect, named “fiber ploughing”, calls for an update in the fiber positions according to the displacement suffered. The electrode manipulation typically performed by an electromyographist while recording concentric needle EMG can also be modeled. The electromyographist tries to bring the core closer to active fibers, producing sharp spikes in MUPs. In the corresponding model, the electrode is allowed to move 1 mm around the original insertion point in the longitudinal direction of the needle, and the position with the shortest distance to an active fiber is selected (Nandedkar et al, 1988; Hamilton-Wright and Stashuk, 2005).

5. Modeling electrical activity of skeletal muscle

In this section, the modeling of the electrical activity of a complete and functional muscle is described. After summarizing the principal anatomical and physiological characteristics of skeletal muscle, the different elements for building a model of a surface EMG signal generated by such a muscle are presented. Together with models of motor unit activity, discussed in the previous section, these elements comprise models for the architecture and geometric organization of the motor units within the muscle, models for the motor neuron pool activity, and models for the potentials recorded at the skin surface.

5.1. Physiological aspects concerning muscle EMG

In order to better understand the anatomical and physiological scenario which the models should recreate, an overview of the arrangement of MUs within the muscle and of the organization and strategies of the motor control is first provided.

5.1.1. Anatomical and architectural organization of the motor units in skeletal muscle

When trying to analyze a muscle's cross-section as a whole, overall statistics for the MUFN, MUTA and MUFD must be provided. We can think of the muscle cross-section as a tightly packed set of muscle fibers that are innervated by a smaller set of motoneurons. The innervation process, which determines the actual set of muscle fibers innervated by each motoneuron, defines the number and spatial distribution of the MUFs, hence the MUFN and the MUTA, and, subsequently, the MUFD.

The MUFNs of the motor units in a muscle tend to distribute exponentially, with fewer motor units with lower MUFN and less motor units with higher MUFN. Indirectly, MUFN can be studied on the basis of the mechanical response of the muscle. Studies of the mechanical response of skeletal muscle support the idea that many factors influence the force production of a single motor unit, including the MUFN, the cross-sectional area of motor unit fibers, and the specific tension output for the different muscle fiber types. However, within motor units of a single type, the main factor determining force production is shown to be the MUFN (Bodine et al, 1987). This makes it reasonable to assume that MUFN is proportional to motor unit maximum tetanic tension (Fuglevand et al, 1993). All these findings support the use of maximum tetanic force estimation techniques in order to investigate the MUFN in humans, where glycogen-depletion techniques obviously cannot be applied. The relative MUFNs can be assessed by measuring the sizes of the electrical and mechanical responses of individual motor units when motor axons are excited by threshold stimuli. Such a strategy involves determining the range of tetanic forces, and estimating the number of muscle fibers required to achieve these forces. These two parameters can be related by linear regression, which enables the estimation of MUFNs in human muscles.

Large ranges of variation of motor unit twitch force are found within single muscles, and the statistical distributions are shown to be highly skewed, with higher number of motor units with small motor unit twitch force, and lower number of motor units with large motor unit twitch force. Fuglevand et al. modeled the distribution of motor unit twitch forces by an exponential function (Fuglevand et al, 1993) that agrees highly with the experimental data (Kernell et al, 1983). Due to the high correlation found between MUFN and motor unit twitch tension, it can be assumed that an exponential function also governs the distribution of MUFN within motor units of a given pool.

The number of overlapping motor units at a given point of the muscle cross-section has been shown to range from 10 to 25 units in some human muscles (McIntosh, 2006). Glycogen-depletion studies have shown that the average MUTA differs for different motor unit types, in the order S<FR<FF, hence showing the same ordering as observed for MUFN. A strong positive correlation (ρ = 0.97) between MUTA and the maximum tension and a strong correlation between MUFN and motor unit maximum tension (ρ = 0.94) can be observed (Bodine et al, 1987). It is also found that the ranges for MUFD are different for each motor unit type and independent of the MUFN within each motor unit type, while MUTA is strongly correlated with MUFN within each motor unit (Kanda and Hashizume, 1992).

5.1.2. Hierarchical organization of motor control

Control of motor function is organized in three levels that act hierarchically and in parallel: the spinal cord, the brain stem and the motor areas of the cortex (Kandel, 1995). Each of these levels receives relevant sensory information from afferent pathways. This organization allows higher centres to give general commands, whilst leaving the control of detailed motor actions to lower level centres. The spinal cord is the lowest level of the hierarchy and contains neuronal circuits responsible for automatic motor patterns and reflexes. It contains a central region of grey matter that is occupied by millions of cell bodies of neurons of two types: interneurons and motor neurons (Guyton, 1994). Interneurons have many interconnections among themselves and with motor neurons, which form interneuronal circuits responsible for integration and feedback-based motor control functions. Motor neurons provide direct control of muscle activity by being directly innervated to muscle fibers. The MNs that innervate individual muscles are arranged in longitudinal columns forming what is known as the ‘motor unit pool’ (Kandel, 1995). The input to the motor neuron pool is the afferent information from peripheral receptors (muscle spindle, Golgi organs, Renshaw cells, etc.) and the efferent information (drive) from higher centres. The output of the motor unit pool consists on the firings of the different motoneurons in response to the different synaptic inputs they receive in the pool. Smooth coordinated movement relies on the subtle interplay of the command orders coming from the higher brain centers and the feedback information obtained by the peripheral receptors. Four interrelated mechanisms constitute the modulators of muscle activity: (1) MU recruitment, (2) motoneuron firing frequency (rate coding), (3) synchronization between pairs of MUs, and (4) the so-called ‘common drive’. The first two are the principal gears of force modulation; while the other two are only secondary mechanisms for control of muscle force output.

5.1.3. Motor unit activation and firing strategies

5.1.3.1. Motor unit recruitment

Motor unit recruitment refers to the way in which the central nervous system selects the specific motor units to come into action as muscle force is required. Henemann and colleagues, in experiments on cats, observed an orderly recruitment of the motor neurons in the pool as the stimuli increased (Henneman, 1965). Specifically, smaller motor neurons were recruited before larger motor neurons. They refer to this behaviour as the ‘size principle’. Many other works have reinforced the ‘size principle’, extending it to other species, different muscles and contraction tasks (Basmajian, 1985). Amplitude and conduction velocity of the MN impulses have been used as indirect indicators of motor neuron size. Twitch tension and MUAP amplitude have also shown strong positive correlation with recruitment order and spike amplitude (Merletti, 2004; Basmajian, 1985). As twitch tension and MUAP amplitude are known to be related to MU size, the ‘size principle’ indirectly links the order of the recruited MUs and their sizes.

The way recruitment is performed varies among muscles. In powerful muscles, such as the biceps brachii or deltoid, recruitment has been observed at least up to 80% maximal voluntary contraction (MVC) (Basmajian, 1985). On the other hand, in small muscles, as those of the hand, the pool of MUs is completely recruited for only 50% MVC. It has also been observed that when the voluntary force decreases, decruitment (deactivation of motor units) is performed in the opposite order to recruitment, in both isometric (De Luca, 1982) and dynamic contractions (Kossev, 1998).

5.1.3.2. Rate coding

Together with recruitment, rate coding is the most important mechanisms to modulate muscle force; the higher the discharge firing rate of a MU in a given task, the higher the force exerted by that MU. However different muscles and different contraction modalities exhibit different characteristics with respect to minimal and maximal firing frequencies and excitation-firing frequency curves (Basmajian, 85). Much of the current knowledge about recruitment and rate coding has been obtained thanks to the availability of reliable techniques for decomposition of MUAP trains from intramuscular EMG signals. One of these techniques is the so-called ‘Precision Decomposition’ technique, which was developed by De Luca and co-workers and included a recording system based on a quadrifilar needle electrode and programs to isolate several MUAP trains from the recorded signals (De Luca, 1993). With this system the researchers were able to track the evolution of several motor units in exercise protocols where force output was intended to follow a trajectory based on a linear increment, a constant level and a linear decrement. Various important paradigms of motor control could be formulated from the results of these experiments: (a) Lower threshold MUs tend to have lower initial firing rates than higher threshold MUs. (b) Earlier recruited MUs have higher firing rates than later recruited MUs. This occurs throughout the duration of the contraction, which gives rise to the typical ‘onion skin’ curves (Fig 6). (c) During sustained contractions the firing rates of MUs tend to decrease (De Luca 82) (Erim, 1996) (Fig 6). (d) The firing rate at decruitment is lower than at recruitment (Erim, 1996) (Fig. 6). (e) The firing rates of MUs with different thresholds tend to converge at maximum contractions (100% MVC).

All these findings indicate a hierarchy of MUs in the pool, which determines the recruitment order and the firing rates of MUs as a function of the required force (Erim, 1996).

Figure 6.

Schematic representation of the output force (red) and MUs mean firing rate evolution (Black) in an increasing ramp-sustained level decreasing ramp exercise.

Even in sustained contractions, the spikes of a MUAP train do not appear with an exact periodicity. The so-called interspike interval (ISI) variability is also an important aspect of rate coding (Merletti, 2004). When a MU is recruited, the ISI variability is relatively high and gets smaller as the firing rate increases (Basmajian, 1985).

5.1.3.3. Synchronization

MN synchronization refers to the “time-locked” spike trains delivered by two MNs in the pool over a certain time interval during a voluntary contraction. The “time-locked” (synchronous) pair of MUAP trains may be simultaneous or have a fixed lag, with a dispersion of a few milliseconds around the average lag (for example, of ±3ms (Sears, 1976) or ±5 ms (Datta, 1990). The number of synchronous spike pairs is more than what would be expected if the two MNs were discharging independently. Synchronization may have different neurological origins, but the most accepted one is the excitatory pre-synaptic potentials coming from brain stem descending neurons that branch into different motor neurons in the pool (Sears, 1976; Datta, 1990).

Using the cross-correlogram technique (Sears, 1976), Datta et. al. analyzed the synchronization of the first dorsal interosseous (FDI) muscle in isometric contractions and observed that the level of synchronization of two MU trains decreased as the difference between the recruitment thresholds of the two MUs increased (Datta, 1990). Kim et. al. studied synchronization in four muscles under the conditions of slight isometric contraction and by means of two simultaneous intramuscular recordings (Kim, 2001). The selection of muscles for study was made in order to investigate proximal-distal and upper-lower dichotomies. Results showed synchronization in the four muscles, with higher degrees in distal relative to proximal muscles, and in muscles of the upper extremity relative to the lower. Analysis of synchronization in the frequency domain with the use of the coherence spectrum revealed coherence peaks in the 1-5 and 25-30 Hz bands, which indicated a common rhythmic input to the MN pool, probably related to oscillating activity in the brain and corticospinal projections (Kim, 2001).

5.1.3.4. Common drive

The tendency of the firing rates of different MU trains to fluctuate together in a low frequency range (1-2 Hz) over the course of contractions was first described by De Luca and his team (De Luca, 1982). They observed this common behaviour of MUs on FDI and deltoid muscles during sustained, linearly increasing and linearly decreasing contractions, and called it ‘common drive’. They pointed to common excitatory inputs to several MUs in the pool as the probable origin of the phenomenon. Several subsequent studies evidenced similar effects in different muscles, with different exercises and in subjects of different ages (De Luca, 2002). Common drive has also been observed in muscles acting simultaneously in a certain task, either synergistically (De Luca, 2002) or in agonist-antagonist pairs (De Luca, 1987).

An interesting question is whether the phenomena of synchronization and common drive share a common origin. Semmler et al. recorded SEMG data from two separate fine-wire electrodes inserted into FDI muscles during slight isometric abduction. A low statistical correlation (<10%) between the indices that measured the extent of these processes was obtained (Semmler, 1997). Using simulated data, Jiang et. al. also studied the relationship between synchronization and common drive (Jiang, 2006). They found no correlation between the phenomena and a relationship of each of them to a different parameter in the proposed model (see Section 5.3.4, below). As Semmler et al., these authors concluded that synchronization and common drive must have a different physiological origin.

5.2. Models for muscle cross section

Several muscle architecture models have been proposed for EMG modeling. One of the objectives of these models is to reproduce the MUFN, MUTA, and MUFD distributions observed in real muscles. This implies modeling of the layout of the muscle fibers that form the muscle volume, sizing and placement of motor unit territories, and recreation of innervation patterns.

All models follow, to some degree, a common simulation scheme, depicted in Fig. 7. In order to compare the different models, we provide a classification of them based on the two components that most significantly affect the resulting properties of the models: the model of the motor unit territory placement, and the model of the innervation pattern. The different solutions proposed for each of the two components are identified. In addition, the parts which are common to most of the approaches are exposed.

The general procedure in all these models follows three main steps: (1) determination of the intended motor unit distributions (MUFN, MUFD, and MUTA); (2) placement of the MUTs; and (3) recreation of the innervation process by identifying the MUFs belonging to each of the motor units. In the following sections we will detail the different approaches available in these three steps.

Figure 7.

A general muscle cross-section model includes three steps in the simulation.

5.2.1. Determination of the motor unit distributions

Several parts of muscle architecture models tend to be common to all the available models. The geometry of the muscle is modeled by a cylinder with cross-sectional radius RMCS. Within the muscle cross-section, muscle fibers are assumed to be densely packed with a constant fiber density. Motor unit territories are modeled as circles on the muscle cross-section. Muscle fibers of different motor units are assumed to be intermingled, hence motor unit territories are allowed to overlap to a variable extent. To explain variation in maximum tetanic forces of motor units, physiological studies indicate a non-uniform distribution, which can be modeled by an exponential distribution (Fuglevand et al, 1993). If we accept that the number of muscle fibers within a motor unit is the main factor affecting the tetanic force variation (Bodine et al, 1997), we can assume, as well, an exponential distribution for the MUFN. A strong positive correlation between the MUFN and the MUTA has also been observed (Bodine et al, 1997). This seems to be consistent with a uniform value of the MUFD over the muscle cross-section. However, there is also evidence to suggest that MUFD distribution can depend on the motor unit type (Kanda and Hashizume, 1992). Hence, we should consider the assumption of constant MUFD as an approximation of the real properties of the entire motor unit pool. Accordingly, MUFD is usually assumed to be constant (20 fibers/mm², taking the value reported in (Burke et al, 1971)), and MUFN is assumed to follow an exponential curve:

n_i = α exp(β (i-1)) ; i = 1, …, N

where i is the motor unit index that ranks the motor units of a muscle from smaller to larger when sorted by the so called “size principle”; n_i is the MUFN of the i-th motor unit; and α and β are calculated to satisfy the MUFN range determined by n_min and n_max (the number of fibers in the smaller and larger motor units, respectively). The former equation is consistent both with the exponential distribution observed for the maximum tetanic forces of motor units, and with the linear relationship observed between the maximum titanic force and fiber number of motor units.

Finally, motor unit territory areas are usually calculated to fit both MUFN and MUFD, defined as MUFN/MUFD. In this way, as long as MUFD is assumed to be constant, MUTAs follow an exponential curve (Stashuk, 1993),

a_i = α d_MUF exp(β (i-1)) ; i = 1, …, N

where a_i is the MUTA of the i-th motor unit, d_MUF is the MUFD, and the rest of quantities are as given in (11). As MUTs are modeled as circles, (12) provides the means to calculate their radii.

5.2.2. MUT location: uniform, optimized, and developmental models

For motor unit territory placement, there are different approaches. The first three are based on independent uniform placement of territory centers within the muscle cross-section. Some authors (Fuglevand et al, 1993; Stashuk, 1993; Duchêne and Hogrel, 2000; Dimitrov et al, 2008) propose a uniform distribution of the territory centers within the muscle cross-section but with the restriction that the territories must lie completely inside the muscle cross-section. Other approach (Shenhav and Gath, 1986; Farina et al 2002; Keenan et al, 2006) is to consider a uniform distribution of the territory centers within the muscle cross-section, allowing the territories to partially exceed the muscle boundary and then cutting off the outlying part. A slight modification (Keenan and Valero-Cuevas, 2007) consists on assuming a uniform distribution of the territory centers within the muscle cross-section, allowing the territories to partially exceed the muscle boundary, cutting off the outlying part as in the previous model, and augmenting the territory radius until the inside region equals the original territory area. In general, these models tend to suffer from severe edge-effects, which lead to much higher MUT overlapping toward the center of the muscle cross-section than toward the edge (Navallas et al, 2009a; b). When the MUTs are cut-off, an additional effect of MUT area loss is present, leading to an increase of the MUFD beyond that intended (Navallas et al, 2009a; b). Although usually neglected, these effects severely affect the desired properties of the simulated muscle, and have to be taken into account.

Another approach, (Hamilton-Wright and Stashuk, 2005), is based on a seed-scattering algorithm. In this model, a grid of possible positions for the motor unit territory centers is created. The MUT center is placed at a seed position which is selected at random from this grid without replacement, and perturbed according to a bivariate normal distribution. When all the seed points have been used, the grid is replaced and the procedure of selecting the seed points without replacement begins again for the remaining MUTs. The resulting positions are refined during the innervation process such that they are the mean of the positions (center of mass) of the currently innervated muscle fibers. Similarly, the MUT radii are recalculated at the end of the innervation process.

The last approach considered here, (Schnetzer et al, 2001), consists in applying optimization algorithms to place the motor unit territories in such a disposition that the final muscle fiber density is as constant as possible. The authors presented two different algorithms. The first one places the motor unit territories at positions where the spatial variance of the muscle fiber density is minimized, and this algorithm can be used either with radius augmentation or without it. The second algorithm places the motor unit territories at positions where the muscle fiber density is minimal, and this approach can also be used with or without radius augmentation. These approaches ensure a uniform overlapping of the MUTs throughout the muscle cross-section, and a uniform muscle fiber density (Navallas et al, 2009b; 2010), which are desired properties for any muscle architecture model. As with the uniform models, cutting-off the MUTs leads to a loss of MUT area with negative effects on the simulation of the MUFD (Navallas et al, 2009b). Therefore, cut-off should be avoided in favor of the “augmented radius” versions of the algorithms.

5.2.3. MUF innervation: scatter, random, and weighted models

In the literature, some models explicitly state the algorithm used to model the innervation process. In these models, prior to innervation, a square or hexagonal grid of muscle fibers is created within the muscle cross-section. Then, for each of the muscle fibers, an innervating motor unit is selected from the set of units whose territories cover the fiber position. This assignment can be done in a completely random manner (Stashuk, 1993) or by weighting the innervation probabilities. The idea behind such assignments is to model non-homogeneous distributions of muscle fibers within the motor unit territory (Cohen et al, 1987; Hamilton-Wright and Stashuk, 2005; Navallas, 2010), or to control the number of innervated fibers (Hamilton-Wright and Stashuk, 2005). A further refineed algorithm (Navallas, 2010) uses an analytical expression of the probability distributions of the outcomes (MUFN, MUFD, and MUTA) to determine the innervation probabilities in order to satisfy certain target distributions.

Other approaches simply state that muscle fibers of a given motor unit are scattered within the motor unit territory (Duchêne and Hogrel, 2000; Stålberg and Karlsson, 2001; Farina et al, 2002; Dimitrov et al, 2008; Keenan et al 2006) with no prior grid of muscle fibers within the muscle cross-section. Often, papers describing scatter models do not provide the algorithms used for the placement of the motor unit fibers within the motor unit territory. The reader might guess that one of two possible approaches have been used: we can assume that the exact number of motor unit fibers is assigned to each motor unit, or that the number of motor unit fibers assigned to each motor unit is calculated to satisfy the expected MUFD. Note that, as long as MUTA remains unchanged, both approaches would be equivalent. In both cases, the final motor unit fiber positions can be obtained at random from a uniform distribution over the motor unit territory. This approach should imply the use of a particular mechanism to avoid collisions between muscle fibers, solving the so called “fiber packing problem”.

All of these models are adequate as long as optimized-augmented MUT placement has been carried out previously (Navallas et al, 2009b). However, only the inverse-model for weighted innervation ensures that exactly the intended MUFN, MUFD, and MUTA distributions will be obtained for the simulated muscle (Navallas et al, 2010).

5.2.4. Motor unit temporal dispersion parameters

If whole-muscle distribution characteristics are used in modeling the distributions of muscle fiber conduction velocities and motor end-plate locations of the motor unit, the result is usually overly complex MUPs. Muscle fiber diameters can be modeled to follow a normal distribution within individual motor units but with an increasing mean and standard deviation as a function of the motor unit index (Hamilton-Wright and Stashuk, 2005). The rationale behind this approach lies in the differences in diameters of fibers of different types. With this approach, the variability in muscle fiber conduction velocities within individual motor units is narrower, leading to more accurate levels of MUP complexity, while the overall distribution for the whole muscle follows a wider and almost normal distribution, as expected.

Generally, the approach to modeling end-plate locations is the same as that used for individual motor units, hence drawing the locations from normal or uniform distributions as described in section 3.3. Thus, end-plate locations are assumed to be uncorrelated between different motor units. However, more complex approaches (Navallas and Stålberg, 2009) can also be used to ensure further realism in the degree of complexity of the simulated MUPs, although end-plate locations are always independent from one motor unit to another.

Axonal delays and the mean neuromuscular transmission delay, as in the case of single motor units, are usually neglected and left unmodeled. However, in more accurate models (Hamilton-Wright and Stashuk, 2005), the neuromuscular junction transmission delay variability is accommodated. This model of “jitter” values, in which individual transmission delay variability can be modeled using the actual values found in real single-fiber EMG recordings, provides simulated jitter values that are highly correlated with those observed in reality.

5.3. Modeling motor unit activation and firing strategies

5.3.1. Models for recruitment

Given a set of motor units of known sizes in the pool, the ‘size principle’ straightforwardly defines the order in which these motor units will be recruited with increasing excitation and decruited with decreasing excitation. The model, explained in (Fuglevand, 1993) and applied in other simulation studies (Yao, 2000; Farina, 2002; Zhou 2004; Gabriel 2009), aims to configure the recruitment threshold excitation (RTE) so that many MUs are assigned low thresholds, while relatively few MUs are assigned high thresholds. An exponential law similar to equation (XX) for the distribution of MUFN in the pool (see Section 5.2.1), was proposed for the RTE values. A motoneuron will remain inactive as long as the excitation target force (Fuglevand, 1993) or torque (Farina, 2002) level in the pool is lower than the motoneuron's threshold, and will start firing when the excitation level reaches that threshold.

A new perspective for recruitment modelling was offered by Wakeling (2009), who introduced two input functions, one excitatory and one inhibitory, for governing recruitment. When the excitation (demanded force) is increased from zero to a maximum level, this mechanism produces the recruitment of motor units of increasingly higher thresholds and the decruitment of motor units with low thresholds.

5.3.2. Models for rate coding

Rate coding models are characterized by three different elements concerning the pool of MUs: the minimum firing rate, the excitation-firing rate curves and the inter-spike interval (ISI) variability (Fuglevand, 1993).

Minimum. firing rate

Although there is a tendency for lower threshold MUs to have lower initial firing rates than higher threshold MUs (Erim, 1996), minimum firing rates are similar for all MUs in the pool (Monster, 1977; Fuglevand, 1993). In (Fuglevand, 1993; Yao 2000; Farina 2002; Zhou 2004) the minimum firing rate was given a constant value of 8 Hz for the whole pool of MUs.

Excitation-firing. rate curves

After the experimental work of Milner Brown (1973), the relationship between the firing rate of the MUs of the pool and the excitation has been modelled by means of a linear function (Fuglevand, 1993; Yao 2000; Farina 2002; Zhou 2004). Different behaviours in relation to the peak firing rates (PFR) of MUs, in combination with the slope of the excitation-firing rate curves (SEFRC), has led to different models for these curves:

1. SEFRC is the same for all MUs, and PFRs are in inverse proportion to the recruitment threshold (Fuglevand, 1993; Zhou, 2004) (Fig 8.A). This model is based on observations made in cats, monkeys and humans (Fuglevand, 1993) and is consistent with the ‘onion skin’ phenomenon (Erim, 1996).

2. SEFRC is the same for all MUs and the PFR is the same for all the MNs in the pool (Fuglevand, 1993; Farina 2002) (Fig 8.B). Several experimental works support this strategy (Fuglevand, 1993).

3. SEFRC is the same for all MUs and PFRs are proportional to MU recruitment thresholds (Fuglevand, 1993; Zhou, 2004) (Fig 8.C). In this case, peak firing rates are related to the mechanical outputs of the MUs in the sense that large force, fast contracting MUs are assigned higher PFRs than small force, slow contracting MUs. This model is also based on experimental observations (Fuglevand, 1993).

4. SEFRC increases with the recruitment threshold, so that all MUs finally reach the same PFR at maximum excitation (Zhou, 2004) (Fig 8.D).

Figure 8.

Linear excitation-firing rate curves.

ISI. variability

In most cases a Gaussian distribution has been confirmed to best fit the experimental data (Merletti, 2004; Fugglevand, 1993). But, a Poisson distribution and a gamma distribution have also been proposed (Merletti, 2004).

Physiological. models for the MU firing rate

A different model for firing rates of MUs was proposed by Matthews (1996). He considered that the repolarization phase of the membrane potential could be modelled by an exponential curve. When this curve crosses a certain threshold, a new action potential is fired. This leads to the repetitive firing of action potentials with a constant firing rate, which depends on the exponential decaying factor and on the threshold. Higher thresholds would lead to higher firing frequencies and vice versa. White Gaussian noise of zero mean was superimposed on the membrane exponential curves, thus introducing variability into the interspike intervals, which was directly related to the power of the noise.

Jiang et al. (2006) proposed a model for the generation of action potential trains in a small set of neurons, which included excitation neurons, motoneurons and synapses. In the particular example developed, one excitation neuron provided common input signals to two different MUs through corresponding synapses, which also received feedback information from the MN outputs (Fig 9). To compare the outputs of the model to data from real experiments, SEMG recordings were obtained from biceps brachii and abductor pollicis brevis muscles during slight isometric contractions. Real and simulated signals showed similar results regarding MN synchronization and common drive (see below).

Figure 9.

Jiang’s model for the generation of action potential trains. I_app: applied currents, responsible of the neuron excitability, I_n: noise currents, V_tm: output potentials.

5.3.3. Models for synchronization

Yao et. al. proposed a model for MU synchronization which basically adjusts the firing instants of MUs in the pool so that they attain a certain degree of time proximity among them (synchronization) (Yao, 2000). This model is used to study the influence of synchronization on SEMG and output force as the neural excitations varies. Simulation results show that both the SEMG level and the variability of output force increase with synchronization, but the level of output force itself is not significantly influenced by it. The same conclusions were reported in the simulation study carried out by Zhou et. al. (2004)using Yao’s model. This model was also used by Gabriel and Kamen (2009) in their inverse modelling study conducted to find out the physiological strategies responsible for elevating the force level in isometric voluntary contractions in biceps brachii. They concluded that either rate coding or synchronization could provide output data fit to the real data and that either of these two strategies, or a combination of the two, could be involved in the motor process.

A different model was proposed by Kleine et al. (2001), who slightly modified Matthews’ firing rate model to introduce synchronization in a controlled way. In essence, the noisy input component is divided into two parts, one which is common to other MUs in the pool and one that is unique and independent of the other MUs. Simulation of SEMG signals correctly predicted the findings observed experimentally in isometric contractions of the trapezius muscle.

5.3.4. Modeling the common drive

Jiangs’ physiological model for modeling the generation of MN firing trains also enables the simulation of common drive (Jiang, 2006). In the example explored (Fig. 9), there are independent MN inputs that determine their excitabilities (I_app1 and I_app2). When these inputs are given a common oscillating signal, emulating interneural or afferent signals reaching the two MNs, their firing rates exhibit a clear correlation (common drive), although synchronization is not appreciably affected.

5.4. Models for the potentials recorded at the skin surface

Models for the potentials recorded at the skin surface usually follow the convolutional approach, that is, they try to find a ‘weighting function’ that provides the potential recorded at any point on the surface of the skin above the active fiber and caused by an elementary current source at the fiber. First of all, a geometric model of the muscle has to be defined, which might accommodate the existence of several tissue layers with different electrical properties (conductivities and capacitances). Together with the ‘weighting function’, a source function of the distribution of current sources in the fiber, such as one of those proposed for the single fiber case (Rodríguez et al, 2012), has to be used. Apart from the infinite homogeneous volume conductors, one of the first approximations for modeling the geometry of the muscle was the semi-infinite structure, which divides the complete (infinite) space into two parts separated by the skin plane: one with finite conductivity representing the muscle tissue, and one with zero conductivity, representing the other side of the skin plane (Merletti et al, 1999).

Muscle tissue presents higher conductivity in the direction longitudinal to the fibers than in the perpendicular direction. To include this behaviour into EMG models a new coefficient was included in the formulation of the ‘weighting function’: the anisotropy ratio (ratio between conductivities in the longitudinal and perpendicular directions) (Merletti et al, 1999). The effects of fat and skin layers, that have different electrical properties (conductivity and capacity), have also been incorporated into some models (Farina, 2001), (Block, 2002; Lowery 2002). The fat layer is normally considered isotropic, with conductivity appreciably lower than the muscle. Skin has a laminar structure with a highly resistant stratum corneoum and a deeper granular tissue with higher conductivity. However, it is normally modelled as a simple layer with homogeneous conductivity, although there is not general agreement about the values of skin conductivity ), (Block, 2002; Lowery 2002). The effects of these layers have been studied through simulation; relative to models which only include one or two layers, multi-layer models generate potentials with peak amplitudes closer to those found in real recordings.

An important step forward in the construction of more elaborated EMG models is the inclusion of finite limb geometries. Cylindrical muscle models have been developed by several authors (Gootzen, 1989; Roeleveld, 1997; Farina, 2004) in which the fibers run parallel to the cylinder axis. But, fibers may also run radially within a cylindrical geometry, for example, in the anal sphincter (Farina, 2004) or have different fiber-pinnation angles (Mesin, 2004). More complicated geometries, which include bones and vessels, have also been included in the models (Mesin, 2008). As the geometrical structure and composition of layers of the EMG model is made more complicated, defining and solving the electrical equations of the problem becomes more difficult. Iterative computational approaches such as the finite-element method (FEM) or boundary element method (BEM) are called for. In (Lowery, 2002) a FEM model with cylindrical geometry was devised for a muscle. This included the muscle tissue, fat and skin layers and a bone, all of them with specific conductivities. Similarly, a FEM model with a realistic geometry taken from magnetic resonance images of a particular subject’s muscle was also modelled (Lowery, 2004). Simulated signals from both models were compared to real EMG data from electrical stimulation of the upper arm. Both models presented similar features with regard to peak amplitude and power spectrum mean frequency as functions of the recording position. However, the more realistic model (Lowery, 2004) provided action potential shapes closer to those actually recorded (Lowery, 2002).

Finally, the effects of the surface EMG electrodes potentials should also be included in the model. In general, placing an electrode on the skin surface does not alter the potential field. This is due to the relatively high impedance “seen” by the conductor tissue, which is, in turn, due to the electrochemical double layer formed between the metal of the electrode and the tissue. The potential recorded by the electrode is then the average potential in the surface covered by the electrode (McGill, 2004). Analytical or numerical procedures may be used to calculate this average either in the spatial domain (Dimitrov, 1998; Merletti, 1999) or in the spatial frequency domain (Farina, 2004).

5.5. Conclusions and open research lines

A general perspective of EMG modeling has been displayed together with a description of the anatomical and functional physiological aspects, in which the described models are grounded. This panoramic view comprises models for the space and size distribution and architectural organization of MUs in the muscle; a view of the hierarchical organization of the motor control; models for the principal MU activation and firing strategies for muscle force production: MU recruitment, ‘rate coding’, MN synchronization and the ‘common drive’; and models for the generation of potentials at electrodes placed on the skin surface. Experimental research works sustaining evidences for the theories and concepts described in the chapter have also been included.

The research effort in modeling and simulating EMG signals in the last three decades has been paramount including both analytical as well as numerical orientations. The degree of complexity and detail has also run parallel to these developments with the aim of recreating the physiological EMG generation system on one hand and the temporal and spectral features of real EMG signals on the other hand. However, there is still room for improvement in several aspects:

1. Most experimental observations related to force modulation mechanisms are referred to isometric contractions. The extent of the validity of existing models of such mechanisms in dynamic contractions is not solidly established yet, deserving more research attention.

2. SMEG generation models with increasingly complex muscle geometry and fibers disposition have been developed (Messin 2004; Messin 2008; Lowery 2002). On the other hand, sophisticated models for the MU geometrical distribution and MUF innervations have been put forward (Navallas 2010). However, these two different and complementary modeling pieces have not been put together yet, as a unified approach of muscle architecture. Their combination could represent even better the complex physiological arrangement of fibers and MU in the muscle, permitting more refined simulation studies.

3. The phenomena of MU substitution (Westgaard, 1999) and replacement (Bawa, 2009), which appears in prolonged contractions, can be considered as a fifth strategy for MU firing rate control. EMG models that include these phenomena are still missing.

4. Fatigue, aging and neuromuscular disease are specific circumstances that degrade muscle performance. EMG models for these situations are scarce (Dimitrov, 2008; Nandedkar, 1989; Stalberg 2001; Enoka, 2003) but may be indeed very useful for better understanding the underlying mechanisms. More research attention should therefore be given to this important area of EMG modeling.