ORIGINAL ARTICLE
Development of a comprehensive musculoskeletal model
of the shoulder and elbow
A. Asadi Nikooyan • H. E. J. Veeger •
E. K. J. Chadwick • M. Praagman •
F. C. T. van der Helm
Received: 7 October 2010 / Accepted: 10 October 2011
 The Author(s) 2011. This article is published with open access at Springerlink.com
Abstract The Delft Shoulder and Elbow Model (DSEM),
a musculoskeletal model of the shoulder and elbow has
been extensively developed since its introduction in 1994.
Extensions cover both model structures and anatomical
data focusing on the addition of an elbow part and muscle
architecture parameters. The model was also extended with
a new inverse-dynamics optimization cost function and
combined inverse-forward-dynamics models. This study is
an update on the developments of the model over the last
decade including a qualitative validation of the different
simulation architectures available in the DSEM. To vali-
date the model, a dynamic forward flexion motion was
performed by one subject, of which the motion data and
surface EMG-signals of 12 superficial muscles were mea-
sured. Patterns of the model-predicted relative muscle
forces were compared with their normalized EMG-signals.
Results showed relatively good agreement between forces
and EMG (mean correlation coefficient of 0.66). However,
for some cases, no force was predicted while EMG activity
had been measured (false-negatives). The DSEM has been
used and has the potential to be used in a variety of clinical
and biomechanical applications.
Keywords Shoulder  Elbow  Musculoskeletal model 
Inverse and forward dynamics  Validation  Muscle force 
EMG
Abbreviations
lf Fiber length
ls Sarcomere length
lopt Optimal fiber length
h Reference joint angle
_h Reference joint angular velocity
€h Reference joint angular acceleration
Lm Muscle length
M Net joint moment
rmax Maximum muscle stress
Fm Predicted muscle force
PCSA Muscle physiological cross sectional area
Fmin Minimum permissible muscle force in the inverse
optimization
Fmax Maximum permissible muscle force in the inverse
optimization
e Excitation dynamics
u Hypothetical neural input of the forward muscle
model
a Active state
Lce Length of contractile element (CE)
Mc Correction moment
hc Calculated joint angle
_hc Calculated joint angular velocity
Electronic supplementary material The online version of this
article (doi:10.1007/s11517-011-0839-7) contains supplementary
material, which is available to authorized users.
A. Asadi Nikooyan (&)  H. E. J. Veeger 
F. C. T. van der Helm
Department of Biomechanical Engineering,
Delft University of Technology, Mekelweg 2,
2628 CD Delft, The Netherlands
e-mail: a.asadinikooyan@tudelft.nl
H. E. J. Veeger  M. Praagman
Research Institute MOVE, VU University Amsterdam,
1081 HV Amsterdam, The Netherlands
E. K. J. Chadwick
Department of Sport and Exercise Science,
Aberystwyth University, Aberystwyth,
Ceredigion, UK
123
Med Biol Eng Comput
DOI 10.1007/s11517-011-0839-7

1 Introduction
Biomechanical models can give insight into the mechanical
basis of musculoskeletal function. In the last few decades, a
variety of models of the entire human musculoskeletal
system, from simple two-dimensional (2D) models [1–5] to
complex three-dimensional (3D) models [6–9], have been
developed. However, of all these models, not many
describe the upper extremity. A major reason for this is the
complex kinematic structure of the upper limb. The many
degrees-of-freedom (DOF) of the shoulder girdle limit the
usefulness of simple 2D models and lead to complex 3D
models.
For clinical applications, sophisticated models are nee-
ded. Such models should be complex enough to realisti-
cally replicate the behavior of the human musculoskeletal
system. Few complex upper extremity models have been
developed such as the Swedish model [10, 11] based on the
model of Hogfors et al. [12, 13], the Newcastle shoulder
model [14], the shoulder part of the AnyBody Modeling
System [6], the Stanford model implemented in SIMM
[15], and the Delft Shoulder Model (DSM) which is the
core of this study.
The Delft Shoulder Model as first described in 1994 [9]
is a comprehensive 3D inverse-dynamic model of the
shoulder complex in which the recorded motions of the
bony segments and external loads are used as input to the
model and muscle and joint contact forces, muscle lengths,
and moment arms are calculated as model outputs through
an inverse-dynamics analysis. To qualitatively validate the
model, estimated force–time curves were compared to
measured EMG signals [16] which showed good agreement
in the timing of muscle activations. Data for the original
model were taken from [17–19]. Later, elbow data were
added based upon a follow-up cadaver study [20]; conse-
quently, the model was renamed to the Delft Shoulder and
Elbow Model (DSEM).
Following a detailed cadaver study on the shoulder [21]
and elbow [22], information about muscle architecture and
optimal fiber length was additionally obtained. It was
expected that this addition would lead to improvements in
the prediction of muscle forces and better insight into the
functioning of specific muscles since force–length and
force–velocity relationships could be implemented. By
including the muscle dynamics, some modifications and
extensions were carried out in the model. First, the inverse-
dynamics model was modified in such a way that the
muscle dynamics were taken into account as constraints on
the maximum permissible muscle force during the inverse
optimization. Second, combined inverse-forward-dynamics
versions of the DSEM [23, 24] could be developed. Third,
a new muscle load sharing cost function for inverse opti-
mization namely the energy-based criterion [25] was
introduced and implemented in the model. This new cost
function is based on the energy-consuming processes in a
muscle needed to produce a contraction.
Although the DSEM has been widely used in a number
of studies, it was not individually addressed in the litera-
ture. The aim of this article is to provide all aspects and
developments of the model since its original introduction in
1994, including the measured elbow data and a qualitative
validation of the three different simulation architectures
available in the DSEM (inverse dynamics optimization—
IDO, inverse-forward-dynamics optimization—IFDO, and
IFDO with controller—IFDOC). The model simulations
will be based on a new anatomical dataset and the appli-
cation of an energy-based load sharing cost function
enabled by the addition of the muscle architecture param-
eters. As such, it fits in the developments sketched in an
Editorial published by Cutti and Veeger [26] in a special
issue of this journal on shoulder biomechanics in which it
was concluded that there is a need for the general biome-
chanical (upper extremity) models to be more thoroughly
validated and tested and in which new cost functions
should be included.
2 Materials and methods
2.1 Anatomical data
The geometrical data for the DSEM were taken from studies
on the shoulder [21] and elbow [22] from the same specimen,
a 57-year-old muscular male cadaver. In these studies a total
number of 31 muscles of the shoulder (23 muscles) and
elbow (8 muscles) were divided into 139 elements. Joint
surfaces and other bony contours were digitized for model-
ing using geometrical forms. Muscle architecture parameters
including tendon length, physiological cross-sectional area
(PCSA), pennation angle, and the fiber length (lf) were
measured. The sarcomere lengths (ls) were also measured
using a laser-diffraction technique [27]. Assuming an opti-
mal sarcomere length of 2.7 lm [28], the optimal fiber length
(lopt) for a muscle was calculated as:
lopt ¼ 2:7
lf
ls
ð1Þ
Since the elbow data have not yet been published other
than in an internal report [22], these data are partly
presented here. These data include the position of bony
landmarks (Fig. 1; Table 1), bony contours for muscle
wrapping (Table 2), and functional axes of rotation of the
elbow (Table 3). For the description of the measurement
methods, the values of the muscle parameters (PCSA, mass
and optimal fiber length), and the relative muscle force-
sarcomere length curves see supplementary materials.
Med Biol Eng Comput
123

