The control of shoulder muscles during goal directed movements, an inverse dynamic analysis

Pergamon

J. Biomechonws,

Vol. 28, No. 10, pp. 1179-1191, 1995 Copyright 0 1995 Eketin Science Ltd Printed in Great Britain. Ail tights resmved 0021-9290/95 19.50 + .oO

0021-9290(94)00181-2

THE CONTROL DIRECTED

OF SHOULDER MUSCLES DURING GOAL MOVEMENTS, AN INVERSE DYNAMIC ANALYSIS R. Happee* and F. C. T. Van der Helm

Man-Machine-Systems Group, Delft University of Technology, Faculty of Mechanical Engineering & Marine Technology, Mekelweg 2, 2628-CD Delft, The Netherlands Abstract-Fast goal directed arm movements in thesag&al plane were analyzed with a three-dimensional shoulder model with 95 muscle elements. Dynamics of the muscle elements were described by a third-order nonlinear muscle model. Muscle forces and activation were estimated using the method of inverse muscular dynamics, an optimization scheme which uses only very limited computational power. Most model results were similar to the EMG but some differences between model results and EMG were found in muscles where the EMG activity was subject dependent. For themovement studied, the thoracoscapular muscles were shown to deliver about 40% of the energy required for the acceleration of the arm during anteflexion and about 22% during retroflexion. Activity of thoracoscapular muscles was also required to ensure contact between the thorax and the scapula which is important for the mechanical stability of the shoulder. The rotator cuff muscles were found to deliver about 19% of the energy required for the acceleration of the arm during anteflexion and about 8% during retrotlexion.

INTRODUCTION

The shoulderis a complex mechanicalstructure containing several joints connecting the humerus, the scapula,the clavicle and finally the sternum(Fig. 1). Furthermore, the scapulacontacts the dorsalpart of the thorax; it can glide over the so-calledscapulothoracic gliding plane. This connection makesthe shouldera closedchain mechanism.As a result, the clavicle and thorax constrain the scapularrotations, and in addition, the moments about the sternoclavicular and acromioclavicularjoint are coupled. The complexity of the shouldermakesit very difficult to understandmuscularcoordination of the shoulder, or even to comprehend the motions of shoulder bones. Until now, the shoulderhasbeenstudiedmainly in static conditions.

Three-dimensional

positions

of the

arm and the bonesof the shouldergirdle have been recordedfor a number of abduction and anteflexion anglesof the humerus(Pronk, 1991;Pronk and Van der Helm, 1991;Van der Helm and Pronk, 1994).The relation betweenthe rotations of humerusand scapula is commonly referred to as the ‘scapulohumeral rhythm’ (H6gfors ef al., 1991; Inman et al., 1944: Pronk, 1989,1991).This relation is reasonablyconsistent and reproducible,although the assessed valuesof the rotations of the scapulaand the humerusshow someindividual variations, and show a dependency on the measurementequipment (roentgen versus Received in final form 14 November 1994. *Current address: TN0 Crash-Safety Research Centre, Biomechanics Section, Schoemakerstraat 97,2628-VK Delft, The Netherlands.

goniometers)and the experimental paradigm (twodimensionalversusthree-dimensional). EMG studies in static conditionshave provided information about the coordinatedactivity of shouldermuscles(Flanders and Soechting,1990,Inman et al., 1944).The mechanical function of suchmuscularactivity hasbeenanalyzed with three-dimensionalmodels(Karlsson and Peterson,1992;Van der Helm, 1994a,b). The coordination of shouldermusclesin dynamic conditionshasbeenstudiedexperimentallyin Happee (1992b).Goal directedarm movementsin the sagittal plane, mainly consistingof an anteflexion or retroflexion in the shoulder (moreformally describedas anteflexion/retroflexion of the glenohumeraljoint) have been studied.The EMG of 13 musclesin the shoulder and the arm has been recorded.Triphasic activation patterns associatedwith accelerationand decelerationof the limb have beenfound in muscles which can be regarded as prime movers of the glenohumeraljoint (deltoideus,pectoralis and latissimusdorsi).The patternsin thoracoscapularmuscles werevery similarto patternsin theseprime moversof the glenohumeraljoint. This indicatesa function of thoracoscapularmusclesin the acceleration of the arm which isalsosupportedby a statisticalanalysisof EMG and kinematics (Happee, 1993a).Antagonist activity was rare in the EMG of thoracoscapular muscles. Thus it is not clearin which way decelerating forcesinducedby the primemoversare transferredto the trunk (Happee,1992b). Such EMG studies have considerable limitations. Recording the EMG of deeperlying musclesyields severetechnicalproblems,especiallyin dynamic conditions and for smallermuscles.EMG provides no information on the forces occurring in joints and

1179

1180

R. Happee and F. C. T. Van der Helm in how far the different muscles contribute mechanical energy for the acceleration/deceleration of the arm. Here the thoracoscapular muscles, which only indirectly affect the arm, are of particular interest. The second question studied is how the scapula and the glenohumeral joint are being stabilized. METHODS

A three-dimensional model of bones and ligaments To analyze displacements of the shoulder, a kinmatic model of the shoulder girdle was developed by Pronk (1989, 1991). To permit analysis of muscle forces, the model was transformed into a dynamic model (Veeger et al., 1991; Van der Helm and Veenbaas, 1991; Van der Helm et al., 1992; Van der Helm, 1994a, b). In the model the thorax, clavicle, scapula, humerus and forearm are represented as rigid bodies. Motion constraints of the scapulothoracic gliding plane are included. This model contains 20 muscles, Fig. 1. Frontal view of the upper part of the body. The bones divided into a total 95 muscle elements as shown in indicated are: the humerus(h), the clavicle(c), the scapula (SC) Table 1. Muscle lines of action curve around underlyand the sternum (s), The articulations indicated are: the ing bony contours. In addition to this existing model, sternoclavicular joint (SC), the acromioclavicular joint (AC), for the current paper, the effects of nonlinear muscular the glenohumeral joint (GH) and the scapulothoracic gliding dynamics are included by using the method of inverse plane (ST); reprinted with permission from Pronk (1989). muscular dynamics (Happee, 1994). The resulting shoulder model will now be described in more detail. Degrees offreedom (DOF): The sternoclavicular, acligaments. EMG can show when musclesare active and glenohumeral joints are debut cannot explain the mechanical function of this romioclavicular scribed as spherical joints each having three DOF. activity. The elbow is described as a joint with only one DOF For these reasons it was decided to analyze the goal which is elbow flexion. The conoid ligament condirected movements studied in Happee (1992b) using strainsthe motionsof the clavicle with respectto the a detailed musculoskeletal model. The goal of this scapula. As the conoid ligament is assumed stiff this simulation study was to provide insight into the function of shoulder muscles. The first question studied is yields one kinematic constraint. The scapulothoracic

