Finite Element Study of Periodontal Ligament Properties for a Maxillary Central Incisor and a Mandibular Second Molar Under Percussion Conditions

The quantitative percussion diagnostics (QPD) response of a mandibular second molar and a maxillary central incisor including their supporting ligament/bone structure was simulated using dynamic 3D finite element analysis (FEA). The focus of the work was on the role of the periodontal ligament (PDL) which acts as a damper in the dental structure and dissipates occlusal forces transmitted from the tooth surface to the surrounding bone. Several FEA models were developed to examine the effects of mechanical characteristics that have been reported for the PDL. Specifically, the effects of changing the PDL’s quasi-static elastic modulus and Rayleigh damping properties were predicted. The present FEA simulations indicate that the PDL can significantly reduce forces for both the incisor and the molar compared to when there is no PDL (i.e. ankylosed tooth) as long as the quasi-static elastic modulus of the PDL is among the lowest reported (~ 0.1 MPa). In addition, the FEA simulations for both the incisor and molar with this lower value of the PDL quasi-static elastic modulus are also in reasonably good agreement with experimental percussion data. A simple approximation for partitioning Rayleigh damping properties between the hard and soft tissues was also found to provide reasonable values of overall damping that are consistent with experimental data. The overall findings indicate that using a quasi-static elastic modulus of approximately 0.1 MPa for the PDL in combination with Rayleigh damping gives realistic predictions of the mechanical response of a tooth under QPD loading conditions.


Introduction
The structure of a natural tooth complex consists of enamel, dentin, ligament tissue and the surrounding bone. The periodontal ligament (PDL) serves as the connective tissue between a tooth root and the bone and has been the most difficult tissue to characterize mechanically. As a result, there is varied and sometimes conflicting information in the scientific literature concerning the mechanical properties of the PDL [1]. The PDL has a densely fibrous structure as well as cellular and vascular structures that provide soft tissue continuity between the hard tissues [2]. As such, the mechanical response of the periodontium to occlusal loading depends significantly on the behavior of the PDL [3]. Previous attempts at determining PDL mechanical properties [4][5][6][7] are briefly discussed in the following.

PDL Constitutive Models
Various constitutive models for the PDL have been developed [8][9][10] and used in FEA simulations [10][11][12][13]. It has long been known that the PDL exhibits non-linear behavior with respect to stress versus strain. However, many researchers have used an isotropic linear elastic law that does not take into account the viscoelastic properties of the PDL [13]. This simplified behavior has been assumed in a number of finite element works involving the PDL [14,15] even though it can result in a significant error for large strains or strain rates. A number of studies have investigated viscoelastic constitutive models for the PDL [11,16]. Others have considered a hyperelastic constitutive model [1,17,18]. Nishihira and coworkers [19] employed a multi-phase constitutive model and Liao et al. [20] considered the PDL as linear elastic but with a piecewise elastic modulus. Finally, other constitutive models have also been developed to describe hyperelastic damage [21], and compressible transversely isotropic hyperelastic behavior [8].
Hyperelastic constitutive models have demonstrated effectiveness in representing the non-linearity of the stress with respect to large strains above about 10% [1,3,17,18,22]. However, linear models have been shown to be adequate for lower strains [1]. Since diagnostic percussion loads are relatively low, it was assumed that PDL strains would be well below those that would justify a hyperelastic model. However, a viscous augmentation to the linear model still needs to be employed in order to capture the damping characteristics of the PDL under percussion loading. Accordingly, linear models augmented with Rayleigh damping were used in the present work.

PDL Thickness
Schroeder [2] has presented a comprehensive review of scientific reports on the thickness of the PDL. These reports have revealed that PDL thickness varies with location in the mouth (maxillary or mandibular), whether it is functioning or non-functioning, age, and other factors. Multiple researchers have reported that the thickness of the PDL is generally greater in function (~ 340 μm) than in non-function (~ 280 μm). By contrast, the minimum thickness was also observed at the tooth apex for an out of function PDL [2]. It was also found that the normal thickness of the PDL was reduced by 50-60% in the absence of an antagonist. Finally, a reduction of PDL thickness with respect to aging was also reported and confirmed by multiple researchers [2]. In general, clinical considerations have shown that the thickness of a healthy PDL generally ranges from 150 to 380 μm [5].

