Scaling and kinematics optimisation of the scapula and thorax in upper limb musculoskeletal models

Accurate representation of individual scapula kinematics and subject geometries is vital in musculoskeletal models applied to upper limb pathology and performance. In applying individual kinematics to a model׳s cadaveric geometry, model constraints are commonly prescriptive. These rely on thorax scaling to effectively define the scapula׳s path but do not consider the area underneath the scapula in scaling, and assume a fixed conoid ligament length. These constraints may not allow continuous solutions or close agreement with directly measured kinematics. A novel method is presented to scale the thorax based on palpated scapula landmarks. The scapula and clavicle kinematics are optimised with the constraint that the scapula medial border does not penetrate the thorax. Conoid ligament length is not used as a constraint. This method is simulated in the UK National Shoulder Model and compared to four other methods, including the standard technique, during three pull-up techniques (n=11). These are high-performance activities covering a large range of motion. Model solutions without substantial jumps in the joint kinematics data were improved from 23% of trials with the standard method, to 100% of trials with the new method. Agreement with measured kinematics was significantly improved (more than 10° closer at p<0.001) when compared to standard methods. The removal of the conoid ligament constraint and the novel thorax scaling correction factor were shown to be key. Separation of the medial border of the scapula from the thorax was large, although this may be physiologically correct due to the high loads and high arm elevation angles.


Introduction
Scaling of musculoskeletal models is important to accurately represent inter-segmental joint moments. Equal scaling of all segments improves accuracy of force predictions in dynamic tasks, but has little effect in static-loaded tasks (Nikooyan et al., 2010): where external loads are important. Once a model's geometry has been scaled and measured kinematics supplied as inputs, it is necessary to test that the model remains physiological; solid structures must not overlap and close agreement with measured kinematics is important: within the measurement errors of dynamic measurement techniques (e.g. Shaheen et al., 2011a).
Commonly, constrained optimisation of measured kinematics is used (e.g. Dickerson et al., 2007). These methods constrain the medial border of the scapula to a fixed offset from the thorax surface, the conoid ligament to a fixed length and do not scale the generic model (Nikooyan et al., 2011). These constraints define the scapula's path based on thorax scaling methods that do not currently consider the area underneath the scapula. This method does not regularly provide continuous solutions (Martelli et al., 2008), and may result in large differences between measured and optimised angles (Bolsterlee et al., 2012).
According to open-MRI (Izadpanah et al., 2012) and CT (Seo et al., 2012) studies, the conoid ligament is not a fixed length during arm abduction. The conoid ligament's length and attachments vary by up to 25% between subjects (Harris et al., 2001;Takase, 2010), and may scale with clavicle length (Rios et al., 2007). It is therefore theorised that relaxation of the standard kinematics optimisation constraints and redefinition of thorax scaling will allow more robust scalable simulations that maintain close agreement with measured kinematics. 0.9-1.6 and 0.9-1.8 s for the front, wide and reverse pull-ups respectively. Legs were kept in a fixed position with posterior-facing heels (Fig. 1). Pull-ups dynamically cover large ranges of motion: 23-1261 humerothoracic elevation and À 56À 101 humerothoracic axial rotation.
The UK National Shoulder Model (UKNSM; Charlton and Johnson, 2006) is the musculoskeletal model used for the simulations: accepting marker positions to calculate local ISB-recommended coordinate systems and Euler rotations (Wu et al., 2002) to then apply inverse dynamics and static optimisation. A nine-camera Vicon motion capture system (200 Hz) recorded marker positions.

