- Email: [email protected]
Prediction of ground reaction forces and moments during various activities of daily living
Journal of Biomechanics 47 (2014) 2321–2329
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Prediction of ground reaction forces and moments during various activities of daily living R. Fluit a,n,1, M.S. Andersen b, S. Kolk c, N. Verdonschot a,d, H.F.J.M. Koopman a a
Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente, Enschede, The Netherlands Department of Mechanical and Manufacturing Engineering, Aalborg University, Aalborg, Denmark c Radboud University Medical Centre, Radboud Institute for Health Sciences, Department of Rehabilitation, Nijmegen, The Netherlands d Radboud University Medical Centre, Radboud Institute for Health Sciences, Orthopaedic Research Laboratory, Nijmegen, The Netherlands b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 18 April 2014
Inverse dynamics based simulations on musculoskeletal models is a commonly used method for the analysis of human movement. Due to inaccuracies in the kinematic and force plate data, and a mismatch between the model and the subject, the equations of motion are violated when solving the inverse dynamics problem. As a result, dynamic inconsistency will exist and lead to residual forces and moments. In this study, we present and evaluate a computational method to perform inverse dynamics-based simulations without force plates, which both improves the dynamic consistency as well as removes the model's dependency on measured external forces. Using the equations of motion and a scaled musculoskeletal model, the ground reaction forces and moments (GRF&Ms) are derived from three-dimensional full-body motion. The method entails a dynamic contact model and optimization techniques to solve the indeterminacy problem during a double contact phase and, in contrast to previously proposed techniques, does not require training or empirical data. The method was applied to nine healthy subjects performing several Activities of Daily Living (ADLs) and evaluated with simultaneously measured force plate data. Except for the transverse ground reaction moment, no significant differences (P40.05) were found between the mean predicted and measured GRF&Ms for almost all ADLs. The mean residual forces and moments, however, were significantly reduced (P40.05) in almost all ADLs using our method compared to conventional inverse dynamic simulations. Hence, the proposed method may be used instead of raw force plate data in human movement analysis using inverse dynamics. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Ground reaction forces and moments Musculoskeletal model Inverse dynamics Dynamic consistency Activities of daily living
1. Introduction Inverse dynamics based simulations on musculoskeletal models are a commonly used method for the analysis of human movement. Despite widespread use, it is well known that solutions obtained with inverse dynamics are sensitive to inaccuracies in the various input variables (Pamies-Vila et al., 2012; Riemer et al., 2008). Errors can stem from estimating body segment parameters (Pearsall and Costigan, 1999; Rao et al., 2006), estimating joint parameters (Schwartz and Rozumalski, 2005), skin movement artifacts (Leardini et al., 2005), noise on skin-mounted marker data (Richards, 1999), estimating the center of pressure (Schmiedmayer and Kastner, 1999) or force plate calibration (Collins et al., 2009). Consequently, when solving the inverse dynamics problem, the equations of motion are violated,
n Correspondence to: Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. Tel.: þ31 53 489 4362; fax: þ31 53 489 2287. E-mail address: r.fl[email protected] (R. Fluit). 1 Visiting address: Campus University, Horstring W213.
http://dx.doi.org/10.1016/j.jbiomech.2014.04.030 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
resulting in dynamic inconsistency, a condition with residual forces and moments (Kuo, 1998). Several algorithms have been proposed that reduce or eliminate the residual forces and moments, such as the least-squared optimization (Kuo, 1998), the residual elimination algorithm (REA) (Thelen and Anderson, 2006) and the residual reduction algorithm (RRA) (Delp et al., 2007). These algorithms adjust the kinematics, ground reaction forces (GRFs) and/or body segment parameters, thereby improving the dynamic consistency. Unfortunately, they have shortcomings too: the REA was shown to dramatically change torso angles for movements longer than 0.5 s (John et al., 2007). For the RRA, an adjustment in the joint angles of up to five degrees is considered reasonable (OpenSim User's Guide). Since these differences are larger than the minimal detectable change for most of the joint angles (Wilken et al., 2012), these adjustments may not be defendable. Alternatively, dynamic consistency can be improved by deriving the ground reaction forces and moments (GRF&Ms) from three-dimensional full-body motion using the equations of motion (Audu et al., 2007; Choi et al., 2013; Oh et al., 2013; Ren et al., 2008; Robert et al., 2013). An additional advantage of this method
2322
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
is that it enables inverse dynamic analysis for studies without force plate data, for example ambulatory measurements with inertial measurements only (e.g. Luinge and Veltink, 2005) or motion capture during treadmill walking (e.g. Hesse et al., 1999). A difficulty of this method is the indeterminacy problem during a double contact phase when the system defines a closed kinetic chain. To overcome this problem, Audu et al. (2003) used optimization techniques to compute the GRF&Ms for different static postures of a standing bipedal model, which were later validated against measured data (Audu et al., 2007). However, it is unknown whether this method is valid for dynamic movements. Ren et al. (2008) introduced a smooth transition assumption to solve the indeterminacy problem. The smoothing functions were based on empirical data and, therefore, the smooth transition assumption may not be applicable for movements other than those present in the empirical data. Oh et al. (2013) and Choi et al. (2013) solved the indeterminacy problem using an artificial neural network. Also their method requires training data, which is not always present. Robert et al. (2013) tested several optimization methods to predict the external contact loads during sit-to-stand movements. Although their method does not require empirical or training data, their contact configurations are simplified and the method was validated for sit-to-stand motion only. Therefore, the purpose of this paper is to demonstrate a universal method for predicting the GRF&Ms based on measured kinematic data only, which is applicable to a variety of Activities of Daily Living (ADLs), and in which the indeterminacy problem during the double contact phase of gait and gait-related ADLs is solved without the use of empirical or training data. The predicted GRF&Ms were evaluated with simultaneously measured force plate data. For the trials where subjects walked at self-selected comfortable walking speed, a sensitivity study was performed to evaluate the effects of the chosen muscle recruitment strategy and parameters of the ground contact model on the accuracy of the predictions.
The local ethics committee approved the study protocol and informed consent was obtained from all subjects prior to the study. 2.2. Instrumentation A six-camera digital optical motion capture system (Vicon MX, Oxford, UK) was used to capture 53 retro-reflective markers placed on the body at 100 Hz sampling rate. Markers were placed according to the Vicons Plug-in-Gait model (Vicons, 2002) with additional markers on the anterior side of the thigh and shank at 1/3 and 2/3 segment length, the medial femoral epicondyle, the medial malleolus and on the fifth metatarsal head of the foot. Two custom-built force plates (AMTI, Watertown, MA, USA), embedded level in the laboratory floor, measured GRF&Ms at 1000 Hz. For the stair negotiation trials, a custom-built four step staircase was used with a riser height of 180 mm and a thread depth of 250 mm. For these trials, GRF&Ms were collected for the initial or last step at floor level, the first step using a box on a force plate and the second step using an additional custom-built portable force plate (ForceLink B.V., Culemborg, the Netherlands). 2.3. Experimental protocol Three-dimensional motion capture with synchronized force plate recordings of several ADLs were collected: level walking at comfortable (CWS), slow (CWS 30%) and fast (CWS þ 30%) walking speed, walking over a 10 (OBS10), 20 (OBS20) and 30 (OBS30) cm obstacle, gait initiation (GI) and termination (GT) starting or ending with the dominant (DL) or non-dominant leg (NDL), deep squatting (DS) and stair ascent (SA) and descent (SD). The dominant leg was determined by asking the participant with which leg they would kick a ball. All ADLs were performed barefooted. For the stair negotiation trials, subjects were instructed to perform the activity without the use of handrails. For the obstacle trials, subjects were allowed to choose their own strategy to step over the obstacle placed on the walkway, as long as the force plate was hit clean. For each ADL, subjects performed several practice trials before the actual trials were recorded. 2.4. Data processing
2. Methods
3D position data of the markers and the force plate data were low-pass filtered using a second order, zero-phase Butterworth filter with a cut-off frequency of 12 and 15 Hz, respectively. For the CWS, CWS 30% and CWS þ 30% trials, three experimental data sets per subject were used for validating the predictions of the GRF&Ms. For all other ADLs, one experimental data set per subject was used, yielding a total of 261 trials. In total, 257 trials were recorded correctly, of which another four trials were excluded due to missing markers.
2.1. Subjects
2.5. Prediction of GRF&Ms
Nine healthy subjects (4 males and 5 females; age: 41.6 7 15.9 yr; height: 1.74 70.12 m; weight: 73.0 7 11.1 kg, Body Mass Index: 23.9 7 2.0 kg/m2) with no history of musculoskeletal disorders volunteered for the study at the Rehabilitation Department of the Radboud University Medical Centre, Nijmegen, the Netherlands.
For all ADLs, the only significant external forces and moments were the GRF&Ms. These can be estimated using the Newton–Euler equations of motion, which state that the sum of all external forces balances the sum of the massacceleration products of all individual body segments (Greenwood, 1988). Hence,
Fig. 1. Left: Visualization of the 12 contact points for each foot. Points were defined at the medial and lateral side of the heel, at the base of the first and fifth metatarsal bone, at the head of each metatarsal bone and at the big, second and fifth toe. Right: Visualization of the characteristic strength function cp;i (Eq. (1)) and the smoothed strength function (Eq. (2)) for the right heel contact node during a gait cycle. Heel contact and toe-off of the right leg are abbreviated as HCR and TOR respectively and, analogously, for the left leg as HCL and TOL.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329 during the single support phase of gait, the GRF&Ms can be obtained directly from the Newton–Euler equations when full-body kinematics is available. During a double contact phase; however, the problem of determining the GRF&Ms under each foot becomes underdetermined and they cannot be resolved from the Newton–Euler equations alone. To overcome this, the computation of the GRF&Ms was made part of the muscle recruitment algorithm by introducing artificial muscle-like actuators at 12 contact points under each foot (Fig. 1, left side). At each contact point p, five artificial muscle-like actuators were added. One actuator was aligned with the vertical direction of the force plate (z) and was able to generate a normal force. The other two pairs of actuators were aligned with the medio-lateral (x) and antero–posterior (y) direction of the force plate and were able to generate positive or negative static friction forces. The force exerted by each actuator was defined as qF max , with q the activation level and F max the strength. A dynamic contact model was used to define the strength profile of each actuator, ensuring they were able to generate a reaction force only if their contact point p was close to the floor and almost without motion. For each contact point p, a characteristic strength function cp;i was defined ( F max if pz o zcrit and jjpjj o vcrit ð1Þ cp;i ¼ 0 otherwise where F max ¼ 0:4 BW, zcrit ¼ 0:03 m (level of the floor was set to z¼ 0), vcrit ¼ 0:8ms and i the index used for the step of the simulation (each step was 0.01 s). These values were based on pilot tests and a sensitivity analysis of the implications of these values was performed (Section 2.8). For the transitions of cp;i at the discrete
2323
time step i ¼t, the characteristic strength function was smoothed (Fig. 1, right side) according to 1 1 cp;i ¼ ðt 3…t þ 8Þ ¼ F max 1 7 cos πn with n ¼ 1…12 ð2Þ 2 12 The smoothing function was based on pilot tests and was introduced to prevent discontinuities in the estimated GRF&Ms during the transitions from being completely determined (single stance) to being underdetermined (double stance) and vice versa. The activation level q of each artificial muscle-like actuators was solved as part of the muscle recruitment problem, thereby determining the magnitude of the GRF&Ms. The solver did not distinguish between single and double contact phases.
2.6. Inverse dynamic model All simulations were performed using a 28 degrees-of-freedom (DOF) full body model (Fig. 2) as available in the AnyBody Modeling System (version 5.3.1, AnyBody Technology A/S, Aalborg, Denmark) (Damsgaard et al., 2006). For the lower limbs the geometry of the new Twente Lower Extremity Model was used (Fluit et al., 2013). The segment masses were linearly scaled using data from Winter (2009). The model marker positions and segment lengths were determined using a parameter optimization algorithm (Andersen et al., 2010) applied to an arbitrary CWS trial. Subsequently, inverse kinematics for all other trials were solved using the optimized parameters. Muscles were included for the lower limbs (Klein Horsman et al., 2007) whereas for the upper extremities, ideal joint torque generators were used. Each leg consisted of 55 muscle-tendon parts described by 166 Hill-type muscle elements. To obtain the muscle forces and thus the predicted GRF&Ms, the muscle recruitment problem was solved through static optimization, minimizing the sum of the cubed muscle activations. To improve numerical stability, additional muscle-like actuators were attached at the origin of the pelvis segment, which could generate small residual forces and moments up to 10 N or Nm. The activation level of these muscles was solved as part of the muscle recruitment optimization problem, so as to minimize their contribution. Skeletal muscles and muscle-like actuators for the GRF&Ms and residuals were weighted equally. However, the strength of the actuators for the GRF&Ms was high compared to the skeletal muscles whereas the actuators for the residuals had a low strength.
2.7. Accuracy analysis
Fig. 2. Visualization of the skeletal model and the 28 degrees of freedom (DOFs).
To assess the quality of the predictions, the predicted GRF&Ms were compared with the force plate data. Except for the DS trials, the timespan over which the GRF&Ms were compared was defined as the time for which the measured vertical GRF was larger than 5% of its maximum. For the DS trials, this timespan was defined from 0.5 s before until 0.5 s after the knee of the dominant leg was flexed at 151 flexion angle. The ground reaction moments (GRMs) were calculated about the point on the force plate surface vertically below the ankle joint. Following the definitions of Ren et al. (2008), the differences were quantified using the Root Mean Squared Difference (RMSD) and, additionally for the CWS and CWS þ30% trials, the relative RMSD (rRMSD). The Pearson correlation coefficient ρ was calculated as well and categorized (in absolute value) as ρ r 0:35, 0:35 o ρ r 0:67, 0:67o ρ r 0:9, 0:9 o ρ to be weak, moderate, strong or excellent correlations respectively (following Taylor, 1990).
Table 1 Comparison of measured and predicted GRFs using the curve magnitude (M) and phase (P) difference (mean (SD)) and Pearson correlation coefficient ρ. Results are averaged over both legs and all subjects. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to the two-tailed Wilcoxon Signed Rank Test (α¼ 0.05). Activity
CWSþ30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
Vertical GRF
Antero-posterior GRF
P (%)
M (%)
3.0 2.4 2.7 2.4 2.3 2.5 3.2 2.3 2.5 2.7 2.7 3.1 4.7
3.4 2.9 3.0 3.4 2.8 2.1 2.0 2.2 2.1 1.4 2.5 4.0 1.2
(0.4) (0.6) (0.4) (0.7) (0.2) (0.6) (0.9) (0.3) (0.2) (0.5) (0.6) (1.3) (0.6)
ρ () (1.5) (0.8) (0.9) (2.2) (2.5) (1.7) (1.1) (0.5) (1.2) (1.0) (1.4) (1.3) (0.9)
0.954 0.967 0.957 0.976 0.980 0.976 0.947 0.962 0.957 0.943 0.621 0.941 0.895
P (%)
n† n†
n†
n† n†
9.0 8.4 9.3 6.0 8.1 6.3 9.5 7.5 8.0 7.9 35.0 26.0 22.3
(1.3) (1.9) (1.2) (1.1) (2.2) (1.5) (1.2) (1.2) (1.3) (1.3) (15.1) (7.6) (2.8)
Medio-lateral GRF
M (%)
ρ ()
P (%)
4.1 26.4 15.8 1.1 13.4 3.3 11.0 6.4 6.1 5.1 71.9 57.5 15.1
0.948 0.965 0.957 0.968 0.944 0.962 0.902 0.966 0.959 0.969 0.202 0.436 0.637
11.8 7.6 9.8 16.7 11.9 17.8 15.2 9.8 11.1 12.4 34.7 23.2 26.2
(8.3) (9.1) (9.3) (11.5) (7.4) (5.8) (6.6) (6.1) (11.3) (7.3) (139.6) (65.3) (12.6)
n†
n†
ρ ()
M (%) (2.6) (1.9) (2.3) (12.4) (2.2) (6.9) (4.9) (3.0) (2.8) (3.5) (14.0) (4.5) (2.1)
1.0 1.9 1.9 8.6 5.1 4.7 6.1 8.6 12.0 14.1 11.2 9.9 10.3
(10.8) (4.0) (6.9) (9.6) (5.5) (17.7) (7.9) (8.7) (12.4) (15.4) (43.8) (14.7) (12.4)
0.753 0.866 0.818 0.898 0.842 0.863 0.726 0.836 0.810 0.713 0.261 0.609 0.455
2324
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
Further, differences were quantified using the Spraque and Geers curve magnitude (M) and phase (P) difference (Schwer, 2007; Sprague and Geers, 2003). The metrics M and P are expressed in percentages and are designed to give zero when the curves are identical. Additionally, two-tailed Wilcoxon signed rank tests (α ¼ 0.05) were performed to test if the medians of the absolute mean and absolute maximum of the measured GRF&Ms were significantly different than those of the predicted GRF&Ms. Similar comparisons were performed for a subset of the computed joint moments. The absolute mean and maximum of the residual forces and moments were calculated for each trial for the timespan in which ground contact was fully defined. One tailed Wilcoxon signed rank tests (α ¼ 0.05) were performed to test if the medians of the absolute mean and absolute maximum of the residual forces and moments using the predicted GRF&Ms were significantly lower than those using the measured GRF&Ms.
The accuracy of the dynamic contact model (Eqs. (1) and (2)) was assessed by comparing heel-strike (HS) and toe-off (TO) events for the CWS trials, defined as the time at which the vertical GRF reached 5% of the maximum.
2.8. Sensitivity analysis For all CWS trials, a sensitivity analysis of the parameters of the dynamic contact model (F max , zcrit and vcrit ) and the muscle recruitment criterion on the prediction accuracy was made. For F max , zcrit and vcrit , simulations were performed using values at 70% and 130% of their original magnitude (Eq. (1)). The sensitivity was quantified as the change in M, P, ρ and rRMSD, where a negative change indicated an improvement in the predicted GRF&Ms and vice versa.
Fig. 3. Calculated GRFs, normalized to body weight (BW), and GRMs, normalized to BW and body height (BH) ( 7 1 SD around mean (shaded area)) compared with force plate data (mean (thick line) 71 SD (thin lines)), averaged over all subjects. For the CWS, OBS30 and SD trials, the GRF&Ms of both legs are shown and the heel contact (HC) and toe-off (TO) events are indicated. For the DS trials, only the GRF&Ms of the dominant leg are shown and at the start (S) and end (E) of the squat movement the knee of the dominant leg was flexed at 15 1 flexion angle.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
3. Results The model showed excellent predictions for the vertical GRF (ρ ranging from 0:621 0:980, median 0:957) and the anteroposterior GRF (ρ ranging from 0:202 0:969, median 0:957) for almost all activities (Table 1, Fig. 3). The magnitude of the vertical GRF was slightly but consistently underestimated (M ranging from 1.2% to 4.0%, median 2:5%), whereas the magnitude of the antero-posterior GRF was consistently overestimated (M ranging from 1:1% to 71:9% median 11:0%). Nevertheless, no significant difference (P 4 0:05) was found between the medians of the absolute mean measured GRFs and those of the predicted GRFs for all activities. Strong correlations were obtained for the medio-lateral GRF for almost all activities (ρ ranging from 0:261 to 0:898, median 0:810). In general, the most extreme activities (DS, SA and SD) showed the weakest correlations, for instance the antero-posterior (ρ ¼ 0:202) and mediolateral (ρ ¼ 0:261) GRFs for the DS trials. The model showed strong correlations for the sagittal GRM (ρ ranging from 0:506 to 0:922, median 0:789) and frontal GRM (ρ ranging from 0:199 to 0:801, median 0:668) for most of the activities (Table 2). The transverse GRM showed the largest differences (ρ ranging from 0:155 to 0:782, median 0:588). A significant difference (P o0:05) was found between the medians of the absolute mean measured transverse GRM and those of the predicted transverse GRM for all activities.
2325
The RMSDs between measured and predicted GRF&Ms were similar across different ADLs (Table 3). Regarding the joint moments, strong correlations were found for almost all activities for hip flexion (ρ ranging from 0:521 to 0:972, median 0:776), hip abduction (ρ ranging from 0.257 to 0.946, median 0:880), hip external rotation (ρ ranging from 0.809 to 0.943, median 0:871), knee flexion (ρ ranging from 0:548 to 0:911, median 0:737) and ankle dorsiflexion (ρ ranging from 0:647 to 0:927, median 0:820). A weak correlation was found for the hip abduction moment for the DS trials only (ρ ¼ 0:258, see Table A1 and Fig. 4). For a subset of activities, Fig. 4 demonstrates the joint moments at the DOFs in the dominant leg only using measured and predicted GRF&Ms (Table A1 contains a complete overview for all activities and both legs). Using our method, the medians of the absolute mean residual forces, transverse and frontal moments were significantly reduced compared to those obtained with a conventional inverse dynamics simulation for almost all activities (P o 0:05). Averaging over all activities and subjects, the absolute mean residual forces and moments applied at the pelvis were about 2 3 N or Nm using the predicted GRF&Ms, thereby improving dynamic consistency (Table A2). The vertical residual force showed the strongest reduction of 90%. Only the residual sagittal moment showed a non-significant reduction of 33%, probably because this moment was already small using traditional inverse dynamics (4:3 Nm, averaging over all activities and subjects).
Table 2 Comparison of measured and predicted GRMs using curve magnitude (M) and phase (P) difference (mean (SD)) and Pearson correlation coefficient ρ. GRMs are expressed at the virtual point on the force plate surface below the ankle joint. Results are averaged over both legs and all subjects. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (α¼ 0.05). Activity
Sagittal GRM P (%)
CWS þ 30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
†
n
6.8 9.6 7.4 9.3 11.4 8.5 10.4 7.1 7.7 8.7 15.6 12.4 8.3
Frontal GRM M (%)
(2.3) (2.3) (1.8) (2.9) (4.7) (2.5) (3.0) (1.8) (2.6) (1.7) (9.8) (3.3) (3.1)
6.6 15.2 8.8 0.4 4.8 14.6 3.7 6.3 6.4 4.4 50.5 14.0 5.5
(28.2) (6.1) (6.3) (7.1) (23.2) (12.1) (10.9) (7.2) (6.0) (7.4) (59.0) (14.7) (9.1)
ρ ()
P (%)
0.929 0.872 0.922 0.727 0.878 0.680 0.777 0.915 0.863 0.801 0.506 0.647 0.830
16.4 10.7 13.4 8.4 10.6 13.0 14.0 12.7 14.6 15.3 24.9 9.7 9.7
Transverse GRM ρ ()
M (%) (7.3) (3.4) (5.8) (2.3) (3.9) (4.8) (4.8) (4.6) (4.3) (6.7) (6.7) (2.0) (2.9)
10.5 7.3 10.8 9.2 7.2 0.7 0.6 3.7 8.3 15.3 37.1 17.4 15.7
(15.6) (11.2) (10.0) (15.3) (9.1) (18.8) (27.9) (20.5) (21.7) (28.6) (23.2) (13.7) (11.7)
0.540 0.777 0.684 0.703 0.801 0.488 0.539 0.694 0.651 0.641 0.199 0.781 0.781
P (%) n† n† n† n† n† n† n† n† n† n† n† n† n†
33.1 46.0 36.8 26.8 30.9 68.7 56.1 33.9 33.2 33.2 42.5 53.9 67.0
ρ ()
M (%) (7.5) (9.5) (6.3) (10.9) (7.6) (10.4) (8.6) (7.2) (8.4) (4.4) (7.2) (8.6) (5.6)
56.7 46.5 48.3 33.3 47.6 79.4 52.0 46.2 52.0 50.6 91.6 33.8 61.0
(20.6) (19.4) (14.6) (27.6) (18.2) (3.4) (41.1) (13.0) (11.7) (13.7) (147.1) (22.5) (15.4)
0.762 0.598 0.704 0.588 0.525 0.128 0.157 0.782 0.759 0.657 0.230 0.155 0.155
Table 3 Absolute differences (mean (SD)) between the measured and predicted GRF&Ms, quantified using the Root Mean Squared Difference (RMSD) following the definitions of Ren et al. (2008). Activity CWS þ 30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
Vertical GRF (N/kg) 0.85 (0.17) 0.64 (0.15) 0.74 (0.13) 0.54 (0.18) 0.63 (0.10) 0.52 (0.12) 0.91 (0.60) 0.64 (0.17) 0.68 (0.12) 0.74 (0.24) 0.81 (0.36) 0.95 (0.25) 0.34 (0.12)
Antero-posterior GRF (N/kg) 0.43 (0.06) 0.30 (0.05) 0.38 (0.07) 0.16 (0.04) 0.29 (0.07) 0.18 (0.06) 0.32 (0.08) 0.32 (0.06) 0.34 (0.07) 0.33 (0.07) 0.30 (0.22) 0.21 (0.06) 0.22 (0.07)
Medio-lateral GRF (N/kg) 0.23 (0.07) 0.12 (0.03) 0.17 (0.04) 0.19 (0.07) 0.17 (0.03) 0.23 (0.07) 0.21 (0.06) 0.18 (0.05) 0.20 (0.07) 0.22 (0.07) 0.24 (0.16) 0.22 (0.04) 0.42 (0.17)
Sagittal GRF (Nm/kg) 0.17 (0.07) 0.22 (0.07) 0.18 (0.05) 0.11 (0.06) 0.19 (0.05) 0.16 (0.06) 0.20 (0.06) 0.23 (0.11) 0.22 (0.08) 0.23 (0.08) 0.22 (0.05) 0.21 (0.05) 0.09 (0.02)
Frontal GRM (Nm/kg) 0.13 (0.03) 0.08 (0.02) 0.11 (0.02) 0.05 (0.03) 0.07 (0.02) 0.06 (0.02) 0.09 (0.03) 0.12 (0.09) 0.11 (0.03) 0.12 (0.04) 0.17 (0.07) 0.14 (0.05) 0.11 (0.05)
Transverse GRM (Nm/kg) 0.28 (0.08) 0.21 (0.07) 0.22 (0.06) 0.08 (0.06) 0.17 (0.06) 0.16 (0.06) 0.24 (0.10) 0.23 (0.06) 0.22 (0.06) 0.21 (0.07) 0.10 (0.03) 0.19 (0.06) 0.13 (0.04)
2326
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
The sensitivity of the recruitment criterion on the values of P and ρ was small and negligible. With respect to M and rRMSD, none of the recruitment criterions performed better than the others (Table 4). The sensitivity of F max , zcrit and vcrit on the values of P and ρ was also small and negligible. vcrit had the largest effects on M and rRMSD, yielding clearly better predictions for the antero-posterio GRF and sagittal GRM for the 30% case (Table 5). For the CWS and CWS þ30% trials, the RMSD and rRMSD showed similar trends as previous literature (Oh et al., 2013; Ren et al., 2008) (Table 6). Our model provided more accurate predictions for the medio-lateral GRF and frontal GRM than Ren et al. (2008) for the CWS trials, although even better results were obtained by Oh et al. (2013). The dynamic contact model was able to predict heel strike (HS) and toe-off (TO) with reasonable accuracy. Averaging over the CWS trials, HS was predicted consistently early with an error of 28 7 13 ms, and TO was predicted consistently late with an error of 16 77 ms.
4. Discussion In this paper, we demonstrated a universal method to predict the GRF&Ms using kinematic data and a scaled musculoskeletal model only, applied to a variety of ADLs. The method entails a dynamic contact model and optimization techniques to solve the indeterminacy problem during a double contact phase. In general, reasonably good results were obtained for all analyzed activities. However, weak or even negative correlations were found for the transverse GRM for multiple activities. Averaging over all subjects and activities, the magnitude of the transverse GRM was underpredicted with about 50%. These differences are likely caused by the hinged knee, which does not allow for tibial rotation. Previous research has reported a range of 16.21 for this DOF during normal walking, in which the internalexternal knee moment was in the same direction as the corresponding rotation (Andriacchi and Dyrby, 2005). Since for normal gait the internal–external knee moment is primarily determined by the transverse GRM (Schache et al., 2007), omitting this DOF
Fig. 4. Computed joint moments, normalized to body weight (BW) and body height (BH), using the predicted GRF&Ms ( 7 1 SD about mean (shaded area)) compared to those obtained using the measured GRF&Ms (mean (thick line) 7 1 SD (thin lines)), averaged over all subjects. For the CWS, OBS30 and SD trials, the joint moments of the dominant leg are shown from heel contact (HC) until toe-off (TO). For the DS trials, the joint moments of the dominant leg are shown and at the start (S) and end (E) of the squat movement the knee of the dominant leg was flexed at 15 1 flexion angle. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (α¼ 0.05).
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
2327
(Ren et al., 2008). Our sensitivity analysis on the contact model showed that consistent over- and underestimation can be reduced, indicating that our method is not necessarily inferior to traditional inverse dynamics in clinical analysis of gait. Compared to previous literature (Table 6), our predictions are more accurate than Ren et al. (2008) for the medio-lateral GRF (rRMSD of 14:9% and 16:6% versus 20:0% and 24:1% for the CWS and CWSþ30% trials, respectively) as well as the frontal GRM (rRMSD of 22:9% and 27:1% versus 32:5% and 37:7% for the CWS and CWSþ30% trials, respectively). When comparing our results with Oh et al. (2013) for the CWS trials, we obtained a much higher rRMSD for the transverse GRM (40:6% versus 25:5%) whereas the rRMSD for the other elements of the GRF&Ms was only roughly 2% higher. Here it should be noted that the method of Oh et al. (2013) involved an algorithm that required training data to relate
will result in smaller magnitudes of the transverse GRM (see Fig. 3 for the CWS activity). Furthermore, weak correlations were found for the antero-posterior, medio-lateral GRF and frontal GRM for the DS trials. A likely explanation is that for this task, the magnitude of these forces and moments are low, in which case noise has a larger influence on the correlation. Additionally, for the DS task, the total external forces and moments need to be distributed over both feet during the complete movement, which may have increased errors. The joint moments computed by our method were similar to those using traditional inverse dynamics. Although our method underestimated the vertical GRF and overestimated the anteroposterior GRF for all activities and thus may have induced errors in the joint moments, it is likely that the joint moments computed with traditional inverse dynamics suffer from similar errors
Table 4 Sensitivity of recruitment criterion expressed as the change in curve magnitude difference M and change in relative Root Mean Squared Difference (rRMSD). Changes are expressed with respect to the results obtained with the simulations using P ¼3 as recruitment criterion. A negative change indicates an improvement in the predicted GRF&Ms and vice versa. The Min/Max Aux criterion is a composite recruitment strategy consisting of two terms: a term which minimizes the maximum muscle stress and an auxiliary quadratic term. P ¼1
Vertical GRF Antero-posterior GRF Medio-lateral GRF Sagittal GRM Frontal GRM Transverse GRM
P¼2
P¼4
Min/Max Aux
M (%)
rRMSD (%)
M (%)
rRMSD (%)
M (%)
rRMSD (%)
M (%)
rRMSD (%)
0.3 5.0 3.7 2.9 0.7 2.6
0.2 0.8 0.7 1.0 2.1 0.0
0.1 2.0 0.9 0.1 0.3 2.9
0.1 0.3 0.6 0.7 0.9 0.9
0.0 1.7 0.7 1.2 0.6 1.4
0.0 0.1 0.2 0.3 0.1 0.1
0.2 1.4 0.0 0.4 0.9 1.1
0.1 0.3 0.1 1.0 0.2 0.6
Table 5 Sensitivity of the parameters of the ground contact model (as defined in Eq. (1)) expressed as the change in curve magnitude difference M and change in relative Root Mean Squared Difference (rRMSD). A negative change indicates an improvement in the predicted GRF&Ms and vice versa. F max
vcrit
30%
Vertical GRF Antero-posterior GRF Medio-lateral GRF Sagittal GRM Frontal GRM Transverse GRM
þ30%
zcrit
30%
þ 30%
30%
þ 30%
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
0.5 1.7 0.1 1.1 2.3 1.2
0.2 0.3 0.2 0.6 0.0 0.1
0.3 1.6 0.1 1.0 1.7 1.0
0.2 0.3 0.1 0.6 0.1 0.2
0.4 5.5 0.5 4.1 1.5 1.5
0.4 0.9 0.2 1.7 0.3 1.2
0.4 3.2 0.6 3.7 1.0 0.6
0.5 0.4 0.2 1.5 0.2 0.6
0.2 1.9 0.0 0.8 0.4 0.9
0.1 0.4 0.1 0.9 0.2 0.1
0.2 1.3 0.0 0.6 0.4 0.7
0.1 0.3 0.1 0.7 0.2 0.1
Table 6 Differences between measured and predicted GRF&Ms for the CWS and CWSþ 30% trials, quantified using the Root Mean Squared Difference (RMSD) and relative RMSD (rRMSD), following the definitions of Ren et al. (2008). Oh et al. (2013) did not investigate CWS þ30% trials. Ren et al. (2008) investigated fast walking trials, where subjects walked about 27% faster than their normal walking speed, which is similar to our CWS þ 30% trials. Method Participants
Smooth transition assumption (Results from Ren et al., 2008) N¼ 3 CWS þ 30%
CWS RMSD (N/kg or Nm/kg)
Artificial neural network Optimization approach (Results from Oh et al., 2013) (Method proposed in this study) N ¼5 N ¼9
rRMSD (%) RMSD (N/kg or Nm/kg)
Vertical GRF 0.710 (0.190) 5.6 (1.5) AP GRF 0.473 (0.068) 10.9 (0.83) ML GRF 0.191 (0.034) 20.0 (2.7) Sagittal GRM 0.199 (0.106) 12.2 (4.8) Frontal GRM 0.148 (0.013) 32.5 (4.3) Transverse GRM 0.039 (0.015) 26.2 (9.4)
1.012 0.507 0.335 0.225 0.235 0.075
(0.223) (0.085) (0.012) (0.029) (0.018) (0.009)
CWS
CWS þ 30%
CWS
rRMSD (%)
RMSD (N/kg or Nm/kg)
rRMSD (%)
RMSD (N/kg rRMSD (%) RMSD (N/kg or Nm/kg) rRMSD (%) or Nm/kg)
6.6 8.9 24.1 14.2 37.7 41.7
0.649 0.154 0.040 0.081 0.052 0.032
5.8 7.3 10.9 9.9 22.8 25.5
0.74 0.38 0.17 0.18 0.11 0.22
( 1.6) (1.3) (1.6) (1.4) (1.0) (1.8)
(0.182) (0.057) (0.022) (0.045) (0.029) (0.018)
(1.0) (0.8) (1.8) (1.9) (4.9) (4.5)
(0.13) (0.07) (0.04) (0.05) (0.02) (0.06)
6.6 9.3 14.9 12.4 22.9 40.6
(1.1) (2.0) (3.4) (3.5) (5.9) (11.3)
0.85 0.43 0.23 0.17 0.13 0.28
(0.17) (0.06) (0.07) (0.07) (0.03) (0.08)
6.9 8.5 16.6 10.4 27.1 38.4
(1.3) (1.6) (4.6) (3.7) (9.0) (10.9)
2328
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
kinematic measures with the GRF&Ms during the double contact phase. In contrast, our method can be used directly without any training data, since optimization techniques were used to solve the indeterminacy problem. Our dynamic contact model was able to predict HS and TO with an error of 28 7 13 ms and 16 7 7 ms respectively for the CWS trials. Oh et al. (2013) used the dynamic contact model of O'Connor et al. (2007), which predicted HS and TO with errors of 16 7 15 ms and 9 7 15 ms, respectively. Our contact model showed larger errors compared to O'Connor et al. (2007), probably because we did not use optimal values for vcrit and zcrit . For example, when vcrit was reduced with 30%, the errors reduced to 22 7 10 ms and 13 7 7 ms for HS and TO, respectively. The magnitude of these errors is similar to O'Connor et al. (2007). Several limitation of this work should be noted. First, our method requires full body kinematics and a corresponding musculoskeletal model, which is not always available. However, Robert et al. (2013) previously showed that a similar optimization approach at the level of joint moments for sit-to-stand motions provided reasonable predictions as well. Second, the knee was modeled as a hinge, which affected the transverse GRM, and the foot was modeled as a single segment, which may have affected the contact dynamics. Third, it has been reported that soft tissue artifacts (STA) are the most significant source of error in human movement analysis (Leardini et al., 2005; Alexander and Andriacchi, 2001). STA may have affected our joint kinematics, predominantly in the non-sagittal plane, and our estimation of the segment lengths. It is likely that these errors have propagated in the predicted GRF&Ms. Fourth, for the skeletal muscles, excitationactivation dynamics were not modeled. Since the actuators for the GRF&Ms, together with the skeletal muscles, were part of the recruitment criterion, inclusion of excitation-activation dynamics might have altered the predicted GRF&Ms. We believe that this may only be the case for activities other than normal walking, since Anderson and Pandy (2001) showed that static and dynamic optimization are practically equivalent for normal walking. Considering the above mentioned possible sources of error, our method could be improved by including STA compensation (Leardini et al., 2005; Alexander and Andriacchi, 2001), imagebased subject-specific geometry (Scheys et al. 2005) and a more advanced and realistic knee joint (Vanheule et al., 2013) or foot model (Oosterwaal et al., 2011). We believe that techniques for creating subject-specific geometry and the development of more advanced joint models should receive priority in future research. In conclusion, reasonably good estimates of GRF&Ms and joint moments for a variety of ADLs can be obtained solely from whole body kinematics, i.e. without the use of training or empirical data. This work may be useful in the development of ambulatory measurement systems using inertial measurement units only (Luinge and Veltink, 2005), motion capture during treadmill walking (e.g. Hesse et al., 1999) or in the field of movement predictions using musculoskeletal optimization (Rasmussen et al., 2000).
Conflict of interest None of the authors have any financial or personal conflict of interest with regard to this study.
Acknowledgments The authors gratefully acknowledge the financial support provided by the European Commission FP7 Programme for the TLEMsafe project (www.tlemsafe.eu) (Grant no. FP7-ICT-247860).
Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jbiomech.2014.04.030.
References Alexander, E.J., Andriacchi, T.P., 2001. Correcting for deformation in skin-based marker systems. J. Biomech. 34, 355–361. Andersen, M.S., Damsgaard, M., MacWilliams, B., Rasmussen, J., 2010. A computationally efficient optimisation-based method for parameter identification of kinematically determinate and over-determinate biomechanical systems. Comput. Methods Biomech. Biomed. Eng. 13, 171–183. Anderson, F.C., Pandy, M.G., 2001. Static and dynamic optimization solutions for gait are practically equivalent. J. Biomech. 34, 153–161. Andriacchi, T.P., Dyrby, C.O., 2005. Interactions between kinematics and loading during walking for the normal and ACL deficient knee. J. Biomech. 38, 293–298. Audu, M.L., Kirsch, R.F., Triolo, R.J., 2003. A computational technique for determining the ground reaction forces in human bipedal stance. J. Appl. Biomech. 19, 361–371. Audu, M.L., Kirsch, R.F., Triolo, R.J., 2007. Experimental verification of a computational technique for determining ground reactions in human bipedal stance. J. Biomech. 40, 1115–1124. Choi, A., Lee, J.M., Mun, J.H., 2013. Ground reaction forces predicted by using artificial neural network during asymmetric movements. Int. J. Precis. Eng. Manuf. 14, 475–483. Collins, S.H., Adamczyk, P.G., Ferris, D.P., Kuo, A.D., 2009. A simple method for calibrating force plates and force treadmills using an instrumented pole. Gait Posture 29, 59–64. Damsgaard, M., Rasmussen, J., Christensen, S.T., Surma, E., de Zee, M., 2006. Analysis of musculoskeletal systems in the anybody modeling system. Simul. Model. Pract. Theory 14, 1100–1111. Delp, S.L., Anderson, F.C., Arnold, A.S., Loan, P., Habib, A., John, C.T., Guendelman, E., Thelen, D.G., 2007. OpenSim: open-source software to create and analyze dynamic Simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. Fluit, R., Pellikaan, P., Carbone, V., Krogt, M.M.v.d., Damsgaard, M., Koopman, H.F.J.M., verdonschot, N., 2013. A new geometrically consistent musculoskeletal model of the lower extremity based on imaging and cadaver measurements. In: 11th International Symposium of Computer Methods in Biomechanics and Biomedical Engineering, Salt Lake City, Utah, USA. Greenwood, D., 1988. Dynamics of a Rigid Body, Principle of Dynamics. Prentice Hall, Englewood Cliffs NJ, pp. 389–468. Hesse, S., Konrad, M., Uhlenbrock, D., 1999. Treadmill walking with partial body weight support versus floor walking in hemiparetic subjects. Arch. Phys. Med. Rehabil. 80, 421–427. John, C.T., Anderson, F.C., Guendelman, E., Higginson, J.S., Delp, S.L., 2007. LongDuration Muscle-Actuated Simulations of Walking at Multiple Speeds. In: American Society of Biomechanics. Stanford, California, USA. Klein Horsman, M.D., Koopman, H.F., van der Helm, F.C., Prose, L.P., Veeger, H.E., 2007. Morphological muscle and joint parameters for musculoskeletal modelling of the lower extremity. Clin. Biomech. (Bristol, Avon) 22, 239–247. Kuo, A.D., 1998. A least-squares estimation approach to improving the precision of inverse dynamics computations. J. Biomech. Eng.Trans. ASME 120, 148–159. Leardini, A., Chiari, L., Della Croce, U., Cappozzo, A., 2005. Human movement analysis using stereophotogrammetry. Part 3 soft tissue artifact assessment and compensation. Gait Posture 21, 212–225. Luinge, H.J., Veltink, P.H., 2005. Measuring orientation of human body segments using miniature gyroscopes and accelerometers. Med. Biol. Eng. Comput. 43, 273–282. O'Connor, C.M., Thorpe, S.K., O'Malley, M.J., Vaughan, C.L., 2007. Automatic detection of gait events using kinematic data. Gait Posture 25, 469–474. Oh, S.E., Choi, A., Mun, J.H., 2013. Prediction of ground reaction forces during gait based on kinematics and a neural network model. J. Biomech. 46, 2372–2380. Oosterwaal, M., Telfer, S., Torholm, S., Carbes, S., van Rhijn, L.W., Macduff, R., Meijer, K., Woodburn, J., 2011. Generation of subject-specific, dynamic, multisegment ankle and foot models to improve orthotic design: a feasibility study. BMC Musculoskelet. Disord., 12. OpenSim User's Guide, Section Residual Reduction Algorithm. Last edit from 10th of June, 2013. 〈http://simtk-confluence.stanford.edu:8080/display/Open Sim/Residual+Reduction+Algorithm〉. Pamies-Vila, R., Font-Llagunes, J.M., Cuadrado, J., Alonso, F.J., 2012. Analysis of different uncertainties in the inverse dynamic analysis of human gait. Mech. Mach. Theory 58, 153–164. Pearsall, D.J., Costigan, P.A., 1999. The effect of segment parameter error on gait analysis results. Gait Posture 9, 173–183. Rao, G., Amarantini, D., Berton, E., Favier, D., 2006. Influence of body segments' parameters estimation models on inverse dynamics solutions during gait. J. Biomech. 39, 1531–1536. Rasmussen, J., Damsgaard, M., Christensen, S.T., 2000. Inverse-inverse dynamics simulation of musculoskeletal systems. In: European Society of Biomechanics Royal Academy of Medicine in Ireland.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
Ren, L., Jones, R.K., Howard, D., 2008. Whole body inverse dynamics over a complete gait cycle based only on measured kinematics. J. Biomech. 41, 2750–2759. Richards, J.G., 1999. The measurement of human motion: a comparison of commercially available systems. Hum. Mov. Sci. 18, 589–602. Riemer, R., Hsiao-Wecksler, E.T., Zhang, X.D., 2008. Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait. Gait Posture 27, 578–588. Robert, T., Causse, J., Monnier, G., 2013. Estimation of external contact loads using an inverse dynamics and optimization approach: general method and application to sit-to-stand maneuvers. J. Biomech. 46, 2220–2227. Schache, A.G., Baker, R., Vaughan, C.L., 2007. Differences in lower limb transverse plane joint moments during gait when expressed in two alternative reference frames. J. Biomech. 40, 9–19. Scheys, L., Jonkers, I., Schutyser, F., Pans, S., Spaepen, A., Suetens, P., 2005. Image based methods to generate subject-specific musculoskeletal models for gait analysis. Int. Congress Ser. 1281, 62–67. Schmiedmayer, H.B., Kastner, J., 1999. Parameters influencing the accuracy of the point of force application determined with piezoelectric force plates. J. Biomech. 32, 1237–1242. Schwartz, M.H., Rozumalski, A., 2005. A new method for estimating joint parameters from motion data. J. Biomech. 38, 107–116.
2329
Schwer, L.E., 2007. Validation metrics for response histories: perspectives and case studies. Eng. Comput. (Germany) 23, 295–309. Sprague, M.A., Geers, T.L., 2003. Spectral elements and field separation for an acoustic fluid subject to cavitation. J. Comput. Phys. 184, 149–162. Taylor, R., 1990. Interpretation of the correlation coefficient: a basic review. J. Diagn. Med. Sonogr. 6, 35–39. Thelen, D.G., Anderson, F.C., 2006. Using computed muscle control to generate forward dynamic simulations of human walking from experimental data. J. Biomech. 39, 1107–1115. Vicons, 2002. Plug-in-Gait modelling instructions. Vicons Manual, Vicons612 Motion Systems. Oxford Metrics Ltd, Oxford, UK. Vanheule, V., Andersen, M.S., Wirix-Speetjens, R., Jonkers, I., Victor, J., Van den Sloten, J., 2013. Modeling of patient-specific knee kinematics and ligament behavior using force-dependent kinematics. In: XXIV Congress of the International Society of Biomechanics. Natal, Brazil. Wilken, J.M., Rodriguez, K.M., Brawner, M., Darter, B.J., 2012. Reliability and minimal detectible change values for gait kinematics and kinetics in healthy adults. Gait Posture 35, 301–307. Winter, D.A., 2009. Biomechanics and Motor Control of Human Movement, Fourth Edition John Wiley & Sons, Inc., Hoboken, NJ, USA.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Prediction of ground reaction forces and moments during various activities of daily living R. Fluit a,n,1, M.S. Andersen b, S. Kolk c, N. Verdonschot a,d, H.F.J.M. Koopman a a
Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente, Enschede, The Netherlands Department of Mechanical and Manufacturing Engineering, Aalborg University, Aalborg, Denmark c Radboud University Medical Centre, Radboud Institute for Health Sciences, Department of Rehabilitation, Nijmegen, The Netherlands d Radboud University Medical Centre, Radboud Institute for Health Sciences, Orthopaedic Research Laboratory, Nijmegen, The Netherlands b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 18 April 2014
Inverse dynamics based simulations on musculoskeletal models is a commonly used method for the analysis of human movement. Due to inaccuracies in the kinematic and force plate data, and a mismatch between the model and the subject, the equations of motion are violated when solving the inverse dynamics problem. As a result, dynamic inconsistency will exist and lead to residual forces and moments. In this study, we present and evaluate a computational method to perform inverse dynamics-based simulations without force plates, which both improves the dynamic consistency as well as removes the model's dependency on measured external forces. Using the equations of motion and a scaled musculoskeletal model, the ground reaction forces and moments (GRF&Ms) are derived from three-dimensional full-body motion. The method entails a dynamic contact model and optimization techniques to solve the indeterminacy problem during a double contact phase and, in contrast to previously proposed techniques, does not require training or empirical data. The method was applied to nine healthy subjects performing several Activities of Daily Living (ADLs) and evaluated with simultaneously measured force plate data. Except for the transverse ground reaction moment, no significant differences (P40.05) were found between the mean predicted and measured GRF&Ms for almost all ADLs. The mean residual forces and moments, however, were significantly reduced (P40.05) in almost all ADLs using our method compared to conventional inverse dynamic simulations. Hence, the proposed method may be used instead of raw force plate data in human movement analysis using inverse dynamics. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Ground reaction forces and moments Musculoskeletal model Inverse dynamics Dynamic consistency Activities of daily living
1. Introduction Inverse dynamics based simulations on musculoskeletal models are a commonly used method for the analysis of human movement. Despite widespread use, it is well known that solutions obtained with inverse dynamics are sensitive to inaccuracies in the various input variables (Pamies-Vila et al., 2012; Riemer et al., 2008). Errors can stem from estimating body segment parameters (Pearsall and Costigan, 1999; Rao et al., 2006), estimating joint parameters (Schwartz and Rozumalski, 2005), skin movement artifacts (Leardini et al., 2005), noise on skin-mounted marker data (Richards, 1999), estimating the center of pressure (Schmiedmayer and Kastner, 1999) or force plate calibration (Collins et al., 2009). Consequently, when solving the inverse dynamics problem, the equations of motion are violated,
n Correspondence to: Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. Tel.: þ31 53 489 4362; fax: þ31 53 489 2287. E-mail address: r.fl[email protected] (R. Fluit). 1 Visiting address: Campus University, Horstring W213.
http://dx.doi.org/10.1016/j.jbiomech.2014.04.030 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
resulting in dynamic inconsistency, a condition with residual forces and moments (Kuo, 1998). Several algorithms have been proposed that reduce or eliminate the residual forces and moments, such as the least-squared optimization (Kuo, 1998), the residual elimination algorithm (REA) (Thelen and Anderson, 2006) and the residual reduction algorithm (RRA) (Delp et al., 2007). These algorithms adjust the kinematics, ground reaction forces (GRFs) and/or body segment parameters, thereby improving the dynamic consistency. Unfortunately, they have shortcomings too: the REA was shown to dramatically change torso angles for movements longer than 0.5 s (John et al., 2007). For the RRA, an adjustment in the joint angles of up to five degrees is considered reasonable (OpenSim User's Guide). Since these differences are larger than the minimal detectable change for most of the joint angles (Wilken et al., 2012), these adjustments may not be defendable. Alternatively, dynamic consistency can be improved by deriving the ground reaction forces and moments (GRF&Ms) from three-dimensional full-body motion using the equations of motion (Audu et al., 2007; Choi et al., 2013; Oh et al., 2013; Ren et al., 2008; Robert et al., 2013). An additional advantage of this method
2322
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
is that it enables inverse dynamic analysis for studies without force plate data, for example ambulatory measurements with inertial measurements only (e.g. Luinge and Veltink, 2005) or motion capture during treadmill walking (e.g. Hesse et al., 1999). A difficulty of this method is the indeterminacy problem during a double contact phase when the system defines a closed kinetic chain. To overcome this problem, Audu et al. (2003) used optimization techniques to compute the GRF&Ms for different static postures of a standing bipedal model, which were later validated against measured data (Audu et al., 2007). However, it is unknown whether this method is valid for dynamic movements. Ren et al. (2008) introduced a smooth transition assumption to solve the indeterminacy problem. The smoothing functions were based on empirical data and, therefore, the smooth transition assumption may not be applicable for movements other than those present in the empirical data. Oh et al. (2013) and Choi et al. (2013) solved the indeterminacy problem using an artificial neural network. Also their method requires training data, which is not always present. Robert et al. (2013) tested several optimization methods to predict the external contact loads during sit-to-stand movements. Although their method does not require empirical or training data, their contact configurations are simplified and the method was validated for sit-to-stand motion only. Therefore, the purpose of this paper is to demonstrate a universal method for predicting the GRF&Ms based on measured kinematic data only, which is applicable to a variety of Activities of Daily Living (ADLs), and in which the indeterminacy problem during the double contact phase of gait and gait-related ADLs is solved without the use of empirical or training data. The predicted GRF&Ms were evaluated with simultaneously measured force plate data. For the trials where subjects walked at self-selected comfortable walking speed, a sensitivity study was performed to evaluate the effects of the chosen muscle recruitment strategy and parameters of the ground contact model on the accuracy of the predictions.
The local ethics committee approved the study protocol and informed consent was obtained from all subjects prior to the study. 2.2. Instrumentation A six-camera digital optical motion capture system (Vicon MX, Oxford, UK) was used to capture 53 retro-reflective markers placed on the body at 100 Hz sampling rate. Markers were placed according to the Vicons Plug-in-Gait model (Vicons, 2002) with additional markers on the anterior side of the thigh and shank at 1/3 and 2/3 segment length, the medial femoral epicondyle, the medial malleolus and on the fifth metatarsal head of the foot. Two custom-built force plates (AMTI, Watertown, MA, USA), embedded level in the laboratory floor, measured GRF&Ms at 1000 Hz. For the stair negotiation trials, a custom-built four step staircase was used with a riser height of 180 mm and a thread depth of 250 mm. For these trials, GRF&Ms were collected for the initial or last step at floor level, the first step using a box on a force plate and the second step using an additional custom-built portable force plate (ForceLink B.V., Culemborg, the Netherlands). 2.3. Experimental protocol Three-dimensional motion capture with synchronized force plate recordings of several ADLs were collected: level walking at comfortable (CWS), slow (CWS 30%) and fast (CWS þ 30%) walking speed, walking over a 10 (OBS10), 20 (OBS20) and 30 (OBS30) cm obstacle, gait initiation (GI) and termination (GT) starting or ending with the dominant (DL) or non-dominant leg (NDL), deep squatting (DS) and stair ascent (SA) and descent (SD). The dominant leg was determined by asking the participant with which leg they would kick a ball. All ADLs were performed barefooted. For the stair negotiation trials, subjects were instructed to perform the activity without the use of handrails. For the obstacle trials, subjects were allowed to choose their own strategy to step over the obstacle placed on the walkway, as long as the force plate was hit clean. For each ADL, subjects performed several practice trials before the actual trials were recorded. 2.4. Data processing
2. Methods
3D position data of the markers and the force plate data were low-pass filtered using a second order, zero-phase Butterworth filter with a cut-off frequency of 12 and 15 Hz, respectively. For the CWS, CWS 30% and CWS þ 30% trials, three experimental data sets per subject were used for validating the predictions of the GRF&Ms. For all other ADLs, one experimental data set per subject was used, yielding a total of 261 trials. In total, 257 trials were recorded correctly, of which another four trials were excluded due to missing markers.
2.1. Subjects
2.5. Prediction of GRF&Ms
Nine healthy subjects (4 males and 5 females; age: 41.6 7 15.9 yr; height: 1.74 70.12 m; weight: 73.0 7 11.1 kg, Body Mass Index: 23.9 7 2.0 kg/m2) with no history of musculoskeletal disorders volunteered for the study at the Rehabilitation Department of the Radboud University Medical Centre, Nijmegen, the Netherlands.
For all ADLs, the only significant external forces and moments were the GRF&Ms. These can be estimated using the Newton–Euler equations of motion, which state that the sum of all external forces balances the sum of the massacceleration products of all individual body segments (Greenwood, 1988). Hence,
Fig. 1. Left: Visualization of the 12 contact points for each foot. Points were defined at the medial and lateral side of the heel, at the base of the first and fifth metatarsal bone, at the head of each metatarsal bone and at the big, second and fifth toe. Right: Visualization of the characteristic strength function cp;i (Eq. (1)) and the smoothed strength function (Eq. (2)) for the right heel contact node during a gait cycle. Heel contact and toe-off of the right leg are abbreviated as HCR and TOR respectively and, analogously, for the left leg as HCL and TOL.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329 during the single support phase of gait, the GRF&Ms can be obtained directly from the Newton–Euler equations when full-body kinematics is available. During a double contact phase; however, the problem of determining the GRF&Ms under each foot becomes underdetermined and they cannot be resolved from the Newton–Euler equations alone. To overcome this, the computation of the GRF&Ms was made part of the muscle recruitment algorithm by introducing artificial muscle-like actuators at 12 contact points under each foot (Fig. 1, left side). At each contact point p, five artificial muscle-like actuators were added. One actuator was aligned with the vertical direction of the force plate (z) and was able to generate a normal force. The other two pairs of actuators were aligned with the medio-lateral (x) and antero–posterior (y) direction of the force plate and were able to generate positive or negative static friction forces. The force exerted by each actuator was defined as qF max , with q the activation level and F max the strength. A dynamic contact model was used to define the strength profile of each actuator, ensuring they were able to generate a reaction force only if their contact point p was close to the floor and almost without motion. For each contact point p, a characteristic strength function cp;i was defined ( F max if pz o zcrit and jjpjj o vcrit ð1Þ cp;i ¼ 0 otherwise where F max ¼ 0:4 BW, zcrit ¼ 0:03 m (level of the floor was set to z¼ 0), vcrit ¼ 0:8ms and i the index used for the step of the simulation (each step was 0.01 s). These values were based on pilot tests and a sensitivity analysis of the implications of these values was performed (Section 2.8). For the transitions of cp;i at the discrete
2323
time step i ¼t, the characteristic strength function was smoothed (Fig. 1, right side) according to 1 1 cp;i ¼ ðt 3…t þ 8Þ ¼ F max 1 7 cos πn with n ¼ 1…12 ð2Þ 2 12 The smoothing function was based on pilot tests and was introduced to prevent discontinuities in the estimated GRF&Ms during the transitions from being completely determined (single stance) to being underdetermined (double stance) and vice versa. The activation level q of each artificial muscle-like actuators was solved as part of the muscle recruitment problem, thereby determining the magnitude of the GRF&Ms. The solver did not distinguish between single and double contact phases.
2.6. Inverse dynamic model All simulations were performed using a 28 degrees-of-freedom (DOF) full body model (Fig. 2) as available in the AnyBody Modeling System (version 5.3.1, AnyBody Technology A/S, Aalborg, Denmark) (Damsgaard et al., 2006). For the lower limbs the geometry of the new Twente Lower Extremity Model was used (Fluit et al., 2013). The segment masses were linearly scaled using data from Winter (2009). The model marker positions and segment lengths were determined using a parameter optimization algorithm (Andersen et al., 2010) applied to an arbitrary CWS trial. Subsequently, inverse kinematics for all other trials were solved using the optimized parameters. Muscles were included for the lower limbs (Klein Horsman et al., 2007) whereas for the upper extremities, ideal joint torque generators were used. Each leg consisted of 55 muscle-tendon parts described by 166 Hill-type muscle elements. To obtain the muscle forces and thus the predicted GRF&Ms, the muscle recruitment problem was solved through static optimization, minimizing the sum of the cubed muscle activations. To improve numerical stability, additional muscle-like actuators were attached at the origin of the pelvis segment, which could generate small residual forces and moments up to 10 N or Nm. The activation level of these muscles was solved as part of the muscle recruitment optimization problem, so as to minimize their contribution. Skeletal muscles and muscle-like actuators for the GRF&Ms and residuals were weighted equally. However, the strength of the actuators for the GRF&Ms was high compared to the skeletal muscles whereas the actuators for the residuals had a low strength.
2.7. Accuracy analysis
Fig. 2. Visualization of the skeletal model and the 28 degrees of freedom (DOFs).
To assess the quality of the predictions, the predicted GRF&Ms were compared with the force plate data. Except for the DS trials, the timespan over which the GRF&Ms were compared was defined as the time for which the measured vertical GRF was larger than 5% of its maximum. For the DS trials, this timespan was defined from 0.5 s before until 0.5 s after the knee of the dominant leg was flexed at 151 flexion angle. The ground reaction moments (GRMs) were calculated about the point on the force plate surface vertically below the ankle joint. Following the definitions of Ren et al. (2008), the differences were quantified using the Root Mean Squared Difference (RMSD) and, additionally for the CWS and CWS þ30% trials, the relative RMSD (rRMSD). The Pearson correlation coefficient ρ was calculated as well and categorized (in absolute value) as ρ r 0:35, 0:35 o ρ r 0:67, 0:67o ρ r 0:9, 0:9 o ρ to be weak, moderate, strong or excellent correlations respectively (following Taylor, 1990).
Table 1 Comparison of measured and predicted GRFs using the curve magnitude (M) and phase (P) difference (mean (SD)) and Pearson correlation coefficient ρ. Results are averaged over both legs and all subjects. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to the two-tailed Wilcoxon Signed Rank Test (α¼ 0.05). Activity
CWSþ30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
Vertical GRF
Antero-posterior GRF
P (%)
M (%)
3.0 2.4 2.7 2.4 2.3 2.5 3.2 2.3 2.5 2.7 2.7 3.1 4.7
3.4 2.9 3.0 3.4 2.8 2.1 2.0 2.2 2.1 1.4 2.5 4.0 1.2
(0.4) (0.6) (0.4) (0.7) (0.2) (0.6) (0.9) (0.3) (0.2) (0.5) (0.6) (1.3) (0.6)
ρ () (1.5) (0.8) (0.9) (2.2) (2.5) (1.7) (1.1) (0.5) (1.2) (1.0) (1.4) (1.3) (0.9)
0.954 0.967 0.957 0.976 0.980 0.976 0.947 0.962 0.957 0.943 0.621 0.941 0.895
P (%)
n† n†
n†
n† n†
9.0 8.4 9.3 6.0 8.1 6.3 9.5 7.5 8.0 7.9 35.0 26.0 22.3
(1.3) (1.9) (1.2) (1.1) (2.2) (1.5) (1.2) (1.2) (1.3) (1.3) (15.1) (7.6) (2.8)
Medio-lateral GRF
M (%)
ρ ()
P (%)
4.1 26.4 15.8 1.1 13.4 3.3 11.0 6.4 6.1 5.1 71.9 57.5 15.1
0.948 0.965 0.957 0.968 0.944 0.962 0.902 0.966 0.959 0.969 0.202 0.436 0.637
11.8 7.6 9.8 16.7 11.9 17.8 15.2 9.8 11.1 12.4 34.7 23.2 26.2
(8.3) (9.1) (9.3) (11.5) (7.4) (5.8) (6.6) (6.1) (11.3) (7.3) (139.6) (65.3) (12.6)
n†
n†
ρ ()
M (%) (2.6) (1.9) (2.3) (12.4) (2.2) (6.9) (4.9) (3.0) (2.8) (3.5) (14.0) (4.5) (2.1)
1.0 1.9 1.9 8.6 5.1 4.7 6.1 8.6 12.0 14.1 11.2 9.9 10.3
(10.8) (4.0) (6.9) (9.6) (5.5) (17.7) (7.9) (8.7) (12.4) (15.4) (43.8) (14.7) (12.4)
0.753 0.866 0.818 0.898 0.842 0.863 0.726 0.836 0.810 0.713 0.261 0.609 0.455
2324
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
Further, differences were quantified using the Spraque and Geers curve magnitude (M) and phase (P) difference (Schwer, 2007; Sprague and Geers, 2003). The metrics M and P are expressed in percentages and are designed to give zero when the curves are identical. Additionally, two-tailed Wilcoxon signed rank tests (α ¼ 0.05) were performed to test if the medians of the absolute mean and absolute maximum of the measured GRF&Ms were significantly different than those of the predicted GRF&Ms. Similar comparisons were performed for a subset of the computed joint moments. The absolute mean and maximum of the residual forces and moments were calculated for each trial for the timespan in which ground contact was fully defined. One tailed Wilcoxon signed rank tests (α ¼ 0.05) were performed to test if the medians of the absolute mean and absolute maximum of the residual forces and moments using the predicted GRF&Ms were significantly lower than those using the measured GRF&Ms.
The accuracy of the dynamic contact model (Eqs. (1) and (2)) was assessed by comparing heel-strike (HS) and toe-off (TO) events for the CWS trials, defined as the time at which the vertical GRF reached 5% of the maximum.
2.8. Sensitivity analysis For all CWS trials, a sensitivity analysis of the parameters of the dynamic contact model (F max , zcrit and vcrit ) and the muscle recruitment criterion on the prediction accuracy was made. For F max , zcrit and vcrit , simulations were performed using values at 70% and 130% of their original magnitude (Eq. (1)). The sensitivity was quantified as the change in M, P, ρ and rRMSD, where a negative change indicated an improvement in the predicted GRF&Ms and vice versa.
Fig. 3. Calculated GRFs, normalized to body weight (BW), and GRMs, normalized to BW and body height (BH) ( 7 1 SD around mean (shaded area)) compared with force plate data (mean (thick line) 71 SD (thin lines)), averaged over all subjects. For the CWS, OBS30 and SD trials, the GRF&Ms of both legs are shown and the heel contact (HC) and toe-off (TO) events are indicated. For the DS trials, only the GRF&Ms of the dominant leg are shown and at the start (S) and end (E) of the squat movement the knee of the dominant leg was flexed at 15 1 flexion angle.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
3. Results The model showed excellent predictions for the vertical GRF (ρ ranging from 0:621 0:980, median 0:957) and the anteroposterior GRF (ρ ranging from 0:202 0:969, median 0:957) for almost all activities (Table 1, Fig. 3). The magnitude of the vertical GRF was slightly but consistently underestimated (M ranging from 1.2% to 4.0%, median 2:5%), whereas the magnitude of the antero-posterior GRF was consistently overestimated (M ranging from 1:1% to 71:9% median 11:0%). Nevertheless, no significant difference (P 4 0:05) was found between the medians of the absolute mean measured GRFs and those of the predicted GRFs for all activities. Strong correlations were obtained for the medio-lateral GRF for almost all activities (ρ ranging from 0:261 to 0:898, median 0:810). In general, the most extreme activities (DS, SA and SD) showed the weakest correlations, for instance the antero-posterior (ρ ¼ 0:202) and mediolateral (ρ ¼ 0:261) GRFs for the DS trials. The model showed strong correlations for the sagittal GRM (ρ ranging from 0:506 to 0:922, median 0:789) and frontal GRM (ρ ranging from 0:199 to 0:801, median 0:668) for most of the activities (Table 2). The transverse GRM showed the largest differences (ρ ranging from 0:155 to 0:782, median 0:588). A significant difference (P o0:05) was found between the medians of the absolute mean measured transverse GRM and those of the predicted transverse GRM for all activities.
2325
The RMSDs between measured and predicted GRF&Ms were similar across different ADLs (Table 3). Regarding the joint moments, strong correlations were found for almost all activities for hip flexion (ρ ranging from 0:521 to 0:972, median 0:776), hip abduction (ρ ranging from 0.257 to 0.946, median 0:880), hip external rotation (ρ ranging from 0.809 to 0.943, median 0:871), knee flexion (ρ ranging from 0:548 to 0:911, median 0:737) and ankle dorsiflexion (ρ ranging from 0:647 to 0:927, median 0:820). A weak correlation was found for the hip abduction moment for the DS trials only (ρ ¼ 0:258, see Table A1 and Fig. 4). For a subset of activities, Fig. 4 demonstrates the joint moments at the DOFs in the dominant leg only using measured and predicted GRF&Ms (Table A1 contains a complete overview for all activities and both legs). Using our method, the medians of the absolute mean residual forces, transverse and frontal moments were significantly reduced compared to those obtained with a conventional inverse dynamics simulation for almost all activities (P o 0:05). Averaging over all activities and subjects, the absolute mean residual forces and moments applied at the pelvis were about 2 3 N or Nm using the predicted GRF&Ms, thereby improving dynamic consistency (Table A2). The vertical residual force showed the strongest reduction of 90%. Only the residual sagittal moment showed a non-significant reduction of 33%, probably because this moment was already small using traditional inverse dynamics (4:3 Nm, averaging over all activities and subjects).
Table 2 Comparison of measured and predicted GRMs using curve magnitude (M) and phase (P) difference (mean (SD)) and Pearson correlation coefficient ρ. GRMs are expressed at the virtual point on the force plate surface below the ankle joint. Results are averaged over both legs and all subjects. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (α¼ 0.05). Activity
Sagittal GRM P (%)
CWS þ 30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
†
n
6.8 9.6 7.4 9.3 11.4 8.5 10.4 7.1 7.7 8.7 15.6 12.4 8.3
Frontal GRM M (%)
(2.3) (2.3) (1.8) (2.9) (4.7) (2.5) (3.0) (1.8) (2.6) (1.7) (9.8) (3.3) (3.1)
6.6 15.2 8.8 0.4 4.8 14.6 3.7 6.3 6.4 4.4 50.5 14.0 5.5
(28.2) (6.1) (6.3) (7.1) (23.2) (12.1) (10.9) (7.2) (6.0) (7.4) (59.0) (14.7) (9.1)
ρ ()
P (%)
0.929 0.872 0.922 0.727 0.878 0.680 0.777 0.915 0.863 0.801 0.506 0.647 0.830
16.4 10.7 13.4 8.4 10.6 13.0 14.0 12.7 14.6 15.3 24.9 9.7 9.7
Transverse GRM ρ ()
M (%) (7.3) (3.4) (5.8) (2.3) (3.9) (4.8) (4.8) (4.6) (4.3) (6.7) (6.7) (2.0) (2.9)
10.5 7.3 10.8 9.2 7.2 0.7 0.6 3.7 8.3 15.3 37.1 17.4 15.7
(15.6) (11.2) (10.0) (15.3) (9.1) (18.8) (27.9) (20.5) (21.7) (28.6) (23.2) (13.7) (11.7)
0.540 0.777 0.684 0.703 0.801 0.488 0.539 0.694 0.651 0.641 0.199 0.781 0.781
P (%) n† n† n† n† n† n† n† n† n† n† n† n† n†
33.1 46.0 36.8 26.8 30.9 68.7 56.1 33.9 33.2 33.2 42.5 53.9 67.0
ρ ()
M (%) (7.5) (9.5) (6.3) (10.9) (7.6) (10.4) (8.6) (7.2) (8.4) (4.4) (7.2) (8.6) (5.6)
56.7 46.5 48.3 33.3 47.6 79.4 52.0 46.2 52.0 50.6 91.6 33.8 61.0
(20.6) (19.4) (14.6) (27.6) (18.2) (3.4) (41.1) (13.0) (11.7) (13.7) (147.1) (22.5) (15.4)
0.762 0.598 0.704 0.588 0.525 0.128 0.157 0.782 0.759 0.657 0.230 0.155 0.155
Table 3 Absolute differences (mean (SD)) between the measured and predicted GRF&Ms, quantified using the Root Mean Squared Difference (RMSD) following the definitions of Ren et al. (2008). Activity CWS þ 30% CWS 30% CWS GI (NDL) GI (DL) GT (NDL) GT (DL) OBS10 OBS20 OBS30 DS SA SD
Vertical GRF (N/kg) 0.85 (0.17) 0.64 (0.15) 0.74 (0.13) 0.54 (0.18) 0.63 (0.10) 0.52 (0.12) 0.91 (0.60) 0.64 (0.17) 0.68 (0.12) 0.74 (0.24) 0.81 (0.36) 0.95 (0.25) 0.34 (0.12)
Antero-posterior GRF (N/kg) 0.43 (0.06) 0.30 (0.05) 0.38 (0.07) 0.16 (0.04) 0.29 (0.07) 0.18 (0.06) 0.32 (0.08) 0.32 (0.06) 0.34 (0.07) 0.33 (0.07) 0.30 (0.22) 0.21 (0.06) 0.22 (0.07)
Medio-lateral GRF (N/kg) 0.23 (0.07) 0.12 (0.03) 0.17 (0.04) 0.19 (0.07) 0.17 (0.03) 0.23 (0.07) 0.21 (0.06) 0.18 (0.05) 0.20 (0.07) 0.22 (0.07) 0.24 (0.16) 0.22 (0.04) 0.42 (0.17)
Sagittal GRF (Nm/kg) 0.17 (0.07) 0.22 (0.07) 0.18 (0.05) 0.11 (0.06) 0.19 (0.05) 0.16 (0.06) 0.20 (0.06) 0.23 (0.11) 0.22 (0.08) 0.23 (0.08) 0.22 (0.05) 0.21 (0.05) 0.09 (0.02)
Frontal GRM (Nm/kg) 0.13 (0.03) 0.08 (0.02) 0.11 (0.02) 0.05 (0.03) 0.07 (0.02) 0.06 (0.02) 0.09 (0.03) 0.12 (0.09) 0.11 (0.03) 0.12 (0.04) 0.17 (0.07) 0.14 (0.05) 0.11 (0.05)
Transverse GRM (Nm/kg) 0.28 (0.08) 0.21 (0.07) 0.22 (0.06) 0.08 (0.06) 0.17 (0.06) 0.16 (0.06) 0.24 (0.10) 0.23 (0.06) 0.22 (0.06) 0.21 (0.07) 0.10 (0.03) 0.19 (0.06) 0.13 (0.04)
2326
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
The sensitivity of the recruitment criterion on the values of P and ρ was small and negligible. With respect to M and rRMSD, none of the recruitment criterions performed better than the others (Table 4). The sensitivity of F max , zcrit and vcrit on the values of P and ρ was also small and negligible. vcrit had the largest effects on M and rRMSD, yielding clearly better predictions for the antero-posterio GRF and sagittal GRM for the 30% case (Table 5). For the CWS and CWS þ30% trials, the RMSD and rRMSD showed similar trends as previous literature (Oh et al., 2013; Ren et al., 2008) (Table 6). Our model provided more accurate predictions for the medio-lateral GRF and frontal GRM than Ren et al. (2008) for the CWS trials, although even better results were obtained by Oh et al. (2013). The dynamic contact model was able to predict heel strike (HS) and toe-off (TO) with reasonable accuracy. Averaging over the CWS trials, HS was predicted consistently early with an error of 28 7 13 ms, and TO was predicted consistently late with an error of 16 77 ms.
4. Discussion In this paper, we demonstrated a universal method to predict the GRF&Ms using kinematic data and a scaled musculoskeletal model only, applied to a variety of ADLs. The method entails a dynamic contact model and optimization techniques to solve the indeterminacy problem during a double contact phase. In general, reasonably good results were obtained for all analyzed activities. However, weak or even negative correlations were found for the transverse GRM for multiple activities. Averaging over all subjects and activities, the magnitude of the transverse GRM was underpredicted with about 50%. These differences are likely caused by the hinged knee, which does not allow for tibial rotation. Previous research has reported a range of 16.21 for this DOF during normal walking, in which the internalexternal knee moment was in the same direction as the corresponding rotation (Andriacchi and Dyrby, 2005). Since for normal gait the internal–external knee moment is primarily determined by the transverse GRM (Schache et al., 2007), omitting this DOF
Fig. 4. Computed joint moments, normalized to body weight (BW) and body height (BH), using the predicted GRF&Ms ( 7 1 SD about mean (shaded area)) compared to those obtained using the measured GRF&Ms (mean (thick line) 7 1 SD (thin lines)), averaged over all subjects. For the CWS, OBS30 and SD trials, the joint moments of the dominant leg are shown from heel contact (HC) until toe-off (TO). For the DS trials, the joint moments of the dominant leg are shown and at the start (S) and end (E) of the squat movement the knee of the dominant leg was flexed at 15 1 flexion angle. n/† indicates a significant difference between the medians of the absolute maximum (n)/absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (α¼ 0.05).
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
2327
(Ren et al., 2008). Our sensitivity analysis on the contact model showed that consistent over- and underestimation can be reduced, indicating that our method is not necessarily inferior to traditional inverse dynamics in clinical analysis of gait. Compared to previous literature (Table 6), our predictions are more accurate than Ren et al. (2008) for the medio-lateral GRF (rRMSD of 14:9% and 16:6% versus 20:0% and 24:1% for the CWS and CWSþ30% trials, respectively) as well as the frontal GRM (rRMSD of 22:9% and 27:1% versus 32:5% and 37:7% for the CWS and CWSþ30% trials, respectively). When comparing our results with Oh et al. (2013) for the CWS trials, we obtained a much higher rRMSD for the transverse GRM (40:6% versus 25:5%) whereas the rRMSD for the other elements of the GRF&Ms was only roughly 2% higher. Here it should be noted that the method of Oh et al. (2013) involved an algorithm that required training data to relate
will result in smaller magnitudes of the transverse GRM (see Fig. 3 for the CWS activity). Furthermore, weak correlations were found for the antero-posterior, medio-lateral GRF and frontal GRM for the DS trials. A likely explanation is that for this task, the magnitude of these forces and moments are low, in which case noise has a larger influence on the correlation. Additionally, for the DS task, the total external forces and moments need to be distributed over both feet during the complete movement, which may have increased errors. The joint moments computed by our method were similar to those using traditional inverse dynamics. Although our method underestimated the vertical GRF and overestimated the anteroposterior GRF for all activities and thus may have induced errors in the joint moments, it is likely that the joint moments computed with traditional inverse dynamics suffer from similar errors
Table 4 Sensitivity of recruitment criterion expressed as the change in curve magnitude difference M and change in relative Root Mean Squared Difference (rRMSD). Changes are expressed with respect to the results obtained with the simulations using P ¼3 as recruitment criterion. A negative change indicates an improvement in the predicted GRF&Ms and vice versa. The Min/Max Aux criterion is a composite recruitment strategy consisting of two terms: a term which minimizes the maximum muscle stress and an auxiliary quadratic term. P ¼1
Vertical GRF Antero-posterior GRF Medio-lateral GRF Sagittal GRM Frontal GRM Transverse GRM
P¼2
P¼4
Min/Max Aux
M (%)
rRMSD (%)
M (%)
rRMSD (%)
M (%)
rRMSD (%)
M (%)
rRMSD (%)
0.3 5.0 3.7 2.9 0.7 2.6
0.2 0.8 0.7 1.0 2.1 0.0
0.1 2.0 0.9 0.1 0.3 2.9
0.1 0.3 0.6 0.7 0.9 0.9
0.0 1.7 0.7 1.2 0.6 1.4
0.0 0.1 0.2 0.3 0.1 0.1
0.2 1.4 0.0 0.4 0.9 1.1
0.1 0.3 0.1 1.0 0.2 0.6
Table 5 Sensitivity of the parameters of the ground contact model (as defined in Eq. (1)) expressed as the change in curve magnitude difference M and change in relative Root Mean Squared Difference (rRMSD). A negative change indicates an improvement in the predicted GRF&Ms and vice versa. F max
vcrit
30%
Vertical GRF Antero-posterior GRF Medio-lateral GRF Sagittal GRM Frontal GRM Transverse GRM
þ30%
zcrit
30%
þ 30%
30%
þ 30%
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
M (%)
rRMSE (%)
0.5 1.7 0.1 1.1 2.3 1.2
0.2 0.3 0.2 0.6 0.0 0.1
0.3 1.6 0.1 1.0 1.7 1.0
0.2 0.3 0.1 0.6 0.1 0.2
0.4 5.5 0.5 4.1 1.5 1.5
0.4 0.9 0.2 1.7 0.3 1.2
0.4 3.2 0.6 3.7 1.0 0.6
0.5 0.4 0.2 1.5 0.2 0.6
0.2 1.9 0.0 0.8 0.4 0.9
0.1 0.4 0.1 0.9 0.2 0.1
0.2 1.3 0.0 0.6 0.4 0.7
0.1 0.3 0.1 0.7 0.2 0.1
Table 6 Differences between measured and predicted GRF&Ms for the CWS and CWSþ 30% trials, quantified using the Root Mean Squared Difference (RMSD) and relative RMSD (rRMSD), following the definitions of Ren et al. (2008). Oh et al. (2013) did not investigate CWS þ30% trials. Ren et al. (2008) investigated fast walking trials, where subjects walked about 27% faster than their normal walking speed, which is similar to our CWS þ 30% trials. Method Participants
Smooth transition assumption (Results from Ren et al., 2008) N¼ 3 CWS þ 30%
CWS RMSD (N/kg or Nm/kg)
Artificial neural network Optimization approach (Results from Oh et al., 2013) (Method proposed in this study) N ¼5 N ¼9
rRMSD (%) RMSD (N/kg or Nm/kg)
Vertical GRF 0.710 (0.190) 5.6 (1.5) AP GRF 0.473 (0.068) 10.9 (0.83) ML GRF 0.191 (0.034) 20.0 (2.7) Sagittal GRM 0.199 (0.106) 12.2 (4.8) Frontal GRM 0.148 (0.013) 32.5 (4.3) Transverse GRM 0.039 (0.015) 26.2 (9.4)
1.012 0.507 0.335 0.225 0.235 0.075
(0.223) (0.085) (0.012) (0.029) (0.018) (0.009)
CWS
CWS þ 30%
CWS
rRMSD (%)
RMSD (N/kg or Nm/kg)
rRMSD (%)
RMSD (N/kg rRMSD (%) RMSD (N/kg or Nm/kg) rRMSD (%) or Nm/kg)
6.6 8.9 24.1 14.2 37.7 41.7
0.649 0.154 0.040 0.081 0.052 0.032
5.8 7.3 10.9 9.9 22.8 25.5
0.74 0.38 0.17 0.18 0.11 0.22
( 1.6) (1.3) (1.6) (1.4) (1.0) (1.8)
(0.182) (0.057) (0.022) (0.045) (0.029) (0.018)
(1.0) (0.8) (1.8) (1.9) (4.9) (4.5)
(0.13) (0.07) (0.04) (0.05) (0.02) (0.06)
6.6 9.3 14.9 12.4 22.9 40.6
(1.1) (2.0) (3.4) (3.5) (5.9) (11.3)
0.85 0.43 0.23 0.17 0.13 0.28
(0.17) (0.06) (0.07) (0.07) (0.03) (0.08)
6.9 8.5 16.6 10.4 27.1 38.4
(1.3) (1.6) (4.6) (3.7) (9.0) (10.9)
2328
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
kinematic measures with the GRF&Ms during the double contact phase. In contrast, our method can be used directly without any training data, since optimization techniques were used to solve the indeterminacy problem. Our dynamic contact model was able to predict HS and TO with an error of 28 7 13 ms and 16 7 7 ms respectively for the CWS trials. Oh et al. (2013) used the dynamic contact model of O'Connor et al. (2007), which predicted HS and TO with errors of 16 7 15 ms and 9 7 15 ms, respectively. Our contact model showed larger errors compared to O'Connor et al. (2007), probably because we did not use optimal values for vcrit and zcrit . For example, when vcrit was reduced with 30%, the errors reduced to 22 7 10 ms and 13 7 7 ms for HS and TO, respectively. The magnitude of these errors is similar to O'Connor et al. (2007). Several limitation of this work should be noted. First, our method requires full body kinematics and a corresponding musculoskeletal model, which is not always available. However, Robert et al. (2013) previously showed that a similar optimization approach at the level of joint moments for sit-to-stand motions provided reasonable predictions as well. Second, the knee was modeled as a hinge, which affected the transverse GRM, and the foot was modeled as a single segment, which may have affected the contact dynamics. Third, it has been reported that soft tissue artifacts (STA) are the most significant source of error in human movement analysis (Leardini et al., 2005; Alexander and Andriacchi, 2001). STA may have affected our joint kinematics, predominantly in the non-sagittal plane, and our estimation of the segment lengths. It is likely that these errors have propagated in the predicted GRF&Ms. Fourth, for the skeletal muscles, excitationactivation dynamics were not modeled. Since the actuators for the GRF&Ms, together with the skeletal muscles, were part of the recruitment criterion, inclusion of excitation-activation dynamics might have altered the predicted GRF&Ms. We believe that this may only be the case for activities other than normal walking, since Anderson and Pandy (2001) showed that static and dynamic optimization are practically equivalent for normal walking. Considering the above mentioned possible sources of error, our method could be improved by including STA compensation (Leardini et al., 2005; Alexander and Andriacchi, 2001), imagebased subject-specific geometry (Scheys et al. 2005) and a more advanced and realistic knee joint (Vanheule et al., 2013) or foot model (Oosterwaal et al., 2011). We believe that techniques for creating subject-specific geometry and the development of more advanced joint models should receive priority in future research. In conclusion, reasonably good estimates of GRF&Ms and joint moments for a variety of ADLs can be obtained solely from whole body kinematics, i.e. without the use of training or empirical data. This work may be useful in the development of ambulatory measurement systems using inertial measurement units only (Luinge and Veltink, 2005), motion capture during treadmill walking (e.g. Hesse et al., 1999) or in the field of movement predictions using musculoskeletal optimization (Rasmussen et al., 2000).
Conflict of interest None of the authors have any financial or personal conflict of interest with regard to this study.
Acknowledgments The authors gratefully acknowledge the financial support provided by the European Commission FP7 Programme for the TLEMsafe project (www.tlemsafe.eu) (Grant no. FP7-ICT-247860).
Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jbiomech.2014.04.030.
References Alexander, E.J., Andriacchi, T.P., 2001. Correcting for deformation in skin-based marker systems. J. Biomech. 34, 355–361. Andersen, M.S., Damsgaard, M., MacWilliams, B., Rasmussen, J., 2010. A computationally efficient optimisation-based method for parameter identification of kinematically determinate and over-determinate biomechanical systems. Comput. Methods Biomech. Biomed. Eng. 13, 171–183. Anderson, F.C., Pandy, M.G., 2001. Static and dynamic optimization solutions for gait are practically equivalent. J. Biomech. 34, 153–161. Andriacchi, T.P., Dyrby, C.O., 2005. Interactions between kinematics and loading during walking for the normal and ACL deficient knee. J. Biomech. 38, 293–298. Audu, M.L., Kirsch, R.F., Triolo, R.J., 2003. A computational technique for determining the ground reaction forces in human bipedal stance. J. Appl. Biomech. 19, 361–371. Audu, M.L., Kirsch, R.F., Triolo, R.J., 2007. Experimental verification of a computational technique for determining ground reactions in human bipedal stance. J. Biomech. 40, 1115–1124. Choi, A., Lee, J.M., Mun, J.H., 2013. Ground reaction forces predicted by using artificial neural network during asymmetric movements. Int. J. Precis. Eng. Manuf. 14, 475–483. Collins, S.H., Adamczyk, P.G., Ferris, D.P., Kuo, A.D., 2009. A simple method for calibrating force plates and force treadmills using an instrumented pole. Gait Posture 29, 59–64. Damsgaard, M., Rasmussen, J., Christensen, S.T., Surma, E., de Zee, M., 2006. Analysis of musculoskeletal systems in the anybody modeling system. Simul. Model. Pract. Theory 14, 1100–1111. Delp, S.L., Anderson, F.C., Arnold, A.S., Loan, P., Habib, A., John, C.T., Guendelman, E., Thelen, D.G., 2007. OpenSim: open-source software to create and analyze dynamic Simulations of movement. IEEE Trans. Biomed. Eng. 54, 1940–1950. Fluit, R., Pellikaan, P., Carbone, V., Krogt, M.M.v.d., Damsgaard, M., Koopman, H.F.J.M., verdonschot, N., 2013. A new geometrically consistent musculoskeletal model of the lower extremity based on imaging and cadaver measurements. In: 11th International Symposium of Computer Methods in Biomechanics and Biomedical Engineering, Salt Lake City, Utah, USA. Greenwood, D., 1988. Dynamics of a Rigid Body, Principle of Dynamics. Prentice Hall, Englewood Cliffs NJ, pp. 389–468. Hesse, S., Konrad, M., Uhlenbrock, D., 1999. Treadmill walking with partial body weight support versus floor walking in hemiparetic subjects. Arch. Phys. Med. Rehabil. 80, 421–427. John, C.T., Anderson, F.C., Guendelman, E., Higginson, J.S., Delp, S.L., 2007. LongDuration Muscle-Actuated Simulations of Walking at Multiple Speeds. In: American Society of Biomechanics. Stanford, California, USA. Klein Horsman, M.D., Koopman, H.F., van der Helm, F.C., Prose, L.P., Veeger, H.E., 2007. Morphological muscle and joint parameters for musculoskeletal modelling of the lower extremity. Clin. Biomech. (Bristol, Avon) 22, 239–247. Kuo, A.D., 1998. A least-squares estimation approach to improving the precision of inverse dynamics computations. J. Biomech. Eng.Trans. ASME 120, 148–159. Leardini, A., Chiari, L., Della Croce, U., Cappozzo, A., 2005. Human movement analysis using stereophotogrammetry. Part 3 soft tissue artifact assessment and compensation. Gait Posture 21, 212–225. Luinge, H.J., Veltink, P.H., 2005. Measuring orientation of human body segments using miniature gyroscopes and accelerometers. Med. Biol. Eng. Comput. 43, 273–282. O'Connor, C.M., Thorpe, S.K., O'Malley, M.J., Vaughan, C.L., 2007. Automatic detection of gait events using kinematic data. Gait Posture 25, 469–474. Oh, S.E., Choi, A., Mun, J.H., 2013. Prediction of ground reaction forces during gait based on kinematics and a neural network model. J. Biomech. 46, 2372–2380. Oosterwaal, M., Telfer, S., Torholm, S., Carbes, S., van Rhijn, L.W., Macduff, R., Meijer, K., Woodburn, J., 2011. Generation of subject-specific, dynamic, multisegment ankle and foot models to improve orthotic design: a feasibility study. BMC Musculoskelet. Disord., 12. OpenSim User's Guide, Section Residual Reduction Algorithm. Last edit from 10th of June, 2013. 〈http://simtk-confluence.stanford.edu:8080/display/Open Sim/Residual+Reduction+Algorithm〉. Pamies-Vila, R., Font-Llagunes, J.M., Cuadrado, J., Alonso, F.J., 2012. Analysis of different uncertainties in the inverse dynamic analysis of human gait. Mech. Mach. Theory 58, 153–164. Pearsall, D.J., Costigan, P.A., 1999. The effect of segment parameter error on gait analysis results. Gait Posture 9, 173–183. Rao, G., Amarantini, D., Berton, E., Favier, D., 2006. Influence of body segments' parameters estimation models on inverse dynamics solutions during gait. J. Biomech. 39, 1531–1536. Rasmussen, J., Damsgaard, M., Christensen, S.T., 2000. Inverse-inverse dynamics simulation of musculoskeletal systems. In: European Society of Biomechanics Royal Academy of Medicine in Ireland.
R. Fluit et al. / Journal of Biomechanics 47 (2014) 2321–2329
Ren, L., Jones, R.K., Howard, D., 2008. Whole body inverse dynamics over a complete gait cycle based only on measured kinematics. J. Biomech. 41, 2750–2759. Richards, J.G., 1999. The measurement of human motion: a comparison of commercially available systems. Hum. Mov. Sci. 18, 589–602. Riemer, R., Hsiao-Wecksler, E.T., Zhang, X.D., 2008. Uncertainties in inverse dynamics solutions: a comprehensive analysis and an application to gait. Gait Posture 27, 578–588. Robert, T., Causse, J., Monnier, G., 2013. Estimation of external contact loads using an inverse dynamics and optimization approach: general method and application to sit-to-stand maneuvers. J. Biomech. 46, 2220–2227. Schache, A.G., Baker, R., Vaughan, C.L., 2007. Differences in lower limb transverse plane joint moments during gait when expressed in two alternative reference frames. J. Biomech. 40, 9–19. Scheys, L., Jonkers, I., Schutyser, F., Pans, S., Spaepen, A., Suetens, P., 2005. Image based methods to generate subject-specific musculoskeletal models for gait analysis. Int. Congress Ser. 1281, 62–67. Schmiedmayer, H.B., Kastner, J., 1999. Parameters influencing the accuracy of the point of force application determined with piezoelectric force plates. J. Biomech. 32, 1237–1242. Schwartz, M.H., Rozumalski, A., 2005. A new method for estimating joint parameters from motion data. J. Biomech. 38, 107–116.
2329
Schwer, L.E., 2007. Validation metrics for response histories: perspectives and case studies. Eng. Comput. (Germany) 23, 295–309. Sprague, M.A., Geers, T.L., 2003. Spectral elements and field separation for an acoustic fluid subject to cavitation. J. Comput. Phys. 184, 149–162. Taylor, R., 1990. Interpretation of the correlation coefficient: a basic review. J. Diagn. Med. Sonogr. 6, 35–39. Thelen, D.G., Anderson, F.C., 2006. Using computed muscle control to generate forward dynamic simulations of human walking from experimental data. J. Biomech. 39, 1107–1115. Vicons, 2002. Plug-in-Gait modelling instructions. Vicons Manual, Vicons612 Motion Systems. Oxford Metrics Ltd, Oxford, UK. Vanheule, V., Andersen, M.S., Wirix-Speetjens, R., Jonkers, I., Victor, J., Van den Sloten, J., 2013. Modeling of patient-specific knee kinematics and ligament behavior using force-dependent kinematics. In: XXIV Congress of the International Society of Biomechanics. Natal, Brazil. Wilken, J.M., Rodriguez, K.M., Brawner, M., Darter, B.J., 2012. Reliability and minimal detectible change values for gait kinematics and kinetics in healthy adults. Gait Posture 35, 301–307. Winter, D.A., 2009. Biomechanics and Motor Control of Human Movement, Fourth Edition John Wiley & Sons, Inc., Hoboken, NJ, USA.