Comparison between kinetic and kinetic-kinematic driven knee joint finite element models

Use of knee joint finite element models for diagnostic purposes is challenging due to their complexity. Therefore, simpler models are needed for studies where a high number of patients need to be analyzed, without compromising the results of the model. In this study, more complex, kinetic (forces and moments) and simpler, kinetic-kinematic (forces and angles) driven finite element models were compared during the stance phase of gait. Patella and tendons were included in the most complex model, while they were absent in the simplest model. The greatest difference between the most complex and simplest models was observed in the internal-external rotation and axial joint reaction force, while all other rotations, translations and joint reaction forces were similar to one another. In terms of cartilage stresses and strains, the simpler models behaved similarly with the more complex models in the lateral joint compartment, while minor differences were observed in the medial compartment at the beginning of the stance phase. We suggest that it is feasible to use kinetic-kinematic driven knee joint models with a simpler geometry in studies with a large cohort size, particularly when analyzing cartilage responses and failures related to potential overloads.

should produce the same results if the kinetic and kinematic parameters can be quantified reliably and necessary joint structures have been considered. However, to our knowledge, it has not been shown whether the simpler kinetic-kinematic driven models yield similar results as the purely kinetic driven models for the tibiofemoral contact.
In the present study, the main aim was to compare the predictions of kinetic and kinetic-kinematic driven knee joint models in terms of contact mechanics and mechanical response of cartilage during the stance phase of gait. For this purpose, four FE models with different levels of complexity were created. The first two more complex models included the patella and surrounding tendons, and were driven using a kinetics 24 and kinetics-kinematics approach, respectively. The third model used a simplified geometry (with no patella or tendons) and was driven by kinetics-kinematics. Lastly, the fourth model used the same geometry as the third model but included the contributions of the patella and tendons through the force input. We hypothesize that the FE models with simplified geometries and inputs could produce similar responses within the knee joint with the more complicated models, and potentially demonstrate the ability of the simplified FE models to be used in studies with a high number of subjects.

