A Hertzian Integrated Contact Model of the Total Knee Replacement Implant for the Estimation of Joint Contact Forces

The prediction of lower limb muscle and contact forces may provide useful knowledge to assist the clinicians in the diagnosis as well as in the development of appropriate treatment for musculoskeletal disorders. Research studies have commonly estimated joint contact forces using model-based muscle force estimation due to the lack of a reliable contact model and material properties. The objective of this present study was to develop a Hertzian integrated contact model. Then, in vivo elastic properties of the Total Knee Replacement (TKR) implant were identified using in vivo contact forces leading to providing reliable material properties for modeling purposes. First, a patient specific rigid musculoskeletal model was built. Second, a STL-based implant model was designed to compute the contact area evolutions during gait motions. Finally, a Hertzian integrated contact model was defined for the in vivo identification of elastic properties (Young’s modulus and Poisson coefficient) of the instrumented TKR implant. Our study showed a potential use of a new approach to predict the contact forces without knowledge of muscle forces. Thus, the outcomes may lead to accurate and reliable prediction of human joint contact forces for new case study.


Introduction
Total Knee Replacement (TKR) implant has been commonly prescribed for patients with severe knee joint damage (e.g., osteoarthritis) associated with progressive pain and impaired function [1,2]. The implant with distal femoral and proximal tibial compartments is usually created with biocompatible materials [3,4]. Besides, TKR implant may include patellafemoral interaction which helps the knee joint to be functionally stabilized. It is important to emphasize that the knee contact force behaviors are activity-specific and the optimal design of the TKR may lead to maximizing the clinical outcomes. However, the implant design may be optimized when its effect on the musculoskeletal tissues and structures could be elucidated [5]. Joint contact forces play an important role in the understanding of the effect of the implant on the joint behavior and musculoskeletal load [6]. Finite element modeling may provide reliable information on the joint contact force behavior. However, the lack of accurate muscle modeling and reliable definition of boundary and loading conditions may affect the contact force computation [7,8].
Research studies demonstrated a deeper knowledge of hip implant effect thanks to computational rigid-body models and available in vivo experimental data measured directly within the instrumented hip implant [9][10][11]. Recently, in vivo knee contact forces can be measured from a new knee device [12,13]. Thus, for the first time, computational rigid musculoskeletal models, used to facilitate the surgical and functional rehabilitation treatments of the musculoskeletal disorders, may be evaluated for their predictive capacities related to muscle force estimation and knee contact load by using these available experimental data [14,15]. Thus, these models may be potentially accepted in a clinical routine procedure in the future [9,16]. Moreover, an open and free access database (http://www.orthoload.com/) of the orthopedic implant loads was also provided leading to promoting extensive research activities in this field of study.
The prediction of in vivo muscle and contact forces, especially within a Total Knee Replacement, is a challenging task for biomechanical engineers/researchers. Model-based 2 Journal of Computational Medicine estimation of these forces has been commonly performed [17][18][19]. Thus, the redundant phenomena of the human neural control (i.e., the number of muscles is greater than the number of degrees of freedom) may be computationally solved using optimization techniques even if the best optimization strategy remains unknown. However, the prediction accuracy of such forces depends on the accurate modeling of musculoskeletal tissues (e.g., muscles, bone, tendons, ligaments, and cartilage) and structures (e.g., joint) as well as appropriate formulation of physical behaviors (e.g., joint contact or multi-body dynamics) [19,20]. Model validation has been a challenging issue. Qualitative pattern of the muscle force was commonly validated against the EMGbased muscle activities.
In vivo knee contact forces may be computationally estimated with or without contact model. On the one hand, model-based efforts related to the use of inverse dynamics approach coupled with an elastic contact model. Within the scope of this vision, Kim et al. [21] developed twostage simulation to compute the muscle forces and then combine them with ground reaction forces and fluoroscopic data within a Hertzian contact model to predict the knee contact forces. Moreover, Hast and Piazza [22] used a coupling between inverse dynamics, forward dynamics, and computed muscle control (CMC) algorithm within a springbased contact model to compute the nodes-to-surface knee contact forces. On the other hand, knee contact forces have been predicted without contact model. Based on the muscle forces estimated from pseudoinverse method and constrained static optimization with Lagrange multipliers formulation and parameter reduction strategy [23], contact forces were computed. Moreover, Electromyography-(EMG-) driven muscle force estimation approach coupled with moment balancing algorithm has been used to predict the medial contact force [24,25]. Furthermore, a parametric numerical model including 9 equilibrium equations reflecting muscle activation levels, passive structure activation, joint moments, external joint forces, maximum physiological muscle moments, and forces has been also developed to predict the knee contact forces [26]. In fact, without contact model, the problem formulation and computational cost may be more advantageous than the use of a contact model. A recent comparison of predicted contact forces in the framework of the Knee Grand Challenge showed that the use of a contact model does not guarantee their accuracy [15]. However, a deformable contact model may provide useful knowledge about the interactions between the implant components and also those at bone-to-implant surface. It is important to note that the Knee Grand Challenge has been created to open a new dimension on the prediction of in vivo muscle and contact forces [14,15]. Almost all available measurements ranging from medical imaging data to gait analysis and experimental contact forces were provided to enhance the research activities related to the prediction of these forces.
Previous studies used only simplified representation (e.g., one point contact or nodes-to-surface contact) of contact surface between implant components. Moreover, material properties of the implants were commonly provided by manufacturer. However, assumptions were commonly performed for musculoskeletal modeling. Thus, these properties need to be calibrated and fitted to the modeling assumption for predicting accurately the knee contact forces. In fact, predictive trends of numerical models are limited for clinical applications. Thus, the objectives of this present study are twofold: (1) the development of a Hertzian integrated contact model based on a global surface-to-surface interaction and (2) the use of in vivo contact forces and the developed contact model to identify the reliable and fitted elastic properties of the implant for modeling purpose.

Computational Workflow.
A computational workflow was developed to identify the elastic properties of the instrumented TKR implant (Figure 1). The workflow consists of the following components: (1) a patient specific rigid musculoskeletal model to compute the knee kinematics from skin-mounted markers' trajectories; (2) a STL-based implant model to compute the contact area evolution during knee motion; and (3) a Hertzian integrated contact model for the in vivo identification of elastic properties (Young's modulus and Poisson coefficient) of the instrumented TKR implant using contact area information and in vivo measured joint contact forces.

