 Research
 Open Access
 Published:
Voluntary EMGtoforce estimation with a multiscale physiological muscle model
BioMedical Engineering OnLine volume 12, Article number: 86 (2013)
Abstract
Background
EMGtoforce estimation based on muscle models, for voluntary contraction has many applications in human motion analysis. The socalled Hill model is recognized as a standard model for this practical use. However, it is a phenomenological model whereby muscle activation, forcelength and forcevelocity properties are considered independently. Perreault reported Hill modeling errors were large for different firing frequencies, level of activation and speed of contraction. It may be due to the lack of coupling between activation and forcevelocity properties. In this paper, we discuss EMGforce estimation with a multiscale physiology based model, which has a link to underlying crossbridge dynamics. Differently from the Hill model, the proposed method provides dual dynamics of recruitment and calcium activation.
Methods
The ankle torque was measured for the plantar flexion along with EMG measurements of the medial gastrocnemius (GAS) and soleus (SOL). In addition to Hill representation of the passive elements, three models of the contractile parts have been compared. Using common EMG signals during isometric contraction in four ablebodied subjects, torque was estimated by the linear Hill model, the nonlinear Hill model and the multiscale physiological model that refers to Huxley theory. The comparison was made in normalized scale versus the case in maximum voluntary contraction.
Results
The estimation results obtained with the multiscale model showed the best performances both in fastshort and slowlong term contraction in randomized tests for all the four subjects. The RMS errors were improved with the nonlinear Hill model compared to linear Hill, however it showed limitations to account for the different speed of contractions. Average error was 16.9% with the linear Hill model, 9.3% with the modified Hill model. In contrast, the error in the multiscale model was 6.1% while maintaining a uniform estimation performance in both fast and slow contractions schemes.
Conclusions
We introduced a novel approach that allows EMGforce estimation based on a multiscale physiology model integrating Hill approach for the passive elements and microscopic crossbridge representations for the contractile element. The experimental evaluation highlights estimation improvements especially a larger range of contraction conditions with integration of the neural activation frequency property and forcevelocity relationship through crossbridge dynamics consideration.
Introduction
Any human movement is produced by muscular and skeletal systems controlled by the nervous system. Mechanism of the human body dynamics has been revealed by many biomechanical researches. The general neuromusculoskeletal model of the whole body and its dynamics computation package have been developed and reported[1]. However, the detailed neuromuscular activation system should still be closely analyzed and modeled from microscopic to macroscopic scales. Macroscopic modeling, such as the Hilltype muscle model, is a phenomenological model. Conversely, microscopic modeling is based on muscle physiology, more specifically the dynamics of actinmyosin, allowing it to display a richer response.
Neuromuscular modeling is important for neuroscience to understand how limb movements are smoothly and effectively controlled[2, 3]. It is also valuable for clinical applications to be used for quantitative analysis on abnormal muscle activation patterns such as spasticity induced by stroke or cerebral palsy[4]. Moreover, neuroprosthetic control such as functional electrical stimulation (FES) for paralyzed muscles[5, 6] can be improved through neuromuscular modeling to optimize the muscle activation control.
Surface electromyography (EMG) is widely used to assess the muscle status in a wide range of clinical fields. The signal generated during a contraction is the summation of the signals of all different recruited motor units (MU) at a given time. Each MU receives information from the central nervous system in the form of action potentials, which are transmitted from the alpha motor neuron to the muscle fibers via the neuromuscular junction. During voluntary contraction, the EMG signal can be seen as the interference pattern of all active MUs, with each MU firing at its own mean frequency with its own timing[7]. The contributions of individual MU cannot be easily recovered from the EMG signal. However, EMG signals still contain precious information about muscle activation and could thus represents the entry of a muscle model. Indeed, the important advantage of EMG usage is that it can account for a subject’s individual activation pattern to estimate muscle force. Moreover, as the EMG signals can be considered for each individual groups of antagonist muscles with their own muscle models, cocontraction can be captured and deeply investigated through the computation of individual force and stiffness. On the contrary, numerical optimization techniques typically do not account for muscular cocontraction.
EMGbased models were already used in many studies to estimate torques around joints[8] or in musculoskeletal models[9, 10]. Most of the muscle models used in such studies are based on phenomenological models derived from Hill’s study[11] and summarized by Zajac[12]. The Hilltype model is attractive in biomechanical simulations because of the computational simplicity and easy to calibrate using experimentally measurable variables. It tends to produce accurate estimation results, however in limited ranges of function.
However, recent work[13] performed the validation of Hill model under functionally wider relevant conditions. The authors concluded that the modeling errors were substantial for different firing frequencies and greatest at low motor unit firing rates that are most relevant in normal movement conditions both in voluntary and FES muscle contractions. They pointed out that this could be due to the Hill model assumption whereby muscle activation, forcelength and forcevelocity properties are considered independently. In their discussion, it was suggested that more physiological coupling between activation and forcevelocity properties can be demonstrated in crossbridge models incorporating dependence between physiologically based activation and crossbridge attachment[13]. The important finding is that Hill modeling errors were large for different firing frequencies and greatest at low motor unit firing rates. In[13], classical linear Hill model was used. Especially for solving the problem of the mismatch at low muscle activation level, a nonlinear conversion was proposed[14] and used by other researchers[9]. This can solve partially the problem, however, there is a still the problem to correspond to the estimations in different contraction speeds. Besides, in Hilltype modeling, the cut off frequency should be carefully chosen depending on the type of task because the envelope of the estimated force is almost defined by the lowpass filtering of the activation. Thus, we aim at proposing multiscale muscle model, which can have more physiological coupling between activation and forcevelocity properties by incorporating microscopic crossbride and dual dynamics pathway for the recruitment / activation level and the calcium triggering the contraction cycle. However, the whole passive structure is based on Hill model as only the contractile element differs.
The first microscopic muscle model was proposed by Huxley[15]. He detailed the interaction of the crossbridges in a sarcomere in order to explain the force generation. The distinctions between microscopic and macroscopic levels are not absolute; thus a sarcomere model can be used to represent a whole muscle, which is assumed to be a homogeneous assembly of identical sarcomeres. Conversely, the Hill model, which was originally developed for whole muscles, has been used to represent the dynamics of individual sarcomeres within a fiber[16]. The distributionmoment approach of Zahalak[17] is a model for sarcomeres or whole muscle, which is extracted via a formal mathematical Gaussian approximation from Huxley crossbridge model. This model fills a gap between microscopic and macroscopic levels. Based on Huxley and Hilltype models, BestelSorine[18] proposed an explanation of how the beating of cardiac muscle may take place, through chemical control input connected to the calcium dynamics of cells in muscles. This chemical input triggers the contractile element model.
Complete Huxley type physiological model have a greater repertoire of response, but are computationally complex. BestelSorine representation provides an efficient tradeoff between complexity and power of the model expression. It consists in an explicit computation of the first two moments reprenting stiffness and force of the contractile element, triggered by the chemical input as regard the contraction start and stop. Starting with this concept, we first adapted it to the striated muscle model to represent muscle responses under FES[19]. The multiscale muscle model was applied for voluntary contraction in a preliminary work[20]. In this study, we aim at developing EMGbased muscle model to study voluntary muscle contraction with this multiscale modeling approach compared to Hill approaches, both linear and nonlinear. Indeed, in Hill approach, the dependency on cutoff frequency choice prevents from good performances in a wide range of contraction types. The multiscale approach should improve the performances in particular regarding fast and slow contraction patterns. In addition, we consider that it is of interest to apply the same type of EMGtoforce estimation with a physiologically detailed model and not only with a phenomenological Hill model because the internal biophysical dynamics could shed new light on neuromuscular activation. Using common data sets of isometric muscle contractions, the force estimation results are compared between Hill model approaches and newly proposed approach.
Methods
Hilltype muscle model
The transformation from EMG to muscle activation is a dominant process because the estimated muscle force is assumed to be proportional to the muscle activation.
EMG processing
The EMG processing method employed in this paper is summarized. For the details, readers can refer to[8, 9].

