A Comprehensive Mathematical Model of Motor Unit Pool Organization, Surface Electromyography, and Force Generation

Neuromuscular physiology is a vibrant research field that has recently seen exciting advances. Previous publications have focused on thorough analyses of particular aspects of neuromuscular physiology, yet an integration of the various novel findings into a single, comprehensive model is missing. In this article, we provide a unified description of a comprehensive mathematical model of surface electromyographic (EMG) measurements and the corresponding force signal in skeletal muscles, both consolidating and extending the results of previous studies regarding various components of the neuromuscular system. The model comprises motor unit (MU) pool organization, recruitment and rate coding, intracellular action potential generation and the resulting EMG measurements, as well as the generated muscular force during voluntary isometric contractions. Mathematically, it consists of a large number of linear PDEs, ODEs, and various stochastic nonlinear relationships, some of which are solved analytically, others numerically. A parameterization of the electrical and mechanical components of the model is proposed that ensures a physiologically meaningful EMG-force relation in the simulated signals, in particular taking the continuous, size-dependent distribution of MU parameters into account. Moreover, a novel nonlinear transformation of the common drive model input is proposed, which ensures that the model force output equals the desired target force. On a physiological level, this corresponds to adjusting the rate coding model to the force generating capabilities of the simulated muscle, while from a control theoretic point of view, this step is equivalent to an exact linearizing transformation of the controlled neuromuscular system. Finally, an alternative analytical formulation of the EMG model is proposed, which renders the physiological meaning of the model more clear and facilitates a mathematical proof that muscle fibers in this model at no point in time represent a net current source or sink. A consistent description of a complete physiological model as presented here, including thorough justification of model component choices, will facilitate the use of these advanced models in future research. Results of a numerical simulation highlight the model's capability to reproduce many physiological effects observed in experimental measurements, and to produce realistic synthetic data that are useful for the validation of signal processing algorithms.

1. Introduction

Electromyography (EMG) denotes the measurement of the electrical fields generated by the electrophysiological processes that lead to muscle fiber contraction. EMG is highly relevant for a number of clinical and scientific applications, since it enables monitoring and analysis of a muscle's electromechanical properties and state, both of which would otherwise remain mostly inaccessible. Surface electromyography (sEMG) denotes the noninvasive measurement of electrical muscle activity by means of electrodes placed on the skin surface, as opposed to the traditional measuring method using needle electrodes. For more background information on sEMG, its analysis and many of its applications, refer to, e.g., Merletti and Parker (2004) and Merletti and Farina (2016).

Mathematical models of sEMG are highly useful, on the one hand to advance understanding of the underlying physiological processes, and on the other hand to analyze the sensitivity of sEMG measurements to various physiological and technical parameters, and to test and validate sEMG signal processing algorithms. Over the past decades, researchers have pursued a number of different approaches for the modeling and simulation of different aspects of sEMG measurements. Phenomenological (Hogan and Mann, 1980; Lo Conte et al., 1994; Sinderby et al., 1995; McGill, 2004) as well as physiologically motivated (Fuglevand et al., 1993; Dimitrov and Dimitrova, 1998; Farina and Merletti, 2001; Farina et al., 2004; Dideriksen et al., 2010a,b; Mordhorst et al., 2015) models have been proposed and analyzed. Overviews can be found in Stegeman et al. (2000), McGill (2004), Rodriguez-Falces et al. (2012), and Merletti and Farina (2016).

Many researchers have worked on modeling the electric signal produced by a single contraction of a single muscle fiber, the so-called single fiber action potential (SFAP). Classically, the propagation of the action potential along a contracting muscle fiber has been modeled using simplified dipole, tripole or quadrupole models Fuglevand et al. (1992); Merletti et al. (1999); Merletti and Parker (2004); Plonsey and Barr (2007). A more general model has been proposed by Dimitrov and Dimitrova (1998), and this model has been successfully employed, modified and combined with various other models for the remaining physiological processes in several publications Farina and Merletti (2001); Farina et al. (2004); Wang et al. (2006); Dideriksen et al. (2010a).

In the present article, the SFAP model originally proposed by Dimitrov and Dimitrova (1998), and subsequently extended by Farina and Merletti (2001), is combined with the well-known motor unit (MU) pool organization model of Fuglevand et al. (1993) and the twitch force parameterization used by Raikova and Aladjov (2002). Care is taken in particular to achieve a consistent parameterization of the electrical and the mechanical components of the model, resulting in a realistic EMG-force relationship of the simulated muscle. Recent results regarding the modeling of MU rate coding and recruitment (De Luca and Hostage, 2010) and the variability of the inter-spike intervals (Moritz et al., 2005) are incorporated, and a new model of motor unit firing rates is proposed. Moreover, we propose a novel nonlinear transformation of the common drive input to the muscle, which ensures that the desired muscle output force can directly be used as a model input. This is achieved by adjusting the rate coding model to the force generating capabilities of the simulated muscle and has interesting consequences for the modeling of physiological force control, which will be discussed as well. Finally, an alternative analytical formulation of the SFAP model of Farina and Merletti (2001) is proposed, which clarifies the physiological meaning of the model. Based on this alternative formulation, a proof is provided that in this model the ingoing and outgoing currents along each fiber sum to zero at all times, which is a physiologically plausible property due to the quasi-static behavior of action potential generation (Plonsey and Barr, 2007).

In section 2, all components of the mathematical model are presented briefly, yet completely, and in a unified way. Several mathematical properties of the model are derived, and the above-mentioned alternative formulation of the model of Farina and Merletti (2001) is proposed. Results of numerical simulations based on the model are presented in section 3 and are assessed with respect to their physiological plausibility. Finally, section 4 concludes the article with a discussion of the various improvements introduced in this article. Note that a preliminary version of the mathematical analysis of the sEMG model of Farina and Merletti (2001) presented in section 2.5 has been the subject of a conference publication (Petersen, 2016).

2. Mathematical Model

The fundamental functional unit of a skeletal muscle is the motor unit (MU), comprising a motor neuron and the muscle fibers innervated by that neuron. The following sections introduce mathematical models of the electrical and mechanical properties of MUs, as well as their organization in a muscle. Figure 1 shows a graphical summary of the main model components and their interactions and may provide a useful reference for the reader while following along the description of the model.

FIGURE 1

Figure 1. Block diagram illustrating the main components of the proposed model of muscular force generation. The model input is the mean normalized desired muscle force ${\widetilde{\stackrel{̄}{F}}}_{\text{ref}}(t)$ with values in [0, 1], the model outputs are the total generated muscle force F(t) and the measured EMG signal. Firing rates of individual MUs are denoted by λ_i(t), the force contribution of each MU by F_i(t), and the EMG contribution of each MU by EMG_i(t). The common drive input CD(t) to the rate coding model is calculated from the desired muscle force by means of a static, nonlinear function, which effectively adjusts the rate coding and recruitment model to the force generation model.

2.1. Motor Unit Pool Structure

Every muscle consists of a number n of MUs. Each MU has various mechanical and electrical properties, most of which have been found to be closely related through the size principle: MU size as measured by the number of fibers contained in the MU is roughly proportional to force twitch amplitude, EMG twitch amplitude, and recruitment threshold (Henneman, 1957; Henneman et al., 1965; Heckman and Enoka, 2012). This means that larger MUs are only activated at higher levels of desired muscle force compared to smaller units, but they also add larger force and EMG contributions to the muscle output once activated. This orderly recruitment appears to be a result of an increase in input resistance with MU size (Powers and Binder, 2001; Heckman and Enoka, 2012) and seems to remain remarkably stable over a wide range of MU conditions (Heckman and Enoka, 2012). Based on these findings, the classical MU pool model of Fuglevand et al. (1993) that has been reused and extended in numerous studies (Stegeman et al., 2000; Zhou and Rymer, 2004; Keenan and Valero-Cuevas, 2007; Dideriksen et al., 2010a; Robertson, 2014; Potvin and Fuglevand, 2017) describes the MU twitch forces and contraction times as a function of the MU's recruitment threshold. Recent results have shown that various other parameters are also related to MU size, and hence can be considered size principle parameters as well (Heckman and Enoka, 2012; Del Vecchio et al., 2017). In the following, we will describe a size-dependent MU pool parameterization model that is based on the model of Fuglevand et al. (1993), but extended in several important aspects. Currently, the most common assumption regarding the neural excitation of MUs is that all MUs in a muscle receive the same net level of neural input, usually termed the common drive to the muscle (De Luca and Erim, 1994; Erim et al., 1996; Piotrkiewicz and Türker, 2017). We therefore formulate our model as a function of the common drive CD(t), with 0 ≤ CD(t) ≤ 1.

The recruitment thresholds of MUs in a muscle appear to follow a continuous distribution with many MUs attaining a small recruitment threshold, and few large MUs only being recruited at high activation levels (Van Cutsem et al., 1997; Raikova et al., 2007; Heckman and Enoka, 2012). This behavior is captured well by the exponential model proposed by Fuglevand et al. (1993), which assigns the recruitment thresholds

$\begin{array}{ll}{\text{CD}}_{\text{rec}(i)}=\frac{e^{ai}}{100}\text{ with }a=\frac{\ln (100·{\text{CD}}_{\text{full}})}{n}& (1)\end{array}$

to MUs i = 1, …, n, where CD_rec(i) denotes the minimum level of common drive at which the MU is recruited, and CD_full denotes the point of full recruitment, i.e., the level of common drive at which all MUs are recruited. Alternatively, following the formulation of De Luca and Contessa (2012), the thresholds can also be modeled as

$\begin{array}{ll}{\text{CD}}_{\text{rec}(i)}=\frac{bi}{n}·\frac{e^{ai}}{100}\text{ with }a=\frac{\ln (\frac{100\text{·}{\text{CD}}_{\text{full}}}{b})}{n},& (2)\end{array}$

where b denotes a scaling factor that influences the shape of the distribution. The latter model results in a more gradual slope compared to the first one and has been used by De Luca and Contessa (2012) for modeling the characteristics of the Vastus Lateralis (VL) muscle. The choice between Equations (1) and (2), in general, should be based on the characteristics of the specific muscle under consideration. With the recruitment thresholds set, peak twitch forces are calculated as a linear function of the recruitment thresholds following

$\begin{array}{ll}{P}_i=P_1+\frac{{\text{CD}}_{\text{rec}(i)}-{\text{CD}}_{\text{rec}(1)}}{{\text{CD}}_{\text{full}}-{\text{CD}}_{\text{rec}(1)}}\text{ · }(P_n-P_1)& (3)\end{array}$

as proposed by Contessa and De Luca (2013), where the twitch peak range (P_n/P₁) is typically large, e.g., P_n/P₁ = 130 for the First Dorsal Interosseus (FDI) muscle (Fleshman et al., 1981; Contessa and De Luca, 2013).

Equivalently to Equation (3), Fuglevand et al. (1993) modeled the number η of innervated muscle fibers, which appears to be the main factor influencing MU twitch force (Totosy de Zepetnek et al., 1992), directly proportional to the peak twitch force and hence also to the recruitment threshold. This direct proportionality, however, neglects the increase in fiber diameter from small to large MUs, which also contributes to the overall rise in MU twitch force (Burke, 1981; Kernell, 2006), since specific fiber tension (tension/area) does not change significantly (Lucas et al., 1987; Heckman and Enoka, 2012). The fiber diameter also determines, among others, action potential amplitude (Hakansson, 1956; Nandedkar and Stålberg, 1983), and electrical twitch conduction velocity (Hakansson, 1956; Nandedkar and Stålberg, 1983; Sadoyama et al., 1988). We thus consider both peak fiber twitch force and single fiber action potential (SFAP) amplitude to be size principle parameters as well, with the former being proportional to the square of the recruitment threshold CD_rec (since it is related to fiber area) and the latter directly proportional to CD_rec (since it is related to fiber diameter). Finally, a fivefold range in contraction speed has been found from the slowest to the fastest MU (Burke, 1981; Fleshman et al., 1981) and there is also a slight negative correlation of time to peak force with the recruitment threshold (Van Cutsem et al., 1997), while the electrical twitch conduction velocity v correlates positively with the recruitment threshold (Del Vecchio et al., 2017). We hence consider the electrical conduction velocity v and the force twitch model parameters T_EMD, T_ri, and T_hr (see section 2.7) to be linearly related to the recruitment threshold as well.