Quantitative Percussion Diagnostics (QPD)
We investigated the mechanical response of a maxillary central incisor and a mandibular second molar under conditions consistent with quantitative percussion diagnostics (QPD), a dynamic transient loading that has a duration of about 0.2 ms. QPD was devised to nondestructively assess the mechanical characteristics of teeth and dental implants. QPD is able to determine overall tooth mobility, and the mechanical integrity of the underlying bone interface for an implant [23][24][25]. It has also been used to detect cracks in teeth and determine their severity [26][27][28][29].
QPD involves lightly percussing the tooth while recording and analyzing the force-time response of the specimen. As opposed to cyclic loading, percussion loading only occurs during the few tenths of a millisecond that the percussion rod is in contact with the tooth. During that period, the QPD percussion rod administers a reproducible level of kinetic energy to the tooth that is governed by the rod's mass and initial velocity prior to impact. Thus, consistent initial kinetic energy is the controlled test condition by which measured mechanical responses to the resulting percussion are compared using QPD.
A primary objective of the present work was to determine mechanical properties for the PDL in an FEA model of a maxillary central incisor and a mandibular molar that lead to realistic simulations of measured force as a function of time under QPD conditions. Another main objective was to determine how much force should normally be reduced by the presence of the PDL. A series of models were performed with various properties for the PDL, and the resulting predicted behaviors were compared to experimental results as well as clinical information on the mechanical function of the PDL. Results for models where the PDL was absent, as in the case of ankylosis, were also compared with results for when the PDL is present to better understand how much it reduces forces on the teeth.