2.2 Kinematics
The DSEM is a finite-element model which has been
implemented by building on the software packet SPACAR
[29] for the analysis of spatial multi-body mechanisms. For
detailed descriptions of model kinematics and implemen-
tation in SPACAR see Ref. [9].
The total number of DOF of the model is 17, namely six
for the thorax (which is considered as the moving base),
three for the shoulder girdle, three for the glenohumeral-
joint, two for the elbow, and three for the wrist.
There are several options to provide input to the model
among which using the joint angles is the most popular. To
calculate the glenohumeral-joint rotation center, which is
necessary for reconstruction of the local coordinate system
of the humerus, the instantaneous helical axes (IHA)
method [30, 31] is mostly used, although alternatives such
as the use of regression equations [32] or the SCoRE
method [33, 34] are also possible in the DSEM.
For the definition of the clavicular orientation only two
landmarks are generally available, thus, the axial rotation
of the clavicle is estimated by minimizing the rotations in
the AC-joint [35].
The shoulder girdle is a closed-chain mechanism and the
motions are constrained by such factors as the shape of the
thorax over which the scapula glides, the length of the
conoid ligament, the length of the clavicle, and the size of
the scapula. As such, the motions of the shoulder girdle of a
measured subject cannot be exactly reproduced by the
model due to differences in the geometry between subject
and model. To ensure that all positions input to the model
can actually be assumed by the model, the measured angles
are adjusted slightly to fit the constraints of the model by
minimization of the following cost function [36]:
J ¼ W1 dCxð Þ2þ dCy
 2þ dCzð Þ2
 
þ W2 dSxð Þ2þ dSy
 2þ dSzð Þ2
 
ð2Þ
where dCx and dSx are the differences between the mea-
sured and optimized angles for the clavicle and scapula
around the x-axis, respectively. A similar definition is
applied for angles around the y- and z-axes. W1 and W2 are
weight factors, and were set at 1 and 2, respectively. For
detailed description of the optimization procedure, see the
supplementary materials.
2.3 The inverse-dynamics optimization (IDO)
In the original model (DSM), the inertial forces and
moments were included but muscle dynamics were not. In
the modified inverse dynamics model (IDO, Fig. 2a), the
muscle dynamics has been taken into account as constraint
on the maximum permissible muscle force in the inverse
optimization process. The joint angles and external loads
are used as the model inputs. The outputs of the model
include muscle and joint reaction forces, muscle and liga-
ment lengths, and muscle power.
The filtering and differentiation routines of Woltring
(the GVC method) [37] have been implemented to calcu-
late velocities and accelerations from the inputs.
The load-sharing problem is solved using a nonlinear
optimization process in which a cost function is minimized.
The stress cost function (SCF) which is based on minimi-
zation of the squared muscle stress [38] was originally
implemented in the DSM as the default objective function,
but recently a new energy-based cost function (ECF) [25]
has been implemented. In the energy-based cost function,
the energy consumption due to calcium pumping and cross-
bridge function is taken into account.
The calculated forces in the optimization process for
each muscle element (m) are bounded by the inclusion of
muscle force–length relation where minimum force is
taken as zero and the maximum force is a function of
maximum muscle stress (rmax), PCSA, and sarcomere
length (lsm):
Fmax lsmð Þ ¼ f lsmð ÞPCSAmrmax ð3Þ
where rmax is taken as 100 N/cm
2 [39]. f(lsm) is the nor-
malized muscle force–length relationship defined as a
Fig. 1 Palpable bony
landmarks on the humerus, ulna,
and radius
Table 1 Positions of palpable bony landmarks (BL) on the humerus,
ulna, and radius in the global coordinate system defined in Ref. [21]
BL X (cm) Y (cm) Z (cm)
Epicondyle medialis (EM) 15.66 -30.79 10.56
Epicondyle lateral (EL) 21.53 -30.21 7.15
Olecranon (OL) 18.97 -30.23 10.58
Styloideus Ulnae (SU) 21.54 -55.71 3.39
Styloideus Radii (SR) 17.29 -55.01 0.01
Med Biol Eng Comput
123