Fuglevand et al. (1993) modeled the relationships between parameter values deterministically, but experimental evidence demonstrates that these relations are of a highly stochastic nature (Heckman and Enoka, 2012). As did several previous studies (Contessa and De Luca, 2013; Robertson, 2014; Al Harrach et al., 2017), we thus draw parameter values for each MU at random from a stochastic distribution. The main novelty of the MU pool part of our model is that the distribution from which each property is drawn depends continuously on the MU's size. Previous modeling studies have followed a different path in assigning each MU a particular MU type or each fiber a fiber type and then drawing parameter values at random from a parameter distribution associated with this type (Robertson, 2014; Al Harrach et al., 2017), or drawing the parameters of all MUs from the same, static distribution (Contessa and De Luca, 2013). It has been argued, however, that both MU and fiber properties follow a continuous distribution that does not support a distinction between discrete types (Fuglevand et al., 1993; Enoka and Fuglevand, 2001; Heckman and Enoka, 2012; Potvin and Fuglevand, 2017). For this reason, we decided to model both MU and fiber properties as continuously distributed quantities. As proposed by Contessa and De Luca (2013), we use Weibull distributions for the parameter values. Summarizing, the MU pool model presented above describes the electrical and mechanical parameters of individual MUs as a stochastic, continuous function of MU size.

One final difference between our model and previous models regards the model input: while previous studies have always described recruitment thresholds as a function of the muscle force F(t) (Fuglevand et al., 1993; De Luca and Contessa, 2012), we define them as a function of the common drive CD(t) instead. This change is necessary since the generated muscle force depends on all properties of the MU pool, and not only on the recruitment thresholds. Finding a model parameterization that fixes the MU recruitment thresholds at particular levels of muscle force is, hence, a hard problem that requires careful adjustment of all model components. It has to the authors' knowledge not been solved so far. The relationship between the common drive and the generated muscle force, as well as an alternative solution to the problem of adjusting model components to each other, will be discussed in more detail in section 2.8.

2.2. Geometrical Distribution of Motor Units and Muscle Fibers

Muscle fibers belonging to the same motor unit spread over a territory that may span a large portion of the muscle cross section (Buchthal et al., 1957; Burke, 1981; Trontelj et al., 2004). The territories of the different MUs overlap, leading to an interweaving of fibers belonging to multiple MUs (Buchthal et al., 1957; Trontelj et al., 2004). Buchthal et al. (1957) have found motor unit territories to attain an irregular round shape, whence we propose the use of an elliptic model for the MU cross sections. With the elliptic axis ratio being fixed, the MU cross-sectional area—and thereby the axes lengths—are calculated by dividing the number of innervated fibers η by the desired MU fiber density ρ (fibers/area):

$\begin{array}{ll}A=\frac{\eta }{\rho }.& (4)\end{array}$

The midpoints of all MUs are then distributed uniformly over the muscle cross section. Note that without further assumptions, the above model directly leads to overlapping regions between MUs, which is a desirable feature, as noted above. Finally, fibers belonging to the MU are then again distributed uniformly inside the elliptic MU cross section. This model is equivalent to the propositions of Fuglevand et al. (1993), except for the means of reducing fiber density variability described next, and the fact that they used circular MU territory shapes as opposed to the more flexible elliptic shape proposed here. It is advisable to divide the muscle cross section into M parts of equal size, and then distribute n/M MUs uniformly in each part, to avoid an unrealistically high variability of the fiber density due to the random MU placement. Alternatively, MU centers could also be placed using an optimization algorithm such as the one proposed by Carriou et al. (2016b), which sequentially places MUs at the position that maximizes the distance to the already placed MUs.

A crucial decision when modeling random MU placement concerns the treatment of muscle boundaries. Those parts of MU territories that exceed the muscle territory must be cut off, and the question remains how to account for this loss in MU territory. Several approaches to solving the problem are conceivable (Rodriguez-Falces et al., 2012; Carriou et al., 2016b):

1. All fibers belonging to the MU are placed in the remaining parts of the MU territory. This approach leads to an increase in the fiber density of boundary MUs, and hence also in the overall fiber density toward the muscle boundaries.

2. The number of fibers innervated by the MU is reduced proportionally. This approach keeps the assigned fiber density constant, but reduces the number of innervated fibers and thus the size of boundary MUs.

3. The axes lengths of the elliptic MU region are adjusted in such a way as to keep the MU area at the desired value in spite of the cut-off. This approach keeps the number of fibers and the fiber density at the desired values but likely leads to strongly increased overall fiber density toward the muscle center, due to many adjusted MU regions overlapping there.

4. MUs with territories exceeding the muscle territory are rejected completely. This approach obviously removes the need for MU property adjustments, but without further modifications leads to reduced overall fiber density close to the muscle boundaries (Carriou et al., 2016b).

There are advantages and disadvantages to each approach, and it does not yet seem to be clear if one of the proposed approaches is generally superior to the others, or which approximates reality best (Rodriguez-Falces et al., 2012). However, muscle fiber diameters appear to be approximately constant throughout a muscle (Johnson et al., 1973; Schnetzer et al., 2001), whence a constant fiber density throughout the muscle cross section seems desirable. To this end, the second of the above approaches has been pursued here.

Reducing only the number of fibers in MUs close to the muscle boundary without adjusting the other MU parameters as well would disturb the relationship between the electrical and mechanical properties of these MUs. Hence, the recruitment threshold CD_rec(i), the peak twitch force P_i and the electrical twitch conduction velocity v were recomputed. Note that this, in turn, distorts the exponential distribution of MU parameters described by Equations (1) and (2). This distortion, however, was considered less grave than a disturbed electromechanical relationship in a significant number of MUs. Finally, to prevent unphysiologically sharp IAP generation and extinction artifacts, we randomize both innervation zone placement as well as fiber end placement, as described previously (Merletti et al., 1999; Carriou et al., 2016b).

An alternative model of MU placement has been proposed by Navallas et al. (2010). They explicitly consider the optimization problem of minimizing the variability of muscle fiber density throughout the muscle, regardless of variances in MU fiber density and while maintaining the exponential relationship (1). Their model has recently been extended (Robertson and Johnston, 2017) to account for the physiologically observed regionalized MU placement (Elder et al., 1982), something that has not been but could be considered in our model. While all of these are desirable properties of a MU placement algorithm, the much simpler division of the muscle region into separate parts and random placement proposed above also reduces the variability of the fiber density. MU regionalization could be implemented using this same division into separate parts; and the distortion to the exponential relationship (1) when using this algorithm has been found to be rather small in practice, see section 3. In summary, both the placement algorithms of Navallas et al. (2010) and the one proposed here represent viable modeling choices, with the former being preferable if the particular influence of different aspects of MU geometry is of interest, and the latter being a much simpler algorithm.

2.3. Firing Rates

The firing rate λ_i(t) of a motor unit denotes the frequency at which its motoneuron discharges, thereby generating both EMG and force twitches in its innervated muscle fibers. Together with the force and EMG twitch amplitudes of each MU, the firing rates are the primary determinant of the generated muscle force and surface EMG. Following the experimental results of Milner-Brown et al. (1973), firing rates have traditionally been modeled as a linear function of excitatory drive (Fuglevand et al., 1993; Erim et al., 1996). Recently, De Luca and Hostage (2010) have proposed a linear-exponential firing rate model based on new experimental findings. In this model, the firing rate of a MU with recruitment threshold CD_rec(i) is given by

$\begin{array}{ll}\text{λ}(\text{CD}(t);{\text{CD}}_{\text{rec}(i)})=C_4\text{·CD}(t)+(C_3-C_1{e}^{-\text{CD}(t)/C_2})\text{·}{\text{CD}}_{\text{rec}(i)}+C_5,& (5)\end{array}$

where 0 ≤ CD(t) ≤ 1 denotes the current level of common drive to the muscle, and the C_j, j = 1, …, 5 are constant shape parameters. For example values of the shape parameters for different muscles, refer to De Luca and Hostage (2010). Figure 2A shows exemplary firing rate characteristics obtained using Equation (5).

FIGURE 2

Figure 2. Firing rates of every tenth MU in one belly of one of the simulated recti, following Equation (5) (A,B) and Equation (6) (C,D). Firing rates are shown as a function of the common drive CD(t) to the muscle (A,C) and as a function of the normalized muscle output force $\tilde{F}$ (B,D). The point of full recruitment was set to CD_full = 0.6. Note that the application of the proposed nonlinear input transformation $\text{CD}={\tilde{F}}^{-1}({\tilde{F}}_{\text{ref}})$ has shifted the point of full recruitment to a much higher force level. This shift and the noticeable upwards bend at high activation levels are a result of adjusting the rate coding model to the force generating capabilities of the muscle.

This model reproduces many phenomena observed experimentally, such as the onion-skin phenomenon: MUs recruited first appear to attain a higher firing rate than those MUs recruited later throughout the whole contraction (De Luca and Erim, 1994; Erim et al., 1996). Interestingly, however, the model also differs from previous results and models in several ways. Firstly, the model does not feature the increased slope of the firing rate characteristics at excitation levels exceeding the point CD_full of full recruitment which has been observed by De Luca and Erim (1994). Secondly, in the observations and the model of De Luca and Hostage (2010), the initial firing rates of motor units are declining with increasing recruitment threshold, while in the earlier study of Erim et al. (1996), they are increasing. Thirdly, the firing rates, especially of the smaller motor units, appear unusually large in the low activity range. These high rates may result from the limitation of the regression data, from which the model was derived, to the range from 20 to 100% MVC, as has been noted by the authors in a subsequent corrigendum (De Luca and Hostage, 2012). While the parameters in Equation (5) of course can be tuned to follow new findings or represent other muscles, the shape of the characteristics cannot be altered arbitrarily due to the small number of model parameters. Finally, other researchers have questioned the validity of the EMG decomposition approach employed by De Luca and Hostage (2010) to obtain the firing rates of individual MUs (De Luca and Nawab, 2011; Farina and Enoka, 2011; Farina et al., 2015), a fact that warrants caution when reusing results obtained using this decomposition approach.

To facilitate the inclusion of the various observations discussed above into a single model, and since the model of De Luca and Hostage (2010) cannot be adjusted in all of these regards, we propose the following, novel model of MU firing rates:

$\begin{array}{l}\text{λ}(\text{CD}(t);{\text{CD}}_{\text{rec}(i)})=-C_1·(C_2-\text{CD}(t))·{\text{CD}}_{\text{rec}(i)}+C_3·\text{CD}(t)\\ +C_4-(C_5-C_6·{\text{CD}}_{\text{rec}(i)})·e^{-\frac{\text{CD}(t)-{\text{CD}}_{\text{rec}(i)}}{C_7}}\text{}.& (6)\end{array}$

Figure 2C shows exemplary firing rate characteristics obtained using this model. This model (with suitably chosen parameter values C_j, j = 1, …, 7) fulfills the following requirements:

• An initial, steep slope of the characteristics is followed by a flatter, linear region, and the transition between the two regions is smooth. Neither the models of Fuglevand et al. (1993) nor Erim et al. (1996) are smooth. Also, the model of Fuglevand et al. (1993) does not feature the initial phase with steeper slopes.

• At each activation level CD(t), the characteristics fulfill the onion skin property dλ(CD; CD_rec(i))/dCD_rec(i) < 0.

• The slope dλ(CD = CD_rec(i); CD_rec(i))/dCD_rec(i) of the initial firing rates can be freely adjusted to either positive or negative values. This is not the case for the model of De Luca and Hostage (2010).

• The degree of convergence of the firing rate characteristics of earlier and later recruited MUs can be freely adjusted by modifying C₂. This is not possible using the model of De Luca and Hostage (2010). The importance of this issue has been discussed in detail by Fuglevand et al. (1993).

Note that as for the recruitment thresholds, we define the firing rates as a function of the common drive, as did Fuglevand et al. (1993) and unlike De Luca and Hostage (2010). The generated muscle force depends on all model parameters, and hence firing rates at a particular level of muscle force must be chosen in accordance with the force twitch parameters of the involved MUs. We will address this issue in section 2.8 by introducing a suitably defined procedure for calculating the common drive input CD(t) to the rate coding model as a function of the desired muscle force level ${\tilde{F}}_{\text{ref}}$, thereby adjusting the rate coding model to the force generating capacity of the model. One benefit of the nonlinear input transformation proposed in section 2.8 is that it effectively results in an increased slope of the firing rate characteristics at excitation levels exceeding the point of full recruitment CD_full, a desirable feature as was mentioned above.