Methods
Physiologically accurate 3D models of a mandibular second molar and a maxillary central incisor were created in Solidworks (Dassault Systèmes) based on computer tomography (CT) data provided by eHuman, Inc. The CT data contained a point-cloud that traced surfaces of enamel, dentin and pulp chamber, as well as the interfaces between enamel and dentin. 3D CAD models were created by connecting the relevant points together and generating surfaces that passed through the said points and surfaces. The tooth models then include enamel and dentin components with a pulp chamber in the dentin. The present tooth geometries and periodontal ligament are shown in Fig. 1.
The solid models of these teeth were meshed for finite element analysis (FEA) in MSC Apex. The pre-processing, processing and post-processing steps were performed in MSC Marc/Mentat software. The finite element mesh that we developed for modeling a mandibular second molar is illustrated in Fig. 2. A percussion rod was also included in the present 3D models to simulate QPD accurately and for predicting the force measured by a sensor in the rod as a function of time. The present FEA meshes have a total number of about one million elements each. The location and angle of the percussion rod to each tooth was based on typical clinical use of the Periometer hand piece. For molars, the rod percussion is generally positioned at the buccal-mesial cusp since this is the easiest cusp on these teeth to access. This rationale was also used to apply similar percussion rod positioning for in vivo studies conducted previously [26][27][28][29].
We used second order isoparametric three-dimensional ten-node tetrahedron elements for the PDL, eight-node, isoparametric, arbitrary hexahedral elements for the percussion probe, and linear isoparametric three-dimensional tetrahedron elements for the remaining structures. Boundary conditions were defined to prevent free body motion in that the elements on the outer surfaces of the supporting bone were constrained. This constraint is reasonable when one considers the mass of the human mandible and maxilla compared to the relatively low forces induced by QPD. The percussion rod tip and at the tooth enamel were specified as touching contact bodies in MSC Marc [30] in to order implement the percussion impact. They are initially separated from each other prior to contact and then are in touching contact during the percussion response. The tooth enamel, dentin, PDL and bone were specified as "glued" contact bodies in MSC Marc [30].
Mesh size sensitivity of the models with various PDL thicknesses was studied. The results showed that increasing the mesh size causes the maximum force to increase slightly and the time duration of the response to decrease slightly. This is expected as reducing the number of elements (increasing the mesh size) in a model generally makes the model stiffer. Regardless, the curves remained bell-shaped which matches experimental results. The present mesh size was then selected to achieve a reasonable combination of accuracy and FEA run time ( Fig. 2).
At least two layers of elements were used in the PDL for each FEA model [7]. To validate the use of two layers, results were compared for two otherwise identical models that had either two or three layers of elements in the PDL. The number of elements in the model with three layers was approximately 50% greater than that in the model with two layers in the PDL. A comparison of the results for two models, one with two layers of elements in the PDL and one with three layers of elements in the PDL, is presented in Fig. 3. The difference in force between the two sets of results is less than about 4% over the duration of the percussion response. Thus, it was concluded that two layers of elements in the PDL are sufficient for reasonably accurate results.
The percussion rod was free to move in only its axial direction with a mass of 9 g and an initial velocity of approximately 60 mm/s, which are consistent with the operation of the Periometer hand piece [31]. The models were run with a time increment of 2 µs. A direct integration method was used to obtain the solution to the equations of motion for our models.
Proportional viscous damping is commonly employed to apply modal analysis of undamped systems to damped systems. A proportional damping model expresses the damping matrix as a linear combination of the mass and stiffness matrices, that is.
where D is the global damping matrix, M i is the mass matrix, and K i is the stiffness matrix, i is the mass matrix multiplier damping coefficient, and i is the stiffness matrix multiplier damping coefficient for the i th element [30]. This damping model is also known as 'Rayleigh damping'. Modes of Rayleigh damped systems preserve the simplicity of the real normal modes as in the undamped case [32]. Thus, i and i can be readily determined from a modal analysis of the undamped structure. In addition, MSC Marc used in  Finite element results for a central incisor for a PDL with two layers of elements and a PDL with three layers of elements for models that were otherwise identical. The PDL elastic modulus and thickness for both models were 0.1 MPa and 300 μm, respectively the present work employs tangent stiffness which avoids the possibility of any artificial Rayliegh damping [30,32]. Accordingly, Rayleigh damping was applied in the present models to simulate the viscoelastic behavior of the PDL. The damping matrix for the present models is defined as a linear combination of the mass and stiffness matrices of the system and damping coefficients are specified on an element-byelement basis.
Orthogonality and relations between the modal equations were used to calculate the and coefficients based on a modal damping factor and the undamped circular natural frequency [20,33]. Accordingly, α and are given by where j and j are the modal damping ratio and the undamped circular natural frequency for the j th mode, respectively.
The damping ratio is a dimensionless parameter associated with the decay in oscillations for a system subjected to an external force and disturbed from its static equilibrium position. Huang and coworkers [34] measured the damping ratios for a maxillary central incisor in vivo, which were determined to be in the range of 0.09 to 0.24. Following the procedure described by Liao and colleagues [20], we assigned these damping ratios to the first and 20th modes, respectively. Interpolation was used to estimate the damping ratios for the second to 19th modes as described in [20,33]. We obtained the first 20 undamped natural frequencies, i.e., i for i = 1 to 20 by performing a modal analysis for the maxillary central incisor model in MSC Marc. We assumed that no significant damping occurs in the tooth enamel and in the stainless steel percussion rod. In addition, the total strain energy of the entire structure and the damping energy dissipated in the entire structure were determined for our FEA maxillary central incisor model. The maximum value of the total strain energy, which is the peak restored energy, was then recorded. The damping ratio was calculated from the present FEA results for the entire tooth complex using the equation: where P is the peak strain energy (stored) and d is the corresponding value of the damping energy (dissipated).
2b) was used to calculate α and for damping ratios determined for several PDL elastic moduli reported by others. The resulting α and properties apply to the entire tooth-ligament-bone complex although each of the different material components of the finite element model should have its own damping properties. Accordingly, a simple Reuss composite model [35] was adapted to estimate the partition of overall α and between these components since they are arranged in series with respect to the loading direction. The Reuss model gives the following relationship for the effective elastic modulus of the tooth structure, E eff : where E E and V E are the elastic modulus and volume fraction of the enamel, E D and V D are the elastic modulus and volume fraction of the dentin, E PDL and V PDL are the elastic modulus and volume fraction of the PDL, and E B and V B are the elastic modulus and volume fraction of the bone. We note that the modulus of the enamel is several times greater than the other tissues in the model while it also has a relatively small volume fraction. Therefore, the first term in Eq. (4) can be ignored. In addition, the moduli for the dentin and cortical bone are roughly the same so that the second and last terms in Eq. (4) may be combined for a single elastic modulus. To a first approximation, Ashby showed that the damping ratio for soft and hard polymers is inversely related to the elastic modulus according to: where constant C has roughly the same value for a range of polymeric materials including biological materials [36]. Assuming that ligament and dentin tissue also follow this relation, Eq. (5) can be substituted into Eq. (4) to give an approximation of the effective damping ratio: where ξ PDL and ξ D are the damping ratios for the PDL and dentin, respectively. In addition, the values of ξ D is approximately equal to 10 − 3 ξ PDL based on Eq. (5) and the ratio of the elastic modulus of the PDL and that for either dentin or bone. For the incisor tooth geometry, the volume of the dentin and adjacent bone is approximately 19 to 20 times greater than that for the PDL. Substituting the corresponding volume fractions into Eq. (6) followed by substitution of 10 − 3 ξ PDL for ξ D gives Thus, the damping ratio for the PDL, ξ PDL , should be about 20 times greater than ξ eff determined for the tooth and (4) PDL together. Since α and are both directly proportional to the corresponding value of ξ, it follows that their values for the PDL should also be about 20 times greater than those determined for the entire tooth complex. Based on the ratio of PDL elastic modulus over dentin elastic modulus and Eq. (4), it follows that the α and values for the dentin and bone should be about 10 − 3 times those for the PDL and about 10 − 2 times those determined for the tooth complex. Accordingly, 20 α and 20 were used in each model for the PDL while 0.001 α and 0.001 were used for the bone and dentin as reasonable damping parameter estimates. Since these partitioned values can be considered material properties, they can also be applied in the molar models even though they were determined using damping data and tissue volume fractions for a central incisor.
The mechanical properties of the hard dental tissues (bone, dentin and enamel) were assumed to be linear-elastic and isotropic. While it is well known that the properties of these hard tissues are somewhat anisotropic, it is generally assumed unnecessary to implement anisotropic properties [4,37,38]. For the present model, this assumption appears reasonable considering the single loading direction imposed by the percussion rod. These properties as well as the density and constitutive model for each dental material are summarized in Table 1.
A large range of elastic moduli for the PDL that have been proposed previously were implemented in the present study. Specifically, we exercised models with a PDL having a typical thickness of 180 μm and each of the different reported linear elastic moduli, which range from 0.01 to 6.9 MPa, as well as Poisson's ratios of 0.45 and 0.49.

