Kriging Surrogate Model for Resonance Frequency Analysis of Dental Implants by a Latin Hypercube-Based Finite Element Method

The dental implantation in clinical operations often encounters difficulties and challenges of failure in osseointegration, bone formulation, and remodeling. The resonance frequency (RF) can effectively describe the stability of the implant in physical experiments or numerical simulations. However, the exact relationship between the design variables of dental implants and RF of the system is correlated, complicated, and dependent. In this study, an appropriate mathematical model is proposed to evaluate and predict the implant stability and performance. The model has merits not only in the prediction reliability and accuracy but also in the compatibility and flexibility, in both experimental data and numerical simulation results. The Kriging surrogate model is proposed to present the numerical relationship between RF and material parameters of dental implants. The Latin Hypercube (LH) sampling method as a competent and sophisticated method is applied and combined with the finite element method (FEM). The methods developed in this paper provide helpful guidance for designers and researchers in the implantation design and surgical plans.


Introduction
The dental implant as a predictable and reliable treatment has been widely applied in the rehabilitation of edentulous patients [1]. The implantation fails sometimes caused by some complicated reasons in oral environment [2]. Fortunately the stability of the implants can be introduced to successfully predict such a failure in most cases by numerical computing or experimental testing [3]. Usually, the stability of implantation can be concluded in two categories: the primary and the secondary stability. The former is largely associated with osseointegration, while the secondary stability is highly corresponding with the bone formulation and remodeling in the process of healing [4][5][6][7]. Bone material quality, geometry characteristics of implants, and cortical bone can affect the primary stability [8][9][10][11][12]. The secondary stability obtains from the bone apposition surrounding the interface of implants and bone [13,14]. The quantification of implant stability offers helpful information to keep reliability in the individual treatment [15,16].
The histological examination is one of the traditional invasive approaches. Noninvasive methods are required for the observation and measurement of implant stability [17,18]. The Periotest and the Osstell Mentor system test are two typical noninvasive methods for the measurement of implant stability in diagnosis. In the Periotest, the interfacial damping characteristics between the implant and the surrounding tissue are evaluated [19,20]. However, the Periotest is often criticized for its lack of prognostic accuracy and poor sensitivity. On the contrary, the Osstell Mentor system is based on resonance frequency analysis (RFA), which appears more competent of assessing the implantation stability [21][22][23][24][25].
The RFA was firstly introduced for dental applications in 1996 and then developed in both experimental and computational aspects [26][27][28][29][30][31]. It provides an effective way to study the relationship between the stiffness of the implant-bone interface and its surrounding local structures during the healing process [32][33][34]. In the numerical computation aspect, RFA can be powerfully implemented by a finite element method (FEM). FEM is capable of simulating puzzling geometry characteristics, material properties, and also boundary conditions, which are often tough to deal with in the laboratory [35][36][37]. Furthermore, it is feasible to allow any independent control of parameters in the implant system by the FEM. Then, performing systematic evaluation for each parameter corresponding to implant stability becomes convenient.
Combining the Latin Hypercube (LH) sampling method with the FEM develops a competent and appropriate stochastic finite element method. The widespread popularity of LH has led to the technical development in various fields, such as the improvement of space filling [38,39], optimization of projective properties [40][41][42], minimization of least square error, maximization of entropy [43,44], and reducing spurious correlations [45]. Meanwhile, LH has been deeply explored in probabilistic analysis, ranging from the assessment of reliability [46][47][48][49] to coefficient evaluation for polynomial chaos and merging with the surrogate models [50,51]. Therefore, the combination of LH and the FEM is a promising method for the RFA of the dental implantation research.
Since the exact relationship between the design variables of dental implants and RF of the system is correlated, complicated, and dependent, traditional regression calculation is difficult to reach a satisfied accuracy. The Kriging model is an interpolation method which finds its roots in geostatistics [52]. As one of the most promising spatial correlation models, the Kriging model is more flexible than the regression model and not as complicated and time-consuming as other metamodels [53]. Recently, there is an increasing interest in applying the Kriging model in the industry, mechanical engineering, and related fields [54][55][56][57][58]. The popularity of Kriging consists in that Kriging is an accurate interpolating approximation model [59]. Kriging model is attractive for its interpolating characteristic, providing predictions with the same values as the observations and reducing the time for the expensive analysis. When there are highly nonlinearity in a great number of factors, polynomial regression modeling becomes insufficient while the Kriging modeling is an alternative choice in spite of the added complexity [60].
The goal of this study is to propose an appropriate mathematical model to evaluate and predict the implant stability and performance. The model is explored to precisely describe the exact relationship between RF and design parameters of dental implants. LH-FEM can reduce the cost of physical experiments, which are expensive and time-consuming. The proposed model in this study has the advantages not only in the prediction reliability and accuracy but also in the compatibility and flexibility in experimental data and numerical simulation results.
This paper consists of four parts. Section 1 provides a brief overview of the Kriging surrogate model and resonance frequency analysis. Parameters of dental implants are presented in Section 2. Besides, the implementation of LH-FEM is also performed in this section. In Section 3, results are demonstrated and discussed comprehensively. A short summary is concluded in the last section. The methods developed in this paper provide useful information and helpful guidance for implant designers and prosthetic professions in the process of implant design and corresponding surgical plans.