2.4. Firing Instants

Given the instant t_i(j−1) of the (j − 1)th firing of MU i and the time course of the MU's firing rate λ_i, the next firing instant t_ij can be calculated. To model the stochastic distribution of the inter-spike intervals, these are assumed to follow a normal distribution with a coefficient of variation that decreases with increasing activation, following

$\begin{array}{ll}{c}_{v,i}(\text{CD}(t))=10+20e^{-(\text{CD}(t)-{\text{CD}}_{\text{rec}(i)})/2.5},& (7)\end{array}$

as proposed by Moritz et al. (2005). The j^th inter-spike interval ISI_j then is drawn from

$\begin{array}{ll}{\text{ISI}}_j~{N}\left({\text{ISI}}_j^{*},c_{v,i}(CD(t_{ij}^{*}))·{\text{ISI}}_j^{*}\right),& (8)\end{array}$

where the mean $t_{ij}^{*}$ of the time of the next firing event and the corresponding mean inter-spike interval ${\text{ISI}}_j^{*}$ are obtained by solving

$\begin{array}{ll}{t}_{ij}^{*}-t_{i(j-1)}=\frac{1}{{\text{λ}}_i(\text{CD}(t_{ij}^{*}))}={\text{ISI}}_j^{*}& (9)\end{array}$

for $t_{ij}^{*}$ as proposed by Fuglevand et al. (1993). The time t_ij of the j^th firing event is finally calculated as

$\begin{array}{ll}{t}_{ij}=t_{i(j-1)}+{\text{ISI}}_j.& (10)\end{array}$

Note that several distributions other than the normal distribution have been proposed for modeling the distribution of the inter-spike intervals (Jones et al., 2002; Barry et al., 2007). For reasons of simplicity, and as the influence of the shape of the distribution of the inter-spike intervals on the overall EMG and force signals found by Barry et al. (2007) appeared to be rather negligible, a normal distribution is used in the simulation described in section 3.

Equation (7) represents a purely phenomenological model of force variability that reproduces the experimentally observed decrease in coefficient of force variation with increasing activation (Moritz et al., 2005). This phenomenon appears to be the result of different physiological processes, with synaptic noise being the main contributor to the coefficient of variation at low activation levels and oscillations in the descending neural drive explaining most of the variation at higher levels of activation (Dideriksen et al., 2012). Since the model described in this article is purely feed-forward, i.e., no physiological feedback or force control are considered, a phenomenological model of force variability is a reasonable choice. If, however, the model presented here was to be included in a more complete model including force control based on sensory feedback, one should strive to reproduce the experimentally observed variations in force variability through synaptic noise and the intrinsic properties of a force feedback controller such as the one recently proposed by Dideriksen et al. (2017).

One aspect of firing instant calculations that has to the authors' knowledge not been addressed before is the fact that firing instants will occur at arbitrary times between two sampled instants. Simulation studies so far have assumed the set of firing instants t_ij to be a subset of the set of sampling instants (Fuglevand et al., 1993), which is a simplifying assumption that may introduce an artificial distortion in the simulated signal. In our model, we hence allow firing instants to occur at arbitrary times, and not just at sampling instants. Note that for implementing this, it is necessary to calculate differently shifted versions of both the force and EMG twitches of each MU at each firing instant. For the force twitch model presented in section 2.7 an analytical functional relationship is available, for which a time shift can be realized by simply evaluating this function at different input values. Regarding the EMG twitch model described in section 2.6, a description of the Fourier spectrum of the EMG twitch is provided which can be evaluated numerically, and which is then inversely transformed to obtain the EMG twitch in time. Here, a time shift can be realized efficiently in the Fourier domain before applying the inverse transform.

2.5. Intracellular Action Potential Propagation

The propagation of an intracellular action potential (IAP) from the neuromuscular junction (NMJ) of a muscle fiber along both directions toward the two fiber ends can be modeled by representing the actively firing fiber by a distributed current source and sink. In the model initially proposed by Dimitrov and Dimitrova (1998), this distributed fiber membrane current source $\widehat{ı}(z,t)$ is composed of two propagating wavefronts and localized contributions at the NMJ and the two fiber ends. These localized contributions model the IAP generation and extinction process. In the formulation of Farina and Merletti (2001), the model reads

$\begin{array}{ll}\widehat{ı}(z,t)=\frac{\text{d}}{\text{d}z}[\psi (z-z_i-vt)p_1(z)-\psi (-z+z_i-vt)p_2(z)].& (11)\end{array}$

Here, z denotes the spatial variable along the muscle fiber, z_i the location of the NMJ, v the IAP's propagation velocity, and

$\begin{array}{l}{p}_1\left(z\right)=H\left(z-z_i\right)-H\left(z-\left(z_i+L_1\right)\right)\text{ and}\\ p_2\left(z\right)=H\left(z-\left(z_i-L_2\right)\right)-H\left(z-z_i\right)& (12)\end{array}$

the characteristic functions of the two fiber halves, where L₁ and L₂ are the distances between the innervation zone and the right and left tendon, respectively. The propagation velocity v, which has been found to depend on the MU's firing rate (Nishizono et al., 1989) among other factors, is assumed constant throughout the whole simulation for each MU in this study to simplify calculations. For an efficient method to simulate different propagation velocities, refer to Dideriksen et al. (2010b). Moreover,

$\begin{array}{ll}\psi (z)=\frac{\text{d}}{\text{d}z}{V}_m(-z),& (13)\end{array}$

denotes the voltage gradient across the fiber membrane along the fiber axis, where the function V_m(z) prescribes a model for the trans-fiber membrane voltage wave shape and can be chosen arbitrarily to match simulated or measured data. For details on the significance of V_m(z), refer to Plonsey and Barr (2007). Here, the analytical model function