Patient Data.
The patient data from the fourth edition of the Knee Grand Challenge [14] was used in this present study. Imaging and gait data were acquired on a patient (male, 88 years old, 168 cm body height, 66.7 kg body mass, and 23.6 kg/m 2 Body Mass Index (BMI)) with knee implant instrumented in his right knee side (Figure 2(a)). The patient had Rockport flat bottom sneakers shoes during the gait data acquisition process. The patient performed an overground normal gait pattern, which was used as reference pattern. Then, 4 overground gait trials reflecting 4 gait modification strategies (bouncy, medial thrust, mild crouch, and mtp) were acquired using an 8-camera Vicon system (Table 1). All skinmounted markers were posed using a modified Cleveland Clinic marker protocol [14]. In vivo knee contact forces were also obtained during the gait data acquisition. Four uniaxial load cell measurements (posterior-medial (PM), anteriormedial (AM), anterior-lateral (AL), and posterior-lateral (PL)) from the instrumented tibial prosthesis were performed ( Figure 2(b)). The first generation of tray design (eKnee) was used. The main implant components were fabricated in CoCr, polyethylene, and titanium. Sampling frequency is 120 Hz. Furthermore, pre-and postsurgery CT data were also acquired on the patient using a Siemens CT machine; the slice thickness was set up as 0.6 mm. The matrix is 512 × 512 and the pixel resolution is 0.9 × 0.9 mm 2 . The number of slices is 1900 for the whole scanned lower limb.