Scaling
Clavicle and scapula segments were homogeneously scaled based on relative segment lengths between model and subject. Clavicle scaling used sternoclavicular (SC) to acromioclavicular (AC) joint distance and scapula scaling AC to scapula inferior angle (AI) distance.
An ellipse represented the scapulothoracic gliding-plane (STGP): described by Charlton (2003). In this study a novel technique is described to improve the thorax STGP scaling. Firstly, a homogeneous scaling factor was defined for the STGP ellipse based on subject height (Eq. (1)). Static measurements of bony landmarks were taken: (1) at rest with the arms by the side, (2) arms horizontal at 451 to the coronal plane and (3) arms at subject's maximal elevation. The following landmarks were measured in these three static poses: Acromial angle, trigonum spinae (TS) and AI (measured using the scapula Palpator; Shaheen et al., 2011b). SC, jugular notch, xiphoid process, 7th cervical vertebra and 8th thoracic vertebra (measured using skin-fixed markers). AC joint measured relative to the scapula Tracker, digitised at 901 elevation and 451 horizontal abduction (Prinold et al., 2011).
where T h is the initial homogeneous scaling factor for the thorax ellipse. The cadaver's height, used for the original UKNSM, was 1.80 m.
An optimisation procedure was then used to define the thorax ellipse size and position. The constraint applied to the optimiser was that the scapula medial border (TS and AI) could not fall within the STGP ellipse. The optimiser objective function kept the non-homogeneous scaling factors as close as possible to the original homogeneous scaling factor in a least squares sense (Eq. (2)). The new nonhomogeneous scaling factors (T x,y,z ) were then applied to the size and centre position of the STGP ellipse and the thorax-body ellipse for the pull-up simulations.
min ½ðTx À T h Þ 2 þðTy À T h Þ 2 þðTz À T h Þ 2 such that : Here T x , T y , T z are the non-homogeneous thorax scaling parameters being optimised; i indicates each of the 3 arm positions; AI x,i , AI y,i , AI z,i are the position of AI in trial i; TS x,i , TS y,i , TS z,i are the position of TS in trial i; M x , M y , M z are the centre of the thorax ellipse; A x , A y , A z are the vectors corresponding to axes lengths. The ellipse thickness of 10 mm is included here through addition to the axes lengths in each dimension. Eq.
(2) is the optimisation and constraints used to find the thorax ellipses scaling factors.

Kinematics optimisation
The fixed closed chain (FCC) method optimises the scapula and clavicle kinematics to reduce the least squares difference to the measured scapula and clavicle kinematics. Meanwhile the scapula medial border is constrained to stay 10 mm from the STGP ellipse and the conoid ligament constrained to a fixed length (Eq. (3); van der Helm, 1994).
The partially closed chain (PCC) method also optimises the scapula and clavicle kinematics to reduce the least squares difference to the measured kinematics. However, the scapula medial border is only constrained to not penetrate the STGP ellipse and the conoid ligament length is not used as a constraint: for PCC method such that : or for FCC method such that : The scapula angles; Sx, Sy, Sz are the scapulothoracic Euler angles being optimised. The clavicle angles; Cx, Cy, Cz are the claviculothoracic Euler angles being optimised. L con is the length of the conoid ligament; computed based on the optimised scapulothoracic and clavicle angles. AI x , AI y , AI z are the coordinates of AI and TS x , TS y , TS z are the coordinates of TS; both of which are computed based on the optimised scapulothoracic and clavicle angles.
The constant values are the measured scapulothoracic Euler rotations (Sx m , Sy m , Sz m ), the claviculothoracic Euler angles (Cx m , Cy m , Cz m ), the rest length of the conoid ligament (L con rest ), the coordinates of the centre of the scapula gliding-plane ellipse (STGP) after the described thorax scaling procedure has been applied (M x , M y , M z ) and the axes lengths of the STGP ellipse after the described thorax scaling procedure has been applied and including an ellipse thickness of 10 mm (A x , A y , A z ).
Eq. (3) is the kinematics optimisation objective function with constraints for the described FCC and PCC methods.
Variations of the PCC and FCC methods are tested to understand the effect of each element of the optimisation (Table 1). Clavicle axial rotation was calculated with regression equations in the PCC-based methods (Charlton and Johnson, 2006) and with a minimisation technique in the FCC-based methods (van der Helm, 1994).
A sequential quadratic programming algorithm implemented in the fmincon function of MATLAB (v2012b) was used to perform the kinematics optimisation. A first-order-optimality measure based on Karush-Kuhn-Tucker (KKT) conditions was calculated, and used as the optimiser's stopping criterion. In frame one of each trial the initial joint angles given to the optimiser were randomly varied five hundred times within physiologically feasible bounds (7 22.51). The most optimal solution to the kinematics optimisation, based on the lowest function value (Eq. (3)), was then found from these five hundred solutions. The most optimal solution was then used as the starting point for the second frame optimisation. For the following frames the previous solution was the starting point to the optimiser.

Data processing
The mean differences to measured rotations were not normally distributed (Shapiro-Wilk). Friedman tests tested for differences across the five optimisationmethods in the three motions (Table 1). A post-hoc Wilcoxon signed-rank test with Holm's correction for multiple comparisons identified significant differences between pairs of optimisation methods within each motion.