$\begin{array}{ll}{V}_m(z[\text{mm}])=\{\begin{array}{ll}{D}_1{z}^3{e}^{-z}+D_2,\hfill & \text{if }z>0\hfill \\ D_2,\hfill & \text{otherwise}\hfill \end{array}& (14)\end{array}$

with D₁ = 96 mV mm⁻³ and D₂ = −90 mV will be used, as originally proposed by Rosenfalck (1969) and as has been done in previous studies (Farina and Merletti, 2001; Carriou et al., 2016a).

The IAP model in Equation (11) can be shown to be equivalent to choosing

$\begin{array}{l}\widehat{ı}(z,t)=\text{GEN}(t)\delta \left(z-z_i\right)+{\psi }^{\prime }\left(z-z_i-vt\right)p_1\left(z\right)\\ +{\text{EOF}}_1(t)\delta \left(z-z_i-L_1\right)+{\psi }^{\prime }\left(-z+z_i-vt\right)p_2\left(z\right)\\ +{\text{EOF}}_2(t)\delta \left(z-z_i+L_2\right),& (15)\end{array}$

with the Dirac distribution δ, the end-of-fiber components

$\begin{array}{ll}{\text{EOF}}_1(t)=-\psi \left(L_1-vt\right)& (16)\end{array}$

and

$\begin{array}{ll}{\text{EOF}}_2(t)=-\psi \left(L_2-vt\right),& (17)\end{array}$

and the potential generation component

$\begin{array}{ll}\text{GEN}(t)=2\psi \left(-vt\right).& (18)\end{array}$

Here, again, the end-of-fiber components describe the IAP extinction process at the fiber ends, and the potential generation component models the influence of IAP generation at the innervation zone on the membrane current.

This formulation renders—to the authors' opinion—the structure of the model more clear, by explicitly distinguishing between propagating and non-propagating signal components, and by revealing the non-smoothness of the resulting distributed current source, the latter following from the presence of the stationary Dirac distributions at the two fiber ends and the location of the innervation zone. The equivalence of the two formulations of the model is summarized in the following lemma, the proof of which is given in the Appendix.

Lemma 1. The expressions given in Equations (11) and (15) to (18) are equivalent, assuming that ψ ∈ ${{C}}^{\infty }$¹.

In order to achieve a smooth current source model as opposed to the discontinuous model proposed by Farina and Merletti (2001) and analyzed above, Carriou et al. (2016a), inspired by earlier works by Merletti et al. (1999), have suggested the use of smooth Tukey window functions w_1/2(z) instead of the characteristic functions p_1/2(z), which results in the progressive generation and extinction of action potentials. In their study, they set the window smoothness parameter to α = 0.1, where α = 0 would result in the characteristic functions specified in Equation (12). However, while appearing physiologically reasonable, the influence of this modification on the resulting EMG signal and the EMG-force relationship has been found to be limited (Carriou et al., 2016a). Note that a formulation similar to Equations (15)– (18) can easily be obtained if smooth window functions are used instead of the characteristic functions in Equation (11), since in this case all derivatives are defined in the classical sense [assuming smoothness of ψ(z)].

Interestingly, one can show that both IAP models ensure

$\begin{array}{ll}\underset{-\infty }{\overset{\infty }{{{\displaystyle \int }}^{\text{}}}}\widehat{ı}(z,t)\text{ d}z=0\forall t,& (19)\end{array}$

which means that the sum of all incoming and outgoing currents along each fiber is zero at all times. This property is well motivated by physiology, considering that the invoked electrodynamical processes can be regarded as quasi-static (Plonsey and Barr, 2007), and it is also in accordance with the predictions of the standard Hodgkin-Huxley model for action potential propagation (Kleinpenning et al., 1990). For the discontinuous model, it is precisely a result of the presence of the three Dirac distributions in Equation (15), as these have the combined effect of collecting all remaining currents exerted by the intermediate fiber sections due to their higher potential. This result is the subject of the following lemma, the proof of which again is deferred to the Appendix.

Lemma 2. For compactly supported ψ(z), the IAP model given in Equations (15) to (18) [or equivalently, Equation (11)] yields a formulation of $\widehat{ı}(z,t)$ that satisfies the condition (19)². This also holds true if the characteristic functions p₁(z), p₂(z) are replaced by smooth window functions w₁(z), w₂(z) which fulfill

$\begin{array}{ll}{w}_1(z_i+L_1)=w_1(z_i)=w_2(z_i)=w_2(z_i-L_2)=0& (20)\end{array}$

and which yield a product ψ(z)w_1/2(z) that is differentiable with an integrable derivative.

The Tukey window function employed in Carriou et al. (2016a) and Al Harrach et al. (2017) falls into the class of window functions supported by the above lemma.

2.6. EMG Measurements

Biological tissues can be considered volume conductors (Plonsey and Barr, 2007). The existence of an electric field implies the existence of electric currents traveling through the tissue, and vice versa³. Due to the comparably low rate of change of physiological systems, it is justified (Plonsey and Barr, 2007) to assume these time-varying electric fields to behave as if they were static at each instant of time, whence they are called quasi-static. This assumption amounts to a neglection of the capacitive properties of the tissues. Accordingly, as for static fields, the electric field in a physiological volume conductor is considered equal to the negative gradient of a scalar potential φ, namely,

$\begin{array}{ll}\overrightarrow{E}=-\nabla \phi .& (21)\end{array}$

By Ohm's law, the current density (current per unit of cross-sectional area) in a volume conductor is proportional to the electric field, that is,

$\begin{array}{ll}\overrightarrow{J}=\sigma \overrightarrow{E}=-\sigma \nabla \phi ,& (22)\end{array}$

where σ denotes the conductivity of the medium. Defining a distributed current density source I throughout the region of interest, the divergence of the current density is constrained by

$\begin{array}{ll}\nabla ·\overrightarrow{J}=I.& (23)\end{array}$

Combining Equations (22) and (23) and assuming a homogeneous, isotropic medium yields Poisson's equation for the diffusion of the potential, namely,

$\begin{array}{ll}\Delta \phi =-\frac{I}{\sigma }.& (24)\end{array}$

In the following, the electric field generated by point sources in planar tissue layers will be considered as a model for flat and large muscles, such as the recti abdominis simulated in section 3. The muscle layer is assumed to be infinitely extended and planar, and to be covered by an infinitely extended planar layer of fat and an infinitely extended planar layer of skin. Muscle tissue is considered anisotropic to reflect the difference in conductivity between currents along the muscle fiber axis and currents across the muscle fiber axis, whereas fat and skin tissue are considered isotropic. Muscle fibers are assumed to run along the z direction, with the x and z dimensions spanning the skin plane, and the y dimension being orthogonal to the skin plane, positive vectors pointing outwards.

The geometrical set-up described above has been analyzed by Farina and Rainoldi (1999). For a point source of strength Î located at (0, y₀, 0), the authors derive the 2-D spatial Fourier transform of the resulting potential distribution at the skin surface to be

$\begin{array}{l}\Phi (\widehat{I},{\omega }_x,{\omega }_z;y_0)=\frac{2\widehat{I}}{{\sigma }_{m,p}}{e}^{-{\omega }_{ya}|y_0|}\\ ·\frac{1}{(1+r_c)\cosh \left({\omega }_y^{+}\right)\nu \left({\omega }_y^{+}\right)+(1-r_c)\cosh \left({\omega }_y^{-}\right)\nu \left({\omega }_y^{-}\right)},& (25)\end{array}$

with the abbreviations

$\begin{array}{l}{\omega }_y^{+}={\omega }_y(d_f+d_s),{\omega }_y^{-}={\omega }_y(d_f-d_s),\\ {\omega }_y=\sqrt{{\omega }_x^2+{\omega }_z^2},{\omega }_{ya}=\sqrt{{\omega }_x^2+r_a{\omega }_z^2}& (26)\end{array}$

and

$\begin{array}{ll}\text{ν}(s)={\omega }_{ya}+s{r}_m\tanh (s),& (27)\end{array}$

where ω_x = 2πf_x and ω_z = 2πf_z denote the spatial angular frequencies in the x and z directions, respectively. The coefficients

$\begin{array}{ll}{r}_c=\frac{{\sigma }_s}{{\sigma }_f},r_m=\frac{{\sigma }_f}{{\sigma }_{m,p}},\text{ and }{r}_a=\frac{{\sigma }_{m,f}}{{\sigma }_{m,p}}& (28)\end{array}$

specify ratios of the different tissue conductivities. Finally, y₀ denotes the depth of the point source in the muscular tissue, d_f the thickness of the fat layer and d_s the thickness of the skin layer.

Equation (25) directly yields an analytic description of the 2-D spatial transfer function of the volume conductor via

$\begin{array}{ll}{H}_{\text{vc}}({\omega }_x,{\omega }_z;y_0)=\frac{1}{\widehat{I}}·\text{Φ}(\widehat{I},{\omega }_x,{\omega }_z;y_0).& (29)\end{array}$

EMG measurements are usually taken differentially between a set of electrodes. Consider a regular grid of R × S electrodes with interelectrode distances d_x and d_z, respectively, where R = R_a + R_b + 1 and S = S_a + S_b + 1. The variables subscripted by a and b denote the number of electrodes on the two sides of an arbitrarily chosen reference electrode. The grid is assumed to be aligned parallel to the z axis. Assigning weights ζ_kℓ to the electrodes and assuming all electrodes to attain the same transfer function, the (spatial) transfer function from a given surface potential distribution to the potential measured by such an electrode configuration at each point on the surface is given by (Farina and Merletti, 2001)

$\begin{array}{ll}{H}_{\text{ec}}\left({\omega }_x,{\omega }_z\right)={\displaystyle \sum _{k=-R_a}^{R_b}\sum _{\text{ℓ}=-S_a}^{S_b}}{\zeta }_{k\text{ℓ}}{e}^{-j{\omega }_{x}k{d}_x}{e}^{-j{\omega }_z\text{ℓ}{d}_z}.& (30)\end{array}$

For the transfer function H_ele of a single electrode, various model assumptions can be made as well. For examples and details, please refer to, e.g., Merletti and Parker (2004).

Concatenating the spatial transfer functions H_vc of the volume conductor, H_ec of the electrode configuration and H_ele of the electrodes themselves, the global transfer function of the combined system emerges as

$\begin{array}{ll}{H}_{\text{glo}}\left({\omega }_x,{\omega }_z;y\right)=H_{\text{vc}}\left({\omega }_x,{\omega }_z;y\right)\cdot H_{\text{ele}}\left({\omega }_x,{\omega }_z\right)\cdot H_{\text{ec}}\left({\omega }_x,{\omega }_z\right).& (31)\end{array}$

From this, the 2-D potential distribution on the skin surface as measured by electrode configurations consisting of electrodes with transfer function H_ec and H_ele, respectively, and positioned at (x, z), can generally be calculated as

$\begin{array}{l}\phi (x,z,t)={\displaystyle \underset{ℝ}{\int }\left(i\left(x,y,z,t\right)\underset{(x,z)}{\ast }{h}_{\text{glo}}\left(x,z;y\right)\right)\text{ d}y}\\ ={\displaystyle \underset{ℝ}{\int }{{F}}_{xz}^{-1}}\left\{i\left({\omega }_x,y,{\omega }_z,t\right)\cdot H_{\text{glo}}\left({\omega }_x,{\omega }_z;y\right)\right\}\text{ d}y& (32)\end{array}$

where i(ω_x, y, ω_z, t) = ${{F}}_{xz}$ {i(x, y, z, t)} is the 2-D Fourier transform of the current density source i(x, y, z, t), and *_{(x, z)} denotes 2-dimensional convolution in the x and z variables. For a particular electrode (configuration) location on the skin surface and a muscle fiber following a straight line parallel to the skin surface, Equation (32) simplifies, and the resulting single-fiber action potential SFAP(t) = φ(t) can be calculated numerically (Farina and Merletti, 2001; Petersen, 2015; Carriou et al., 2016a). One can prove that in this case the integration kernel only has removable singularities, which ensures the convergence of a numerical integration scheme (Petersen, 2015). Figure 3 shows exemplary SFAPs resulting from the evaluation of Equation (32) for such fibers using nested numerical integration schemes. Note that while the above derivation has been performed for the case of planar volume conductors, Farina et al. (2004) have derived a similar model for cylindrical volume conductors, which is much more appropriate for the simulation of limb muscles.

FIGURE 3

Figure 3. Simulated SFAPs evoked by a single firing muscle fiber as detected by three surface electrodes positioned close to the NMJ (solid), in between the NMJ and the fiber end (short dashes) and above the fiber end (long dashes). The constant k_EMG represents an arbitrary, static EMG scaling factor.

The numerical solution of the (simplified version of) Equation (32) is computationally moderately expensive (Carriou et al., 2016a), refer to section 3 for a discussion of computation times. Fortunately, this only has to be done once for each fiber before the actual simulation, to calculate the SFAPs of all fibers. During the simulation, multiple shifted versions of these SFAPs are then superposed to generate the actual EMG measurement

$\begin{array}{ll}\text{EMG}\left(t\right)={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{N_i}{\text{MUAP}}_i}}\left(t-t_{ij}\right),& (33)\end{array}$

where n denotes the number of MUs, N_i the number of firing events of that MU, and

$\begin{array}{ll}{\text{MUAP}}_i\left(t\right)={\displaystyle \sum _{j=1}^{{\eta }_i}{\text{SFAP}}_{ij}}\left(t\right)& (34)\end{array}$

the motor unit action potential, which is obtained by summing over the SFAPs of all muscle fibers j belonging to MU i. For a discussion of how to evaluate Equation (32) efficiently especially in the HD-sEMG case, refer to Carriou et al. (2016a). The simulated SFAPs correctly reproduce the dependency of the SFAP shape on the relative position of the recording electrode and the depth of the muscle fiber, as well as the distinction between propagating and localized signal components at the NMJ and the two fiber ends. In particular, the SFAPs display the experimentally observed end-of-fiber extinction components.

2.7. Force Twitches

Each MUAP generates a corresponding force contribution, denoted as a force twitch f_i(t). Most previously proposed muscle models (Fuglevand et al., 1993; Dideriksen et al., 2010a) have employed the force twitch parameterization of Milner-Brown et al. (1973), in which a force twitch is entirely described by its twitch rise time T_ri and its peak twitch force P. This model, however, does not allow setting the half relaxation time T_hr independently of T_ri, which is essential for modeling different muscle fiber types (refer to Figure 1.2 in Merletti and Parker, 2004, for example). For this reason, we propose the use of a different model. Raikova and Aladjov (2002) introduced an additional degree of freedom to obtain the following force twitch model of increased expressiveness:

$\begin{array}{ll}{f}_i\left(t\right)=\{\begin{array}{ll}0,\hfill & t<T_{\text{EMD},i},\hfill \\ p_i\cdot {\left(t-T_{\text{EMD},i}\right)}^m\cdot e^{-\kappa \left(t-T_{\text{EMD},i}\right)},\hfill & t\ge T_{\text{EMD},i}\hfill \end{array}& (35)\end{array}$

with

$\begin{array}{ll}{p}_i=P_i\cdot e^{-\kappa T_{\text{ri,}i}\left(\log T_{\text{ri,}i}-1\right)},& (36)\end{array}$

$\begin{array}{ll}m=\kappa \cdot T_{\text{ri},i},& (37)\end{array}$

and

$\begin{array}{ll}\kappa =\frac{\log 2}{T_{\text{hr,}i}-T_{\text{ri},i}\log \left(T_{\text{hr},i}/T_{\text{ri},i}\right)-T_{\text{ri},i}}.& (38)\end{array}$

This model satisfies f_i(T_ri,i) = P_i, and f_i(T_hr,i) = P_i/2. The parameter T_EMD denotes the electromechanical delay between the onset of electrical and mechanical activity of the motor unit. The three parameters T_ri, T_hr and T_EMD are sampled from a Weibull distribution for each MU, as proposed (for the former two) by Contessa and De Luca (2013). Figure 4 shows some exemplary force twitches generated by this model. Note that fixing m = 1 and freely selecting p_i and κ results in the force model of Milner-Brown et al. (1973).

FIGURE 4

Figure 4. Force twitches generated by every tenth MU of the muscle simulated in section 3. These are the same MUs for which the rate coding characteristics are shown in Figure 2. The constant k_f represents an arbitrary, static force scaling factor.

Kernell et al. (1983) have found experimentally that there is a nonlinear relationship between generated isometric muscle force and firing rates of MUs: At high firing rates, the force twitch amplitude decreases. The proposed model includes this nonlinear relationship by scaling individual force twitches in the impulse train of a particular MU by a factor

$\begin{array}{ll}{g}_{ij}\left({\widetilde{\text{λ}}}_{ij}\right)=\{\begin{array}{ll}1\hfill & 0\le {\widetilde{\text{λ}}}_{ij}\le 0.4\hfill \\ \frac{0.4}{{\widetilde{\text{λ}}}_{ij}\left(1-{\gamma }_2\right)}\left(1-{\gamma }_2{e}^{\frac{\left(0.4-{\widetilde{\text{λ}}}_{ij}\right)}{{\gamma }_1}}\right)\hfill & {\widetilde{\text{λ}}}_{ij}>0.4,\hfill \end{array}& (39)\end{array}$

with γ₂ and γ₁ constant muscle parameters, as proposed by Contessa and De Luca (2013). Here, g_ij denotes the gain factor assigned to the jth firing of MU i, and ${\tilde{\lambda }}_{ij}$ is the normalized instantaneous firing rate at that firing event:

$\begin{array}{ll}{\widetilde{\text{λ}}}_{ij}=T_{\text{ri},i}/{\text{ISI}}_j.& (40)\end{array}$

Finally, the total force generated by a muscle is calculated as the superposition

$\begin{array}{ll}F\left(t\right)={\displaystyle \sum _{i=1}^n{F}_i}\left(t\right)={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^{N_i}{f}_i}}\left(t-t_{ij}\right)& (41)\end{array}$

of the individual force twitches of all MUs, with n the number of MUs in the muscle, F_i(t) the force contribution of MU i over time, N_i the number of firing events of MU i, and t_ij the j^th firing instant of MU i, calculated following Equation (10).

2.8. Excitation-Force Relationship and Physiological Force Control

Motor unit firing rate models such as the model of De Luca and Hostage (2010) described in section 2.3 usually define the firing rate as a function of the desired normalized muscle force ${\tilde{F}}_{\text{ref}}=F_{\text{ref}}/F_{\text{max}}$, i.e., they choose the common drive input as $\text{CD}\equiv {\tilde{F}}_{\text{ref}}$. This choice is a consequence of the fact that these models are derived from experimental measurements of MU firing rates at different muscle force levels. Defining the firing rates λ_i of all MUs (e.g., as described in section 2.3) as a function of the force target ${\tilde{F}}_{\text{ref}}$ and also defining the force generating properties of these MUs (e.g., as described in section 2.7) uniquely determines the generated muscle force F, see Figure 1. It is, however, by no means guaranteed that the generated muscle force F will be equal to the force reference F_ref that has been used as an input to the firing rate model. In other words, the generated muscle force output does not match the desired muscle force. This discrepancy is a consequence of the fact that the rate coding/recruitment model is not adjusted to the force generating capacity of the simulated muscle (or vice versa). Both have been defined individually, in parts based on physiological measurements, but they need to be adapted to each other. Note that this is a general problem that every EMG-force model needs to solve.

Dideriksen et al. (2010b) and Venugopal et al. (2017) solved this issue by introducing a simple PID controller for adjusting the model input, such that the error between desired and actual force output is minimized. This method corresponds to the introduction of an artificial feedback loop in Figure 1. In their model, parameters change over time (simulating fatigue), rendering the introduction of such a feedback loop an elegant solution to the problem of time-varying input-output consistency. Contessa and De Luca (2013) introduced a similar force feedback loop, although they implemented the feedback in an offline instead of online fashion, updating the input signal in hindsight and re-running the simulation if the force output deviated too strongly from the desired output. While all three articles (Dideriksen et al., 2010b; Contessa and De Luca, 2013; Venugopal et al., 2017) include some remarks on physiological feedback processes, neither model was meant to replicate properties of actual physiological feedback control, but rather to account for the fact that the rate coding and force generation model components had not been adjusted to each other. One significant drawback of this method is that the feedback loop unpredictably distorts the characteristics of the assumed rate coding model, effectively leading to a different model being used in simulation than the one that has been initially described. Finally, note that many researchers have worked on understanding and modeling actual physiological feedback control (Wolpert and Ghahramani, 2000; Todorov and Jordan, 2002; Dideriksen et al., 2017).

While there certainly is a feedback element in physiological force control, we consider it reasonable that this feedback is mainly required to compensate for external disturbances and changing muscle properties, not to account for a static input-output inconsistency of the neuromuscular system. For these reasons, we propose to model physiological force control as the combined action of two separate components:

1. a feed-forward component that adjusts the rate coding and recruitment model to the force generation model such that in the undisturbed and static case, the muscle's force output equals the desired target force, and

2. a feedback component that processes sensory information to counteract external disturbances and changing muscles properties due to muscular fatigue and changes in muscle geometry.

Figure 5 illustrates the proposed force control scheme, as opposed to the scheme used traditionally (Dideriksen et al., 2010b; Contessa and De Luca, 2013; Venugopal et al., 2017). In this article, we focus on describing a comprehensive feed-forward model of muscular EMG and force generation, and hence we only consider the first of the two components here. Our approach to the modeling of the feed-forward component is based on a simple, static input nonlinearity that is applied equally to all MUs in the muscle.

FIGURE 5

Figure 5. Illustration of two physiological force control schemes. Traditionally, authors have followed approach (A), in which the controller needs to compensate for the nonlinearity of the muscular force generation model in addition to external disturbances and time-varying parameters (Dideriksen et al., 2010b; Contessa and De Luca, 2013; Venugopal et al., 2017). This nonlinearity is not a physiological property but rather results from the fact that the rate coding model is not adjusted to the force generating capacity of the simulated muscle (or vice versa). We propose the use of control scheme (B), in which a static nonlinearity is used to adjust the rate coding model to the force generation model, yielding an adjusted, exactly linearized model of muscular force generation and leaving the feedback controller only with the task of compensating for external disturbances and time-varying model parameters.

The general idea is as follows: denoting the desired normalized muscle force by ${\tilde{F}}_{\text{ref}}\in \left[0,1\right]$, our aim is to choose the common drive $\text{CD}({\tilde{F}}_{\text{ref}})\in \left[0,1\right]$ to the MU pool such that

$\begin{array}{ll}\tilde{F}\left(\text{CD}\left({\tilde{F}}_{\text{ref}}\right)\right)={\tilde{F}}_{\text{ref}}\forall {\tilde{F}}_{\text{ref}}\in [0,1],& (42)\end{array}$

where $\tilde{F}(\text{CD})$ is the normalized generated force output. To achieve this aim, we evaluate $\tilde{F}(\text{CD})$ for all CD ∈ [0, 1], i.e., the normalized generated muscle force as a (nonlinear) function of the common drive. Utilizing the inverse

$\begin{array}{ll}\text{CD}={\tilde{F}}^{-1}\left({\tilde{F}}_{\text{ref}}\right)& (43)\end{array}$

as the common drive input to the firing rate model then yields a simulation model that satisfies condition (42). Next, we will derive an analytical expression for the nonlinear function $\tilde{F}(\text{CD})$ that can be evaluated efficiently.

Averaging over individual firing events, the mean generated muscle force as a function of the common drive is given by

$\begin{array}{ll}\overline{F}\left(\text{CD}\right)={\displaystyle \sum _{i=1}^n}{g}_i\left({\text{λ}}_i\left(\text{CD}\right)\right)\cdot {\Omega }_i·{\text{λ}}_i\left(\text{CD}\right),& (44)\end{array}$

where

$\begin{array}{ll}{\Omega }_i={\displaystyle \underset{0}{\overset{\infty }{\int }}{f}_i}\left(t\right)\text{ d}t=p_i\cdot {\kappa }^{-m-1}·\Gamma \left(m+1\right),& (45)\end{array}$

with Γ(x) the Gamma function, denotes the total impulse generated by a single force twitch f_i(t) of MU i, g_i denotes the nonlinear force gain factor defined in Equation (39), and $\bar{x}$ in general denotes the mean over time of a variable x. By evaluating

$\begin{array}{ll}\tilde{\overline{F}}=\frac{1}{{\overline{F}}_{\max }}\cdot {\displaystyle \sum _{i=1}^n}g\left({\text{λ}}_i\left(\text{CD}\right)\right)\cdot {\Omega }_i·{\text{λ}}_i\left(\text{CD}\right)& (46)\end{array}$

for different values of the common drive CD, one can determine $\widetilde{\stackrel{̄}{F}}(\text{CD})$. Employing $\text{CD}={\widetilde{\stackrel{̄}{F}}}^{-1}({\widetilde{\stackrel{̄}{F}}}_{\text{ref}})$ as the common drive input to the firing rate model then yields a simulation model that satisfies

$\begin{array}{ll}\frac{\overline{F}(\text{CD}({\tilde{\overline{F}}}_{\text{ref}}))}{{\overline{F}}_{\max }}={\tilde{F}}_{\text{ref}}\forall {\tilde{\overline{F}}}_{\text{ref}}\in [0,1].& (47)\end{array}$

Condition (47) ensures that while there may be differences between desired and generated muscle force at individual time instants due to the stochastic nature of the firing instants (see section 2.4), the two forces agree on average.

The model described above for the feed-forward component of physiological force control corresponds to the static input nonlinearity also depicted in Figures 1 and 5 determining the common drive input to the rate coding and recruitment model as a function of the target force level. The introduction of this static input nonlinearity represents an adjustment of the firing rate model to the force generation model, since the rate coding model now no longer receives ${\tilde{F}}_{\text{ref}}$ as an input, but rather its nonlinear transformation ${\tilde{F}}^{-1}({\tilde{F}}_{\text{ref}})$. The consequences of this model adjustment will be assessed in the following section, and we will discuss its plausibility. Note that from a control theoretic point of view, the approach presented here represents an exact linearizing transformation of the control plant, as indicated in Figure 5B.

3. Simulation Results

Using the mathematical model presented in the previous section, a numerical simulation of the rectus abdominis muscle has been implemented as a test scenario, using the R programming language (R Core Team, 2017). In total, 300 MUs have been simulated, each consisting of between 30 and 150 muscle fibers, organized into three separate muscle bellies of each of the two recti (left and right). Single differential detection was assumed for the EMG signals, with the two electrodes placed on the linea alba between the second and third belly. Parameter values obtained from physiological measurements of the rectus abdominis have been chosen for those parameters where they were available (Delp et al., 2001; Rankin et al., 2006; Teyhen et al., 2012). A series of constant-force isometric contractions has been simulated at progressively increasing force levels in steps of 5%MVC up to 100%MVC. Each step lasted 4 s, adding to a total simulation time of 80 s at a simulated sampling frequency of 1024 Hz. Using parallel processing on an Intel Core i7-6700K 4.0 GHz quad core processor with 32 GB of memory, the generation of one such muscle and the calculation of the corresponding sEMG twitch shapes took approximately 3.5 h, while the subsequent simulation of the described contraction using this muscle model took (on a single core) approximately 1 d. In order to assess the effect of the inherent randomness of the MU placement algorithm, the muscle generation procedure and the subsequent simulation each have been executed five times. The simulated signals are available online at the Dryad data repository (Petersen and Rostalski, 2019).

Figure 2 shows a comparison of the two discussed rate coding models with and without the static input nonlinearity introduced in section 2.8. It can be observed that the input nonlinearity introduces an increase in the slope of the firing rate curves of all MUs at high activation levels in both models. This is in perfect agreement with the experimental findings of Erim et al. (1996) and also makes sense from an intuitive point of view: Once most MUs are recruited, firing rates must increase faster than before in order to achieve an increase in muscle force output. It hence appears that the introduction of this static input nonlinearity, calculated by considering all model parameters simultaneously, increases the degree of similarity between experimental observations and models of rate coding and recruitment. The non-smoothness of the adjusted rate coding characteristics is a consequence of the same non-smoothness in the original activation-force relationship, due to new MUs being recruited at discrete levels of excitation. Accordingly, the non-smoothness of the adjusted rate coding characteristics is required to obtain a smooth input-output force relationship. Note that for lack of better data, we employed rate coding model parameter values similar to the ones described by De Luca and Hostage (2010) for the vastus lateralis (VL) muscle, values that have been selected to match isometric measurements between 20%MVC and 100%MVC (De Luca and Hostage, 2010). This may be a reason for the relatively high initial firing rates observable in Figure 2, as compared to physiological measurements (Erim et al., 1996). The following further results have all been obtained using the model of De Luca and Hostage (2010), since no experimentally justified parameters values for our newly proposed rate coding model are currently available.

Figure 6 shows the amplitude distribution of the simulated EMG signal, which resembles a smoothed Laplacian distribution. This type of distribution has been reported previously for the amplitude distributions of real measurement signals (Clancy and Hogan, 1999). Figure 7 displays the simulated series of isometric contractions; Figure 8 shows the coefficient of variation and the standard deviation of the simulated force signal at different levels of muscle activation. Both force variability graphs resemble those reported by Barry et al. (2007) for experimental measurements of index finger force steadiness. Finally, Figure 9 shows the sEMG-Force relationship calculated over the simulated series of isometric contractions.

FIGURE 6

Figure 6. Amplitude distribution of the simulated sEMG signal over the course of a series of isometric contractions at progressively increasing activation levels (solid line). Superimposed are a Gaussian (dashed line) and a Laplacian (dotted line) distribution. The constant k_EMG represents an arbitrary, static EMG scaling factor.

FIGURE 7

Figure 7. Simulated force signal of the simulated recti during a series of isometric contractions. The constant k_f represents an arbitrary, static force scaling factor.

FIGURE 8

Figure 8. Coefficient of variation (A) and normalized standard deviation (B) of the simulated force signal of one belly of the simulated recti during a series of isometric contractions as a function of muscle activation. Shown are mean ± standard deviation of these values over two simulation runs (including repeated MU placement).

FIGURE 9

Figure 9. Steady-state sEMG-Force relation of the simulated recti abdominis over the course of the simulated constant-force isometric contractions. The first 250 ms of each contraction were discarded, and mean force and mean absolute value (MAV) of the EMG signal were calculated over the remainder of the contraction. Both values were normalized to their minima and maxima over all contractions. Shown are the means of these values over five simulation runs (including repeated MU placement).

4. Discussion and Conclusion

In this article, a comprehensive, feed-forward model of surface EMG and force generation during voluntary isometric contractions in skeletal muscles has been described. The model consists of a large number of linear PDEs, ODEs, and various static, sometimes stochastic nonlinear relationships, some of which are solved analytically, while others are calculated numerically. Particular emphasis is placed on choosing electrical and mechanical MU properties so as to achieve realistic EMG-force relationships, and physiological justification is provided for every component of the model. The model consolidates and extends several previously presented results regarding the different components of the physiological system, incorporating recent advances in understanding of physiology. In addition to combining previously isolated results in a unified model and rigorously deriving a more intuitive formulation of the EMG model of Farina and Merletti (2001), new models for several subcomponents of the muscular system have been proposed.

Firstly, a novel nonlinear input transformation is proposed which ensures that the generated muscle forces match the target force level in the isometric, undisturbed case, and which adjusts the rate-coding and recruitment characteristics of the model in a physiologically meaningful way to the force generating capabilities of the model. While previous studies (Dideriksen et al., 2010b, 2017; Contessa and De Luca, 2013; Venugopal et al., 2017) have used PID feedback control to solve this problem, we propose to separate the tasks of

1. adjusting the different model components (rate coding/recruitment and force generation) to each other, and

2. compensating for external disturbances and time-varying parameters.

While the second task inevitably must be (and physiologically is) solved using feedback control, the first one can be solved by means of a static nonlinear transformation between the desired output force level and the common drive input to the rate coding model, as proposed here. This solution entails several useful consequences. On the one hand, it allows for an analytic derivation of the rate coding scheme adjusted to the force generating capabilities of the muscle, as illustrated in Figure 2. It can be observed that the adjustment introduces physiologically reasonable modifications to the model, such as an upwards bend in the firing rate curves at high activation levels (Erim et al., 1996). This is especially interesting since the properties of sEMG-force models have been found to depend strongly on the assumed rate coding model and parameters (Keenan and Valero-Cuevas, 2007). On the other hand, from a control theoretic perspective, the introduction of this static nonlinearity represents an exact linearization of the control plant, which significantly simplifies the feedback control task as well as a subsequent analysis of the closed-loop system. Note that in this article, only the feed-forward component of the model has been considered; combining the presented model with an actual feedback controller is an obvious next step (refer to Figure 5 for an illustration of the two discussed force control schemes).

Secondly, the classical Fuglevand model (Fuglevand et al., 1993) has been extended to reflect findings regarding the dependence of MU and muscle fiber properties on MU and fiber size. In particular, not only the recruitment threshold and the peak twitch force are considered size principle parameters (Henneman, 1957; Henneman et al., 1965), but also the electrical conduction velocity, the time to peak force, the single fiber action potential (SFAP) amplitude and the single fiber peak twitch force. By describing the latter two quantities as a function of MU size, we have introduced an implicit measure of fiber size (fiber diameter) that depends on MU size. This model reflects the differences between small and large MUs and the different fiber types, while still maintaining a continuous distribution of MU and fiber properties that has been found to represent physiological findings better than discrete types with different properties (Enoka and Fuglevand, 2001; Heckman and Enoka, 2012).

Thirdly, a new model of individual MU firing rates has been proposed, which combines several desirable features: A steep, initial slope of the characteristics transitions smoothly into a flatter, linear region at higher activation levels. The onion skin principle is followed, i.e., earlier recruited MUs retain higher firing rates at all points. The initial firing rate trend, i.e., the development of the firing rates of each MU at the point of recruitment, can be freely adjusted to be either increasing or decreasing. And, finally, the convergence behavior of the firing rates of all MUs for high activation levels can also be adjusted freely. By combination with the aforementioned nonlinear input transformation, the model also displays the experimentally observed (Erim et al., 1996) increased slopes starting from the point of full recruitment. Currently, this model has not been fit to any physiological data; this is of course an obvious avenue for future work.

The described model of muscular force generation is in no way to be considered fixed: any component of the model can and should be exchanged for other models of that particular physiological subsystem, e.g., to increase or reduce model complexity, or to take future physiological insights into account. If, for example, more complex and realistic geometries are to be taken into account, the analytical, planar volume conductor and action potential propagation model used here can be exchanged against much more detailed models proposed in the literature (Lowery et al., 2004; Mordhorst et al., 2015), while still using the MU pool model, rate coding and recruitment models and force control concepts described in this article. To further increase physiological realism, regionalized MU placement could be modeled, e.g., as proposed by Robertson and Johnston (2017). Currently, effects due to the presence of muscle fatigue are not modeled; these could further be added to the model, e.g., by employing the metabolic model of Dideriksen et al. (2010a) and a force feedback loop (Dideriksen et al., 2010a; Venugopal et al., 2017). For the simulation of limb muscles, the employed model of the volume conductor should be exchanged for the cylindrical model presented by Farina et al. (2004) or a more flexible finite-element model (Lowery et al., 2004; Mordhorst et al., 2015). The restriction to isometric contractions could be resolved by implementing time-varying muscle geometry (and hence time-varying MU-electrode transmission paths), and by considering the force-length and force-velocity characteristics of skeletal muscles (Yamaguchi, 2001). And finally, an important challenge is raised by recent insights into the effect of neuromodulation on the neuromuscular system (Heckman et al., 2009), an effect that has to the authors' knowledge not been taken into account in any computational neuromuscular models so far.

The proposed model may be useful, among others, for the simulation of realistic sEMG and force signals for the validation of signal processing algorithms, for analyzing the sensitivity of the recorded signals to various physiological parameters, as well as to enhance the general understanding of physiology by evaluating the impact of model or parameter modifications on the system behavior. Preliminary versions of this model have already been used successfully for the validation of algorithms for the separation of inspiratory and expiratory activity in sEMG measurements of the respiratory muscles (Buchner et al., 2016; Petersen et al., 2017), as well as for the derivation of a novel algorithm for sEMG-based muscular activity quantification (Olbrich et al., 2018).

Data Availability

The simulated EMG, force, and impulse train signals of all six rectus bellies, as well as the parameter values used to simulate these signals, are available for one simulated rectus abdominis muscle at the Dryad data repository (Petersen and Rostalski, 2019).

Author Contributions

EP and PR contributed to the conception of the work. The resulting mathematical model was derived and implemented by EP under the supervision of PR. EP wrote the first draft of the manuscript. All authors contributed to manuscript revision and approved the submitted version.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Part of the research leading to the presented mathematical model was financially supported by Drägerwerk AG & Co. KGaA, Lübeck, Germany.

Acknowledgments

The authors would like to thank Marcus Eger, Drägerwerk AG & Co. KGaA, for many valuable discussions on the subject of this article, Jakob L. Dideriksen for a brief but inspiring discussion on the subject of physiological feedback control, and the reviewers for many valuable comments that helped to greatly improve this article. We furthermore acknowledge financial support by Land Schleswig-Holstein within the funding programme Open Access Publikationsfonds.

Footnotes

1. ^Note that – strictly speaking – the expression for V_m(z) given in Equation (14) does not result in a compactly supported ψ(z) ∈ ${{C}}^{\infty }$, but Equation (14) could trivially be modified to comply with this requirement, e.g., by convolution with a smooth and compactly supported mollifier. As this is a purely theoretical operation with no practical relevance at all, this has not been pursued here.

2. ^See previous footnote.

3. ^Note that in biological tissues, the charge carriers are ions, as opposed to electrons in electric wires (Plonsey and Barr, 2007, p. 25).

References

Al Harrach, M., Carriou, V., Boudaoud, S., Laforet, J., and Marin, F. (2017). Analysis of the sEMG/force relationship using HD-sEMG technique and data fusion: a simulation study. Comput. Biol. Med. 83, 34–47. doi: 10.1016/j.compbiomed.2017.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Barry, B. K., Pascoe, M. A., Jesunathadas, M., and Enoka, R. M. (2007). Rate coding is compressed but variability is unaltered for motor units in a hand muscle of old adults. J. Neurophysiol. 97, 3206–3218. doi: 10.1152/jn.01280.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchner, H., Petersen, E., Eger, M., and Rostalski, P. (2016). “Convolutive blind source separation on surface EMG signals for respiratory diagnostics and medical ventilation control,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Orlando, FL).