LH-FEM for Dental Implants
2.1. Parameters of Dental Implants. The geometrical parameters of dental implants are expressed in Figure 1 and Table 1. They include height, length, and width of dental implants, cortical bone, cancellous bone, and ceramic dent. There are 8 parameters (H1, H2, H3, W3, h1, h, w, and d) for dental implants and 2 parameters (W2 and H4) for ceramic dent. H, W, and W1 for the cortical bone are 26 mm, 16 mm,

H1
The height of the bottom part in dental implants 14.2 H2 The height of the middle part in dental implants 2 H3 The height of the top part in dental implants 5

H
The height of the cortical bone 26 h1 The height (in the front end) of the thread in dental implants 0.1 h The distance (in the front end) between the edges of the thread in dental implants 0.9

W
The width of whole FEM 16 W1 The thickness of the cortical bone 1.3 W2 The diameter of the ceramic crown 5.1 W3 The diameter of dental implants 3.1 w The width of teeth in the thread part of dental implants 0.4

H4
The height of the ceramic crown 2 Besides, d represents the degree of the included angle in the front end of thread and is settled as50 ∘ . and 3.1 mm, respectively. With the development of material science, stainless steel, titanium, gold, and even fiber can reach a biomedical grade and can be applied in dental implants. For material properties, Young's modulus, Poisson's ratio, and physical density of dental implants (E1, R1, and D1), ceramic crown (E2, R2, and D2), cortical bone (E3, R3, and D3), and cancellous bone (E4, R4, and D4) are 12 input variables in the finite element model. In order to perform the Latin Hypercube-based finite element method, the specific intervals of parameters corresponding with material properties are given in Table 2.
In terms of the material properties, on the one hand, the work in this study provides the appropriate and wide interval ranges for the material parameters instead of the certain and specific values. The specific values of material parameters are included in the interval ranges. On the other hand, the material properties are analyzed and discussed in the aspects of stiffness and mass matrix, which are more comprehensive to analyze the effects of material properties in RFs of dental implants.
According to the FEM, (1) all materials are homogeneous, isotropic, and linear elastic; (2) the different components exhibit different physical properties; the Young's modulus and Poisson's ratio are given as input variables according to common values in the certain range; (3) perfect bonding is involved between implant-abutment-screw, abutment-crown, and implant-bone. In the contact pair, the objective surface and contact surface are supposed to be bonded. Perfect bonding makes the computation process avoid the problem of convergence, speeds up the calculation, and is proper to have the analysis in large deformation and nonlinear problems; (4) there are no flaws in any components; and (5) for the boundary condition in the FEM of the dental implant system, three surfaces (in the bottom and left and right sides) of the cortical bone surface are chosen and completely fixed; the displacement in X, Y, and Z directions of these three surfaces is constrained to be zero.
Besides, in the previous work of Chu et al. in reference [49], it has been proven that a sufficient number of samples can guarantee a satisfied accuracy of LH and MCS. Enlarging the sampling pool can effectively improve the result accuracy when the number of samples is small; however, when the amount of samples reaches a certain number, the improvement is not evident. Therefore, the number of samples for each variable is supposed to be 500 to obtain an appropriate accuracy while reducing computation costs.