Model Development.
To compute the knee kinematics, a patient specific musculoskeletal model, which is provided   from the fourth edition of the Knee Grand Challenge [14], was used. This model has 8 segments (pelvis, right femur, femoral implant component, patella, tibial implant component, right tibia, right fibula, and right foot). Image segmentation using SliceOmatic (Montreal, Canada) was performed on the post-and presurgery CT data and laserscan-based implant component to extract bony segments and implant components. STL-based surface geometries were created and imported into OpenSIM engine [27] to develop corresponding bony segments. The model has 24 degrees of freedom (DOF) (6-DOF ground-to-pelvis joint, 3-DOF hip joint, 6-DOF tibiofemoral joint, 6-DOF patellofemoral joint, 2-DOF ankle joint, and 1-DOF toes joint). The interaction between implant components and its respective bone was modeled using an OpenSIM weld joint. The interaction between implant components (femoral component, tibial tray, and patella/femoral component) was modeled using an OpenSIM reverse joint. All segment coordinate systems were based on the ISB recommendations [28]. The model has 22 skin-mounted virtual markers' cluster according to the real market cloud used in the motion capture acquisition. This model has 44 lower limb muscles. Each muscle was modeled as a Schutte muscle model [29]. Muscle attachment points and wrapping surface were created using an available lower limb model [30]. A transformation process using information from bone-to-bone alignment was performed to morph the patient specific musculature. Some attachments points and unrealistic intersection of muscle lines of action were manually adjusted.

Joint Angle
Computing. Based on skin-mounted marker's trajectories, inverse kinematics was performed to compute the knee joint angles of each gait pattern (normal, bouncy, medial thrust, mild crouch, and mtp). At each time step (i.e., frame of motion), generalized coordinate values were computed to put the model in a best-fitting pose in minimizing the deviation between the experimental markers' and estimated coordinate values. The deviation is expressed as a weighted least squares problem as follows: where is the vector of generalized coordinates; is the marker weight; exp is the experimental position of the marker ; ( ) is the position of the corresponding marker on the model.

STL-Based Implant Model and Contact Area Computing.
Based on the reconstructed implant components, a geometrical computing process was applied to compute the surface-tosurface contact area evolution according to each series of knee kinematics (e.g., flexion/extension). Threshold-based minimal distance principle was used to compute the implant-toimplant contact area: minimal distance map was computed and area under the threshold of 4 mm was calculated [31]. The threshold was calibrated using experimental data of the knee contact area [32]. An illustration of contact area map during knee flexion/extension from 0 to 90 degrees is shown in Figure 3.

Contact Model.
To describe the physical relationships between the applied force and the contact behavior, a Hertzian linearly elastic contact model was formulated. This model describes conservative interaction forces [33]. Constitutive equation of this model is expressed as follows: where Hertz and Hertz are applied force and its respective contact radius. is the radius of the contact cylinder (e.g., femoral condyle) ( Figure 4). and V are Young's modulus and Poisson coefficient, respectively. It is assumed that is constant for different knee flexion/extension angles.

