This article was downloaded by: [Cukurova Universitesi] On: 05 November 2014, At: 06:02 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Computer Methods in Biomechanics and Biomedical Engineering Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/gcmb20

Surface EMG-force modelling for the biceps brachii and its experimental evaluation during isometric isotonic contractions a

b

b

b

Hua Cao , Sofiane Boudaoud , Frédéric Marin & Catherine Marque a

Département EEA, Haute Etude de l'Ingénieur (HEI), Lille, France

b

UMR CNRS 7338: Biomécanique et Bioingénierie, Centre de Recherche Royallieu, Université de Technologie de Compiègne, BP 20529, 60205, Compiègne, France Published online: 27 Jan 2014.

To cite this article: Hua Cao, Sofiane Boudaoud, Frédéric Marin & Catherine Marque (2014): Surface EMG-force modelling for the biceps brachii and its experimental evaluation during isometric isotonic contractions, Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255842.2013.867952 To link to this article: http://dx.doi.org/10.1080/10255842.2013.867952

PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions

Computer Methods in Biomechanics and Biomedical Engineering, 2014 http://dx.doi.org/10.1080/10255842.2013.867952

Surface EMG-force modelling for the biceps brachii and its experimental evaluation during isometric isotonic contractions Hua Caoa, Sofiane Boudaoudb*, Fre´de´ric Marinb and Catherine Marqueb a

De´partement EEA, Haute Etude de l’Inge´nieur (HEI), Lille, France; bUMR CNRS 7338: Biome´canique et Bioinge´nierie, Centre de Recherche Royallieu, Universite´ de Technologie de Compie`gne, BP 20529, 60205, Compie`gne, France

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

(Received 15 October 2012; accepted 18 November 2013) The aim of this study was to evaluate a surface electromyography (sEMG) signal and force model for the biceps brachii muscle during isotonic isometric contractions for an experimental set-up as well as for a simulation. The proposed model includes a new rate coding scheme and a new analytical formulation of the muscle force generation. The proposed rate coding scheme supposes varying minimum and peak firing frequencies according to motor unit (MU) type (I or II). Practically, the proposed analytical mechanogram allows us to tune the force contribution of each active MU according to its type and instantaneous firing rate. A subsequent sensitivity analysis using a Monte Carlo simulation allows deducing optimised input parameter ranges that guarantee a realistic behaviour of the proposed model according to two existing criteria and an additional one. In fact, this proposed new criterion evaluates the force generation efficiency according to neural intent. Experiments and simulations, at varying force levels and using the optimised parameter ranges, were performed to evaluate the proposed model. As a result, our study showed that the proposed sEMG– force modelling can emulate the biceps brachii behaviour during isotonic isometric contractions. Keywords: biomechanics; sEMG; muscle force; modelling; motor unit recruitment scheme; isometric contractions; Monte-Carlo simulations; sensitivity analysis

1. Introduction Skeletal muscles generate active forces to move or stabilise the skeleton in association with human activities (Basmajian and De Luca 1985). The force exerted by a skeletal muscle cannot be directly measured noninvasively (Disselhorst-Klug et al. 2009; Staudenmann et al. 2010). Estimations based on mechanical Newton’s equations are also controversial due to the redundancy of the muscle function, which generates an indeterminate equation system (Disselhorst-Klug et al. 2009; Staudenmann et al. 2010). Consequently, any attempt for muscle force estimation based on neuro-musculoskeletal modelling cannot be validated using this paradigm (Lloyd and Besier 2003; Buchanan et al. 2004). From a physiological point of view, the production of mechanical force is always associated with an electrical activity which can be measured by surface electromyography (sEMG) (De Luca et al. 1996). The sEMG signal is an objective and non-invasive technique to assess the muscle activation in human voluntary contractions (Hagg et al. 2000). Several methods are proposed to correlate the sEMG signal and the level of external force produced by muscle contractions (Staudenmann et al. 2010). However, all methods are based on the hypothesis that sEMG and muscle force are directly related, notwithstanding the fact that this hypothesis has not yet been clearly demonstrated

*Corresponding author. Email: [email protected] q 2014 Taylor & Francis

(Fuglevand et al. 1993; Lloyd and Besier 2003; Disselhorst-Klug et al. 2009). From the motor control theory point of view, the force generated by a muscle is controlled by the central nervous system (CNS), depending on two parameters: the recruitment of motor units (MUs) and the firing rate of active MUs (Erim et al. 1996). These two parameters are directly connected with the generation of electrical activity inside the muscle and also influence the sEMG signal (Fuglevand et al. 1993). Consequently on the sEMG/force relationship, on the one hand side, the amplitude of the sEMG signal is affected by the shape of action potentials and the position of active MUs that do not affect the production of force (Farina et al. 2004), and on the other, instantaneous muscle length, rate of change in length (Hof 1997), and contraction history (Welter and Bobbert 2001) affects the force at a given level of muscle activity that is not reflected in the sEMG signal. Consequently, the relationship between the sEMG and the muscle force production is not trivial and highly dependent on the studied muscle and the measurement protocol (Woods and Bigland-Ritchie 1983; Fuglevand et al. 1993; Farina et al. 2002; Zhou and Rymer 2004; Disselhorst-Klug et al. 2009). Numerical modelling is useful for linking the internal physiological parameters – e.g. MU recruitment scheme, conduction velocity (CV), and generated force – with the

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

2

H. Cao et al.

