A comparison among different Hill-type contraction dynamics formulations for muscle force estimation

Muscle is a type of tissue able to contract and, thus, shorten, producing a pulling force able to generate movement. The analysis of its activity is essential to understand how the force is generated to perform a movement and how that force can be estimated from direct or indirect measurements. Hill-type muscle model is one of the most used models to describe the mechanism of force production. It is composed by different elements that describe the behaviour of the muscle (contractile, series elastic and parallel elastic element) and tendon. In this work we analyze the differences between different formulations found in the literature for these elements. To evaluate the differences, a flexo-extension movement of the arm was performed, using as input to the different models the surface electromyography signal recorded and the muscle-tendon lengths and contraction velocities obtained by means of inverse dynamic analysis. The results show that the force predicted by the different models is similar and the main differences in muscle force prediction were observed at full-flexion. The results are expected to contribute in the selection of the different formulations of Hill-type muscle model to solve a specific problem.


Introduction
The study of mechanical muscle models is one of the major topics in Biomechanics of human movement.The analysis of the forces that produce a given movement (inverse dynamics, ID) or the movement induced by a set of muscle forces or activations (forward dynamics, FD) are typical problems that need the description of muscle mechanical properties.Since the classic work of Hill (1938), the number of models proposed has been grown up, ranging from simple models to the most complex (as e.g. the model developed by Hatze (1977) with up to 50 parameters required to describe the motion of a simple joint).Muscle models are commonly categorized into three groups according to Winters and Stark (1987): (i) simple second-order models: it is a "black box" approximation where the inputs are either the neural signal or the external load and the output corresponds to either the joint position or torque; (ii) Hill-based lumped-parameter model, which is the most widely used and will be described in Sect.2.1; and (iii) Huxley-based distributed-parameter models that attempts to explain correctly the mechanism of contraction with great accuracy but at a high computational effort.Some recent works, not described nor discussed here, deal with fractional order models of different kinds of muscles, such as the gastrocnemius muscle in Sommacal et al. (2007Sommacal et al. ( , 2008) ) or the hamstring muscle group in Grahovac and Žigić (2010).This type of muscle attempts to describe the viscoelastic properties of muscle tissue as a whole (HosseinNia et al., 2012).
The mechanical behaviour of muscle tissue can be described by means of passive elements such as springs and damping elements (SE and DE, respectively).These elements, combined properly, allow to understand the response of muscle tissue under compressive and tensile loads.In the literature (see, e.g., Yamaguchi, 2001) it is possible to find different models combining the properties of those mechanical components: the Maxwell model (Fig. 1a) uses both elements attached in series.Contrariwise, in the Voight model (Fig. 1b), those elements are used in parallel.Lastly, the Kelvin model (Fig. 1c), modifies the Voight model to include an additional spring in series with the DE.The dif- ferent combinations of springs and dashpots are intended to improve the physiological response, however, as these models are composed by passive elements, they are not able to reproduce properly the active muscle contraction.This issue was solved by introducing the contractile element (Hill, 1938).This model (Fig. 1d) combines the passive properties of the Kelvin model with the active properties given by the contractile element (CE).Nowadays, the Hill-type muscle model is the most used in biomechanical studies involving muscular coordination (Zajac, 1989;Van Soest and Bobbert, 1993;Van Den Bogert et al., 1998;Thelen, 2003;Silva, 2003;Ackermann and Schiehlen, 2006;Ackermann, 2007;García-Vallejo, 2010;Alonso et al., 2012).
Figure 1e shows an schematic representation of the muscle-tendon unit using the Hill-type muscle model.In this sense, the elastic properties of the tendon are represented by a spring attached in series with the Hill-type muscle model.Regarding the muscle, the CE is responsible of the active force generated in the muscle, and two non-linear passive springs describe the properties of muscle tissue: on the one hand, the series elastic element (SE) represents the elasticity of the actin-miosyn crossbridges (Yamaguchi, 2001) and, on the other hand, the parallel elastic element (PE) describes the passive elastic properties of the muscle fibers.The SE element can be neglected with little inaccuracy if the study does not involve short-tendon actuators (Silva, 2003;Zajac, 1989).
The contribution of tendon is important for physiological purposes, but is often neglected for simplicity, as the process to distinguish between tendon and muscle length increases the difficulty of the problem.However, as pointed out by Yamaguchi (2001), tendon elasticity is important if the tendon stretches an amount approaching the fiber length of a particular muscle.This fact is relevant in some muscles such as the soleus or the gastrocnemius at the ankle joint or the rectus femoris at the knee.In the other cases, or as a first approximation, it is possible to consider the tendon as a rigid element, and include as the tendon length its slack length (l T slack ), that is, the length for which the tendon just begin to resist lengthening.
The models available in the literature, based on the Hilltype description, differ basically in the amount of parameters to define the muscle mechanical behaviour.An extensive review of Hill-based muscle models was performed by Winters (1990b).In the present work, we analyse two parameterbased models derived from the Hill-type description of muscle tissue.Specifically, the model proposed by Van Soest and Bobbert (1993), the one by Thelen (2003), and a mathematical adjustment proposed by Kaplan (2000) and extended by Silva (2003) are presented and discussed here.The objective is to have a reference collecting the most widely used muscle models in biomechanics, highlighting their characteristics and discussing their use for specific problems.To do so, we analyse the differences observed in muscle force production by using the proposed models.Moreover, the variability of each muscle element between models is also studied.The order in which the different elements will be addressed is summarized in Table 1.