Latin Hypercube-Based Finite Element Method.
A finite element model of the dental implant system is created as shown in Figure 2 by ANSYS (Mechanical Parameter Design Language, Version 14.5, USA). Solid 285 is the chosen finite element, which is a tetrahedron solid element with 4 nodes in each element, and each node has 3 degrees of freedom (displacement in X, Y, and Z directions). The  tetrahedron shape of Solid 285 is flexible and convenient to mesh nonlinear and complicated geometry components, such as thread of dental implants. There are a total of 133961 elements and 22168 nodes. The thread components of dental implants have been fine meshed. By performing a finite element (FE) procedure, the RF can be obtained by solving the eigenvalue problem in the govern equation through the block Lanczos method.
The study of the relationship between the various corresponding parameters of the dental implant system and the values of RF requires fabrication of a huge sample set for experiments, which is expensive and time-consuming. LH as an advanced Monte Carlo method is especially efficient, which divides the sample space into a number of subspaces, then samples from subspaces, thereby perfectly avoiding sample clustering or variation in the boundary [49,50]. Combing the traditional finite element model of dental implants with LH is a sophisticated method to overcome the disadvantages of physical experiments.
The input variables of the finite element model were assigned using LH sampling to effectively avoid sample repetition or clustering. For each variable, there are five hundred randomly samples, which are uniformly distributed in the specific interval range as shown in Table 2. The stochastic sampling process of the input variable creates a reliable database for FEM computation. The flowchart of LH-FEM in Figure 3 presents the programing process. The process marked by blue color represents the deterministic finite element model of the dental implant system. The convergence and accuracy of the deterministic finite element model for dental implants are verified before the next steps. After the validation of the deterministic finite element model, the original codes are assigned to LH-FEM. The loop of LH-FEM does not stop until all of the samples are computed by the finite element model. Then, the results of RF are captured and stored in the output database, as shown in the right side of the flowchart in Figure 3 which are marked by red color.
In order to find the availability of LH-FEM, a parameter sampling record (E1) of input variables is presented in Figure 4. The difference between samples is not regular because of the stochastic sampling process, which makes the sampling set including different situations as far as possible. The probability distribution in statistical mathematics is following a uniform distribution as in Figure 4. The LH sampling method is successfully implemented in the whole sampling spaces.

Results and Discussion
3.1. Results of the Deterministic Finite Element Model. RFA (resonance frequency analysis) as a noninvasive and nondestructive technique is proposed to assess the implant stability. A deterministic finite element model of dental implants is performed to evaluate the feasibility of this method. The resonant displacements of the first five order resonance frequencies are plotted in Figure 5, for the X, Y, and Z directions and vector sums, respectively. The axes in Figure 5 are the same as that in Figure 2. The contour results allow visualizing the displacement and mode shape of the dental implant system. The contours present that the displacement happened in the direction of Y axis in RFA is a more dominant vector that influences the whole system in the first-order RFA than others. While in the second-order RFA, the more important displacement is observed in the direction of X axis. In addition, the displacement vector sums in the third-order RFA are more approximated to the displacement in the direction of Z axis. For the fourth-and fifthorder RFA, the displacements in vibration modes are more complicated and the natural frequencies are larger than the lower-order vibration modes. Besides, the largest deformation in the firstand secondorder vibration modes happens in the top of the system, which well agreed with the results in the work of Li et al. [36].
In order to be more convenient, the results of dental implants are demonstrated independently. Figure 6 presents the displacement vector sum and Von Mises stress of dental implants in the first five order RFA. The associated displacement corresponding to the first bending modes of a cantilever beam has been experimentally observed in the literature [61]. The results of the deterministic finite element model are consistent with those in experiments of the literature. The accordance can be confirmed in both the displacement vector sums and Von Mises stress. The failures and risks most likely occur at the beginning of the threads and the junction of implant collar. This fact is in a good agreement with numerical and experimental investigation [37,62]. Therefore, the results in this study are validated. The finite element model of dental implants is feasible and appropriate for further studies.
In this study, not only the resonant frequencies of dental implants are calculated by the FEM but also the vibration modes of the dental implants are provided. As the direct measurement of the vibration performance of dental implants higher than the third mode is difficult, the displacements of dental implants under the vibration mode in this numerical study are an important supplement of the experimental measurements. Furthermore, the high mode vibration is an appropriate reference for the safety and reliability of the dental implants in the real operating environment.

