**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.

EMG Modeling 19

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

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'

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

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

output.

*5.1.3. Motor unit activation and firing strategies* 

indirectly links the order of the recruited MUs and their sizes.

*5.1.3.1. Motor unit recruitment* 

dynamic contractions (Kossev, 1998).

*5.1.3.2. Rate coding* 

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). Glycogendepletion 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*

Computational Intelligence in Electromyography Analysis – 18 A Perspective on Current Applications and Future Challenges

of MUFN within motor units of a given pool.

*5.1.2. Hierarchical organization of motor control* 

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

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). Glycogendepletion 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).

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 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).

EMG Modeling 21

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

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,

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

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

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

either synergistically (De Luca, 2002) or in agonist-antagonist pairs (De Luca, 1987).

synchronization and common drive must have a different physiological origin.

parts which are common to most of the approaches are exposed.

and corticospinal projections (Kim, 2001).

**5.2. Models for muscle cross section** 

innervation patterns.

*5.1.3.4. Common drive* 

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*

Computational Intelligence in Electromyography Analysis – 20 A Perspective on Current Applications and Future Challenges

contractions (100% MVC).

*5.1.3.3. Synchronization* 

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

All these findings indicate a hierarchy of MUs in the pool, which determines the recruitment

Force (% MVC)

**Figure 6.** Schematic representation of the output force (red) and MUs mean firing rate evolution (Black)

Time (s)

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

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

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.*

in an increasing ramp-sustained level decreasing ramp exercise.

Mean firing rate (Hz)

gets smaller as the firing rate increases (Basmajian, 1985).

neurons in the pool (Sears, 1976; Datta, 1990).

order and the firing rates of MUs as a function of the required force (Erim, 1996).

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.

Computational Intelligence in Electromyography Analysis – 22 A Perspective on Current Applications and Future Challenges

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.

EMG Modeling 23

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

ai = α dMUF exp(β (i-1)) ; i = 1, …, N (12)

where ai is the MUTA of the i-th motor unit, dMUF 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.

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 crosssection 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

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

desired properties of the simulated muscle, and have to be taken into account.

are recalculated at the end of the innervation process.

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

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

follow an exponential curve (Stashuk, 1993),

**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 crosssection. 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/mm2, taking the value reported in (Burke et al, 1971)), and MUFN is assumed to follow an exponential curve:

$$\mathbf{n} = \alpha \exp(\beta \text{ (i-1)}) \quad \text{; } \mathbf{i} = 1, \dots, \mathbf{N} \tag{11}$$

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"; ni is the MUFN of the i-th motor unit; and α and β are calculated to satisfy the MUFN range determined by nmin and nmax (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),

$$\mathbf{a} = \mathbf{a} \text{ dɔɔɔ } \exp(\mathbf{\beta} \text{ (i-1)}) \text{ : } \mathbf{i} = 1, \dots, \mathbf{N} \tag{12}$$

where ai is the MUTA of the i-th motor unit, dMUF 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*

Computational Intelligence in Electromyography Analysis – 22 A Perspective on Current Applications and Future Challenges

*5.2.1. Determination of the motor unit distributions* 

assumed to follow an exponential curve:

fiber number of motor units.

these three steps.

muscle sizing and motor unit realizations

*muscle parameters and motor unit distributions*

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

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 crosssection. 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/mm2, taking the value reported in (Burke et al, 1971)), and MUFN is

motor unit sizing and placement

*motor units placement algorithm*

> muscle fiber sizing and innervation

*muscle fibers innervation algorithm*

 ni = α exp(β (i-1)) ; i = 1, …, N (11) 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"; ni is the MUFN of the i-th motor unit; and α and β are calculated to satisfy the MUFN range determined by nmin and nmax (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

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

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 crosssection 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.

EMG Modeling 25

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

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,

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

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

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

although end-plate locations are always independent from one motor unit to another.

**5.3. Modeling motor unit activation and firing strategies** 

expected.

in reality.

threshold.

*5.3.1. Models for recruitment* 

decruitment of motor units with low thresholds.