Modeling PDL Thicknesses
The percussion behaviors of the present teeth were simulated with a PDL thickness of 180 μm and one of three different PDL elastic moduli: 0.01, 0.1, and 6.9 MPa. Thicknesses in the range of 180 μm have been commonly reported in the literature for healthy teeth [2].

Results
A plot of predicted force as a function of time for the sensor region of the percussion rod is shown in Fig. 4 for the present tooth models with a PDL thickness of 180 μm and a range of different PDL linear elastic properties. Results for models in which the PDL is not present, such as in the case of ankylosis, are also shown. As expected, this figure shows that the peak force increases as elastic modulus of the PDL increases and that taking out the PDL results in the greatest percussion force. It can also be seen that the present FEA models predicted similar responses for the two values of Poisson's ratio reported in the literature, ν = 0.45 and ν = 0.49 (E PDL = 0.01 and 6.9 MPa). The models also  Fig. 5 for both an in vivo maxillary central incisor (a), and in vivo mandibular second molar (b). These data were measured following protocols described previously [26][27][28][29] and are the only experimental percussion data available for unrestored teeth that were determined to be intact based on tooth disassembly examination [26]. We note that the force amplitude of about 17 N for the molar is in reasonably good agreement with the predicted behavior for the higher values of the PDL elastic modulus (i.e. E = 0.1 to 6.9 MPa). For the incisor, a PDL elastic modulus of 0.1 MPa gives the closest approximation of maximum force for the experimental data in Fig. 5. However, it is difficult to closely correlate model results to experimental measurements since factors such as the density of the tooth-supporting bone are not easily measurable and can make a significant difference in the results. Accordingly, the data in Fig. 5 as well as any experimental results should only be used as a rough guide as to the accuracy of the present model predictions. Taking this into account, it appears that a elastic modulus of about 0.1 MPa is reasonably accurate for both the incisor and the molar.  Figure 6 contains plots of total strain energy and damping energy as a function of time for the entire tooth structure for the maxillary central incisor (a) and mandibular second molar (b) FEA models with a PDL thickness of 180 μm and elastic modulus (E = 0.1 MPa). The dashed vertical line in these plots indicates the peak strain energy for which the damping ratio is determined. Based on these results, Eq. (3) gives a damping ratio of approximately 0.1 for the maxillary central incisor model with a PDL elastic modulus of 0.1 MPa. This damping ratio falls within the range measured in vivo by Huang and coworkers [32] for a central incisor (0.09 to 0.24). The damping ratio calculated for the mandibular second molar is lower at about 0.03, which is expected due to the higher stability of the molar compared to the incisor tooth. This is evident in Fig. 6 by the substantially higher strain energy for the molar while there is lower damping energy where the strain energy is a maximum for this tooth.
The maximum principal strain in the PDL for an incisor with a PDL was surveyed for an elastic modulus of 0.1 MPa. This maximum strain was found to be less than 1.1 × 10 − 3 . As noted above, strains exceeding 0.1 are typically imposed when hyperelastic models are deemed useful [1,3,17,18,22]. Thus, the relatively low strain predicted in the present model supports the assertion that a more involved hyperelastic model is not required for the PDL under the present percussion conditions.
Force as a function of time for multiple values of the damping partition for the PDL (10, 20 and 30) are shown in Fig. 7 for both the incisor and molar models. These results show that this three-fold variation in the PDL partition factor has a relatively small effect on the maximum force predicted as well as the duration of the percussion response in time. Thus, it is reasonable to assume that the measured percussion response predicted for the rod sensor is not very sensitive to the value of this parameter and that its rough estimate of 20 should not lead to significant errors in predicting the present behavior. Figure 4 demonstrates that the incisor models with higher PDL elastic modulus values can result in force amplitudes that are nearly as high as that when this tooth is fused to the bone without a PDL. This outcome suggests that the PDL would have little effect on the maximum force with the higher elastic modulus values reported in the literature. As the elastic modulus of the PDL approaches that for the dentin and bone, the PDL becomes less distinguishable from these other tissues in terms of how it affects overall stiffness of the tooth-PDL-bone complex. Thus, it appears that a PDL elastic modulus greater than about 0.1 MPa goes somewhat contrary to the established notion that the PDL acts as a shock absorber that should significantly reduce percussive loads. Consideration of the equation of motion used in the present models to predict resultant force, F, can provide insight as to the effect quasi-static elastic modulus has on this parameter:

Discussion
where M is the mass matrix, D is the damping matrix, K is the static stiffness matrix that depends on the quasi-static value of the elastic modulus assigned to a given material, and x is the displacement. We note that the value of D depends on both M and K according to Eq. (1). Accordingly, an amplification of the effect of K on force can be substantial Fig. 7 Force as a function of time for three PDL partition factors of 10, 20 and 30 in the present FEA models of an a maxillary central incisor and b mandibular second molar under QPD conditions. The PDL elastic modulus used in these models was 0.1 MPa through the value of D when the displacement rate, ̇x , is relatively high as in the case of percussion. Therefore, even a slightly excessive elastic modulus, represented in Eq. (8) by K, can result in a force that is uncharacteristically high in the present models. It also suggests that experimental measurement of the elastic modulus of the PDL could be overestimated unless the PDL is tested at a very low displacement rate. This latter point provides some rationale for the wide range of values reported for the elastic modulus of the PDL and why our results appear more realistic when a relatively low value for the PDL elastic modulus is used. In the present work, elastic modulus is a quasi-static property that can only be measured accurately at zero or relatively low strain rates. Measurements of stiffness at higher strain rates is more reflective of a dynamic modulus which is typically greater than the corresponding quasi-static elastic modulus, as implied in Eq. (8).
The Rayleigh damping properties assigned to the PDL, dentin and bone appear to be appropriate for these tissues based on a comparison of the overall damping ratio measured for a central incisor and that generated with the present FEA models (Fig. 6). However, these damping properties are only approximations based on a simple partitioning method and a rough approximation of the relationship between elastic modulus and damping ratio. While a more comprehensive study of damping properties for these individual tissues is beyond the scope of the present work, it could lead to greater insight regarding the role of the PDL for damping impact loads on teeth in the future.
The present FE models have demonstrated reasonably good agreement with available experimental results in terms of force amplitude from a comparison of Figs. 4 and 5. Another indicator of accuracy for the present models is demonstrated by the ability of the incisor model to predict an overall damping ratio 0.1 that matches damping ratio data reported for experiments [32]. It should be noted that the time duration of percussion for the molar is in very good agreement with experimental data as demonstrated by comparing Figs. 4b with 5b. The experimental percussion duration is somewhat longer for the incisor in in Fig. 4a than that determined experimentally in Fig. 5a. However, it is possible that the incisor measured experimentally may have a thicker ligament than assumed in the present work. Thicker ligaments are often observed on anterior teeth [39]. While exercising the present models with thicker PDLs is beyond the scope of the present work, this exercise should be helpful for determining how this variable affects the percussion response in the future. Another potential source of error could be due to the damping partition approach introduced in the present work. However, Fig. 7 shows that varying the damping partition for the PDL by as much as 50% results in relatively small changes in predicted results. Thus, any errors due to the estimation of damping partitions for the present models is likely to be small.
The validity of the present models is likely limited to relatively small forces and strains (~ 0.1%) that are induced during quantitative percussion diagnostics (QPD). For higher loads and strains more representative of trauma for example, hyperelastic models may be more appropriate for predicting PDL behavior. In the case of nondestructive QPD, such stresses and strains are not feasible due to the relatively low value of the initial kinetic energy of the percussion rod [25].