Probabilistic
Results of LH-FEM. The results of LH-FEM are a large database, and the effective and useful information is extracted and presented in Table 3. The maximum, minimum, mean value, and variance of the first five order RFs are all shown in Table 3. The results of the first-order and second-order RFs are very close, which is because of the geometrical symmetry in the finite element model. However, the vibration modes of the first and second orders are different in Figure 5. In the first-order vibration mode, the dominant displacement happens in the direction of Y axis, while in the second-order vibration mode, the important displacement occurs in the direction of Z axis.
Different with the probability distribution of input variables, the output results of the first five order RFs for dental implants do not have uniform distributions as shown in Figures 7 and 8. From Figure 7, the results of firstand second-order RFs are more concentrated in the smaller interval than the third, fourth, and fifth orders, which means the difference in material property of dental implants causes a large deviation of RF in higher-order vibration modes. Besides, in Figures 7 and 8, the probability distribution and cumulative probability of the firstand second-order RFs are approximated, but when compared with the third-order RF, the difference is very evident. These characteristics of the probabilistic results provide helpful guidance in experimental designs and research.
Furthermore, the primary (first order) RF of LH-FEM is fluctuated in the interval range from 6100 to 10000 Hz approximately, and the average is 7800 Hz. In regard to RFA, tremendous substantial studies have been accomplished in experiments and clinical research. In the work of Glauser et al. [26], the implant stability based on the RFA of dental implants under an early functional loading is explored.
During a 12-month interval, the average RF ranges from 6100 Hz at the 1st month to approximately 6600 Hz at the 12th month. Besides, the results of the experiment in the tibia  of three groups of guinea pigs in a period of 4 weeks [34] also indicated that the average RF is roughly 5650 Hz amongst the three groups. Compared with these literature data, the results of LH-FEM are in a reasonable range. Therefore, the database of LH-FEM for RFA of dental implants is available and feasible. The work of this study provides the maximum, minimum, mean values, and the variances of the resonant frequencies of dental implants. The primary (first order) resonant frequencies of dental implants corresponding to the implant stability in the literatures fall within the interval range of LH-FEM, which strongly proves the feasibility of the proposed model in this study. Besides, the probability density distribution results of the first five order RFs for dental implants in Figures 7 and 8 provide important information corresponding to the safety and reliability of dental implants' stability.

