To what extent is joint and muscle mechanics predicted by musculoskeletal models sensitive to soft tissue artefacts?

Musculoskeletal models are widely used to estimate joint kinematics, intersegmental loads, and muscle and joint contact forces during movement. These estimates can be heavily affected by the soft tissue artefact (STA) when input positional data are obtained using stereophotogrammetry, but this aspect has not yet been fully characterised for muscle and joint forces. This study aims to assess the sensitivity to the STA of three open-source musculoskeletal models, implemented in OpenSim. A baseline dataset of marker trajectories was created for each model from experimental data of one healthy volunteer. Five hundred STA realizations were then statistically generated using a marker- dependent model of the pelvis and lower limb artefact and added to the baseline data. The STA's impact on the musculoskeletal model estimates was ﬁ nally quanti ﬁ ed using a Monte Carlo analysis. The modelled STA distributions were in line with the literature. Observed output variations were comparable across the three models, and sensitivity to the STA was evident for most investigated quantities. Shape, magnitude and timing of the joint angle and moment time histories were not signi ﬁ cantly affected throughout the entire gait cycle, whereas magnitude variations were observed for muscle and joint forces. Ranges of contact force variations differed between joints, with hip variations up to 1.8 times body weight observed. Variations of more than 30% were observed for some of the muscle forces. In conclusion, musculoskeletal simulations using stereophotogrammetry may be safely run when only interested in overall output patterns. Caution should be paid when more accurate estimated values are needed. & 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Stereophotogrammetric recordings of skin-mounted marker trajectories and ground reactions are fed to musculoskeletal models (MSMs) with the aim of estimating joint angles, intersegmental loads, and muscle and joint contact forces during movement Delp et al., 2007). Unfortunately, the skin-mounted markers move over the underlying bone generating the so-called "soft tissue artefact" (STA) which makes the estimation of the instantaneous skeletal pose awkward . Normally, MSMs cope with this problem by using a multibody kinematics optimization method which embeds a least squares approach and articular constraints (Delp et al., 2007;Lu and O'Connor, 1999). The residual artefact, however, might still propagate to MSM estimates, with an effect that is still unclear, especially as far as muscle and joint forces are concerned.
Recent studies attempted to address the aforementioned problem by quantifying the sensitivity of MSMs estimates to the STA. El Habachi et al. (2015), using a global probabilistic approach and, contrary to the available evidence Peters et al., 2010), modelling the STA with the same statistics for all markers independently from their location on the body, showed that the STA may cause joint angle variations of up to 36°. The variations of muscle and joint contact forces were not investigated. Myers et al. (2015) investigated the effects of the propagation of the STA for the MSM proposed by Delp et al. (1990) through a Monte Carlo analysis and showed that the STA can induce variations in the joint angles that are 1.8 times higher than the uncertainties due to anatomical landmark identification. Myers et al. (2015) also investigated the variations Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com induced by the STA on the joint moments, and found that these were 2.3 to 4 times higher than those induced by improper positioning of skin markers on the anatomical landmarks and uncertainties in estimating the inertial parameters (i.e., mass, moment of inertia and centre of mass). The same authors also reported an impact on muscle forces, with variations due to the STA that, for gluteus medius and medial gastrocnemius, reached 50%. These effects, however, were about half of those generated by the inaccuracies affecting musculotendon parameters such as pennation angle, maximum isometric force, and tendon slack length. In this study, the STA model embedded marker-specific parameters which were also gait-phase dependent. However, STAs were constrained to have a maximal amplitude of 15 mm, in contrast with the values reported in the literature, for example the 40 mm observed at the thigh Peters et al., 2010). Finally, the effects on joint contact forces were not investigated. It therefore appears that the available information is limited to particular types of MSMs, not all of which are publically available, to a specific subset of model outputs, and to simplified STA designs. Thus, a conclusive quantification of the sensitivity of the estimates of different MSMs to a realistic and comprehensive STA representation is still lacking.
The aim of the present study was thus to investigate the sensitivity of joint angles, joint moments, and muscle and joint contact forces to a STA consistent with the best knowledge available in the literature using three different open-source MSMs and relevant tools, which are commonly used in research contexts (Arnold et al., 2010;Delp et al., 1990;Modenese et al., 2011). A probabilistic approach and published STA models were used to design a realistic set of artefact-affected marker trajectories and, through a Monte Carlo analysis, assess the statistical impact of the artefact on the outputs of the selected MSMs when studying the gait of a representative subject.