Google Scholar

Buchthal, F., Guld, C., and Rosenfalck, P. (1957). Multielectrode study of the territory of a motor unit. Acta Physiol. 39, 83–104. doi: 10.1111/j.1748-1716.1957.tb01411.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, R. E. (1981) “Motor units: anatomy, physiology, and functional organization,” in Handbook of Physiology: The Nervous System, Motor Control, ed V. B. Brooks (Bethesda, MD: American Physiological Society), 345–422. doi: 10.1002/cphy.cp010210

CrossRef Full Text | Google Scholar

Carriou, V., Boudaoud, S., Laforet, J., and Ayachi, F. S. (2016a). Fast generation model of high density surface EMG signals in a cylindrical conductor volume. Comput. Biol. Med. 74, 54–68. doi: 10.1016/j.compbiomed.2016.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Carriou, V., Laforet, J., Boudaoud, S., and Harrach, M. A. (2016b). “Realistic motor unit placement in a cylindrical HD-sEMG generation model,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Lake Buena Vista, FL: IEEE).

Google Scholar

Clancy, E. A., and Hogan, N. (1999). Probability density of the surface electromyogram and its relation to amplitude detectors. IEEE Trans. Biomed. Eng. 46, 730–739. doi: 10.1109/10.764949

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, F. (2013). Functional Analysis, Calculus of Variations and Optimal Control. Graduate Texts in Mathematics. London: Springer.

Google Scholar

Contessa, P., and De Luca, C. J. (2013). Neural control of muscle force-indications from a simulation model. J. Neurophysiol. 109, 1548–1570. doi: 10.1152/jn.00237.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., and Contessa, P. (2012). Hierarchical control of motor units in voluntary contractions. J. Neurophysiol. 107, 178–195. doi: 10.1152/jn.00961.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., and Erim, Z. (1994). Common drive of motor units in regulation of muscle force. Trends Neurosci. 17, 299–305. doi: 10.1016/0166-2236(94)90064-7

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., and Hostage, E. C. (2010). Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J. Neurophysiol. 104, 1034–1046. doi: 10.1152/jn.01018.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., Hostage, E. C. (2012). Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J. Neurophysiol. 107:1544. doi: 10.1152/jn.z9k-1251-corr.2011. Corrigendum.

CrossRef Full Text

De Luca, C. J., and Nawab, S. H. (2011). Reply to farina and enoka: The reconstruct-and-test approach is the most appropriate validation for surface EMG signal decomposition to date. J. Neurophysiol. 105, 983–984. doi: 10.1152/jn.01060.2010

CrossRef Full Text | Google Scholar

Del Vecchio, A., Negro, F., Felici, F., and Farina, D. (2017). Associations between motor unit action potential parameters and surface EMG features. J. Appl. Physiol. 123, 835–843. doi: 10.1152/japplphysiol.00482.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Delp, S. L., Suryanarayanan, S., Murray, W. M., Uhlir, J., and Triolo, R. J. (2001). Architecture of the rectus abdominis, quadratus lumborum, and erector spinae. J. Biomech. 34, 371–375. doi: 10.1016/S0021-9290(00)00202-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dideriksen, J. L., Farina, D., Baekgaard, M., and Enoka, R. M. (2010a). An integrative model of motor unit activity during sustained submaximal contractions. J. Appl. Physiol. 108, 1550–1562. doi: 10.1152/japplphysiol.01017.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Dideriksen, J. L., Farina, D., and Enoka, R. M. (2010b). Influence of fatigue on the simulated relation between the amplitude of the surface electromyogram and muscle force. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 368, 2765–2781. doi: 10.1098/rsta.2010.0094

PubMed Abstract | CrossRef Full Text | Google Scholar