Prediction of the Kriging Surrogate Model.
The resonance frequency is directly attributed to the stiffness matrix and mass matrix. The material stiffness and the stiffness of implant-bone interfaces and surrounding tissues have positive effects on RF, while RF usually has the inverse proportion to the mass matrix. However, the Young's modulus of the cancellous bone is crucially related with the physical density, and the increase of physical density of cancellous bones (mass matrix) can contribute to the improvement of the Young's modulus (stiffness matrix). On the other hand, cancellous bones have a much larger contact area to the dental implant and can hugely affect the whole dental implant system. Therefore, the relationship between RF and parameters of material properties in the dental implant system is complicated and correlated. It is difficult to be described by a simple linear interpolation or regression function.    Based on the reliable database of LH-FEM, a huge amount of samples for parameters of material properties in the dental implant system and corresponding RF is provided. By performing the Kriging surrogate model, the numerical relationship between RF and parameters of material property is created. The accuracy and convergence of the prediction by the Kriging surrogate model are demonstrated in Figure 9. The black spots in Figure 9 are prediction results of the Kriging surrogate model; the mesh surface is the results obtained     from LH-FEM. Satisfied level of accuracy is reached by the Kriging surrogate model. Figure 9 well confirms the prediction accuracy of the proposed Kriging surrogate model in this study. In order to have quantitative comparison, the prediction results of the Kriging surrogate model and the results in the literature [36] are presented in Figures 10 and 11. In clinical and experimental observation, RF increases and changes in the beginning of the first several months and reaches a certain steady state after a period. In the healing process, the Young's modulus of the cortical region and cancellous region increases. Therefore, the other input variables are settled as certain values in the Kriging surrogate model; the Young's modulus of cortical (E3) and cancellous (E4) bones are input variables.
From Figure 10, it can be found that in general, according to the increase of the Young's modulus cortical region, the first three order RFs are all amplified in the beginning period in the results of the literature [36] and the prediction of the Kriging surrogate model, where β is the ratio of the RF related to the specific Young's modulus with the steady RF in the final phase. Due to the database of LH-FEM, the first-order and the second-order prediction results are very close, while the prediction results of the Kriging surrogate model in the third RF are approximated with the results in the work of Glauser et al. [26], especially in the beginning and end parts of the result curve.
Furthermore, the prediction results of RF in the cancellous region by the Kriging surrogate model are also compared with the results in the literature [36] as shown in Figure 11. A good agreement is achieved when the Young's modulus in the cancellous region is large. However, the prediction results of the Kriging surrogate model are larger than the deterministic results [36] when the Young's modulus in the cancellous region is large. The reasons causing this phenomenon can be the original database of LH-FEM, the low sensitivity of the Kriging surrogate model for the cancellous region, or computational relative errors of the method applied in this study. On the other hand, the Kriging surrogate model provides continuous result prediction of RF for the dental implant system, which is more convenient, comprehensive, and time-saving for RFA than the traditional finite element method calculation and clinical test.
In addition, based on the accuracy and reliability of the Kriging surrogate model prediction, it is also useful in the analysis process of dental implant stability. However, the dental stability cannot be directly measured or simulated with an in vivo or vitro model under specific clinical situations. The Kriging surrogate model proposed in this study can play important roles and provide believable prediction results of RF corresponding with the dental stability. Besides, the Kriging surrogate model is compatible to combine the numerical simulation results with clinical and experimental results in the original database, which makes it a more sophisticated surrogate model for the dental implant system.
Through LH-FEM proposed in this paper, a series of meaningful results are obtained for RFA of dental implants. However, there are still some limitations caused by simplification and assumptions involving in the numerical simulation. Firstly, material properties of each component in the dental implant system are supposed to be homogenous and isotropic. Secondly, the damping effect is totally neglected during the RFA, while damping may have a certain extent influence on the accuracy of RF. Thirdly, the interfaces in the system are all treated as perfect bonding without considering the local specialty. In the real operation situation, the bone-implant interface is a dynamic living surface that evolves from a debonded interface to a bonded interface. Li et al. [60] applied the bonded model to calculate the RF of the bone-implant interface in a dental implant. The results were in a quantitative agreement with some experimental measurements, because the bone-implant interface was fully osseointegrated after several weeks [4,61]. Despite such simplification in the simulation process, the computational results in RFA are fairly reasonable, which could be useful and helpful for implant designers and prosthetic researchers. Furthermore, the prediction results of the Kriging surrogate model for RFA of dental implants by LH-FEM are reliable and believable.