sEMG signal generation process (Farina and Merletti 2001). Moreover, the influence of the physiological parameters on the EMG signal can be analysed by simulation (Farina et al. 2002; Zhou and Rymer 2004). Fuglevand et al. (1993) have developed a model to estimate the EMG signal and the force generated by a muscle, with a given target isometric contraction according to recruitment /rate coding of the active MUs. The force is generated by accumulating individual MU twitches (Fuglevand et al. 1993). Several activation patterns were tested to analyse the influence of the activation pattern on the sEMG – force relationship. This model was then applied in physiological and pathological conditions (Zhou and Rymer 2004; Zhou et al. 2007). The ranges of input parameters of this model were optimised by using a Monte Carlo method to test its sensitivity to parameter variations (Keenan and Valero-Cuevas 2007). However, in these works, the used sEMG – force relationship model deals with isometric contractions generated by the same Fuglevand’s MU recruitment scheme (Fuglevand et al. 1993). Based on this approach, the important role played by the neural parameters in driving the sEMG-force relationship has been observed (Keenan and Valero-Cuevas 2007). In addition, to better describe the muscle force generation, another model has been proposed. It included a ‘twitch’ force model associated with a basic MU recruitment scheme, based on needle experimental results, simple conduction volume and constant MU firing rates (Wheeler et al. 2010). However, at this stage, the realism of the MU recruitment schemes remains questionable. The objective of this study was to define and evaluate an sEMG – force model during isometric isotonic contractions of the biceps brachii. The novelty of the proposed model concerns the proposition of a new MU rate coding scheme and a new analytical mechanogram. It is assessed by an evaluation of optimal model parameter ranges coupled with a sensitivity analysis. Finally, the model will be compared to an in vivo experimental set-up recording external force and sEMG activities of the biceps brachii muscle. 2. 2.1

Figure 1. This schematic diagram provides the number of MUs activated and their firing rates for a given muscle force command P(t), expressed in percentage of the maximal voluntary contraction (MVC). It includes two MU types [slow motor units (SMUs) and fast motor units (FMUs)], with different number of muscle fibres. Each MU type has its specific minimum firing rate (MFRF and MFRs) and peak firing rate (PFRF and PFRs). According to the fact that motoneurons innervating fast twitch muscle fibres exhibit shorter after-hyperpolarisations than those innervating slow twitch muscle fibres (Gardiner 1993; Christie and Kamen 2010), we defined MFRF . MFRs and PFRF . PFRs. The firing rate of an active MU increases linearly from its MFR to its PFR with muscle force increase (Erim et al. 1996). For human muscles, the MFR ranges from 7 to 23 Hz, and the PFR from 14 to 50 Hz (Erim et al. 1996; Moritz et al. 2005). All MUs of the same type reach their PFR at 100% MVC. MUs are progressively recruited in an orderly sequence. According to the size principle, SMUs are activated first and then the FMUs (Henneman 1980). The whole MUs are recruited at a given force (RR) following the Fuglevand model (Fuglevand et al. 1993). RR ranges between 30% and 90% of MVC according to the muscle type (Kukulka and Clamann 1981). The recruitment threshold (RT) is linked to the requested neural excitatory current at a specific force level and needs to recruit a new MU. RTi of the ith MU is modelled as follows (Fuglevand et al. 1993): RTi ¼ eðlnRR=NÞi ;

ð1Þ

where RR is the recruitment range, i.e. the RT of the last activated MU in a muscle (30 – 90% MVC in our test), N is the number of MUs and i the MU index. The second model component deals with the sEMG signal generation. The sEMG is the result of the electrical activity propagating along the fibre inside the muscle,

Materials and methods sEMG – force modelling

sEMG – force modelling can be decomposed in three main components (Fuglevand et al. 1993; Keenan and ValeroCuevas 2007): (1) the spatio-temporal recruitment pattern of the active MUs, (2) the sEMG signal generation and (3) the muscle force simulation according to MU recruitment. All these components are driven by the input force command profile originating from the CNS (Fuglevand et al. 1993). The first component is the spatio-temporal MU recruitment pattern (Cao et al. 2008). In our model, it is generated according to the proposed relation depicted in

Figure 1. MU recruitment pattern (PFRF, MFRF, PFRS, MFRS: peak and MFRs of FMUs and SMUs, respectively; RT1 ¼ 0). In this schematic diagram, four SMUs (solid lines) and three FMUs (dashed lines) are recruited for the given muscle input force level P(t).

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

Computer Methods in Biomechanics and Biomedical Engineering measured at the skin surface, through the tissues located between the muscle and the electrodes (conducting volume). Consequently, in order to simulate the sEMG signal, it is necessary to have a geometrical/anatomical model of the muscle, a propagation model for the electrical activity of the active fibres through the soft tissues (conducting volume), then a source generation model including the generation, propagation and extinction of the propagating electrical sources along the fibres (action potentials). For the geometrical model of the anatomy of biceps brachii, the MUs are uniformly distributed within the muscle cross section. This assumption is a compromise due to the fact that some studies point out a possible spatial regionalisation of the FMUs near the muscle surface and the SMUs deeper in the biceps muscle suffering from important variability (age, sex and physical condition) (Jennekens et al. 1971). The muscle fibres are also located within each MU with a uniform spatial distribution. Due to the variability of the human muscle anatomy (Maclntosh et al. 2006), the number of FMUs and of SMUs in a muscle can vary from 250 to 600 units. The number of fibres per MU can vary from 30 to 100 units (Guyton and Hall 2000). Moreover, the muscle fibre length is supposed to vary from 4 to 16 cm (Keenan and Valero-Cuevas 2007). For the source modelling, the Rosenfalck model (Andreassen and Rosenfalck 1978) has been used. The centre of the innervation zones is located in the middle of the fibres. However, to take into account the neuromuscular junction (NMJ) geometric dispersion, for a given MU, we distributed a uniform time delay (from 0 to 5 ms) for the starting time of action potentials among the MU fibres, at the innervation zone that corresponds to an uniform spatial distribution of the fibre innervation admitted in previous models (Farina and Merletti 2001; Keenan and Valero-Cuevas 2007). The muscle fibre CV is considered as a Gaussian distribution with a mean value of 3 –4 m/s and a standard deviation of 0.5 m/s (Farina and Merletti 2001). For a given muscle, the fast twitch fibre CV is greater than for the slow twitch fibre (Farina and Merletti 2001; Keenan and Valero-Cuevas 2007). The inter pulse interval (IPI) distribution of the MUs firing is modelled as a Gaussian distribution function with a standard deviation between 10% and 30% of the mean IPI. We compute the motor unit action potential train (MUAPT) of each MU, according to the MU recruitment/ firing pattern. Concerning the propagation through the conducting volume, the potential distribution on the skin surface, named output plan, is the sum of the MUAPT of each MU filtered by a spatial transfer function depending on the thickness of the volume conductor (muscle, fat and skin layers) (Farina and Merletti 2001). The simulated sEMG signal is finally obtained by filtering the output potential plan, on the skin, with an electrode transfer function (Farina and Merletti 2001). In this study, single differential