Gaussian-type shape function (see [40] for a detailed
description).
The model stability is defined as being maintained when
the joint reaction force vector is directed inside the rim of
the glenoid fossa, modeled by an ellipse.
2.4 The inverse-forward-dynamics optimization
(IFDO)
The combination of an inverse dynamic optimization
approach with inclusion of the muscle dynamics by a for-
ward dynamic muscle model (IFDO) is an efficient way to
obtain dynamically feasible muscle forces and neural
inputs. Happee and van der Helm [41, 42] showed that
inclusion of the muscle dynamics in the inverse optimiza-
tion had considerable effects on the model-estimations of
the neural inputs and individual muscle forces.
The muscle models are required to account for the effect
of muscle electromechanical delays and force–velocity
relationship in an inverse-dynamic optimization. In the
IFDO (Fig. 2b) both forward and inverse muscle models
are used. As for muscle model, a three-component Hill type
model [40, 43] consisting of a second-order activation
dynamics part and a first-order contractile dynamics part is
being used (for a detailed mathematical formulation of the
muscle model see supplementary materials).
At each time-step (i, Fig. 2b), the calculated optimal
muscle forces are constrained by maximum (Fmax,i) and
minimum (Fmin,i) permissible values of the muscle forces
estimated by a forward muscle model with use of the muscle
states of the previous time-step (ei-2, ai,-1, Lce,i-1). At the
same time-step (i), an inverse muscle model is used to
estimate the neural inputs (ei-1, ai, Lce,i) that will be used as
the inputs to the forward muscle model in the next step
(i ? 1). The starting position is assumed to be quasi-static
in which the initial vales of u, e, a, and Lce are estimated
through a steady state equilibrium condition. During the
next time steps, the states are updated in a dynamic opti-
mization procedure. The values of e, a, and Lce are updated
iteratively, while u is estimated analytically.
2.5 The IFDOC
Due to discretization errors and possibly an unstable sys-
tem, the motions calculated by the forward-dynamic part in
the IFDO will not be exactly the same as the recorded
motions which were input to the inverse dynamics part.
Therefore, the IFDO was modified in such a way that the
difference in position and velocity will be fed back to the
inverse-dynamic model. The modified model is called the
IFDOC or simply the IFDOC (Fig. 2c) [23]. At each time
step, the feedback controller will adjust the neural input
signal in the next time step by calculating a correction
moment (Mc) using the errors in angle and angular veloc-
ity. Therefore, a forward-dynamic simulation will be
obtained which should result in exactly the same motion as
the recorded motion. In the forward dynamics part of the
model (Fig. 2c), the forward musculoskeletal model
developed by van der Helm and Chadwick [24] was
implemented. For integration of the motion equations in
forward dynamics analysis, two different algorithms
namely the Adams–Moulton and Euler algorithms have
been implemented in the model. The forward dynamics
simulation based on the Euler method is up to four times
faster than the Adams–Moulton algorithm, but less stable.
In contrast to the IDO model in which each time step is
considered to be independent of the preceding time steps, in
the IFDO and IFDOC analyses each time-step is coupled to the
following time-steps through sets of differential equations.
2.6 Biomechanical applications of the model
The DSEM has frequently been applied. It was used to
study goal-directed movements [42], wheelchair propul-
sion [44, 45], rotator cuff tears [46], tendon transfers [47,
48], loads on the arm [49, 50], rotator cuff changed fol-
lowing scapular neck fracture [51], effect of rotator cuff
dysfunctions on wheelchair propulsion [52], weight trans-
fers in wheelchair users [53], effect of including the neural
activity in the modeling process [54], and stability of ce-
mentless glenoid prostheses [55].
Table 2 Bony contours. [dx dy dz] is the direction of the central axis of the cylinder; [Px Py Pz] is the coordinate of an arbitrary point on the
central axis of the cylinder; R is the radius
Wrapping object Bony structure Px Py Pz dx dy dz R
Cylinder 1 Radius 19.92 -35.32 5.79 0.0993 0.9003 0.4239 0.92
Cylinder 2 Ulna 19.36 -30.80 9.02 0.8531 0.0183 -0.5108 1.90
Cylinder 3 Ulna 19.36 -30.80 9.02 -0.8531 -0.0183 0.5108 1.50
Cylinder 4 Radius 19.79 -43.87 3.87 0.0186 0.9767 0.2136 0.90
Cylinder 5 Ulna 20.77 -51.88 3.88 -0.1157 0.9547 0.274 0.70
Cylinder 6 Radius 19.28 -40.60 3.49 -0.1481 -0.8839 -0.4436 0.71
All values are in cm
Med Biol Eng Comput
123