Mathematical description of muscle dynamics
The categories described previously are used depending on the type of problem to solve.From the mechanical point of view, although Huxley's muscle model describes precisely the chemical and mechanical process that take place in the muscle contraction, its use is not recommended in coordination studies that involve several muscle actuators as the complexity of the problem increases considerably.In this type of studies, the mechanical behaviour of muscle tissue can be modelled by means of passive elements such as springs and damping elements (SE and DE, respectively), and therefore, by using Hill-type models.The concepts introduced in this section are deeply described in Winters (1990a) and Zajac (1989).In a Hill-type muscle model, the dynamic behaviour of the active element, i.e., the contractile element can be expressed in terms of two cascaded differential equations (2a), the excitation-to-activation dynamics (Eq. 1, Fig. 2b) and the activation-to-force (contraction) dynamics (Eq.2, Fig. 2c):

Formulation of muscle's dynamic equations
(1) The expression in Eq. (1) transforms the muscular excitation signal, u, from the central nervous system into an activation signal, a, representing muscle recruitment level.This excitation-to-activation dynamics corresponds to the time lag due to the electrochemical process produced by the action potential that leads to muscle contraction.Equation (2) defines the force exerted by the muscle as a function of the physiological state, that is, dependent of the activation level (a), the muscle length (l MT ), the contraction velocity (v MT = − lMT ) and the actual force (f MT ).The following sections will describe the different dynamics that govern muscle activity and the contribution of the different tissues to the force development.

Excitation-to-activation dynamics
Excitation-to-activation dynamics or, simply, activation dynamics (Eq.1), represents the muscle fibers recruitment state.The most widely used expression to obtain muscle activation from a set of neural excitation is (Nagano and Gerritsen, 2001): where t 2 = 1/t d and t 1 = 1/(t a −t 2 ) are time constants, u(t) is the neural excitation 0 ≤ u ≤ 1, a(t) is the muscle activation, where 0 ≤ u ≤ 1, t a is a time activation constant and t d is the deactivation constant.The values for those constants are taken from Umberger et al. (2003).
F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation As depicted on Fig. 2b, activation dynamics represents a delay between input and output, that is, the delay between the reception of the signal (neural excitation) and the transformation of those signals into action potentials to activate the muscle fibers due to electrochemical process that take place in muscle tissue.
The direct acquisition of neural excitations represents a major drawback.Instead, a set of activations can be obtained from indirect measurements such as electromyography (EMG) signals.Buchanan et al. (2004) proposed an EMG-to-activation relationship, as EMG signals can be measured easily.This relationship, namely the A-model, can be written as: where a j (t) is the activation of the j th muscle, u j (t) the processed sEMG signal of the j th muscle, and A a nonlinear shape factor constrained to −3 < A < 0 according to Buchanan et al. (2004) or −5 < A < 0 as reported by Sartori et al. (2009) (in this work A = −4).

