UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 SURROGATE BASED MODELING OF JOINT CONTACT MECHANICS By YICHUNG LIN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 PAGE 2 2 2008 YiChung Lin PAGE 3 3 To my parents PAGE 4 4 ACKNOWLEDGMENTS First, I would like to sin cerely thank Dr. Fregly and Dr. Haftka for the direction they have given me. Their valuable suggestions and dedica tions contributed greatly. I would also like to thank Dr. Banks, Dr. Kim, and Dr. Jiang for serv ing on my doctoral committee. It was a great honor to have them. I would also like to thank all the Computational Biomechanics Lab members for their support in my past and present work. Last but not least, I thank my parents and my wife HsiuJen for being so supportive. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........7 LIST OF FIGURES.........................................................................................................................8 ABSTRACT...................................................................................................................................10 CHAP TER 1 INTRODUCTION..................................................................................................................12 Need for Efficient Contact Model Evaluation........................................................................ 12 Specific Aims..........................................................................................................................14 2 TWODIMENSIONAL SURROGATE CONTACT MODELING FOR COMPUT ATIONALLYEFFICIENT DYNAM IC SIMULATION OF TOTAL KNEE REPLACEMENTS.................................................................................................................15 Introduction................................................................................................................... ..........15 Methods..................................................................................................................................17 Surrogate Contact Model Development..........................................................................17 Surrogate Contact Model Evaluation.............................................................................. 21 Surrogate Contact Model Application............................................................................. 23 Results.....................................................................................................................................25 Discussion...............................................................................................................................26 3 THREEDIMENSIONAL SURROGATE CONTACT MODELS FOR COMPUT ATIONALLYEFFICIENT MULTI BODY DYNAMIC SIMULATIONS.......... 37 Introduction................................................................................................................... ..........37 Methods..................................................................................................................................39 Surrogate Contact Model Development..........................................................................39 Surrogate Contact Model Evaluation.............................................................................. 45 Results.....................................................................................................................................50 Discussion...............................................................................................................................51 4 SIMULTANEOUS PREDICTION OF MUSC L E AND CONTACT FORCES IN THE KNEE DURING GAIT........................................................................................................... 63 Introduction................................................................................................................... ..........63 Methods..................................................................................................................................65 Data Collection................................................................................................................65 Inverse Dynamics Loads.................................................................................................66 PAGE 6 6 Musculoskeletal Model................................................................................................... 66 Surrogate Contact Model................................................................................................. 68 Muscle and Contact Force Prediction.............................................................................. 70 Results.....................................................................................................................................72 Discussion...............................................................................................................................73 5 CONCLUSION..................................................................................................................... ..83 APPENDIX COMPUTATIONAL WEAR PREDICTION.................................................... 85 LIST OF REFERENCES...............................................................................................................86 BIOGRAPHICAL SKETCH.........................................................................................................96 PAGE 7 7 LIST OF TABLES T able page 21 Comparison between the wear volume predictions (mm3) made with the elastic foundation contact model and the surrogate contact model for each of the nine nominal component placements......................................................................................... 33 31 Comparison between the wear volume predictions (mm3) made with the elastic foundation contact model and the surrogate contact model for dynamic simulations with different input curves................................................................................................. 57 41 Eight different optimization problem fo r mulations. C1 to C4 represent four cheating formulations while H1 to H4 represent four honest formulations............... 78 42 Summary of root mean square and ma xim um absolute joint kinematics errors produced by surrogate contact model for problem formulation C1................................... 79 43 Summary of root mean square and ma xim um absolute joint dynamics errors produced by surrogate contact model for problem formulation C1................................... 80 44 Summary of root mean square and maxi m um absolute muscle force errors produced by surrogate contact model fo r problem formulation C1................................................... 81 PAGE 8 8 LIST OF FIGURES Figure page 21 Overview of surrogate contact mode l creation and im plementation process.................... 32 22 Overview of global and local sensitiv ity analyses of predicted wear volume perform ed with the surr ogate contact model..................................................................... 33 23 Rootmeansquare errors in joint moti ons and co ntact forces for each nominal component placement........................................................................................................ 34 24 Predicted wear volume as a function of la rge variations in component placem ents as calculated from 400 dynamic simulations pe rformed with the surrogate contact model..................................................................................................................................34 25 Box plot distribution of predicted w ear volume generated by four Monte Carlo analyses using the surrogate contact model....................................................................... 35 26 Predicted wear volume as a function of variations in loads and fe moral flexion angle.... 36 31 A 6 degree of freedom total knee replacement contact model........................................... 56 32 Overview of the surrogate modeling approach.................................................................. 58 33 Box plot distribution of predicted w ear volume generated by four Monte Carlo analyses using the surrogate contact model....................................................................... 59 34 Comparison between the six pose parame ters m ade with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation....... 60 35 Comparison between the mediallateral contact fo rces made with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation..................................................................................................................... ......61 36 Comparison between the mediallateral contact torques m ade with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation..................................................................................................................... ......62 41 Overview of the process of using muscul oske letal model to predict muscle forces and contact loads simultaneously over the entire gait cycle.............................................. 78 42 A 12 degree of freedom tota l knee replacem ent contact model......................................... 79 43 Simultaneous approache for predictin g m uscle forces and contact loads.......................... 80 44 Comparison between the in vivo contact m easurements (color in grey) and predictions from eight optimi zation problem formulations............................................... 81 PAGE 9 9 45 Predictions of muscle forces from eight optim ization problem formulations................... 82 PAGE 10 10 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SURROGATE BASED MODELING OF JOINT CONTACT MECHANICS By YiChung Lin December 2008 Chair: Benjamin J. Fregly Cochair: Raphael T. Haftka Major: Mechanical Engineering Musculoskeletal computer models are useful for estimating internal quantities that cannot be measured experimentally, designing new medical devices and rehabil itation approaches, and predicting the outcome of surgical procedures. Unfortunately, use of articular contact in such models makes computational speed a limiting f actor, rendering dynamic simulations either completely intractable or else so slow that optimization is impossible. Lacking of articular contact will significantly affect many contact quantities predicted from musculoskeletal computer models (e.g., muscle forces, ligament strains, bone load s, and cartilage and implant wear). Computational limitations in other engineering disciplines have been overcome through the use of surrogate modeling. Surrogatebased modeling involves fitting a model to a model, where the original model is computationa lly expensive and the surrogate model is computationally cheap. The sample data points us ed to fit the surrogate model are generated by running the original model repeatedly with different combinations of input variables. Once fitted, the surrogate model can be used in place of the original model in simulations or optimizations to eliminate computational cost as a limiting factor. Though surroga te modeling techniques have PAGE 11 11 successfully eliminated computational bottlenecks in other fields, they have not yet been applied to articular contact problems. Primary objectives were twofold. First, we present a computational evaluation and practical application of the proposed surrogatebased modeli ng approach using dynamic wear simulation of a total knee replacement constrained to both twoand threedimensional motions in a Stanmore machine. The sample points needed for surrogate modeling fitting are generated by an elastic foundation contact model. Second, we apply the surrogate contact model into the musculoskeletal model to simultaneously calcul ate medial and lateral tibiofemoral contact forces, patellofemoral contact forces, and the muscle forces. This project provides a computationally efficient way to apply diffe rent types of contact models to predict physiologically significant contact quantities. PAGE 12 12 CHAPTER 1 INTRODUCTION Need for Efficient Cont act Model Evaluation Musculoskeletal com puter models are useful for estimating physiol ogical quantities that cannot be measured experimentally (Buchana n and Shreeve, 1996), designing new medical devices and rehabilitation approaches (Neptune et al., 2000), and pred icting the outcome of surgical procedures (Delp et al., 1996). In such models, muscles, ligaments, and bones interact (as in the real musculoskeletal system) to de termine the motions and loads experienced by the anatomic structures. Much of this interaction occu rs at the joints, where forces from each type of structure influence the forces exer ted on and by the other structures. Lack of articular (i.e., surfacesurface) contac t in musculoskeletal computer models can lead to inaccurate prediction of quantities influenced by muscle, ligament, and bone interactions. At least four types of predicti ons can be adversely affected by omission of articular contact models: 1) Muscle force predictions the amount of constraint present in traditional engineering joint models (e.g., hinge versus ball and sock et) can significantly influence muscle force predictions during optimization analysis (Challis and Kerwin, 1993; Glitsch and Baumann, 1997; Li et al., 1999; Pierce and Li, 2005) 2) Ligament strain predictions articular geometry, via its influence on joint kinematics, can have a significant eff ect on the strains calc ulated in individual ligaments (Blankevoort et al., 1991; Wilson et al ., 1998). 3) Bone remodeling predictions bone remodeling simulations are sensitive to the appl ied joint contact and muscle forces, both of which are influenced by articula r contact conditions (Bitsakos et al., 2005). 4) Cartilage and implant wear predictions contact forces, pressu res, and sliding conditi ons within a joint are dependent on articular surface geometry as well as muscle cocontractions and ligament forces (Kwak et al., 2000; Fregly et al., 2005). PAGE 13 13 Unfortunately, use of articular contact in musculoskeletal computer models makes computational speed a limiting factor, rendering dynamic simulations or optimization analysis too slow to perform. Unlike engi neering joint models, articular contact models require repeated evaluation of surface geometry, and these geometry evaluations consume the vast majority of the CPU time in a dynamic contact simulation (Bei and Fregly, 2004). A slow dynamic simulation requiring hours or days of CPU time (Giddings et al., 2001; G odest et al., 2002) means painfully slow repeated analyses and impractical optimization studies, which typi cally require hundreds or even thousands of repeated analyses. While se veral labs have developed their own articular contact models to address computational speed issues (Blankevoort et al., 1991; Pandy et al., 1997; Kwak et al., 2000; Cohen et al., 2001; Piazza and Delp, 2001; Chao, 2003; Elias et al., 2003; Bei and Fregly, 2004; Car untu and Hefzy, 2004), only three m odels have been used to perform dynamic simulations (Piazza and Del p, 2001; Caruntu and Hefz y, 2004; Bei and Fregly, 2004), only one functions within a general dynam ic simulation environment (Bei and Fregly, 2004), and none can be labeled as fast for pe rforming repeated dynamic simulations involving articular contact at one or more joints. Similar computational challenges are being tackled in other engineering fields using surrogatebased modeling (Giunta et al., 1997; Kurt aran et al., 2001; Queipo et al., 2002). This approach involves fitting a computationally cheap model to data points generated by a computationally expensive model. The surroga te model (also sometimes called a response surface) is then used in place of the original model to elim inate a computational bottleneck either in some aspect of a simulation (e.g., the contact calculations) or in an entire simulation (e.g., by providing fast evaluation of the cost function and constraints). PAGE 14 14 Unlike traditional engineering applications, join t contact analyses pos e unique challenges to existing surrogatebased modeling approach es. Sampling data points within a hypercube defined by the upper and lower bounds of the relative pose parameters will result in few physically realistic data points fo r surrogate model creation. Most sample points will be either out of contact (i.e., zero contact force) or in excessive penetration (i .e., unrealistically high contact force) since contact forces are highly se nsitive to the small displacements in the contact normal direction. Consequently, tr aditional surrogatebas ed modeling approaches fail to provide the accuracy needed to perform dynamic simulations. A special surrogatebased modeling approach adapted to the specific needs of jo int contact analyses is therefore required. Specific Aims In short, the goal of present work is to develop and evaluate a novel surrogate contact modeling ap proach to permit rapid dynamic simulation, sensitivity analysis, and optimization of musculoskeletal models with articular contact. The objectives of the current work can be summarized as follows: 1. Evaluate the effect of varia tions in component placements a nd input profiles on wear volume generated by a knee simulator machine. 2. Develop and evaluate an extended surrogate contact modeling approach capable of simulating threedimensional situations. 3. Predict muscle and contact forces simultaneously at the knee using a musculoskeletal model incorporating threedimensional surrogate contact models. 4. Evaluate the muscle and contac t forces prediction using in vivo contact force measurements provided by an instrumented knee replacement. PAGE 15 15 CHAPTER 2 TWODIMENSIONAL SURROGATE CONTACT MODELING FOR COMPUT ATIONALLYEFFICIENT DYNAM IC SIMULATION OF TOTAL KNEE REPLACEMENTS Introduction For m ore than three decades, total knee repl acement (TKR) surgery has been performed on patients with severe osteoarthritis. Though mode rn TKR surgery has a high success rate, some patients still experience substantial functional impairment postsurgery compared to healthy individuals (Noble et al., 2005). Furthermore, patients today desire more than just pain relief, seeking a high level and broad variety of daily activities following TKR surgery (e.g., tennis, hiking, gardening, even jogging) (M ancuso et al., 2001; Weiss et al., 2002). Thus, maximization of durability and minimization of functional limitations have become two primary design goals. One way to pursue these goals is through th e development of com putational technology that permits sensitivity and optimization studies of new TKR designs. A recent step in this direction has been the development of computat ional models of knee wear simulator machines (Giddings et al., 2001; Halloran et al., 2005a; Rawlinson et al., 2006; Fregly et al., 2008; Zhao et al., 2008). Such models permit dynamic contact and wear simulations of knee implant designs to be performed in a timeand costefficient ma nner, with model predictions being validated against experimental motion and wear measuremen ts (Halloran et al., 20 05a; Knight et al., 2007; Fregly et al., 2008; Zhao et al., 2008). Though fi nite element models can be used to perform contact calculations for these simulations, recent work has shown that elastic foundation contact models produce similar contact forces and pressu res in a fraction of the CPU time (Halloran et al., 2005b). Even so, CPU times on the order of 5 to 10 minutes for a onecycle dynamic contact simulation can still be prohibitive for design sensitivity and optimiza tion studies requiring thousands of simulations. PAGE 16 16 Computational limitations in other engineering disciplines have been overcome through the use of surrogate modeling approaches. Surrogate modeling (also known as metamodeling or response surface approximation) involves repl acing a computationally costly model with a computationally cheap model constructed using data points sampled from the original model. Once the surrogate model is constructed, it is used in place of the original computationally costly model when subsequent engineering analyses are performed. A vari ety of surrogate model fitting methods have been proposed, including polynomial response surfaces (Box and Draper, 1987; Myers and Montgomery, 1995; Khuri and Cornell, 1996), Kriging (Sacks et al., 1989; Jones et al., 1998), radial basis functions (Wendland, 1995; Chen et al., 1996), splines (Wahba, 1987), and support vector machines (Girosi, 1998; Va pnik, 1998). While surrogate models have been utilized successfully fo r a variety of engineering applications (Koch et al., 1999; Liu et al., 2000; Cox et al., 2001; Queipo et al., 2002; Queipo et al., 2005; Wang and Shan, 2007), only a few studies have used surrogate models to fit inputo utput relationships from a computational contact model (Bouzid et al., 1998; Chang et al., 1999; Lin et al., 2006). Furthermore, none of these studies has proposed a surrogate modeling technique to improve the computational efficiency of dynamic contact simulations. This study proposes a surrogate contact modeling approach for performing dynamic contact and wear simulations of total knee replac ements in a computationally efficient manner. The approach addresses the unique challenges involved in applying surrogate modeling techniques to joint contact problems. The primar y objectives of this chapter are two fold. The first is to evaluate the computational speed and accuracy of the proposed surrogate contact modeling approach by comparing its results with those generated by an elastic foundation contact model. Sample points fo r the surrogate contact model are generated by the same elastic PAGE 17 17 foundation contact model used in the comparative dynamic wear simulations. The second objective is to demonstrate the practical applicab ility of the proposed approach by analyzing the sensitivity of knee replacement wear predictions to realistic varia tions in input motions and loads and component placements. A Monte Carlo appr oach requiring thousands of dynamic wear simulations is used to perform the sensitivity study. Both objectives are pursued using a dynamic contact model of a total knee repl acement constrained to planar motion in a Stanmore machine, and both demonstrate the significa nt computational benefits that surrogate contact modeling can provide for analysis of TKR designs. Methods Surrogate Contact Model Development We have developed a novel surrogate m odeling approach to replace computationally costly contact calculations in dynamic simulations of total knee replacements. To develop the approach, we utilized a threedimensional (3D) elastic foundation contact model (An et al., 1990; Blankevoort et al., 1991; Li et al., 1997; Pandy et al., 1998; Fregly et al., 2003) of a cruciateretaining commercial kn ee implant (Osteonics 7000, Str yker Howmedica Osteonics, Inc, Allendale, NJ) constrained to sagittal plane moti on in a Stanmore simulator machine (Walker et al., 1997). The model possessed three degrees of freedom (DOFs) rela tive to ground: tibial anteriorposterior translation, femoral s uperiorinferior translation, and femoral flexionextension. Similar to a Stanmore mach ine, femoral flexionextension rotation was prescribed while the remaining two DOFs were loaded using ISO standard input curves (DesJardins et al., 2000 ). A pair of parallel spring bumper lo ads was also applied to the tibia in the anteriorposterior direction to approximate the effect of ligam ent forces (Walker et al., 1997). The EF contact model was treated as an applied load on the femoral component and tibial insert and used linear elastic material properties (Fregly et al., 2005 ). The dynamic equations for the PAGE 18 18 model were derived via Kanes method using Autolev symbolic manipulation software (OnLine Dynamics, Sunnyvale, CA). The equations were incorporated into a Matlab program (The Mathworks, Natick, MA) that was used to perfor m forward dynamic simulations with the stiff solver ode15s. To develop a surrogate contact model to re place the EF contact model in the larger dynamic model described above, we followed a f ourstep process: 1) design of experiments (DOE), 2) computational experiments, 3) surr ogate model selection, and 4) surrogate model implementation (Figure 21). The first step, DOE is a statistical method for determining at which locations in design space computational experiments should be performed with the computationally costly model to sample the outputs of interest (Figure 21A). For our sagittal plane knee implant model, locations in design space (or sample points) were defined by the following three relative pose parameters: tibi al anteriorposterior (AP) translation X femoral superiorinferior (SI) translation Y and femoral flexionextension rotation The associated outputs of interest calculat ed by the EF contact model were the AP contact force Fx and the SI contact force Fy. To minimize the total number of sample points while maximizing the quality of the resulting surrogate model fit, we chose a Hammersley qua sirandom (HQ) sampling method (Hammersley, 1960). This method uses an optimal design scheme to uniformly place sample points within a multidimensional hypercube (K alagnanam and Diwekar, 1997; Diwekar, 2003). The second step, computational experiments, involved performing repeated simulations with the original contact model at the sample points defined in the first step. Because contact forces are highly sensitive to some pose paramete rs but not others, we made one modification to the computational experiments performed at the HQ sample points (Figure 21B). A traditional HQ sampling method would sample values of the three inputs X Y and within their specified PAGE 19 19 bounds and use the EF contact model to calculate the two outputs Fx and Fy at each sample point. However, since SI contact force Fy is highly sensitive to small variations in SI translation Y (Lin et al., 2006; Fregly et al., 2008) this approach produced a large number of sample points that were either out of contact or deeply penetrating. Inclusion of these unrealistic points in the surrogate model fitting pr ocess reduced the accuracy of the re sulting surrogate contact model. To resolve this issue, we inverted the HQ samp ling method for superiorinferior translation by sampling forces in the sensitive direction (i.e., Fy) and displacements in the insensitive directions (i.e., X and ). Specifically, we used the HQ method to sample 300 pairs of X and values. For each of these pairs, we generated five sample poi nts by performing five static analyses with the EF contact model using Fy values ranging from 600 to 2400 N (approximately 3x body weight) in 600 N increments plus a 10 N value for dete rmining contact initia tion. The outputs of each static analysis were Fx and Y. Evaluating 1500 sample points (X Fy, ) using static analyses required approximately 3 hours of CP U time on a 3.4 GHz Pentium IV PC. The third step, surrogate mode l selection, involved determinin g the most appropriate type of surrogate model for fitting the inputoutput relationships defined by the computational experiments in the second step. Based on its ability to interpolate multidimensional nonuniformly spaced sample points, Kriging (Sacks et al., 1989; Cressie, 1992) was selected as the surrogate modeling method to fit the inputo utput relationships defined by the 1500 sample points. Originally developed for geostatistics a nd spatial statistics, Kriging has been widely applied in many different fields (Chang et al., 1999; Gupta et al., 2006). In mathematical terms, Kriging utilizes a combination of a polynomial trend model and a systematic departure model: y( x) = f( x) + z(x) (21) where x is a vector of design variable inputs, y( x ) is the output to be fitted, f( x) is a polynomial trend model, and z( x) is a systematic departure model whos e correlation structure is a function of PAGE 20 20 the distances between sample points. The f( x) term approximates the design space globally while the z(x ) term interpolates local deviations from f( x). The necessary parameters for z(x) were obtained through maximum likelihood estimatio n (Sacks et al., 1989). The Matlab toolbox DACE (Lophaven et al., 2002) was us ed to construct all necessary Kriging models using a cubic polynomial trend function and a cubic spline corre lation function based on preliminary tests. The fourth step, surrogate model implementa tion, required the development of a special model inversion process to accommodate our m odified HQ sampling approach (Figure 21C). During a dynamic simulation, the EF contact model requires X Y and as inputs and provides Fx and Fy as outputs, whereas the surrogate contact model requires X Fy, and as inputs and predicts Fx and Y as outputs. To address this inconsistency, we fitted Fy as a function of Y using the general form of a Hertzian contact model (Equation 22). 10 ,,ynX yyNFkXYXY (22) In this equation, yF is the predicted SI contact force, ky is the contact stiffness, 10 NY is the SI translation at contact initiation, and ny is the contact exponent. Use of Equation 22 to calculate yF required fitting ky, 10 NY and ny as a function of X and Each of the 300 (X ) sample pairs had five sampled values of yF (i.e., 10N, 600N, 1200N, 1800N, 2400N) and five associated output values of Y. Thus, 10NY was already known for each (X, ) sample pair. To find ky and ny for each (X, ) sample pair, we linearized Equation 22 by taking 10log of both sides and solved for 10log()yk (and hence ky) and ny using linear least square s (5 equations in 2 unknowns). This process yielded one value of ky, 10 NY and ny for each ( X ) sample pair. To calculate the AP contact force x F, we used Equation 23 based on the observation that Fx is close to a constant fraction of yF at each ( X ) sample pair. This fraction ( ratio ) is the PAGE 21 21 average of the five ( Fx, Fy) quotients obtained for each sample pair, and yF is the SI contact force predicted by Equation 22. Thus, the final surrogate contact model utilized four separate Kriging models to fit ky, 10NY ny, and ratio as a function of X and (Figure 21D), permitting the calculation of (,,)yFfXY and (,, xFfX )yF from Equation 22 and 23 at any point during a dynamic simulation (Figure 21E). ,xFratioX yF (23) Surrogate Contact Model Evaluation To evaluate the computational speed and accuracy of our surrogate contact modeling approach, we performed nine dynamic wear simu lations with both the su rrogate contact model and the EF contact model used to generate it. Each wear simulation util ized different nominal component placements that were variations away from the original placements (Figure 22A). The specific variations simulated were all possible combinations of three SI locations of the femoral flexion axis (original 10 mm) and thr ee AP slopes of the tibial insert (original 10 deg). To cover all possible dynamic simulation pa ths, we specified the ranges of the 300 HQS points using the limits of motion obtained from dynamic simulations performed using the EF contact model when the SI locations of the femoral flexion axis and AP slope of the tibial insert were set to their extreme values. We quantified differences between the sets of results by calculating root mean square erro rs (RMSE) in predicted contact forces, component translations, and wear volumes. Generation of wear volume predictions with the surrogate contact model required the development of a separate surroga te modeling approach based on pr ediction of medial and lateral center of pressure (CoP) locations. This approach used Archards wear law (Archard and Hirst, 1956, and Appendix A) to predict wear volume for any specified number of cycles from the time PAGE 22 22 history of contact forces and CoP slip ve locities on each side over a onecycle dynamic simulation (Equation 24). 22 1111nn ijijijij ijijVNkFdNkFt v (24) In this equation, V is the total wear volume of th e medial and lateral sides over N cycles, N is the assumed total number of cycles (5 million), i is a discrete time instant during the onecycle dynamic simulation (1 through n ), j is the side (medial or lateral), k is the assumed material wear rate (1 107 mm3/Nm; Fisher et al ., 1994), ijF is the total contact force at time instant i on side j and ijd is the sliding distance of the CoP at time instant i on side j. Due to symmetry of the implant geometry, the predicted value ofijF on each side was taken as half the total contact force obtained from Equation 22 and 23. Prediction of ijd on each side required determining the magnitude of the slip velocity ijv of the CoP and multiplying by the time increment t (Fregly et al., 2005). To find ijv we used the kinematic concept of two points fixed on a rigid body (Equation 25). AOABOCoP ij vv p (25) In this equation, A represents the tibial insert, B represents the femoral component, O is the femoral coordinate system origin located at the intersection of the flexionextension and superiorinferior translation axes, CoP is the medial or lateral center of pressure location treated as a point fixed in body B, AOv is the velocity of point O with respect to reference A as determined from the superiorinferior translation of the femoral component, AB is the angular velocity of body B with respect to reference frame A as determined from the flexionextension of the femoral component, and OCoPp is the position vector from point O to point CoP. While the values of AOv and AB were obtained directly from the dynamic simulations, calculation of PAGE 23 23 the CoP location on each side re quired six additional Kriging models (i.e., 3 components per side). Each Kriging model was created by fitting the average of the 5 CoP locations obtained from each (X, ) sample pair. Wear volumes calculated by the surrogate contact model using the CoPbased approach were compared with those calculated by the EF contact model using both the CoPbased approach and a pr eviously published elementbased approach (Fregly et al., 2005; Zhao et al., 2008) to verify that the surrogate model wear volume predictions were accurate. Surrogate Contact Model Application As a practical application of the proposed surrogate co ntact modeling approach, we performed two types of sensitivity analyses to investigate how machine setup issues affect predicted wear volume. The first was a global sensitivity analysis to investigate the effect of large variations in component placements within the ranges defined in Figure 22A (i.e., mm/deg). Specifically, we generated 400 combin ations of component placements using 20 evenly spaced values of femoral flexion axis SI location and tibial insert AP slope. For each combination, the surrogatebased contact model was used to perform a dynamic simulation and predict wear volume. Percent changes in predicte d wear volume were calculated relative to the maximum predicted value from all global sensitiv ity analysis results. The second category was a local sensitivity analysis to investigate the effect of small variations in nominal component placements and motion and load inputs by using M onte Carlo techniques (Figure 22B). For the local sensitivity analysis, allowa ble variations in motion and load input curves were defined using experimental observations from eight different implant designs tested in a Stanmore simulator machine (DesJard ins et al., 2000). Deviations of th e experimental input curves from the ISO standard curves were calculated and analyzed using principal component analysis (Daffertshofer et al., 2004), w ith the first two principal comp onents being sufficient to account PAGE 24 24 for over 98% of the variability in each type of curve. These two components were used to generate new deviation curves by selecting com ponent weights as uniformly distributed random numbers within the bounds of the experimental cu rves. These deviation curves were added to the ISO standard curves to create new input motion and load curves as needed for the Monte Carlo analyses. To distinguish between the effects of small errors in motion and load inputs and small errors in component placements, the local sensit ivity analysis performed four types of Monte Carlo analyses for each of the nine nominal comp onent placements for a total of 36 Monte Carlo analyses in all. Each type of Monte Carlo analysis varied input profiles and component placements separately or together. The first type varied motion and load profiles within 100% of their allowable ranges and also va ried component placements within 1 mm/deg. The second type was similar but used smaller variations of only 10% of the allowable range for motion and load profiles and 0.1 mm/deg for component placements The third type varied only motion and load profiles within 100% of their allowable ranges and impos ed no variations on component placements. Finally, the fourth type va ried component placements within 1 mm/deg but imposed no variations on motion and load profiles. For all types, percent changes in predicted wear volume were calculated relative to values obtained from the nominal component placements using the ISO standard inputs. Each Monte Carlo analysis was performed on a 3.4 GHz Pentium IV PC and involved at least 1000 dynamic contact simulations. The convergence criterion for each analysis was met when the mean and coefficient of variance (i.e., 100*standard deviation/mean) for the last 10% of the wear predictions were w ithin 2% of the final mean and coefficient of variance (Val eroCuevas et al., 2003). PAGE 25 25 Results For the computational evaluation, the dynamic simulations performed with the surrogate contact model closely reproduced the planar moti ons, contact forces, and wear volumes predicted with the EF contact model for all nine component placements. On average, the RMSE for planar motions and contact forces were less than 0.125 mm and 4 N, respectively (Figure 23). While each dynamic simulation performed with the EF contact model required approximately 13 minutes of CPU time, each simulation performed using the surrogate contact model required 13 seconds. For both the surrogate contact model and the EF contact model, the CoPbased approach for predicting wear vol umes reproduced the elementbase d EF results to within 0.5% error (Table 21). For the practical application, the global and lo cal sensitivity analyses revealed that large variations in component placements relative to the original locations and small variations in input motion relative to the ISO standard had a significant affect on predicted wear volume. For the global sensitivity analysis, predicted wear volum e was sensitive to large variations in the SI location of the femoral flexion axis (maximum ch ange of 26% when AP slope of tibial insert held fixed), the AP slope of the tibial insert (maximum change of 13% when SI location of femoral flexion axis held fixed), or both para meters simultaneously (maximum change of 33%) (Figure 24). For the local sensitivity analysis wear volume ranges were larger (mean of 18 mm3 compared to 2 mm3) and 50% percentile wear volume re sults were lower (mean of 23 mm3 compared to 31 mm3) when motion and load inputs and component placements were varied within 100% rather than 10% of their allowable ranges (Figure 25A and 25B). When inputs and placements were varied separately within 100% of their allowable ranges, predicted wear volume was much more sensitive to motion and load i nputs than to component placements (maximum change of 64% compared to 5%; Figure 25C and 25D). Finally, when an additional global PAGE 26 26 sensitivity analysis was performed with the or iginal component placements to separate the influence of the different inputs, predicted wear volume varied approximately linearly with small changes in each input curve but was more sensitive to changes in flexionextension motion (maximum change of 47%; Figure 26) than to ch anges in input loads (maximum change of 20% for Fx and 14% for Fy; Figure 26A and B). Each Monte Carlo analysis performed with the surrogate contact model required 4 hours of CP U time compared to an estimated 230 hours with the EF contact model. Discussion Computational speed is a major limiting factor for performing sensitivity and optimization studies of total knee replacement designs. Unlik e constraintbased joint models, knee joint contact models require repeated geometry evaluati ons that consume the vast majority of the CPU time in a dynamic contact simulation (Bei a nd Fregly, 2004). This paper has presented a surrogate modeling approach for performing effi cient dynamic contact simulations of human joints. The surrogate contact model was developed to replace the original EF contact model to avoid timeconsuming repeated evaluation of the surface geometry during dynamic simulation. The surrogate contact model produced dynamic simulation and wear results that accurately matched those from the EF contact model but with an order of magnitude improvement in computational speed. This speed improvement would have been even more significant if a finite element contact model had been replaced. The 36 Monte Carlo analyses demonstrated how analyses that would have been impractical (o r at a minimum would have required extensive parallel processing) can be eas ily achieved on a single computer with the surrogate contact model. The sensitivity results also have practic al value, suggesting that wear volume generated by simulator machines may be sensitive to large variations in component placements within the machine as well as small variations in the flexionextension motion input to the machine. PAGE 27 27 The primary benefit of using surrogate contact models in dynamic simulations is improved computational efficiency. For elastic foundation or finite element contact models, the computation time per dynamic simulation is largel y determined by the number of deformable contacts in the model. Thus, a musculoskeleta l computer model utilizing deformable contact models for multiple joints could easily require hours or days of CPU time to complete a single dynamic simulation. Repeated dynamic simulations as part of a sensitiv ity or optimization study would become impractical. Although it took 3 hours of CPU time on generating 1500 sample points, the use of the surrogate contact model in place of the EF model reduced Monte Carlo analysis time from 230 hours to 4 hours, a near ly 99% reduction. By using surrogate contact models instead, one could eliminate the computa tional cost of the contact solver as the limiting factor. The larger the number of deformable contacts in a multibody dynamic model, the greater the anticipated benefit from us ing surrogate contact models. While we have shown that surrogate contact models can be beneficial for sensitivity studies, the greatest computational benefit is likely to occur for optimization studies. For stochastic sensitivity studies, statisticallybased alternatives to Monte Carlo analyses exist that require fewer repeated simulations. For example, Laz et al. (2006) used the mean value method (Wu et al., 1990) to perform probabilistic elas tic foundation simulations of a TKR design and reduced the number of repeated simulations by a factor of four compared to a traditional Monte Carlo approach. Since the CPU time per simulati on was 6 minutes, use of a surrogate contact model could still reduce the overall computation time by an additional factor of 28, assuming we could achieve comparable speed improvements for the 3D situation. In contrast, no good methods exist for reducing the number of repeated simulations in optimization studies. For those PAGE 28 28 situations, the factor of 60 reduction in computation time (i.e ., 13 minutes to 13 seconds) could mean the difference between an impossibl e and an achievable optimization study. Other than computational efficiency, surr ogate contact models possess the additional benefits of being adaptable and flexible. First, the surrogateb ased contact modeling approach presented here is not limited to the knee joint. Theoretically, the surrogatebased approach could be applied to any joint contact model whose i nputoutput relationships can be sampled. The modified HQ sampling method should be suitable for any joint contact m odel once the sensitive directions are identified. Therefor e, neither type of joint (e.g., hi p or knee or ankle) nor type of contact model (e.g., EF or finite element) nor t ype of material model (e .g., linear or nonlinear) is critical for developing a surrogate contact model. The benefit of using an EF contact model to generate the required sample points was that it is relatively fast computationally for performing repeated static analyses. Finite element models would require more time for each static analysis, plus the inputs are usually applied displacements rather than loads. For the finite element situation, it could be necessary to curve fit the forcedisplacement curve at each sample point to produce new points at fixed force increments. Sec ond, the surrogate contact model represents the real surface geometry implicitly within the desi gn space. Since the surrogate modeling approach uses the real surface geometry to generate the sampling points, no idealization of complex surface geometry is required. Third, the approa ch is not limited to a single surrogate model fitting method. While we used Kriging in the pres ent study, any of the f itting methods mentioned previously could be investigated as well (Jin et al., 2001; Simp son et al., 2001; Queipo et al., 2005; Wang and Shan, 2007). Though polynominal re sponse surfaces are the most commonly used surrogate model, Kriging has the advantage of interpolating the sample point outputs, which is appealing for the determinis tic computer models (Sacks et al., 1989; Simpson et al., 2001). PAGE 29 29 The global sensitivity results demonstrate th e importance of devel oping a standardized method for defining the location of the femoral fl exion axis. For a given tibial slope, predicted wear volume varied approximately linearly with changes in the superiori nferior location of the femoral flexion axis, increasing as the axis was translated superiorly (Figure 24). Superior translation of the flexion axis increases the dist ance to the CoP, thereby increasing the magnitude of the slip velocity in Equation 24. Thus, to minimize wear for any given implant design, one can simply locate the femoral flexion axis clos e to the articular surfaces. In contrast, the predicted wear volume varied approximately quadratically with changes in tibial slope, increasing as the insert was tilted anteriorly or posteriorly away from its neutral position (Figure 24). Increasing the tibial anterior or posterior slope increases AP tibial sliding distance (Bai et al., 2000), thereby leading to an increase in predicted wear volu me (Blunn et al., 1991; Kawanabe et al., 2001). In contrast, the local sensitivity results dem onstrate the importance of controlling machine inputs closely, especially the femoral flexi on angle. The Monte Carlos analyses and the additional global sensitivity analysis revealed wh en small variations were imposed on the input motion, loads, and component placements, only vari ations in the femoral flexion angle resulted in large changes in predicted wear volume. When the flexion angle was varied within 100% of its allowable range, the drop in 50th percentile wear volume resu lts (Figure 25A and 25C) relative to the ISO standa rd inputs (i.e., 50% percentile results in Figure 25B and 25D) was due to the use of imposed variations that always reduced the amplit ude of the curve (DesJardins et al., 2000). This finding is consistent with previo us studies showing that predicted wear volume decreases when the amplitude of the input curves is reduced (Barnett et al., 2001; D'Lima et al., 2001; Johnson et al., 2001). Thus, another way to minimize wear for any given implant design is PAGE 30 30 to use a control system that systematically unde rshoots the peaks in the input flexionextension curve. Combining this approach with an inferior ly located femoral flexion axis could cut wear volume by more than 50% compared to the orig inal component placement with ISO standard inputs. The development of a standard for accep table variations in input flexionextension motion should be c onsidered as well. Despite its computational advantages, the pr oposed surrogate contac t modeling approach possesses several important limita tions. The biggest one is that new surrogate contact models must be generated if implant geometry, insert th ickness, or material prope rties are changed, at least for the current formulation. The geometry limitation makes it impossible to use surrogate contact models for progressive wear simula tions, where the surface geometry is changed gradually over an iterative sequence of wear si mulations. However, this limitation is not serious if only wear volume is of interest, as recent progressive wear simulation studies have reported that predicted wear volume (but not wear area or depth) is relatively insensitive to whether or not the surface geometry is changed gradually (Knigh t et al., 2007; Zhao et al., 2008). Theoretically, changes in insert thickness or material m odel parameter values could be accommodated by adding more design variables to the surrogate mo del fitting process (e.g ., make insert thickness an additional surrogate model input). For contact forces, a changing value of Youngs modulus (E) in a linear elastic material model can be accommodated without adding another design variable. If all sample points are generated using a value of E = 1, then the contact forces output by the surrogate model can be scaled by the desired value of E (Lin et al., 2006). A number of other limitations exist as well. Fi rst, the current surrogate contact model is restricted to sagittal plane motion, where cont act force is dominated by SI translation of the femoral component. For 3D simulations, contact force will be highly sensitive to changes in PAGE 31 31 varusvalgus rotation as well (Fregly et al., 2008). In addition to an increased number of sensitive directions, the number of sample points will also need to increase. However, the increased computational cost for evaluating more sample poi nts is paid up front since the sampling process is only performed once. Second, the calculation of Fx is dependent on the value of Fy computed from Equation 22. Thus, the quality of the Fy prediction will directly impact the Fx calculation. The current approach works fine for sagittal pl ane models with a single sensitive direction, but an independent calculation of Fx may be needed for 3D situations where more than one sensitive direction is present. Third, the CoP approach for calculating w ear volume may not work well for some types of 3D motion. For example, for a sp here spinning about a vertical axis while in contact with a plane, the CoP would be located at the bottom of the sphere where the slip velocity is zero. Thus, no wear volume would be calculated, even t hough the elementbased approach would calculate a nonzero wear volume based on the pressures and slip velocities throughout the contact area. Consequently, the CoP approach will work well for analysis of TKR designs undergoing 3D motion only if the amount of internalexternal rotation is small. Fourth, our Monte Carlo results are limited to cruciate retaining designs with a nonconformal tibial insert. Whether or not our conclusions are genera lizable to other implant designs and simulator machines will require further investigation. In summary, this chapter has demonstrat ed that surrogate contact modeling can significantly improve the computational speed of dynamic contact simulations and is appropriate for 2D sensitivity and optimization studies incorporating deformable contact. Refinement of the current 2D approach and extension to 3D problems are topics of ongoing research. Once surrogate contact models are deve loped for 3D situations, it shoul d be possible to create dynamic musculoskeletal models with multiple deformable joints that can be simulated in a short amount PAGE 32 32 of CPU time. Such models should improve our ability to analyze new knee implant designs rapidly and to predict muscle forces across join ts where contact forces may have an important stabilizing effect. Figure 21. Overview of surroga te contact model creation and implementation process. A) 300 points (black dots) are sampled in the (X, ) design space. B) Five static analyses using five axial loads are performed at each sample point location using an elastic foundation contact model. C) Generalized Hert zian contact theory is used to model Fy as a function of Y while Fx is modeled as a linear function of Fy. D) The values of 10NY, ky, ny, and ratio from the 300 sample points are fitted as functions of X and using Kriging. E) During a dynamic si mulation, the surrogate contact model calculates Fy and Fx from four Kriging models given the current values X, Y, and PAGE 33 33 Table 21. Comparison between the wear volume predictions (mm3) made with the elastic foundation contact model and the surrogate contact model for each of the nine nominal component placements. Each pr ediction was made for 5 million motion cycles. Both contact models can predict wear volume using a simple center of pressurebased approach, while only the elastic foundation contact model can make predictions using a more detailed elementbased approach. Elastic foundation contact m odel Surrogate contact model Nominal Placement Elementbased CoPbased CoP based A 37.07 37.05 36.98 B 34.03 34.02 33.91 C 39.04 39.03 38.87 D 32.91 32.89 32.88 E 30.78 30.76 30.77 F 34.20 34.19 34.14 G 29.37 29.36 29.32 H 26.45 26.44 26.38 I 29.37 29.35 29.31 Figure 22. Overview of globa l and local sensitivity analys es of predicted wear volume performed with the surrogate contact model. A) Large variations in nominal component placements (labeled A through I) us ed for the global sensitivity analysis. B) Small variations in component placemen ts, input loads, and input motion applied about each nominal component placement us ed for the local sensitivity analysis involving Monte Carlo methods PAGE 34 34 Figure 23. Rootmeansquare errors in joint motions and contact forces for each nominal component placement. Errors are computed by comparing dynamic simulation results generated with the elastic foundation contac t model and the surrogate contact model. Figure 24. Predicted wear volume as a function of large variations in component placements as calculated from 400 dynamic simulations performed with the surrogate contact model. Wear increases linearly when the femo ral flexion axis is translated superiorly and quadratically when the tibial slope is increased anteriorly or posteriorly. PAGE 35 35 Figure 25. Box plot distributi on of predicted wear volume ge nerated by four Monte Carlo analyses using the surrogate contact mode l. Each box has (from bottom to top) one whisker at the 10th percentile, three lines at the 25th, 50th, and 75th percentile, and another whisker at the 90th percentile. Outliers are indicated by black crosses located beyond the ends of the whiskers. For the first and second Monte Carlo analyses, motion and load inputs and component placem ents are varied together within A) 100% or B) 10% of their maximum ranges, respectively. For the third and fourth Monte Carlo analyses, C) motion and load inputs or D) component placements are varied separately within 100% of their maximum ranges, respectively. PAGE 36 36 Figure 26. Predicted wear volume as a function of variations in A) anteriorposterior load and femoral flexion angle and B) axial load a nd femoral flexion angle. Both plots were created by performing dynamic simulations with the surrogate model and changing the weight of the first principal component for each input profile variation away from the ISO standard. Each weight was normali zed to be between 0 and 1. The cross marker (x) represents the ISO standard input with no superimposed input profile variations (i.e., weights of 0) while the circle marker (o) represents the ISO standard input with maximum superimposed input pr ofile variations (i .e., weights of 1) PAGE 37 37 CHAPTER 3 THREEDIMENSIONAL SURROGATE CONTACT MODELS FOR COMPUT ATIONALLYEFFICIENT MU LTIBODY DYNAMIC SIMULATIONS Introduction Multibody dynamic simulation of systems experi encing one or more contacts is valuable for engineering applications ranging from the de sign of industrial machines and mechanisms to the analysis of human joints. Computational cont act models used to perform such simulations generally fall into two categories: rigid body contact models and deformable body contact models (Sharf and Zhang, 2006). Rigid body cont act models use unilateral constraints to maintain contact between opposing surfaces and are highly efficient for predicting motion (Glocker, 1999; Piazza and Delp, 1999; Stewart, 2000; Porta et al., 2005; Brogliato et al., 2002). However, they do not permit calculation of unique contact forces under sta tically indeterminate conditions, such as when contact occurs in tw o or more places between the same pair of contacting bodies (Cheng et al. 1990). Furthermore, they normally require different computational schemes to simulate continuous contact versus intermittent impact. In contrast, deformable body contact models, su ch as finite element models (Gradin and Cardona, 2001; Halloran et al., 2005a; Gerstmayr and Schoberl, 2006) and elastic foundation models (Fregly et al., 2003; Hippmann et al., 2004; Santamaria et al., 2006; Moran et al., 2008), provide a unified approach for simulating continuous contact versus intermittent impact and can calculate unique contact forces in statically indeterminate situations. Furthermore, they permit calculation of distributed contact pressures and strains across the opposing surfaces. These benefits are achieved by using elas tic contact theory rather than constraints to determine the forces and torques exerted by the contacting bodies on each other. Unfortunately, the high computational cost of deformable body contact anal yses, particularly due to extensive geometry evaluations, is a significant lim iting factor to their use in multibody dynamic simulations (Bei PAGE 38 38 and Fregly, 2004). Thus, no contact m odeling approach curren tly exists that prov ides the benefits of deformable body contact methods and the com putational speed of rigid body contact methods. In other engineering discip lines, surrogate modeling approaches have been used successfully to overcome similar com putational challenges (e.g., Roux et al., 1996; Liu et al., 2000; Cox et al., 2001; Queipo et al., 2002; Wang et al., 2007). In general terms, these approaches involve replacing a computationally costly model with a computationally cheap model constructed using data points sampled from the original model. Once constructed, the cheap surrogate model is used in place of the costly original model to eliminate computational bottlenecks. Until recently, only a small number of studies have applied surrogate modeling techniques to cont act problems (Chang et al., 1999; Bouzid and Champliaud, 2003; Lin et al., 2006; Lin et al., 2008). In Chapter 2, we have presente d a 2D surrogatebased dynamic contact simulation by using the knowledge of phys ics. However, such technique was suitable only when contact force is dominated by the tran slation normal to the co ntact surface. For 3D knee contact simulations, contact force will be highly sensitive to changes in varusvalgus rotation as well (Fregly et al., 2008). None of the existing techniques is suitable for threedimensional elastic contact problems wh ere the contacting bodies undergo large relative displacements. This study proposes a novel surrogate modeli ng approach for performing computationallyefficient 3D elastic contact analyses within multibody dynamic simulations. The approach uses two key concepts sensitive directions and a reasonable design space to address the unique challenges involved in applying surrogate modeli ng techniques to elastic contact problems. Sensitive directions define which degrees of fr eedom in the original model should be sampled using applied loads rather than displacements while a reasonable design space screens out PAGE 39 39 infeasible sample points before they are evaluate d with the original model. The performance and accuracy of the proposed approach are evaluated by performing thousands of computational wear simulations for a TKR tested in a Stanmore knee simulator machine. Methods Surrogate Contact Model Development Special surrogate model creation techniques are required to apply surrogate modeling methods to elastic contact problems. From a big picture perspective, we treat surrogate contact model creation as a nested process. First, a coar se surrogate model is c onstructed by evaluating a limited number of sample points with the original contact model. Next, the coarse surrogate model is used to rapidly evaluate the complete set of sample points and screen out infeasible ones. Finally, a fine surrogate contact model is constructed by evaluating all feasible sample points with the original contact model. This approach minimizes the number computationally costly analyses performed by the original co ntact model during the surrogate model creation process. The remainder of this section describes our special surrogate mode l creation process in detail. Eight key assumptions pr ovide the foundation for the subsequent theo retical development: 1. Threedimensional elastic cont act is being modeled between a master body and a slave body possessing general surface geometry and material properties. 2. The pose of the slave body relativ e to the master body is defined by six pose parameters (i.e., three translations x, y, z along with three rotations using a specified rotation sequence). 3. For each combination of six pose parameters there exists a unique combination of six contact loads (i.e., three forces Fx, Fy, Fz and three torques T, T, T) experienced by the two contacting bodies in an equal and opposite sense. 4. For any pose, the 6 x 6 matrix defining the sens itivity of the contact loads to pose parameter changes is diagonally dominant (e.g., Fy is most sensitive to changes in y). PAGE 40 40 5. A deformable body contact model (henceforth called the original model) is the computationally costly contact model to be replaced. 6. A sample point is any combination of six pose and load inputs to the original model (e.g., x, Fy, z, ) for which the six corresponding load and pose outputs are desired (e.g., Fx, y, Fz, T, T, T). 7. Given a large number of sample points, repeated static analys es can be performed with the original model to generate inputoutput rela tionships to be fitted by the computationally cheap surrogate contact model. 8. A surrogate contact model is actually a collectio n of six separate surrogate models one for each of the six contact loads. Based on these foundational assumptions, we ha ve devised a threestep surrogate model development process to address the unique characteristics of elastic contact models. Step 1) Identify a nonstandard samplin g method using the concept of sensitive directions. As noted above, the surrogate model creati on process involves performing repeated computational experiments with the original cont act model to generate in putoutput pairs to be fitted by the surrogate contact model. The combinations of inputs to be analyzed are selected using design of experiments (DOE), which employs statistical methods to scatter sample points uniformly throughout the bounded sixdimensional design space. Within the context of a multibody dynamic simulation, the desired inputs to a 3D elastic contact model are the six pose parameters and the desired outputs are the six co ntact loads (Bei and Fr egly, 2004). Thus, using a traditional sampling method, the sample points would be combinations of the six pose parameters x, y, z, with an upper and lower bound place d on each one. The problem with this approach is that physically realistic contact occurs over a thin hypersurface in sixdimensional pose parameter design space. Cons equently, many of the selected sample points will correspond to situations where the opposing surfaces are either out of contact or deeply penetrating. Inclusion of physically unrealisti c sample points in the surrogate model fitting process is undesirable since it reduces the accu racy of the resulting su rrogate contact model. PAGE 41 41 To resolve this issue, we propose the concep t of sensitive direct ions to modify the traditional sample point definition so as to avoid physically unrealistic situations. When the master and slave body are in contact, some contact loads will be highly sensitive to changes in their corresponding pose parameters, while othe rs will be insensitive. To quantify these sensitivities, we calculate six central difference de rivatives and collect the results in a sensitivity vector s (Equation 31) //////xyzFxFyFzTTT s (31) Due to differences in units, sensitive directions are determined separately for translations (first three entries of s ) and rotations (last three entries of s ). If one translational (or rotational) derivative in s is significantly larger th an the other two, the corres ponding direction is deemed a sensitive direction. If two transl ational (or rotational) derivatives are significantly larger than the third, then two sensitive directions exist. For so me situations (e.g., conformal contact between a sphere and a spherical cup of slightly larger radius), all three translational (or rotational) derivatives may be of comparable magnitude, in which case knowledge of the physical situation must be used to determine whether three or zero sensitive directions exist. At least one translational sensitive direction will always exist corresponding to an approximate contact normal direction, and two or more sensitive directions may exist depending on the geometry of the contacting bodies. Once sensitive directions have been identified, the definition of a sample point is modified such that pose parameters for sensitive direc tions are replaced by thei r corresponding contact loads. For example, if y is identified as a sensitive direction, then y would be replaced with yF in the sample point definition. With this nontraditional sampling method, sample points become PAGE 42 42 combinations of pose parameters and contact loads (e.g., x, yF z ), with an upper and lower bound placed on each pose parameter and contact load. Given this modified sample point definition, we generate a specified number of sample points n using a DOE approach. Common approaches include optimal Latin hypercube sampling (Mckay, 1979), Hammersley quasirandom (HQ) sampling (Hammersley, 1960), and facecentered central composite design (Myers and Montgomery, 1995). We choose to use HQ sampling for two reasons. First, HQ sampling provi des better uniformity properties than do other sampling techniques for sample points generated within a multidimensional hypercube (Kalagnanam and Diwekar, 1997; Diwekar 2003). Second, HQ sample points are generated sequentially rather than simultaneously. Thus, the first m sample points from a larger set of n sample points ( n > m ) will always be approximately equidistant from one another and can be used as a subsample. Furthermore, all existing sa mple points can be kept if new sample points need to be added. The repeated static analyses are performed with the original contact model to find the physically realistic poses for the given sample points. During each static analysis, sensitive directions are free to equilibrate under the specif ied applied loads, while insensitive directions are constrained to the specified pose paramete r values. For each static analysis, the resulting loads in sensitive directions are compared to th eir corresponding applied values. A sample point is chosen only when the maximum percent error of 0.1% is achieved. By the end of static analyses, we ensure that all chosen points fr ee from situations like out of contact or deeply penetrating. Step 2) Filter out infeasible sample points before evaluation using a coarse surrogate contact model and the concept of a reasonable design space Though physically unrealistic PAGE 43 43 sample points are avoided by using the modified sample point definiti on described above, the resulting sample points may still be outside the envelope of a realistic dynamic simulation (i.e., infeasible sample points). In addition to placing physically re alistic upper and lower bounds on the pose parameter and contact load inputs (e.g., x yF z ), we place physically realistic upper and lower bounds on the corresponding co ntact load and pose parameter outputs (e.g., x F y zF T T T ). These bounds are estimated by performing a nominal dynamic simulation with the original contact model and expanding the observed range of each pose parameter and contact load by some specified percentage p. Sample points where the corresponding contact model outputs are ou tside their allowable bound s are deemed infeasible. Similar to the inclusion of unrealistic sample points, inclusion of infeasib le sample points in the surrogate model fitting process adversely affects the accuracy of the resulting surrogate contact model. Furthermore, evaluation of infeasible sample points by the original contact model wastes a large amount of CPU time performing computationally costly sta tic analyses whose results will be discarded. To resolve this issue, we utilize the concept of a reasonable design space to screen out infeasible sample points before evaluating them w ith the original contact model. Generally, the design space can be refined by either directly reducing the rang e of design variables or by using some computationally inexpensive tool to define irregular design space bo undaries (Mack et al., 2007). For the latter case, the reasonable design space approach has been successfully applied to different studies for developing efficient and accurate surrogatebased optimization (Kaufman et al ., 1996; Roux et a l., 1999; Hosder et al., 2001). To estimate the reasonable design space, we select the first m sample points from the complete set of n sample points generated by HQ sampling. These m sample points are the foundations of reasonable design space approach. First, a static analysis mentioned in previous step is performed with the origin al contact model for each PAGE 44 44 of the m sample points, and physically unrealistic po ints are identified an d eliminated. Second, six coarse surrogate models are created based on fitting six static analysis outputs (e.g., x F y zF T T T ) of the remaining points (< m ) as functions static analysis inputs (e.g., x yF z ). Third, six coarse surrogate models are eval uated repeatedly to determine the feasibility of the nm sample points not used to construct it. Fourth, the final set of all feasible sample points (<< n ) is used to construct a fi ne surrogate contact model. Since a large number of the n m sample points will turn out to be infeasible, sc reening these points with the coarse surrogate contact model provides significant computational savings. Step 3) Construct a fine surrogate contac t model using the feasible subset of the original sample points Fitting a fine surrogate contact mode l to the feasible sample points is complicated by the fact that contact model outp uts and inputs are different for fitting versus sampling. To illustrate this issue, consider an elastic contact model whose only sensitive direction is y translation. During sampling of the original contact model, outputs would be x F y zF T T T while inputs would be x yF z During fitting of the surrogate contact model, outputs would be x F yF zF T T T and inputs would be x y z , as required by the dynamic simulation. Altered in terpretation of outputs and inputs between sampling and fitting is not problematic, since at the end of the sampling process, we are left with a table containing six pose parameters a nd six contact loads for each sample point. With this issue in mind, the final surrogate model creation process proceeds as follows. Repeated static analyses are pe rformed with the original contact model for all feasible sample points not yet evaluated. Since th e coarse surrogate models may have wrongly identified some infeasible sample points as feasible, sample points whose outputs are outside their allowable bounds are discarded. The final set of all feasible sample points (<< n ) is used to construct the PAGE 45 45 final surrogate contact model. C ontact loads corresponding to sensit ive directions are fitted as a function of the six pose parameters (e.g., (,,,,,)yFfxyz ), while contact loads corresponding to insensitive directions are fitt ed as a function of contact loads in sensitive directions and pose parameters in insensitive ones (e.g., (,,,,,)xyFfxFz ). This approach keeps the fitting process as close to the sampli ng process as possible wi thout solving a nonlinear rootfinding problem for the sensitiv e directions (e.g., given (,,,,,)yyfxFz find yF such that (,,,,,)0ycfxFzy where c y is the current value of y in the dynamic simulation). Actual fitting of the six inputoutput relationships is performed using Kriging (Krige, 1951), as our preliminary studies reveal ed that Kriging produces more accurate dynamic simulation results than do a variety of other surrogate model fitting methods (e.g., polynomial response surface and support vect or regression). Once a Krigingbased surrogate contact model has been constructed, it is used to calculate six contact loads at each time instant during a dynamic simulation given the six pose parameters for the contacting bodies. Surrogate Contact Model Evaluation The proposed surrogate modeling approach was evaluated by pe rforming Monte Carlo analyses and forward dynamic simulati ons using a multibody dynamic model of a cruciateretaining commercial knee implant (Depuy Orthopedics Warsaw, IN) mounted in a Stanmore knee simulator machine (Walker et al., 1997). The model possessed six degrees of freedom (DOFs) relative to ground: tibial anterio rposterior translation x, tibial medi allateral translation z, tibial internalexternal rotation femoral superiorinferior translation y, femoral varusvalgus rotation and femoral flexionextension (Figure 31) Similar to an actual Stanmore machine, was motion controlled, x, y, and were load controlled, and z and were left free, where all controlled quantities used ISO standard input curves (DesJardins et al. 2000). PAGE 46 46 A 3D elastic foundation (EF) contact m odel (An et al., 1990; Blankevoort et al., 1991; Li et al., 1997; Pandy et al., 1998; Fregly et al., 2003) of th e knee implant was used to calculate contact loads between the femoral component and tibial insert given the six pose parameters for the contacting bodies. The EF model utilized line ar material properties (Youngs modulus = 463 MPa, Poissons ratio = 0.46; Bei and Freg ly, 2004) and surface geometry taken from manufacturer computeraided de sign models, with contact occu rring on the medial and lateral sides of the tibial inse rt. Symbolic dynamics equations for th e system were derived using Kanes method (Kane and Levinson, 1985) and Autolev (OnLine Dynamics, Sunnyvale, CA). These equations and the EF contact model were incorpor ated into a Matlab program (The Mathworks, Natick, MA) that was used to perform forwar d dynamic simulations with Matlabs stiff numerical integrator ode15s. The surrogate modeling approach presented is a general algorithm, and it can be applied to any number of elastic contacts. To demonstr ate this ability during the surrogate model evaluation, two respective surroga te contact models were constr ucted to replace the medial and lateral contacts occurring on the tibial insert. These two models were created using threestep process outlined above. Step 1) Identify a nonstandard sampling meth od using the concept of sensitive directions. With the implant component s in a nominal anatomic pose and both sides barely touching, we calculated the sensitivity vector s defined in Equation 31 and identified two sensitive directions. Specifically, contact force yF and contact torque T were found to be highly sensitive to small variations in superiorinf erior translation y and varusvalgus rotation respectively. Thus, sample points were defined as x, yF, z, T Upper and lower bounds for the sample points were determined by performing 16 dynamic contact simulations with the EF PAGE 47 47 contact model using motion and load input curves at the extremes of thei r allowable variations for the subsequent Monte Carlo analyses After expanding the resulting bounds by p = 50%, we generated n = 2000 sample points using the HQ samp ling method (Figure 32A), where the values for p and n were chosen based on previous experien ce with surrogate contact models (Lin et al., 2008). Step 2) Filter out infeasible sample points before evaluation using the concept of a reasonable design space To estimate which sample points were within the reasonable design space, we selected m = 500 sample points from the complete set of 2000. We then performed 500 static analyses with the EF contact model to calcula te the corresponding outputs x F, y, zF, T T where each static analysis required appr oximately 24 seconds of CPU time on a 3 GHz Pentium IV PC (Figure 32B). 18 sample points were identif ied as unrealistic, leaving 482 sample points for constructing coarse su rrogate models by fitting six outputs (i.e., x F, y, zF, T and T ) as functions of sample inputs (Figure 32C). After performing a cross validation analysis with these 482 sample points, we sele cted a cubic polynomial and Gaussian correlation function to construct all surrogate models using the DACE Toolbox fo r Matlab (Lophaven et al., 2002). Since each sample points contained informati on of medial and latera l side contact loads, two sets of six coarse surrogate models were created to predict their corresponding values for the remaining 1500 sample points (Figure 32D). 93 4 (937) of these points were found to be infeasible for medial (lateral) contact, leaving 56 6 (563) additional sample points to be evaluated with the EF contact model (Figure 32E). Step 3) Construct the final surrogate contac t model using the feas ible subset of the original sample points Considering the overlap between th ree groups of sample points (566, 563, and first 500 sample points), only 560 additional static analyses were performed with the EF PAGE 48 48 contact model for the feasible sample points identified in the previous step (Figure 32F). The CPU time required for 1060 (=500+560) static analyses was 6 hours. The results of these static analyses showed that 500 (560) out of 566 (563) sample points were feasible for creating the medial (lateral) part of the final surrogate contact model. Each part fitted six contact loads as a function of the following inputs (Equation 32, Figure 32G): (,,,,,) (,,,,,) (,,,,,) (,,,,,) (,,,,,) (,,,,,)xy y zy y yFfxFzT Ffxyz FfxFzT Tfxyz TfxFzT TfxFzT (32) Since x, y, z, are the inputs, Equation 32 indicates that contact loads in the sensitive directions (i.e., yF and T ) must be calculated before contact loads in the insensitive directions. The final surrogate contact model wa s composed of two sets of these six Kriging models and was incorporated into the multi body dynamic model in place of the EF contact model (Figure 32H). During the dynamic simulation, the total contact loads are generated by summarizing Krigingbased medial and lateral contact loads. Six additional Kriging models were created to fit medial and lateral center of pressure (CoP) location (3 components per side) as functions of the samp le inputs. The CoP predictions and previously calculated contact forces were us ed to calculate medial and lateral wear volume. Details of the process are given by Lin et al ( 2008). Briefly, Archards wear law (Equation 24) with a wear factor of 1 107 mm3/Nm (Fisher et al., 1994) was used to predict wear volume on each side from the time history of contact forc e and CoP slip velocity along with the time increment for numerical integration. Once the Co P location was predicted for a given pose, the first time derivative of the pose parameters a nd the kinematic concept of two points fixed on a PAGE 49 49 rigid body were used to calculate the correspond ing CoP slip velocity (Kane and Levinson, 1985). To evaluate the performance of surrogateb ased dynamic simulation, we performed five Monte Carlo analyses to investigate how realisti c variations in motion a nd load inputs affect predicted wear volume. The first four Monte Carlo analyses varied the three load inputs and one motion input separately while the fifth analysis varied all input profiles together. Principal component analysis was used to create realistic variations of the input curves based on experimentally observed variations for eight di fferent implant designs tested in a Stanmore simulator machine (DesJardins et al. 2000; Lin et al., 2008). Each new input curve was generated by selecting weights between 0 and 1 for the fi rst two principal components, which captured 98% of the variability in each type of curve. Each Monte Carlo analysis was performed on a 3 GHz Pentium IV PC and involved at least 1000 forward dynamic simulations incorporating the surrogate contact model. The convergence criteri on for each analysis was met when the mean and coefficient of variation (i.e., 100*standard deviation/mean) for the last 10% of the wear predictions were within 2% of the final mean and coefficient of variation (Fishman 1996; ValeroCuevas et al., 2003). To evaluate surrogate contact model accu racy, we performed eleven dynamic wear simulations with both the surrogate contact model and the EF cont act model used to create it. A nominal dynamic simulation was first performed using the ISO standard input curves and the simulation results from two models compared qua ntitatively. We then performed ten additional EFbased and surrogatebased dynamic simulations to evaluate the extremes in input curve variations used during the Monte Carlo analyses Specifically, for each M onte Carlo analysis, we performed two dynamic simulations with each contact model where each input curve was PAGE 50 50 generated using either the maximum or minimum sum of the component weights. The ten pairs of results were compared quantitatively by calcu lating root mean square errors (RMSE) and maximum absolute errors (MAE) for predicted pose parameters, medial and lateral contact loads, and wear volumes. The CoPbased approach fo r calculating wear volume was also evaluated by calculating wear volume with the EF model using both an element ba sed approach (Fregly et al., 2005; Zhao et al., 2008) and the CoPbased approach (Lin et al., 2008). Results For the Monte Carlo analyses, the results reveal ed small variations in input motion relative to the ISO standard had different effects on predicted wear volu me (Figure 33). The predicted wear volume was more sensitive to changes in fl exionextension motion than to changes in input loads. The most variation was obt ained when all of the input curv es were varied simultaneously. Each Monte Carlo analysis perf ormed with the surrogate contac t model required 1.4 hours of CPU time compared to an estimated 284 hours w ith the EF contact model. As for a single dynamic simulation, the CPU time required for EF and surrogate contact model was 17 minutes and 5 seconds, respectively. The comparison between multiple EFbased and surrogatebased dynamic simulations demonstrated that the surrogate contact mode l accurately reproduced the contact kinematics, contact loads, and wear volumes predicted with the EF contact model. The nominal dynamic simulation results showed similar behavior betw een the surrogate and EF contact model (Figure 34 to 36). Similar to the nominal dynamic simulation, the results of the other ten surrogatebased dynamic simulations were close to the EFbased simulation results. On average, the RMSE and MAE for 3D translations/rotations were less than 0.2 mm/0.6 deg and 0.3 mm/ 1.3 deg, respectively. For the medial contact forces/torques, RMSE and MAE were less than 20 N/0.540 Nm and 46 N/1.11 Nm, respectively. For the lateral contact forc es/torques, RMSE and PAGE 51 51 MAE were less than 20 N/0.561 Nm and 41 N/1.12 Nm, respectively. For both the surrogate contact model and the EF contact model, the Co Pbased approach for predicting wear volumes reproduced the elementbased EF results to within 2% error (Table 31). Discussion Computational speed is a major limiting factor for performing sensitivity and optimization studies of dynamic contact problems. This paper has presented a novel surrogate modeling approach for performing efficient dynamic contact simulations. The approach requires use of the contact surfaces and the knowledge of sensitive dire ctions to generate the sample points while no computationally expensive geometry evaluation is required during the dynamic simulation. The methodology was successfully applied to replace the EFbased artifici al knee contact model created from the manufactured CAD data. The evaluation demonstrated the computational efficiency of the approach by performing 5 Monte Carlo analyses and highlighted that the wear predicted form a simulator machine was sensitive to the changes in machine inputs. The results of 11 dynamic simulation comparisons demonstrated that the presented surrogate contact model accurately reproduced the EFbased dynamic simulati on results even when the input curves were varied from their standard values during the Monte Carlo analyses. The proposed surrogate contact model is base d on two key concepts. First, knowledge of sensitive directions is used to generate pre liminary sample points w ithout produci ng physically unrealistic outputs. This key concept is based on the assumption that th e value of each contact load is dominated by the value of the correspond ing pose parameter. This assumption implicitly assumed that the derivatives of contact for ces/torques with respect to the noncorresponding translations/rotations were not significant. To validate this assumption, the derivatives with respect to the noncorresponding pos e parameters were calculated for the proposed TKR design and the results proved that Fy and Tx were truly dominated by the values of Y and respectively. PAGE 52 52 The results were also consistent with previous study done on the same TKR design (Fregly et al., 2008) reported that the variations in superiorinferior translation and varusvalgus rotation have much more influence on the contact quantities than did variations in the remaining DOFs. Although this assumption is valid for this study, there are other cases in which it fails. For example, a hip joint (a sphere in a spherical cup) would have three sensitiv e translation directions and each sensitive direction might have less infl uence on the corresponding force than other two forces. Under this circumstance, we may no long er sample the corresponding force but instead sample one of the other two forces. Future data collected from a hip joint should be used to investigate this issue further. Second, a reasonable design space approach is developed to refine the design space and generate an accurate surrogate contact model. It screened out infeasible points from the design space and provided more accurate surrogate models. This study used the initial surrogate models developed from a subset of the samp le points to screen all sample poi nts before they were used in the computational experiments. These initial su rrogate models may not provide highly accurate predictions, but they efficiently identified most of the appropri ate sample points even based on their limited prediction abilities. About 88% (i .e., 500 out of 563 and 496 out of 566) of the accepted sample points turned out to be within the reasonable desi gn space. The additional static analyses performed on the rest of the 2000 sample points further indicated that only 5% of the bad points (including the unrealistic and infeasible points) were not captured by the initial surrogate models. A hybrid approach was utilized to develop the final surrogate contact model for couple reasons. Theoretically, two fitting approaches could be used to develop the final surrogate contact model: samplebased fitting and usebased fitting. For samplebased fitting, we fitted Fx, Fz, Ty, PAGE 53 53 Tz, and two sensitive directions (i.e., Y and ) as functions of sample inputs and used a rootfinding techni que to solve for Fy and Tx during the forward dynamic simulation. However, the extra computational time spe nding on the rootfinding made this approach less appealing in practice. For usebased fitting, we fitted all contact loads as func tions of six pose parameters and then performed the forward dynamic simulatio n by manipulating all six pose parameters simultaneously. Although the preliminary results showed that this approach successfully completed a dynamic simulation within 5 seconds of CPU time, the accuracy was not as good as the current approach. Thus, given the weakness of both fitting approaches, we developed a hybrid approach to combine these two approaches for performing fast and accurate dynamic simulation. Use of the surrogate contact model to replace repeated contact analyses is worthwhile for several reasons. The primary benefit is improved computational efficiency. Although the computational cost of performing repeated static analyses is inevita ble, this cost is pay up front cost since the sampling process is only performed once. In addition, such cost would be recovered quickly when performing sensitivity or optimization studies. In the presented application, the CPU time spent on 1060 static analyses was approximately 6 hours. However, the use of the surrogate contact model in place of the EF model reduced Monte Carlo analysis ti me from 284 hours to 1.4 hours, a nearly 99% reduction. While the applic ation has demonstrated that surrogate contact models can be beneficial for sens itivity studies, the greatest computational benefit is likely to occur for optimization studies. For stochastic se nsitivity studies, the mean value method (Wu et al., 1990) has been used to repr oduce the Monte Carlo wear predictions while requiring only 6% of computational time (Pal et al., 2008). In contrast, no good me thods exist for reducing the number of repeated simulations in optimization studies. For thos e situations, the factor of 200 reduction in computation time (i.e., 17 minutes to 5 seconds) could mean the difference between PAGE 54 54 an impossible and an achievable optimization stu dy. Other benefits of using surrogate contact models were provided in Chapter 2. The sensitivity analysis results demonstrate the importance of controlling machine inputs closely, especially the femoral flexion angle. It should be noted that th e current surrogatebased wear prediction did not include creep information. However, the absence of creep calculation should not affect the wear volume prediction significantly (Zhao et al., 2008). The Monte Carlo analyses revealed when small variations were imposed on the input motions and loads, only variations in the femoral flexion angle resulted in large changes in predic ted wear volume (Figure 35). This finding was consistent with results from our previous study where the similar Monte Carlo analyses were applied to a similar TKR design constrained to planar motion in the Stanmore machine (Lin et al., 2008). Pal et al. (2008) also discovered a similar result where they developed a probabilistic wear prediction model to identify that femo ral flexionextension alignment was one of the parameters most affecting predicted wear. The current modeling approach is limited in at least three as pects. First, the results of current central finite difference scheme were in fluenced by the coordinate system definition of the master body. For instance, X translation might become as sensitive as Y translation once we rotated our tibial insert coor dinate system 45 degree about it s zaxis. However, there is a redundancy between these two sensitive directions (i.e., X and Y translation) since they both are highly related to contact force Fy in original yaxis directi on. To minimize the redundancy between the sensitive directions, the technique likes principal com ponent analysis could be used to determine a principal sensitive direction from any highly correlated sensitive directions. More extensive data would be required to dr aw any conclusions about the effectiveness of principal component analysis for reducing the redundancy. PAGE 55 55 Second, new surrogate contact models must be generated if geometries or material properties of the contact bodies are changed. Ho wever, this limitation is not serious for the presented application since recent studies have reported that predicte d wear volume (but not wear area or depth) is relatively insensitive to whether or not the surface geometry is changed gradually (Knight et al., 2007; Zh ao et al., 2008). Furthermore, the additional design variables process (e.g., make material property an additional surrogate model input) could be added into the surrogate modeling process. Third, dynamic simulations obtained from the computationally expensive model must be performed to evaluate the accuracy of the surroga te contact model. This limitation prevents us from evaluating all surrogatebased dynamic simu lations performed in Monte Carlo analysis. In general, the predictive capability of a surr ogate model is evaluate d either by using the nontraining points (i .e., points not used to develop a surrogate model) or by performing the prediction error sum of squares analysis (M artens and Naes, 1989) on the training points. However, these two methods may be insufficien t for dynamic contact problems since it is difficult to draw firm conclusions about the co ntributions that the predictive capability of a surrogate model make to accura tely reproduce the original dynamic simulation. One of the possible alternatives to evaluate the surrogate based dynamic simulation without relying on the computationally expensive model is to perfor m a probability analysis. Such analysis would establish likelihood measures of a new pose parameter given by the forward dynamic simulation being consistent with the training points. Higher probability indicates greater confidence that the surrogate model will produce a more accurate pr ediction, and consequently improve the quality of dynamic simulation results. PAGE 56 56 In summary, this study has presented a deta iled surrogate contact modeling approach for significantly improving the computational speed of 3D dynamic contact simulations. The evaluation demonstrated the advantage of using th e surrogate contact model to perform sensitivity and optimization studies incorporating multiple de formable contacts. In the biomechanics field, such computationally efficient model would be id eal for the simultaneous predictions of contact forces and muscle forces since optimization requiring thousands of contact simulations is the essential part of the prediction process. Exte nsion of the current approach to a fullleg musculoskeletal model incorporating multiple su rrogate contact models is one of the ongoing research topics. Once surrogate contact models ar e developed for all major lower extremity joints, it should be possible to perform dynamic simu lations in a short amount of CPU time. Figure 31. A 6 degree of freedom total knee replacement contact model. The tibial insert (light blue object) poses three DOFs (a nteriorposterior translation x mediallateral translation z internalexternal rotation ) and the femoral component (dark blue object) poses another three DOFs (superiorinferio r translation y varusvalgus rotation and flexionextension ). PAGE 57 57 Table 31. Comparison between the wear volume predictions (mm3) made with the elastic foundation contact model and the surrogate contact model for dynamic simulations with different input curves. Quantity varied Sum of component weights Elastic foundation contact model Surrogate contact model Difference (%) Fx 28.44 28.51 0.24% Fy 30.31 30.27 0.14% Ty 18.14 18.24 0.56% Flex 17.10 17.23 0.72% All Maximum 29.76 29.52 0.80% Fx 29.49 29.20 1.02% Fy 30.19 29.73 1.53% Ty 27.63 27.37 0.96% Flex 27.08 26.90 0.66% All Minimum 30.11 30.08 0.10% PAGE 58 58 Figure 32. Overview of the surrogate modeling approach. A) N points are sampled in the ( X Z Fy, Tx) design space where the first M points are used to generate the initial surrogate models. B) Static an alysis is performed at the selected sample points using an elastic foundation (EF) c ontact model. C) The init ial surrogate models are developed based on the first M sample point s. D) The initial surrogate models are used to predict the static analyses results for N sample points. E) The predictions are screened based on the ranges defined by the EFbased dynamic simulation results. F) Additional static analyses are performed to obtain the outputs of interest for all remaining sample points. G) The final su rrogate models are developed based on the remaining sample points. H) Given the current values of six pose parameters, the surrogate contact model calculates the contact lo ads throughout the dynamic simulation. PAGE 59 59 Fx Fy Ty Flex All 10 15 20 25 30 35 Wear Volume (mm 3 )Quantity Varied 10th Percentile 25th Percentile 50thPercentile 75th Percentile 90th Percentile Outlier Figure 33. Box plot distribution of predicted wear volum e generated by four Monte Carlo analyses using the surrogate contact model. PAGE 60 60 4 2 8 X (mm) 42 44 46 Y ( mm ) 0 0.5 1 2 0 2 Z (mm) 10 0 10 (deg) 10 0 10 (deg) 0 0.5 1 60 30 0 (deg)Kriging Elastic Foundation Figure 34. Comparison between the six pose parameters made with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation. PAGE 61 61 200 0 200 Medial SideFx (N) 0 750 1500 Fy (N) 0 0.5 1 200 0 200 Fz (N) 200 0 200 Lateral Side 0 750 1500 0 0.5 1 200 0 200 Figure 35. Comparison between the mediallateral contact fo rces made with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation. PAGE 62 62 0 20 40 Medial SideTx (Nm) 5 0 5 Ty (Nm) 0 0.5 1 20 0 20 Tz (Nm) 40 20 0 Lateral Side 5 0 5 0 0.5 1 20 0 20 Figure 36. Comparison between the mediallateral contact to rques made with the elastic foundation contact model and the surrogate contact model for the nominal dynamic simulation. PAGE 63 63 CHAPTER 4 SIMULTANEOUS PREDICTION OF MUSCLE AND CONTAC T FORCES IN THE KNEE DURING GAIT Introduction Knowledge of muscle forces during walki ng is necessary to characterize muscle coordination and function, which ar e factors in neurological disorders such as cerebral palsy and stroke (Schmidt et al., 1999; Higginson et al., 2006; Gibson et al., 2007; Sujith, 2008), and to characterize joint and softtissue loading, which are factors in bone and joint diseases such as osteoarthritis and patellofemoral pain syndrome (H urwitz et al., 2000; Fregly et al., 2007; Milner et al., 2007). For all of these disorders, the abili ty to determine muscle and joint loads accurately in vivo could lead to significant improvements in clinical outcome. Unfortunately, direct measurement of muscle forces in vivo is possible only under limited special circumstances (Schuind et al., 1992; Dennerlein et al., 1998; Finni et al., 1998; Ishikawa et al., 2005), and similarly for joint contact loads (Stansfied et al., 2003; D'Lima et al., 2006). Due to these limitations, musculoskeletal computer models have become the primary scientific approach for predicting muscle forces and contact loads during human movement (Anderson and Pandy, 2001a; Buchanan et al., 2004; Higginson et al., 2006; Erdemir et al., 2007). A primary challenge to developi ng such predictions is the nonuniqueness of the calculated muscle forces, often referred to as the muscl e redundancy problem. Since more muscles act on the skeleton than the number of degrees of fr eedom in the skeleton, an infinite number of possible muscle force solutions ex ist. Consequently, optimization is commonly used to solve the indeterminate problem (Kaufman et al., 1991; Anderson and Pandy, 2001; Erdemir et al., 2007). This approach starts by calculati ng the net reaction loads (i.e., three forces and three torques) at each joint using a dynamic musculoskeletal model and movement data collected from a patient. A common constrained optimizatio n problem is then used to estimate muscle forces by PAGE 64 64 minimizing the sum of squares of muscle activati ons with constraints that muscles balance some (but not all) of the net reacti on loads. In physiological terms, the net loads are produced by a combination of muscle and contact forces. Thus, if a contact model is not included explicitly in the musculoskeletal model, then assumptions need to be made about which loads are generated primarily by muscles (i.e., little contribution from contact loads), and these loads are the only ones that can be used as constraints. In contrast if a contact model is included explicitly in the musculoskeletal model, then all six net loads can theoretically be used as constraints in the optimization problem formulation, thereby elimin ating the need for assumptions about the lack of contact contributions to particular loads. To date, no study has included a contact model explicitly in the musculoskeletal model to predict knee joint muscle forces and contac t loads simultaneously during the optimization process. Previous musculoskeletal models of th e knee generally used c onstraintbased joints (e.g., pin joint, ballandsocket joint) with limite d degrees of freedom (DOFs) when in reality the knee joint possesses 12 DOFs (6 for the tibiofemoral joint and 6 for the patellofemoral joint). Despite this limitation, some studies have reported that increasing the DOFs in a constraintbased knee model from 1 to 3 significantly affects the predicted muscle forces (Buchanan et al., 1996; Glitsch et al., 1997; Li et al., 1998). Though this issue can be eliminated by replacing each constraintbased joint with a geometrybased joint (i.e., a contact model), the computational cost of solving repeated contact proble ms currently makes such an appr oach infeasible (Fregly et al., 2005). This chapter compares simultaneous muscle and contact force pred ictions in the knee generated using different optimization problem fo rmulations with and without matching of in vivo medial and lateral contact fo rce measurements as constraints. The goals of this chapter are PAGE 65 65 twofold. The first is to evaluate the ability of a musculoskeletal model incorporating surrogate contact models to reproduce optimization predic tions of muscle and co ntact forces generated using the same musculoskeletal model but incorporating elastic foundatio n contact models. The second is to utilize the musculos keletal model with surrogate cont act to evaluate the ability of eight optimization problem formul ations to reproduce in vivo contact force measurements provided by an instrumented knee implant. A 12 DOF surrogate contact model of the tibiofemoral (TF) and patellofemoral (PF) jo ints is created and incorporated into the musculoskeletal model for performing the optim izations. The predicted muscle forces are evaluated using in vivo medial and lateral contac t force measurements from a single patient with an instrumented knee implant (D'Lima et al., 2006; D'Lima et al., 2007; Zhao et al., 2007b). Methods Data Collection Experimental gait data were collected from a single patient (male, age 80 years; mass 68 kg) implanted with an instrumented total knee replacement (TKR) in his right leg. Institutional review board approval and patie nt informed consent were obtained. The instrumented knee implant consisted of four uniaxial force transducers, a micro transmitter, and an antenna ((D'Lima et al., 2005). In vivo tibial force data we re recorded simultaneously with video motion (Motion Analysis Corporation, Santa Rosa, CA) and ground reaction data when the patient performed different walking trials. For each trial, the distribution of axial contact forces between the medial and lateral compartmen ts of the knee was calculated fr om the four force transducer measurements using validated regression equatio ns (Zhao et al., 2007a). A modified Cleveland Clinic marker set with additional markers placed on the feet was used for the video motion data collection. Static markers placed over the medi al and lateral epicondyles of the knee and the PAGE 66 66 medial and lateral malleoli of the ankle were us ed along with dynamic markers to define segment coordinate systems and provide threedimensional (3D) movement data. Inverse Dynamics Loads Net forces and torques (henceforth referred to as inverse dynamics loads) at the knee joint were calculated using a publis hed bottomup inverse dynamics approach (Reinbolt et al., 2008). In brief, a foursegment leg model was deve loped specifically for th e patient based on the collected surface marker data. The four segments c onsisted of the foot, shank, thigh, and pelvis. The foot was treated as the root segment connected to the laboratoryfixed coordinate system via a six DOF joint. Two nonintersecting pin joints were used to model the ankle (van den Bogert et al., 1994), a single pin joint was used to model the knee (Reinbolt et al., 2005),and a ballandsocket joint was used to model the hip. A series of optimizations were applied to the model to calibrate its inertia parameter values (body segment ma sses, mass centers, and central principal moments of inertia) and joint parameter values (positions and orientations of ankle, knee, and hip joint axes in the body segments ) by minimizing the errors between model and experimental marker lo cations (Reinbolt et al., 2005). The optimized tibiofemoral (TF) joint center was defined as the point on the optimized knee axis at the midpoint between the medial and lateral femoral epicondyle markers. Six inverse dynamics loads at the knee for one cycle (right heel strike to subsequent ri ght heel strike) of a selected gait trial were then calculated about the optimized TF joint center a nd expressed in the sh ank reference frame. Corresponding medial, lateral, and total tibial contact forces for the same gait trial were provi ded by the instrumented knee implant. Musculoskeletal Model A composite implant/bone/muscle knee model was constructed based on medical imaging data obtained from two subjects. The model cons isted of computeraided design (CAD) models PAGE 67 67 of the knee implant components, geometric models of four bones (femur, patella, tibia, and fibula), and origin and insertion locations for eleven muscles and the patellar ligament (Figure 41). Coarse bone models were constructed fo r the patient using his postsurgery computed tomography (CT) data. Both the bones and the metallic implant components were segmented from the CT images, allowing the implant CAD models to be aligned to the segmented outlines of the coarse bone models. Fine bone models with corresponding muscle and ligament attachment points were constructed for a health y subject of similar stature using magnetic resonance imaging (MRI) data. Th e healthy subjects fine bone models were aligned to the segmented outlines the patients coarse bone mode ls, resulting in a composite geometric model consisting of the implant CAD models and the fine bone models with corresponding muscle and ligament attachment points. All alignment tasks were performed to minimize distances between aligned objects using the commercial reverse e ngineering software Ge omagic Studio (Raindrop Geomagic, Research Triangle Park, NC). Once the registration process was finished, the kinematic structure of the inverse dynamics leg model was surperimposed onto the composite geometric model by aligning the ankle, kn ee, and hip centers of the two models. Eleven major muscles spanning the knee were included in the mode l: Vastus medialis (VM), vastus lateralis (VL), vastus intermedius (VI), rectus femoris (RF), semimembranosus (SM), semitendinosus (ST), biceps femoris long head (BFLH), bicpes femoris short head (BFSH), tensor fascia latae (TFL), gastrocnemiu s medial head (GM), and gastrocnemius lateral head (GL). Each muscle was activ ated independently except for medial hamstrings (group of SM and ST) and vasti (group of VM, VI, and VL), which were each controlled by a single muscle activation, with muscle activation being used as an indicato r of neural control input. Muscle force was generated by multiplying each muscle activation by its corresponding peak isometric PAGE 68 68 force, as a previous study has shown that this simple muscle model produces nearly identical results to those of more complex muscle mode ls when gait is being simulated (Anderson et al., 2001b). The patellar ligament was m odeled as three parallel springs with a total stiffness of 2000 N/mm (Reeves et al., 2003). The TF and patellofemoral (PF) joints were each modeled as 6 DOF joints with elastic contact (Figure 42). The femoral component was fixed to gr ound and the tibial insert was connected to the femoral component via a 6 DOF joint (generalized coordinates TTTTTTX,Y,Z,,,and for the three translations and rotations, respectively). The patella was also connected to the femoral component vi a a 6 DOF joint (gener alized coordinates PPPPPPX,Y,Z,,,and ). These two 6 DOF joints were used to quantify the relative position and orientation of each pair of articulating im plant components for contact calculations. For the TF joint, anteriorposterior translation, intern alexternal rotation, and flexionextension were prescribed to match the kinematics measured from video motion, since c ontact force calculations are insensitive to small errors in these DOFs (Fre gly et al., 2008). In contrast, for the PF joint, none of the 6 generalized coordinates was pres cribed. No relative motion was allowed between the thigh and femoral component or between the shank and tibial insert. For each time frame, the values of the 9 free generalized coordinates were determined based on sta tic equilibrium of the shank and patella produced by the applied inverse dynamics loads, contact loads, muscle forces, and patellar ligament force. Surrogate Contact Model Two surrogate contact models, one for the TF joint and one for the PF joint, were developed to perform repeated contact analyses efficiently. For each join t, a surrogate contact model was created using the corresponding im plant CAD geometry, a validated elastic PAGE 69 69 foundation (EF) contact model (Fre gly et al., 2003; Bei et al., 2004; Zhao et al., 2008), and the methods presented in Chapter 3. In brief, we used a threedimensional EF contact model with nonlinear material properties (Fre gly et al., 2007) to generate sa mple points for creating TF and PF surrogate contact models. For both contact models, superiorinfer ior translation and varusvalgus rotation were identifi ed as sensitive directions sin ce the contact loads were highly sensitive to small variations in these generalized coordinate values. We therefore sampled force and torque (Fy and Tx) in these two sensitive directions and pose (X, Z, and ) in the remaining four insensitive directions. All sampling was again performed using a Hammersley quasirandom sequence. Physically realistic upper and lower b ounds on the sample point inputs (i.e., X, Z, Fy and Tx) were determined by performing a sequence of static analyses across the gait cycle using a knee model incorporating EF contact mode ls. A simplified version of model was created for this purpose where vasti was the only muscle and possessed a constant force of 200 N. For each time frame, the pose of the TF joint was prescribed to reproduce the in vivo medial and lateral tibial contact force measurements while the pose of the PF joint was determined based on static equilibrium of the patella. A reasonable design space approach as descri bed in Chapter 3 was used to identify a feasible set of sample points within the complete set of sample points. Sample points where the corresponding contact m odel outputs were outside thei r allowable bounds were deemed infeasible. First, initial surrogate contact mode ls developed from a subset of the sample points were used to screen the complete set of sample points. Second, repeated static analyses were performed with the EF contact model to find physically realistic poses for the sample points passing the screening process. For each static an alysis, the static contact loads in sensitive directions were compared to their corresponding applied values. A sample point was chosen only PAGE 70 70 when the maximum percent error was below 0.1%. Once all static analyses were completed and screened, contact loads corresponding to sensitive directions were fitted as a function of the six pose parameters, while contact lo ads corresponding to in sensitive directions were fitted as a function of the contact loads in sensitive directions and pose pa rameters in insensitive ones. Muscle and Contact Force Prediction Once the surrogate contact models were inco rporated into the knee model (Figure 43), muscle forces and associated TF and PF contact loads were predicted simultaneously using a twolevel optimization a pproach. An outer level optimization (i.e., muscle force optimization) modified the muscle activations to minimize the cost function calculated by the inner level optimization (i.e., pose optimization) over all ev aluated time frames. Given the initial guess for muscle activations, the pose optimization modified the nine free pose parameters (q) to to determine a static configurat ion that minimized the corres ponding nine residual loads: 3 222 22 1minResidual Residual Residual ResidualResidual tibiaSI tibiaML tibiaVV patella patella i iFFTFT q (41) whereResidual tibiaSIF, Residual tibiaMLF, and Residual tibiaVVTare the three tibia l residual loads corresponding to superiorinferior translation, medi allateral translation, and varu svalgus rotation, respectively, while Residual patellaFand Residual patellaT are patellar residual loads corresponding to three mutually perpendicular directions (i=1, 2, and 3). The tibial residual loads quantify the imbalance in inverse dynamics loads, muscle forces, and cont act loads acting on the tibia, while the patellar residual loads quantify the imbalance in muscle forces and contact load s acting on the patella. The muscle force optimization modified the muscle activations to minimize one of two cost functions: 1) the sum of squares of muscle activ ations (Crowninshield and Brand, 1981; Forster et al., 1998; Anderson and Pandy, 2001) or 2) the sum of absolute values of the three PAGE 71 71 compressive loads (medial TF, lateral TF, a nd PF) (Schultz and A nderson, 1981), as shown below: 11 2 1Minimize cost function 1:min or cost function 2:min Subjectto constraints 1:0 0 or constraints 2:0i i medcontactlatcontactcontact tibiaSItibiaSIpatellaSI i residual tibiaFE i residua tibiaFEa FFF a T a T a a 0 0 0 by varying muscle activation design variables1,...,11l residual tibiaAP residual tibiaIE iF T ai (42) Each cost function was employed wi th two different constraint sets Both constraint sets limited muscle activations to be positive and required the flexionextension torque applied to the tibia by muscle and contact loads to match the flexionex tension torque from inverse dynamics. The first constraint set represents the approach commonly used in the lite rature. The second constraint set added matching of the anteriorposterior for ce and internalexternal torque from inverse dynamics, thereby narrowing the feasible solution space further. Muscle activations were not constrained to be less than 1 (i .e., the physiological s ituation) to allow the optimizer to find feasible solutions that we re not limited by upper bounds on muscle force production. Eight optimization problems were formulat ed based on Equation 42 (Table 41). The problems consisted of all possible combinations of the two cost f unctions and two constraint sets, with each problem being solved with and without additional constraints (not shown in Equation PAGE 72 72 42) to match the in vivo medial and lateral contact force measurements. Problems that forced the solution to match the in vivo contact fo rce measurements were termed "cheating" formulations (i.e., problems C1 through C4 in Table 41) whil e those that did not force the solution to match the in vivo measurements were termed "honest" formulations (i.e., problems H1 through H4 in Table 41). Problem C1 was used to evaluate the ability of a musculoskeletal model with surrogate contact to reproduce optimization predictions of muscle and contact forces generated using the same musculos keletal model but with elastic foundation contact. Differences between the solutions generated by the two mo dels were quantified by calculating maximum absolute errors (MAE) and root mean square er rors (RMSE) in bone poses, contact loads, and muscle forces. The "cheating" formulations al so allowed us to verify that the proposed musculoskeletal model was theoretically capab le of matching the in vivo contact force measurements and to generate muscle force and PF contact force solutions that are likely to be within a physiological range. Th e "honest" formulations allowed us to evaluate how well different optimization problem formulations could predict the in vivo contact force measurements without knowing them. To k eep computation time tractable, optimization solutions were generated at 5% increments across the gait cycle. Results The optimization performed with the surrogate contact models closely repr oduced the 3D motions, contact forces, and muscle forces pred icted with the EF contact models for problem formulation C1. On average, for both TF and PF jo ints, the RMSE for translations/rotations and forces/torques were less than 0.6 mm/ 0.8deg an d 7 N/ 0.5 Nm, respectively (Table 42 and Table 43). The RMSE between surrogatebased and EFbased muscle force predictions were below 15 N (Table 44). While each optimization pe rformed with the EF c ontact model required PAGE 73 73 approximately 90 minutes of CPU time, each optim ization performed using the surrogate contact model required 2 minutes. The medial and lateral tibial contact forces obtained from all cheati ng formulations but not all honest formulations were in good ag reements with in vivo contact measurements (Figure 44). Overall, honest formulations predicted medial tibi al contact forces that were comparable in magnitude and lateral tibial cont act forces that were too low compared to measured tibial contact forces. Among other hon est problem formulations, formulation H1 was able to obtain a better agreement between predicted and measured medial contact force throughout whole gait cycle while formulation H3 achieved the best agreement for the first 25% of the gait. For PF joint, the contact forces Fy obtained form eight problem formulations had the similar shape. Muscle forces calculated from eight problem formulations were significantly different from each other (Figure 45). Compared to honest formulations, cheating formulations generated higher central and medial muscle forc es and higher lateral mu scle forces over the midstance and most of stance phase, respectively. None of the muscles exceed their maximum isometric forces except for BFSH and TFL. Discussion This study used a novel surrogate contact mode ling approach to faci litate simultaneous prediction of muscle and contact forces in the knee during gait. The surrogate modeling approach involves replacing computationally costly contact models with computationally cheap contact models constructed using data points sampled from the original contact models. Once the surrogate models are constructed, they are used in place of the original models to eliminate a computational bottleneck caused by repeated cont act analyses. Though the computational cost of constructing the three surrogate contact models is high (on the order of 30 hours of CPU time), PAGE 74 74 this time is quickly redeemed by the greatly re duced computational cost of each one cycle gait optimization (about 42 minutes of CPU time co mpared to 32 hours of CPU time for the EF contact models), especially when multiple gait optimizations are performed. The cheating problem formulations demonstrated that the proposed musculoskele tal model provided enough dimensionality in the magnitudes and directions of muscle forces to reproduce the in vivo contact force measurements. In addition, problem C1 de monstrated that the surrogate contact models could accurately reproduce the EFbased optimiza tion results in a much shorter amount of CPU time. The "honest" formulations evaluated four different optimization problems and demonstrated that contact loads and muscle forces were significan tly influenced by the optimization formulations, with different proble m formulations matching medial contact force the best at different points in the gait cycle. Simultaneous prediction of muscle and contact forces provides at least three advantages over sequential prediction, where muscle forces calculated without cons ideration of contact forces are subsequently used to calculate cont act forces (Taylor et al., 2004; Kim et al., 2008). First, a simultaneous approach eliminates the need to make assumptions about the inverse dynamics loads to which contact forces do not co ntribute. By incorporat ing contact geometry into the model, one can satisfy all six inverse dynamics loads simultaneously at the tibiofemoral joint, and more loads than just the flexionexte nsion torque can be used as constraints. This benefit was demonstrated by problem formulati ons C2, C4, H2, and H4, where three inverse dynamics loads were used as equality constraint s. The results from th ese four problems were different from the corresponding problems where onl y the flexionextension torque served as a constraint, consistent with prev ious studies reporting that musc le force predictions are greatly PAGE 75 75 influenced by the number of inverse dynamics lo ads used as constraints (Buchanan et al., 1996; Glitsch et al., 1997; Li et al., 1998). Second, a simultaneous approach makes it po ssible to utilize cost functions that are unrelated to muscles. This benefit was demons trated by problems C3, C4, H3, and H4, where the cost function minimized compressive force in th e three compartments of the knee rather than minimizing muscle activations, forces, stresses, or other musclerelated quantities. Minimization of compressive contact force predicted better me dial contact force results in early stance phase than did minimization of muscle activations. This observation suggest s that the body may use different cost functions (or differe nt combinations of cost functions ) at different points in the gait cycle to produce muscle and contact forces (Figure 44). Third, a simultaneous approach allows one to match in vivo contact force measurements as additional constraints. Without these constraints, it would have been extremely difficult to evaluate whether our musculoske letal model was theore tically capable of re producing the in vivo medial and lateral contact force measurements. The addition of these co nstraints also made it possible to hypothesize whet her the cheating solutions are represen tative of the in vivo situation. Whereas the predicted compressive forces on the patella were similar for all four cheating solutions (Figure 44), the predicted muscle forces were quite di fferent (Figure 45). Thus, the predicted PF contact forces may provide reasonable estimates for the in vivo situation, while it is unclear which cheating formulation (if any) pr oduced the most physiologi cally realistic muscle force predictions. This finding is not surprising gi ven the amount of indeterminacy in the muscle force prediction problem. To gain confidence in the muscle force predictions, one would need to identify an honest problem formulation that consistently predicted the correct medial and lateral contact forces for a vari ety of gait motions. Muscle electromyographic data from the same PAGE 76 76 patient would also be valuable for evaluati ng the reasonableness of the predicted muscle activations. For the PF contact forces the fact that the shapes (but not magnitudes) of the curves are similar to published studies th at did not utilize contact models provides additional evaluation of these results (Ward and Power, 2004). Our results highlight at least four possible musculoskeletal model inadequacies that may have contributed to the high errors observed in the predicted lateral contact forces. First, control of frontal plane stability may be critical for pr edicting lateral contact force accurately. Honest formulations may have generated low lateral co ntact force since the need for frontal plane stability was not accounted for in the problem formulation. It may be possible to increase frontal plane stability, and thus lateral contact force, by preferentially selecting muscles that contribute the most to varusvalgus stiffness (Nathan et al., 2008). Second, our mu sculoskeletal model only included a knee joint model, though some muscles were biarticular (e.g., TFL). Consequently, forceproducing constraints impos ed on biarticular muscles by neighboring joints were not account for in our problem formulations. In particular, for lateral biarticul ar muscles crossing the knee and hip (e.g., TFL), the need to balance the frontal plane moment at the hip to keep the pelvis level could contribute to increased lateral muscle force and hence lateral knee contact force. Third, the lack of lateral collateral ligam ent (LCL) in our musculoskeletal model may have affected the lateral contact for ce predictions, as the LCL is the primary structure resisting varus angulation of the knee (Gollehon et al., 1987; LaPrade et al., 2004). Cl early, adding an LCL with high stiffness to the model could increase lateral compressive load. To reproduce the in vivo lateral contact force measurements, the LCL wo uld need to provide approximately 400 N of force throughout most of the gait cycle (F igure 413). However, this requirement is physiologically unrealistic as the mean failure strength for the LCL is approximately 450 N PAGE 77 77 (Ciccone et al., 2005). Furthermore, in vitro measuremen ts indicate that la teral compartment contact force due to ligaments alone is on th e order of only 50 to 100 N (Hawary et al., 2003). Fourth, no muscle physiology (i.e., forcelength and forcevelocity properties) and no activation or contraction dynamics were in corporated into our muscle m odel. Anderson et al. (2001b) reported that inclusion of those properties had only a small effect on muscle force predictions during gait, suggesting that the accu racy of the muscle forces pr edicted by optimization methods depends mainly on the accuracy of the inverse dynamics loads. In summary, this study presented a way to predict muscle and contact forces simultaneously at the knee us ing surrogate contact mode ling methods. The cheating optimization problem formulations revealed that our musculoskeletal model with surrogate contact was capable of reproduci ng in vivo medial and lateral c ontact force measurements and that the surrogate contact solutions were reliable. The honest formulations revealed that our musculoskeletal model could predict in vivo medial contact force measurements well but underestimated lateral contact force measurements. Demonstra ting the feasibility of this simultaneous approach is a crucial step toward its application to a large number of trials or patients. Future studies will increase musculoskeletal model complexity by including control of frontal plane stability in the optimization form ulation and by adding ligaments and other lower extremity joints. PAGE 78 78 Figure 41. Overview of the pro cess of using musculoskeletal mode l to predict muscle forces and contact loads simultaneously over the entire gait cycle. Surrogate contact model is created to provide computa tionally efficient contact an alysis for musculoskeletal model during optimization process. The optimization is solv ed subject to constraints that the inverse dynamics loads were bala nced out by the combination of predicted contact and muscle forces. The predictions are evaluated by comparing measured and predicted tibial contact forces. Table 41. Eight different optimization problem fo rmulations. C1 to C4 represent four cheating formulations while H1 to H4 repr esent four honest formulations. Constraints 1 Constraints 2 Cost Function 1 C1 (H1) C2 (H2) Cost Function 2 C3 (H3) C4 (H4) PAGE 79 79 TXTYTZPZPXPYP P P T T T Figure 42. A 12 degree of freedom total kne e replacement contact model. The femoral component is fixed to the ground while tibial insert and patella are connected to the femoral component via two 6 DOF joints (anteriorposterior translation X, superiorinferior translation Y, mediallat eral translation Z, varusvalgus rotation internalexternal rotation and flexionextension ). Table 42. Summary of root mean square a nd maximum absolute joint kinematics errors produced by surrogate contact model for problem formulation C1. TF joint PF joint RMSE MAE RMSE MAE X (mm) 0.001 0.002 0.3000.882 Y (mm) 0.005 0.020 0.5681.183 Z (mm) 0.017 0.042 0.4520.785 (deg) 0.007 0.024 0.7922.113 (deg) 0.001 0.003 0.1330.356 (deg) 0.000 0.001 0.3300.765 PAGE 80 80 Figure 43. Simultaneous approache for predictin g muscle forces and contact loads. The pose optimization performs repeated static anal yses to equilibrate the free DOFs of the musculoskeletal model given the current muscle activations. The muscle force optimization adjusts the muscle activati ons to evaluate the cost function and constraints calculated by the pose optimization. Table 43. Summary of root mean square and maximum absolute joint dynamics errors produced by surrogate contact model fo r problem formulation C1. TF joint PF joint Medial contact Lateral contact RMSE MAE RMSE MAE RMSE MAE Fx (N) 0.75 2.63 2.32 7.83 0.82 2.58 Fy (N) 0.49 1.88 1.91 8.54 6.42 16.50 Fz (N) 3.42 12.69 2.73 5.90 2.71 6.68 Tx (Nm) 0.14 0.40 0.35 1.30 0.05 0.15 Ty (Nm) 0.02 0.04 0.05 0.17 0.00 0.01 Tz (Nm) 0.03 0.06 0.08 0.25 0.04 0.11 PAGE 81 81 Figure 44. Comparison between the in vivo contact measurements (color in grey) and predictions from eight optimi zation problem formulations. Table 44. Summary of root mean square and maximum absolute muscle force errors produced by surrogate contact model fo r problem formulation C1. Peak Isometric Force (N) RMSE (N) MAE (N) RF 1319.8411.5346.92 VAS 6857.8913.5053.85 MHAMS 1941.467.0719.45 MGAS 825.477.6120.27 LGAS 825.4710.9425.16 BFLH 872.2512.1433.71 BFSH 228.267.9726.90 TFL 262.167.5418.56 PAGE 82 82 Figure 45. Predictions of mu scle forces from eight optim ization problem formulations. PAGE 83 83 CHAPTER 5 CONCLUSION The m ain motivation of bringing surrogate m odeling approach into biomechanics field is to efficiently predict the quantities that are difficult to measure in vivo (e.g., joint contact force, contact pressure, muscle force). Those quantities would be valuable for preventing and treating joint injuries and for improving the longevity of joint replacements. Although a variety of computational models have been used for this pu rpose, none of them is computationally efficient enough for studies involving repe ated contact analysis. Unlike other computational models, a surroga te modeling approach provides a way to eliminate contactrelated computational bottleneck s for any type of simulations possessing one or more articular contacts. Our results first de monstrated that the surrogate modeling approach can accurately and efficiently reproduce either 2D or 3D dynamic simulations results originally generated by the elastic foundation contact mode l. The fact that a Stanmore knee simulator machine is highly sensitive to small variations in motion and load inpu ts is presented as a practical application of usi ng surrogate contact model. Given what we have learned from chapter 3, sa me technique could be applied to the similar analysis in chapter 2. Instead of using Hertzi an contact theory and lin ear relationship between Fx and Fy, the concept of sensitive directions and reasonable design space approach could be used to generate a new accurate surrogate contact mode l for 2D dynamic contact simulation. The new surrogate should perform better than the original one since new approach will minimize the number of surrogates needed to calculate the contact forces. Followed by the successful single joint simulation, we then demonstrated that same surrogate modeling approach can be applied to a multijoint contact analysis for predicting muscle and contact forces simultaneously. The ta sk was done by integrating a surrogate contact PAGE 84 84 model into a musculoskeletal model of the knee. The result s showed that the proposed surrogatebased simultaneous approach could finish an optimization analysis within a practical CPU time to investigate the relationship between contact loads and musc le forces during gait. Although we have shown that the simultaneous a pproach could reproduce the in vivo contact measurements when they were used as additiona l constraints, further re search is required to resolve the underestimation of lateral contact force. Given the agreements of the surrogatebased results and EFbased results, a possible extension of current study could be to develop multiple surrogate contact models for all lower extremity joints other than just knee joint. Once finished, it would be possible to create a fullleg musculoskeletal model with multiple surrogate contact models that can be simulated in a short amount of CPU time. Another futu re topic involves the developm ent of the surrogate modeling software. Once developed, the software along w ith a brief tutorial will be freely available throughout the internet. These reso urces will allow researchers to learn how the surrogatebased modeling approach can be used for their projects. PAGE 85 85 APPENDIX A COMPUTATIONAL WEAR PREDICTION The wear an alysis calculated the depth of material removed from each element over N cycles based on Archards wear law (Archard and Hirst, 1956) 11 nn iiii iiNkpdNkptv (A1) where k is the assumed material wear rate (1 107 mm3/Nm; Fisher et al., 1994), i is a discrete time instant during the onecycl e dynamic simulation (1 through n), pi is the contact pressure on the element at that instant, and di is the sliding distance experien ced by the element, calculated as the product of magnitude vi of the elements relative slip velocity and time increment t. The wear volume was calculated by multiplying each element wear depth by the corresponding element A 1 n ii iVANkptA v (A2) Since i p A is just the normal contact force Fi applied to the element, the wear volume can be expressed concisely as 1 n ii iVNkFtv (A3) PAGE 86 86 LIST OF REFERENCES An, K.N., Him enso, S., Tsumura, H., Kawai, T., Chao, E.Y.S., 1990. Pressure distribution on articular surfaces: applic ation to joint stability analysis. Journal of Biomechanics 23, 10131020. Anderson, F.C., Pandy, M.G., 2001a. Dynamic op timization of human walking. Journal of Biomechanical EngineeringTransact ions of the Asme 123, 381390. Anderson, F.C., Pandy, M.G., 2001b. Static a nd dynamic optimization solutions for gait are practically equivalent. Journa l of Biomechanics 34, 153161. Archard, J.F., Hirst, W., 1956. The wear of me tals under unlubricated conditions. Proceedings of the Royal Society of London Se ries aMathematical and P hysical Sciences 236, 397&. Bai, B., Baez, J., Testa, N.N., Kummer, F.J., 2000. Effect of posteri or cut angle on tibial component loading. Journal of Arthroplasty 15, 916920. Balabanov, V.O., Giunta, A.A., Golovidov, O., Grossman, B., Mason, W.H., Watson, L.T., Haftka, R.T., 1999. Reasonable design space approach to response surface approximation. Journal of Aircraft 36, 308315. Barnett, P.I., Fisher, J., Auger, D.D., Stone, M.H., Ingham, E., 2001. Comparison of wear in a total knee replacement under different kine matic conditions. Jour nal of Materials ScienceMaterials in Medicine 12, 10391042. Bei, Y.H., Fregly, B.J., 2004. Multibody dynamic simulation of knee contact mechanics. Medical Engineering & Physics 26, 777789. Bitsakos, C., Kerner, J., Fisher, I., Amis, A.A ., 2005. The effect of muscle loading on the simulation of bone remodeling in the proximal femur. Journal of Biomechanics 38, 133139. Blankevoort, L., Kuiper, J.H., Huiskes, R., Grootenboer, H.J., 1991. Articular contact in a threedimensional model of the knee. Journal of Biomechanics 24, 10191031. Blunn, G.W., Walker, P.S., Joshi, A., Hardinge K., 1991. The dominance of cyclic sliding in producing wear in total knee replacements. Clin ical orthopaedics and related research, 253260. Bottasso, C. L., Prilutsky, B. I., Croce, A., Imberti, E., Sartirana, S., 2006. A numerical procedure for inferring from expe rimental data the optimization cost functions using a multibody model of the neuromusculoskeletal sy stem. Multibody System Dynamics 16, 123. Bouzid, A.H., de Technologie Superieure, E., Ch ampliaud, H., 1998. On the use of dual kriging interpolation for the evaluation of the gasket stress distribution in bolted joints. In Proceeding of the ASME Pressure Vessels and Piping Conference, 457, 7783. Box, G.E.P., Draper, N.R., 1987. Empirical mo delbuilding and response surfaces. Wiley, New York. PAGE 87 87 Brogliato, B., ten Dam, A.A., Paoli, L., Genot F., Abadie, M., 2002. Numerical simulation of finite dimensional multibody nonsmooth mechanical systems. Applied Mechanics Reviews 55, 107. Buchanan, T.S., Shreeve, D.A., 1996. An ev aluation of optimization techniques for the prediction of muscle activati on patterns during isometric task s. Journal of Biomechanical Engineering 118, 565574. Buchanan, T.S., Lloyd, D.G., Manal, K., Besier T.F., 2004. Neuromuscu loskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. Journal of Applied Biomechanics 20, 36795. Caruntu, D.I., Hefzy, M.S., 2004. 3D anatomically based dynamic modeling of the human knee to include tibiofemoral and patellofemoral jo ints. Journal of Biomechanical Engineering 126, 4453. Challis, J.H., Kerwin, D.G., 1993. An analytical examination of muscle force estimations using optimization techniques. Proceedings of the Inst itution of Mechanical Engineers Part H 207, 139148. Chang, P.B., Williams, B.J., Santner, T.J., Notz W.I., Bartel, D.L., 1999. Robust optimization of total joint replacements incor porating environmental variable s. Journal of Biomechanical Engineering 121, 304310. Chao, E.Y.S., 2003. Graphicbase d musculoskeletal model for biomechanical analyses and animation. Medical Engin eering & Physics 25, 201212. Chen, S., Chng, E.S., Alkadhimi, K., 1996. Regul arized orthogonal leas t squares algorithm for constructing radial basis func tion networks. International Jo urnal of Control 64, 829837. Cheng, B., Titterington, D.M., 1994. Neural Networ ks a Review from a Statistical Perspective. Statistical Science 9, 230. Ciccone, W.J., Bratton, D.R., Weinstein, D.M ., Walden, D.L., Elias, J.J., 2005. Structural properties of lateral collateral ligament reconstruction at the fi bular head. The American Journal of Sports Medicine 34, 2428. Cohen, Z.A., Roglic, H., Grelsamer, R.P., Henr y, J.H., Levine, W.N., Mow, V.C., Ateshian, G.A., 2001. Patellofemoral stress es during open and closed chain exercises an analysis using computer simulation. The American Jour nal of Sports Medicine 29, 480487. Cox, S.E., Haftka, R.T., Baker, C.A., Gro ssman, B., Mason, W.H., Watson, L.T., 2001. A comparison of global optimization methods for the design of a highspeed ci vil transport. Journal of Global Optimization 21, 415433. Cressie, N., 1992. STATIS TICS FOR SPATIAL DATA. Terra Nova 4, 613617. PAGE 88 88 Crowninshield, R.D., Brand, R. A., 1981. A physiologically base d criterion of muscle force prediction in locomotion. Journa l of Biomechanics 14, 793801. Daffertshofer, A., Lamoth, C.J.C., Meijer, O.G ., Beek, P.J., 2004. PCA in studying coordination and variability: a tutorial. C linical Biomechanics 19, 415428. Delp, S.L., Arnold, A.S., Speers, R.S., Moore, C.A., 1996. Hamstrings a nd psoas lengths during normal and crouch gait: implications for muscletendon surgery. Journal of Orthopaedic Research 14, 144151. Dennerlein, J.T., Diao, E., Mote, C.D., Rempel, D.M., 1998. Tensions of the flexor digitorum superficialis are higher than a current model predicts. Journal of Bi omechanics 31, 295301. DesJardins, J.D., Walker, P.S., Haider, H., Perry J., 2000. The use of a forcecontrolled dynamic knee simulator to quantify the mechanical perfor mance of total knee replacement designs during functional activity. Journal of Biomechanics 33, 12311242. Diwekar, U.M., 2003. Introduction to Applied Optimization. Springer. D'Lima, D.D., Hermida, J.C., Chen, P.C., Colwell, C.W., 2001. Polyethylene wear and variations in knee kinematics. Clinical Orthop aedics and Related Research, 124130. D'Lima, D.D., Townsend, C.P., Arms, S.W., Morri s, B.A., Colwell, C.W., 2005. An implantable telemetry device to measure intr aarticular tibial forces. Journa l of Biomechanics 38, 299304. D'Lima, D.D., Patil, S., Steklov, N., Slamin, J.E., Colwell, C.W., 2006. Tibial forces measured in vivo after total kne e arthroplasty. Journal of Arthroplasty 21, 255262. D'Lima, D.D., Patil, S., Steklov, N., Chien, S ., Colwell, C.W., 2007. In vivo knee moments and shear after total knee arthroplasty. Jo urnal of Biomechanics 40, S11S17. Elias, J.J., Wilson, D.R., Adamson, R., Cosgarea, A.J., 2004. Evaluati on of a computational model used to predict the patellofemoral contac t pressure distribution. Journal of Biomechanics 37, 295302. Erdemir, A., McLean, S., Herzog, W., van den Bogert, A.J., 2007. Modelbased estimation of muscle forces exerted during movements. Clin Biomech (Bristol, Avon) 22, 13154. Fang, H., RaisRohani, M., Liu, Z., Horstemeyer, M.F., 2005. A comparative study of metamodeling methods for multiobjective cras hworthiness optimization. Computers & Structures 83, 21212136. Finni, T., Komi, P.V., Lukkariniemi, J., 1998. Achilles tendon loading during walking: application of a novel optic fibe r technique. European Journal of Applied Physiology 77, 28991. Fishman, G.S., 1996. Monte Carlo: Concepts Algorithms, and Applications. Springer. PAGE 89 89 Forster, E., Simon, U., Augat, P., and Claes, L., 2004. Extension of a stateoftheart optimization criterion to predict cocontra ction. Journal of Bi omechanics 37, 577581. Fregly, B.J., Bei, Y.H., Sylvester, M.E., 2003. Experimental evaluation of an elastic foundation model to predict contact pressures in knee replacements. Journal of Biomechanics 36, 16591668. Fregly, B.J., Sawyer, W.G., Harman, M.K., Banks, S.A., 2005. Computational wear prediction of a total knee replacement from in vivo kinema tics. Journal of Biomechanics 38, 305314. Fregly, B.J., Reinbolt, J.A., Rooney, K.L., Mitc hell, K.H., Chmielewski, T.L., 2007. Design of patientspecific gait modifications for knee osteoa rthritis rehabilitation. IEEE Transactions on Biomedical Engineering 54, 16871695. Fregly, B.J., Banks, S.A., D'Lima, D.D., Colw ell, C.W., Jr., 2008. Sensitivity of knee replacement contact calculations to kinematic measurement errors. Journal of Orthopaedic Resarch 26, 11739. Gradin, M., Cardona, A., 2001. Flexible multibody dynamics : a finite element approach. John Wiley, Chichester ; New York. Gerstmayr, J., Schoberl, J., 2006. A 3D finite element method for flexible multibody systems. Multibody System Dynamics 15, 309324. Gibson, N., Graham, H.K., Love, S., 2007. Botu linum toxin A in the management of focal muscle overactivity in children with cerebral palsy. Disability and Rehabilitation 29, 18131822. Giddings, V.L., Kurtz, S.M., Edidin, A.A., 2001. Total knee replacement polyethylene stresses during loading in a knee simulator. Journal of Tribology 123, 842847. Girosi, F., 1998. An equivalence between spar se approximation and support vector machines. Neural Computation 10, 14551480. Giunta, A.A., Balabanov, V., Burgee, S., Grossm an, B. Haftka, R.T., Mason, W.H., Watson, L.T., 1997. Multidisciplinary optimisation of a s upersonic transport using design of experiments theory and response surface modeling. Aeronautical Journal 101, 347356. Glitsch, U., Baumann, W., 1997. The threedimensi onal determination of internal loads in the lower extremity. Journal of Biomechanics 30, 11231131. Glocker, C., 1999. Formulation of spatial c ontact situations in rigid multibody systems. Computer Methods in Applied Mech anics and Engineering 177, 199214. Godest, A.C., Meaugonin, M., Haug, E., Taylor, M., Gregson, P.J., 2002. Simulation of a knee joint replacement during a gait cycle using explicit finite el ement analysis. Journal of Biomechanics 35, 267276. PAGE 90 90 Gollehon, D.L., Tozilli, P.A., Warren, R.F., 1987. The role of the poste rolateral and cruciate ligaments in the stabil ity of the human knee: a biomechani cal study. The Journal of Bone and Joint Surgery 69, 233242. Gupta, A., Ding, Y., Xu, L., Reinikainen, T., 2006. Optimal parameter selection for electronic packaging using sequential computer simula tions. Journal of Manufacturing Science and Engineering 128, 705715. Halloran, J.P., Easley, S.K., Petrella, A.J., Rullkoetter, P.J., 2005a. Comparison of deformable and elastic foundation finite element simulations for predicting knee re placement mechanics. Journal of Biomechanical Engineering 127, 813818. Halloran, J.P., Petrella, A.J., Rullkoetter, P.J., 2005b. Explicit finite el ement modeling of total knee replacement mechanics. Journal of Biomechanics 38, 323331. Hammersley, J.M., 1960. Related Problems .3. MonteCarlo Methods fo r Solving Multivariable Problems. Annals of the New York Academy of Sciences 86, 844874. Hawary R.E., Roth, S.E., Harwood, J.C., Johnson, J.A., King, G.J.W., Chess, D.G., 2003. Ligamentous balancing in total kn ee arthroplasty: an invitro load cell analysis. Journal of Bone and Joint Surgery British Volume 90. Higginson, J.S., Zajac, F.E., Neptune, R.R., Kaut z, S.A., Delp, S.L., 2006. Muscle contributions to support during gait in an indi vidual with poststroke hemiparesis. J Biomech 39, 176977. Hippmann, G., 2004. An algorithm for compliant contact between comp lexly shaped bodies. Multibody System Dynamics 12, 345362. Hosder, S., Watson, L.T., Grossman, B., Mason, W.H., Kim, H., Haftka, R.T., Cox, S.E., 2001. Polynomial Response Surface Approximations for th e Multidisciplinary Design Optimization of a High Speed Civil Transport. Optim ization and Engineering 22, 431452. Hu, N., 1997. A solution method for dynamic cont act problems. Computers & Structures 63, 10531063. Hurwitz, D.E., Ryals, A.R., Block, J.A., Shar ma, L., Schnitzer, T.J., Andriacchi, T.P., 2000. Knee pain and joint loading in subjects with oste oarthritis of the knee. Journal of Orthopaedic Research 18, 572579. Ishikawa, M., Komi, P.V., Grey, M.J., Lepola, V., Bruggemann, G.P., 2005. Muscletendon interaction and elastic energy usage in human walking. Journal of Applied Physiology 99, 603608. Jin, R., Chen, W., Simpson, T.W., 2001. Comp arative studies of me tamodelling techniques under multiple modelling criteria. Structural and Multidisciplinary Optimization 23, 113. Jinha, A., AitHaddou, R., Herzog, W., 2006. Predic tions of cocontraction depend critically on degreesoffreedom in the musculoskeletal model. Journal of Biomechanics 39, 11451152. PAGE 91 91 Johnson, T.S., Laurent, M.P., Yao, J.Q., Gilber tson, L.N., 2001. The e ffect of displacement control input parameters on tibiofemoral prosthetic knee wear. Wear 250, 222226. Jones, D.R., Schonlau, M., Welch, W.J., 1998. Efficient global optimi zation of expensive blackbox functions. Journal of Global Optimization 13, 455492. Kalagnanam, J.R., Diwekar, U.M., 1997. An effi cient sampling technique for offline quality control. Technometrics 39, 308319. Kane, T.R., Levinson, D.A., 1985. Dynamics, theo ry and applications. McGrawHill, New York. Kaufman, K.R., An, K.N., Litchy, W.J., Chao, E. Y.S., 1991. Physiological Prediction of Muscle Forces .2. Application to Isokinetic Exercise. Neuroscience 40, 793804. Kaufman, M.D., Balabanov, V., Burgee, S.L., Giunt a, A.A., Grossman, B., Haftka, R.T., Mason, W.H., and Watson, L.T., 1996. Variablecompl exity response surface approximations for wing structural weight in HS CT design. Computationa l Mechanics 18, 112126. Kawanabe, K., Clarke, I.C., Tamura, J., Akagi, M., Good, V.D., Williams P.A., Yamamoto, K., 2001. Effects of AP translati on and rotation on the wear of UHMWPE in a total knee joint simulator. Journal of Biomedical Materials Research 54, 400406. Khuri, A.I., Cornell, J.A., 1996. Response surface s : designs and analyses. Marcel Dekker, New York. Kim, H.J, Fernandez, J.W., Akbarshahi, M., Walter, J.P., Fregly, B.J., Pandy. M.G., 2008. Evaluation of predicted kneejoint muscle forces during gait using an instrumented knee implant. Journal of Orthopaedic Research (in press). Knight, L.A., Pal, S., Coleman, J.C., Bronson, F., Haider, H., Levine, D.L., Taylor, M., Rullkoetter, P.J., 2007. Comparison of longterm numerical and experimental total knee replacement wear during simulated gait load ing. Journal of Biomechanics 40, 15501558. Koch, P.N., Simpson, T.W., Allen, J.K., Mist ree, F., 1999. Statistical approximations for multidisciplinary design optimization: The problem of size. Journal of Aircraft 36, 275286. Kurtaran, H., Eskandarian, A., Marzougui, D ., Bedewi, N.E., 2001. Crashworthiness design optimization using successive response surface approximations. Computational Mechanics 29, 409421. Kwak, S.D., Blankevoort, L., Ateshian, G.A., 2000. A mathematical formulation for 3D quasistatic multibody models of diarthrodial joints. Computer Methods in Biomechanics and Biomedical Engineering 3, 4164. LaPrade, R.F., Tso, A., Wentorf, F.A., 2004. Force measurements on the fibular collateral ligament, popliteofibular ligament, and popliteus tendon to applied loads. American Journal of Sports Medicine 32, 16951701. PAGE 92 92 Li, B., Melkote, S.N., 1999. An elastic contact model for the prediction of workpiecefixture contact forces in clamping. Journal of Ma nufacturing Science and Engineering 121, 485493. Li, G., Sakamoto, M., Chao, E.Y.S., 1997. A compar ison of different methods in predicting static pressure distribution in articulating join ts. Journal of Biomechanics 30, 635638. Li, G., Kaufman, K.R., Chao, E.Y.S., Rubash, H.E., 1999. Pred iction of antagonistic muscle forces using inverse dynamic opt imization during flexion extension of the knee. Journal of Biomechanical Engineering 121, 316322. Lin, Y.C., Farr, J., Carter, K., Fregly, B.J., 2006. Response surface optimization for joint contact model evaluation. Journal of A pplied Biomechanics 22, 120130. Liu, B., Haftka, R.T., Akgun, M.A., 2000. Twolevel composite wing structural optimization using response surfaces. Structural and Multidisciplinary Optimization 20, 8796. Lophaven, S.N., Nielsen, H.B., Sndergaard, J., 2002. DAC EA Matlab Kriging Toolbox. Technical Report IMMTR200213, Informatics and Mathematical Modelling, Technical University of Denmark. Mack, Y., Goel, T., Shyy, W., Haftka, R.T ., 2007. Surrogate ModelBased Optimization Framework: A Case Study in Aerospace Design. Studies in Computational Intelligence 51, 323342. Mancuso, C.A., Sculco, T.P., Wickiewicz, T. L., Jones, E.C., Robbins, L., Warren, R.F., WilliamsRusso, P., 2001. Patients' expectations of knee surgery. Journa l of Bone and Joint SurgeryAmerican Volume 83A, 10051012. Martens, H., Naes, T., 1984. Multi variate Calibration. Chemometrics, Mathematics and Statistics in Chemistry. Milner, C.E., Hamill, J., Davis, I., 2007. Are knee mechanics during early stance related to tibial stress fracture in runners? Clinical Biomechanics 22, 697703. Moran, M.F., Bhimji, S., Racanelli, J., Pi azza, S.J., 2008. Computational assessment of constraint in total knee replacement. Journal of Biomechanics 41, 201320. Myers, R.H., Montgomery, D.C., 1995. Respons e surface methodology : process and product in optimization using designed e xperiments. Wiley, New York. Bunderson, N.E., Burkholder, T.J., Ting, L. H., 2008. Reduction of neuromuscular redundancy for postural force generation using an intrinsic stability criterion. Journal of Biomechanics 41,15371544. Neptune, R.R., Wright, I.C., van de n Bogert, A.J., 2000. Influence of orthotic devices and vastus medialis strength and timing on pa tellofemoral loads during runni ng. Clinical Biomechanics 15, 611618. PAGE 93 93 Noble, P.C., Gordon, M.J., Weiss, J.M., Reddix, R.N., Conditt, M.A., Mathis, K.B., 2005. Does total knee replacement restore normal knee f unction? Clinical Orthopaedics and Related Research, 157165. Pal, S., Haider, H., Laz, P.J., Knight, L.A., Rullkoetter, P.J., 2008. Probabilistic computational modeling of total knee replacem ent wear. Wear 264, 701707. Pandy, M.G., Sasaki, K., and Kim, S., 1997. A th reedimensional musculoskeletal model of the human knee joint. Part 1: Theoretical construction. Computer Methods in Biomechanics and Biomedical Engineering 1, 87108. Pandy, M.G., Sasaki, K., Kim, S., 1998. A thr eedimensional musculoskeletal model of the human knee joint. part 1: theoretical constr uct. Computer Methods in Biomechanics and Biomedical Engineering 1, 87108. Piazza, S.J., Delp, S.L., 2001. Threedimensiona l dynamic simulation of total knee replacement motion during a stepup task. Journal of Biomechanical Engineering 123, 599606. Pierce, J.E., Li, G., 2005. Muscle forces predic ted using optimization methods are coordinate system dependent. Journal of Biomechanics 38, 695702. Potra, F.A., Anitescu, M., Gavreal, B., Trinkle, J., 2006. A linearly imp licit trapezoidal method for integrating stiff multibody dynami cs with contact, joints, and fr iction. International Journal for Numerical Methods in Engineering 66, 10791124. Queipo, N.V., Pintos, S., Rincon, N., Cont reras, N., Colmenares, J., 2002. Surrogate modelingbased optimization for th e integration of static and dynamic data into a reservoir description. Journal of Petroleum Science and Engineering 35, 167181. Queipo, N.V., Haftka, R.T., Shyy, W., Goel T., Vaidyanathan, R., Tucker, P.K., 2005. Surrogatebased analysis and optimization. Progress in Aerospace Sciences 41, 128. Rawlinson, J.J., Furman, B.D., Li, S., Wright, T. M., Bartel, D.L., 2006. Retrieval, experimental, and computational assessment of the perfor mance of total knee replacements. Journal of Orthopaedic Research 24, 13841394. Reeves, N.D., Maganaris, C.N., Narici, M.V., 2003. Effect of strength training on human patella tendon mechanical properties of older indivi duals. Journal of PhysiologyLondon 548, 971981. Reinbolt, J.A., Haftka, R.T., Chmielewski, T.L. Fregly, B.J., 2008. A co mputational framework to predict posttreatment outcome for gaitrelate d disorders. Medical E ngineering & Physics 30, 434443. Reinbolt, J.A., Schutte, J.F., Freg ly, B.J., Koh, B.I., Haftka, R.T ., George, A.D., Mitchell, K. H., 2005. Determination of patientspecific multijoint kinematic models through twolevel optimization. Journal of Biomechanics 38, 621. PAGE 94 94 Roux, W.J., Stander, N., Haft ka, R.T., 1998. Response surface approximations for structural optimization. International Journal for Nu merical Methods in E ngineering 42, 517534. Sacks, J., Welch, W.J., Mitchell, T.J., Wynn, H.P., 1989. Design and analysis of computer experiments. Statistical Science 4, 409435. Santamaria, J., Vadillo, E.G., Gomez, J., 2006. A comprehensive method for the elastic calculation of the twopoint wheelrail contact. Vehicle Sy stem Dynamics 44, 240250. Schmidt, D.J., Arnold, A.S., Carroll, N.C., Delp S.L., 1999. Length changes of the hamstrings and adductors resulting from derotational osteotomies of the femur. Journal of Orthopaedic Research 17, 279285. Schuind, F., GarciaElias, M., C ooney, W.P., 3rd, An, K.N., 1992. Flexor tendon forces: in vivo measurements. Journal of Hand Surgery 17, 291298. Schultz, A.B., Anderson, G.B.J., 1981. Analysis of loads on the lumbar spine. Spine 6, 76. Sharf, I., Zhang, Y.N., 2006. A contact for ce solution for noncol liding contact dynamics simulation. Multibody System Dynamics 16, 263290. Simpson, T.W., Peplinski, J.D., Koch, P.N., A llen, J.K., 2001. Metamodels for computerbased engineering design: survey and recommendations. Engineeri ng with Computers 17, 129150. Smith, M., 1993. Neural Networks for Statistica l Modeling. John Wiley & Sons, Inc. New York, NY, USA. Stansfield, B.W., Nicol, A.C., Paul, J.P., Kelly I.G., Graichen, F., Be rgmann, G., 2003, Direct comparison of calculated hip join t contact forces with those measured using instrumented implants. An evaluation of a threedimensional ma thematical model of the lower limb. Journal of biomechanics 36, 929936. Stewart, D.E., 2000. Rigidbody dynamics with friction and impact. Si am Review 42, 339. Streator, J.L., 2003. Dynamic contact of a rigid sp here with an elastic halfspace: A numerical simulation. Journal of Tribology 125, 2532. Sujith, O.K., 2008. Functional electrical stimulati on in neurological disorders. European Journal of Neurology 15, 437444. Taylor, W., Heller, M., Bergmann, G., Duda, G. 2004. Tibiofemoral loading during human gait and stair climbing. Journal of Or thopaedic Research 22, 625632. ValeroCuevas, F.J., Johanson, M.E., Towles, J.D., 2003. Towards a realistic biomechanical model of the thumb: the choice of kinematic desc ription may be more critical than the solution method or the variability/uncertainty of musculoskeletal parameters. Journal of Biomechanics 36, 10191030. PAGE 95 95 van den Bogert, A.J., Smith, G.D., Nigg, B.M., 1994. In vivo determination of the anatomical axes of the ankle joint complex: an optim ization approach. Journal of Biomechanics 27:1477. Vapnik, V.N., 1998. Statistical l earning theory. John Wiley & Sons. Wahba, G., 1987. Spline models for observationa l data. CBMSNSF Regional Conference Series in Applied Mathematics, Based on a series of 10 lectures at Ohio State University at Columbus, March 2327, 1987, Philadelphia: Society for Industrial and Applied Mathematics, 1990. Walker, P.S., Blunn, G.W., Broome, D.R., Perry, J ., Watkins, A., Sathasivam, S., Dewar, M.E., Paul, J.P., 1997. A knee simulating machine for performance evaluation of total knee replacements. Journal of Biomechanics 30, 8389. Wang, G.G., Shan, S., 2007. Review of metamo deling techniques in s upport of engineering design optimization. Journal of Mechanical Design 129, 370380. Ward, S.R., Powers, C.M., 2004. The influence of patella alta on patellofemoral joint stress during normal and fast walking. Cl inical Biomechanics 19, 10401047. Weiss, J.M., Noble, P.C., Conditt, M.A., Kohl, H.W., Roberts, S., Cook, K.F., Gordon, M.J., Mathis, K.B., 2002. What functional activities ar e important to patients with knee replacements? Clinical Orthopaedics and Related Research, 172188. Wendland, H., 1995. Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree. Advances in Computational Ma thematics 4, 389396. Wilson, D.R., Feikes, J.D., O'Connor, J.J., 1998. Ligaments and articular contact guide passive knee flexion. Journal of Biomechanics 31, 112736. Wu, Y.T., Millwater, H.R., Cruse, T.A., 1990. Advanced probabilistic structuralanalysis method for implicit performance functions. AIAA Journal 28, 16631669. Zhao, D., Banks, S.A., D'Lima, D.D., Colwell, C.W., Fregly, B.J., 2007. In vivo medial and lateral tibial loads during dynamic and high flexion activities. Journal of Orthopaedic Research 25, 593602. Zhao, D., Banks, S.A., Mitchell, K.H., D'Lima D.D., Colwell, C.W., Fregly, B.J., 2007. Correlation between the knee adduc tion torque and medial contact force for a variety of gait patterns. Journal of Orthopaedic Research 25, 789797. Zhao, D., Sakoda, H., Sawyer, W.G., Banks, S.A., Fregly, B.J., 2008. Predicting knee replacement damage in a simulator machine using a computational model with a consistent wear factor. Journal of Biomecha nical Engineering 130, 011004. PAGE 96 96 BIOGRAPHICAL SKETCH YiChung Lin was born in Tainan, a city in southern Taiwan, in 1976. During his childhood, he enjoyed assem bling plastic models. He has assembled many different varieties of plastic models (aircraft, automobile, robot), and gradually found himself a bi g fan of robots. This hobby eventually led him into the Department of Mechanical Engineering of the National Cheng Kung University, Taiwan. He received his Bachelor of Science in 1998 and joined the Army to fulfill the military service obligation. After rece iving a Honorable Discharge, he went overseas and enrolled in the masters program in the Department of Mechanical and Aerospace Engineering at the University of Florida. During his masters program, he extended his interest from robots to the human body and joinined the Computational and Biom echanics Laboratory of Associate Professor B. J. Fregly in fall 2002. In 2004, YiChung received a Master of Science degree in mechanical engineering with focus on biom echanics. At that time, he decided to stay in the same lab to pursue the Doctor of Philos ophy degree. YiChung would like to continue his research in biomechanics and to improve the quality of human daily life in the future. 