Table 1. Number of elements, physiological cross-sectional area (PCSA) per element and location order of the elements for future reference in figures, for the muscles of the shoulder mechanism Muscle (part) Trapezius pars scapularis Trapezius pars clavicularis Levator scapulae Pectoralis minor Rhomboideus Serratus anterior Deltoideus pars scapularis Deltoideus pars clavicularis Coracobrachialis Infraspinatus Teres minor Teres major Supraspinatus Subscapularis Biceps brachii caput longum Biceps bra&ii caput breve Triceps bra&ii caput longum Latissimus dorsi Pectoralis major pars thoracaIis Pectoralis major pars davicularis

No. of elements

PCSA/element (cm*)

6 6

2.39 0.52

3 4 3 6 6 6 6 6 6 6 6 6

1 1 2 5 5 5

1.15 0.86 2.52 1.90 2.76 1.35 0.53 1.36 0.52 2.09 0.78 2.50 3.12 3.12 3.12 1.70 1.74 0.71

Element order Caudal-cranial Caudal-cranial Caudal-cranial Medial-lateral Caudal-cranial Caudal-cranial Medial-lateral Lateral-medial Lateral-medial Cranial-caudal Cranial-caudal Cranial-caudai Ventral-dorsal Cranial-caudal Medial-lateral Cranial-caudal

Caudal-cranial Medial-lateral

Shoulder model

gliding plane is modeledby two surfaceelementsat the angulusinferior (AI) and the trigonum spinae(TS) of the scapula,respectively,It is assumedthat during the active arm movementsstudied,contact at these two surface elementsis maintained. This contact yields two kinematic constraints.How this contact will be maintained will be treated in the following section.The sternoclavicular,acromioclavicular,glenohumeral and elbow joints yield 10DOF, the conoid ligament and the scapulothoracic gliding plane impose3 kinematic constraintswhich yieldsa model with 7 kinematicDOF. This specifies 7 equality constraints relating muscleforcesto joint torques. Furthermore,the model contains4 inequality constraint on the muscularforceswhich will now be treated. Linear inequality constraints: The conoid ligament can practically transfer only tension forces,which is taken into account by one inequality constraint (Van der Helm, 1994a).The forces at the scapulothoracic gliding plane have to be compressionforces which yields two inequality constraintsfor AI and TS, respectively. Glenohumeral constraint: The glenohumeraljoint hasbeenimplementedasa sphericaljoint with three rotational DOF, not allowing translation (e.g.dislocation). A constraint on muscleforces, preventing dislocationwill now be described.The passivestructures in and around the glenohumeraljoint do not guaranteejoint stability (Lippitt and Matsen, 1993).It is assumedthat the glenohumeraljoint is actively stabilized;the surroundingmuscleshaveto becontrolled so as to prevent dislocation. Stability of the glenohumeraljoint has beenincluded asa nonlinear constraint as describedin Van der Helm (1994a). Stability is obtained if the reactionforcespoint inside the glenoid. The rim of the glenoidjoint surfacehas been describedas an ellipse, Muscular forces have beenconstrainedsothat the net joint reactionforce is directed from the joint rotation center to any point within this ellipse.If the reaction force vector would

1181

setof elbow muscles.For the movementstudiedhere, for accelerationand decelerationof the limb, large torques are required in the shoulderwhereasin the elbow a small nearly constant torque is required (Happee, 1992b).For these reasonsthe movement equationprescribingthe elbow torque hasbeendisregarded.The kinematic degreeof freedomis left in the model.Thus 6 out of the 7 linear equality constraints representingthe movement equations will be considered.Three linear inequality constraintsrepresent both capulothoracic forcesand the conoid ligament force. The stability of the glenohumeraljoint yields one nonlinear inequality constraint. Muscle forces will be chosensuch that theseconstraintsare being met. This will ensurethat the assumptionsmadedefining the DOF arevalid. Muscleforceswill alwaysbe positive or zero, yet musculardynamicsimposefurther constraintson muscleforcesaswill be described in a following section. Muscle model

A nonlinear,lumped-parametermodelof muscular dynamicsis employed which has been presentedin Happee(1994)and which is largely basedon Winters and Stark (1985, 1988).The model contains a firstorder description of neural-excitation dynamicsand a first-order active-statedynamics. The third order emergesfrom the force-velocity relation (Hill curve) combinedwith a serieselasticelement(SE). An active force-length relation is not yet implementedasmajor problemsoccurred in the choiceof musclelengthsat which maximal active force is delivered. This is acceptableasfor the current paper only relatively small changesof musclelength are considered.For similar reasons the model does not contain a passive force-length relation describinga parallel elasticelement (PE) which isjustified becausefor this paperno extreme joint anglesare considered(Engin, 1980). Information on musclepropertieslike the slow/fast fiber composition was not available for all muscles point outside the glenoid, it could not be counteracted considered,thereforemuscleparameterswerechosen by the glenoid reaction force and dislocationwould according to an averagefiber composition(Winters result. and Stark, 1985).The minimal active state Aminwas Elbowfunction: T&e model includes the bi-articular taken as the resting active state from Hatze (1981). bicepsbrachii caput longumand tricepsbrachii caput The parametersare given in Table 2, wherethe muslongumand caput breve. This providesonly a limited cular physiologicalcross-sectionalarea (PCSA) was Table 2. Normalized muscular parameters as a function of muscle size, PCSA = physiological crosssectional area in cm2, L,, = total rest length, L, = muscle fiber rest length Description Maximal isometric force Minimal or resting active state Series elastic stretch with maximal isometric force Series elastic shape parameter Hill, maximal shortening velocity Hill, shape parameter Time constant increasing active state Time constant decreasing active state Time constant neural excitation

Parameter

Value

F max Aili” SE,, SE,,, MVvm MVs,, TX T da T “S

100 PCSA 0.005 0.05 L,o 3 5 Lo 0.25

10 50 40

R. Happee andF. C. T. Vander Helm

1182

chosenaccordingto Table 1 and the rest length (L,,) wastaken from Van der Helmet al. (1992).In accordancewith dissectionstudiesthe meanfiber restlength Lo was taken as 33% of LtO for the distal muscles (bicepsbrachii and triceps brachii) and 75% for the other muscles(Karlsson, 1992). Dynamic constraint dynamics

