Using Design of Experiments Methods for Assessing Peak Contact Pressure to Material Properties of Soft Tissue in Human Knee

Contact pressure in the knee joint is a key element in the mechanisms of knee pain and osteoarthritis. Assessing the contact pressure in tibiofemoral joint is a challenging mechanical problem due to uncertainty in material properties. In this study, a sensitivity analysis of tibiofemoral peak contact pressure to the material properties of the soft tissue was carried out through fractional factorial and Box-Behnken designs. The cartilage was modeled as linear elastic material, and in addition to its elastic modulus, interaction effects of soft tissue material properties were added compared to previous research. The results indicated that elastic modulus of the cartilage is the most effective factor. Interaction effects of axial/radial modulus with elastic modulus of cartilage, circumferential and axial/radial moduli of meniscus were other influential factors. Furthermore this study showed how design of experiment methods can help designers to reduce the number of finite element analyses and to better interpret the results.


Introduction
Knee joint contact pressure is of critical importance in the mechanisms of knee pain and osteoarthritis [1,2]. Computational models and finite element analyses (FEA) have been utilized to study contact characteristics of normal and injured knees, as well as total knee replacements (TKR) [3][4][5][6][7][8].
The purpose of these studies was to determine peak contact pressure in order to predict either tissue degradation of the knee or wear of ultra-high molecular weight polyethylene (UHMWPE) in TKR. Some biomechanical factors, such as material properties and geometries of tissues [9,10], and knee kinematic [11] can affect the contact behavior of the knee and consequently the design of TKR. Impacts of horn attachments stiffness and meniscal material properties on tibiofemoral contact pressure using "semiautomatic" optimization method were investigated by Haut Donahue et al. [9], who set tolerances on the variables to restore the contact pressure to within a specified error. The authors, however, performed more than 60 analyses to determine whether an individual factor is of importance. Meanwhile, interaction effects between different factors were not considered in their study. In order to better interpret the effects of variations in the material properties of soft tissue, a powerful statistical approach is required to design computational experiments.
Design of experiments (DOE) is a formal mathematical method that helps to solve complicated problems and to save time and resources (cost) by reducing the number of required experiments (runs) while obtaining all the necessary information. However, reducing runs associate with decrease in resolution. Usually in an experiment, one or more factors are deliberately changed in order to observe the effect of these changes on one or more response variables. The statistical design of experiments is an efficient method for planning experiments so that it can be analyzed to yield valid and objective conclusions that can be obtained for a given amount of experimental efforts. Recently, there has been an increasing interest in use of DOE for sensitivity analysis based on FEA in biomedical applications [11][12][13][14][15]. So far, this method has only been applied in few studies related to the human knee joint [11,14,15]. Yao et al. [15] focused on the medial compartment of the knee and investigated the sensitivities of medial meniscal motion and deformation to material properties of soft tissues. They used Taguchi approach and central composite design to fit the finite element model (FEM) to the experimental data in the anterior cruciate ligament-deficient knee. Furthermore, Julkunen et al. [14] used a three-level fractional factorial design in combination with compositionbased finite element model to determine the effect of different cartilage constituents on the mechanical response of the tissue. Due to the uncertainty in material properties [16], finite element analysis of the tibiofemoral joint becomes a very challenging mechanical problem. Therefore, the aims of this paper are to explore the most important parameters related to the material properties of meniscus and cartilage affecting the tibiofemoral joint contact pressure and to make a regression model based on main interaction and quadratic effects of variables to understand how they influence and minimize the error of FEA output. In this regard, fractional factorial design was applied in screening step and Box-Behnken method was used in response to surface method and optimization process.

