Author's Accepted Manuscript

Prediction of ground reaction forces and moments during various activities of daily living R. Fluit, M.S. Andersen, S. Kolk, N. Verdonschot, H.F.J.M. Koopman

www.elsevier.com/locate/jbiomech

PII: DOI: Reference:

S0021-9290(14)00251-6 http://dx.doi.org/10.1016/j.jbiomech.2014.04.030 BM6631

To appear in:

Journal of Biomechanics

Accepted date: 18 April 2014 Cite this article as: R. Fluit, M.S. Andersen, S. Kolk, N. Verdonschot, H.F.J.M. Koopman, Prediction of ground reaction forces and moments during various activities of daily living, Journal of Biomechanics, http://dx.doi.org/10.1016/j. jbiomech.2014.04.030 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Original Article PREDICTION OF GROUND REACTION FORCES AND MOMENTS DURING VARIOUS ACTIVITIES OF DAILY LIVING R. Fluit1*, M.S. Andersen2, S.Kolk3, N. Verdonschot1,4, H.F.J.M. Koopman1 1

Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente, Enschede, The Netherlands 2

3

Department of Mechanical and Manufacturing Engineering, Aalborg University, Aalborg, Denmark

Radboud University Medical Centre, Radboud Institute for Health Sciences, Department of

Rehabilitation, Nijmegen, The Netherlands 4

Radboud University Medical Centre, Radboud Institute for Health Sciences, Orthopaedic Research Laboratory, Nijmegen, The Netherlands

*

Corresponding author:

René Fluit, Laboratory of Biomechanical Engineering, Faculty of Engineering Technology, University of Twente P.O. Box 217, 7500 AE Enschede, The Netherlands Visiting address: Campus University, Horstring W213 Tel: +31 53 489 4362 Fax: +31 53 489 2287 Email: [email protected]

Keywords: Ground reaction forces and moments, Musculoskeletal model, Inverse dynamics, Dynamic consistency, Activities of Daily Living Word count: Abstract: 248+2 Manuscript: 605+25 (introduction), 1470+106 (methods), 613+10 (results), 809+150 (discussion): 3497+291= 3788 (total) Number of Figures and Tables: 4 Figures, 6 Tables