derived with inverse muscular

A new method for estimatingmuscularforce and activation from experimentalkinematicdata waspresentedin Happee(1994).The methodcombinesconventional inversedynamicswith optimization utilizing a dynamic musclemodel.The method usesonly very limited computational power, which makes it a usefultool especiallyfor morecomplexsystems.The net torques/forcesare calculated by using conventional inversedynamics.Then, a solution of the load sharingproblem is determinedby optimization (minimization or maximization) of a cost criterion. This yields an inversedynamic optimization. The solution determinedis suchthat it satisfiesa dynamic physiological constraint on muscularforce. This dynamic constraint is determinedwith the inverseof the nonlinear musclemodeldescribedin the previoussection and thus takes into account the muscleproperties describedin this section. With this model muscular statesand neural inputs are alsoestimated. The criterion for the load sharing problem

The load sharing problem is solved by minimization of a metabolic cost criterion. Such a criterion can representthe metabolic cost of musclecontraction but can also be task dependent.Severalcriteria have beenproposedin the literature (Pedotti, 1978; CrowninshieldandBrand, 1981;Patriarco et al., 1981; Dul et al., 1984a,b; Pedersenet al., 1987).In Happee and van der Helm (1994) experimental support is provided for a quadratic criterion, physiologicalsupport is provided for a weightingof musclevolume (V) and it is shownthat sucha weighting contributes to similar stressesin synergistmuscles.This yields the following criterion:

&Vj (+, >

of the hand. The arm washeld in the sagittal plane, the forearm stayedroughly horizontal and the upper arm ranged betweenabout 20” anteflexion and 15’ retroflexion. As the modelincludesmusculardynamics it is suitablefor analysisof such fast movements, Such experimentswere performed with 18 healthy male subjectsfrom the student population (Happee, 1992b,1993band unpublishedresults).A very tight relation existsbetweenthe kinematic and the EMG timing. As this relation would be obscuredby averaging,individual movementswere analyzed with the shouldermodel.One maximally fast anteflexion(forward) and one retroflexion (backward) movement wereanalyzed.The subjectwasfrom Happee(1993b) and wasselectedbecausehe moved relatively fast as comparedto other subjects(but not extremely fast). The movementswere selectedas to have no abnormalitiesin the kinematicsignals.This meansthat they had a small overshoot and no subsequentfluctuations.The muscleactivity predictedwith the model will be comparedwith surfaceEMG recordedduring thesemovements.SomeEMG resultsfrom the other subjectswill be also discussed.The musclesstudied are: deltoideuspars anterior, media and posterior, latissimusdorsi,pectoralismajor parsclavicularisand thoracalis,serratusanterior, trapeziuspars transversalis,infraspinatusand bicepsbrachii. Bi-polar surfaceEMG wasrecordedwith disposable Ag/AgCl electrodeswith an inter-electrodedistanceof 23mm directly attachedto a pre-amplifier.The signal wasbandpassfiltered from 20-1500Hz, rectified and low-passfiltered (100Hz). TheseEMG signalswere sampledat 250Hz, together with the position and velocity of the hand and the forceexertedby the hand. As input variablesthe methodof inversedynamics requiresthe position,velocity and accelerationsignals for every degreeof freedom.The positionof the hand was recorded. This noisy signal was processedas describedin Happee(1994),yielding kinematic data with a time-stepof 12ms shownin Fig. 2. This processingis essentialbecausethe method of inverse

2

j=l

max j

(1)

Here j is the musclenumber and F/F,,, represents the active state which is related to the calcium ion concentrationand henceto the muscularenergyconsumption. In dynamic conditions, F,,, will be time dependentdue to the physiologicalforce-velocity relation and, whenimplemented,the force-length relation. Kinematic

" ....(._.__, __....-'.A'. /-x_:::;:i

and EMG data

Maximally fast goal directedmovementswereper.l .2 .!i 0 .* .6 tirn*‘tS) formed with a set-up describedin Ruitenbeekand Jansen(1984)and Happee(1992b).With their right Fig.2. Kinematicdata of the anteflexion(---) and the arm the subjectscontrolled a manipulator with one retroflexion( .) movement,P = handposition,V = velocity, A = acceleration. degreeof freedom;a forward or backward movement

Shoulder model dynamics involves multiple differentiations which would give unacceptable effects of measurement noise. The set-up does not constrain adduction/abduction. If necessary the experimenter instructed subjects to keep the arm in the sagittal plane. Therefore the position of the elbow was computed assuming that the arm was held in the sagittal plane. The movements of the scapula and the clavicle could not be measured in these fast movements. Therefore shoulder bone positions were chosen according to static measurements. These measurements were performed with the palpator, an instrument for measuring 3-D positions of bony landmarks (Pronk and Van der Helm, 1991). Shoulder bone positions at zero degrees anteflexion were chosen according to Van der Helm and Pronk (1994). As discussed in the introduction, in static experiments a reasonably consistent relation has been found between the rotations of humerus and scapula. This relation is commonly referred to as the ‘scapulohumeral rhythm’. Based on the results of Inman et al. (1944) and Pronk (1991) we assumedthe medial/lateral-rotation of the scapulato be one-third of the rotation of the humerus.As describedin the section Degrees offreedom the assumedmotions of the scapulaare suchthat contact with the thorax is maintained at the angulusinferior (AI) and the trigonum spinae(TS). This derivation of scapularmovementis describedin detail in Van der Helm (1994b).

RESULTS

Muscular

dynamics