Materials and methods
A single healthy participant (male, age: 28 years, stature: 1.90 m, mass: 82 kg) was enrolled in the study after providing informed consent. Ethical approval for the study was obtained from the University Research Ethics Committee at the University of Sheffield.
Overall, twenty-eight 8mm-diameter reflective skin-markers were attached using double-sided tape to the feet (8), shanks (8), thighs (8), and pelvis (4). They were placed on the following anatomical landmarks (anatomical markers): anterior and posterior superior iliac spines (ASIS and PSIS), lateral femoral condyle (LE), tibial tuberosity (TT), lateral malleolus (LM), posterior distal aspect of the heel (HEE), forefoot (midpoint between second and third metatarsal heads; FF), heads of first and fifth metatarsals (MT1 and MT5). Furthermore, additional markers were placed in the following positions (technical markers): laterally and equidistant along the length of the thigh (TH1, TH2 and TH3), and anterior and lateral to the mid-shank (SH1 and SH2). Marker trajectories were recorded using an 8-camera stereophotogrammetric system (Vicon MX, Vicon Motion Systems Ltd, Oxford, UK, 100 frames per second) with synchronized measurement of the ground reaction forces obtained using two strain-gauge force plates (Bertec Corp., Columbus, OH, USA, 1000 samples per second). Motion tasks included a static standing posture with each foot on the two separate force platforms and five acquisitions of level walking at self-selected speed.

Musculoskeletal models
Three lower limb MSMs, named ALLM, G2392, and LLLM respectively were downloaded from www.simtk.org. ALLM and G2392 were chosen for being widely adopted and cited. LLLM was chosen as being the one that most differed from them in terms of bone geometries, joint constraints, muscular attachment sites and linesof-action, number of muscle bundles, and for being a single lower limb model (Table 1). This last characteristic influences the model estimates because a multibody kinematics optimization is employed.
Each generic MSM, which includes the above-mentioned anatomical markers, was scaled to match the volunteer's anthropometry estimated using the ratio between the lengths of the model segments and those computed from the experimental data. The pelvis was scaled using the distance between the right and left anterior superior iliac spines, and the distance between the mid-points of the anterior and posterior superior iliac spines. The joint centres were located using the marker positions as acquired in a static trial and the Harrington regression equations (Harrington et al., 2007) for the hip joint, the mid-point between the femoral epicondyles for the knee joint, and the mid-point between the malleoli for the ankle joint. The size of the thighs, shanks and feet was scaled using the distances between the hip and knee centres, knee and ankle centres, and heel and second metatarsal head markers, respectively. The technical markers were finally embedded in the scaled MSMs by registering, using the multibody kinematics optimization method, the anatomical markers of each model with the corresponding anatomical markers placed on the volunteer as recorded during the static trial. The segment masses in the model were uniformly scaled to match the total body mass of the participant.
The maximal isometric forces of the muscles represented in the MSMs, which are parameters needed to solve the myoskeletal indeterminacy problem (Viceconti et al., 2006), were uniformly scaled following criteria described in previous studies (Arnold et al., 2013;Laughlin et al., 2011;Mokhtarzadeh et al., 2014). In particular, a scaling factor equal to the ratio between the volunteer lower limb mass, estimated as a percentage of the total mass (De Leva, 1996), and the corresponding generic MSM lower limb mass was used. However, when using ALLM and LLLM during gait, some muscles resulted fully activated, reaching the maximal force values permitted. Given the nature of walking as a sub-maximal motor act, this is an unlikely outcome, so the affected maximal forces defined in the MSMs, were increased by up to a factor of three, confident in the fact that this would not significantly influence the sensitivity analysis of the present study.
One gait cycle was simulated for the participant's dominant lower limb using the standard OpenSim pipeline (Delp et al., 2007). First run was the "inverse kinematics" analysis which uses a multibody kinematics optimization algorithm to determine the joint angles that best fit the experimental trajectories collected during one selected walking trial (Lu and O'Connor, 1999). The RMS difference between the virtual and experimental markers was on average 1.3 cm, 1.2 cm and 0.9 cm for ALLM, G2392 and LLLM, respectively, with a maximum tracking error lower than 4.1 cm, 3.6 cm and 3.6 cm, respectively. The joint moments were calculated through inverse dynamics and decomposed into muscle forces by minimizing the sum of the squared muscle activations while neglecting the forcelength-velocity relationships of muscles (Anderson and Pandy, 2001). The residuals at the hip, knee and ankle were all below 0.06 Nm and hence far less than 1% of the COM height times the magnitude of the measured net external force, which is the limit suggested by Hicks et al. (2015). Finally, joint contact forces were calculated by solving the static equilibrium conditions for each segment. The estimation of the knee contact force was only performed for G2392. This was due to the fact that in both ALLM and LLLM the pose of the patella is defined as a function of the tibiofemoral joint flexion-extension angle, which has been proven to lead to inaccurate estimates of the overall tibio-femoral contact force when computed using the available OpenSim tools (Koehle and Hull, 2008;Wagner et al., 2013). Since implementing ad-hoc tools to perform this calculation was beyond the scope of this study, relevant data will not be reported for these models. All analyses were conducted using OpenSim 3.1 (Delp et al., 2007) and MATLAB scripts (The MathWorks Inc., USA, version 2015a), including the publically available libraries (Barre and Armand, 2014;Mantoan et al., 2015).
All estimated joint angles, joint moments, and muscle and joint contact forces showed good agreement with the literature (Kadaba et al., 1989;Martelli et al., 2014;Modenese and Phillips, 2012;Prinold et al., 2016;Valente et al., 2014; Table 1 Musculoskeletal models used to perform the sensitivity analysis. a Single lower limb model. Wesseling et al., 2015), with slightly increased contact forces observed at the hip, consistent with larger measured ground reaction forces.