Dideriksen, J. L., Feeney, D. F., Almuklass, A. M., and Enoka, R. M. (2017). Control of force during rapid visuomotor force-matching tasks can be described by discrete time PID control algorithms. Exp. Brain Res. 235, 2561–2573. doi: 10.1007/s00221-017-4995-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Dideriksen, J. L., Negro, F., Enoka, R. M., and Farina, D. (2012). Motor unit recruitment strategies and muscle properties determine the influence of synaptic noise on force steadiness. J. Neurophysiol. 107, 3357–3369. doi: 10.1152/jn.00938.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimitrov, G. V., and Dimitrova, N. A. (1998). Precise and fast calculation of the motor unit potentials detected by a point and rectangular plate electrode. Med. Eng. Phys. 20, 374–381. doi: 10.1016/S1350-4533(09)00014-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Elder, G. C., Bradbury, K., and Roberts, R. (1982). Variability of fiber type distributions within human muscles. J. Appl. Physiol. 53, 1473–1480. doi: 10.1152/jappl.1982.53.6.1473

PubMed Abstract | CrossRef Full Text | Google Scholar

Enoka, R. M., and Fuglevand, A. J. (2001). Motor unit physiology: Some unresolved issues. Muscle Nerve 42, 4–17. doi: 10.1002/1097-4598(200101)24:1<4::AID-MUS13>3.0.CO;2-F

CrossRef Full Text | Google Scholar

Erim, Z., De Luca, C. J., Mineo, K., and Aoki, T. (1996). Rank-ordered regulation of motor units. Muscle Nerve 19, 563–573.

PubMed Abstract | Google Scholar

Farina, D., and Enoka, R. M. (2011). Surface EMG decomposition requires and appropriate validation. J. Neurophysiol. 105, 981–982; author reply 983-4. doi: 10.1152/jn.00855.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., and Merletti, R. (2001). A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Trans. Biomed. Eng. 48, 637–646. doi: 10.1109/10.923782

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Merletti, R., and Enoka, R. M. (2015). Reply to De Luca, Nawab, and Kline: the proposed method to validate surface EMG signal decomposition remains problematic. J. Appl. Physiol. 118:1085. doi: 10.1152/japplphysiol.00107.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Mesin, L., Martina, S., and Merletti, R. (2004). A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Trans. Biomed. Eng. 51, 415–426. doi: 10.1109/TBME.2003.820998

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., and Rainoldi, A. (1999). Compensation of the effect of sub-cutaneous tissue layers on surface EMG: a simulation study. Med. Eng. Phys. 21, 487–497. doi: 10.1016/S1350-4533(99)00075-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleshman, J. W., Munson, J. B., Sypert, G. W., and Friedman, W. A. (1981). Rheobase, input resistance, and motor-unit type in medial gastrocnemius motoneurons in the cat. J. Neurophysiol. 46, 1326–1338. doi: 10.1152/jn.1981.46.6.1326

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuglevand, A. J., Winter, D. A., and Patla, A. E. (1993). Models of recruitment and rate coding organization in motor-unit pools. J. Neurophysiol. 70, 2470–2488. doi: 10.1152/jn.1993.70.6.2470

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuglevand, A. J., Winter, D. A., Patla, A. E., and Stashuk, D. (1992). Detection of motor unit action potentials with surface electrodes: influence of electrode size and spacing. Biol. Cybern. 67, 143–153. doi: 10.1007/BF00201021

PubMed Abstract | CrossRef Full Text | Google Scholar

Hakansson, C. H. (1956). Conduction velocity and amplitude of the action potential as related to circumference in the isolated fibre of frog muscle. Acta Physiol. Scand. 37, 14–34. doi: 10.1111/j.1748-1716.1956.tb01338.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Heckman, C., and Enoka, R. M. (2012). “Motor unit,” in Comprehensive Physiology, ed R. Terjung. Available online at: https://onlinelibrary.wiley.com/doi/abs/10.1002/cphy.c100087

Google Scholar

Heckman, C. J., Mottram, C., Quinlan, K., Theiss, R., and Schuster, J. (2009). Motoneuron excitability: The importance of neuromodulatory inputs. Clin. Neurophysiol. 120, 2040–2054. doi: 10.1016/j.clinph.2009.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Henneman, E. (1957). Relation between size of neurons and their susceptibility to discharge. Science 126, 1345–1347. doi: 10.1126/science.126.3287.1345

PubMed Abstract | CrossRef Full Text | Google Scholar

Henneman, E., Somjen, G., and Carpenter, D. O. (1965). Excitability and inhibitibility of motoneurons of different sizes. J. Neurophysiol. 28, 599–620. doi: 10.1152/jn.1965.28.3.599

PubMed Abstract | CrossRef Full Text | Google Scholar

Hogan, N., and Mann, R. W. (1980). Myoelectric signal processing: Optimal estimation applied to electromyography - part I: derivation of the optimal myoprocessor. IEEE Trans. Biomed. Eng. 27, 382–395.

PubMed Abstract | Google Scholar

Johnson, M. A., Polgar, J., Weightman, D., and Appleton, D. (1973). Data on the distribution of fibre types in thirty-six human muscles – an autopsy study. J. Neurol. Sci. 18, 111–129.

PubMed Abstract | Google Scholar

Jones, K. E., Hamilton, A. F., and Wolpert, D. M. (2002). Sources of signal-dependent noise during isometric force production. J. Neurophysiol. 88, 1533–1544. doi: 10.1152/jn.2002.88.3.1533

PubMed Abstract | CrossRef Full Text | Google Scholar

Keenan, K. G., and Valero-Cuevas, F. J. (2007). Experimentally valid predictions of muscle force and EMG in models of motor-unit function are most sensitive to neural properties. J. Neurophysiol. 98, 1581–1590. doi: 10.1152/jn.00577.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Kernell, D. (2006). The Motoneurone and Its Muscle Fibres. Oxford: Oxford University Press UK.

Google Scholar

Kernell, D., Eerbeek, O., and Verhey, B. A. (1983). Relation between isometric force and stimulus rate in cat's hindlimb motor units of different twitch contraction time. Exp. Brain Res. 50, 220–227.

PubMed Abstract | Google Scholar

Kleinpenning, P. H., Gootzen, T. H., Van Oosterom, A., and Stegeman, D. F. (1990). The equivalent source description representing the extinction of an action potential at a muscle fiber ending. Math. Biosci. 101, 41–61.

PubMed Abstract | Google Scholar

Lo Conte, L. R., Merletti, R., and Sandri, G. V. (1994). Hermite expansions of compact support waveforms: applications to myoelectric signals. IEEE Trans. Biomed. Eng. 41, 1147–1159.

PubMed Abstract | Google Scholar

Lowery, M. M., Stoykov, N. S., Dewald, J. P., and Kuiken, T. A. (2004). Volume conduction in an anatomically based surface EMG model. IEEE Trans. Biomed. Eng. 51, 2138–2147. doi: 10.1109/TBME.2004.836494

PubMed Abstract | CrossRef Full Text | Google Scholar

Lucas, S. M., Ruff, R. L., and Binder, M. (1987). Specific tension measurements in single soleus and medial gastrocnemius muscle fibers of the cat. Exp. Neurol. 95, 142–154.

PubMed Abstract | Google Scholar

McGill, K. C. (2004). Surface electromyogram signal modelling. Med. Biol. Eng. Comput. 42, 446–454. doi: 10.1007/BF02350985

PubMed Abstract | CrossRef Full Text | Google Scholar

Merletti, R., and Farina, D., editors (2016). Surface Electromyography: Physiology, Engineering and Applications. Hoboken, NJ: John Wiley and Sons, Inc.

Google Scholar

Merletti, R., Lo Conte, L., Avignone, E., and Guglielminotti, P. (1999). Modeling of surface myoelectric signals. I. Model implementation. IEEE Trans. Biomed. Eng. 46, 810–820.

PubMed Abstract | Google Scholar

Merletti, R., and Parker, P. A., editors (2004). Electromyography. Hoboken, NJ: John Wiley & Sons, Inc.

Google Scholar

Milner-Brown, H. S., Stein, R. B., and Yemm, R. (1973). Changes in firing rate of human motor units during linearly changing voluntary contractions. J. Physiol. 230, 371–390.

PubMed Abstract | Google Scholar

Mordhorst, M., Heidlauf, T., and Röhrle, O. (2015). Predicting electromyographic signals under realistic conditions using a multiscale chemo-electro-mechanical finite element model. Interface Focus 5:20140076. doi: 10.1098/rsfs.2014.0076

PubMed Abstract | CrossRef Full Text | Google Scholar

Moritz, C. T., Barry, B. K., Pascoe, M. A., and Enoka, R. M. (2005). Discharge rate variability influences the variation in force fluctuations across the working range of a hand muscle. J. Neurophysiol. 93, 2449–2459. doi: 10.1152/jn.01122.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Nandedkar, S. D., and Stålberg, E. (1983). Simulation of single muscle fibre action potentials. Med. Biol. Eng. Comput. 21, 158–165.

PubMed Abstract | Google Scholar

Navallas, J., Malanda, A., Gila, L., Rodríguez, J., and Rodríguez, I. (2010). A muscle architecture model offering control over motor unit fiber density distributions. Med. Biol. Eng. Comput. 48, 875–886. doi: 10.1007/s11517-010-0642-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishizono, H., Kurata, H., and Miyashita, M. (1989). Muscle fiber conduction velocity related to stimulation rate. Electroencephalogr. Clin. Neurophysiol. 72, 529–534.

PubMed Abstract | Google Scholar

Olbrich, M., Petersen, E., Hoffmann, C., and Rostalski, P. (2018). Sparse estimation for the assessment of muscular activity based on sEMG measurements. IFAC-PapersOnLine 51, 305–310. doi: 10.1016/j.ifacol.2018.09.152

CrossRef Full Text | Google Scholar

Petersen, E. (2015). Integrative Mathematical Modeling of Surface Electromyographic Modeling. Master's thesis, University of Hamburg.

Petersen, E. (2016). “A mathematical model of surface electromyographic measurements,” in Proceedings of the Workshop Biosignal Processing (Berlin).

Google Scholar

Petersen, E., Buchner, H., Eger, M., and Rostalski, P. (2017). Convolutive blind source separation of surface EMG measurements of the respiratory muscles. Biomed. Eng. Biomed. Technik 62, 171–181. doi: 10.1515/bmt-2016-0092

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, E., and Rostalski, P. (2019). Data from: A comprehensive mathematical model of motor unit pool organization, surface electromyography and force generation. Dryad Digit. Repository. doi: 10.5061/dryad.326qs26

CrossRef Full Text

Piotrkiewicz, M., and Türker, K. S. (2017). Onion skin or common drive? Front. Cell. Neurosci. 11:2. doi: 10.3389/fncel.2017.00002

CrossRef Full Text | Google Scholar

Plonsey, R., and Barr, R. C. (2007). Bioelectricity, 3rd Edn. New York, NY: Springer.

Google Scholar

Potvin, J. R., and Fuglevand, A. J. (2017). A motor unit-based model of muscle fatigue. PLoS Comput. Biol. 13:e1005581. doi: 10.1371/journal.pcbi.1005581

PubMed Abstract | CrossRef Full Text | Google Scholar

Powers, R. K., and Binder, M. D. (2001). “Input-output functions of mammalian motoneurons,” in Reviews of Physiology, Biochemistry and Pharmacology, eds E. Bamberg, F. M. P. Blaustein, B. R. Greger, F. H. Grunicke, I. R. Jahn, G. W. J. Lederer, B. L. M. MendeU, S. B. A. Miyajima, T. N. Pfanner, F. G. Schultz, B. M. Schweiger (Berlin; Heidelberg: Springer), 137–263.

Google Scholar

R Core Team (2017). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna.

PubMed Abstract

Raikova, R., Krutki, P., Aladjov, H., and Celichowski, J. (2007). Variability of the twitch parameters of the rat medial gastrocnemius motor units - experimental and modeling study. Comput. Biol. Med. 37, 1572–1581. doi: 10.1016/j.compbiomed.2007.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Raikova, R. T., and Aladjov, H. T. (2002). Hierarchical genetic algorithm versus static optimization and investigation of elbow flexion and extension movements. J. Biomech. 35, 1123–1135. doi: 10.1016/S0021-9290(02)00031-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rankin, G., Stokes, M., and Newham, D. J. (2006). Abdominal muscle size and symmetry in normal subjects. Muscle Nerve 34, 320–326. doi: 10.1002/mus.20589

PubMed Abstract | CrossRef Full Text | Google Scholar

Robertson, J. W. (2014). Modifying Muscle Properties in a Leading Neuromuscular Model: The Fuglevand Model Revisited. Ph.D. thesis. Alberta: University of Calgary.

