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 needed. Such models should be complex enough to realistically 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]; consequently, 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 optimization 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 literature. 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 application of an energy-based load sharing cost function enabled by the addition of the muscle architecture parameters. 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 biomechanical (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 modeling using geometrical forms. Muscle architecture parameters including tendon length, physiological cross-sectional area (PCSA), pennation angle, and the fiber length (l f) were measured. The sarcomere lengths (l s) were also measured using a laser-diffraction technique [27]. Assuming an optimal sarcomere length of 2.7 μm [28], the optimal fiber length (l opt) for a muscle was calculated as:

$$ l_{\text{opt}} = 2.7\frac{{l_{\text{f}} }}{{l_{\text{s}} }} $$
(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.

Fig. 1
figure 1

Palpable bony landmarks on the humerus, ulna, and radius

Full size image
Table 1 Positions of palpable bony landmarks (BL) on the humerus, ulna, and radius in the global coordinate system defined in Ref. [21]
Full size table
Table 2 Bony contours. [d x d y d z ] is the direction of the central axis of the cylinder; [P x P y P z ] is the coordinate of an arbitrary point on the central axis of the cylinder; R is the radius
Full size table
Table 3 The position ([P x P y P z ]) and orientation ([d x d y d z ]) of the functional axes of rotation of the elbow
Full size table

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 implementation 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 = W_{1} \left( {\left( {{\text{d}}C_{x} } \right)^{2} + \left( {{\text{d}}C_{y} } \right)^{2} + \left( {{\text{d}}C_{z} } \right)^{2} } \right) + W_{2} \left( {\left( {{\text{d}}S_{x} } \right)^{2} + \left( {{\text{d}}S_{y} } \right)^{2} + \left( {{\text{d}}S_{z} } \right)^{2} } \right) $$
(2)

where dC x and dS x are the differences between the measured 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. W 1 and W 2 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 ligament lengths, and muscle power.

Fig. 2
figure 2

Schematic of the a IDO, b IFDO, and c IFDOC models

Full size image

The filtering and differentiation routines of Woltring (the GVC method) [37] have been implemented to calculate 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 minimization 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 (σ max), PCSA, and sarcomere length (l sm):

$$ F_{\max } \left( {l_{\text{sm}} } \right) = f\left( {l_{\text{sm}} } \right){\text{PCSA}}_{\text{m}} \sigma_{\max } $$
(3)

where σ max is taken as 100 N/cm2 [39]. f(l sm) is the normalized muscle force–length relationship defined as a 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 forward 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 optimization 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 (F max,i ) and minimum (F min,i ) permissible values of the muscle forces estimated by a forward muscle model with use of the muscle states of the previous time-step (e i−2, a i,−1, L ce,i−1). At the same time-step (i), an inverse muscle model is used to estimate the neural inputs (e i−1, a i , L ce,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 L ce are estimated through a steady state equilibrium condition. During the next time steps, the states are updated in a dynamic optimization procedure. The values of e, a, and L ce are updated iteratively, while u is estimated analytically.

2.5 The IFDOC

Due to discretization errors and possibly an unstable system, 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 (M c) using the errors in angle and angular velocity. 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 propulsion [44, 45], rotator cuff tears [46], tendon transfers [47, 48], loads on the arm [49, 50], rotator cuff changed following scapular neck fracture [51], effect of rotator cuff dysfunctions on wheelchair propulsion [52], weight transfers in wheelchair users [53], effect of including the neural activity in the modeling process [54], and stability of cementless glenoid prostheses [55].

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 segments, 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 range of motion of the patient we used an acromion sensor [56] for scapular motion tracking. Local coordinate systems 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 International, Enschede, The Netherlands) at the sampling frequency of 1000 Hz. The SENIAM recommendations [58] were followed for the EMG sensor positioning. We visually 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 clavicular and thoracic parts, the biceps short head, the triceps medialis, and the brachioradialis. To determine the maximum 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 simulations 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 calculated 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 frequencies 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., trapezius 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 negatives 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 relatively high (R > 0.80) for deltoid anterior, pectoralis major clavicular part, and triceps medialis.

Fig. 3
figure 3figure 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

Full size image

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 ≤100°). A very similar pattern was observed in two studies for the predicted muscle forces of infraspinatus and deltoid anterior muscles.

3.2 Glenohumeral joint reaction force comparison

The IFDOC predicted considerably higher (~24%) reaction force in the glenohumeral joint at the peak elevation angle as compared with the IDO and IFDO (Fig. 4).

Fig. 4
figure 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

Full size image

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 patterns 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 study, the new anatomical dataset and optimization criterion (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 relationship 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 relative muscle forces as well as the glenohumeral joint 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 forward-dynamics optimization and/or the feedback controller 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 calculated 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 differed up to ~20° between IDO and IFDOC. Such difference 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 optimized 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 acromioclavicular (AC) joint and the x-coordinate of the trigonum spinae (TS) are chosen for generalized coordinates. In IFDOC, 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 flexion 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 rigorous model validation is still required for a decisive conclusion. 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 differences 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.

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 significant 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 original 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 criterion, 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-forward 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.