Figure 3 showsthe effect of the dynamic constraint imposedby musculardynamicson the neural inputs and stresses predicted.Both the inversedynamicsolution including muscular dynamics and the conventional inversedynamic solution which ignoresmuscular dynamicsare given. An agonist force is required betweenabout t = 0.05s and t = 0.20s ascan beseen from the unconstrainedsimulation(Fig. 3, right upper

dvnomic

constraint(-),

1183

plot, dotted trace).However, musculardynamicsconstrain the rate of force decrease(Fig. 3, right upper plot, continuous trace. The slowly decaying agonist forceshave to be compensated for by increasedantagonistforces(Fig. 3, lower right plot, continuoustrace). The predictedneural muscularinputs in the left plots of Fig. 3 show the well-known triphasic pattern. A comparisonof the inputs predictedwith and without the dynamic constraint, once more shows the relevanceof musculardynamics. Comparison

of model predictions

Figure 4 showsneural inputs estimatedwith the model and the rectified surfaceEMG of the corresponding muscleparts. With the model,neural inputs wereestimatedfrom the kinematic data of the movementsfrom which theseEMG signalswererecorded. The modelinputs and EMG signalsin Fig. 4 will be compared, also taking into account the EMG recorded in identical experimentsin 18 other subjects (Happee,1992b,1993b;and unpublishedresults). Deltoideus: In Fig. 4a both the simulatedand the EMG activity of the deltoideuspars anterior, media and posterior showthe well-knowndifferencein function betweenthesemuscleparts (e.g. Bassettet al., 1990,Flandersand Soechting,1990).In the simulations the other elementsof the deltoideusshowedan intermediate activity. The agonist EMG activity in the deltoideusparsanterior during anteflexionand in the deltoideusparsposterior during retroflexion correspondswell with the predictedmodelinputs. In the EMG of thesemuscles,antagonistactivity is not very clear in Fig. 4a. However, in other subjectssuch antagonist activity was more pronounced. In the deltoideuspars media no clear activity wasfound in both the model predictionsand the EMG. Latissimus dorsi: Figure 4b showsdifferencesbetween the predictedactivity for three elementsof the latissimusdorsi.The other elementsshowedan intermediatebehavior. The EMG was recordedin a position matching the middle element(lati3). Figure 4b

dynamic

static(.)

with EMG

constraint(

-),

static(.)

$1

1

0

.l

js) time neural inputs .2

.4

.5

.6

0

I

.l

.2

‘3s)

.4

.5

.6

time stress

(N

cm-2)

Fig. 3. Theeffectof thedynamicconstraintimposed by muscular dynamics for twopartsof thedeltoideus showingagonist(uppertraces)and antagonistactivity (lower traces)during retroflexion,left column= neuralinputs,right column= stress,(-) simulatedwith musculardynamics,c . .) simulated withoutmuscular dynamics where theneuralinputsimplydescribes theforcerelativeto themaximalforce.

R. Happee and F. C. T. Van der Helm neural

inputs

:

neural

EMG

-. 1

inputs

:

(mV)

:

0

.l

A 0

.l

.2

.4 timi

.5

-. 1

.6

(s)

.2

.3 time

(s)

timi

(s)

timi

(s)

tim;?

(s)

.4

.5

.6

.4

.5

.6

.4

.5

.6

.4

.5

.6

(4 neurol

-.-

inputs

:

1

zi

0 .l

!2 z

0 Ln .1 5 0

.-; .07vi

z ,I-_ 1

.e %

.07

5 o

.I

.2

.3. time

.4

.5

EMG

(s)

:

0

.l

o. -.l

.6

CmV)

.2

(b) neural

inouts

:

neurol 2.2

c1

4

3 EMG

inputs

:

(mV)

:

0

.l

ia. -.l

0

.l

.2 timi

-. 1 (4

.4

.5

.6

-.l

.2

Cs)

Cc)

0

.l

.2 timi

(s)

.4

.5

.6

Fig. 4(a-d).

-.l

0

.l

.2

Shoulder

1185

model

neural

EMG I

-_1

t-3

0

neural

.1

inputs

.2

time

.4

.3

.5

-_1

.6

ts)

:

z E

inputs

:

(mV)

:

0

.l

neural

inputs

.5 r

.2

timl

(s)

.4

.5

.6

:

0 EMG

-.I

(mV)

:

0

.l

.2

neural

-0 .o n 2 .u n u) & .;

.3 time

(f) inputs

EMG

0

-.l

.5

.6

(mV) I

:

0

.l

I

.2

.4

.5

.6

.4

.5

.6

time3 (s)

:

.6 0 H-L .a 0 .2

.4 (s)

.E 4 IEMG a.’ u1 e I 0 -.l

n.-

(mV)

:

0

.l

.2 EMG r

.2

.4 timi

kid

.5

.6

(s)

-.l

(mV)

:

0

.l

.2 tim3

(s)

Fig. 4. Validation: comparison of neuralinputsestimated with the model(uppertraces), andtherectified EMG in mV(lowertraces),for anteflexion(left)andretroflexion(right):(a) deltoideus, pars anterior (delcl), media(dels6) and posterior (dels2); (b) latissimus dorsi; (c) pectoralis major, pars clavicularis (pccl3 and EMG) andthoracalis(pcth2andpcth4);(d) serratus anterior; (e) trapezius pars transversalis (within pars scapularis);

(f) infraspinatus;

(g) biceps brachii, biclo = caput longum, bicbr representing both caput longum and breve.

showsa good correspondencebetweenlati3 and the EMG. The small EMG burst around t = 0.05s during antetlexion is not reproducedby the model but such activity was not found either in most other subjects. Pectoralis major: Figure 4c shows the predicted inputs for three representativeelementsof the pectoralis major. Apparently, in the simulation,the pars clavicularishardly showsagonistactivity during anteflexion wherethe EMG doesshowsuchactivity. Resultsfrom other subjectsindicate a similar activity in the pars clavicularis and thoracalis of the pectoralis major which correspondswell to the averageof the neural inputs in Fig. 4c. Serratus anterior: Figure 4d demonstratesa good correspondence of the predictedinputswith the EMG of the serratusanterior. Both contain a clear triphasic pattern with agonist activity during anteflexion and antagonist activity during retroflexion. In Happee

= caput

breve, biceps

= EMG

(1992b),in the thoracoscapularmuscles,antagonist patterns were rarely found. However, in later measurements (Happee, 1993b;and unpublishedresults)in mostsubjectsantagonistactivity wasfound in the serratusanterior during retroflexion. Trapezius: Figure 4e showsthe predicted input of a part of the trapeziuswhich shouldmatch the EMG of the trapezius pars transversalis.Theseare in reasonableagreementduring anteflexion. During retroflexion the first EMG burst is well predicted by the model but the secondis not. As a whole, the predictions for the trapezius are hard to validate as large differencesbetween subjectswere observed in the EMG of the three parts recorded(trapeziuspars descendens, transversalisand ascendens). Infiaspinatus: Figure 4f showsthat during both anteflexion and retroflexion activity of the infraspinatusis predicted nearthe movementonset.Such activity is also found in the EMG in Fig. 4f and in

1186

R. HappeeandF. C. T. Vander Helm

