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

Abstract

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.

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/ doi:10.5194/ms-7-19-2016 Author(s) 2016. CC Attribution 3.0 License.

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

F. Romero and F. J. Alonso

Department of Mechanical, Energy and Materials Engineering, University of Extremadura, Avda. de Elvas s/n, 06006 Badajoz, Spain

Correspondence to: F. Romero ([email protected])

Received: 3 September 2015 Revised: 4 December 2015 Accepted: 19 December 2015 Published: 18 January 2016

Abstract. 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 exo-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-exion. The results are expected to contribute in the selection of the different formulations of Hill-type muscle model to solve a specic problem.

1 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. (2007, 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 nd 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), modies the Voight model to include an additional spring in series with the DE. The dif-

Published by Copernicus Publications.

20 F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation

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.

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 Bob-bert, 1993; Van Den Bogert et al., 1998; Thelen, 2003; Silva, 2003; Ackermann and Schiehlen, 2006; Ackermann, 2007; Garca-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 bers. 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 difculty of the problem. However, as pointed out by Yamaguchi (2001), tendon elasticity is important if the tendon stretches an amount approaching the ber 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 rst approximation, it is possible to consider the tendon as a rigid element, and include as the tendon length its slack length (lTslack), that

is, the length for which the tendon just begin to resist lengthening.

The models available in the literature, based on the Hill-type description, differ basically in the amount of parameters to dene the muscle mechanical behaviour. An extensive review of Hill-based muscle models was performed by Winters (1990b). In the present work, we analyse two parameter-based models derived from the Hill-type description of muscle tissue. Specically, 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 specic 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.

2 Methods

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

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/

F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation 21

Table 1. Schematic summary of the different sections addressed.

Formulation of muscles dynamic equations

Activation dynamics Contraction dynamics

Mechanical expressions of muscle elements

Tendon Muscle

Parallel element Contractile element

Forcevelocity Forcelength

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.

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

a = f (u,a) (1) f MT = g a,lMT,vMT,f MT

[parenrightbig]

will describe the different dynamics that govern muscle activity and the contribution of the different tissues to the force development.

2.1.1 Excitation-to-activation dynamics

Excitation-to-activation dynamics or, simply, activation dynamics (Eq. 1), represents the muscle bers recruitment state.The most widely used expression to obtain muscle activation from a set of neural excitation is (Nagano and Gerritsen, 2001):

a(t) = (u(t) a(t)) (t1u(t) + t2), (3)

where t2 = 1/td and t1 = 1/(tat2) are time constants, u(t) is

the neural excitation 0 u 1, a(t) is the muscle activation,

where 0 u 1, ta is a time activation constant and td is

the deactivation constant. The values for those constants are taken from Umberger et al. (2003).

www.mech-sci.net/7/19/2016/ Mech. Sci., 7, 1929, 2016

(2)

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) denes the force exerted by the muscle as a function of the physiological state, that is, dependent of the activation level (a), the muscle length (lMT), the contraction velocity (vMT =

lMT) and the actual force (f MT). The following sections

22 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 bers 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:

aj (t) =

eAuj(t) 1 eA 1

, (4)

where aj (t) is the activation of the jth muscle, uj (t) the processed sEMG signal of the jth 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).

2.1.2 Activation-to-force dynamics

The force exerted by a muscle depends on the number of muscle bers recruited to perform a given movement (activation level), as well as on the actual ber 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:

F MT = F T = F M, (5)

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:

F MT = F M0 f T lT

[parenrightbig]

evM = vM/vmax, where vmax = lM0/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, specically the models proposed by Van Soest and Bobbert (1993), Thelen (2003) and Kaplan (2000)Silva (2003).

2.2 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:

F T "T

[parenrightbig]

= F TM0 fT "T

[parenrightbig]

F M

0 0.10377

=

[braceleftBigg]

e91"T 1[parenrightBig]if 0 "

T < 0.01516,

F M

0 37.526"

if 0.01516 "

T < 0.1 (8)

T

0.26029

[parenrightbig]

where "T = l

T

lslack

lslack .

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 dened by:

fT lT

0 if lT < lslack

= kT lT lslack

2 if lT lslack

(9)

2.

KaplanSilvas 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 dened 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.

2.3 Muscle

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

2.3.1 Parallel element

The parallel element represents the elasticity of the tissue attached in series to the muscle. According to Thelen (2003),

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/

where kT = F M0/ "T lslack

= F PE + F M

[parenrightbig]

cos p (6)

where F M0 is the maximum isometric force, f T(lT) is the normalized force exerted by the tendon that depends on its length (lT). 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):

F M = F M0 a fL