Robertson, J. W., and Johnston, J. A. (2017). Modifying motor unit territory placement in the Fuglevand model. Med. Biol. Eng. Comput. 55, 2015–2025. doi: 10.1007/s11517-017-1645-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Falces, J., Navallas, J., and Malanda, A. (2012). “EMG modeling,” in Computational Intelligence in Electromyography Analysis-A Perspective on Current Applications and Future Challenges, ed G. R. Naik, InTech, 3–36.

Google Scholar

Rosenfalck, P. (1969). Intra and extracellular fields of active nerve and muscle fibers. A physico-mathematical analysis of different models. Acta Physiol. Scand. 321, 1–168.

Sadoyama, T., Masuda, T., Miyata, H., and Katsuta, S. (1988). Fibre conduction velocity and fibre composition in human vastus lateralis. Eur. J. Appl. Physiol. Occupat. Physiol. 57, 767–771.

PubMed Abstract | Google Scholar

Schnetzer, M.-A., Rüegg, D. G., Baltensperger, R., and Gabriel, J.-P. (2001). “Three-dimensional model of a muscle and simulation of its surface emg,” in 23rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Istanbul).

Google Scholar

Sinderby, C., Lindström, L., and Grassino, A. E. (1995). Automatic assessment of electromyogram quality. J. Appl. Physiol. 79, 1803–1815.

PubMed Abstract | Google Scholar

Stegeman, D. F., Blok, J. H., Hermens, H. J., and Roeleveld, K. (2000). Surface EMG models: properties and applications. J. Electromyogr. Kinesiol. 10, 313–326. doi: 10.1016/S1050-6411(00)00023-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Teyhen, D. S., Childs, J. D., Stokes, M. J., Wright, A. C., Dugan, J. L., and George, S. Z. (2012). Abdominal and lumbar multifidus muscle size and symmetry at rest and during contracted states. J. Ultrasound Med. 31, 1099–1110. doi: 10.7863/jum.2012.31.7.1099

PubMed Abstract | CrossRef Full Text | Google Scholar

Todorov, E., and Jordan, M. I. (2002). Optimal feedback control as a theory of motor coordination. Nat. Neurosci. 5, 1226–1235. doi: 10.1038/nn963

PubMed Abstract | CrossRef Full Text | Google Scholar

Totosy de Zepetnek, J. E., Zung, H. V., Erdebil, S., and Gordon, T. (1992). Innervation ratio is an important determinant of force in normal and reinnervated rat tibialis anterior muscles. J. Neurophysiol. 67, 1385–1403. doi: 10.1152/jn.1992.67.5.1385

PubMed Abstract | CrossRef Full Text | Google Scholar

Trontelj, J. V., Jabre, J., and Mihelin, M. (2004). “Needle and wire detection techniques,” in Electromyography, eds R. Merletti and P. A. Parker (Hoboken, NJ: John Wiley & Sons, Inc.), 27–46.

Google Scholar

Van Cutsem, M., Feiereisen, P., Duchateau, J., and Hainaut, K. (1997). Mechanical properties and behaviour of motor units in the tibialis anterior during voluntary contractions. Can. J. Appl. Physiol. 22, 585–597. doi: 10.1139/h97-038

PubMed Abstract | CrossRef Full Text | Google Scholar

Venugopal, G., Deepak, P., Ghosh, D. M., and Ramakrishnan, S. (2017). Generation of synthetic surface electromyography signals under fatigue conditions for varying force inputs using feedback control algorithm. Proc. Inst. Mech. Eng. H 231, 1025–1033. doi: 10.1177/0954411917727307

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., De Stefano, A., and Allen, R. (2006). A simulation model of the surface EMG signal for analysis of muscle activity during the gait cycle. Comput. Biol. Med. 36, 601–618. doi: 10.1016/j.compbiomed.2005.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolpert, D. M., and Ghahramani, Z. (2000). Computational principles of movement neuroscience. Nat. Neurosci. 3(Suppl.):1212–1217. doi: 10.1038/81497

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, G. T. (2001). Dynamic Modeling of Musculoskeletal Motion. 1st Edn. Springer Science+Business Media Dordrecht.

Google Scholar

Zhou, P., and Rymer, W. Z. (2004). Factors governing the form of the relation between muscle force and the EMG: a simulation study. J. Neurophysiol. 92, 2878–2886. doi: 10.1152/jn.00367.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix A: Proofs

Lemma 1. The expressions given in Equations (11) and (15) to (18) are equivalent, assuming that ψ ∈ ${{C}}^{\infty }$.

Proof: Expanding all terms in Equation (11) yields

$\begin{array}{l}\widehat{ı}(z,t)=\frac{\text{d}}{\text{d}z}\left[\psi \left(z-z_i-vt\right)H\left(z-z_i\right)\right]\\ -\frac{\text{d}}{\text{d}z}\left[\psi \left(z-z_i-vt\right)H\left(z-z_i-L_1\right)\right]\\ -\frac{\text{d}}{\text{d}z}\left[\psi \left(-z+z_i-vt\right)H\left(z-z_i+L_2\right)\right]\\ +\frac{\text{d}}{\text{d}z}\left[\psi \left(-z+z_i-vt\right)H\left(z-z_i\right)\right].& (\mathrm{A1})\end{array}$

Note that the derivatives in Equation (A1) can only be understood in the sense of distributions, since the derivative of the Heaviside function H(z) can not be defined in the classical sense at z = 0. For some basic properties of distributions, refer to Appendix B.

Assuming ψ ∈ ${{C}}^{\infty }$ and recalling that 〈(ψH)′, ζ〉 = 〈ψ′H + ψδ, ζ〉 (refer to Appendix B), the duality $\langle \widehat{ı}(t),\zeta \rangle$ can hence be formulated as

$\begin{array}{l}\langle \widehat{ı}(t),\zeta \rangle =\langle {\psi }_{z_i+vt}^{\prime }{H}_{z_i},\zeta \rangle +\langle {\psi }_{z_i+vt}{\delta }_{z_i},\zeta \rangle \\ -\langle {\psi }_{z_i+vt}^{\prime }{H}_{z_i+L_1},\zeta \rangle -\langle {\psi }_{z_i+vt}{\delta }_{z_i+L_1},\zeta \rangle \\ +\langle {\overline{\psi }}_{z_i-vt}^{\prime }{H}_{z_i-L_2},\zeta \rangle -\langle {\overline{\psi }}_{z_i-vt}{\delta }_{z_i-L_2},\zeta \rangle \\ -\langle {\overline{\psi }}_{z_i-vt}^{\prime }{H}_{z_i},\zeta \rangle +\langle {\overline{\psi }}_{z_i-vt}{\delta }_{z_i},\zeta \rangle \\ =\langle {\psi }_{z_i+vt}^{\prime }(H_{z_i}-H_{z_i+L_1})+{\overline{\psi }}_{z_i-vt}^{\prime }(H_{z_i-L_2}-H_{z_i}),\zeta \rangle \\ +\langle ({\psi }_{z_i+vt}+{\overline{\psi }}_{z_i-vt}){\delta }_{z_i},\zeta \rangle \\ -\langle {\psi }_{z_i+vt}{\delta }_{z_i+L_1}+{\overline{\psi }}_{z_i-vt}{\delta }_{z_i-L_2},\zeta \rangle & (\mathrm{A2})\end{array}$

with notations $\overline{f}(x)=f(-x)$, f_a(x) = f(x − a) and ${\overline{f}}_a(x)=f(-x+a)$. This is exactly equivalent to Equations (15) to (18). □

Lemma 2. For compactly supported ψ(z), the IAP model given in Equations (15) to (18) [or equivalently, Equation (11)] yields a formulation of $\widehat{ı}(z,t)$ that satisfies the condition (19). This also holds true if the characteristic functions p_1/2(z) are replaced by smooth window functions w_1/2(z) which fulfill

$\begin{array}{ll}{w}_1\left(z_i+L_1\right)=w_1\left(z_i\right)=w_2\left(z_i\right)=w_2\left(z_i-L_2\right)=0& (\mathrm{A3})\end{array}$

and which yield products ψ(z)w_1/2(z) that are differentiable with an integrable derivative.

Proof: First, we consider the case of $\widehat{ı}(z,t)$ as described by Equations (15) to (18), i.e., using the characteristic functions p₁(z), p₂(z). For compactly supported ψ(z),

$\begin{array}{ll}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\psi }^{\prime }(z)\text{ d}z=0}& (\mathrm{A4})\end{array}$

generally holds. Furthermore,

$\begin{array}{ll}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\text{GEN}(t)\delta \left(z-z_i\right)\text{ d}z=2\psi (-vt)=2}{\displaystyle \underset{-\infty }{\overset{-vt}{\int }}{\psi }^{\prime }(z)\text{ d}z},& (\mathrm{A5})\end{array}$

$\begin{array}{l}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\text{EOF}}_1(t)\delta \left(z-z_i-L_1\right)\text{ d}z=-\psi \left(L_1-vt\right)}\\ =-{\displaystyle \underset{-\infty }{\overset{L_1-vt}{\int }}{\psi }^{\prime }(z)\text{ d}z,}& (\mathrm{A6})\end{array}$

and equivalently for EOF₂. Combining everything yields

$\begin{array}{l}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\widehat{ı}(z,t)\text{ d}z=2}{\displaystyle \underset{-\infty }{\overset{-vt}{\int }}{\psi }^{\prime }(z)\text{ d}z}-{\displaystyle \underset{-\infty }{\overset{L_1-vt}{\int }}{\psi }^{\prime }(z)\text{ d}z}-{\displaystyle \underset{-\infty }{\overset{L_2-vt}{\int }}{\psi }^{\prime }(z)\text{ d}z}\\ +{\displaystyle \underset{z_i}{\overset{z_i+L_1}{\int }}{\psi }^{\prime }\left(z-z_i-vt\right)\text{ d}z}+{\displaystyle \underset{z_i-L_2}{\overset{z_i}{\int }}{\psi }^{\prime }\left(-z+z_i-vt\right)\text{ d}z}\\ =0.& (\mathrm{A7})\end{array}$

Secondly, considering the case of smooth window functions w_1/2(z), we have

$\begin{array}{l}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\widehat{ı}(z,t)\text{ d}z}={\displaystyle \underset{z_i}{\overset{\infty }{\int }}\frac{\text{d}}{\text{ d}z}\left[\psi \left(z-z_i-vt\right)w_1(z)\right]\text{ d}z}\\ -{\displaystyle \underset{-\infty }{\overset{z_i}{\int }}\frac{\text{d}}{\text{ d}z}\left[\psi \left(-z+z_i-vt\right)w_2(z)\right]\text{ d}z}\\ =0& (\mathrm{A8})\end{array}$

by the fundamental theorem of calculus if the conditions on w_1/2(z) stated in the lemma hold. □

Appendix B: Distributions

A distribution T is a continuous linear mapping T: D → ℝ, where D is a given set of so-called test functions. In the following, the value of a distribution T acting on a test function ζ shall be denoted by the duality 〈T, ζ〉, with test functions chosen from the set D of compactly supported, smooth functions. Furthermore, two basic properties of distributions are

$\begin{array}{l}\langle \eta T,\zeta \rangle =\langle T,\eta \zeta \rangle \forall \eta \in {{C}}^{\infty },\zeta \in D\end{array}$

and

$\begin{array}{l}\langle T^{\prime },\zeta \rangle =-\langle T,{\zeta }^{\prime }\rangle \forall \zeta \in D.\end{array}$

Moreover (following from the above), considering η ∈ ${{C}}^{\infty }$,

$\begin{array}{l}\langle {(\eta H)}^{\prime },\zeta \rangle =-\langle \eta H,{\zeta }^{\prime }\rangle =-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}(\eta H)(x){\zeta }^{\prime }(x)\text{ d}x}\\ =-{[\eta \zeta ]}_0^{\infty }+{\displaystyle \underset{0}{\overset{\infty }{\int }}{\eta }^{\prime }\zeta \text{ d}x=\langle \delta ,\eta \zeta \rangle +\langle {\eta }^{\prime }H,\zeta \rangle }\\ =\langle {\eta }^{\prime }H+\eta \delta ,\zeta \rangle .\end{array}$

For more background on distributions (i.e., continuous linear functionals), refer to, e.g., Clarke (2013).