Materials and Methods
The study workflow is shown in Fig. 1. Experiments. Magnetic  echo train length 32, matrix 384 × 384 and field of view 16 mm). The insertion points of cruciate and collateral ligaments, as well as patellar and quadriceps tendons were determined from the MR images. The segmentation was done using Seg3D software (NIH/NIGMS CIBC, University of Utah, UT, USA), and the accuracy of the process was ensured by an experienced musculoskeletal radiologist. Surface meshes were generated after post-processing in Mimics (Materialise, Leuven, Belgium) and converted into solid geometries using a custom Matlab script (MathWorks, Natick, MA, USA). The resulting 3D solid parts were then imported in Abaqus (v6.13-3, Dassault Systémes, Providence, RI, USA). To reduce the computation time, bones were considered as rigid bodies 24,26,27 .
Motion capture and musculoskeletal modeling. The gait of the subject was obtained using a previously established protocol 28 and is briefly described below. A 10-camera motion capture system (Vicon, Oxford, UK) and 41 retro-reflective markers were used to collect segment position data. Ground reaction force data were collected simultaneously using two in-ground force plates (AMTI, Watertown, MA, USA). The joint rotations (using Cardan sequence) and joint moments (standard inverse dynamics) were computed from a seven-segment subject-specific musculoskeletal model (Visual 3D, C-Motion, Rockville, MD). The average of three gait trials was computed and used as inputs for the FE models (Fig. 1d).
For the models that included patellar contributions directly, the quadriceps muscle force was needed and computed using OpenSim (SimTK, Stanford, CA, USA). The generic Gait2392 model 29 was scaled in order to match the generic marker positions with those from the Visual3D (C-Motion, Germantown, MD, USA) musculoskeletal model. Inverse kinematics data was exported using Visual3D and used as inputs for the musculoskeletal simulation within OpenSim. Residual Reduction Algorithm (RRA), a form of forward dynamics, was coupled with Computed Muscle Control (CMC). RRA minimizes large non-physical compensatory forces, called residuals, by altering the torso mass center of the model. This makes inverse kinematics dynamically more consistent with ground reaction force data 30,31 . CMC is a tracking controller and was used to compute lower extremity muscle forces using a static optimization algorithm. The algorithm reduces muscle force redundancies by minimizing the sum of the squares of the muscle excitations and includes muscle activation and contraction dynamics 32,33 . Several studies have shown that OpenSim-CMC produces similar muscle excitation patterns as electromyography (EMG) measurements [34][35][36][37] . The quadriceps force (Fig. 1e) was calculated as the algebraic sum of rectus femoris, vastus lateralis and intermedius, and medialis muscles, and used in the FE simulation 24 . We also want to emphasize that the contribution of other muscles was indirectly considered in the kinetic driven model by assuming muscles to absorb most of the moments 24 . Similarly, the kinetic-kinematic driven model took indirectly muscles and other structures into account by directly implemented rotations in the model (see more below). The average knee joint angles, moments and ground reaction forces were also compared against literature values to ensure that they were within physiologically relevant ranges (Supplementary material).
Finite element models. Geometries, boundary and loading conditions. Two geometries were created from the segmented solid geometries. The complex geometry (Geometry A, Fig. 1b) included tibial, femoral and patellar cartilage, menisci, collateral and cruciate ligaments, patello-femoral ligaments, quadriceps and patellar tendons. The simplified geometry (Geometry B, Fig. 1c) only contained tibial cartilage, femoral cartilage, menisci, cruciate and collateral ligaments. In studies with a high number of patients, the tissues included in the simplified geometry could represent the minimum representation needed for high accuracy in the results, while still maintaining a low degree of complexity. To compare different knee joint driving methods, four models were created: Model A uses kinetic actuation and Geometry A (Fig. 1f). This model is the most complex in this study and similar to the one detailed by Halonen et al. 24 . The implemented moments included internal-external and varus-valgus, with flexion-extension as rotation. The moment scaling was 50% of the measured values, similarly as before 24 , indicating that muscles absorb half of the measured moments. By this assumption, Halonen et al. 24 obtained a good match between modeled and experimental femoral rotations. Translational forces (anterior-posterior, medial-lateral and distal-proximal) and quadriceps muscle forces (divided in anterior-posterior and distal-proximal components) were implemented in this model (Fig. 1e). From this on, we call this model as a kinetic driven model, even though the flexion-extension angle was implemented with kinematics.
Model B uses kinetic-kinematic actuation and Geometry A (Fig. 1g). The implemented translational and quadriceps forces were the same as in Model A, while flexion-extension and internal-external motions were implemented as rotations instead of moments. The varus-valgus rotation was kept free, since the variation in motion capture for varus-valgus rotation is known to be higher than with the other two rotation directions [38][39][40] . This assumption has also been applied in several earlier studies 24,41 .
Models C and D use kinetic-kinematic actuation and Geometry B (Fig. 1h,i). For Model C, only the translational forces and rotations were implemented, similarly as for Model B. For Model D, in addition to the translational forces and rotations, the anterior-posterior component of the quadriceps force was added to the anterior-posterior component of the translational forces. This was done in order to consider the quadriceps force without a need for more complex geometry with patella and tendons. A parametric study was performed to estimate the scaling of this parameter so that knee joint translations and contact forces as well as cartilage responses match with those of the more complex Model A. The scaling factor was varied between 0.5 and 1.5 times the anterior-posterior quadriceps force (QF). The best match between this model and Model A was obtained when no scaling was applied to the QF. Thus, scaling factor of 1 was used in further analyses.
The knee joint translational forces, rotations or moments (Fig. 1d) were implemented to the femoral reference point (RP Femur , Fig. 1b,c), which was taken as the midpoint of the distance between medial and lateral femoral epicondyles 41,42 . The anterior-posterior and distal-proximal components of the quadriceps force were implemented into another reference point (RP Quad , Fig. 1b), similarly as was done before 24 . The rotations were implemented as time-dependent boundary conditions, while forces and moments were implemented as time-dependent loads. For all models, the cartilage-tibia bone interface, and the meniscus and ligament insertion points in the tibia bone were fixed for all six degrees of freedom, similarly as before 41 . Additionally for Model A and B, the patellar tendon insertion point into the tibia bone was kept fixed 24 . As for the ligament-femur insertion points, they were coupled kinematically with the femoral reference point (Fig. 1b, RP Femur ). The tendon-patella insertion point was modeled by kinematically coupling it with a reference point in the patella (Fig. 1b, RP Patella ). The contact between soft tissues was modeled using a hard pressure-overclosure and frictionless surface-to-surface discretization, similarly as before 20,26 .
Meshing and Material properties. For all models, 8-node hexahedral poroelastic (C3D8P) elements were used for cartilage and corresponding elastic elements (C3D8) for menisci. The characteristic element length of femoral, tibial and patellar cartilage, and menisci were 1.5, 0.7, 0.7 and 1 mm, respectively 24 , resulting in a three depth-wise element layers for cartilage and six depth-wise element layers for menisci.
Cartilage was modeled as a transversely isotropic poro-elastic (TIPE) material 27 , while the menisci were transversely isotropic elastic (TIE) 43,44 . Details on the material properties are shown in Table 1. A more detailed description of the TIPE and TIE materials is shown in Supplementary Materials. For this study, more complicated fibril reinforced materials were not implemented, since Klets et al. 27 showed that the present TIPE material exhibits similar behavior as the fibril reinforced poroviscoelastic material, in terms of both knee joint stress and strain magnitudes and distributions. Other motivations for using the TIPE material are its simplicity in implementation and speed, making it more suitable for clinical purposes. Similarly, the TIE formulation for menisci adequately describes the behavior of the tissue during short-term loading and can produce similar results as the corresponding poroelastic model 45,46 . This also speeds up the model solution when no fluid pressure is investigated in meniscus, only the load-carrying capacity of the tissue is considered.
The cruciate and collateral ligaments were modeled as bi-linear springs, resisting tension but not compression. For the ACL a tensile stiffness of k = 380 N/mm was used 47 , while for the PCL it was k = 200 N/mm 48 . Both MCL and LCL had a tensile stiffness of k = 100 N/mm 48 . QT and PT were implemented also with bi-linear springs having tensile stiffnesses of k = 475 N/mm and k = 545 N/mm, respectively 49 . According to Gantoi et al. 50 , a pre-strain of 5% was implemented in the ACL and PCL, while for MCL, LCL, QT and PT it was 4%. Similar to Halonen et al. 24 , the MPFL and LPFL were modeled as linearly elastic truss elements with no compressive stiffness; E = 19.1 MPa for MPFL, E = 17 MPa for LPFL and ν = 0.499 51 . The bi-linear elastic spring model for ligaments and tendons is acceptable 52-55 due to the relatively short toe-region in the ligament stress-strain response 56,57 , indicating that ligaments were mostly at the linear region during loading 54,58 . A recent study 55 also showed that the knee joint model with ligaments modeled as springs produce comparable joint contact forces with the model with fibril-reinforced poroelastic ligaments. Moreover, other studies have used the spring formulation for ligaments with similar pre-strain values as used here 19,25,59,60 . The meniscal horn attachments were modeled using springs with a total spring constant for each horn of k = 350 N/mm 61 .
Simulations. The simulations were done in three steps, similarly as before 41,62 . Briefly, in the first step all soft-tissues were brought in contact. In the second step, the knee joint was rotated to its initial position and the initial forces were applied. Here, we considered the initial orientation of the knee joint from MR images. In the third step, the kinetic and kinetic-kinematic inputs were applied (Fig. 1d-i) and soils consolidation analyses were simulated in Abaqus.
Joint reaction forces, rotations and translations were compared between the models. For evaluating cartilage responses, average values over the cartilage-cartilage contact area were computed as a function of time. Further, we visually evaluated the distributions of maximum principal stresses and maximum logarithmic strains at 20%, 50% and 80% of the stance phase. These correspond to the first peak load, mid-stance and second peak load. The computation time for each model was estimated. It includes all tasks needed for the FE model generation and simulation: MRI segmentation, meshing, implementation of the material properties, extraction of the motion data, musculoskeletal modeling (needed in Models A, B, D), implementation of the boundary conditions, and adjusting the model parameters (such as contact interactions, mesh, convergence parameters) to enable a converged solution.