2.7 Evaluation of the DSEM
To qualitatively validate the three different simulation
architectures available in the DSEM (i.e., IDO, IFDO, and
IFDOC), we compared EMG signals with model-predicted
muscle forces. To this end, one patient (male, 64 years,
163 cm, 85 kg) with shoulder hemi-arthroplasty was
measured after giving informed consent. Measurements
included the recording of pose and EMG. The subject was
asked to perform the standard shoulder dynamic tasks
including forward flexion motions up to maximum possible
arm elevation angle. The speed of movement was about
0.1 Hz.
For motion recordings, marker clusters on bony seg-
ments, including thorax, scapula, upper arm, and forearm,
were measured using four Optotrak camera bars (Northern
Digital Inc., Canada, nominal accuracy 0.3 mm) at a
sampling frequency of 50 Hz. Considering the limited
Inverse skeleton
dynamics
Mi Muscle load sharing
(inverse optimization)
Fm,i, ,i i i
F
x
a
m
,i
Muscle force-
length relationshipLm,i
Inverse skeleton
dynamics
Mi Muscle load sharing
(inverse optimization)
Inverse muscle
model
Fm,i
ei-1 , ai , Lce,i
, ,i i i
Forward muscle
model
F m
in
,i
F m
a
x ,
i
Lm,i
ei-2 , ai-1 , Lce,i-1
u
i-2
=
0
u
i-2
=
1
Inverse skeleton
dynamics
Mi Muscle load sharing
(inverse optimization)
Inverse muscle
model
Fm,i
i
Forward muscle
model
F m
in
,
i
F m
a
x,
i
Lm,i ai-1 , Lce,i-1
u
i-2
=
0
u
i-2
=
1
Mc,i
+
-
ai , Lce,i
Forward
musculoskeletal
model
,
_
_
1ci
ei-2
i
ai-1 , Lce,i-1
ei-1
i
,1ci
,
,
c i
c i
Feedback
controller
Mci+1
(a)
(b)
(c)
Fig. 2 Schematic of the a IDO,
b IFDO, and c IFDOC models
Med Biol Eng Comput
123