Results
The FCC optimisation method leads to solutions with substantial jumps in the data for most subjects. The PCC optimisation method gives data that could be fitted to a continuous function, without substantial jumps, in all trials and subjects ( Table 2).
The effect of different optimisation parameters are compared through mean differences to measured rotations and 95% confidence intervals (C.I.) of those optimised rotations (Figs. 2 and 3). The PCC optimisation method is closest to the measured valueswith significantly smaller errors in many cases. All the other methods fall outside the scapulothoracic measurement errors (Fig. 2).
The distance between the bony landmarks on the scapula medial border and the ellipse representing the STGP are of a similar order to the landmarks' resting distances, although the TS distance can be up to twice the resting distance (Fig. 4).

Discussion
The PCC optimisation and novel thorax scaling provide continuous solutions in all pull-up trials. The FCC method allows a continuous Table 1 Details of the five kinematics optimisation methods simulated. The scaling methods used are also described. Note that the PCC and FCC methods have been described in detail in Section 2.3; the other three methods are variations of these. 'Correction factor' refers to the optimisation-based thorax scaling method described in Section 2.2. 'Homogeneous scaling based on height' refers to scaling of the STGP ellipse according to Eq. (1). 'Segment homogeneous' refers to homogeneous scaling of each segment individually, as described in Section 2.2.  solution in less than 25% of trials ( Table 2). The errors compared to measured kinematics are substantially and significantly (po0.001) reduced with the PCC optimisation and scaling correction factor, when compared to the FCC method. The errors only fall within measurement errors (at the 95% C.I.) with the PCC method (Figs. 2 and 3). The FCC method facilitates a comparison with modelling standards (Nikooyan et al., 2011), but some models are homogeneously scaled before optimisation with the FCC method (e.g., Charlton, 2003).

