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

#### *5.3.1. Models for recruitment*

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

"augmented radius" versions of the algorithms.

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

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

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

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

be obtained for the simulated muscle (Navallas et al, 2010).

*5.2.4. Motor unit temporal dispersion parameters* 

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.

26 A Perspective on Current Applications and Future Challenges

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

EMG Modeling 27

*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

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

**Figure 9.** Jiang's model for the generation of action potential trains. *I*app: applied currents, responsible of

Vtm1 Vtm2

Excitation neuron

Iapp0 In0

Synapse Synapse

n1 Iapp2 I

Motor neuron 2 n2

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

interspike intervals, which was directly related to the power of the noise.

results regarding MN synchronization and common drive (see below).

Motor neuron 1

Iapp1 I

the neuron excitability, *I*n: noise currents, *V*tm: output potentials.

*5.3.3. Models for synchronization* 

#### *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:


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

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

curves (SEFRC), has led to different models for these curves:

on experimental observations (Fuglevand, 1993).

PFR at maximum excitation (Zhou, 2004) (Fig 8.D).

B

Excitation (% max)

skin' phenomenon (Erim, 1996).

strategy (Fuglevand, 1993).

**Figure 8.** Linear excitation-firing rate curves.

Excitation (% max)

have also been proposed (Merletti, 2004).

*ISI. variability* 

A

Mean firing rate (Hz)

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

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.

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

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

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

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

d. SEFRC increases with the recruitment threshold, so that all MUs finally reach the same

C

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

Excitation (% max)

Excitation (% max)

D

*5.3.2. Models for rate coding* 

*Minimum. firing rate* 

*Excitation-firing. rate curves* 

(ISI) variability (Fuglevand, 1993).

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.

EMG Modeling 29

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

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

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

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.

those found in real recordings.

those actually recorded (Lowery, 2002).

in the spatial frequency domain (Farina, 2004).

**5.5. Conclusions and open research lines** 

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.