Results
All four models showed similar trends in varus-valgus rotation (Fig. 2a). The average difference between the kinetic (Model A) and kinetic-kinematic (Models B, C and D) driven models over the stance phase was ~0.15°. Differences between the kinetic and kinetic-kinematic driven models were found in the internal-external rotation (Fig. 2b), the average difference over the stance phase being ~4°. In terms of translations, all models behaved quite similarly for the distal-proximal (Fig. 2c) and medial-lateral (Fig. 2d) directions, with the average differences of ~0.12 mm and ~0.11 mm between Models A and C. For the anterior-posterior direction (Fig. 2e), the maximum  difference between Model A and the rest of the models was ~2 mm or less at ~20% and ~80% of the stance. For the rest of the stance, this difference was on average ~0.2 mm. Models A and B showed almost identical behavior in joint reaction forces through the stance phase. Model C showed the smallest distal-proximal (Fig. 2f) and anterior-posterior (Fig. 2h) reaction forces, with the maximum differences of ~0.6 BW and ~0.25 BW, respectively, at ~20% of the stance when compared to Model A. These differences reduced to ~0.05 BW for the rest of the stance phase. Even though Model D only included the anterior-posterior contribution of the quadriceps force, the maximum difference in the distal-proximal reaction force, as compared to Model A, was reduced to ~0.25 BW (Fig. 2f). On the other hand, this force affected only slightly the anterior-posterior reaction force (Fig. 2h). For all models, the medial-lateral reaction force (Fig. 2g) was almost the same, with the average difference of ~0.1 BW between Models A and C throughout the stance phase.
At 20% of the stance, the kinetic driven model (Model A) showed a different stress distribution on the cartilage surface than the kinetic-kinematic driven models (Models B, C and D). However, at 80% of the stance, all models showed almost identical distribution of stresses (Fig. 3a). The main difference between the models was between 10% and 30% of the stance, where the medial contact area of Model A was different than that in Models B, C and D. After 30% of the stance, the contact areas for all models were almost the same for both medial and lateral compartments (Fig. 3b).
To compare the FE models more quantitatively, the average values of maximum principal stresses, logarithmic strains and pore pressures over the cartilage-cartilage contact area were calculated as a function of stance (Fig. 4). All kinetic-kinematic driven models (Models B, C and D) showed higher values in the lateral compartment compared to the medial compartment (~60% of the total values) from 0% to 50% of the stance and almost even distribution between the lateral and medial compartment from 50% to 100% of the stance. The kinetic driven model (Model A) revealed nearly even distribution of stresses and strains in the medial and lateral compartments throughout the entire stance phase.
In the medial compartment, the differences were up to two-fold between the kinetic driven model (Model A) and the kinetic-kinematic driven models (Models B, C and D) for all parameters at 10-30% of the stance. Onwards from ~30% of the stance, the differences in these parameters became very small, with the average differences of ~0.5 MPa, ~1% and ~1 MPa for maximum principal stresses (Fig. 4a), logarithmic strains (Fig. 4b) and pore pressure (Fig. 4c), respectively, between Models A and C. In the lateral compartment, all models behaved similarly, with the maximum differences of ~0.5 MPa, ~0.5% and ~1 MPa in maximum principal stresses (Fig. 4d), logarithmic strains (Fig. 4e) and pore pressure (Fig. 4f), respectively, between Model A and the rest of the models. Estimates of the total time needed for each model generation and simulation are shown in Table 2. The most complex, kinetic driven model (Model A) required ~10 times longer to get a converged solution, when compared to the simplest, kinetic-kinematic driven model (Model C).