Creation of Finite Element
Analysis. Geometries of bony structures and soft tissues were taken from a healthy human knee of a 24-year-old man. Solid models of the femur and tibia and geometries of soft tissues, including articular cartilages and menisci, were developed from the magnetic resonance images (MRI). Each image was taken at 3.2 mm interval in a sagittal plane. The obtained data, subsequently, was used to create a three dimensional computer aided design (3D CAD) model in order to import into ABAQUS 6.8 software (Dassault Systèmes Simulia Corp., Providence, RI, USA). The model consisted of two bony structures (femur and tibia), both the femoral and tibial articular cartilages, and both the medial and lateral menisci. Figure 1 shows the generated 3D model in details. The model did not include ligaments. The finite element mesh generation was performed leading to 41709 linear 4-noded tetrahedron elements for articular cartilage and menisci (25293 for femoral cartilage, 9130 for tibial cartilage, 3866 for medial meniscus, and 3420 for lateral meniscus). Contact was defined between the femoral cartilage and meniscus, the meniscus and tibial cartilage, and femoral cartilage and tibial cartilage for both lateral and medial compartments, resulting in six contactsurface pairs. Completely general contact condition involving small sliding of pairs was applied on the model and all contact surfaces were modeled as frictionless. The cartilage in the knee is a complex structure, composed mainly of networks of collagen fibrils that embed water and a nonfibrillar matrix. The cartilage is known to be inhomogeneous and anisotropic material, but considering that the loading time of interest is related to a single leg stance and that the viscoelastic time constant for cartilage is approximately 1500 seconds from biphasic theory [3,9], the elastic solution does not diverge from the biphasic solution [17]. The cartilage, therefore, was assumed to behave as a homogeneous linearly isotropic elastic material for contact pressure computations, similar to the previous studies [18,19]. The meniscus, also, has similar structure to that of cartilage and it is also known to be inhomogeneous and anisotropic material, but various material property definitions can be found in the literature for this component [20][21][22]. Furthermore, the meniscus has a time constant, as large as 3300 seconds [9], and can also be considered as an elastic material for compression of the joint during the short loading times (single leg stance). In this study, the menisci were treated as linearly elastic, transversely isotropic material to represent the circumferential fiber arrangement. Femur and tibia were represented by rigid bodies because this is time efficient in a nonlinear analysis and accurate due to their much larger stiffness compared to that of soft tissues. Meanwhile, the previous study [3] confirmed that this simplification has no substantive effect on contact variables. Horn attachments, in the current model, were defined by 10 linear springs. For boundary conditions, the tibia was fully constrained and femur was constrained from rotation and free to translate in anterior-posterior, medial-lateral, and inferior-superior axes.
For validation of the model, static loads equivalent to 0, 500, 734, 800, 1000, 1500, 2000, and 2500 N were applied on the model at 0 ∘ flexion angle in order to compare with previously reported measurements and predictions [3,[29][30][31]. In this regard, the initial cartilage elastic modulus and Poisson's ratio were considered as 15 MPa and 0.475, respectively [25], and for menisci the primary values were moduli of 20 MPa in axial/radial directions and 140 MPa in circumferential direction. The values used for in-plane and out-of-plane Poisson's ratios and shear modulus were 0.2, 0.3, and 50 MPa, respectively [20,21,23,24]. The stiffness of horn attachments was considered 200 N/mm, which resulted in 2000 N/mm total stiffness. Figure 2 shows the finite element representation of the joint.

Verifying the Results of FEA.
The results of peak contact pressure for different magnitudes of force are shown in Figure 3. The applied force is transferred through the femurmeniscus, femur-tibia, and meniscus-tibia at the contact regions. The stresses were computed and it was seen that the total stress multiples by area equilibrate the total load in the knee joint. The predicted reaction forces at each loading condition, also, were in equilibrium with the applied load. Although, the finite element solution may have satisfied the equilibrium, indicating that the finite element solution was accurate to some extent, confidence in the validity of the model itself were obtained by comparing the computed values of the peak contact pressure with the previously reported measurements and predictions. Among the various researches that have measured the peak contact pressure on the tibial plateau [3,4,[8][9][10][29][30][31], studies of Brown and Shaw [31], Ahmed et al. [29], Fukubayashi and Kurosawa [30] and Donahue et al. [3] were chosen because they used a load application system with various compressive loads (734, 800, 1000, 1500, and 2500 N) at 0 flexion angle and computed the peak contact pressure on the tibial cartilage. It can be seen that, the results of present study fall well within the ranges provided by the literature. Hence the present results are verified.

Design of Experiments.
DOE starts with determining the objectives of an experiment. These objectives are as follows: comparative, screening, and modeling [32][33][34]. Objective of comparative designs is to find a suitable method for an initial comparison. Screening designs identify which factors are important and help to screen out unimportant factors. Response surface modeling seeks for one or more of the following objectives: hit a target, maximize or minimize a response or make it robust.
In this research, seven factors including axial/radial and circumferential elastic moduli of meniscus ( 2,3 and 1 ),   stiffness of meniscus horn attachment ( ), in-plane and outof-plane Poisson's ratios ( 12 , 23 ), shear modulus ( 12 ), and elastic modulus of cartilage ( ) were considered as initial variables for sensitivity analysis. Due to the large number of factors and levels, in the first step fractional factorial design was applied to screen out less significant factors. It is useful and efficient when full factorial design becomes unpractical [35]. Two levels for each factor impose (2 7 = 128) treatment combinations for full factorial design, but the 1/8 fractional factorial design suggested 16, so it can be used as a rational way for choosing the treatment combinations of experiments. The 1/8 fractional design corresponds to resolution IV in which the main effects are not confounded with two-way interactions. However, a limitation of fractional factorial design is the use of only two levels for each factor and the responses are assumed to be approximately linear over the range of the factor levels chosen. More detailed discussions can be found in Montgomery [36]. In the next step, after  screening out the less significant factors, Box-Behnken design was applied to do more investigations and hit the value of experiment in the FE model. The Box-Behnken is a good design in response surface methodology due to estimation of the parameters in the quadratic model. Furthermore, it is slightly more efficient than the central composite design [37], which was used by Yao et al. [15] in the sensitivity analysis of the knee joint. In this study, sensitivity analysis was performed under 2500 (N) static load at full extension and the reference value of the peak contact pressure was taken from the experimental work of Brown and Shaw [31], which was equivalent to 6.5 (MPa). At this load, the optimum values of parameters were obtained by the estimated model based on Box-Behnken design. The predicted optimum values were subsequently tested for other applied loads at 0 flexion angle. The study of Brown and Shaw [31] was chosen because it measured the peak contact pressure on tibial cartilage under the same loading condition (static load of 2500 N at 0 flexion angle). The considered factors and their investigated ranges based on the literature are demonstrated in Table 1. The combination of parameters was generated and analyzed using Minitab software [38], and the design generators were E = ABC, F = BCD, and G = ACD. Table 2 shows the treatment combinations and the results of peak contact pressure, according to fractional factorial design. The ranges in this step were chosen according to the most prevalent values used in the literature. For example, however the given range for elastic modulus of cartilage is 5-20 MPa, the majority of studies have considered values equal or more than 12 MPa [3,4,8,9,11,39,40].
The other parameters are not significant at the 5% level. It can be seen that, however, other factors including 1 and are not significant, interaction of 1 and is significant. Interaction is the variation among the differences between means for different levels of one factor over different levels of the other factor; and since in resolution IV designs, two-factor interaction effects may be confounded with other two-factor interactions, screening out of factors should be done carefully, because confounding pattern makes it difficult to determine which factors are the most important ones. In this regard, Table 3 also shows alias structure up to order 3. The alias structure indicates which effects are confounded with each other. As shown in Table 3, effect of 1 * aliased to 2,3 * and 23 * 12 . Since factors 2,3 and have been significant, so probably the main reason for significance of 1 * is due to the interaction of these two significant factors ( 2,3 * ); but for more confidence, factors 1 and are kept for further investigation due to their higher effect (0.0126, 0.0131, resp.) comparing to 23 and 12 (0.0121, 0.0091, resp.). Moreover, this result is in agreement with research of Haut Donahue et al. [9] that showed the importance of 1 and for contact variables of the tibial plateau. Furthermore, the importance of axial/radial modulus ( 2,3 ) and circumferential moduli of meniscus ( 1 ), as a result of the present study, is consistent with the findings of Yao et al. [15], who revealed that the meniscal motion and deformation are most sensitive to the circumferential and radial/axial moduli of menisci.
According to screening experiments, the following results can be estimated. (1) Peak contact pressure for = 12 (Mpa) is, on average, 1.1049 (Mpa) more than that for = 20 (Mpa).
(2) Peak contact pressure at 2,3 = 15 (Mpa) is, on average, 0.0496 (Mpa) more than that at 2,3 = 50 (Mpa). obvious that the main effect of factor is much larger than the other factors and overshadows them. Interaction effect of 1 and is shown in Figure 6. From the results of the screening experiment, factors of , 1 , , and 2,3 were collected for more investigation by response surface method (RSM).

Response Surface Method .
In this section, the response surface method was used as a statistical design of experiment tool, in order to produce precise maps based on mathematical models leading to optimum performance [41]. The regression model was built in two phases. First started with linear model but due to the lack of linear fit, quadratic model was applied subsequently. For choosing the level of factors which screened out in the last section, less interactions with other factors were considered ( 23 = 0.2, 12 = 60 MPa, and 12 = 0.2). The Box-Behnken method was used at three levels of each factor and a single center point was considered because the FEM experiments include no actual experimental error; thus, duplication of center point was not necessary [42,43]. The results of simulation runs and Box-Behnken design in RSM are given in Table 4. In this step, in order to assess more levels of each factor, other ranges of variables were chosen. Results of ANOVA revealed that with 99% confidence, increase of one unit of ( value <0.01) will result in decrease of peak contact pressure by 0.13871 MPa.
Although, the adjusted 2 demonstrates that 95.8% of variation in peak contact pressure can be explained by According to the outputs of ANOVA the following can be concluded.  confidence ( * * * means value <0.01), as −0.7265 * + 0.0182 * 2 , and with 95% confidence ( * * means value <0.05) as 0.0001 * 2,3 , if the effects of factor 2,3 are held constant.
(2) 2,3 affects the peak contact pressure with 95% confidence ( value <0.05) as 0.0001 * * 2,3 , if the effects of factor are held constant. (3) There is no reason to believe the importance, of other factors, interactions and any other quadratic effects. Residuals of full quadratic regression model are shown in Figure 7(b). According to the results of RSM, surface plot of peak contact pressure versus 2,3 and is shown in Figure 8. It demonstrates the negative effect of on peak contact pressure. Contour plot of peak contact pressure in Figure 9 shows that 2,3 has more effect on contact pressure than 1 . Figure 10 represents that the peak contact pressure does not change when is greater than 4000 (N/mm).

Optimization.
The optimal region to run a process is typically determined after a sequence of experiments and developing empirical models. From a mathematical viewpoint, the objective is to find the operating conditions that maximize, minimize, or close the system's response to the true one. Therefore, the goal of this section is minimizing difference between estimated quadratic model obtained in the last section and experimental data of Brown and Shaw [31]. By considering 99% confidence, the estimated quadratic model will only include factor . Figure 11 shows the behavior of peak contact pressure with respect to . In regression analysis, usually developing a model includes the fewest numbers of explanatory variables which permit an adequate interpretation: Min = PCP − 6.5 = (13.4733 − 0.7265 * + 0.0182 * 2 ) − 6.5 The value of was obtained to be 16.059 (Mpa) by solving the above nonlinear problem using generalized reduced gradient (GRG) algorithm on Microsoft Excel software. Furthermore, according to Figure 11, it is obvious that the minimum contact pressure is at = 20 (MPa).   Further analyses were carried out after optimizing the and removing its strong shadow on the other factors. Table 5 shows the Box-Behnken design with one center point, three factors, and peak contact pressure. In this stage, wider ranges of parameters were considered to investigate the maximum effects of factors.
Estimated full quadratic regression coefficients for peak contact pressure by considering optimized value of , factors According to the outputs of ANOVA, it can be concluded that with 99% confidence, 2,3 affects peak contact pressure as 0.00099 * 2,3 − 0.00001 * 2 2,3 , 1 as 0.00019 * 1 and there is no reason to believe the importance of factor , other interactions, and any other quadratic effect. Figure 12 shows normality of residual in the above estimated regression model The optimum values for 1 and 2,3 were 100 and 20 MPa, respectively. The suitable value for can be considered ≥2000 N/mm. This value supports the idea that meniscal replacement surgery should attach the horns through a technique providing high stiffness. Using the values obtained from the optimization process, finally the results of the proposed model were tested by running another two simulations with two compressive loads of 1000 N [30] and 3000 N [31] at 0 degrees of flexion. The errors between FEA and experimental data of peak contact pressure decreased to be less than 10% and 5% for 1000 and 3000 N, respectively.
It should be pointed out that in the FEM, some anatomical geometries are missing or simplified depending on the complexity of the problem. Ligaments were not included in our FE model. The posterior cruciate ligament (PCL) and lateral collateral ligament (LCL) are slack under axial compressive loading at 0 flexion angle [8], but the anterior cruciate ligament (ACL) and the medial collateral ligament (MCL) both contribute in the axial compression experienced by the joint. Under no external load, the joint is primarily compressed due to the prestress in these two ligaments and the axial compression sustained by the joint is, thus, greater than the applied external load. Therefore, the influence of missing these two ligaments in this FEA is that the results of peak contact pressure only correspond to the external load. Future studies will consider including ACL and MCL. However, according to Haut Donahue et al. [9], contact characteristics are not sensitive to the nonlinear material properties of 6.55 6.555 6.56 6.565 6.57 6.575 MCL during axial compression. Furthermore, cartilage can be assumed to have linear elastic material property in contact analyses from the biphasic solution, but the huge influence of its elastic modulus on contact pressure might indicate the requirement for more precise material model in this component. This is supported by the study of Julkunen et al. [14] which showed that the mechanical responses of the cartilage under different loading conditions are dependent on tissue composition and structure. Therefore, future investigations will focus on the effect of anisotropic nonlinear behavior of cartilage on contact outputs. Moreover, it will be interesting to investigate the sensitivity of contact pressure to material properties under different degrees of flexion, which was not considered in this research.

Conclusion
A sensitivity analysis of tibiofemoral peak contact pressure to the material properties of soft tissue was performed and design of experiments methods was used to reduce the number of program runs and to minimize the contact pressure error. The present study evaluated the effect of cartilage elastic modulus and interaction effects of the parameters in addition to previous research. It was demonstrated that elastic modulus of the cartilage is the most influential factor. Another important finding was that after cartilage elastic modulus, interaction of axial/radial modulus with elastic modulus of cartilage, circumferential and axial/radial moduli of meniscus are significant factors. The importance of circumferential and axial/radial moduli of meniscus as a result of this study is in agreement with the past predictions. Furthermore, this research demonstrated the complex relations between material properties of tissue and contact pressure of tibiofemoral joint. The result of sensitivity analyses can be used as a guideline for experimental efforts intended at determining material properties of soft tissue, because estimating the most sensitive parameters should be done precisely. However, this analysis is only valid under full extension loading mode and with elastic assumptions of soft tissues. Further biomaterial studies may reveal more factors or more realistic form of material properties of human tissues.
However, more investigation in this regard based on DOE techniques will provide a remarkably versatile strategy for analysis of knee joint biomechanics and help researchers with faster and more reliable analysis.