[tildewide]lM[parenrightbig]

, (7)

where a is the muscle activation and fL(

elM) and fV(

evM) the forcelength and forcevelocity relationships, respectively. These relationships depend on the normalized muscle length

elM = lM/lM0 and normalized contraction velocity,

fV

[tildewide]vM[parenrightbig]

F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation 23

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: KaplanSilva. Red solid line: Thelen.

the forcelength relationship for the parallel element is dened by:

FPE = F M0 fPE

[tildewide]lM[parenrightbig]

e kPE

[tildewide]lM[parenrightbig]/"M

[parenrightbig]

fV

(13)

The mathematical expressions for both relationships are detailed below.

Forcelength relationship

According to Thelen (2003), the forcelength relationship is described by:

fL

[tildewide]lM[parenrightbig]

= F M0

[tildewide]vM[parenrightbig]

0 1 , (10)

where kPE is a shape factor and "M0 is the parallel element strain. In this work, kPE = 3.0 and "M0 = 0.5.

Silva (2003), based on the work of Kaplan (2000), proposed an analytical expression for the PE:

ekPE 1

8 >

>

>

>

<

>

>

>

>

:

0 for lMT lM0

0 for

elM 1

FPE =

8F M0 lM30

3 for lM0 lM 1.63lM0

2F M0 for lM > 1.63lM0

lM lM0

=

8 <

:

e

[tildewide]lM1[parenrightbig]2/

for

(11)

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:

FPE = F M0 fPE(

elM) = [braceleftBigg]

0 for

elM > 1 , (14)

where represents the half-width of the curve for fL = 1/e.

The value of is set to 0.45 in this work.

Silva (2003), again, proposes an analytical expression for the forcelength relationship fL:

fL

[tildewide]lM[parenrightbig]

= e

[bracketleftbigg][bracketleftBig] 94 [parenleftBig][tildewide]

lM 1920 [parenrightBig][bracketrightBig]4 14 [bracketleftBig] 94 [parenleftBig][tildewide]

elM < 1

F M0 4

[tildewide]lM

elM 1(12)

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

elM 2 [0,1.6]. Values above elM = 1.6 present signicant divergences, however values over this limit are rarely reached in normal movements.

2.3.2 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 forcelength and forcevelocity relationships, that is:

F M = F M0 a fL

[tildewide]lM[parenrightbig]

lM 1920 [parenrightBig][bracketrightBig]2[bracketrightbigg] (15)

Contrariwise, Van Soest and Bobbert (1993) use directly a forcelengthvelocity curve. The expressions for the force length and forcevelocity relationships can be derived from it by setting

evM = 0 and

elM = 1 respectively (for a detailed explanation see Ackermann, 2007). Both curves are depicted in Fig. 4a and b. The proposed relationship is:

Concentric contraction vM < 0

f M

[tildewide]lM,[tildewide]vM[parenrightbig]

1

2 for

Br (fiso + Ar) Ar

: Br [tildewide]vMfac [parenrightBig]

F M0 =

, (16)

where a is the muscle activation, fac = min(1,3.33 a)

and, Ar and Br are Hills constants corrections according to Van Soest and Bobbert (1993). The values for

www.mech-sci.net/7/19/2016/ Mech. Sci., 7, 1929, 2016

a

Br [tildewide]vMfac

24 F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation

Figure 4. Comparison of the different formulations. (a) Forcelength relationship. (b) Forcevelocity relationship. Green solid line: Van Soest and Bobbert. Blue solid line: KaplanSilva. Red solid line: Thelen (c) Surface plots of the forcelengthvelocity relationship for the different formulations in the CE (Van Soest and Bobbert, KaplanSilva and Thelen, respectively.

those constants are set to Ar = 0.41 and Br = 5.2. Fi

nally, the isometric force relative to the maximum force, fiso can be obtained as:

fiso =

1 elM2 +

1

2 elM

width2 + 1, (17)

were width is the maximum range of force production relative to lM0.

Eccentric contraction vM > 0

:

f M

[tildewide]lM,[tildewide]vM[parenrightbig]

where f maxv is the maximum normalized achievable muscle force when the muscle is lengthening (in this work this value is set to 1.6), fac and Sf are shape factors usually set to 1 to reduce computational effort.

The comparison of the different forcelength relationship formulations are represented in Fig. 4b.

Forcevelocity relationship

According to Thelen (2003), the forcevelocity relationship can be written as:

evM

evM , (18)

where shape factors b1, b2 and b3 can be obtained as:

b2 = fisof maxv

b1 =

facBr(fiso + b2)2 (fiso + Ar)slope factor

F M0 =

b1 b2 b3

0 for

evM 1,

[tildewide]vM[parenrightbig]

8 >

>

>

>

>

<

>

>

>

>

>

:

a

b 3

evM 1

evM/kCE1

1 +

for 1 <

evM 0,

fV

=

1 +

evMf maxv /kCE2 1 +

evm/kCE2 for

evM > 0

, (20)

where kCE1 and kCE2 are forcevelocity shape factors (0.25 and 0.06 respectively in this work)and f maxv is the maximum normalized achievable muscle force when the muscle is lengthening (in this work this value is set to 1.6).

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/

(19)

b3 =

b1 fiso + b2

,

F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation 25

The analytical expression given by Silva (2003) is dened as:

8 >

>

>

>

<

>

>

>

>

:

0 for

evM < 1

[tildewide]vM[parenrightbig]

arctan 5

evM

arctan(5) + 1 for 1

evM 0.2

. (21)

The expression given by Van Soest and Bobbert (1993) was already addressed on the previous section. The comparison of the different curves for the forcevelocity relationship is represented in Fig. 4b. Besides, it is possible to represent in a surface plot the forcelengthvelocity 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 forcelength relationship described by Silva (2003) is not centred on

elM = 1 and the bell shape is slenderer than the other approaches. However, the forcevelocity relationship is quite similar to the model proposed by Thelen (2003). In this case, the model proposed by Van Soest and Bobbert (1993) presents signicant differences in the interval

evM 2 [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 (lM0) the passive force generated by the PE becomes signicant.

3 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 exo extension movement under load was performed. The recordings consist of the acquisition of the kinematics of reective 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 exoextension 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 subjects dimensions by using the recorded positions of the reective markers. Parameters such as lM0 or lslack were obtained after the scale. Others, as F M0, 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 subjects 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 exion 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 reective 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 supercial 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 ltered by means of singular spectrum analysis (see Romero et al., 2015 for further details) and the rst 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:

M =

[radicalbigg]

fV

=

4arctan(5) + 1 for

evM > 0.2

vcc

vmm 1,

P =

1[pi1] cos1

vmc

pvmmvcc , (22)

pM2 + P 2,

where the values vmm, vcc and vmc can be dened as:

vmm =

1t2 t1

C =

t2

[integraldisplay]

t1

m(t)2dt,

vcc =

1t2 t1

t2

[integraldisplay]

t1

c(t)2dt, (23)

www.mech-sci.net/7/19/2016/ Mech. Sci., 7, 1929, 2016

26 F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation

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 lM/lM

0 . The sEMG signals were normalized to the maximum voluntary contraction

(MVC).

Regarding to muscle elements (see Figs. 3b, 4), the three proposed models show qualitative similarities. For the parallel element, a notable difference between SilvaKaplans model and the rest for values below

elM = 1.5 (see Fig. 3b).

As the range of normalized muscle length of the recorded movement corresponds to

elM 2 [0.5,1.2] (Fig. 5a) signi-cant 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 forcelength 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 forcelength relationship (Fig. 4a), the differences in the range 3560 % are the ones expected. In fact, similar values are observed in the neighbouring of the optimal length, where all the forcelength curves present similar values. The major differences are observed at 50 % of the cycle, corresponding to full exion. These variations correspond to the differences in the bell-shaped curve of the forcelength relationship, as for full-exion values (

elM 0.5) there are considerable differences in the forcelength relationship. SilvaKaplans 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 SilvaKaplans models present major similarities between them compared to Thelen and Soest and Bobberts. These differences are more evident during the extension phase with deviations of 40.68 and 28.38 % between Soest and Bobberts models and those of Thelen or SilvaKaplan, respectively.Differences between curves are lower during exion in both

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/

vmc =

1t2 t1

t2

[integraldisplay]

t1

m(t)c(t)dt,

where m(t) and c(t) are the measured and the computed signals and [t1,t2] 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 dened as the interest signal, which is compared with the signal c(t). Major similarities between signals are related with indices near to zero.

4 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 gure, 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 bers 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 ber length of a particular muscle, and therefore it must be considered in these cases (Zatsiorsky and Prilutsky, 2012).

F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation 27

(a)Magnitudecomparison.Completecycle.(b)Phasecomparison.Completecycle.(c)Combinederror.Completecycle.

SoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelen

SoestandBobbert00.20350.154500.12520.066900.23890.1683

SilvaKaplan0.255500.06150.125200.06020.284600.0861

Thelen0.18270.058000.06690.060200.19460.08360

(d)Magnitudecomparison.Flexion.(e)Phasecomparison.Flexion.(f)Combinederror.Flexion

SoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelen

SoestandBobbert00.17640.133000.11270.062000.20930.1467

SilvaKaplan0.214200.05260.112700.05260.242000.0745

Thelen0.15340.050000.06200.052600.16550.07260

(g)Magnitudecomparison.Extension.(h)Phasecomparison.Extension.(i)Combinederror.Extension.

SoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelenSoestandBobbertSilvaKaplanThelen

SoestandBobbert00.28910.221000.15560.075200.32840.2335

SilvaKaplan0.406800.09580.155600.08210.435500.1262

Thelen0.28380.087400.07520.082100.29360.11990

Figure 6. Normalized muscle forces of biceps brachii in the exo extension 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 .

Table2.ComparisonofdifferentformulationsbyusingGeersindex.

www.mech-sci.net/7/19/2016/ Mech. Sci., 7, 1929, 2016

tion. Thelens 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, SilvaKaplans model can be used in experiments in which computing time and computational effort are critical but not the use of physiological parameters whereas Thelens 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 specic conditions of the problem.

Acknowledgements. This work was supported by the Spanish Ministry of Economy and Competitiveness under project DPI2012-38331-C03, co-nanced by the European Union through EFRD funds.

Edited by: P. FloresReviewed by: two anonymous referees

References

Ackermann, M.: Dynamics and energetics of walking with pros-theses, PhD thesis, University of Stuttgart, Stuttgart, Germany, 2007.

Ackermann, M. and Schiehlen, W.: Dynamic analysis of human gait disorder and metabolical cost estimation, Arch. Appl. Mech., 75, 569594, 2006.

Alonso, J., Romero, F., Pmies-Vil, R., Lugrs, U., and Font-

Llagunes, J.: A simple approach to estimate muscle forces and orthosis actuation in powered assisted walking of spinal cord-injured subjects, Multibody Syst. Dyn., 28, 109124, 2012.Buchanan, T., Lloyd, D., Manal, K., and Besier, T.: Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command,J. Appl. Biomech., 20, 367395, 2004.

Buchanan, T., Lloyd, D., Manal, K., and Besier, T.: Estimation of muscle forces and joint moments using a forward-inverse dynamics model, Med. Sci. Sport. Exer., 37, 19111916, 2005.Criswell, E.: Crams introduction to surface electromyography,

Jones & Bartlett Publishers, Sudbury, USA, 2010.

Delp, S., Anderson, F., Arnold, A., Loan, P., Habib, A., John, C.,

Guendelman, E., and Thelen, D.: OpenSim: open-source software to create and analyze dynamic simulations of movement, IEEE T. Bio.-Med. Eng., 54, 19401950, 2007.Garca-Vallejo, D.: Simulacin de la marcha humana mediante optimizacin paramtrica, XVIII Congreso Nacional de Ingeniera Mecnica, Ciudad Real, 35 November 2010.

Geers, T. L.: An Objective Error Measure for the Comparison of Calculated and Measured Transient Response Histories, The Shock and Vibration Bulletin, The Shock and Vibration Information Center, Naval Research Laboratory, Washington, D.C., Bulletin 54, Part 2, 99107, June 1984.

Grahovac, N. and igi, M.: Modelling of the hamstring muscle group by use of fractional derivatives, Comput. Math. Appl., 59, 16951700, 2010.

Mech. Sci., 7, 1929, 2016 www.mech-sci.net/7/19/2016/

28 F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation

magnitude and phase as expressed by Geers indices. On the contrary, the major differences observed in the models are related to the full-exion state, where the models show certain divergences. At full-exion, Van Soest and Bobebrts model overestimates the force whereas KaplanSilva 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. KaplanSilvas model is based in the adjustment of Hills 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 Bobberts 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, KaplanSilva (0.021811 s) outperforms Thelen(0.023208 s) and Van Soest and Bobberts (0.026339 s) muscle descriptions. The differences may not be signicant 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 ltered properly. Nevertheless, the results presented here can be used as a reference to muscle force production during arm exo-extension movement and to evaluate differences between models prior their selection for biomechanical studies.

5 Conclusions

In this paper, the mechanical expressions to represent muscle tissue have been presented. To compute muscle efforts, the most signicant 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, SilvaKaplans model is based in analytical expressions, and from the computational point of view, their use is more efcient as the expressions are simpler than the other cases, however, for some specic applications it may be interesting to include physiological information which is not provided with this formula-

F. Romero and F. J. Alonso: Different Hill-type formulations for muscle force estimation 29

Hatze, H.: A myocybernetic control model of skeletal muscle, Biol.Cybern., 25, 103119, 1977.

Hermens, H., Freriks, B., Disselhorst-Klug, C., and Rau, G.: Development of recommendations for sEMG sensors and sensor placement procedures, J. Electromyogr. Kines., 10, 361374, 2000.Hill, A.: The Heat of Shortening and the Dynamic Constants of

Muscle, P. Roy. Soc. B-Biol. Sci., 126, 136195, 1938.

Holzbaur, K., Murray, W., and Delp, S.: A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control, Ann. Biomed. Eng., 33, 829840, 2005.HosseinNia, S., Romero, F., Tejado, I., Vinagre, B., and Alonso,J.: Effects of Introducing Fractional Dynamics in Hills Model for Muscle Contraction, in: System Identication, International Federation of Automatic Control, 16, 17431748, 2012.Kaplan, M.: Efcient optimal contro of large-scale biomechanical systems, PhD thesis, Stanford University, USA, 2000.

Konrad, P.: The ABC of EMG, in: A Practical Introduction to Kinesiological Electromyography, Version 1.4, Noraxon INC, USA, 2643, 2006.

Lund, M. E., de Zee, M., and Rasmussen, J.: Comparing calculated and measured curves in validation of musculoskeletal models, in: XIII International Symposium on Computer Simulation in Biomechanics, Leuven, Belgium, 30 June2 July 2011.Martins, J., Pires, E., Salvado, R., and Dinis, P.: A numerical model of passive and active behavior of skeletal muscles, Comput.Method. Appl. M., 151, 419433, 1998.

Nagano, A. and Gerritsen, K.: Effects of neuromuscular strength training on vertical jumping performance A computer simulation study, J. Appl. Biomech., 17, 113128, 2001.

Nigg, B. and Herzog, W.: Biomechanics of the musculo-skeletal system, vol. 192, Wiley New York, 1999.

Romero, F., Alonso, F., Cubero, J., and Galn-Marn, G.: An automatic SSA-based de-noising and smoothing technique for surface electromyography signals, Biomed. Signal Proces., 18, 317 324, 2015.

Sartori, M., Reggiani, M., Mezzato, C., and Pagello, E.: A lower limb EMG-driven biomechanical model for applications in rehabilitation robotics, in: IEEE International Conference on Advanced Robotics, ICAR 2009, Munich, Germany, 2226 June 2009, 17, 2009.

Silva, M.: Human motion analysis using multibody dynamics and optimiation tools, PhD thesis, Instituto Superior Tcnico, Unversidade Tcnica de Lisboa, Lisbon, Portugal, 2003.

Sommacal, L., Melchior, P., Cabelguen, J.-M., Oustaloup, A., and Ijspeert, A.: Fractional multimodels of the gastrocnemius muscle for tetanus pattern, in: Advances in fractional calculus, Springer, 271285, 2007.

Sommacal, L., Melchior, P., Oustaloup, A., Cabelguen, J.-M., and Ijspeert, A. J.: Fractional multi-models of the frog gastrocnemius muscle, J. Vib. Control, 14, 14151430, 2008.

Thelen, D.: Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults, J. Biomech. Eng.-T. ASME, 125, 7077, 2003.

Umberger, B., Gerritsen, K., and Martin, P.: A model of human muscle energy expenditure, Comput. Method. Biomec., 6, 99111, 2003.

Van Den Bogert, A., Gerritsen, K., and Cole, G.: Human muscle modelling from a users perspective, J. Electromyogr. Kines., 8, 119124, 1998.

Van Soest, A. and Bobbert, M.: The contribution of muscle properties in the control of explosive movements, Biol. Cybern., 69, 195204, 1993.

Winters, J. M.: Hill-based muscle models: a systems engineering perspective, Springer-Verlag, 6993, 1990a.

Winters, J. M.: Hill-based muscle models: a systems engineering perspective, in: Multiple muscle systems, Springer, 6993, 1990b.

Winters, J. M. and Stark, L.: Muscle models: what is gained and what is lost by varying model complexity, Biol. Cybern., 55, 403420, 1987.

Yamaguchi, G.: Dynamic modeling of musculoskeletal motion: a vectorized approach for biomechanical analysis in three dimensions, Kluwer Academic Publishers Norwell, MA, 2001.

Zajac, F.: Muscle and tendon: Properties, models, scaling and applications to biomechanics and motor control, Crit. Rev. Biomed. Eng., 17, 359411, 1989.

Zatsiorsky, V. and Prilutsky, B.: Biomechanics of skeletal muscles, Human Kinetics, Human Kinetics Publisher, Champaign, USA, 2012.

www.mech-sci.net/7/19/2016/ Mech. Sci., 7, 1929, 2016