Conclusion
A wide range of elastic moduli has been reported in the scientific literature for periodontal ligament (PDL) tissue. The present study examined the effect of this property in conjunction with Rayleigh damping on QPD response for a maxillary central incisor and a mandibular second molar using FEA simulations. The model predictions demonstrate that the PDL exhibits reasonable force levels under QPD conditions if the assumed quasi-static elastic modulus is relatively small (~ 0.1 MPa). However, a lower elastic modulus of 0.01 MPa resulted in unstable behavior for the central incisor. The present work also suggests that the variation in elastic modulus values reported for the PDL could be largely the result of experimental measurements that were made under different strain rates or peak strain conditions. Overall, it appears that using a quasi-static elastic modulus of approximately 0.1 MPa for the PDL in combination with Rayleigh damping leads to reasonably realistic predictions of the mechanical response of a tooth under the present percussion loading conditions. Author contributions Drs. James Earthman and Cherilyn Sheets are the inventors of quantitative percussion diagnostics. They formed Perimetrics, Inc., as a result of their research efforts with this technology and currently maintain a minority stock ownership in the company. Numerous patents have been issued and are pending. Dr. Aboozar Mapar currently works for Perimetrics, Inc. as Lead Scientist. Drs.Taheri-Nassaj and Mapar contributed essentially equally in the present work Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.