PCC
Homogeneous scaling of all model segments (including the thorax) with the constraints of the FCC method (Table 1) improved solution continuity and accuracy relative to the FCC method, as previously observed . However, the errors were larger than with the PCC method and significantly larger than measurement errors (Figs. 2 and 3 caption text). Additionally, 37% of trials were noncontinuous with a constant conoid length and 3% non-continuous without a constant length.  Table 1), including the 95% confidence intervals of the differences. Scapulothoracic measurement errors are shown as a dashed horizontal line. These values are the average absolute error plus three SDs (note not RMS) found for each rotation in a previous study (Prinold et al., 2011). The measurement error calculated in upward rotation was 5.61 (compared to 8.41 RMS error in abduction or forward flexion in a bone-pin study: Karduna et al., 2001), the value in internal rotation was 5.81 (compared to 3.81 RMS error in Karduna et al. (2001)) and for posterior tilt 4.91 (compared to 6.21 RMS error in Karduna et al. (2001)). n indicates po 0.05, nn p o 0.01, nnn p o0.0001. All trials were included in the statistical analysis. The abbreviations used for the optimisation methods are described in Table 1. Additional simulations using homogeneous scaling of all segments based on segment length (including the thorax) and the FCC method constraints, found average differences across the three pull-up configurations of: 9. The PCC optimisation shows large separations of the scapula medial border from the thorax (maximum mean 5.3 cm; Fig. 4). The presented wide pull-up is representative of the other pull-up's pattern and magnitude of separation. Separation primarily occurs at TS, with AI generally staying below the rest separation distance. This pattern is likely to be physiological due to the loading condition; the scapula is levered backwards (Fig. 1), increasing separation of TS and pressing AI against the thorax. Also, the large hand forces will cause the arms to translate superiorly, potentially lifting TS above the thorax surface. Cadaver values for rest separation (TS: 3.72 cm and AI: 2.63 cm; Klein Breteler et al., 1999) are similar to the mean values found here (Fig. 4). Although reasonable explanations for the large separations, the true values are not known in highly loaded overhead activities like pull-ups. Further work could utilise bi-planar fluoroscopy (Bey et al., 2006) or open-MRI (Sahara et al., 2007) in quantifying these separations. A maximum separation bounding may be sensible, particularly using a weighted optimisation approach (Bolsterlee et al., 2013). Fig. 3. Average clavicle rotation differences to measured kinematics for (a) front, (b) wide and (c) reverse configuration pull-ups using various scaling and kinematics optimisation strategies (described in Table 1), including the 95% confidence intervals of the kinematics. The average RMS clavicle measurement errors found for the Scapula Tracker in a bone pin validation study (Karduna et al., 2001) are included as a dashed horizontal line. These values are 1.21 and 1.61 for clavicle protraction and upward rotation respectively. n indicates p o0.05, nn p o 0.01, nnn p o 0.0001. Note all trials were included in the statistical analysis. The abbreviations used for the optimisation methods are described in Table 1. Additional simulations using homogeneous scaling of all segments based on segment length (including the thorax) and the FCC method constraints, found average differences across the three pull-up configurations of: 10.9 7 3.31 (CT internal) and 9.0 7 2.81 (CT upward). The same values, with homogeneous scaling of all segments, but without a constrained conoid ligament length were smaller: 6.1 7 2.21 (CT internal) and 7.9 7 3.21 (CT upward).
A study using an optimisation approach to restrict the scapula separation  found smaller separation ( À 1 to 1.5 cm) during a quasi-static movement, where separation is expected to be small compared to pull-ups. The errors from measured scapula rotations are of a similar order, although the range is smaller with the PCC methods here. The range of errors from measured scapula rotations with the FCC method is similar between the two studies.
Future work could utilise non-homogeneous bone morphing techniques (Kaptein and van der Helm, 2004;Yang et al., 2008). The large variation in the externally palpatable landmarks makes scaling difficult (MacGillivray et al., 1998;Pappas et al., 2006). Recent work used coracoid process palpation to improve this with a significant improvement in biceps origin scaling (Bolsterlee and Zadpoor, 2013).
The kinematics optimisation objective function is a convex programming problem, but with the constraints taken into account the feasible regions of both methods are non-convex. Testing five hundred variations of the first frame joint angles increased the likelihood of finding a global minimum within the physiologically feasible region. The mean first-order optimality measure of 3.8 Â 10 À 5 for the FCC method shows that KKT points are generally reached, and a jump in the kinematics data corresponded to a within-trial maximum firstorder-optimality measure above 1 Â 10 À 3 on only two occasions, with values of 1.1 Â 10 À 3 and 1.2 Â 10 À 3 observed. These values are low and imply KKT values were reached, suggesting that the jumps observed are not the results of the optimiser failing to reach a minimum.
The constraint that the conoid ligament is a fixed length is shown to have the largest effect on the solution continuity (Table 2) and the errors relative to measured rotations (Figs. 2 and 3). This constraint was originally designed to define the clavicle axial rotation, given the ligament's high stiffness and large moment arm around the clavicle long-axis (van der Helm, 1994). Model sensitivity to clavicle axial rotation is probably small compared with sensitivity to the scapula position in relation to the thorax (determined by clavicle protraction and elevation). The ligament is probably not a fixed length during arm elevation tasks (Izadpanah et al., 2012;Seo et al., 2012) and the ligament attachments may be inconsistent between subjects. Therefore, given the serious detrimental effect that the constraint has on solution continuity and the large and significant effects on accuracy, it is recommended that the constraint not be used in kinematics optimisation.
Soft-constraint of the conoid ligament length is not considered appropriate because of high model sensitivity to a parameter that cannot be accurately determined due to inter-subject anatomical variations and a lack of non-invasive clavicle axial rotation measurement techniques. Inclusion in later stages of musculoskeletal modelling (e.g. load-sharing optimisation) is, however, likely to be required -probably with constant or regression-predicted strain values.
The thorax scaling correction factor had little effect on solution continuity ( Table 2). The effects on errors compared to measured rotations are generally significant (po0.01) but small (Figs. 2 and 3). Therefore the limiting factor for solution continuity is seemingly the conoid ligament constraint. Standard thorax scaling uses thorax width, depth and height for scaling (Charlton, 2003), with no consideration of the area under the scapula. RMS errors of 5.9 mm have been found between thorax surface points and a fitted ellipse . For future work other shapes may be more appropriate. Fig. 4. Separation between ellipse representing the STGP and the (a) inferior angle of the scapula and (b) trigonum spinae during the wide motion: using the PCC method of kinematics optimization. The rest position of the two landmarks is shown as the measured separation distances at a position of rest, with the arms by the side and no load on the hands. The grey area represents the positive and negative standard deviation around the mean (black line). The separations presented are the mean values during the wide pull-ups, but should be considered representative of those seen in the two other pull-up tasks.