Feasibility of using EMG driven neuromusculoskeletal model for prediction

Feasibility of using EMG driven neuromusculoskeletal modelfor prediction of dynamic movement of the elbowTerry K.K.Koo,Arthur F.T.Mak*Jockey Club Rehabilitation Engineering Centre,The Hong Kong Polytechnic University,Hung Hom,Kowloon,Hong KongAbstractNeuromusculoskeletal (NMS)modeling is a valuable tool in orthopaedic biomechanics and motor control research.To evaluate the feasibility of using electromyographic (EMG)signals with NMS modeling to estimate individual muscle force during dynamic movement,an EMG driven NMS model of the elbow was developed.The model incorporates dynamical equation of motion of the forearm,musculoskeletal geometry and musculotendon modeling of four prime elbow flexors and three prime elbow extensors.It was first calibrated to two normal subjects by determining the subject-specific musculotendon parameters using computational opti-mization to minimize the root mean square difference between the predicted and measured maximum isometric flexion and extension torque at nine elbow positions (0–120°of flexion with an increment of 15°).Once calibrated,the model was used to predict the elbow joint trajectories for three flexion/extension tasks by processing the EMG signals picked up by both surface and fine electrodes using two different EMG-to-activation processing schemes reported in the literature without involving any trajectory fitting procedures.It appeared that both schemes interpreted the EMG somewhat consistently but their prediction accuracy varied among testing pro-tocols.In general,the model succeeded in predicting the elbow flexion trajectory in the moderate loading condition but over-drove the flexion trajectory under unloaded condition.The predicted trajectories of the elbow extension were noted to be continuous but the general shape did not fit very well with the measured one.Estimation of muscle activation based on EMG was believed to be the major source of uncertainty within the EMG driven model.It was especially so apparently when fine wire EMG signal is involved primarily.In spite of such limitation,we demonstrated the potential of using EMG driven neuromusculoskeletal modeling for non-invasive prediction of individual muscle forces during dynamic movement under certain conditions.Ó2004Elsevier Ltd.All rights reserved.Keywords:Electromyography;Neuromusculoskeletal modeling;Muscle mechanics;Elbow1.IntroductionIn vivo estimation of individual muscle forces is important to understanding the mechanisms of injuries,and predicting the outcome of surgical interventions.Unfortunately,direct measurement of muscle force in vivo is beyond the current state-of-the-art.Therefore,indirect methods based on mathematical models are widely used.Four major approaches including static optimization [3],dynamic optimization [50],neural net-work [54],and EMG assisted method [18,38]have been used to estimate individual muscle forces.For the static optimization,an objective function such as minimizing the muscle stress or muscle activation is employed at each instant of a motor task to partition the net joint moments calculated by inverse dynamics into respective muscular components.The drawbacks of this method are that the results are very sensitive to the choice of the objective function and it is very difficult to determine an appropriate objective function,which results in com-parable temporal patterns between the estimated muscle forces and measured EMG.Consequently,the muscle forces predicted by this method are questionable.1050-6411/$-see front matter Ó2004Elsevier Ltd.All rights reserved.doi:10.1016/j.jelekin.2004.06.007*Corresponding author.Tel.:+852********;fax:+852********.E-mail address:rcafmak@.hk (A.F.T.Mak)./locate/jelekinJournal of Electromyography and Kinesiology 15(2005)12–26Dynamic optimization is a process of choosing neural excitations by optimizing a well-defined performance criterion such as maximizing the jumping height in a ver-tical jumping task[50]and combining them with a neu-romusculoskeletal model to predict muscle forces and resulting kinematics.The temporal EMG patterns and the movement kinematics are usually used to validate the model indirectly.Theoretically,this approach pro-vides the preferred avenue for understanding human motor control strategy and muscle excitation signals. However,it is only applicable to normal subjects prima-rily because the so-called‘‘optimal performance’’may not be readily achieved by pathological subjects and is limited to movements in which the performance crite-rion is well defined.Obviously,the lack of versatility over a wide range of movements and its demanding computational cost limited its potential of being a useful tool for muscle force estimation.The use of artificial neural network to predict individual muscle forces is a relatively new approach.Although the method has potential advantages in terms of data collection,adapt-ability,and prediction;and the outputs can be validated indirectly by comparing the predicted net joint moments or kinematics with the measured one,its reliability is yet to be validated as the neural network algorithm acts in essence as a‘‘black box’’and is often not explicitly transparent.Using EMG signals as inputs to a musculoskeletal model to estimate individual muscle forces have been used by several groups for either isometric or dynamic tasks at different anatomical locations such as elbow [7,16,18,31,55],shoulder[35],knee[38–40,49,56],ankle [25–29,49],jaw[52],lower back[12,22,42,43,48],and wrist[9].The theoretical basis behind EMG driven model is that if the EMG signals can be measured precisely and processed adequately to reflect the acti-vation of each muscle crossing the joint and if the activation can be modulated properly by anatomical and musculotendon models,it is possible to accurately estimate individual muscle forces over a wide range of tasks and contraction modes.This makes EMG driven model to be a very attractive choice if individual mus-cle force estimation,rather than others like under-standing the motor control strategy,is the primary concern of the researchers.Implicit to the EMG dri-ven model is the EMG-to-activation processing,the musculotendon model,and the geometrical model of the joint,which entail a set of parameters that could affect the behavior of the model.Among many param-eters,some of them are regarded as calibration param-eters,which required to be tuned by some kinds of calibration procedure and the other parameters are re-garded asfixed parameters,whose values are either extracted from the literature[42]or assigned based on some subjective criteria[55].Depending on the model structure and assumptions made,the calibra-tion parameters can range from parameters with no physiological meaning such as gain or error factor [29],and EMG coefficients[9]to parameters with physiological meaning such as muscle rest length [11,33],tendon slack length[33,38],maximum isomet-ric stress[22,33],shape factor of the EMG–force rela-tionship[38,48],and electromechanical delay[18,48]. Different methods have also been used to determine the calibration parameters.For instance,some researchers use many calibration trials that manipulate certain conditions while attempting to hold others constants so that contribution of individual muscle can be isolated[27,56].Others attempt to determine multiple calibration parameters simultaneously using a set of calibration trials by systematically vary the parameters using non-liner least square optimization until the difference between the measured and pre-dicted net joint moment is minimized[7,11,33,40].An-other method is to use electrical stimulation to selectively activate a particular muscle and determine the parameters specific to that muscle[33,34].One of the limitations of most of the EMG driven models reported in the literature is that they tend to as-sume that once the model was calibrated by a limited set of calibration trials,the same set of model parameters could still yield reasonable prediction in other tasks and testing configurations.That is especially so for the previous elbow joint models[7,16,18,55].Before an EMG driven model can be used as a reliable tool for the estimation of muscle forces,additional validation across other tasks and testing configurations must be done.This study sought to develop an EMG driven neuro-musculoskeletal(NMS)model of the elbow that predicts individual muscle forces using surface andfine wire EMG recorded during voluntary elbowflexion–exten-sion movements.A set of subject-specific musculoten-don parameters within the model wasfirst determined based on a series of maximum isometricflexion and extension trials at nine elbow positions(0–120°offlex-ion with an increment of15°).Once calibrated,we sought to use the same set of parameters and the EMG signals as inputs to test the validity and versatility of the model at different dynamic movement activities by comparing the accuracy and consistency of the joint tra-jectories predicted by the model.2.Methods2.1.SubjectsTwo normal subjects(subject A:male,25years,174 cm,64.5kg;subject B:female,30years,156cm,44.5 kg)were recruited.Subjects gave their informed consent to the investigation,which was approved by the HongT.K.K.Koo,A.F.T.Mak/Journal of Electromyography and Kinesiology15(2005)12–2613Kong Polytechnic University human subject ethics committee.2.2.Experimental protocolThe subject sat on an assessment chair with the upper arm vertical and supported by an adjustable wooden plate.The shoulder was orientated at0°offlexion,10°of abduction,and neutral rotation.The forearm was in a supinated position and the wrist was in a neutral position with the hand completely relaxed.Aflexible electrogoniometer(Penny&Giles,UK),which was at-tached on the lateral side of the arm when the forearm was in a fully supinated position,was used to measure the angular displacement of the elbow joint.The elbow angle was defined as0°when the elbow was fully ex-tended.An eight channel telemetric surface EMG sur-face system with a bandwidth of10–500Hz and a gain of2000(Noraxon U.S.A.Inc.,Scottsdale,AZ)was used to capture the surface EMG signals of the long and short head of biceps(LHB and SHB),the brachioradi-alis(BRD),and the medial,lateral,and long heads of triceps(MHT,LatHT,and LngHT)picked up by pre-gelled Ag/AgCl6mm diameter bipolar surface elec-trodes(Blue Sensor,Medicotest,Denmark).Only a pair of electrodes was attached over the muscle belly of the biceps because the LHB was assumed to have the same activation as the SHB.The electrodes were placed paral-lel to the presumed direction of the musclefibers and at the locations suggested by Cram et al.[13]with inter-electrode distance set to2cm.A pre-sterilized Teflon coated nickel alloy pairedfine wire electrode(Nicolet Biomedical,Madison,WI),which was ready-loaded in a27Ga.30mm long hypodermic needle was used to measure the EMG signals of the brachialis(BRA). Diameter of each wire was50l m with2mm of metal exposed at the tips of the hooked ends.The hook ends were positioned so that2mm of one wire and5mm of the second wire exited at the end of the hypodermic needle,which ensured that there was no overlap of the bared tips.With the forearm positioned to the extended and supinated position,thefine wire electrode was in-serted about twofingerbreadths proximal to elbow crease along and just lateral to the biceps tendon to a depth about20mm.Correct placement of the electrode was confirmed by asking the subjects toflex the forearm in pronation against minimal resistance and see whether the EMG signals augmented significantly[51].After insertion,thefine wire electrode was connected to a dif-ferential amplifier(B&L Engineering,Tustin,CA)with a CMRR of95dB,an input impedance>100M X,a bandwidth of10–1000Hz,and a gain of330via a pair of screw-springs adaptors.A2cm loop was formed be-tween the insertion site and the pre-amplifier to avoid pulling on thefine wire electrode during movement.In addition,a pre-gelled Ag/AgCl10mm diameter surface ground electrode(Blue Sensor,Medicotest,Denmark) was placed over the olceranon.The skin around the elec-trode attachment or insertion sites were shaved and rubbed lightly by a cotton wool thatfilled with abrasive skin prepping gel(Nuprep e,USA)so that skin imped-ance could be greatly reduced.A cotton wool pad that soaked with alcohol was then used to clean the skin.Prior to the test,each subject was instructed to be completely relaxed with his/her forearm in a fully ex-tended position(i.e.0°)for theflexion test or in a fully flexed position($130°)for the extension test.The sub-jects were then instructed toflex/extend the elbow at his/her preferred speed(mean speed:151.6°±31.1°). Three experimental tasks were conducted and each pro-tocol was repeated three times.Thefirst task involved voluntary elbowflexion from a fully extended position to a fullyflexed position with a1.82-kg dead weight strapped around the arm at the location just proximal to the wrist.The second and the third tasks involved vol-untary elbowflexion and extension,respectively,with-out any added weight.The subjects had several practice trials prior to data collection.During the tests, both the EMG signals and elbow joint angles were re-corded simultaneously at2000Hz and stored in a PC computer via a12bit A–D converter(Data Translation, USA)for further off-line processing using Matlab(The Mathworks,Inc.,Natick,MA).A series of maximum voluntary contraction(MVC) trials was done right after completion of the experimen-tal tasks.The goals were to determine the subject specific musculotendon parameters and collect the MVC EMG amplitude of each muscle to normalize the EMG signals recorded during the experimental tasks.The MVC test-ing procedure has been reported elsewhere[33,34]. Briefly speaking,a Cybex dynamometer(Cybex Interna-tional Inc.,Medway,MA)was used to measure the tor-que during elbowflexion and extension while the subjects exerted maximum isometric voluntaryflexion (MIVF)and maximum isometric voluntary extension (MIVE)torques,respectively.Nine elbow positions (i.e.0–120°offlexion in increments of15°)were tested. At each position,the subject wasfirst instructed to relax and then gradually to increase theflexion or extension torque to maximum and to maintain at that level for approximately2s.Visual feedback of the torque was gi-ven to the subjects.The MIVF and MIVE were con-ducted alternatively,and a2-min recovery period was allowed between contractions to avoid muscle fatigue. Before the test,warm-up trials were allowed until the subjects became familiar with the testing protocol.2.3.Model developmentAn EMG driven NMS model of the elbow was devel-oped using Matlab(The Mathworks,Inc.,Natick,MA), which consisted of the linkage dynamics of the forearm14T.K.K.Koo,A.F.T.Mak/Journal of Electromyography and Kinesiology15(2005)12–26(including the hand)(Section2.5),the musculoskeletal geometry(Section2.7),the EMG-activation processing (Section2.4)and the musculotendon modeling(Section 2.6)of four prime elbowflexors(i.e.,LHB,SHB,BRA, BRD)and three prime elbow extensors(i.e.,MHT, LatHT,LngHT).Selection of these muscles was based on computing and comparing the potential moment contribution among the24muscles wrapping across the elbow joint[2].Inputs to the model included the EMG signals of each targeted muscle as well as the ini-tial elbow joint position and ing the subject-specific musculotendon parameters estimated from a series of maximum isometricflexion and extension trials (Section2.9),the model computed the elbow joint trajec-tories without trajectoryfitting procedure involved.The overall structure of the EMG driven model is illustrated in Fig.1.2.4.EMG-activation processingThefiltered EMG signals were converted to muscle activations using two different EMG-activation process-ing schemes reported in the literature–EMG to active state converter[26]and linear envelope processing[49].A brief description of each of them is as follows.2.4.1.EMG to active state converterBriefly speaking,the EMG to active state converter can be described mathematically asaðtÞ¼MaxMax½EðtÀtÃÞ 8tÃ:3:06tÃ<s2;Max½EðtÀtÃÞeÀðtÃÀs2Þs3 8tÃ:3:tÃP s2;(ð1Þwhere.’.stands for‘‘such that’’;a(t)is the active state, E(t)is processed EMG(i.e.bandpassfiltered,full-wave rectified,moving averaged with a50ms window and normalized by the corresponding MVC EMG).The ‘‘Max’’operator denotes that the maximal value of the respective functions is to be taken for the interval indi-cated.The working is such that a(t)follows E(t)when E(t)is rising.When a relative maximum of E(t)occurs, this maximum is held for a time period of s2and fol-lowed by an exponential decay with a time constant s3.This plateau or decay of a(t)goes on until E(t)or the plateau resulting from a subsequent relative maxi-mum E(t)becomes larger than a(t).The plateau dura-tion s2and decay constant s3were set to be30and60 ms,respectively,according to[28].This EMG-active state converter captures some well known properties of the muscle activation.For example, the response to a short EMG burst has a form of muscle twitch.Moreover,in quasi-static state,the smoothed rectified EMG is linearly proportional to the active state.Furthermore,the model shows summation effects in response to a series of impulses,to a step function, and to some combinations of steps and impulses.2.4.2.Linear envelope processingLinear envelope processing can be regarded as an EMG processing scheme that converts the raw EMG signals to a linear envelope(LE)profile that mimics the muscle tension waveform during dynamic changes of isometric tension.The LE profile is obtained by full-wave rectification of the bandpass EMG signals fol-lowed by forward and backward lowpassfiltering using a second-order Butterworthfilter with cut-offfrequency of3Hz andfinally normalized by the corresponding MVC EMG.The theoretical basis of LE processing is that the frequency characteristics of muscle twitches could be regarded as a second order,critically damped, low-passed system[44].Since each twitch can be consid-ered to be the impulse response of an active motor unit associated with a motor unit action potential(MUAP) impulse,the full-wave-rectified EMG can be considered as a summation of MUAPs of various amplitudes, whereas,muscle tension can be considered as a graded summation of the twitches resulting from all active mo-tor units.It has been reported that the contracting timeT.K.K.Koo,A.F.T.Mak/Journal of Electromyography and Kinesiology15(2005)12–2615of a twitch for different muscles ranged between20and 90ms[10,44].Since the cut-offfrequency of3Hz of a second-order Butterworthfilter corresponds approxi-mately to a contracting time of50ms,which can be re-garded as the averaged value reported in the literature,it represents a reasonable selection for the elbow muscles.2.5.Linkage dynamicsThe linkage dynamics of the EMG driven model de-scribed theflexion/extension of the elbow in the sagittal plane.The upper extremity was treated as consisting of two rigid bodies–the upper arm and the forearm,which included the hand.Theflexion/extension movement of the elbow joint was modeled as a frictionless hinge joint rotation of the ulna with respect to the humerus,with the axis of rotation passing through the centers of capi-tellum and trochlear sulcus[41].Carrying angle of the elbow joint was assumed to be unchanged as the fore-armflexed and extended about the elbow joint.For each subject,the segment mass,the moment of inertia,and the location of center of gravity of the fore-arm with respect to theflexion/extension axis were esti-mated using the formulas suggested by Winter[57]and Dempster[15].For thefirst protocol,the dead weight was treated as a hollow cylinder so that its moment of inertia about theflexion/extension axis could be approx-imated.Therefore,the second-order differential equa-tions that govern the motion of the forearm with and without the dead weight attached can be described, respectively,as follows:€q¼1I CX7j¼1F tjðtÞÃMA jðqðtÞÞÀm B gLÃBþm A gL AI CsinðqðtÞÞ;ð2aÞ€q¼1I BX7j¼1F tjðtÞÃMA jðqðtÞÞÀm B gLÃBI BsinðqðtÞÞ;ð2bÞwhere q;€q are the elbow angle and angular acceleration,respectively,F tj are tendon force transients for the j thmuscle,MA j are theflexion–extension moment arms ofF tj relative to the elbow joint center,g is the gravita-tional acceleration,m A and m B are the masses of the dead weight and the forearm,respectively.I B and I C are the moment of inertia of the forearm and the com-bined moment of inertia of the forearm and the dead weight about theflexion–extension axis,respectively.The parameter LÃB indicates the distance between thecenter of elbow joint and the mass center of the forearm and L A indicates the distance between the dead weight and the elbow joint center.2.6.Musculotendon modelingThe musculotendon model used in the current study was based on Giat et al.[20].Briefly speaking,each mus-cle was biomechanically modeled as a unit of three par-allel components(i.e.,CE,PE,and VE)where CE represents the active contractile element,and PE and VE represent the non-linear elastic component and the viscous component of thefluid-filled connective tissues of the muscles,respectively(Fig.2).It was assumed that all the musclefibers are parallel,equal in length and ori-ented at the same pennation angle(a).T represents the tendon,which was modeled as a stiffand non-linear elastic spring.M is the muscle mass and l m,l t,and l mt represent the musclefiber length,tendon length,and musculotendon length,respectively.The tendon length value represents the total of proximal and distal tendon lengths within a muscle.The empirical relationship be-tween muscle,tendon,and musculotendon lengths can be stated as follows:l mt¼l m cos aþl tð3Þand the equation of motion of each muscle along the tendon axis of pull can be given asMÁ€l m cos a¼F tÀF m cos a;ð4Þwhere F m=F p+F d+F a and F a¼r mÁPCSAÁf lðl mÞÁf vð_l mÞÁaðtÞ.Herein,F t is the tendon force,which was modeled to have an exponential force–strain relation-ship for strain<2%and a linear relationship afterwards[59].The muscle force(F m)equals the sum of the passive(F p),viscous(F v),and active(F a)force components of a muscle.F a is a function of the maximum muscle stress (r m),the physiological cross-section area(PCSA),the generic force-length function(f l),the generic force-ve-locity function(f v),and the level of muscle activation (a).F p was modeled to increase exponentially with re-spect to the musclefiber length[20]and F v was assumed to be proportional to the rate of change of musclefiber length[20].f l and f v were those employed by Giat et al.[20].Values of thefixed parameters such as PCSA and pennation angle within the model were either extracted or scaled from the literature[1,2,46,59],leaving three subject-specific musculotendon parameters per muscle16T.K.K.Koo,A.F.T.Mak/Journal of Electromyography and Kinesiology15(2005)12–26(i.e.,maximum muscle stress(r m),muscle optimal length(l mo),and tendon slack length(l to)),which were used to specify the force generating characteristics of each musculotendon unit.r m was assumed to be the same for muscles that belonged to the same func-tional group.Hence,the maximum muscle stress of the flexors muscles(i.e.,LHB,SHB,BRA,BRD)was de-noted as r f and that of the extensor muscles(i.e., MHT,LatHT,LngHT)was denoted as r e.All together, the EMG driven model featured16musculotendon parameters.2.7.Musculoskeletal geometryMoment arms(MA j)and musculotendon lengthsðl mtj Þof each targeted muscle over the elbow range ofmotion were estimated using a subject-specific interac-tive graphic-based geometrical model of the elbow joint, which was scaled from a generic geometrical model of the elbow joint[33,34]developed in SIMM(Musculo-Graphics,Chicago,USA)environment[14].The generic model defined the bone surfaces of the humerus,ulna, radius,hand,rib cage,scapula and clavicle,the joint kinematics of the elbow,and the musculotendon paths of those selected muscles.The bone surfaces were pro-vided by the SIMM,which were obtained by digitizing a male skeleton with known anthropometric dimen-sions.The musculotendon paths of those selected mus-cles were modeled as a series of points connected by line segments,which included the muscle origin and insertion points as well as the additional intermediate points that were defined when the muscle wrapped over the joint surface.Muscle origin and insertion points of each selected muscle were defined as the centroids of muscle attachment areas marked on a commercially available A11muscle skeleton(Kappa Medical,Pres-cott,USA).For each bone segment,three anatomical landmarks were identified to define a local coordinate system and the coordinates of those muscle origin and insertion points were digitized with respect to the afore-said local coordinate systems using a mechanical3-D digitizer(MicroScribe-3D,USA).These coordinates were then transformed and scaled to match the local coordinate systems of the generic model.The generic model was then scaled to specific subject model based on the ratio of segment lengths.2.8.Solution procedurePutting together,the dynamics of the NMS system can be represented as a system of second-order ordinary differential equations,which consists of the linkage dynamics(Eq.(2))and the musculotendon dynamics (Eq.(4))of each muscle.The system can also be ex-pressed as afirst-order vector–matrix differential equa-tion(Eq.(5))by defining a set of state variables:x1¼q;x2¼_q;x2jþ1¼l mj;x2jþ2¼_l mjfor j¼1;2; (7)_xðtÞ¼fðxðtÞ;aðtÀEMD emÞÞ;ð5Þwhere x=[x1,x2,...,x16]T is the state vector,a=[a1, a2,...,a7]T is the input vector,whose elements represent the muscle activation of each targeted muscle,and EM-D em is the electromechanical delay,which was set to50 ms for all muscles.The initial elbow position,i.e.x1(0), was set to0°forflexion tasks and120°for extension tasks.The initial elbow angular velocity,i.e.,x2(0), and the initial muscle velocity,i.e.,x2j+2(0)(for j=1,2,...,7)were all set to zero,and the initial muscle lengths,i.e.,x2j+1(0)(for j=1,2,...,7)were estimated from Eq.(3)by setting muscle velocity and acceleration as zero.A variable-order solver that based on the numerical differentiation formulas was used to solve the initial value problem owing to its high efficiency in solving stiffproblem[23].2.9.Determination of musculotendon parametersThe musculotendon parameters of theflexor(exten-sor)muscles were estimated by minimizing the root mean square difference between the measured and pre-dicted maximum isometricflexion(extension)torques using Nelder–Mead simplex method[47].A number of constraints were imposed on the optimization process. First,the total moment produced by the summation of passive elastic component of each muscle at10°offlex-ion was set to be less than2N m.This constraint is con-sistent with in vivo measurement of passive moment [21].Second,the tendon slack lengths must be greater than zero.Third,the muscle optimal lengths(l mo)fol-lowed the sequence:l mo(BRA)<l mo(LHB)<l mo (SHB)<l mo(BRD)for theflexors and l mo(MHT)<l mo (LatHT)<l mo(LngHT)for the extensors.This con-straint is consistent with cadaver data[1,2,46].Muscle optimal angles(/o)were either estimated using electrical stimulation[33]or adapted from the cadaver data[46]to act as constraints(i.e.,l mt(/o)=l to+l mo cos a)to guide the parameter searching process.Initial estimates of the length-related parameters of each muscle were based on anthropometric scaling of the cadaver data[46]and an initial estimate of maximum muscle stress was set to100N/cm2,which is within the range reported in the literature.2.10.Data analysisThe accuracy of the prediction was quantified by cal-culating the root mean square(RMS)difference between the predicted and measured joint trajectories.One-way ANOVA was used to test whether the prediction accu-racy was task dependent and paired t test was used toT.K.K.Koo,A.F.T.Mak/Journal of Electromyography and Kinesiology15(2005)12–2617。