Baseline dataset and probabilistic design of the parametric STA
A probabilistic approach, based on a Monte Carlo analysis, was used to evaluate the impact of the STA on the three selected MSMs. For each MSM, a set of marker trajectories was synthetically created as a baseline dataset for this sensitivity analysis. This was achieved by running the "point kinematics" tool which, given the calculated joint angles, provides the global coordinates of the markers that are rigidly attached to the model body segments for each instant in time. As such, these coordinates are time-invariant when observed from their respective anatomical frames.
The STAs for the feet, shanks, lateral femoral epicondyles and pelvis markers were modelled as sinusoidal functions of time described by nine parameters representing amplitude, frequency and phase of each marker's spatial coordinates (Chèze et al., 1995). Their statistical representation was obtained using their mean range7 three standard deviations. The pelvis marker amplitudes were varied nonuniformly for the three spatial coordinates using the mean values and 95% confidence intervals reported by Rozumalski et al. (2007). The shank STA amplitudes were computed using the values suggested by Dumas and Cheze (2009) with the foot amplitudes similarly determined using the values reported by Tranberg and Karlsson (1998).
The STAs for the lateral-thigh markers were modelled as a linear function of the three hip angular rotations and the knee flexion angle (Bonci et al., 2014). Each coordinate of these STAs was described by four coefficients, which were used to multiply the reference hip flexion, abduction, rotation and knee flexion angles respectively, and by one constant (h°, Table 2). The mean value for the statistical distribution of the four coefficients was set to be equal to the values of the ex-vivo dataset of Subject 1 reported in Bonci et al. (2014). The standard deviation was computed using the ratio between the root mean square values of the STA components of the same subject and the average value of the corresponding joint angles over the gait cycle. The mean and standard deviation values for h°were set using the standing joint angle statistics reported in Hemmerich et al. (2006). As a result of the above calculations, 22 sinusoidal STAs and 6 kinematics-dependent STAs were defined for G2392 and ALLM, resulting in a total of 324 stochastic input variables for the Monte Carlo analysis. For the single-leg LLLM, only 13 sinusoidal STAs and 3 kinematics-dependent STAs were used, resulting in 162 input variables for the statistical analysis.
A Latin Hypercube Sampling (LHS) method was then used to generate 500 samples for each of the stochastic variables, reflecting the mean and standard deviation of each variable. The distributions generated were then checked for normality using the Lilliefors test (Lilliefors, 1967). This process produced 500 STA realizations in the local anatomical frames. A coordinate transformation to the laboratory frame was then performed, in order to sum the STA realizations to the reference marker trajectories and create the artefact-affected trajectories. Finally, the artefact-affected trajectories were then iteratively fed to the corresponding MSM. Joint angles, joint moments and muscle and joint contact forces were estimated using the generated artefact-affected trajectories while keeping the same measured ground reaction forces. Joint moments were normalized to the volunteer's mass and muscle and contact forces were expressed as multiples of body weight (BW).
The appropriateness of using 500 as the sample size was assessed via convergence analysis of the entire set of generated input and output variable distributions, where changes observed for samples higher than 300 were found to be below a convergence threshold of 2% of the 50th-85th percentile (Martelli et al., 2015). No discrepancies were observed among the investigated MSMs.