Elastic Property Identification.
To identify the appropriate values of the implant elastic material properties ( and V), an optimization problem was formulated as follows: Subject to: > 0, 0 < V < 0.5 where eKnee is the in vivo measured contact force. is a set of knee angles. is the total contact area computed by the method described in Section 2.4. To compute the value of , a cylinder-fitting process was applied. The values of in vivo measured contact forces were calculated from the load cell measurements ( AM , PM , AL , PL ) (Figure 2 The constrained optimization process was performed using a quasi-Newton approach implemented in the Optimization Toolbox in Matlab R2010b (Mathworks, USA).

Determination of Contact Cylinder
Radius . The fitted cylinder used to compute the value of the radius is shown in Figure 5. It is important to note that two physiological condyles were assumed to be symmetric. Thus, the knee joint behavior functions like a hinge joint. A value of 37 mm was found with a mean error of 16.52 mm from 1024 points between the fitted cylinder and the external femoral surface of the implant component. Note this error is relatively high due to the assumption of constant for different knee flexion/extension angles. Figure 6 shows the evolution of knee flexion/extension angles during the normal gait with instrumented TKR implant. The total squared error between real markers' coordinates and virtual makers' coordinates computed by inverse kinematics process ranges from 0.00072 to 0.00215 mm for 326 frames. The same pattern was obtained for other gait patterns with a first peak-to-peak absolute deviation ranging from 5.6 to 18.5 degrees as well as a second peak-to-peak absolute deviation ranging from 1.8 to 8.52 degrees. It is important to note that these deviations are due to temporal shifts when modifying the gait patterns.

Implant-to-Implant Contact Area Evolution.
The evolution of implant-to-implant contact area during the stance phase (60%) of the normal and modified gaits is illustrated in Figure 7. The peak contact area values range from 1016 to 1117 mm 2 for all gait patterns. According to the normal gait, the patterns of the contact area are different for all modified strategies, except for the Medthrust. These differences are directly related to the differences in the kinematics patterns of these modified strategies.

In Vivo Contact Forces and Implant's Elastic Properties.
The in vivo medial contact forces during the normal and 4 modified gaits are shown in Figure 8. Root mean square errors (RMSE) between the normal gait and bouncy, medial thrust, mild crouch, mtp gaits are 222. 5, 192.8, 132.3, and 1777.7 N, respectively, for the medial side. On the other hand,  lateral side revealed a RMSE of 120.2, 122.6, 91.1, and 89.6 N for comparison between the normal and the bouncy, medial thrust, mild crouch, and mtp gaits, respectively. The maximal RMSE of the medial contact force is 222.5 N between the normal and the bouncy gaits. The maximal RMSE of the lateral contact force is 122.6 N between the normal and the medial thrust gaits. Using the total contact forces with the integrated contact model, the optimization process was converged after 6 iterations (3,15,30,33,36,39, and 42 function calls, resp.) at the minimal objective function value of 1.16437 − 009. Thus, optimal Young's modulus was 4.67 GPa and its respective Poisson coefficient was 0.25. Figure 9 shows the evolutions of predicted contact forces using integrated contact model and experimental eKnee data for the normal gait. A good convergence (i.e., error between predicted results and experimental data leads nearly to 0) was  noted over the 54% of the gait cycle. However, the amplitudes and trends of the predicted values were largely different in early phase of the gait cycle. This is due to the modeling assumption of the ideal symmetric contact model.

Discussion
The prediction of in vivo muscle and contact forces plays an important role in the objective and precise diagnosis and treatment of the musculoskeletal disorders in the future [9,10,14,35,36]. In fact, the use of computational models may allow the clinicians to perform a series of virtual treatment trials without any dangerous risk for the involved patients. Moreover, optimal treatment plan may be elucidated with objective data related to its effect on the musculoskeletal tissues and structures. However, the validation of model behavior/prediction against critical experimental data is extremely a challenging task, especially for in vivo estimated muscle and contact forces. Skeletal muscle with its multiscale architecture revealed a complex (e.g., anisotropic or heterogeneous) living material. Hill-based macroscopic rheological model has been commonly used to define the physiological force-length and force-velocity relationships [37]. Despite many efforts during the last several decades to estimate the in vivo muscle forces using inverse dynamics and optimization or EMG-driven approach [17,19,20,24,25,29], the lack of in vivo experimental data for direct and quantitative validation purpose may be one of the obstacles to promote the use of computational musculoskeletal models in the clinical setting. Recently, technological progress allows measuring in vivo contact forces within instrumented implants (e.g., at the hip and knee levels) [9,12,13], which could be used to validate indirectly the muscle force estimation. However, these data are only available for patients with instrumented implant. Thus, the accurate and reliable prediction and its quantitative validation of these forces on the healthy subjects to prevent the musculoskeletal diseases remain a challenge for musculoskeletal modeling community.
Regarding the Knee Grand Challenge, a series of editions has been organized and many efforts have been investigated to predict the in vivo contact forces using different approaches [15,[21][22][23][24][25][26]. In fact, contact forces may be predicted with or without a contact model formulation. Even if the approach without a contact model revealed interesting results [26], this approach is based purely on the numerical power to find the statistical correlation pattern between the contact forces and other physically measurable or computable quantities such as EMG-based muscle activation levels or maximum physiological muscle moments. Besides, the use of a contact model may provide more descriptive and meaningful information inside the implant-to-implant interaction under musculoskeletal load. However, current approach based on Hertzian contact model has two disadvantages. The first one relates to the strong constraint of the use of muscle forces to predict the contact forces due to the lack of contact area information. The second one deals with the simplified representation of the contact interaction. Thus, nodes-to-surface approach is commonly used. However, the number and location of nodes may be sensitive factors influencing the prediction of the contact forces. This present study developed a Hertzian contact model integrating contact area information. Indeed, muscle forces are not necessary to predict the contact forces. Moreover, thanks to available in vivo contact forces, elastic properties (Young's modulus and Poisson coefficient) of the implant were numerically identified and then may be used for the further prediction of the new patient case. Furthermore, when the reliable contact forces are known, muscle forces may be further estimated in a reliable manner. In addition, the knowledge of the contact area may allow other contact properties such as contact pressure and surface energy to be computed, leading to better understanding of the effect of TKR implant on the musculoskeletal tissues and structures.
From methodological point of view, the fact that muscle forces are not required to predict the contact forces leads to avoiding the limitations of the muscle force estimation approaches. For example, inverse dynamics and static optimization showed high residual results due to dynamic inconsistencies between ground reaction forces (GFR) and captured kinematics [17,19]. Kinematical errors coming from anatomical representation, marker registration, and soft tissue artifacts may alter the accuracy of the results [12]. Moreover, the inaccuracy may come from the noise in GRF or drift effect over time. Furthermore, the choice of the objective function may create nonuniqueness solutions. EMG-driven approach has computational complexity. In addition, the accuracy of estimated muscle forces depends on the accuracy of EMG signals related to electrode type and location. Furthermore, EMG signal processing approach may alter significantly the simulation results [38].
Regarding our identification process, the value of cylinder radius ( ) fitted to the femoral implant component is very close to the value used in the literature ( = 40 mm) [21]. Moreover, a good agreement was also noted between our contact area and the literature value acquired from cadaveric testing with 1000 N of load ( Figure 8) [32]. The contact evolution maps revealed a more descriptive surface-to-surface interaction between implant components. 8

Journal of Computational Medicine
The accuracy of our prediction depends strongly on the contact area estimation. Thus, the input data such as joint angle needs to be accurately computed. The use of a patient specific musculoskeletal model without scaling step may guarantee the minimal errors between virtual marker coordinates and the real marker ones during inverse kinematics process [19]. However, the lack of knowledge about the muscles forces may lead to the neglect of the contact force balance ratio between medial and lateral sides [14]. Moreover, the present definition of the contact model is global and its behavior is symmetric. Furthermore, from computational point of view, this contact area process is dissociated from the inverse kinematics analysis. Indeed, further work will be investigated for the cosimulation of these physics with asymmetric definition of the contact model integrated side-dependent kinematic data such as joint abduction/adduction patterns. In this case, advanced computational paradigm needs to be developed [18,23]. Thus, available fluoroscopic data will provide precise and local information on the translation and rotation of the implant leading to possibly achieving this challenging goal.
In addition, the use of a Hertzian contact model allows only the description of conservative interaction forces via a linearly elastic behavior. More advanced contact models need to be investigated. For example, Derjaguin-Muller-Toporov (DMT) model [39] may provide elastic behavior and van der Waals forces outside the contact region. Moreover, Johnson-Kendall-Robert (JKR) model [40] allows defining the adhesion (i.e., short-range forces in the contact area) phenomena and this model may be used to predict infinite stress at edge of contact circle and nonconservative interaction forces. Furthermore, the model needs to be refined for other types of implants with different designs.
In conclusion, reliable elastic properties of the TKR implants were identified using a Hertzian integrated contact model and in vivo contact forces. The availability of contact area information allows the in vivo prediction of contact forces without knowledge of the muscle forces. Thus, our present study demonstrated a potential use of such approach in the accurate and reliable prediction of human joint contact forces leading to better diagnosis and a more suitable treatment prescription for patients with instrumented implants to recover locomotion function and to improve the quality of daily activities of life.