Activation-to-force dynamics
The force exerted by a muscle depends on the number of muscle fibers recruited to perform a given movement (activation level), as well as on the actual fiber lengths and the contraction velocity, as expressed in Eq. ( 2).This mathematical description can be derived from the muscle-tendon components.
According to the distribution of the muscle elements shown in Fig. 1e, and considering tendon and muscle as springs attached in series, it is possible to write: where superscript MT is referred to muscle-tendon unit, T is related to tendon and M corresponds to muscle.If one expand those equations for tendon and muscle (CE) force, the following expression is obtained: where F M 0 is the maximum isometric force, f T (l T ) is the normalized force exerted by the tendon that depends on its length (l T ).The terms F PE and F M are, respectively, the forces exerted by the parallel element and the contractile element.The force exerted by this last component can be expressed as (Buchanan et al., 2005): where a is the muscle activation and f L ( l M ) and f V ( v M ) the force-length and force-velocity relationships, respectively.These relationships depend on the normalized muscle length l M = l M /l M 0 and normalized contraction velocity, where v max = l M 0 /τ c , being τ c a time constant (set to 0.1 s in this work).
In the next sections, the mathematical models for each element of the Hill-type muscle model, i.e., tendon, SE, PE and CE will be described, specifically the models proposed by Van Soest and Bobbert (1993), Thelen (2003) and Kaplan (2000)-Silva (2003).

Tendon
The force exerted by the tendon element is expressed in terms of the tendon strain, ε T , however, other authors use parabolic relationships that depends directly on the tendon length.
The relationship proposed by Thelen (2003) to obtain the force exerted by the tendon as a function of the stiffness is: where ε T = l T −l slack l slack .Van Soest and Bobbert (1993) use a simple relationship to obtain the force developed in the tendon.In its work the tendon is presented as a quadratic spring and the expression to obtain the force is defined by: where Kaplan-Silva's description of tendon is not included in this work as they considered in their works that this element could be neglected in slow movements to reduce the dimensionality of the problem.The comparison of the models presented here is shown in Fig. 3b.The relationship given by Thelen is composed by an exponential function followed by a linear expression.On the contrary, the relationship defined by Van Soest and Bobbert is a quadratic curve.Both expressions show slight differences in the region where the muscle works.However, the curves diverge for greater values.

Muscle
Muscle tissue transforms the muscle fiber recruitment level into muscle contraction.Muscle tissue is defined by the contractile element (CE), that describes the contraction process, and the parallel element (PE), which defines the passive force exerted by the muscle when lengthen over the optimal length.Both elements are described right after.

Parallel element
The parallel element represents the elasticity of the tissue attached in series to the muscle.According to Thelen (2003), the force-length relationship for the parallel element is defined by: where k PE is a shape factor and ε M 0 is the parallel element strain.In this work, k PE = 3.0 and ε M 0 = 0.5.Silva (2003), based on the work of Kaplan (2000), proposed an analytical expression for the PE: Van Soest and Bobbert (1993) does not include an expression for this element, however, a quadratic curve is described in their work.Martins et al. (1998) present a relationship of this type for the PE: The comparison between different formulations is shown on Fig. 3a.The results show a qualitative similarity between curves, specially in the normalized muscle length interval l M ∈ [0, 1.6].Values above l M = 1.6 present significant divergences, however values over this limit are rarely reached in normal movements.

Contractile element
The contractile element is responsible of the active force generation in muscle tissue.As expressed by Eq. ( 7), written below for a better understanding, the force developed by the CE depends on the force-length and force-velocity relationships, that is: The mathematical expressions for both relationships are detailed below.