Data analysis
The distribution of the STA realizations of the pelvis and right lower limb was calculated and compared with published STA measurements excluding those used to generate the artefact-affected trajectories (Akbarshahi et al., 2010;Cappozzo et al., 1996;Hara et al., 2014;Maslen and Ackland, 1994;Stagni et al., 2005;Tsai et al., 2009;Wrbaskić and Dowling, 2007).
The sensitivity of each of the three MSMs was determined by calculating the 5th, 50th and 95th percentiles of their output for the right lower limb over the entire gait cycle. The difference between the 95th and 5th percentile variation of each output of interest (hereinafter referred to as the variation interval) was described using maximum, mean and standard deviation values. Relative variations were also quantified by calculating for each output the ratio between the mean variation interval and the range of the 50th percentile.

Results
The marker-depended STA distribution showed good agreement to published STA measurements (Table 3). The STA for the markers in the thigh segment exhibited the largest range of values with mean and standard deviation reaching 39.77 17.6 mm (peak STA value: 46.6 721.7 mm found for the RTH2 marker). The mean STA for the markers on the pelvis was 28.7 76.7 mm (peak STA value: 36.97 8.2 mm found for the ASIS markers), was for those on the shank 11.973.8 mm (peak STA value: 14.5 74.6 mm found for  the RSH1 marker) and was for those on the foot 1.97 0.6 mm (peak STA value: 2.5 7 0.8 mm found for the RHEE marker).
The shape, magnitude and timing of the joint angle and moment time histories were not significantly affected throughout the entire gait cycle (Fig. 1-4). Minor magnitude variations across the different gait phases were observed in the muscle and joint forces. These variations were consistently found in the estimates of the three models, with the highest percentage values occurring for the peak values of the quantities involved (Figs. 5 and 6).
The time histories of the muscle forces estimated using the three MSMs showed similar patterns but different magnitudes. The largest difference was found for the soleus and the gastrocnemius muscles, where the peak values of the 95th percentile were 3.1 BW and 2.1 BW as calculated by the ALLM, 3.3 BW and 1.5 BW as calculated by the G2392 and 1.4 BW and 3.2 BW as calculated by the LLLM, respectively. The STA effect on the estimation of the joint reaction forces was jointdependent, showing the highest effect at the hip and a reduced impact at the ankle. This was consistent across the three MSMs analysed. The maximum force variation intervals at the hip were 1.8 BW, 1.5 BW and 1.6 BW for ALLM, G2392 and LLLM, respectively, while the maximum knee force variation was 0.9 BW for G2392 and the maximum ankle force variations were of about 0.6 BW for all three MSMs. (Fig. 7).
For the pelvis and hip angles, relative variations higher than 30% occurred consistently for the three MSMs, whereas values below 20% were observed for both the knee and ankle angles. The relative variations found for the joint moments and muscle forces ranged between 5% and 25% and were to a great extent consistent across MSMs. For LLLM only, slightly higher values were found for the soleus (38%), the gluteus maximus (31%) and the lateral gastrocnemius (31%). Finally, relative variations ranging from 5% to 15% were consistently found across the three MSMs for the joint contact forces.