3

arrangement of the electrodes (electrode diameter, 10 mm and inter-electrode distance, 20 mm) will be considered for both simulation and experimentation. They will be positioned near the belly of the muscle at a middle distance between the NMJ and the elbow (Barbero et al. 2012). The last model component is also used to compute the muscle force resulting from the MUs recruitment. For a given P(t) force input command defined as the requested force driven by the CNS trough electrical neural excitation, the force generated by the muscle Fm(t) is calculated by summing the A active MU forces (FMUi). To minimise the time computation, for the ith active MU, FMUi is supposed to increase linearly with a positive slope ai according to the firing rate FRi (Cao et al. 2009) obtained, as depicted on Figure 1. In addition, the force variability is simulated by a zero-mean Gaussian noise N(0, sP(t)) with varying standard deviation according to P(t).

F m ðtÞ ¼

A X i¼1

F MUi;PðtÞ ¼

A X

ai ·FRi;PðtÞ þ Nð0; sPðtÞ Þ: ð2Þ

i¼1

Let ai represent the slope associated with the ith MU, FRij (i [ [1, N ], j [ [1, N ]) is its firing rate at the force level RTj and A is the total number of MUs. To have a determinist identification of ai, let A ¼ N and sP(t) ¼ 0. Consequently, the N MU of the muscle generated the force Fm(t). Based on Equation (2) at each recruitment threshold (RTi), the values of Fm(t) are related with the FR of i muscles. The following equations system related to the MU recruitment pattern could be deduced as follows: 8 > > > > > > > > > > > >
> > > > > > > > > > > : PFR1 ·a1 þ PFR2 ·a2 þ · · · þ PFRN ·aN ¼ 100%MVC ð3Þ with FRij ¼ ððPFRi 2 MFRi Þ=ð100 2 RTi ÞÞ·RTj . Solving this equation system, a unique set of ai values is determined. Finally, the force Fm(t), corresponding to the input command P(t), is given by the following equation from (2): F M ðtÞ ¼ FR1;PðtÞ ·a1 þ FR2;PðtÞ ·a2 þ · · · þ FRA;PðtÞ ·aA ð4Þ

4

H. Cao et al.

with

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

FRi;PðtÞ ¼

PFRi 2 MFRi ·PðtÞ; 100 2 RTi

where A is the index of the last recruited MU for the input command P(t) as described in (2). The signal P(t) can be considered as the attended force profile by the CNS. This injunction is physiologically translated by the synaptic neural current intensity that flows through motor neurons and tuned both MU recruitment and rate coding. Since the generated force Fm(t), given by (4), is constant for a constant given P(t) input command, a Gaussian noise has been added to mimic the force variability during constant isometric contraction included in (2). In fact, its standard deviation has been chosen to increase linearly with increasing contraction level. At this stage, the novelty of our model, when compared to previous works (Fuglevand et al. 1993; Cheng et al. 2000) is the possibility to associate to each active MU a different positive slope for the linear regression force/firing rate in relation with type (slow or rapid) and the number N of MU during isometric contraction. This approach is supported by the Monster and Chan’s study on the extensor digitorum communis muscle that evidenced different slopes from experimental data (Monster and Chan 1977). Furthermore, including more realistic MU firing behaviour may offer a higher neural drive realism to the sEMG/force simulation (Keenan and Valero-Cuevas 2007). A schematic graph of the proposed sEMG/force generation model is shown in Figure 2. 2.2 Model sensitivity evaluation and confidence intervals estimation Since many parameters can interact in the sEMG-force relationship, it is necessary to evaluate their sensitivity (Keenan and Valero-Cuevas 2007). The proposed sEMG –

Figure 2.

force model requires 12 physiological parameter values coming from the literature (Fuglevand et al. 1993; Keenan and Valero-Cuevas 2007). However, the variability and uncertainty of the parameters used in these models have increased the difficulty of their assessment (Keenan and Valero-Cuevas 2007). Another critical point is the estimation of the confidence intervals for these parameters that will permit a realistic behaviour of the simulated sEMG and force signals. To analyse multivariate interactions of these parameters, the Monte Carlo method demonstrated the ability to evaluate the sensitivity of a sEMG – force model, based on a classical MU recruitment model during isometric contraction (Keenan and Valero-Cuevas 2007). The same method is performed here. For each Monte Carlo simulation, a set of parameter values has been defined by random sampling from uniform distributions, for each one of the 12 parameter ranges. The simulated outputs (sEMG and muscle force Fm(t)) were computed at 1000 samples/s. They used two criteria, namely sEMG/force regression slope and force variability/ force (Keenan and Valero-Cuevas 2007), to evaluate the confidence interval of each parameter. The sEMG signal was rectified and smoothed by a four-order low-pass Butterworth filter with a cut-off frequency of 8 Hz to obtain its envelope (Buchanan et al. 2004). For the sensitivity analysis, the convergence is achieved if the running coefficients of multiple correlation of the EMG envelope and of the force are lower than 2% for the last 20% of the simulations at all four command levels (Santos and Valero-Cuevas 2006; Keenan and Valero-Cuevas 2007). The sensitivity of the model has been analysed by counting the number of Monte Carlo simulations required to achieve the convergence test (Keenan and ValeroCuevas 2007). For our sensitivity analysis, the same approach has been used adding a new criterion to compute the confidence intervals of the various parameters. This criterion is obtained by estimating the similarity between the input (or target) force P(t) and the output muscle force

A schematic graph of the proposed sEMG/force generation model with the description of the principal parameter settings.

Computer Methods in Biomechanics and Biomedical Engineering

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

FM(t) obtained. It is estimated by the root mean square error (ErrorRMS) between P(t) and FM(t). This criterion will be compared (see Section 3.1) with the ones used by Keenan and Valero-Cuevas (2007). Furthermore, we also analysed the obtained confident intervals. These confident intervals are the results of a refining process. In fact, two series (Series 1 and Series 2) of simulations were launched following the methodology described before. The Series 1 roughly delimits the confident intervals while the Series 2, initialised by the results of the Series 1, refines these confident intervals. Only the sets of parameters that lay within the confidence intervals obtained from the Series 2 were used for the comparison with the experimental data. 2.3 Experimental protocol To assess the degree of realism of the model, we compared the simulated sEMG and force signals with the experimental data. We expected the best correlation between biceps brachii force and the hand pull force, for an elbow flexion at 908 in an isometric situation. The elbow angle has been controlled by a goniometer (Figure 3). In these conditions, the intrinsic Biceps force behaviour should be linearly linked to the measured hand pull force. Ten healthy subjects have been involved in this experimental evaluation (five males and five females; age ¼ 25.8 ^ 2.0 years, height ¼ 172.6 ^ 9.0 cm, weight ¼ 68.1 ^ 13.6 kg). The study has been carried out in accordance with the Declaration of Helsinki and after approval of the local ethic comity. None of the volunteers were regular exercise practitioners. They did not present any symptoms of neuromuscular or ligament disorders, and were not using anti-inflammatory medication or muscle relaxants during the experimental sessions. Written informed consent had been given by all subjects prior to the sessions. After preparing the skin of the subject to lower the inter-electrode impedance (less than 10 kV), the sEMG activity was recorded using a telemetry acquisition system. We used bipolar electrode configuration including two Ag/AgCl electrodes (Noraxon dual electrodes, Noraxon Inc., Scottsdale, AZ, USA), 10 mm in diameter, with an interelectrode distance of 20 mm, placed along the longitudinal axis of the biceps of the subject’s right arm between the middle of the muscle belly and the elbow as in simulated

Figure 3.

Description of the experimental protocol.

5

configurations. A reference electrode was also placed over the bony part of the elbow (olecranon). The sEMG signals of the biceps brachii were recorded via a portable transmitter system (TeleMyo 2400T G2 transmitter, Noraxon), combined with a receiver (TeleMyo 2400R G2 receiver, Noraxon). The subject is in a supine position with an elbow angle of 908 (see Figure 3). By a handle-cablepulley set-up, connected to a force transducer, the isometric contraction mechanical action was recorded. An oscilloscope connected to this force transducer was used as visual feedback for the subject. To control the isometric behaviour of the muscle contraction, a flexible goniometer (Biometrics LTD, Newport, UK) was attached to the surface of the arm and the forearm along the bone permitting to measure the elbow angle. An elbow angle of 908 ^ 5% was accepted in our exercise. During the measurement sessions, subjects were asked to maintain the wrist static, hand closed aligned with the forearm (Figure 3) and to keep their triceps and shoulder in contact with the ergometer bed. These postural constraints were visually checked all through the experiments. No control of cross-muscle force was included in the protocol since it was supposed that, if the wrist and shoulder joints remained static in an optimal position, the measured force was to be majorly generated by the biceps brachii contraction. All signals were simultaneously recorded at a sampling frequency of 1000 Hz and recorded by the acquisition system (CMS-HS system, Zebris GmbH, Isny im Allga¨u, Germany). The experimental protocol included two successive exercises. First, the subject performed three maximal arm pulls of 5 s, separated by 3-min rest periods. The highest force measured defined the maximal voluntary isometric contractions (MVC). After a short training, and following the same protocol for obtaining MVC values, three trials were performed on each level of force (20%, 50%, 80%, and 100% MVC) and a rest interval of 3 min was allowed between trials or force levels. To compare experimental and simulated data, all parameters were relative to the force normalised by the MVC of each subject. For analysis purposes, a selective window has been applied to sEMG recordings to extract a sequence (1 s duration) in the middle of the recordings where the measured force is stabilised on the requested value as depicted in Figure 4. In parallel with the experimental data, the simulated sEMG and the muscle force of the biceps brachii were computed by using the proposed sEMG–force model for isotonic isometric contraction of 2 s with a maximum at 20%, 50%, 80%, and 100% MVC. To simulate the biceps brachii of healthy subjects, we respected the following hypotheses: all MUs are recruited at a force level of approximately 80% MVC (Fuglevand et al. 1993); the average CV of action potentials along the fibre is 3.62 m/s (Gydikov and Gantchev 1989); the average length of fibres is about 14 cm (An et al. 1989); the biceps muscle is composed of around 770 MUs

6

H. Cao et al.

P(t) (%MVC)

and the median frequency (MDF) obtained from the power spectrum density of the sEMG signals estimated by the Welch’s method (Welch 1967). A statistical analysis is then performed to compare the experimental results of the 10 subjects with data obtained from the 10 simulated muscles. Due to the studied sample population smallness, neither normality nor homoscedasticity has been supposed to be present among analysed data. For these reasons, nonparametric ANOVA tests have been preferred.

1s

150 100 50 0

0

500

1000

1500

2000

0

500

1000

1500

2000

0

500

1000 Time (ms)

1500

2000

Fm(t) (%MVC)

150 100 50 0

EMG(t) (AU)

1000 0

3.

Results

3.1 Model sensitivity and fitness analysis Following the methodology described Section 2.2, the first obtained Monte Carlo simulation series (Series 1) converged after 61 simulations. Eighteen simulations were selected based on the two selection criterions proposed by Keenan and Valero-Cuevas, and the whole 18 simulations also met the new proposed criterion, the RMS error between the input force and the output force (success rate: 18/18 ¼ 100%) (Figure 5 left). We then performed a second simulation series (Series 2) by using as initial intervals the ones refined from Series 1, with the muscular parameter ranges optimised by the three criteria. The convergence was achieved after 169 simulations. A total of 138 simulations were selected by using only the new ErrorRMS criterion on the force. Of these, 133 simulations satisfy the three criteria simultaneously (Figure 5 right). The success rate of the ErrorRMS alone is 96.4%. The 12 parameter ranges selected from the Series 2 are given in Table 1. We observed that the main parameters that interact on the three studied criteria are the neural parameters and especially the PFRs and the recruitment range (RR). In fact, the optimised interval for PFRs of fast MUs is higher

Figure 4. Constant force command (upper: force command P(t); middle: simulated force Fm(t); bottom: simulated sEMG(t)).

(Christensen 1959; Buchtal 1961) with 40% of slow fibres (type I) and 60% of fast fibres (type II) (Dahmane et al. 2005) and the average number of fibres per MU is considered equal to 72 (Gath and Stalberg 1981). To take into account the normal variability of healthy subjects, 10 muscular geometries are generated with parameters randomly chosen within a range of 5% of the average value of healthy subjects. All these parameters were chosen within the range that guaranteed the robustness of our model for isotonic contractions. In total, 120 simulations (4 contraction levels £ 3 trials £ 10 muscles) were made with a sampling frequency of 1000 Hz. The signal processing of the simulated sEMG signals and output normalised muscle force were similar to the experimental data.

Data analysis To qualify and quantify the measured sEMG signals, we calculated the root mean square (RMS) amplitude value 2

2

1.8

1.8

1.6

1.6

1.4

1.4

EMG/force (slope)

EMG/force (slope)

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

-1000

1.2 1 0.8 0.6

1.2 1 0.8 0.6

0.4

0.4

0.2

0.2

0

0

0.5

1 Force/force variability (slope)

1.5

2

0

0

0.5

1

1.5

2

Force/force variability (slope)

Figure 5. Selection of the simulations by the two experimental criteria (EMG/force and force/force variability represented by a blue frame) and by the RMS error criterion (red points) for the two series (left: Series 1; right: Series 2).

Computer Methods in Biomechanics and Biomedical Engineering

7

Table 1. Initial and optimised input parameter ranges for the sEMG – force model. Parameters RR (%MVC) Slow fibre CV (m/s) Fast fibre CV (m/s) SMUs number FMUs number Fibres number per MU VC of IPI Fibre length (cm) MFR of SMU (Hz) MFR of FMU (Hz) PFR of SMU (Hz) PFR of FMU (Hz)

Initial

Optimised (Series 2)

30–90 3–4 3–4 250–600 250–600 30–100 0.1–0.3 4–16 7–23 7–23 14–50 14–50

71– 89 3.0– 3.8 3.0– 4.0 250– 590 250– 570 34– 100 0.1– 0.3 6– 16 7– 20 8– 23 14– 41 27– 48

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

Note: VC, variation coefficient.

than the one obtained for the slow MUs. In addition, these criteria seem to be also sensitive to CV of slow MUs as depicted in Table 1. 3.2 Comparison with experimental data The experimental data obtained by the protocol described in Section 2.3 were compared with the simulated data obtained by means of the proposed model. According the previous section, 10 muscle configurations with the characteristics of a biceps brachii (Table 1) and the optimal parameter ranges obtained from the sensitivity analysis were used (previous study, Section 3.1). The force used as input command of the model (2 s duration) presented three phases: an increase of 0.3 s duration, a constant force maintained during 1.4 s and a decrease of 0.3 s duration. Four constant force levels (20%, 50%, 80% and 100% MVC) were simulated. We performed three simulations for each simulation condition (3 simulations of 2 s duration £ 4 force levels £ 10 muscles). In total, we generated 120 simulations at a sampling frequency of 1000 Hz. The simulated sEMG and force signals selected during a 1-s window in the middle of the whole recording (between 0.5 and 1.5 s) were taken into account in our analysis. We processed these simulated sEMG and force signals in the same way as for the measured signals (Figure 6). The RMS error between the simulated force and the command, in all the simulations, was less than 3%, which fits with the new validity criterion that we have proposed for the optimisation of our model. In the experimentation, the force level had a significant influence on the two parameters extracted from sEMG data (normalised RMS, MDF) (Friedman test, p , 0.05). During simulation, we noticed significant influence of the force level on the RMS parameter (Friedman test, p , 0.05) (Figure 7). However, no significant influence is observed for MDF parameter in both simulation and experiment (Figure 8).

Figure 6. Recorded sEMG signal, output force during isometric contraction at 50% MVC and the processing window (red, 1 s duration).

4. Discussions and conclusion In this paper, we propose a new MU rate coding scheme and analytical mechanogram for simulation of both sEMG and force data from the biceps brachii. Using the proposed sEMG/Force model, the surface EMG signal and the force generated by a muscle can be simulated according to a target force input supposed to be requested by the CNS. The new rate coding scheme supposes different firing/ excitation slopes and firing extreme values (MFR and PFR) between two MU types (slow and rapid). An analytical formalism for force generation has been presented. The proposed analytical mechanogram allows the tuning of MU force contributions by changing the slope of the force/firing rate relationship which is supposed to be linear. First, the proposed rate coding strategy supposed varying firing rate/excitation slope (Cheng et al. 2000) and varying MFR, PFR according to

Figure 7. Normalised sEMG/force relationship for simulated (SIM) and experimental (EXP) data.

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

8

H. Cao et al.

Figure 8. Evolution of the sEMG MDF according to the output force for simulated (SIM) and experimental (EXP) data.

the MU type (I: slow and II: fast). This strategy can be considered as a compromise in comparison with previous works with fixed slopes and firing rate intervals (Fuglevand et al. 1993) and random firing rate intervals like in Keenan and Valero-Cuevas (2007). Second, a linear behaviour is supposed concerning the produced MU force according to its firing rate and can be considered as fitting with the sigmoid relationship that exists between MU firing rate and force output (Fuglevand et al. 1999; McNulty et al. 2000). Regardless of the final tetanic phase, it can be considered that MU force increases quasi-linearly with increasing its firing rate as depicted in the study of Thomas et al. (1991) on human thenar muscle. Consequently, to simplify this relationship in the proposed model, it has been assumed as linear with a positive slope. The main advantage is also that the obtained force profile is deterministic and can be analytically obtained. A random force component has been added to the proposed force model to take into account force variability during sustained contractions. To evaluate the robustness/sensitivity of the proposed model against the variability of an input model parameter set, we have used a Monte Carlo method. For this purpose, we also proposed new simulation convergence and criterion. This new simulation criterion consists in measuring the error between the normalised input (target) force linked to neural command and the normalised output force delivered by the proposed model. The new convergence criterion named ErrorRMS is the limitation of this error to a specific threshold of 3%. The obtained results demonstrate the accuracy of the new criterion in comparison with the existing ones. Practically, it seems that the fitting accuracy of the realised force to the requested force by the CNS is highly correlated to the sEMG/force slope and force/force variability ratio. In fact, this new criterion speeded up the optimisation procedure

of the input parameters, since the simulations satisfying the new criterion met almost all the three criteria in the simulation Series 2 in 95% of the cases. This new criterion expressed the efficiency of the simulated neural drive and the modelled muscle to produce a specific force profile requested by the CNS. In comparison to two criteria based on sEMG – force slope and force variability (Keenan and Valero-Cuevas 2007) previously used, the new criterion has an equivalent performance with less computation time. As previously described, we also noticed that our model seems to be sensitive to rate coding and spatial recruitment parameters (Keenan and Valero-Cuevas 2007). In order to experimentally evaluate our model, muscle force and sEMG from the model and external force and sEMG from data recorded during isometric isotonic contractions of the biceps for four levels of contraction were compared. For the sEMG amplitude comparison, we noticed that the absolute or normalised RMS of the recorded sEMG signal increased almost linearly with the increase in force, which is consistent with simulated sEMG results (Figure 7). In the case of biceps isometric efforts, previous studies demonstrated that the relations between the amplitude of the sEMG amplitude (RMS value) and the associated force are either linear (Beck et al. 2005) or curvilinear (Lawrence and De Luca 1983). This contradiction may be due to the methodology of sEMG recording (electrode type and shape, electrodes placement, unipolar or bipolar configurations, etc.) (Moritani and deVries 1978) and differences of joint angles in which isometric muscle efforts were performed (Solomonow et al. 1991). In our experiment, for the MDF, the force level had a significant influence, while it had no significant influence in the simulation (Figure 8). The previous experimental studies on the effect of force on the MDF have also been inconclusive. When the force increases, the MDF can increase (Doheny et al. 2008) as in our experimental results (before 80% MVC), and decrease (Kaplanis et al. 2009), or vary non-significantly (Ravier et al. 2005; Talebinejad et al. 2009), as in our simulated results. Moreover, the value of MDF ranges between 60 and 80 Hz for the experiment, as in the studies of Talebinejad et al. (2009) and Doheny et al. (2008), and between 80 and 100 Hz for the simulation, as in the studies of Krivickas et al. (1998) and Kaplanis et al. (2009). These differences could be multiple and caused by (listed in a no exhaustive way): the limited number of subjects examined in the experiments, the difference in the duration of the experimentation, since fatigue can affect the results with longer durations (Lariviere et al. 2001), the methods used for statistical analysis, the used force level, the position of electrodes, the changes in experimental conditions (Shankar et al. 1989), the orientation of muscle fibres compared to surface electrodes, fibre length, thickness of individual subcutaneous layers (Farina et al. 2002),

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

Computer Methods in Biomechanics and Biomedical Engineering the number of active MUs for each level and each studied muscle and the anatomical simplifications. The proposed model also presents some limitations. First, personalised geometrical and morphological variations of the volume conductor are not taken into account during contraction. However, such parameterisation requires complex procedures and data processing (Marin 2013). In addition to that, the proportion of PFRs and MFRs according to the principal MU types (slow and rapid) is still a conjecture. Experimentally, PFR and MFR can be determined locally only invasively (Farina et al. 2002). Some promising works based on sEMG decomposition techniques pointed out possible varying MFRs and PFRs (De Luca and Hostage 2010) but could not be validated. Second, a rate coding scheme is presented and tested in simulation. Evaluating experimentally the exact rate coding strategy for a specific muscle also is a complex task and no consensus could be found in the literature (Fluglevand et al. 1993; Keenan and Valero-Cuevas 2007; De Luca and Hostage 2010). Third, the electrode positioning is an important parameter relative to the muscle NMJ. The modelling procedure placed the electrode at a theoretical position which cannot be exactly located in the experimental setup. However, as large interindividual differences are demonstrated (Barbero et al. 2012), this theoretical position could be considered as a compromise. To conclude, according to the obtained results and despite the discussed limitations, the proposed model seems precise enough to provide a realistic and promising behaviour of the sEMG –force relationship for the biceps brachii muscle. Preliminary investigation with isometric anisotonic contractions showed possible extension of the present model to various muscle contraction types (Cao et al. 2009). Further work will focus on a better understanding of the MUs recruitment and de-recruitment strategies at varying force level and pattern and possible MU activation turnover, during sustained contractions, which could be easily implemented in the present modelling to increase its physiological relevance. Funding This research was sponsored by the BDI grant supported by the Centre National de Recherche Scientifique (CNRS) and the Re´gion Picardie. Conflict of interest: There is no conflict of interest.

References An KN, Kaufman KR, Chao EY. 1989. Physiological considerations of muscle force through the elbow joint. J Biomech. 22:1249– 1256. Andreassen S, Rosenfalck A. 1978. Recording from a single motor unit during strong effort. IEEE Trans Biomed Eng. 25:501 – 508.

9

Basmajian JV, De Luca CJ. 1985. Muscles alive: their functions revealed by electromyography. Baltimore: Williams and Wilkins Ed. Barbero M, Merletti R, Rainoldi A. 2012. Atlas of Muscle Innervation Zones: understanding surface electromyography and its applications. Berlin: Springer-Verlag. ISBN: 978-88470-2462-5. Beck TW, Housh TJ, Johnson GO, Weir JP, Cramer JT, Coburn JW, Malek MH. 2005. The effects of interelectrode distance on electromyographic amplitude and mean power frequency during isokinetic and isometric muscle actions of the biceps brachii. J Electromyogr Kinesiol. 15:482 – 495. Buchanan TS, Lloyd DG, Manal K, Besier TF. 2004. Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. J Appl Biomech. 20:367 – 395. Buchtal F. 1961. The general concept of the motor unit. Res Publ Ass Nerv Ment Dis. 38:1– 30. Cao H, Boudaoud S, Marin F, Marque C. 2009. Optimization of input parameters of an emg-force model in constant and sinusoidal force contractions. Paper presented at: Proceeding of the 31th IEEE EMBS Conference; Minneapolis, MN, USA. Cao H, Marin F, Boudaoud S, Marque C. 2008. Muscle force simulation by using two motor-unit recruitment strategies. Comput Methods Biomech Biomed Eng. 11(Suppl. 1): 51 – 52. Cheng EJ, Brown IE, Loeb GE. 2000. Virtual muscle: a computational approach to understanding the effects of muscle properties on motor control. J Neurosci Methods. 101:117 – 130. Christensen E. 1959. Topography of terminal motor innervation in striated muscles from stillborn infants. Am J Phys Med. 38:65 –78. Christie A, Kamen G. 2010. Short-term training adaptations in maximal motor unit firing rates and after hyperpolarization duration. Muscle Nerve. 41(5):651– 660. Dahmane R, Djordjevic S, Simunic B, Valencic V. 2005. Spatial fiber type distribution in normal human muscle: histochemical and tensiomyographical evaluation. J Biomech. 38:2451 – 2459. De Luca CJ, Foley PJ, Erim Z. 1996. Motor unit control properties in constant-force isometric contractions. J Neurophysiol. 76:1503 –1516. De Luca CJ, Hostage EC. 2010. Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J Neurophysiol. 104:1034 – 1046. Disselhorst-Klug C, Schmitz-Rode T, Rau G. 2009. Surface electromyography and muscle force: limits in sEMG – force relationship and new approaches for applications. Clin Biomech. 24:225 – 235. Doheny EP, Lowery MM, Fitzpatrick DP, O’Malley MJ. 2008. Effect of elbow joint angle on force – EMG relationships in human elbow flexor and extensor muscles. J Electromyogr Kinesiol. 18:760 –770. Erim Z, De Luca CJ, Mineo K, Aoki T. 1996. Rank-ordered regulation of motor units. Muscle Nerve. 19:563– 573. Farina D, Fosci M, Merletti R. 2002. Motor unit recruitment strategies investigated by surface EMG variables. J Appl Physiol. 92:235 – 247. Farina D, Merletti R. 2001. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Trans Biomed Eng. 48:637 – 646. Farina D, Mesin L, Martina S, Merletti R. 2004. A surface EMG generation model with multilayer cylindrical description of

Downloaded by [Cukurova Universitesi] at 06:02 05 November 2014

10

H. Cao et al.

the volume conductor. IEEE Trans Biomed Eng. 51:415 – 426. Fuglevand AJ, Macefield VG, Bigland-Ritchie B. 1999. Forcefrequency and fatigue properties of motor units in muscles that control digits of the human hand. J Neurophysiol. 81:1718 – 1729. Fuglevand AJ, Winter DA, Patla AE. 1993. Models of recruitment and rate coding organization in motor-unit pools. J Neurophysiol. 70:2470– 2488. Gardiner PF. 1993. Physiological properties of motoneurons innervating different muscle unit types in rat gastrocnemius. J Neurophysiol. 69:1160 – 1170. Gath I, Stalberg E. 1981. In situ measurement of the innervation ratio of motor units in human muscles. Exp Brain Res. 43:377 – 382. Guyton C, Hall JE. 2000. Textbook of medical physiology. 10th ed. Philadelphia, PA: W.B. Saunders. Gydikov A, Gantchev N. 1989. Velocity of spreading of the excitation along the muscle fibres of human motor units measured by superficial electrodes. Electromyogr Clin Neurophysiol. 29:131– 138. Hagg GM, Luttmann A, Jager M. 2000. Methodologies for evaluating electromyographic field data in ergonomics. J Electromyogr Kinesiol. 10:301 –312. Henneman E. 1980. Skeletal muscle: the servant of the nervous system. In: Mountcastle VB, editor. Medical physiology. St Louis, MO: Mosby; p. 674– 702. Hof AL. 1997. The relationship between electromyogram and muscle force. Sportverletz Sportschaden. 11:79– 86. Jennekens FGI, Tomlinson BE, Walton JN. 1971. Data on the distribution of fibre types in five human limb muscles: an autopsy study. J Neurol Sci. 14:245 –257. Kaplanis PA, Pattichis CS, Hadjileontiadis LJ, Roberts VC. 2009. Surface EMG analysis on normal subjects based on isometric voluntary contraction. J Electromyogr Kinesiol. 19:157 – 171. Keenan KG, Valero-Cuevas FJ. 2007. Experimentally valid predictions of muscle force and EMG in models of motorunit function are most sensitive to neural properties. J Neurophysiol. 98:1581 – 1590. Krivickas LS, Taylor A, Maniar RM, Mascha E, Reisman SS. 1998. Is spectral analysis of the surface electromyographic signal a clinically useful tool for evaluation of skeletal muscle fatigue? J Clin Neurophysiol. 15:138– 145. Kukulka CG, Clamann HP. 1981. Comparison of the recruitment and discharge properties of motor units in human brachial biceps and adductor pollicis during isometric contractions. Brain Res. 219:45 – 55. Lariviere C, Arsenault AB, Gravel D, Gagnon D, Loisel P. 2001. Median frequency of the electromyographic signal: effect of time-window location on brief step contractions. J Electromyogr Kinesiol. 11:65 –71. Lawrence JH, De Luca CJ. 1983. Myoelectric signal versus force relationship in different human muscles. J Appl Physiol. 54:1653 – 1659. Lloyd DG, Besier TF. 2003. An EMG-driven musculoskeletal model to estimate muscle forces and elbow joint moments in vivo. J Biomech. 36:765 – 776. Maclntosh BR, Gardiner PF, McComas AJ. 2006. Skeletal muscle: form and function. 2nd ed. Champaign, IL: Human Kinetics Publishers Inc. Marin F. 2013. Influence of the mechanical definition of the functional spinal unit on the musculoskeletal modeling of the cervical spine. Movement Sport Sci, in press.

McNulty PA, Falland KJ, Macefield VG. 2000. Comparison of contractile properties of single motor units in human intrinsic and extrinsic finger muscles. J Physiol. 526:445– 456. Monster AW, Chan H. 1977. Isometric force production by motor units of extensor digitorum communis muscle in man. J Neurophysiol. 40:1432 –1443. Moritani T, deVries HA. 1978. Reexamination of the relationship between the surface integrated electromyogram (IEMG) and force of isometric contraction. Am J Phys Med. 57:263– 277. Moritz CT, Barry BK, Pascoe MA, Enoka RM. 2005. Discharge rate variability influences the variation in force fluctuations across the working range of a hand muscle. J Neurophysiol. 93:2449 – 2459. Ravier P, Buttelli O, Jennane R, Couratier P. 2005. An EMG fractal indicator having different sensitivities to changes in force and muscle fatigue during voluntary static muscle contractions. J Electromyogr Kinesiol. 15:210– 221. Santos VJ, Valero-Cuevas FJ. 2006. Reported anatomical variability naturally leads to multimodal distributions of Denavit – Hartenberg parameters for the human thumb. IEEE Trans Biomed Eng. 53:155 – 163. Shankar S, Gander RE, Brandell BR. 1989. Changes in the myoelectric signal (MES) power spectra during dynamic contractions. Electroencephalogr Clin Neurophysiol. 73:142 – 150. Solomonow M, Baratta RV, D’Ambrosia R. 1991. EMG – force relations of a single skeletal muscle acting across a joint: dependence on joint angle. J Electromyogr Kinesiol. 1:58 –67. Staudenmann D, Roeleveld K, Stegeman DF, van Dieen JH. 2010. Methodological aspects of SEMG recordings for force estimation – a tutorial and review. J Electromyogr Kinesiol. 20:375 – 387. Talebinejad M, Chan AD, Miri A, Dansereau RM. 2009. Fractal analysis of surface electromyography signals: a novel power spectrum-based method. J Electromyogr Kinesiol. 19:840 – 850. Thomas CK, Bigland-Richie B, Johansson RS. 1991. Force– frequency relationships of human thenar motor units. J Neurophysiol. 65:1509 –1516. Welch PD. 1967. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans Audio Electroacoust. 15:70– 73. Welter TG, Bobbert MF. 2001. During slow wrist movements, distance covered affects EMG at a given external force. Motor Control. 5:50 –60. Wheeler KA, Kumar DK, Shimada H. 2010. An Accurate bicep muscle model with sEMG and muscle force outputs. J Med Biol Eng. 30(6):393– 398. Woods JJ, Bigland-Ritchie B. 1983. Linear and non-linear surface EMG/force relationships in human muscles. An anatomical/functional argument for the existence of both. Am J Phys Med. 62:287– 299. Zhou P, Rymer WZ. 2004. Factors governing the form of the relation between muscle force and the EMG: a simulation study. J Neurophysiol. 92:2878 – 2886. Zhou P, Suresh NL, Rymer WZ. 2007. Model based sensitivity analysis of EMG –force relation with respect to motor unit properties: applications to muscle paresis in stroke. Ann Biomed Eng. 35:1521 – 1531.