Force-length relationship
According to Thelen (2003), the force-length relationship is described by: where γ represents the half-width of the curve for f L = 1/e.The value of γ is set to 0.45 in this work.Silva (2003), again, proposes an analytical expression for the force-length relationship f L : Contrariwise, Van Soest and Bobbert (1993) use directly a force-length-velocity curve.The expressions for the forcelength and force-velocity relationships can be derived from it by setting v M = 0 and l M = 1 respectively (for a detailed explanation see Ackermann, 2007).Both curves are depicted in Fig. 4a and b.The proposed relationship is: where a is the muscle activation, f ac = min (1, 3.33 • a) and, A r and B r are Hill's constants corrections according to Van Soest and Bobbert (1993).The values for those constants are set to A r = 0.41 and B r = 5.2.Finally, the isometric force relative to the maximum force, f iso can be obtained as: were width is the maximum range of force production relative to l M 0 .
-Eccentric contraction v M > 0 : where shape factors b 1 , b 2 and b 3 can be obtained as: where f max v is the maximum normalized achievable muscle force when the muscle is lengthening (in this work this value is set to 1.6), f ac and S f are shape factors usually set to 1 to reduce computational effort.
The comparison of the different force-length relationship formulations are represented in Fig. 4b.

Force-velocity relationship
According to Thelen (2003), the force-velocity relationship can be written as: where k CE1 and k CE2 are force-velocity shape factors (0.25 and 0.06 respectively in this work)and f max v is the maximum normalized achievable muscle force when the muscle is lengthening (in this work this value is set to 1.6).
The analytical expression given by Silva (2003) is defined as: The expression given by Van Soest and Bobbert (1993) was already addressed on the previous section.The comparison of the different curves for the force-velocity relationship is represented in Fig. 4b.Besides, it is possible to represent in a surface plot the force-length-velocity relationship (Fig. 4c).Results are similar in qualitative terms, and therefore, any of the proposed formulations can be used to obtain muscle forces.However, there are some slight differences that must be mentioned.The force-length relationship described by Silva (2003) is not centred on l M = 1 and the bell shape is slenderer than the other approaches.However, the force-velocity relationship is quite similar to the model proposed by Thelen (2003).In this case, the model proposed by Van Soest and Bobbert (1993) presents significant differences in the interval v M ∈ [−0.5, 0.5], that may lead to variations of 20 % in F CE compared to the other ones in the same region.
Lastly, it is possible to represent the combined actuation of the PE and CE.As it can be shown in Fig. 2c, if the muscle tissue is stretched beyond its optimal length (l M 0 ) the passive force generated by the PE becomes significant.

Experimental setup
In order to test the different muscle models analysed in this work, the following experiment was carried out.The idea is to obtain a set of muscle forces to study the variability of the results by using different muscle models in which all the inputs are known.To do so, a simple movement, an arm flexoextension movement under load was performed.The recordings consist of the acquisition of the kinematics of reflective markers attached to anatomical landmarks and muscle activity (EMG signal) of the long head of the biceps brachii.On the one hand, the EMG signal is used as input in the A-model to obtain a set of activations.On the other hand, the kinematics of the flexo-extension movement under load was used to obtain the muscle lengths and contraction velocities to be used with the activations to obtain the muscle forces.This process was performed in OpenSim (Delp et al., 2007) using the upper extremity model from Holzbaur et al. (2005).The model was scaled to subject's dimensions by using the recorded positions of the reflective markers.Parameters such as l M 0 or l slack were obtained after the scale.Others, as F M 0 , were taken directly from literature (Holzbaur et al., 2005).
The experimental test was carried out by a voluntary subject (26 years, healthy with no muscular nor neurological disorders).The test was performed under subject's consent and approved by the local Ethics Committee at University of Extremadura.The acquired movement consists of a dumbbell weightlifting with the right arm, from full extension at anatomical position to full flexion and return to the initial position.A load of 2.5 kg was used.A total of seven recordings were carried out for the experiment.The voluntary was trained previously to maintain constant the velocity of the cycle in order to obtain an adequate repeatability.Seven reflective markers were placed according to Nigg and Herzog protocol (Nigg and Herzog, 1999) in the right arm.The motion was recorded with 12 infra-red light cameras Opti-Track V100:R2 at 100 Hz.The motion capture system was fully synchronized with the Trigno ™ Wireless System from Delsys ® .The muscle activity was recorded on the superficial long head of the biceps brachii.The sEMG electrode was placed following the recommended standard of the Surface EMG for a Non-Invasive Assessment of Muscles (SENIAM) project (Hermens et al., 2000).The skin was abraded and cleaned with isopropyl alcohol.Then, a thin layer of conductive gel was extended along the point of application of the electrode, that is, in the middle zone of the muscle, far from innervated an tendinosus zones (Criswell, 2010).The sEMG signals were filtered by means of singular spectrum analysis (see Romero et al., 2015 for further details) and the first component was retrieved as the trend of the signal, i.e., the input to the A-model.The values were normalized to the maximum voluntary contraction (MVC) following the recommendation given in Konrad (2006).The processing of kinematic and sEMG signals was performed in MATLAB ® , running on an Intel ® Core ™ i5 CPU at 3.20 GHz.
In order to assess the differences in muscle force production between models, a validation metric for curve comparison proposed by Geers (Geers, 1984;Lund et al., 2011) is used.This metric allows to quantify independently differences in magnitude (M) and phase (P ).A combined error (C) based on the previous differences is also used.The expressions can be written as: where the values v mm , v cc and v mc can be defined as: www.mech-sci.net/7/19/2016/Mech.Sci., 7, 19-29, 2016 where m(t) and c(t) are the measured and the computed signals and [t 1 , t 2 ] is the time interval in which the analysis is performed.In this work we perform a pairwise comparison (see Table 2), therefore m(t) can be defined as the interest signal, which is compared with the signal c(t).Major similarities between signals are related with indices near to zero.