Discussion
The aim of this study was to quantify the sensitivity of three different musculoskeletal models to the soft tissue artefact affecting their input positional measurements. This was achieved through a probabilistic analysis, which overall showed that the output variations increased from the ankle to the hip, while the shape and magnitude of the outputs of interest were mostly preserved throughout the entire gait cycle. The observed effects were similar across MSMs.
The STA realizations generated in this study were found to be in line with measured STAs reported in the literature. The magnitude of the STA estimated for the pelvis markers during walking was on average 20 mm (Table 3) which, as expected, was higher than the 17 mm found by Hara et al. (2014) for multiple static standing postures. STA magnitude was higher at the thigh than at the shank (Table 3), as per in previous studies (Akbarshahi et al., 2010;Stagni et al., 2005;Tsai et al., 2009). The average thigh STA peak-to-peak value was higher than those previously reported (23 mm vs 9 mm in Akbarshahi et al., 2010, and14 mm in Tsai et al., 2009), but its 40 mm peak value was in line with that reported by Cappozzo et al. (1996) and Sati et al. (1996). Average values at the shank, were similar to the 8 mm reported by Akbarshahi et al. (2010), with lateral malleus values comparable to the 2 73 mm observed in static positions by Maslen and Ackland, (1994). Finally, the lowmagnitude STA generated at the foot (peak values of 1.97 0.6 mm) confirmed the fluoroscopy-based results of Wrbaskić and Dowling (2007), who found strongly correlated patterns for skin and bone mounted marker trajectories. The sensitivity of the three selected MSMs to the STA was evident for most of the investigated output quantities, with different amplitudes but similar patterns observed for the three models, despite the differences in their bone and joint definitions, muscular-tendon parameters and even number of limbs. Maximum, mean and standard deviation of the output variation intervals increased from the ankle to the hip for most of the variables. When investigating the probabilistic effect of the STA on joint kinematics for different models, El Habachi et al. (2015) found a different pattern, with the highest values reported for the ankle. This disagreement is likely due to the fact that, in contrast with the literature, they considered the pelvis and foot segments affected by a STA modelled with the same amplitude. Our results partially differ also from those of Myers et al. (2015), who also investigated G2392. They observed mean joint angle variation  intervals of 5°, 2°and 6°for hip flexion, knee flexion and ankle dorsiflexion, respectively, whereas we found values of 15°, 8°and 6°, respectively. However, Myers et al. (2015) set a maximum constraint of 15 mm for their probabilistically generated STA. Although this constraint seems to be plausible for the markers used to estimate the foot kinematics, much higher values would be expected for the pelvis and thigh markers, and this may explain the divergence from our results. This may also explain the similar variation of the joint moments at the ankle (0.03 Nm/kg) and the almost doubled variation at the hip and the knee as compared to corresponding variations reported by Myers et al.,(hip: 0.09 Nm/ kg vs 0.24 Nm/kg; knee: 0.07 Nm/kg vs 0.14 Nm/kg).
The present study was affected by some limitations. Firstly, only gait was investigated while different motor tasks exhibiting a larger range of motion such as squatting or running may have shown different sensitivities to the STA. Further studies are however needed to prove this prediction. Secondly, we limited the analysis to data from one representative subject of an adult healthy population and caution should therefore be used when considering the reported results in association with data from pathological or paediatric populations. Their gait kinematics may in fact be characterised by different baseline datasets with the corresponding sensitivity analysis leading to different output variations. Thirdly, we investigated the model sensitivity to STA alone while the interaction between this aspect and other parameters and assumptions (e.g. model anatomy, inertial parameters, and muscle function) may have altered the model sensitivity further. More research is needed to fully address this aspect. Lastly, the STA design used in this study might be limited by the fact that the experimental STA measurements used for the different body segments were obtained with different techniques and this might have affected the different variations observed among the investigated joints.
Despite the above limitations, the present study provides useful information to deal with the unsolved problem of the STA that inevitably propagates to the estimates of MSMs. The reported results  suggest that current MSMs, driven by stereophotogrammetric recordings of skin-mounted marker trajectories, might be effectively used when an overall pattern is more important than an accurate quantitative estimation, such as in comparative cohort-based studies. The amount of observed variation, higher than 30% in some cases, suggests that caution should be exercised in interpreting the results  when MSMs are used for applications requiring a very accurate level of estimation, and in particular when they are used for subject-specific estimation of joint kinetics and bone strains. More research is required to optimize marker sets and their placement, the inverse kinematics algorithm and to develop STA compensation techniques.