other subjects.Thus the model supportsthe experimentalfinding that the infraspinatuscannot besimply classifiedas an anteflexor or a retroflexor. The predicted activity around t = 0.2 during retroflexion is not well matchedby the EMG in Fig. 4f and in other subjects. Bicepsbrachii: Figure 4g showsthe predicted activity in the biceps brachii which correspondswell with the EMG only during anteflexion.This indicates a limited value of the simulationresultsfor the biceps. The predictions in the biceps brachii and triceps brachii alsoare of a limited value sincethe elbow is incompletelydescribedin the current model. Overview of stresses and forces

Figure 5 showsall musclestresses. Most stresses are relatively low compared to the maximal isometric stressof 100Ncm-‘. This is partially due to the force-velocity relation which causesa reduction of the maximal stressto about 40 Ncm-* in several muscleparts. Furthermore, the predictedmuscleinputs are largely submaximalasis illustrated in Fig. 4. The largest stressis predicted in the bicepsbrachii caput longum which is the only musclein which maximal inputs are estimated(Fig. 4g). Reaction forces and constraints: Figure 6 showsthe reactionforcesat the scapulothoracicgliding plane,in the conoid ligamentand in the glenohumeraljoint. As can be seenin Fig. 6, the contact force at the angulus inferior (AI) is positive throughout both the anteflexion and the retroflexion movement.This indicates that the inequality constraint that this force must be at least zero is not active; the constraint does not affect the solution.The reaction force at the trigonum spinae(TS) is positive in the first phaseof the anteflexion movementand in a later phaseof the retroflexion movement.The inequality constraint is active for the rest of the time.Thus in the antagonistphaseof the anteflexion movement and the agonist phaseof the retroflexion movement,muscleforceshave to be adaptedto let the force at the TS be at leastzero and thus to maintain contact betweenthe thorax and the scapula.In these phases,retroflexor activity is required for either decelerationor acceleration.To assessthe muscularcontributions to the TS force, the modelwasalso simulatedwithout the TS constraint. Comparingthe resultswith the normal simulations,it was found that considerablyincreasedactivity was required to satisfy the TS constraint, mainly in the trapeziusparsscapularisand in the levator scapulae. Furthermore, without the TS constraint, the pectoralis minor remainedcompletely inactive. The TS constraint had negligibleeffectson the estimatedactivity of the first 5 (lower) parts of the serratusanterior, which showanteflexor activity asillustrated in Fig. 4d. The upper part (part 6) showedagonistactivity during retroflexion. This activity wasfound only whenthe TS constraintwasactive, and thusapparently contributesto the TS force. Such a contribution of

the upper part (parts 5,6) of the serratusanterior to the TS stability constraint was also found during static abduction (Van der Helm, 1994b). The tensionforce in the conoid ligamentwaspositive throughout both the anteflexion and the retroflexion movement.Only in the rest position (15’ retroflexion) before anteflexion wasa zero conoid force found. Glenohumeral constraint: The trace GH in Fig. 6 describesthe glenohumeralreactionforce. Large contact forcesare estimatedthroughout both the anteflexion and the retroflexion movement. The traces theta and phi in Fig. 6 describethe direction of the glenohumeralreaction force. This direction is describedmoreclearly in Fig. 7 which alsoillustratesthe stability of the glenohumeraljoint. This stability was included as a nonlinear inequality constraint in the model. Figure 7 shows the behavior of the glenohumeral constraint during anteflexion. The reaction force is first directedcloseto the ellipsedescribingthe rim of the glenoidcavity and then movesto a more central position. During retroflexion the force describesa similar path in the oppositedirection. The constraint requiresthat the joint reaction force stays within the ellipse.As illustratedin Fig. 7 the ellipseis nevermet and thus the constraint is not active; it does not affect the solution. In preliminary simulations with different model parametersthe constraint becameactive during severalphasesof the movements. Omitting the constraint yieldeda glenohumeralforce in a direction almost tangential to the glenoidjoint surface.Thesepotentially dangerousconditions occurred at timeswhenthe overall force levelswerevery low. In these conditions only slight modifications of muscle forces were required to satisfy the glenohumeralconstraint. When this also holds for other conditions it will be very difficult to detect such modificationsexperimentallyfrom the EMG. Thus it will be practically impossibleto discriminatebetween the muscularcoordination of the patient with habitual luxation and the normal coordination. Energy and scapular motion

The mechanicalenergydeliveredand dissipatedby musclesis presentedin Table 3. Until now (Figs 3-7) a scapular movement according to the scapulohumeral rhythm was assumed.Table 3 also describesresultsof simulationsperformedwith a fixed scapula;its position waskept constant but the movement equationsspecifyingequilibrium of the scapula werestill to be satisfied.The scapularmotion strongly affects the mechanicalenergy delivered by muscles (Table 3).Energy is mainly deliveredto acceleratethe limb and dissipatedto deceleratethe limb. This was verified from plots of the mechanicalpower delivered by individual muscleelementsas a function of time, where it emergedthat during accelerationvery little power is dissipated.In the simulation with scapular motion, the thoracoscapularmusclesdeliver a considerableamount of energy. Of course,with a fixed

Shoulder model

G

20

tropezius

scopul.

serrotus F

-I’” ;

10

t

~~~

onterior

1

‘1

20

clovic.

0