Results and discussion
The differences between models have already been highlighted along the methods section.The variability of the acquired signals during the different contractions is represented in Fig. 5.As shown in the figure, the standard deviation is qualitatively low for the reconstructed normalized muscle length (represented by the limits of the shaded region).In the same sense, the measurements of sEMG signals present certain variability, mainly due to differences in muscle fibers recruitment by the central nervous system to prevent fatigue or damage.
The results related to the different elements of the muscle will be discussed in the order they were described in Sect. 2. First, regarding to tendon (see Fig. 3a), the presented models show similarities for typical values of muscle lengths in normal activities, however for greater values the differences are important.However, due to the high slope of the tendon force curve it is usual to consider this component as a stiff element.This fact reduces the dimensionality of the problem, however, as commented before, tendon elasticity is important if the tendon stretches an amount approaching the fiber length of a particular muscle, and therefore it must be considered in these cases (Zatsiorsky and Prilutsky, 2012).
Regarding to muscle elements (see Figs. 3b, 4), the three proposed models show qualitative similarities.For the parallel element, a notable difference between Silva-Kaplan's model and the rest for values below l M = 1.5 (see Fig. 3b).As the range of normalized muscle length of the recorded movement corresponds to l M ∈ [0.5, 1.2] (Fig. 5a) significant differences will appear in the normalized PE force, as shown in Fig. 6a, and therefore in the total normalized muscle force (Fig. 6c).In the case of the contractile element, it was pointed out previously that the main variability was found on the force-length relationship, and therefore, differences in muscle force production will be extended to the total muscle force.These differences in the normalized CE force are depicted in Fig. 6b.Considering the range of normalized muscle force and the shape of the force-length relationship (Fig. 4a), the differences in the range 35-60 % are the ones expected.In fact, similar values are observed in the neighbouring of the optimal length, where all the force-length curves present similar values.The major differences are observed at 50 % of the cycle, corresponding to full flexion.These variations correspond to the differences in the bell-shaped curve of the force-length relationship, as for full-flexion values ( l M ≈ 0.5) there are considerable differences in the force-length relationship.Silva-Kaplan's model presents major variation in this point.Moreover, when computing the total force (PE + CE) the differences are increased in the limits of the cycle as a consequence of the contribution of the PE.
Regarding to the results shown in Table 2 the main differences between models can be observed in magnitude rather than phase.According to the results, Thelen and Silva-Kaplan's models present major similarities between them compared to Thelen and Soest and Bobbert's.These differences are more evident during the extension phase with deviations of 40.68 and 28.38 % between Soest and Bobbert's models and those of Thelen or Silva-Kaplan, respectively.Differences between curves are lower during flexion in both