1.
highpass filtering of raw EMG using zerolag 4th order Butterworth filter (30 Hz) to remove motion artifacts

2.
wave rectification

3.
lowpass filtering with a 2 Hz cutoff frequency using zerolag 2nd order Butterworth filter

4.
normalization with the peak of maximum voluntary contraction (MVC)
There is an alternative method instead of the lowpass filtering of the third step. It is referred to as the mean absolute value (MAV) method. A moving average filter is applied within a fixed time window, which is equivalent to lowpass filtering. The recommendations of the SENIAM project indicate a time window of 0.25 to 0.5 seconds[21]. It corresponds to the choice of cutoff frequency of the lowpass filter of 2 to 4 Hz. This critical parameter choice depends on the task dynamics. Indeed, this choice of cut off frequency is a trade off between how fast the fluctuations in amplitude can be and how reliable the estimate of the amplitude is. A low cut off frequency gives a more reliable estimate of the amplitude, but cannot capture fast changes, while higher cut off frequency results in noisier output but captures fast contractions. The cut off frequency mentioned in the SENIAM recommendations is 2 Hz for slow motions, and 6 Hz for fast motions[21]. In Hilltype model, this cut off frequency should be carefully set depending on the type of task because the envelope of the estimated force is almost determined by this lowpass filtering step. The normalized EMG is processed with the following activation dynamics, which mainly captures the delay from muscle activation to mechanical output.
Activation dynamics
The process of transforming the normalized, rectified and filtered EMG, referred as e(t) to neural activation p(t), is called the activation dynamics as shown in the Figure1 flowchart. When a muscle fiber is activated by a single action potential (AP), the muscle generates a twitch response. This delayed response can be approximated by a critically damped linear secondorder differential system[8]. Its recursive discrete form is as follows in Eq. 1 to calculate p(t_{ k }).
where d is the electromechanical delay and γ, β_{1} and β_{2} are the coefficients that define the secondorder dynamics. To obtain a positive stable solution, a set of constraints are employed, i.e. β_{1}=C_{1}+C_{2},β_{2}=C_{1}C_{2} where C_{1}<1,C_{2}<1. In addition, the gain of this filter should be maintained by ensuring γ−β_{1}−β_{2}=1 as described in[8].
Nonlinearization of muscle activation for nonlinear Hill model
Many researchers assume that the above p(t) is a reasonable approximation of muscle activation. However, a nonlinear relationship has been reported between individual muscle EMG and the joint moment for some muscles, especially at lower forces[14]. In studies on single motor units, multiple APs cause succesive twitch responses. If the time between APs decreases, twitches start to merge and the muscle force increases steadily. However, at high frequency we reach maximal contraction, where force development saturates even if the frequency increases. This induces a nonlinear relationship between the frequency and force for single motor units[8].
This modification was thus recently introduced to correspond to the mismatch of Hill model estimation at low activation levels. Indeed, in the classical linear Hill model, neural activation p(t) is directly treated as muscle activation a(t) without the nonlinear conversion. When we refer to the modified nonlinear Hill model, the following conversion from neural activation p(t) to muscle activation a(t) is performed as in Figure1. As a simple and adequate solution, Lloyd and Besier[9] proposed the following formulation:
where A is a constant parameter for the nonlinear shape factor allowed to vary between 3 and 0, with A=−3 being highly exponential and A=0 being linear.
Hilltype contraction dynamics
The muscletendon unit is modeled as a contractile element in series with an elastic tendon, as shown in Figure2(A) as same as the general Hill model[12]. The tendon element in series was modeled as having a linear stiffness k_{ t } of 180 N/mm[22] as the parameter in human triceps surae muscletendon complex. ϕ is the pennation angle between the tendon and muscle fibers. Moreover, in this paper, we assess the EMGforce relationship in isometric conditions. Then, only concentric contraction is considered where a contractile element is always shortening during contraction. The parallel elastic element is not included, since its length does not change during isometric contraction, which is assumed in this study. In addition, we would include the parallel element effect in joint level modeling rather than the muscle level when we apply the proposed approach in nonisometric condition. This common macroscopic structure as shown in Figure2(A) is used for the three types of modeling, only the contractile element dynamics change.
The Hilltype muscle model is used to estimate the force F_{ c }(t) generated by the contractile element with the general form:
where ε_{ c } is the strain of the contractile element, f_{ l }(ε_{ c }) and{f}_{v}\left({\stackrel{\u0307}{\epsilon}}_{c}\right) are the normalized forcelength and forcevelocity relationships, respectively. F^{m} is the maximum isometric muscle force.
The forcelength relationship shows a Gaussian curve around the optimal length for which the muscle generates the maximum force[12]:
where b is a constant parameter.{f}_{v}\left({\stackrel{\u0307}{\epsilon}}_{c}\right) represents the relationship between the velocity and the normalized force. The muscle contracts at its maximum velocity v_{ max } without load and slows down as the load increases. In the case of concentric contraction, this relationship can be formulated as follows:
where L_{c 0} is the length of contractile element at rest position and V_{ sh } is a constant parameter. The b and V_{ sh } parameters were obtained from[23, 24]. At each time step, the fiber velocity should be solved and the muscle fiber length is computed by forward integration using RungeKutta algorithm[8]. Since the value of ε_{ c } changed, the calculation should continue iteratively until the end of the input time series of a(t).
The muscle tendon parameters were adopted from Delp[25]. For gastrocnemius, the parameters are averages for two heads (med/lat). Only EMG of the medial head is measured, but the force of both heads is estimated using the sum of maximum force of two heads through F^{m}. F^{m} and pennation angle ϕ are also obtained from[25]. However, in the following, the final result is normalized by the maximum contraction, so the effect of these parameters is not relevant. Parameters of muscle model are detailed in Table1.
Physiological muscle model
We breifly describe the main points of the multiscale physiological muscle model. It integrates the characteristics both of macroscopic and microscopic muscle dynamics. The details of this modeling are given in[19]. This model was originally designed to represent the muscle response under FES and is applied to voluntary muscular force estimation based on EMG in this paper.
Generation of chemical input
The physiological model is composed of two inputs as shown in Figure1. One is the recruitment rate α which determines the percentage of recruited muscle fibers. The recruitment rate α has the same physiological meaning as the neural activation p(t) in Hilltype model. Then normalized p(t) without the nonlinearization was directly used for α. The second input for physiological muscle model is the chemical signal u. It represents the underlying physiological processes between the contraction and relaxation states independent from recruitment. Muscle contraction is initiated by an AP along the muscle fiber membrane, which goes deeply into the cell through Ttubules. It causes calcium release which induces the contraction process when the concentration goes above a threshold and is sustained until the concentration decreases below this threshold again[26]. Hatze[27] gives an example of calcium dynamics [ C a^{2+}] modeling. The contractionrelaxation cycle is then triggered by the [ C a^{2+}] to be defined, and associated with two phases: i) contraction and ii) relaxation. We use a delayed model to take into account the propagation time of the AP and an average delay due to the calcium dynamics. It corresponds to the electromechanical delay for e(t) into p(t) conversion in Hilltype approach. For chemical input generation, the rectified EMG was lowpass filtered with a 30 Hz cutoff frequency. Then, the chemical input u(t) was created by thresholding the extracted EMG signals as shown in Figure3. The signal measured in EMG is the summation of the APs of all different MUs. However, as we are trying to deal with one muscle model for the whole muscle force estimation, the EMG signal was considered as a representative average triggering signal. The thresholding can be assumed to reflect a muscle cell’s allornothing response corresponding to the above chemical nature of the muscle contraction process.
Thus, the chemical threshold consists basically in extracting the contraction time duration of the muscle and to indicating if the muscle is in contraction or not. The level of activation itself is managed by recruitment rate which is classically used in EMG usage for Hill model. Thus, especially for the low activation state, even if the chemical input is inaccurately determined, its effect is low since the recruitment activation p(t) is very low during this period. And for the higher activation state, the contraction period is accurately captured by this method. As long as the threshold is taken little above the baseline, the threshold is less sensitive as long as it gives the picture of contraction/relaxation state. The parameters for chemical input generation was obtained from[18, 28].
Sarcomere scale
All the sarcomeres are assumed to be identical, and the deformation of both sarcomere and muscle is proportional. If S is the sarcomere length, we can write (S−S_{0})/S_{0}=(L_{ c }−L_{c 0})/L_{c 0}=ε_{ c }.
Huxley proposed that a crossbridge between actin filaments and myosin heads could exist in two biochemical states, i.e. attached and detached states as in Figure2(B). Filaments sliding is the result of interactions between the myosin crossbridges and the thin actin filaments. The crossbridges reversibly bind to actin and produce a mechanical impulse, which in turn results in force transmission along the filaments. M and A in the figure represents the myosin head and actin binding site respectively.
One myosin head can only attach to one actin site. Then the dynamics of the fraction n(y,t) of the attached cross bridges is given by
where h is the maximum elongation of the myosin spring, x is the distance and y the normalized distance between the actin binding site and the myosin head:y=\frac{x}{h}. n(y,t) is a distribution function representing the fraction of attached cross bridges relative to the normalized relative position y. To ensure a one way displacement, this attachment is considered possible only when y is between 0 and 1.{S}_{0}{\stackrel{\u0307}{\epsilon}}_{C} represents the velocity of the actin filament relative to the myosin filament. f and g denote the rate functions of attachment and detachment, respectively.
Several f and g functions have been defined and recently a chemical input was introduced by BestelSorine[18] to modify the ability of the cross bridge to attach or not. They also proposed that these rates depend on the relative velocity between actin and myosin. Indeed, the higher the velocity is, the greater the probability to break bridges is. f and g can thus be defined by: During the contraction phase
During the relaxation phase
U_{ c } and U_{ r } are the chemical kinetics levels under contraction and relaxation phases respectively. That can be resumed as below
To complete the description at the sarcomere scale, we assume that the force generated by one attached cross bridge is modeled by a linear spring with constant stiffness. All cross bridges are parallel so the global stiffness and force generated by the whole sarcomere is proportional to the number of formed bridges. Let k_{0} (Nm ^{−1}) denote the maximum stiffness obtained when all the available bridges are attached. Let ξ(y,t) denote the elongation of a cross bridge due to the contribution of the global extension of the sarcomere, where ε_{ c }(0) is the initial value, and to the local distribution of elongations y:
The stiffness and force generated by a muscle sarcomere is obtained by computing the first and second moment of the distribution n(ξ(y,t),t)[17]. The two first moments are defined as:
From Eq. 6, it can be rewritten as follows:
Contrary to Zahalak approximation, we do not make any assumption on the distribution, rather the chosen f and g functions allows for a straightforward computation of the stiffness and the force generated by the contractile element. It provides a computational effective model that however relies on Huxley theory.
Myofiber and muscle scale
This set of differential equations can be easily extended to the whole muscle fiber considering that each fiber is composed of identical sarcormeres in series. Then{k}_{f}={k}_{s}\frac{{S}_{0}}{{L}_{c0}} and F_{ f }=F_{ s }. In addition, the maximum available cross bridges varied depending on the relative length of the contractile element. This is known as the forcelength relationship (f_{ l }(ε_{ c })). Contrary to previous studies where this effect is introduced at the macroscopic level, we take into account this relation at the microscopic scale[17]. Indeed, this relation is directly linked to the maximum available actin and myosin sites[16]. For the stiffness k_{ f } and the force F_{ f } at the fiber scale including Eq. 9, we get:
Next, we introduced the recruitment process at muscle scale. At each contraction phase, the recruitment ratio is updated but remains constant during its phase itself. Let k_{ c } and F_{ c } denote the stiffness and the force for whole contractile element, and N the number of all MUs. The recruited number is written as α N using recruitment ratio α. Finally the complete model of contractile element can be written as a set of differential equations:
where k_{ m }=S_{0}N k_{0}f_{ l }(ε_{ c })/L_{c 0},F_{ m }=N k_{0}h f_{ l }(ε_{ c })/2. k_{0} (Nm ^{−1}) is the maximum stiffness obtained when all the available bridges are attached.
Recruitment rate is thus mixed with the dynamics of the contraction relaxation cycle providing a wider range of mechanical responses than Hill approaches. Compared to Zahalak approximation, we provide a two input model with intricate dynamics due to both internal muscle dynamics and recruitment rate dynamics in a straightforward way.
Computation of the whole dynamics
For the macroscopic representation, the same configuration with Hill models as in Figure2(A) is used including the muscle tendon parameters. The contractile element is replaced with the above nonlinear differential equations.
The dynamics of the contractile element coupled with the tendon in series should also comply with the following equation:
In order to have one solution for two equations of{\stackrel{\u0307}{F}}_{c}, we must verify that: k_{ t }L_{c 0}/c o s ϕ+k_{ c }L_{c 0}−F_{ c }>0. The differential of elongation ε_{ c } is computed with the following equation:
{S}_{{\stackrel{\u0307}{\epsilon}}_{c}} is the sign of{\stackrel{\u0307}{\epsilon}}_{c}. From the condition: k_{ t }L_{c 0}/c o s ϕ+k_{ c }L_{c 0}−F_{ c }>0,{S}_{{\stackrel{\u0307}{\epsilon}}_{c}} can be obtained from the sign of these terms:{k}_{t}{L}_{0}\stackrel{\u0307}{\epsilon}/\mathrm{cos\varphi}+{F}_{c}u\alpha {F}_{m}{\Pi}_{c}\left(t\right){U}_{c}. Then we can compute{\stackrel{\u0307}{k}}_{c} and{\stackrel{\u0307}{F}}_{c} with (14). The internal state vector of this system should be set as\mathbf{x}=\left[\phantom{\rule{0.3em}{0ex}}\begin{array}{lll}{k}_{c}& {F}_{c}& {\epsilon}_{c}\end{array}\right].
Experimental measurement and data processing
Four healthy volunteers (2 males and 2 females, age =30.8±1.3 years, mass =66.5±16.0 kg, and stature =1.7±0.13 m) participated in the study after signing an informed consent form. The study was approved by the Agence française de sécurité sanitaire des produits de santé (Afssaps) committee for persons’ protection managed by CHRU Montpellier (2011A0103338), dedicated to the study of the lower limb biomechanics within ANR SoHuSim project.
The subjects were seated on a chair with their right foot attached to a Biodex dynamometer (Biodex Medical Systems Inc., New York, USA) as shown in Figure4. The setting was 90° for the ankle joint and 110° for the knee joint, while straps were used on the pelvis and shoulders to secure the subjects’ position on the chair. The axis of the dynamometer was aligned with the ankle’s rotation axis. The torque around the ankle joint was measured when it was voluntarily generated for the plantar flexion. For EMG measurements, bipolar surface Ag/AgClelectrodes were placed on the muscle belly of the medial gastrocnemius (GAS) and soleus (SOL) with an interelectrode distance of 20 m m. The reference electrode was placed on the left patella. Synchronous acquisition of the force and differential EMG signal was performed with a sampling frequency of 2048 H z by the EMG100C amplifier and Biopac MP100 system (Biopac Systems, Inc., Santa Barbara, USA).
In this preliminary test, we focused on the isometric condition and tried to compare the output force on a normalized scale against the output during MVC. The common parameters between three approaches were set at the same value as summarized in Table1. In order to maintain similar conditions, the common macroscopic model as in Figure2(A) was used for all the three approaches. Only the dynamics representation of the contractile element was replaced. A was the only parameter identified for each subject. A was obtained by minimizing the root mean square difference between the measured and the model estimation at MVC. Other muscle parameters except A are taken from literatures as explained in the model section.
The estimated muscle force is multiplied by the moment arm in order to obtain the torque. The moment arm was estimated from the Hawkins data[29] from the joint angle in the measured condition. The moment arms for MG and SOL are 0.0515 [m] and 0.0464 [m], respectively. The contribution ratio was calculated using the values reported by Delp[25]. This was obtained as the maximum force by the moment arm considering the pennation angle. The resulting ratio is MG 0.41 vs SOL 0.59. The SUM in the estimation result is plotted using the ratio as the sum of the two muscles.
The generated input command u(t) and activation rate α were given to the contractile element of the physiological model and the active stiffness k_{ c } and the muscle force F_{ c } were computed using Eq. 14. The same value for F_{ m } for the three approaches were used. k_{ m } is proportional to F_{ m } from the condition in Eq. 14. It was multipied by 20 times, as shown in Table1, based on the knowledge that the muscle active stiffness reaches the tendon stiffness at approximately maximum contraction[30]. However, this muscle stiffness parameter does not influence much the force output in isometric condition.
Results
The predicted joint torque based on EMG signals was compared with the directly measured torque of the ankle joint for the plantar flexion. The comparison was made against the case in MVC in normalized scale. Thus, we avoided the situation effected by the subjectspecific difference regarding muscle strength, which is significant for torque estimation in absolute scale.
First, the normalized torques for 30% and 70% of MVC estimated by the nonlinear Hill model and multiscale model were computed as in Figure5. The normalized data of measured torque (red), normalized estimated torque of SOL (green), normalized estimated torque of GAS (magenta) and the sum of two muscles (blue) are plotted. The obtained RMS errors for the first two subjects are given in Table2 for both contraction levels and by both estimation approaches. The estimation results of physiology based model show the smooth transition compared to the estimation by Hill model keeping the realistic torque development as in the measured torque. In addition, especially for low activation levels, RMS errors in the physiology based model showed better accuracy. Note that the modified Hill model has already the nonlinearization conversion to correspond to a low activation level. The proposed physiological model doesn’t include such empirical modifications.
Next, in order to confirm the estimation ability for both fastshort and slowlong contraction, random contractions including these two types of contraction were requested to all subjects. The generated chemical inputs are shown in Figure3. Normalized torques estimated by linear Hill model (A), nonlinear Hill model (B) and multiscale physiology based model (C.1) are depicted in Figure6. In the physiological model, active stiffness and strain of the contractile element are explicitly defined in differential equations, then it can be estimated along with the muscle force. The active muscle stiffness k_{ c } is shown at Figure6(C.2). The results of linear Hill model is computed to show the contribution of the nonlinear conversion. As we observe the difference between them in Figure6(A,B), nonlinear Hill model estimation is improved especially for low level contractions. However, the errors for fastshort contraction remained high. In contrast, the estimation result as in Figure6(C.1) with the proposed method demonstrates the nice estimation for all conditions of contractions.
The obtained RMS errors for the four subjects are given in Table3 for both contraction types and by three estimation approaches. In addition, mean error of the peak torque for fastshort contraction are given. Average results is also reported as m e a n±s t a n d a r d deviation over all the subjects. Average error was 16.9% with the linear Hill model, 9.3% with the nonlinear Hill model. In contrast, the error in the multiscale model was 6.1% while maintaining the uniform estimation performance in both fast and slow contraction speeds. Mean peak error for fastshort contraction are 21.9% with the linear Hill model, 17.6% with the nonlinear Hill model. In contrast, the error in the multiscale model was 5.8%. In general, results estimated by physiology based model show the best performance both in fastshort and slowlong contractions in these random tests, as reported in Table3.
Discussion
Hilltype muscle is a phenomenological model based on experimental facts with no link to the microscopic physiology. First, even if the estimation performances are equivalent, it is meaningful to understand and capture the muscle dynamics with a more detailed representation. Indeed, the same kind of force estimation could be obtained with the newly proposed physiological model using voluntary EMG signals.
Second, previous work by Perreault[13] reported that Hill modeling errors were large for different firing frequencies and greatest at low MU firing rates. The classical linear Hill model was used. They suggested that more physiological coupling between activation and forcevelocity properties may help to deal with these issues[13]. Thus, the aim of introducing new approach based on a physiologically detailed model[19] was to develop EMGforce estimation with natural velocity dependency incorporating microscopic crossbridge dynamics on top of recruitment dynamics.
In general, the estimation accuracy of the proposed method is better than that of the linear and nonlinear Hilltype approach as in the obtained result regarding RMS errors. In addition, uniform accuracy was obtained for different types of contractions in the proposed method. In particular, the proposed model is able to properly deal with both slow and fast contraction speeds. When considering the Hill model results, Figure6(A,B) shows largest errors for the estimation of fastshort contractions. The reason of this error could be explained as follows: the signal measured in EMG is the summation of the APs of all different MUs. Even in the fastshort contraction, the amplitude of EMG itself does not reflect differences except duration with the same level as for slowlong contraction. However, the resulting fastshort contraction force is actually much less than that of slowlong contraction. This means that there is sometimes hysteresis regarding the neural command. In Hill approach, muscle activation a(t) is proportional to the resulting force. It therefore does not include the effect of the time hysteresis in contraction[13], thus Hill model can not accurately represent both short and longterm contractions simultaneously with only one choice of cutoff frequency in EMG lowpass filtering. In the proposed approach, the derivative of the contraction force is directly given by the neural command and it brings time hysteresis in force generation. It is interesting to note that this kind of uniform accuracy for different types of contractions was appeared along with the introduction of crossbridge dynamics.
In addition, it is known that there is a nonlinear relationship between frequency during contraction and force for single MUs[14]. In modified Hill model, this frequency dependency is offset only by the nonlinear conversion. This nonlinearization was recently proposed to modify the classical Hill model. This process was originally not introduced in the socalled Hill model. In fact, this modification brings much better estimation, especially at lower forces. However, even with the modification, it is not a time function so it still can not correspond to the different muscle contraction speeds. In Hilltype approach, neural activation p(t) is dominantly imposed by lowpass filtering with a 2–10 Hz cutoff frequency. The choice of cutoff frequency is very sensitive to the dynamics of p(t). p(t) is more or less proportional to the resulting force. Therefore, it is advised that the choice of cut off frequency should be carefully selected depending on the type of task in Hilltype modeling. In this study, a 2 Hz cutoff frequency was used according to the SENIAM recommendation. Sometimes by changing the cutoff frequency for different types of contraction, better results can be obtained with Hilltype model. However, in any case, the expected motion by subject is unknown in advance for practical use of EMGforce estimation, so it is difficult to have an appropriate cutoff frequency in advance. Moreover it assumes that the subject stays in the same kind of motion dynamics, that is not always the case. In the proposed model, the nonlinear activation property coming from forcevelocity dependency is internally integrated in the model dynamics.
Furthermore, the physiological model satisfies the well established properties observed for the muscle’s behavior: i) the forcelength relation is included in the definition of Eq. 13, ii) the forcevelocity relationship can be expressed in isotonic and tetanic conditions. Before isotonic contraction could occur, the muscle contracts in isometric conditions until the force generated by the muscle balances the imposed one. Isotonic contraction then becomes possible. Assuming{\stackrel{\u0307}{F}}_{c}=0 in Eq. 14, the following equation may be formulated:{F}_{c}=(\alpha {F}_{m}{U}_{c}+{k}_{c}{L}_{c0}{\stackrel{\u0307}{\epsilon}}_{c})/({U}_{c}+\left{\stackrel{\u0307}{\epsilon}}_{c}\right).
At the beginning of this phase (t = 0), k_{ c }(0)≠0 but it depends on the initial isometric contraction phase. If we defineA\left(0\right)=\frac{{L}_{c0}{k}_{c}\left(0\right)}{{U}_{c}\alpha {F}_{m}},B=\frac{1}{{U}_{c}}, in concentric contraction where{\stackrel{\u0307}{\epsilon}}_{c}<0, we get a Hill type forcevelocity relation:
This correspondence can be found in Eq. 5. It can be verified that our physiological muscle model integrates the forcevelocity relation naturally when considering the actinmyosin cross bridge. Finally, this improvement in estimation accuracy is in line with the suggestion put forward in a previous work[13] on introducing crossbridge model.
Future study will be focused on simulations with multiple sets of the proposed model corresponding to multiple MUs in one muscle considering size principal, fast/slow fibers. However, for the application of voluntary EMGforce estimation, we should keep the simplicity as Hilltype model. Future work will also focus on increasing the number of experimental tests and we should develop an algorithm to statistically determine if a muscle is turned on/off regarding the chemical thresholding. In addition, the further interpretation of neuromuscular system both in voluntary and FES activations should be pursued. We have already performed the study on the identification and validation of the physiological model in invivo rabbit experiment[31] and in paraplegic subjects including nonisometric situation under FES[32]. In order to be applied to a broader range of clinical situations, further investigations in isokinetic and isotonic cases should be carried out.
Conclusion
In this paper we presented a method that allows estimation of muscle force from EMG signals with a multiscale physiology based model with a link to underlying microscopic filament dynamics. The experimental results highlight the feasibility of the torque estimation and its comparison with Hilltype models using the same EMG signals. This is the first report on EMGforce estimation based on a multiscale physiology model integrating Hilltype and microscopic cross bridge representations. The proposed method features:

a novel physiologically detailed model for EMGforce estimation in the place of a phenomenological Hilltype muscle model,

the estimation improvement especially for different types of contraction incorporating dual dynamics of recruitment and microscopic crossbridge formation toward coupled activationvelocity relationship.
References
Delp SL, Anderson FC, al E: OpenSim: Opensource software to create and analyze dynamic simulations of movement. IEEE Trans Biomed Eng 2007, 54(11):1940–1950.
Flash T, Hogan N: The coordination of arm movements: an experimentally confirmed mathematical model. J Neurosci 1985, 5(7):1688–1703.
Kawato M, Maeda Y, Uno Y, Suzuki R: Trajectory formation of arm movement by cascade neural network model based on minimum torquechange criterion. Biol Cybern 1990, 62(4):275–288. 10.1007/BF00201442
Smaga S: Tremor. Am Fam Physician 2003, 68: 1545–1552.
Kobetic R, Triolo RJ, Uhlir JP, Bieri C, Wibowo M, Polando G, Marsolais EB, Davis JA: Implanted functional electrical stimulation system for mobility in paraplegia:a followup case report. IEEE Trans Rehabil Eng 1999, 7(4):390–398. 10.1109/86.808942
Guiraud D, Stieglitz T, Koch KP, Divoux JL, Rabischong P: An implantable neuroprosthesis for standing and walking in paraplegia: 5year patient followup. J Neural Eng 2006, 3: 268–275. 10.1088/17412560/3/4/003
Stegeman DF, Blok JH, Hermens HJ, Roeleveld K: Surface EMG models: properties and applications. J Electromyogr Kinesiol 2000, 10: 313–326. 10.1016/S10506411(00)000237
Buchanan TS, Lloyd DG, Manal K, Bessier TF: Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. J Appl Biomech 2004, 20: 367–395.
Lloyd DG, Bessier TF: An emgdriven musculoskeletal model to estimate muscle forces and knee joint moment in vivo. J Biomech 2003, 36: 765–776. 10.1016/S00219290(03)000101
Koo TKK, Mak AFT: Feasibility of using emg driven neuromusculoskeletal model for prediction of dynamic movement of the elbow. J Electromyogr Kinesiol 2005, 15: 12–26. 10.1016/j.jelekin.2004.06.007
Hill AV: The heat of shortening and the dynamic constants in muscle. Proc R S Ser B; London 1938, 126: 136–195. 10.1098/rspb.1938.0050
Zajac FE: Muscle and tendon: properties, models, scaling and application to biomechanics and motor control. CRC Critic Rev Biomed Eng 1989, 17: 359–411.
Perreault E, Heckman C, Sandercock T: Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J Biomech 2003, 36: 211–218. 10.1016/S00219290(02)003329
Manal K, Buchanan TS: A oneparameter neural activation to muscle activation model. J Biomech 2003, 36: 1197–1202. 10.1016/S00219290(03)001520
Huxley AF: Muscle structure and theories of contraction. Prog Biophys Biophys Chem 1957, 7: 255–318.
Leeuwen JL: Optimum power output and structural design of sarcomeres. J Theor Biol 1991, 19: 229–256.
Zahalak GI: A distributionmoment approximation for kinetic theories of muscular contraction. Math Biosci 1981, 55: 89–114. 10.1016/00255564(81)900146
Bestel J, Clement F, Sorine M: A biomechanical model of muscle contraction. Proc Int Conf Med Image Comput ComputAssist Intervent MICCAI 2001, Lect Notes Comput Sci 2001, 2208: 1159–1161.
Makssoud HE, Guiraud D, Poignet P, Hayashibe M, Wieber PB, Yoshida K, Azevedo C: Multiscale modeling of skeletal muscle propertiesexperimental validations in isometric conditions. Biol Cybern 2011, 105(2):121–138. 10.1007/s0042201104457
Hayashibe M, Guiraud D, Poignet P: EMGtoforce estimation with fullscale physiology based muscle model. Proc IEEE/Int Conf Intell Robots Syst 2009, 1621–1626.
Hermens HJ, Freriks B, Merletti R, Stegeman D, Blok J, Rau G, DisselhorstKlug C, Hagg G: European recommendations for surface electromyography, results of the seniam project. Enschede, The Netherlands: Roessingh Research and Development; 1999.
Lichtwark GA, Wilson AM: Optimal muscle fascicle length and tendon stiffness for maximising gastrocnemius efficiency during human walking and running. J Theor Biol 2008, 252(4):662–673. http://www.ncbi.nlm.nih.gov/pubmed/18374362 10.1016/j.jtbi.2008.01.018
Stroeve S: Learning combined feedback and feedforward control of a musculoskeletal system. Biol Cybern 1996, 75: 73–83. 10.1007/BF00238741
Stroeve S: Impedance characteristics of a neuromusculoskeletal model of the human arm I: posture control. Biol Cybern 1999, 81: 475–494. 10.1007/s004220050577
Delp SL: Surgery simulation: a computer graphics system to analyze and design musculoskeletal reconstructions of the lower limb. PhD thesis, Dissertation, Stanford University, CA, USA; 1990
Keener J, Sneyd J: Mathematical Physiology, II: Systems Physiology. Series 8, Vol. 2. New York: Springer; 2009.
Hatze H: A general myocybernetic control model of skeletal muscle. Biol Cybern 1978, 28: 143–157. 10.1007/BF00337136
Bestel J: Modele Differentiel de la Contraction Musculaire controlee application au system cardiovasculaire. PhD thesis, Dissertation, Universite de Paris IX, France; 2000
Hawkins D, Hull M: A method for determining lower extremity muscletendon lengths during flexion/extension movements. J Biomech 1990, 23: 487–494. 10.1016/00219290(90)90304L
Cook C, McDonagh M: Measurement of muscle and tendon stiffness in man. Eur J Appl Physiol Occup Physiol 1996, 72(4):380–382. 10.1007/BF00599700
Hayashibe M, Guiraud D, Poignet P: In vivo identification of skeletal muscle dynamics with nonlinear Kalman filter: comparison between EKF and SPKF. ISRN Rehabil 2013. 2013: Article ID 610709
Benoussaad M, Poignet P, Hayashibe M, Azevedo C, Fattal C, Guiraud D: Experimental parameter identification of a multiscale musculoskeletal model controlled by electrical stimulation: application to patients with spinal cord injury. Med Biol Eng Comput 2013, 51(6):617–631. 10.1007/s115170131032y
Acknowledgements
This work was supported by French National ANR SoHuSim project (http://sohusim.gforge.inria.fr/).
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
MH proposed and built the structure of EMGforce estimation with a multiscale model and wrote the first draft. DG was in charge of the physiological muscle modeling and proof reading. MH performed the scientific computation and experiments. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Hayashibe, M., Guiraud, D. Voluntary EMGtoforce estimation with a multiscale physiological muscle model. BioMed Eng OnLine 12, 86 (2013). https://doi.org/10.1186/1475925X1286
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1475925X1286
Keywords
 Muscle model
 EMG
 Muscle force estimation
 Hill model
 Multiscale physiology