c-r (v I E

20

,:-.. ‘..

s

lzif!Gl levotor

40

;

10

pect

mojor

thor.

pect

mojor

clov.

10

0 a subscopuloris 30

corocobrochiolis t

0 20

20

k

1

10 0

infrospinotus

m

20

.:; ‘A. ;‘\:. ._.. .i ‘i. .-.Q.,

20

rhomboideus

20

suprospinotus

60

10

0

G

clovic

20

0

dorsi

2o

0 deltoideus

10

-...._

10

0

10

0 iIizl lotissimus

major

t

scopuloe

10

t

teres

20

10

I uL

scopul

20

*.. “.:....:. i:::;;-+. .. _,_ . ... _. 0 m 0 . minor. pectorolls 2. 20

3

triceps

10

deltoideus

10

t

minor

0 tropezius 3

” fL

teres

20

0

G

1187

1

4o

biceos

60

Ion-

brv.

40 20 e

0

20

h hl

0

.2

.4

time

(s)

tropezius

200 .6

scopul.

0

.2 time

serrotus

20

.4 (s)

0 ~ 0

.6

onterior

-

20

.2

.4

time

(s)

teres

.6

minor

20

triceps

1

;

10

10

10

10

....


0

h

20

;

10

z

0

c-4

0 m teres 20

0 tropezius

clovic.

deltoideus

20

scopul

10

“.__, 0 major

20

10

Lriizl lotissimus

dorsi

10 .g\

/-.

t-4

20

0 levotor

scapulae

0 deltoideus I

20

clovic

I 40

iiz_l suprospinotus

0

I

lIii!zl oect

*. moior

thor.

2o

1 ;

10

3

0

cs

20

r-4

;

10

3

0

A

r-4

k ” z -

20

10

pectorolls

minor,

2.

20

;orocobrochiolis

10

,

; ;K

0 I-AA ’ infrospinotus 30 c

rhomboideus I

I so

biceps -

Ion-

brv.

10 10 0

0

.2

.4

time

(s)

.6

0 0

.2 time

.4 (s)

.6

0

.2 time

Fig. 5. Muscle stresses predicted, (a) anteflexion, (b) retroflexion: (-) element(s); (-. -. -) last element(s).

.4 (s)

.6

first element(s); (- - -) middle

R. Happee and F. C. T. Van der Helm

1188

9 2

-40

I .l

0

1 2

I .3

time

/ .4

= .ca .5

0 -40

.6

0

L .1

I .2

(f)

I .3

time

I .4

I .5

.6

(s)

Fig. 6. Reaction forces, (TS) trigonum spinae; (AI) angulus inferior; (conoid) conoid ligament tension force; (GH) glenohumeral resultant force; theta and phi describe the direction of the glenohumeral force as illustrated in Fig. 7.

change and therefore deliver more energy. These effects of scapular motion have alsobeenevaluatedwith respect to the shortening velocities. With scapular movement, the lower element of the rhomboideus

40

0

attained about 65% of its maximal shorteningvelocity which is much more than the 20-30% at which muscular energy generation is most efficient. The shorteningvelocity of most other muscleswasabout 25%. The scapularmotionsaffect the metaboliccost

20

E s

0

0

z s

-20

required -40

-50 phi

0 (degrees)

50

to obtain

the desired forces. This cost is

representedby the criterion optimized [Equation (111. The summedcriterion over the whole anteflexion movementis increased13% when no scapularmovement is allowed.

Similarly

during

retroflexion

the cri-

terion is increased6%. This indicatesthat the asFig. 7. The glenohnmeral constraint during antellexion. The sumed scapular movement reducesthe metaboliccost. curve describes the intersection of the joint reaction force Besidesthe shorteningvelocitiesthe movementsof vector with the articular surface of the glenoid. The ellipse the scapulaalsoaffect the muscularmomentarmsand describes the border of the glenoid cavity, as viewed from a lateral position. Axes are in degrees, (right) anterior, (up) the stability of the glenohumeraljoint. From separate cranial. (0) the initial position, (x) the final position. simulationsit emergedthat the moment arms affect the summedcost by lessthan 0.7%. The glenohumeral stability is positively affectedby the scapularmotion: In the retroflexed position, marked 0 in Fig. 7, with scapula the thoracoscapular muscles cannot deliver scapularmotion no force modificationsare required any energy as their length is constant. With a fixed to obtain a stableglenohumeraljoint. Without scapscapula, the muscles with an origin on the scapulaand ular motion the force movesto the ellipseand minor insertion on the humerus have an increased length force modificationsare required to stabilizethe joint.

Shoulder model

1189

Table 3. Mechanical energy in Nm, delivered and dissipated during antegexion (left) and retroflexion (right), the value in brackets ( ) is computed without scapular motion, the energy is computed as the product of force and shortening velocity integrated over the whole movement Anteflexion Muscle (part)

Delivered

Dissipated

Delivered

Dissipated

Trapezius pars scapularis Trapezius pars clavicularis Levator scapulae Pectoralis minor Rhomboideus Serratus anterior Subtotal scapular muscles

0.35 (0) 0.15 (0)

0.16(O) ,

0.18 (0)

0.20 (0) 0.13 (0)

0.24(O) 0.01 (0) 0.57(O)

0.30(O) 0.01 (0) 0.36(O) 0.03 (0)

1.45 (0)

0.98 (0)

0.88 (0)

1.78 (0)

Infraspinatus Supraspinatus Subscapularis Teres minor Subtotal rotator cuff muscles

0.47 (0.74) O.Zl(O.56) O.OS(O.18)

0.02 (0.06) 0.03 (0.03) O.OS(O.08) 0.06(0.09)

0.03 (0.06) 0.04 (0.04) 0.22(0.24) 0.04(0.13)

O.SO(O.78) 0.17(0.53) 0.05 (0.11)

0.76(1.48)

0.19 (0.26)

0.33 (0.47)

0.62 (1.42)

Deltoideus pars scapularis Deltoideus pars clavicularis Coracobrachialis Teres major

0.03 (0.08) 0.18 (0.50) 0.15 (0.09)

1.15 (1.26)

1.12(1.32)

0.08 (0.15) 0.16(0.38) 0.22(0.14)

1.06(1.95)

0.89 (1.05)

Biceps brachii caput longum Biceps brachii caput breve Triceps brachii caput longum Latissimus dorsi Pectoralis major pars thoracahs Pectoralis major pars clavicularis Total of all muscles

O.Sl(1.24) 0.02(0.02)

1.09 (0) 1.59 (0)

0.72(1.15) 0.04 (0.06) 0.20(0.21) 0.30(0.17)

3.95 (3.83)

The model

Three-dimensionalmodelsof shoulderkinematics have beenpresentedby Engin and Peindl (1987),by Pronk (1989,199l) and by Hogfors et al. (1991). Muscle forceshave been included by Hdgfors et al. (1987)and by Karlssonand Peterson(1992).A threedimensionaldynamic shouldermodel,includinginertia and muscleforceswasrecently developed( Van der Helm, 1994a,b;Van der Helm and Veenbaas,1991; Van der Helm et al., 1992;Veegeret al., 1991).In the current paper,musculardynamicshavebeenincluded shoulder

model.

This

is done

0.31(0.39) 0.02 (0.08)

0.01 (0.02)

Fast goal directed arm movementsin the sagittal planewere analyzed with a three-dimensionalshoulder model.The simulationsprovided estimatesof the activity of shouldermusclesincluding the activity of deeperlying muscles.Table 3 for instanceshowsthe largeamount of energycontributed by the teresmajor to retroflexion.

this

0.24(0.36) 0.54(0.60)

0.40 (0.40)

DISCUSSION

into

Retroflexion

with

the

method of inverse dynamic optimization including musculardynamicsdevelopedby Happee(1994). The activity estimatedwith the model has been comparedwith surfaceEMG measurements (Fig. 4). EMG measurements are informative with respectto timing rather than amplitude of muscular activity.

3.88 (3.85)

4.00(3.80)

4.05 (3.77)

Most model resultswere similar to the EMG but somedifferencesbetween model results and EMG were found in muscleswhere the EMG activity was subjectdependent. As describedin the section Kinematic and EMG data, the assumedmovementsof the scapulaand the claviclewerebasedon static measurements performed on other subjects.For dynamicmeasurement of these movements, the possibility of three-dimensional roentgenis currently assessed in our group. The consequence of the assumed scapularmovementhasbeen evaluated by performing an additional simulation without scapular movement. A marked effect was found on the muscularshorteningvelocitiesand on the mechanicalenergy contributions (Table 3). Furthermore, in the simulation without scapularmovement the summedmetaboliccostcriterion [Equation (l)] wasincreased13% for the anteflexionmovement and 6% for the retroflexion movement.This indicates that the scapularmovementassumedin the normal simulationreducesthe metaboliccost. Muscle elements: Becausemany shouldermuscles have a complex three-dimensionalgeometry, muscle fibershave different linesof action and thereforehave different mechanicaleffects.Van der Helm and Veenbaas(1991)analyzedthe mechanicaleffect of shoulder musclesand showedthat, dependingon the muscle geometry,at most 6 elementsper muscleare required

1190

R. Happee

and F. C. T. Van der Helm

to obtain a sufficiently accurate description of the mechanical effect of muscles (Table 1). Figures 4 and 5 show large differences between the predicted activity in different elements of the same muscle. As these forces are predicted by optimization of a mechanical criterion, this shows the mechanical relevance of modelling so many muscle elements. Another point is to what extent the nervous system can selectively activate muscle elements. The motor-neuron pool and the motor-unit distribution in muscles may constrain this selectivity. In the solution of the load sharing problem it may be relevant to include neural constraints on the differences between activations of muscle elements. On the other hand, the current results contain a rather gradual change of the activation through muscles: No extreme differences were predicted between the activity of adjacent muscle elements. Muscular dynamics: The inclusionof musculardynamicsinto the modelhad a considerableeffecton the neural inputs and forces predicted (Fig. 3). In the EMG a large differencein timing was observedbetweentwo synergistmuscles,the latissimusdorsi and the deltoideuspars posterior. Activity of the latissimusdorsi precededthe deltoideusparsposterior by up to 62 ms (Table 3 in Happee, 1992b).In Happee (1992b)it wassuggested that this phenomenoncan be explained by the fact that the latissimushas longer musclefibersand a longertendon than the deltoideus pars posterior and therefore (probably) exhibits a slower responseof force to activation. The earlier activation of the latissimuswassuggested to bea compensation,generatedin the nervous system,for the slowerresponseof the latissimus.In the modelresults no relevant differencein timing hasbeenfound. This can be explained by the fact that all musclesare modelledas having equal active state dynamicsand neuralactivation dynamics(Table 2). The mechanical dynamicsaredeterminedby the ratio of the maximum shorteningvelocity and the serieselasticstretch. As both these parameterswere taken proportional to musclelength, in the current model, musclelength does not affect muscular dynamics, which seems a point of interest in further modelling. muscles The mechanicalenergy contributions in Table 3 showthat thoracoscapularmusclescontribute strongly to the accelerationand decelerationof the limb. The serratusanterior delivers 28% of the positive energy during anteflexion and dissipates35% of the energyduring retroflexion. The rhomboideusdelivers 9% of the positive energy during retroflexion and dissipates15% of the energy during anteflexion. Furthermore, considerable activity of thoracoscapularmusclesis shownto be requiredto satisfythe TS constraint. This constraint representsthat the trigonum spinaeof the scapulahasto stay in contact with the thorax. This constraint wasalsoshownto be active during abduction (Van der Helm, 1991).Contact betweenthe thorax and the scapulais very relThoracoscapular

evant, as this contact makesthe shouldera closedchain mechanism (thorax-clavicle-scapula-thorax). Due to this closed-chainthe clavicle and the scapula provide a mechanically stable base for the glenohumeraljoint. In Happee(1992b)in the EMG of thoracoscapular musclesantagonistactivity wasrare. This left unclear how deceleratingforcesinducedby the prime movers are transferred to the trunk. Meanwhile, measurementsin other subjectsshowedantagonistactivity in the serratusanterior during retroflexion asis alsopredicted by the model (Fig. 4d). During anteflexion in thoracoscapularmusclesno clear antagonistpatterns were simulated. Yet the deceleratingforces can be transferred to the trunk by the more complicated activity of, for instancethe trapeziuspars scapularis, wherethis activity is alsorelatedto the TS constraint. Stability

of the glenohumeral

joint

Van der Helm(1994b)showedthat modificationsof muscle forces are required to stabilize the glenohumeraljoint during posturemaintenance.He found increasedforcesin the so-calledrotator cuff muscles, four musclesinserting as a half circle close to the glenohumeraljoint. These results support the conclusionthat rotator cuff musclesactively stabilizethe glenohumeraljoint. The current paper presentsthe first evaluation of glenohumeralstability in dynamic conditions.In the goal directedmovementsanalyzed, considerablemuscleforcesare required to accelerate and deceleratethe limb. However, it was found that the glenohumeralconstraint is not active during acceleration and deceleration;the constraint does not affect the solution. Apparently, the force distribution resulting from the optimization of the criterion [Equation (l)] provides a stablejoint. Thus the current resultsprovide no indication that voluntary accelerationsare critical with respectto humerusdislocation. This observationis,however,basedon simulations in a smallmovementrangeand certainly cannot be extended to extreme movements like baseball pitching. A relatively largestresswasfound in 3 of the 4 rotator cuff muscles:the infraspinatus,supraspinatusand subscapularis(Fig. 5). This large stressis a result of the costcriterion optimized[Equation (l)] containing a weighting of musclevolume (motivated in Happee and van der Helm, 1994).The volume is proportional to muscle fiber length, thus a larger cost will be attributed to stressin longermuscles.Minimization of the cost criterion will therefore attribute a relatively largestressto shortermuscleslike thoseof the rotator cuff. This accountsfor agonistactivity in rotator cuff muscles.Here agonistactivity is definedascontributing somehowto the net torquesrequiredat a certain instant. It hasbeensuggestedthat even antagonistic cuff activity should be required to stabilize the joint (Van der Helm, 1994a).Antagonist activity will never be predictedby any criterion of a form like equation (1). Antagonist activity can be enforced by the

1191

Shoulder model

glenohumeralconstraint but this wasnot the casein the current simulations as the constraint was not active; the constraint did not affect the solution. According to Table 3 rotator cuff musclesdeliver considerable positive energy: The infraspinatus and supraspinatuscontribute to anteflexion and the subscapularisto retroflexion. This once more indicates agonistactivity of rotator cuff muscles.

Inman, V. T., Saunders, J. B. and Abbot, L. C. (1944) Observations on the function of the shoulder joint. J. bone Jt surg. 26a, l-30.

Karlsson, D. (1992) Force distributions in the human shoulder. Ph.D. Thesis, Giiteborg. Karlsson, D. and Peterson, B. (1992) Towards a model for force predictions in the human shoulder. J. Biomechanics 25, 189-199. Lippitt, S. and Mattsen, F. (1993) Mechanism of glenohumeral joint stability. Clini. Orthop. Related Res. 291, 20-28.

REFERENCES

Basset, R. W., Browne, A. O., Morrey, B. F. and An, K. N. (1990).Glenohumeral muscle force and moment dynamics in a position of shoulder instability. J. Biomechanics 23, 40-415 Crowninshield, R. D. and Brand, R. A. (1981) A physiologically based criterion of muscle force prediction in locomotion. J. Biomechanics 14, 793-801. Dul, J., Townsend, M. A., Shiavi, R. and Johnson, G. E. (1984a) Muscular synergism--I: On criteria for load sharing between synergistic muscles. J. Biomechanics 17, 663-673. Dul, J., Johnson, G. E., Shiavi, R. and Townsend, M. A. (1984b) Muscular synergism-II. A minimum-fatigue criterion for load sharing between synergistic muscles. J. Biomechanics

17. 675-684.

Engin, A. E. (1980) On the biomechanics of the shoulder complex. J. Biomechanics 13-7; 575-590. Engin, A. E. and Peindl, R. D. (1987). On the biomechanics of human shoulder complex--I. Kinematics for determination of the shoulder complex sinus. J. Biomechanics 20, 103-l 17. Flanders, M. and Soechting, J. F. (1990) Arm muscle activation for static forces in three dimensional space. J. Neurophysiology 64, 1818-1837. Happee, R. (1992a) Time optimality in the control of human movements. Biological Cybernet. 66, 357-366. Happee, R. (1992b) Goal-directed arm movements. I. Analysis of EMG records in shoulder and elbow muscles. J. Electromyogr. Kinesiol. 2, 165-178. Happee, R. (1993a) Goal-directed arm movements. II: A kinematicmodel and its relation to EMG. J. Electromyogr. Kinesiol. 3, 13-23. Happee, R. (1993b) Goal-directed arm movements. III: Feedback and adaptation in response to inertia perturbations. J. Electromyogr. Kinesiol. 3, 112-122. Happee, R. (1994) Inverse dynamic optimization including muscular dynamics, a new simulation method applied to goal directed movements. J. Biomechanics 27, 953-960. Happee, R. and Van der Helm, F. C. T. (1994) Criteria representing the cost of muscular contraction, used to solve the load sharing problem in inverse dynamic optimizations (submitted). Hatze, H. (1981) Myocybernetic control models of skeletal muscle, characteristics and applications. University of South Africa, ISBN 0 86981 216 5. Hogfors, C., Sigholm, G. and Herberts, P. (1987) Biomechanical model of the human shoulder-I: Elements. J. Biomechanics 20, 157-166. Hiigfors, C., Peterson, B., Sigholm, G and Herberts, P. (1991) Biomechanical model of the human shoulder joint-II: The shoulder rhythm. J. Biomechnnics 24, 699-709.

Patriarco, A. G.. Mann, R. W., Simon, S. R. and Mansour, J. M. (1981) An evaluation of the approaches of optimization models in the prediction of muscles forces during human gait. J. Biomechanics 14, 513-525. Pedersen, D. R., Brand, R. A., Cheng, C. and Aurora, J. S. (1987) Direct comparison of muscle force predictions using linear and nonlinear programming. J. biomech. Engng 109, 192-198. Pedotti, A., Krishan, V. V. and Stark, L. (1978) Optimization of muscle-force sequencing in human locomotion. Math. Biosci. 38, 57-76. Pronk, G. M. (1989) A kinematic model of the shoulder girdle: a resume. J. Med. Eng. Technol. 13, 119-123. Pronk, G. M. (1991) The shoulder girdle: analysed and modelled kinematically. Ph.D. Thesis, Delft University of Technology, The Netherlands, ISBN 90-370-0053-3. Pronk, G. M. and Van der Helm, F. C. T. (1991) The palpator, an instrument developed for measuring the 3-D positions of bony landmarks in a fast and easy way. J. Med.

Engng

Technol.

15, 15-20.

Ruitenbeek, J. C. and Jansen, R. J. (1984) Computer controlled manipulator display system for human movement studies. Med. Biological Engng Comput. 22, 304-308. Van der Helm, F. C. T. (1994a) A finite element musculoskeletal model of the shoulder mechanism. J. Biomechanics 27, 551-569. Van der Helm, F. C. T. (1994b) Analysis of the kinematic and dynamic behavior of the shoulder mechanism. J. Biomechanics

27, 527-550.

Van der Helm, F. C. T. and Veenbaas, R. (1991) Modelling the mechanical effect of muscles with large attachment sites: application to the shoulder mechanism. J. Biomechanics 24, 1151-1163. Van der Helm, F. C. T. and Pronk, G. M. (1994) Threedimensional recording and description of the motions of the shoulder mechanism. J. biomech. Engng. (accepted). Van der Helm, F. C. T., Veeger. H. E. J., Pronk, G. M., Van der, Woude L.H.V. and Rozendal, R.H. (1992) Geometry parameters for musculoskeletal modelling of the shoulder mechanism. J. Biomechanics 25, 129-144. Veeger, H. E. J., Van der Helm, F. C. T., Van Der Woude, L. H. V.. Pronk, G. M. and Rozendal, R. H. (1991) Inertia and muscle contraction parameters for musculoskeletal modelling of the shoulder mechanism. J. Biomechanics 24, 615-629.

Winters, J. M. and Stark, L. (1985) Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models. IEEE Trans Biomed Engng. BME-32-10,

826-839.

Winters, J. M. and Stark, L. (1988) Estimated properties of synergistic muscles involved in movements of a variety of human joints. J. Biomechanics 21, 1027-1041.