range of motion of the patient we used an acromion sensor
[56] for scapular motion tracking. Local coordinate sys-
tems of the segments were defined according the ISB
standardization protocol [57].
EMG signals of 12 superficial muscles were measured
using Ambu N-00-S ECG bipolar surface EMG electrodes
and recorded by a 16-channels Porti system (TMS Inter-
national, Enschede, The Netherlands) at the sampling fre-
quency of 1000 Hz. The SENIAM recommendations [58]
were followed for the EMG sensor positioning. We visu-
ally checked the measured signals for possible crosstalk.
The measured muscles included the trapezius ascendens,
transversum, and descendens, the infraspinatus, the deltoid
anterior, medialis, and posterior, the pectoralis major cla-
vicular and thoracic parts, the biceps short head, the triceps
medialis, and the brachioradialis. To determine the maxi-
mum EMG values, maximum voluntary contractions
(MVCs) were also performed.
From the segment poses, joint angles were calculated
based upon the ISB-standard and were used to run the three
models. The ECF was used in all models to solve the muscle
load sharing problem in the inverse dynamic optimization.
To guarantee the stability of the forward dynamics simula-
tions in the IFDOC, the Adams–Moulton algorithm with an
integration time-step of 0.005 s was used. The individual
muscle forces as well as the glenohumeral joint reaction
forces were estimated as the outputs of the model. The cal-
culated muscle forces were normalized (relative muscle
force) to the maximum muscle force (Eq. 3).
Measured EMGs were high-pass filtered, rectified, and
subsequently low-pass filtered. For high- and low-pass
filtering, second-order Butterworth filter with cut-off fre-
quencies of, respectively, 25 and 2 Hz were used. For each
muscle, the measured EMG was normalized with respect to
the maximum value found for that muscle during MVCs.
To evaluate the model, the time series of the relative
forces and normalized EMG were compared. For each
muscle, the comparison was carried out for the muscle
element which was the closest to the position of the EMG
electrodes on the subject body. Since the time series of
forces and EMG were compared, we used the bivariate
two-tailed Pearson correlation coefficient (R) as indicator
of goodness of fit. Moreover, the resultant glenohumeral
joint reaction force was compared between different
modeling architectures.
3 Results
The simulation times for IDO, IFDO, and IFDO were 43.3,
68.2, and 423.12 s, respectively.
3.1 EMG–force comparison
For most conditions the predicted forces followed the
pattern of the EMG signals (for IDO and IFDO mean
R * 0.71, for IFDOC mean R * 0.60). The IDO and
IFDO showed almost the same but in few cases (e.g., tra-
pezius ascendens and pectoralis major clavicular) different
results from IFDOC. In a few cases false-negatives were
found in which no muscle force was calculated by the
model while EMG showed activity for that muscle (Fig. 3).
For IDO and IFDO models, the false negatives occurred for
trapezius descendens, deltoid posterior, and pectoralis
major thoracic muscles. For IFDOC model, the false neg-
atives were related to pectoralis major thoracic part. Except
for false-negatives, three models followed the pattern of
EMG signals. A very high correlation (R * 0.97) was
found between the IDO and IFDO estimated force–time
curves with the EMG signal of trapezius transversum,
deltoid anterior, and triceps medialis. The correlation
between the IFDOC predictions and the EMG was rela-
tively high (R [ 0.80) for deltoid anterior, pectoralis major
clavicular part, and triceps medialis.
The results of comparing the estimated muscle force–
time curves to measured EMG signals during forward
flexion motion in the current study are comparable to the
ones in the study by van der Helm [16]. In the study by van
der Helm, the comparison was performed for the inverse
dynamics model and using the DSM original anatomical
dataset and the SCF for IDO. In both the studies the false-
negatives occurred for deltoid posterior and pectoralis
major thoracic part (for humeral elevation B100). A very
similar pattern was observed in two studies for the pre-
dicted muscle forces of infraspinatus and deltoid anterior
muscles.
3.2 Glenohumeral joint reaction force comparison
The IFDOC predicted considerably higher (*24%) reac-
tion force in the glenohumeral joint at the peak elevation
angle as compared with the IDO and IFDO (Fig. 4).
Table 3 The position ([Px Py Pz]) and orientation ([dx dy dz]) of the functional axes of rotation of the elbow
Rotation axis Bony structure Px Py Pz dx dy dz
Flexion–extension Humerus-ulna 19.36 -30.80 9.02 0.8531 0.0183 -0.5108
Pronation–supination Ulna-radius 20.11 -32.27 6.87 -0.0604 0.9874 0.1465
All values are in cm
Med Biol Eng Comput
123

4 Discussion
This study aimed to, first, report about all developments of
the DSEM since its early introduction and, second, to
compare the force–time curves with the EMG signals to
qualitatively evaluate the three simulation architectures
available in the DSEM.
The qualitative validation was carried out for the generic
model for one subject performing a typical dynamic
shoulder task (i.e., forward flexion). The model was not
scaled to subject’s geometry. Main reasons for this are the
difficulties related to scaling in general, but also the choice
to present the model as currently and up till now mostly
used, which is not scaled, but with scaled (optimized)
kinematics. Results (Fig. 3) showed a relatively good
agreement between the model-predicted normalized forces
and measured EMG of individual muscles, although a few
cases of false-negatives were observed. The IDO and IFDO
showed very similar patterns but somewhat different pat-
terns from IFDOC. The predicted glenohumeral joint
reaction force by the IFDOC was also higher in comparison
to the other two models (Fig. 4).
As discussed earlier, the original DSM was previously
validated by comparing the estimated muscle force–time
curves to measured EMG signals [16]. That comparison
was performed using the DSM original anatomical dataset
and the SCF for IDO, while the specific muscle force–
length relationship was not included in the model. In this
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
trap. ascend.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
trap. transv.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
trap. descen.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
delt. ant.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
delt. med.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
delt. post.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
infrasp.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
biceps
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
pect. maj. thor.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
pect. maj. clav.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
triceps med.
0 50 100
0
0.5
Arm elevation
re
la
tiv
e
fo
rc
e
brachrad.
IDO IFDO IFDOC
Fig. 3 The relative muscle
forces versus normalized EMG
for 12 muscles (muscle parts)
during forward flexion motion.
The simulations from three
modeling architectures (IDO,
IFDO, and IFDOC) are
compared. The force and EMG
are plotted against arm elevation
angles (in degrees). trap ascen
trapezius ascendens, trap transv
trapezius transversum, trap
descen trapezius descendens,
infrasp infraspinatus, delt ant
deltoid anterior, delt med
deltoid medialis, delt post
deltoid posterior, pect maj thor
pectoralis major thoracic, pect
maj clav pectoralis major
clavicular, biceps biceps short
head, triceps med triceps
medialis, brachrad
brachioradialis
Med Biol Eng Comput
123

study, the new anatomical dataset and optimization crite-
rion (the ECF) were used and validation was performed for
the three modeling architectures (IDO, IFDO, and IFDOC).
In the modified IDO model used for the current study, the
muscle force–length relationships were also considered in
the inverse optimization.
Praagman et al. [25] compared the model predicted
muscle forces to measured muscle oxygen consumption.
They used comparable elbow isometric contractions and
applied both the SCF and the ECF. They concluded that the
ECF led to fewer false-negatives and a higher correlation
between predicted muscle forces and measured oxygen
consumption. In a more recent study [50], it was shown
that comparing to the SCF using the ECF makes a better
consistency between the experimentally measured and
model estimated so-called principal action. The results of
these studies suggest that the ECF would be the preferred
optimization criterion for the DSEM. Therefore, in the
current study, the ECF was used in the process of model
evaluation.
The IDO and IFDO predicted similar forces. Therefore,
the effects of considering the muscle force–velocity rela-
tionship in case of low-speed motions (in our case
*0.1 Hz) are not remarkable. One would expect the
muscle force–velocity relationship to be more of influence
during high-speed shoulder movements like throwing a ball
in baseball. However, the IFDOC predictions of the rela-
tive muscle forces as well as the glenohumeral joint
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
trap. ascend.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
trap. transv.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
trap. descen.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
delt. ant.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
delt. med.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
delt. post.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
infrasp.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
biceps
0 50 10
0
0.5
Arm elevation
n
o
rm
.
EM
G
pect. maj. thor.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
pect. maj. clav.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
triceps med.
0 50 100
0
0.5
Arm elevation
n
o
rm
.
EM
G
brachrad.
Fig. 3 continued
Med Biol Eng Comput
123

reaction force were somewhat different from those of the
other two models. Such differences can be due to one of (or
a combination of) the following reasons:
The differences between models can relate to the for-
ward-dynamics optimization and/or the feedback control-
ler which may lead to calculation of noticeably different
neural inputs and/or large additional moment (i.e., Mc,
Fig. 2c). For the motion simulated in this study, the
additional moment at the peak elevation angle is calcu-
lated to be *22% of the net shoulder joint moment (i.e.,
M, Fig. 2c). The additional moment which is added to the
net joint moment in the IFDOC model, can explain the
difference (*24%) between the glenohumeral joint
reaction force estimations from the IDO and IFDOC.
Considerably large additional moment can be related to
either large differences between the optimized angles as
input to IFDOC and the resulting angles from the forward
simulation or to large feedback gains in the feedback
controller.
Another reason could be the differences between the
optimized angles as input to IDO (and IFDO) and IFDOC.
For the motion simulated in this study, all optimized angles
in different versions of the models were almost the same
except the optimized clavicular axial rotation which dif-
fered up to *20 between IDO and IFDOC. Such differ-
ence may explain the different relative muscle force
estimations from IDO and IFDOC in the case of muscles
with an origin/insertion on the clavicle such as pectoralis
major clavicular part. The difference between the opti-
mized angles in IDO and IFDOC may relate to the different
choices of generalized coordinates (for detailed description
about the generalized coordinates see Ref. [9]). Although
different versions of the model have the same number of
generalized coordinates (equal to the total number of
DOF), the choice of these coordinates to describe the
orientation of the scapula and clavicle in IDO is different
form that in IFDOC. In IDO and IFDO, the position
coordinates of the bony landmarks including the y- and
z-coordinates of the most dorsal point on the acromiocla-
vicular (AC) joint and the x-coordinate of the trigonum
spinae (TS) are chosen for generalized coordinates. In IF-
DOC, these coordinates include the rotations around the y-
axis (pro/retraction) and z-axis (elevation/depression) in
the sternoclavicular joint and the rotation around the y-axis
(pro/retraction) in the AC joint.
In a recent study and in an attempt to quantitatively
validate the DSEM [59], the glenohumeral-joint reaction
forces estimated by the IDO model were compared to those
measured by the instrumented shoulder implant [60]. The
results of that study showed that the generic IDO model
generally underestimates the glenohumeral-joint reaction
forces during standard dynamic tasks such as abduction and
forward flexion. According to the results of the current
study, the IFDOC predicts higher glenohumeral joint
reaction forces during dynamic motions like forward flex-
ion compared with the IDO. One may, therefore, conclude
that the IFDOC can potentially be a better candidate for
modeling dynamic tasks. Less number of false-negatives
predicted by the IFDOC compared with the IDO and IFDO
(Fig. 3) supports this line of reasoning. However, a rigor-
ous model validation is still required for a decisive con-
clusion. To this end, the model needs to be modified to
have exactly the same input optimized angles as in the
IDO, IFDO, and IFDOC. Scaling the model to the subject-
specific geometry can also minimize the effect of differ-
ences between the optimized angles (as inputs to IFDOC)
and the resulting (corrected) angles from the forward
simulation on the magnitude of additional moment.
Moreover, the feedback gains in the feedback controller of
the IFDOC need to be optimized. This should be done in a
separate study in which considerable numbers of subjects
are present and both the measured muscle activities as well
as in vivo measured joint reaction forces are used to tune
these parameters.
When the DSEM is compared to other existing upper
extremity models some differences can be discerned:
In contrast to the DSEM in which the motions of the
scapula and clavicle are used as inputs, the Swedish
Shoulder Model (SSM) and the SIMM model use the
‘‘shoulder rhythm’’ as input for scapular motions. While
this is practical since kinematic data collection is highly
simplified, the downside of that option is the limitation of
their use to applications where scapular motion is not
disturbed.
For optimization most models use the quadratic cost
function, although the SSM also uses the so-called soft-
saturation criterion [61]. The AnyBody model uses a min/
max criterion [62] as cost function.
0 20 40 60 80 100
100
300
500
Arm elevation (deg.)
R
es
ul
ta
nt
G
HJ
RF
(N
)
IDO
IFDO
IFDOC
Fig. 4 The predicted resultant reaction force in the glenohumeral
joint by the IDO, IFDO, and IFDOC models versus arm elevation
angle during forward flexion motion. GHJRF glenohumeral joint
reaction force
Med Biol Eng Comput
123

For the SSM, the predicted forces and the normalized
EMG patterns of four muscles of the shoulder have been
compared [63, 64]. Most importantly, results showed sig-
nificant differences above 60–90 degrees humeral elevation
during abduction.
The Newcastle shoulder model (NSM) is based to a
large extent on the same data as the DSM, but also includes
data from [65]. Although the muscle force predictions from
NSM have been compared with those of DSM and SSM,
there is no individual report on validation of the model.
The AnyBody shoulder and elbow model uses the ori-
ginal anatomical dataset of the DSEM but its structure is
slightly different: the scapulothoracic-gliding plane and
wrapping contours for the deltoid muscle have differently
been modeled. The validation of the model was performed
for wheelchair propulsion [66].
The recently developed model by Blana et al. [67] uses
the new DSEM anatomical dataset and optimization cri-
terion, but the SIMM algorithms for calculating muscle
wrapping paths rather than those of SPACAR. The model
can be used in both inverse and forward dynamics analyses,
and was evaluated by force-EMG comparison for standard
dynamic and activity daily living tasks using a similar
method to that described here.
A sophisticated musculoskeletal model of the entire
shoulder and elbow was represented and qualitatively
validated. The developed model has capability of inverse-
dynamics, forward-dynamics, and combined inverse-for-
ward dynamics analysis and has potential to be recruited in
different clinical and biomechanical applications. For more
realistic predictions, it is recommended to use the new
anatomical dataset and the new developed muscle load
sharing cost function.
Open Access This article is distributed under the terms of the
Creative Commons Attribution Noncommercial License which per-
mits any noncommercial use, distribution, and reproduction in any
medium, provided the original author(s) and source are credited.
References
1. Nikooyan AA, Zadpoor AA (2011) An improved cost function
for modeling of muscle activity during running. J Biomech 44(5):
984–987
2. Zadpoor AA, Nikooyan AA (2006) A mechanical model to
determine the influence of masses and mass distribution on the
impact force during running—a discussion. J Biomech 39(2):
388–390
3. Zadpoor AA, Nikooyan AA (2010) Modeling muscle activity to
study the effects of footwear on the impact forces and vibrations
of the human body during running. J Biomech 43(2):186–193
4. Zadpoor AA, Nikooyan AA, Arshi AR (2007) A model-based
parametric study of impact force during running. J Biomech 40(9):
2012–2021
5. Nikooyan AA, Zadpoor AA. Mass-spring-damper modeling of
the human body to study running and hopping: an overview. Proc
Inst Mech Eng Part H J Eng Med (in press)
6. Damsgaard M, Rasmussen J, Christensen ST et al (2006) Anal-
ysis of musculoskeletal systems in the AnyBody Modeling Sys-
tem. Simul Model Pract Theory 14(8):1100–1111
7. Glitsch U, Baumann W (1997) The three-dimensional determi-
nation of internal loads in the lower extremity. J Biomech
30(11–12):1123–1131
8. McLean SG, Su A, van den Bogert AJ (2003) Development and
validation of a 3-D model to predict knee joint loading during
dynamic movement. J Biomech Eng Trans ASME 125:864–874
9. van der Helm FCT (1994) A finite element musculoskeletal
model of the shoulder mechanism. J Biomech 27(5):551–553
10. Karlsson D, Peterson B (1992) Towards a model for force pre-
dictions in the human shoulder. J Biomech 25(2):189–199
11. Makhsous M, Hogfors C, Siemien’ski A et al (1999) Total
shoulder and relative muscle strength in the scapular plane.
J Biomech 32(11):1213–1220
12. Hogfors C, Peterson B, Sigholm G et al (1991) Biomechanical
model of the human shoulder joint—II. The shoulder rhythm.
J Biomech 24(8):699–709
13. Hogfors C, Sigholm G, Herberts P (1987) Biomechanical model
of the human shoulder—I. Elements. J Biomech 20(2):157–166
14. Charlton IW, Johnson GR (2006) A model for the prediction of
the forces at the glenohumeral joint. Proc Inst Mech Eng H
220(H8):801–812
15. Holzbaur KRS, Murray WM, Delp SL (2005) A model of the upper
extremity for simulating musculoskeletal surgery and analyzing
neuromuscular control. Ann Biomed Eng 33(6):829–840
16. van der Helm FCT (1994) Analysis of the kinematic and dynamic
behavior of the shoulder mechanism. J Biomech 27(5):527–550
17. van der Helm FCT, Veeger HEJ, Pronk GM et al (1992) Geometry
parameters for musculoskeletal modelling of the shoulder system.
J Biomech 25(2):129–144
18. van der Helm FCT, Veenbaas R (1991) Modelling the mechanical
effect of muscles with large attachment sites: application to the
shoulder mechanism. J Biomech 24(12):1151–1163
19. Veeger HEJ, Van Der Helm FCT, Van Der Woude LHV et al
(1991) Inertia and muscle contraction parameters for musculo-
skeletal modelling of the shoulder mechanism. J Biomech 24(7):
615–629
20. Veeger HEJ, Yu B, An K-N et al (1997) Parameters for modeling
the upper extremity. J Biomech 30(6):647–652
21. Klein Breteler MD, Spoor CW, Van der Helm FCT (1999)
Measuring muscle and joint geometry parameters of a shoulder
for modeling purposes. J Biomech 32(11):1191–1197
22. Minekus J (1997) Kinematical model of the human elbow. Leiden
University, Leiden
23. Chadwick EKJ, Van Der Helm FCT (2003) Combined inverse/
forward dynamic optimization of a large-scale biomechanical
model of the shoulder and elbow. In: IX international symposium
on computer simulation in biomechanics. Sydney, Australia
24. van der Helm FCT, Chadwick EKJ (2002) A forward-dynamic
shoulder and elbow model. In: 4th meeting of the international
shoulder group. Cleveland, OH
25. Praagman M, Chadwick EKJ, van der Helm FCT et al (2006) The
relationship between two different mechanical cost functions and
muscle oxygen consumption. J Biomech 39(4):758–765
26. Cutti AG, Veeger HEJ (2009) Shoulder biomechanics: today’s
consensus and tomorrow’s perspectives. Med Biol Eng Comput
47(5):463–466
27. Young LL, Papa CM, Lyon CE et al (1990) Comparison of
microscopic and laser diffraction methods for measuring sarco-
mere lengths of contracted muscle fibers of chicken pectoralis
major muscle. Poult Sci 69:1800–1802
Med Biol Eng Comput
123

28. Walker SM, Schrodt GR (1974) I segmental lengths and thin
filament periods in skeletal muscle fibres of the rhesus monkey
and human. Anat Rec 178:63–82
29. van Soest AJ, Schwab AL, Bobbert MF et al (1992) SPACAR: a
software subroutine package for simulation of the behavior of
biomechanical systems. J Biomech 25(10):1219–1226
30. Woltring HJ, Long K, Osterbauer PJ et al (1994) Instantaneous
helical axis estimation from 3-D video data in neck kinematics
for whiplash diagnostics. J Biomech 27(12):1415–1425
31. Nikooyan AA, van der Helm FCT, Westerhoff P et al (2011)
Comparison of two methods for in vivo estimation of the gle-
nohumeral joint rotation center of the patients with shoulder
hemiarthroplasty. PLoS One 6(3):e18488
32. Meskers CGM, van der Helm FCT, Rozendaal LA et al (1998) In vivo
estimation of the glenohumeral joint rotation center from scapular
bony landmarks by linear regression. J Biomech 31(1):93–96
33. Ehrig RM, Taylor WR, Duda GN et al (2006) A survey of formal
methods for determining the centre of rotation of ball joints.
J Biomech 39(15):2798–2809
34. Monnet T, Desailly E, Begon M et al (2007) Comparison of the
SCoRE and HA methods for locating in vivo the glenohumeral
joint centre. J Biomech 40(15):3487–3492
35. van der Helm FCT, Pronk GM (1995) Three-dimensional
recording and description of motions of the shoulder mechanism.
J Biomech Eng T ASME 117:27–40
36. de Groot JH (1998) The shoulder: a kinematic and dynamic
analysis of motion and loading. Delft University of Technology,
Delft
37. Woltring H (1986) A FORTRAN package for generalized, cross-
validatory spline smoothing and differentiation. Adv Eng Softw
8:104–113
38. Crowninshield RD, Brand RA (1981) A physiologically based
criterion of muscle force prediction in locomotion. J Biomech
14(11):793–801
39. An KN, Kaufman KR, Chao EYS (1989) Physiological consid-
erations of muscle force through the elbow joint. J Biomech
22(11–12):1249–1256
40. Winters JM, Stark L (1985) Analysis of fundamental human
movement patterns through the use of in-depth antagonistic
muscle models. IEEE Trans Bio-Med Eng 32:826–839
41. Happee R (1994) Inverse dynamic optimization including mus-
cular dynamics, a new simulation method applied to goal directed
movements. J Biomech 27(7):953–960
42. Happee R, Van der Helm FCT (1995) The control of shoulder
muscles during goal directed movements, an inverse dynamic
analysis. J Biomech 28(10):1179–1191
43. Winters JM (1990) Hill based muscle models: a systems engi-
neering perspective. In: Winters JM, Woo SL-Y (eds) Multiple
muscle systems: biomechanics and movement organization.
Springer, New York, pp 70–93
44. Van Drongelen S, Van Der Woude LH, Janssen TW et al (2005)
Mechanical load on the upper extremity during wheelchair
activities. Arch Phys Med Rehabil 86(6):1214–1220
45. Veeger HEJ, Rozendaal LA, van der Helm FCT (2002) Load on
the shoulder in low intensity wheelchair propulsion. Clin Bio-
mech 17(3):211–218
46. Steenbrink F, de Groot JH, Veeger HEJ et al (2009) Glenohu-
meral stability in simulated rotator cuff tears. J Biomech 42(11):
1740–1745
47. Magermans DJ, Chadwick EKJ, Veeger HEJ et al (2004) Effec-
tiveness of tendon transfers for massive rotator cuff tears: a
simulation study. Clin Biomech 19(2):116–122
48. Magermans DJ, Chadwick EKJ, Veeger HEJ et al (2004) Bio-
mechanical analysis of tendon transfers for massive rotator cuff
tears. Clin Biomech 19(4):350–357
49. Praagman M, Stokdijk M, Veeger HEJ et al (2000) Predicting
mechanical load of the glenohumeral joint, using net joint
moments. Clin Biomech 15(5):315–321
50. Steenbrink F, Meskers CGM, van Vliet B et al (2009) Arm load
magnitude affects selective shoulder muscle activation. Med Biol
Eng Comput 47(5):565–572
51. Chadwick EKJ, Van Noort A, Van Der Helm FCT (2004) Bio-
mechanical analysis of scapular neck malunion—a simulation
study. Clin Biomech 19(9):906–912
52. Van Drongelen S, Van Der Woude LH, Janssen TW et al (2005)
Glenohumeral contact forces and muscle forces evaluated in
wheelchair-related activities of daily living in able-bodied sub-
jects versus subjects with paraplegia and tetraplegia. Arch Phys
Med Rehabil 86(7):1434–1440
53. Van Drongelen S, Van Der Woude LHV, Janssen TWJ et al
(2006) Glenohumeral joint loading in tetraplegia during weight
relief lifting: a simulation study. Clin Biomech 21(2):128–137
54. Nikooyan AA, Veeger HEJ, Westerhoff P et al. An EMG-driven
musculoskeletal model of the shoulder. Hum Mov Sci (in press)
55. Suarez DR, Valstar ER, Linden JC et al (2009) Effect of rotator
cuff dysfunction on the initial mechanical stability of cementless
glenoid components. Med Biol Eng Comput 47(5 SPEC. ISS):
507–514
56. Karduna AR, McClure PW, Michener LA et al (2001) Dynamic
measurements of three-dimensional scapular kinematics: a vali-
dation study. J Biomech Eng Trans ASME 123:184–190
57. Wu G, van der Helm FCT, Veeger HEJ et al (2005) ISB rec-
ommendation on definitions of joint coordinate systems of vari-
ous joints for the reporting of human joint motion—part ii:
shoulder, elbow, wrist and hand. J Biomech 38(5):981–992
58. Merletti R, Wallinga W, Hermens H et al (1999) Guidelines for
reporting SEMG data. Roessingh Research and Development
59. Nikooyan AA, Veeger HEJ, Westerhoff P et al (2010) Validation
of the Delft Shoulder and Elbow Model using in vivo glenohu-
meral joint contact forces. J Biomech 43(15):3007–3014
60. Westerhoff P, Graichen F, Bender A et al (2009) An instrumented
implant for in vivo measurement of contact forces and contact
moments in the shoulder joint. Med Eng Phys 31(2):207–213
61. Siemienski A (1992) Soft saturation—an idea for load sharing
between muscles. Application to the study of human locomotion.
In: Cappozzo A, Marchetti M, Tosi V (eds) Biolocomotion: a
century of research using moving pictures. Promograph, Rome,
pp 293–303
62. Rasmussen J, Damsgaard M, Voigt M (2001) Muscle recruitment
by the min/max criterion—a comparative numerical study.
J Biomech 34(3):409–415
63. Jarvholm U, Palmerud G, Herberts P et al (1989) Intramuscular
pressure and electromyography in the supraspinatus muscle at
shoulder abduction. Clin Orthop Relat Res 245:102–109
64. Jarvholm U, Palmerud G, Karlsson D et al (1991) Intramuscular
pressure and electromyography in four shoulder muscle. J Orthop
Res 9(4):609–619
65. Johnson G, Spalding D, Nowitzke A et al (1996) Modelling the
muscles of the scapula morphometric and coordinate data and
functional implications. J Biomech 29(8):1039–1051
66. Dubowsky SR, Rasmussen J, Sisto SA et al (2008) Validation of
a musculoskeletal model of wheelchair propulsion and its appli-
cation to minimizing shoulder joint forces. J Biomech 41(14):
2981–2988
67. Blana D, Hincapie JG, Chadwick EK et al (2008) A musculo-
skeletal model of the upper extremity for use in the development
of neuroprosthetic systems. J Biomech 41(8):1714–1721
Med Biol Eng Comput
123