Discussion
In the present study, more complex kinetic driven and simpler kinetic-kinematic driven finite element models were created. The most complex kinetic driven model included patella and tendons, in addition to other major joint structures, while the simplest kinetic-kinematic driven model was generated only to characterize the tibiofemoral contact. The effects of different implementations of gait cycle, obtained from motion analysis, on knee joint motion and tibiofemoral joint reaction forces were compared. Moreover, pore pressures, stresses and strains at the tibial cartilage contact surface were compared. In most of the analyzed parameters, differences between the models were small. Particularly, varus-valgus rotation, translations and cartilage responses in each model were highly similar. Differences in the joint reaction forces between kinetic (with patella) and kinetic-kinematic (without patella) driven models could be minimized by adding an additional force component to estimate the effect of patella and quadriceps forces. These results suggest that simpler models, in terms of motion and geometry, could be applied when analyzing specific kinetic and kinematic parameters and especially cartilage stresses and strains for large patient cohorts. Even though different amount of joint structures was included in the models and motion was different, in all models the obtained values were within the range of values reported in the literature. Here we only present the results for varus-valgus rotation (Fig. 5a), joint reaction force (Fig. 5b) and contact pressure (Fig. 5c).
The internal-external rotation depended on the method used to prescribe its motion. For the moment driven models, the moment scaling (i.e. the assumed absorption of forces by muscles) was based on a previous study 24 and it may influence the results 19,63 . For the rotation driven models, the contributions of the muscles and other joint structures are already considered indirectly in the motion input. However, measurement uncertainties from motion analysis, such as skin marker movement and cross-talk, may affect the accuracy of the measurement [38][39][40] . For these reasons, we cannot conclude whether moment or rotation is better to drive the FE model, especially as both methods produced results that are within the range of values presented in the literature (Supplementary Material and Fig. 5). For varus-valgus rotation, all models behaved almost identically. This was due to the fact that all models had the same ligament prestrain and stiffness, and ligaments have been suggested to have a significant role in varus-valgus rotation [64][65][66] .
It was not a surprise that the inclusion of patella and quadriceps force (Models A and B) increased both distal-proximal and anterior-posterior translational forces. The greatest differences between the models can be seen between 10% and 30% of the stance, which coincides with the first peak of the quadriceps force. However, by just adding the anterior-posterior component to the anterior-posterior translational force (Model D), almost these same joint reaction forces could be obtained. This could offer a simple and fast solution for modeling to account for the contribution of patella. A better solution, though more time-consuming, might be to estimate  the patello-femoral joint force, using for example musculoskeletal modeling, and implement it directly into the femoral reference point. On the other hand, in terms of medial-lateral translations and joint reaction forces there were only minor differences between the models. Even though the axial (distal-proximal) joint reaction force was at maximum ~0.6 BW different between Models A and C, the average maximum principal stress, logarithmic strain and pore pressure of cartilage at the lateral tibial contact area were similar in all models, independent of the motion implementation. In addition, differences in these parameters between the models were observed in the medial joint compartment only at around 10-30% of the stance, while during the last ~70% of the stance the differences between the models were nearly negligible. The differences were mainly caused by the different internal-external rotation. This leads to the conclusion that the internal-external rotation determines more the distribution of stresses and strains between the lateral and medial joint compartments, and less the magnitude.
Generation and simulation of the complex models required ~10 times more time than the simpler models. Adding the patella, tendons and patello-femoral ligaments increased mainly the time needed for the FE models to converge. This is because the added contact interactions, constraints and boundary conditions forced us to refine meshes several times and to find optimal convergence parameters, introducing additional work. It is also worth emphasizing that the simplest model (Model C) can be generated without musculoskeletal modeling, speeding up the model generation. Obviously, the time estimates presented in Table 2 apply only to the presented model. However, the general trend in time needed to generate and simulate this kind of models would remain the same.
This study has a few limitations. First, only one subject was included in the study. However, this was a methodological study showing that it is possible to obtain similar results from kinetic and kinetic-kinematic finite element models. In the future more patients will be added, however the main conclusions from the study should not change. A second limitation is the force implementation. Here the total joint force was produced primarily by the translational forces and ligament forces, due to their pre-strains and stiffnesses, together with quadriceps forces pulling the patella in the most complex model, similarly as has been done before 59,67 . Other muscle groups, such as hamstrings, contribute also to the total joint contact force 68 . However, for the kinetic-kinematic driven models (Models B, C and D), addition of other muscle groups should not influence the results, as their contribution was already considered in the motion capture measurement 41 . On the other hand, scaling of the measured moments in the kinetic driven model means that majority of these forces are absorbed by muscles and the final, scaled values reflect the moments caused by passive structures 24,63 . Therefore, muscle contributions were indirectly considered in the models by either implemented rotations in Models B, C and D or by moment scaling in Model A. These implementations of loading are identical to earlier studies 24,41,55 . However, we acknowledge that musculoskeletal modeling could be used to estimate all muscle forces and the FE model could include all these muscles 19,25 . This would, however, increase the complexity of the model further away from a clinical application and it is unclear if these would produce any different results.
Lack of direct experimental validation is also a limitation of the study. Real experimental validation of contact forces or pressures is possible with patients having instrumented knee replacements [69][70][71] . This is not possible in normal, healthy patients, since this would entail unnecessary surgical interventions. However, as indicated above, we compared our FE model results against literature reported values [72][73][74][75][76] , including experimental studies, to ensure the model outputs are within physiologically relevant ranges (Fig. 5). We would also like to mention that Model A is based on Halonen et al. 24 , in which the rotations obtained from the FE model matched quite swell with the experimental values. One can also compare the results of the FE model against those obtained from musculoskeletal modeling 77,78 . Obviously, this is not optimal comparison since both methods are computational. Another verification method is to compare the FE model results against experimental follow-up information of patients, if available. For example, T 1ρ or T 2 relaxation times may reveal degenerative changes in the knee joint cartilage [79][80][81] . This could be correlated with stresses and strains obtained from the FE models, similar to Liukkonen et al. 8 If we assume that the most complex model (Model A) is also the most realistic, then from the clinical point of view, it would be important to recognize how much a simple and fast model can differentiate from that model. Here, maximum values of the maximum principal stresses differed by ~1 MPa between Models A and C in the medial compartment, while the difference in the lateral compartment was close to zero. For the biomechanical analysis to predict joint and cartilage failure, 1 MPa is really small, especially as cartilage failure depends on many factors such as joint location, age and OA [82][83][84] . In the literature, cartilage stresses that lead to OA show a large range from 9 to 36 MPa, with fibrillar crosslink damage from 9 to 15 MPa and proteoglycan loss and nonvisual collagen damage from 11 to 36 MPa [85][86][87] . Based on these studies, 1 MPa uncertainty in the analysis makes a difference only if the result is close to this threshold value, i.e. 8-10 MPa.
In clinical practice, FE models could be used to identify areas susceptible to OA development and potentially indicate interventions to avoid or delay OA 6,[16][17][18]88 . Recent studies used FE modelling to simulate progressive collagen degeneration and/or proteoglycan loss in articular cartilage 8,11,15,88 . When applying these kind of models in the knee joint 8,9 , they would be particularly useful when studying OA development due to biomechanical alterations in the knee, such as due to ACL injury and reconstruction. However, most of the above-mentioned studies use FE models with complex material properties, which are extremely time-consuming, especially in a clinical setting. Simpler and faster FE models, as presented here, coupled with recent developments in semi-automatic segmentation methods [89][90][91][92][93][94][95] , might provide a better option.
In conclusion, the results showed that the kinetic-kinematic driven models with a simple geometry, including only cartilage, menisci, and ligaments, can produce mostly similar cartilage responses as the more complicated FE models, which include more tissue structures of the joint. They are also much faster and easier to implement and might therefore be more applicable for clinical practice or in large cohort studies, where the cartilage responses and possible failures are modeled.

Data Availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.