Abstract 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 ( P ! 0.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 (

P  0.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.

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, 2009). Consequently, when solving the inverse dynamics problem, the equations of motion are violated, 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 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, 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 gaitrelated 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 selfselected 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.

2.

Methods

2.1. Subjects Nine healthy subjects (4 males and 5 females; age: 41.6 ± 15.9 yr; height: 1.74 ± 0.12 m; weight: 73.0 ± 11.1 kg, Body Mass Index: 23.9 ± 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. 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 Vicon® Plug-in-Gait model (Vicon®, 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 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.5. Prediction of GRF&Ms 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 mass-acceleration products of all individual body segments (Greenwood, 1988). Hence, 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 (Figure 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

q ˜ Fmax , with q the activation level and Fmax 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

­° F ® max °¯ 0

c p ,i

where

p , a characteristic strength function c p ,i was defined 

if pz  zcrit and p  vcrit otherwise

Fmax 0.4 BW , zcrit

(1)

0.03 m (level of the floor was set to z 0 ), vcrit

0.8

m and i the s

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

c p ,i at the discrete time step i t , the characteristic strength function was smoothed

(Figure 1, right side) according to

c p ,i

 t 3}t 8

§ 1 §1 ·· Fmax ¨ 1 r cos ¨ S n ¸ ¸ with n 1}12 2 © 12 ¹ ¹ ©

(2)

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 (Figure 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 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

5% of its maximum. For the DS trials, 0.5 s before until 0.5 s after the knee of the dominant leg was flexed

the time for which the measured vertical GRF was larger than this timespan was defined from

at 15 q 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  d 0.35 , 0.35   d 0.67 , 0.67   d 0.9 ,

0.9   to be weak, moderate, strong or excellent correlations respectively (following Taylor, 1990). 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 ( Fmax ,

zcrit

and

vcrit ) and the muscle recruitment criterion on the prediction accuracy was made. For Fmax , zcrit

and

vcrit , simulations were performed using values at 70 and 130 percent 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.

3. Results 0.621  0.980 , median 0.957 ) and the antero-posterior GRF ( ranging from 0.202  0.969 , median 0.957 ) for almost all

The model showed excellent predictions for the vertical GRF ( ranging from

activities (Table 1, Fig. 3). The magnitude of the vertical GRF was slightly but consistently

1.2%  4.0 %, median 2.5% ), whereas the magnitude of the antero-posterior GRF was consistently overestimated (M ranging from 1.1%  71.9% median 11.0% ). Nevertheless, no significant difference ( P ! 0.05 ) was found between the medians of the absolute underestimated (M ranging from

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  0.898 , median

0.810 ). In general, the most extreme activities (DS, SA and SD) showed the weakest correlations, for 0.202 ) and medio-lateral (  0.261 ) GRFs for the DS trials.

instance the antero-posterior ( 

The model showed strong correlations for the sagittal GRM ( ranging from

0.506  0.922 , median

0.789 ) and frontal GRM ( ranging from 0.199  0.801 , median 0.668 ) for most of the activities 0.155  0.782 , median 0.588 ). A significant difference ( P  0.05 ) was found between the medians of the absolute (Table 2). The transverse GRM showed the largest differences ( ranging from

mean measured transverse GRM and those of the predicted transverse GRM for all activities. 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 (

0.521  0.972 , median 0.776 ), hip abduction ( ranging from 0.257  0.946 , median 0.880 ), hip external rotation ( ranging from 0.809  0.943 , median 0.871 ), knee flexion ( ranging from 0.548  0.911 , median 0.737 ) and ankle dorsiflexion ( ranging from 0.647  0.927 , median 0.820 ). A weak correlation was found for the hip abduction moment for the DS trials only (  0.258 ,

ranging from

see Table A1 and Figure 4). For a subset of activities, Figure 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  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

33% , probably because this moment was already small using traditional inverse dynamics ( 4.3 Nm , averaging over all activities and subjects). reduction of

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

Fmax , 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 r 13 ms , and TO was predicted consistently late with an error of 16 r 7 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.2q for this DOF during normal walking, in which the internal-external 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 will result in smaller magnitudes of the transverse GRM (see Figure 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 antero-posterior 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 (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

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

the medio-lateral GRF (rRMSD of

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 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 r 13 ms and 16 r 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 r 15 ms and 9 r 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 r 10 ms and 13 r 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, 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, excitation-activation 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, 2005; Alexander and Andriacchi, 2001), image-based 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, 1999) or in the field of movement predictions using musculoskeletal optimization (Rasmussen et al., 2000).

5. Conflict of interest None of the authors have any financial or personal conflict of interest with regard to this study.

6. Acknowledgements The authors gratefully acknowledge the financial support provided by the European Commission FP7 Programme for the TLEMsafe project (www.tlemsafe.eu).

7. References Alexander, E.J., Andriacchi, T.P., 2001. Correcting for deformation in skin-based marker systems. Journal of Biomechanics 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 overdeterminate biomechanical systems. Computer Methods in Biomechanics and Biomedical Engineering 13, 171-183. Anderson, F.C., Pandy, M.G., 2001. Static and dynamic optimization solutions for gait are practically equivalent. Journal of Biomechanics 34, 153-161. Andriacchi, T.P., Dyrby, C.O., 2005. Interactions between kinematics and loading during walking for the normal and ACL deficient knee. Journal of Biomechics 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. Journal of Applied Biomechanics 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. International Journal of Precision Engineering and Manufacturing 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 and 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. Simulation Modelling Practice and 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 Transactions on Biomedical Engineering 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, 11th International Symposium of Computer Methods in Biomechanics and Biomedical Engineering, Salt Lake City, Utah, USA. Greenwood, D., 1988. Dynamics of a Rigid Body. In: (Eds.), 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. Archives of Physical Medicine and Rehabilitation 80, 421427. John, C.T., Anderson, F.C., Guendelman, E., Higginson, J.S., Delp, S.L., 2007. Long-Duration 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. Clinical Biomechanics (Bristol, Avon) 22, 239-247. Kuo, A.D., 1998. A least-squares estimation approach to improving the precision of inverse dynamics computations. Journal of Biomechanical Engineering-Transactions of the ASME 120, 148159. Leardini, A., Chiari, L., Della Croce, U., Cappozzo, A., 2005. Human movement analysis using stereophotogrammetry. Part 3. Soft tissue artifact assessment and compensation. Gait and Posture 21, 212-225. Luinge, H.J., Veltink, P.H., 2005. Measuring orientation of human body segments using miniature gyroscopes and accelerometers. Medical & Biological Engineering & Computing 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 and 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. Journal of Biomechanics 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 Musculoskeletal Disorders 12. OpenSim

User's

Guide,

section

Residual

Reduction

Algorithm.

http://simtk-

confluence.stanford.edu:8080/display/OpenSim/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. Mechanism and Machine Theory 58, 153164. Pearsall, D.J., Costigan, P.A., 1999. The effect of segment parameter error on gait analysis results. Gait and 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. Journal of Biomechanics 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. Ren, L., Jones, R.K., Howard, D., 2008. Whole body inverse dynamics over a complete gait cycle based only on measured kinematics. Journal of Biomechanics 41, 2750-2759. Richards, J.G., 1999. The measurement of human motion: A comparison of commercially available systems. Human Movement Science 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 and 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. Journal of Biomechanics 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. International Congress Series 1281, 62-67. Schmiedmayer, H.B., Kastner, J., 1999. Parameters influencing the accuracy of the point of force application determined with piezoelectric force plates. Journal of Biomechanics 32, 1237-1242. Schwartz, M.H., Rozumalski, A., 2005. A new method for estimating joint parameters from motion data. Journal of Biomechanics 38, 107-116. Schwer, L.E., 2007. Validation metrics for response histories: perspectives and case studies. Engineering with Computers (Germany) 23, 295-309. Sprague, M.A., Geers, T.L., 2003. Spectral elements and field separation for an acoustic fluid subject to cavitation. Journal of Computational Physics 184, 149-162. Taylor, R., 1990. Interpretation of the correlation coefficient: a basic review. Journal of Diagnostic Medical Sonography 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. Journal of Biomechanics 39, 1107-1115. Vicon®, 2002. Plug-in-Gait modelling instructions. Vicon® Manual, Vicon®612 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 and Posture 35, 301-307. Winter, D.A., 2009. Biomechanics and Motor Control of Human Movement, Fourth Edition ed. John Wiley & Sons, Inc., Hoboken, NJ, USA.

Figure 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

c p ,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.

Figure 2:Visualization of the skeletal model and the 28 degrees of freedom (DOFs).

Figure 3:Calculated GRFs, normalized to body weight (BW), and GRMs, normalized to BW and body height (BH) (± 1 SD around mean (shaded area)) compared with force plate data (mean (thick line) ± 1 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 q flexion angle.

Figure. 4:Computed joint moments, normalized to body weight (BW) and body height (BH), using the predicted GRF&Ms (± 1 SD about mean (shaded area)) compared to those obtained using the measured GRF&Ms (mean (thick line) ± 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 q flexion angle. */†

indicates a significant difference between the medians of the absolute maximum (*)/ absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (=0.05).

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. */† indicates a significant difference between the medians of the absolute maximum (*)/ absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (=0.05).

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. */† indicates a significant difference between the medians of the absolute maximum (*)/ absolute mean (†) of the measured and predicted value according to a two-tailed Wilcoxon Signed Rank Test (=0.05).

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).

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.

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.

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.

Activity

Vertical GRF

Antero-posterior GRF

P (%)

M (%)

 (-)

CWS +30%

3.0 (0.4)

-3.4 (1.5)

CWS -30%

2.4 (0.6)

-2.9 (0.8)

CWS

2.7 (0.4)

GI (NDL)

2.4 (0.7)

GI (DL) GT (NDL)

Medio-lateral GRF

P (%)

M (%)

 (-)

P (%)

M (%)

 (-)

0.954

9.0 (1.3)

4.1 (8.3)

0.948

11.8 (2.6)

1.0 (10.8)

0.753

0.967

* 8.4 (1.9)

26.4 (9.1)

0.965

7.6 (1.9)

-1.9 (4.0)

0.866

-3.0 (0.9)

0.957

* 9.3 (1.2)

15.8 (9.3)

0.957

* 9.8 (2.3)

1.9 (6.9)

0.818

-3.4 (2.2)

0.976

6.0 (1.1)

1.1 (11.5)

0.968

16.7 (12.4)

-8.6 (9.6)

0.898

2.3 (0.2)

-2.8 (2.5)

0.980

* 8.1 (2.2)

13.4 (7.4)

0.944

11.9 (2.2)

-5.1 (5.5)

0.842

2.5 (0.6)

-2.1 (1.7)

0.976

6.3 (1.5)

3.3 (5.8)

0.962

17.8 (6.9)

-4.7 (17.7)

0.863

GT (DL)

3.2 (0.9)

-2.0 (1.1)

0.947

9.5 (1.2)

11.0 (6.6)

0.902

15.2 (4.9)

-6.1 (7.9)

0.726

OBS10

2.3 (0.3)

-2.2 (0.5)

0.962

7.5 (1.2)

6.4 (6.1)

0.966

9.8 (3.0)

8.6 (8.7)

0.836

OBS20

2.5 (0.2)

-2.1 (1.2)

0.957

8.0 (1.3)

6.1 (11.3)

0.959

11.1 (2.8)

12.0 (12.4)

0.810

OBS30

2.7 (0.5)

-1.4 (1.0)

0.943

7.9 (1.3)

5.1 (7.3)

0.969

* 12.4 (3.5)

14.1 (15.4)

0.713

DS

2.7 (0.6)

-2.5 (1.4)

0.621

* 35.0 (15.1)

71.9 (139.6)

0.202

34.7 (14.0)

11.2 (43.8)

0.261

SA

3.1 (1.3)

-4.0 (1.3)

0.941

* 26.0 (7.6)

57.5 (65.3)

0.436

23.2 (4.5)

-9.9 (14.7)

0.609

SD

4.7 (0.6)

-1.2 (0.9)

0.895

22.3 (2.8)

15.1 (12.6)

0.637

26.2 (2.1)

-10.3 (12.4)

0.455



Activity

Sagittal GRM

Frontal GRM

Transverse GRM

P (%)

M (%)

 (-)

P (%)

M (%)

 (-)

CWS +30%

6.8 (2.3)

6.6 (28.2)

0.929

16.4 (7.3)

-10.5 (15.6)

0.540

CWS -30%

† 9.6 (2.3)

-15.2 (6.1)

0.872

10.7 (3.4)

-7.3 (11.2)

0.777

P (%)

M (%)

 (-)

*† 33.1 (7.5)

-56.7 (20.6)

0.762

*† 46.0 (9.5)

-46.5 (19.4)

0.598

CWS

7.4 (1.8)

-8.8 (6.3)

0.922

13.4 (5.8)

-10.8 (10.0)

0.684

*† 36.8 (6.3)

-48.3 (14.6)

0.704

GI (NDL)

9.3 (2.9)

0.4 (7.1)

0.727

8.4 (2.3)

-9.2 (15.3)

0.703

† 26.8 (10.9)

-33.3 (27.6)

0.588

11.4 (4.7)

-4.8 (23.2)

0.878

10.6 (3.9)

-7.2 (9.1)

0.801

*† 30.9 (7.6)

-47.6 (18.2)

0.525

14.6 (12.1)

0.680

* 13.0 (4.8)

0.7 (18.8)

0.488

*† 68.7 (10.4)

-79.4 (3.4)

-0.128

GI (DL) GT (NDL)

* 8.5 (2.5)

GT (DL)

10.4 (3.0)

3.7 (10.9)

0.777

14.0 (4.8)

0.6 (27.9)

0.539

*† 56.1 (8.6)

-52.0 (41.1)

0.157

OBS10

7.1 (1.8)

-6.3 (7.2)

0.915

12.7 (4.6)

3.7 (20.5)

0.694

*† 33.9 (7.2)

-46.2 (13.0)

0.782

OBS20

7.7 (2.6)

-6.4 (6.0)

0.863

14.6 (4.3)

8.3 (21.7)

0.651

*† 33.2 (8.4)

-52.0 (11.7)

0.759

OBS30

8.7 (1.7)

-4.4 (7.4)

0.801

15.3 (6.7)

15.3 (28.6)

0.641

*† 33.2 (4.4)

-50.6 (13.7)

0.657

DS

15.6 (9.8)

50.5 (59.0)

0.506

24.9 (6.7)

37.1 (23.2)

0.199

*† 42.5 (7.2)

91.6 (147.1)

0.230

SA

12.4 (3.3)

-14.0 (14.7)

0.647

9.7 (2.0)

17.4 (13.7)

0.781

*† 53.9 (8.6)

-33.8 (22.5)

-0.155

SD

8.3 (3.1)

5.5 (9.1)

0.830

9.7 (2.9)

15.7 (11.7)

0.781

*† 67.0 (5.6)

-61.0 (15.4)

-0.155



Activity



Vertical GRF

Antero-posterior

Medio-lateral

Sagittal GRF

Frontal GRM

Transverse GRM

(N/kg)

GRF (N/kg)

GRF (N/kg)

(Nm/kg)

(Nm/kg)

(Nm/kg)

CWS +30%

0.85 (0.17)

0.43 (0.06)

0.23 (0.07)

0.17 (0.07)

0.13 (0.03)

0.28 (0.08)

CWS -30%

0.64 (0.15)

0.30 (0.05)

0.12 (0.03)

0.22 (0.07)

0.08 (0.02)

0.21 (0.07)

CWS

0.74 (0.13)

0.38 (0.07)

0.17 (0.04)

0.18 (0.05)

0.11 (0.02)

0.22 (0.06)

GI (NDL)

0.54 (0.18)

0.16 (0.04)

0.19 (0.07)

0.11 (0.06)

0.05 (0.03)

0.08 (0.06)

GI (DL)

0.63 (0.10)

0.29 (0.07)

0.17 (0.03)

0.19 (0.05)

0.07 (0.02)

0.17 (0.06)

GT (NDL)

0.52 (0.12)

0.18 (0.06)

0.23 (0.07)

0.16 (0.06)

0.06 (0.02)

0.16 (0.06)

GT (DL)

0.91 (0.60)

0.32 (0.08)

0.21 (0.06)

0.20 (0.06)

0.09 (0.03)

0.24 (0.10)

OBS10

0.64 (0.17)

0.32 (0.06)

0.18 (0.05)

0.23 (0.11)

0.12 (0.09)

0.23 (0.06)

OBS20

0.68 (0.12)

0.34 (0.07)

0.20 (0.07)

0.22 (0.08)

0.11 (0.03)

0.22 (0.06)

OBS30

0.74 (0.24)

0.33 (0.07)

0.22 (0.07)

0.23 (0.08)

0.12 (0.04)

0.21 (0.07)

DS

0.81 (0.36)

0.30 (0.22)

0.24 (0.16)

0.22 (0.05)

0.17 (0.07)

0.10 (0.03)

SA

0.95 (0.25)

0.21 (0.06)

0.22 (0.04)

0.21 (0.05)

0.14 (0.05)

0.19 (0.06)

SD

0.34 (0.12)

0.22 (0.07)

0.42 (0.17)

0.09 (0.02)

0.11 (0.05)

0.13 (0.04)

P=1



P=2

P=4

Min/Max Aux

M (%)

rRMSD (%)

M (%)

rRMSD (%)

M (%)

rRMSD (%)

M (%)

rRMSD (%)

Vertical GRF

-0.3

-0.2

-0.1

-0.1

0.0

0.0

-0.2

-0.1

Antero-posterior GRF

5.0

0.8

2.0

0.3

-1.7

-0.1

1.4

0.3

Medio-lateral GRF

3.7

-0.7

-0.9

-0.6

0.7

0.2

0.0

-0.1

Sagittal GRM

-2.9

-1.0

0.1

-0.7

1.2

0.3

0.4

-1.0

Frontal GRM

-0.7

-2.1

0.3

-0.9

-0.6

0.1

0.9

-0.2

Transverse GRM

2.6

0.0

2.9

0.9

-1.4

-0.1

1.1

0.6





Fmaz -30%

+30%



vcrit -30%

+30%

zcrit -30%

+30%

M

rRMSE

M

rRMSE

M

rRMSE

M

rRMSE

M

rRMSE

M

rRMSE

(%)

(%)

(%)

(%)

(%)

(%)

(%)

(%)

(%)

(%)

(%)

(%)

Vertical GRF

-0.5

-0.2

0.3

0.2

0.4

-0.4

-0.4

0.5

-0.2

-0.1

0.2

0.1

Antero-posterior GRF

-1.7

-0.3

1.6

0.3

-5.5

-0.9

3.2

0.4

-1.9

-0.4

1.3

0.3

Medio-lateral GRF

-0.1

-0.2

0.1

0.1

-0.5

-0.2

0.6

0.2

-0.0

-0.1

0.0

0.1

Sagittal GRM

-1.1

-0.6

1.0

0.6

-4.1

-1.7

3.7

1.5

-0.8

-0.9

0.6

0.7

Frontal GRM

2.3

-0.0

-1.7

0.1

1.5

-0.3

-1.0

0.2

0.4

-0.2

-0.4

0.2

Transverse GRM

1.2

-0.1

-1.0

0.2

1.5

-1.2

-0.6

0.6

0.9

-0.1

-0.7

-0.1

Method

Smooth Transition Assumption

Artificial Neural Network

Optimization Approach

(Results from Ren et al., 2008)

(Results from Oh et al., 2013)

(Method proposed in this study)

N=3

N=5

Participants CWS

CWS +30%

CWS +30%

RMSD (N/kg rRMSD (%)

or Nm/kg)

or Nm/kg)

or Nm/kg)

6.6 ( 1.6)

0.649 (0.182) 5.8 (1.0)

0.74 (0.13)

6.6 (1.1)

0.85 (0.17)

6.9 (1.3)

8.9 (1.3)

0.154 (0.057) 7.3 (0.8)

0.38 (0.07)

9.3 (2.0)

0.43 (0.06)

8.5 (1.6)

0.335 (0.012)

24.1 (1.6)

0.040 (0.022) 10.9 (1.8)

0.17 (0.04)

14.9 (3.4)

0.23 (0.07)

16.6 (4.6)

0.225 (0.029)

14.2 (1.4)

0.081 (0.045) 9.9 (1.9)

0.18 (0.05)

12.4 (3.5)

0.17 (0.07)

10.4 (3.7)

0.148 (0.013) 32.5 (4.3)

0.235 (0.018)

37.7 (1.0)

0.052 (0.029) 22.8 (4.9)

0.11 (0.02)

22.9 (5.9)

0.13 (0.03

27.1 (9.0)

0.039 (0.015) 26.2 (9.4)

0.075 (0.009)

41.7 (1.8)

0.032 (0.018) 25.5 (4.5)

0.22 (0.06)

40.6 (11.3)

0.28 (0.08)

38.4 (10.9)

RMSD (N/kg

rRMSD

RMSD (N/kg

or Nm/kg)

or Nm/kg)

(%)

Vertical GRF

0.710 (0.190) 5.6 (1.5)

1.012 (0.223)

AP GRF

0.473 (0.068) 10.9 (0.83)

0.507 (0.085)

ML GRF

0.191 (0.034) 20.0 (2.7)

Sagittal GRM

0.199 (0.106) 12.2 (4.8)

Frontal GRM Transverse GRM

rRMSD (%)

CWS RMSD (N/kg

RMSD (N/kg rRMSD (%)



N=9

CWS

rRMSD (%)

Figure 1

Figure 2

Figure 3

Figure 4