Table 2.
Comparison of different formulations by using Geers index.F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation magnitude and phase as expressed by Geer's indices.On the contrary, the major differences observed in the models are related to the full-flexion state, where the models show certain divergences.At full-flexion, Van Soest and Bobebrt's model overestimates the force whereas Kaplan-Silva does the opposite (see Fig. 6c and Table 2).Moreover this last model underestimates also the force at full-extension compared to the other models.The reason of these deviations arises from the muscle description itself.Kaplan-Silva's model is based in the adjustment of Hill's model for multibody formulation purposes and the mathematical description depends only on the normalized muscle length and contraction velocity.On the contrary, Thelen and Van Soest and Bobbert's models contains enough muscle parameters to adjust the mechanical properties of muscle tissue adequately.
On the other hand, if one attends to the computational effort, Kaplan-Silva (0.021811 s) outperforms Thelen (0.023208 s) and Van Soest and Bobbert's (0.026339 s) muscle descriptions.The differences may not be significant in this simple off-line experiment but can make the difference in real time applications.Only the comparison of these results with in vivo measurements will provide the better solution for each case, as a trade-off between accuracy and execution time is required.Finally, this study presents some limitations that must be mentioned.The muscle parameters considered in this study have been taken from literature.Moreover, the methodology proposed to analyze the differences between the presented models is based in a simple exercise for a single muscle.Complex exercises involving several muscles may report different results as the EMG signal used as input to the A-model may contain noise or crosstalk if not filtered properly.Nevertheless, the results presented here can be used as a reference to muscle force production during arm flexo-extension movement and to evaluate differences between models prior their selection for biomechanical studies.

Conclusions
In this paper, the mechanical expressions to represent muscle tissue have been presented.To compute muscle efforts, the most significant models available in the literature have been presented.As it has been shown, there are slight differences between the obtained muscle efforts.However, there are some relevant aspects that should be mentioned.The model proposed by Soest and Bobbert gives a detailed description of muscle behaviour, however, this model involves the measurement of too many parameters and most of them are not available in practice.By contrast, Silva-Kaplan's model is based in analytical expressions, and from the computational point of view, their use is more efficient as the expressions are simpler than the other cases, however, for some specific applications it may be interesting to include physiological information which is not provided with this formula-tion.Thelen's model contains enough parameters to describe muscle force production and it has been implemented in a widely used and validated open source software: OpenSim.The use of each model depends mainly on the computational effort and the control of the parameters involved in the experiments.In this way, Silva-Kaplan's model can be used in experiments in which computing time and computational effort are critical but not the use of physiological parameters whereas Thelen's model provides enough control of physiological parameters to perform simulations in a reasonable time.Therefore, the selection of the most appropriate model depends on the specific conditions of the problem.

Figure 1 .
Figure 1.Muscle models (tendon not included in a-d models).(a) Maxwell model.(b) Voigth model.(c) Kelvin model.(d) Hill model.(e) Nomenclature used in the Hill-type muscle model.

Figure 2 .
Figure 2. (a) Schematic representation of muscle dynamics: inputs and outputs for activation and contraction dynamics.(b) Activation dynamics.(c) Contraction dynamics including the effects of the parallel elastic element.

Figure 3 .
Figure 3.Comparison of the different formulations for (a) Tendon and (b) parallel elastic element.Green solid line: Van Soest and Bobbert.Blue solid line: Kaplan-Silva.Red solid line: Thelen.

Figure 4 .
Figure 4. Comparison of the different formulations.(a) Force-length relationship.(b) Force-velocity relationship.Green solid line: Van Soest and Bobbert.Blue solid line: Kaplan-Silva.Red solid line: Thelen (c) Surface plots of the force-length-velocity relationship for the different formulations in the CE (Van Soest and Bobbert, Kaplan-Silva and Thelen, respectively.

Figure 5 .
Figure 5. (a) Variability in muscle length during tests.(b) Variability in the recorded sEMG during tests.Mean ± standard deviation.The normalized muscle lengths were obtained as the ratio l M /l M 0 .The sEMG signals were normalized to the maximum voluntary contraction (MVC).

Figure 6 .
Figure 6.Normalized muscle forces of biceps brachii in the flexoextension movement.Mean ± standard deviation.(a) PE force.(b) CE force.(c) CE + PE force.Green solid line: Van Soest and Bobbert.Red dotted line: Thelen.Blue dashed line: Silva-Kaplan.The normalized muscle forces were obtained as the ratio F PE/CE /F M 0 .

Table 1 .
Schematic summary of the different sections addressed.