Conclusion
In conclusion, the computational model proposed in this paper is a successful numerical tool for noninvasive RFA in dental implant research. The response surface method is a typical local model and can be written in polynomial series as follows [49]: Equation (A.1) and equation (A.2) are the local firstand second-order response surface model, respectively, where F β x is a deterministic regression model and β i is the corresponding regression coefficient.
In addition, a global surrogate model generally has global searching space. Kriging models fit a spatial correlation function as follows: where z x is assumed to have mean zero and covariance.
In the first-order linear Kriging surrogate model, the predictor can be expressed as follows [63]: where c = c x ∈ R m . The relative error iŝ In order to make sure the predictor is unbiased, F T c is equal to f x . The mean squared error of the predictor can be written as follows: The objective of the optimization program is to have the minimum value of φ; the constraint is The computation of the gradient of the constraint function with respect to c can be expressed as follows: Suppose λ = −λ/2σ 2 , the following system of equations is The solution can be obtained: where R −1 is the inverse matrix of the correlation matrix R.
The generalized least squares solution is Substitute β * into the predictor of the Kriging surrogate model, Besides, the corresponding maximum likelihood estimate of the variance is written as follows: If the errors are uncorrelated and have different variances, E e i e i = σ 2 i and E e i e j = 0 for i ≠ j. It is logical to find that R is the diagonal matrix, Besides, the weight matrix W is given as follows: Therefore, And, the below equations are satisfied: Replacing F and Y by weighted function, the results can be depicted as follows:

A.2. Resonance Frequency
Modal analysis is efficient to identify the RF of objects in the fields of engineering, industry, and medicine. Resonance happens when a structural system vibrates at its RF with a tendency to oscillate in higher amplitudes. Since the value of RF is associated with material stiffness, structural damping, physical density, and boundary condition, RF can be employed to symbolize a structural system. The calculation of RF of a system is essentially a generalized eigenvalue problem. The free vibration without damping is governed by the following equation: where K is the structure stiffness matrix and M is the structure mass matrix.
For the linear system, free vibrations will be harmonic of the form: In the above equation, ϕ i is eigenvector representing the mode shape of the ith natural frequency, ω i is the ith natural circular frequency, and t is time.
Thus, the free vibration can be presented as follows: This equality is satisfied if either ϕ i = 0 or −ω 2 i M + K equals to zero. Then, This is an eigenvalue problem which may be solved for up to n values of ω 2 and n eigenvectors ϕ i , where n is the number of degree of freedom.
The eigenvalue and eigenvector problem needs to be solved for mode-frequency and buckling analyses. It also has the form as follows: where λ i is an eigenvalue. The block Lanczos method uses the sparse direct solver to perform Lanczos iterations and extracts the requested eigenvalues. Rather than the natural circular frequencies ω , the natural frequency f is Therefore, the natural frequency f is not only corresponding with structure stiffness matrix K but also be affected by M structure mass matrix. It is an effective index to observe the changes or deviation of the structure geometry and material property.

A.3. Latin Hypercube Sampling Method
The Latin Hypercube sampling method is one kind of advanced Monte Carlo simulation (MCS). This approach divides the range of each variable into disjoint intervals of equal probability, and one value is randomly selected from each interval. It improves the stability of MCS and keeps satisfied accuracy and good convergence [49,50]. Consider a statistic system described by the function, Y = F X , X = X 1 , X 2 ,⋯,X n , A 26 where X is the random vector and represents the independent input random variables. F is the operator and performs computer simulation, such as finite element computation. Traditional MCS relies on simple random sampling, in which realization of X, denoted by x k , k = 1, … , N, where N is the amount of samples. In the sample space, where U i are the uniformly distributed samples on [0, 1] and P is the cumulative distribution function. Besides, Nataf or Rosenblatt transformations can be used to produce a set of uncorrelated random variables from correlated variables. Latin Hypercube sampling method divides the range of each vector component into disjoint subsets of equal probability. Samples of each vector component are captured from the respective subsets according to equation as follows: x j ki = P −1 X i U ij , A 28 where i = 1, … , n and j = 1, … , m, where n refers to the total number of vector components or dimensions of the vector and m is the number of the subset in a design. k is the subscript which denotes a specific sample. Besides, U ij are uniformly distributed samples on ξ j , ξ j ′ , By the Latin Hypercube sampling method, the ranges of all random input variables are all divided into intervals with equal probability, the efficiency of sampling is improved, and the disadvantage of clustering points is well overcome, as presented in Figure 12.
In Figure 12, the data are the examples to compare the efficiency of the Latin Hypercube sampling method and Monte Carlo method. Actually, it can be seen as the example of sampling space of parameters in the model of dental implants.

Data Availability
The original data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.