Dynamic Modelling of Tooth Deformation Using Occlusal Kinematics and Finite Element Analysis

Background Dental biomechanics based on finite element (FE) analysis is attracting enormous interest in dentistry, biology, anthropology and palaeontology. Nonetheless, several shortcomings in FE modeling exist, mainly due to unrealistic loading conditions. In this contribution we used kinematics information recorded in a virtual environment derived from occlusal contact detection between high resolution models of an upper and lower human first molar pair (M1 and M1, respectively) to run a non-linear dynamic FE crash colliding test. Methodology MicroCT image data of a modern human skull were segmented to reconstruct digital models of the antagonistic right M1 and M1 and the dental supporting structures. We used the Occlusal Fingerprint Analyser software to reconstruct the individual occlusal pathway trajectory during the power stroke of the chewing cycle, which was applied in a FE simulation to guide the M1 3D-path for the crash colliding test. Results FE analysis results showed that the stress pattern changes considerably during the power stroke, demonstrating that knowledge about chewing kinematics in conjunction with a morphologically detailed FE model is crucial for understanding tooth form and function under physiological conditions. Conclusions/Significance Results from such advanced dynamic approaches will be applicable to evaluate and avoid mechanical failure in prosthodontics/endodontic treatments, and to test material behavior for modern tooth restoration in dentistry. This approach will also allow us to improve our knowledge in chewing-related biomechanics for functional diagnosis and therapy, and it will help paleoanthropologists to illuminate dental adaptive processes and morphological modifications in human evolution.


Introduction
Dental biomechanics is a pivotal research area in dentistry to test materials for tooth restoration and to analyse the causes of dental failures [1,2]. In biology and palaeontology, dental biomechanics is the basis for understanding the relationship between tooth morphology, functional demands and constraints, and evolutionary adaptive processes in the masticatory apparatus of living and extinct animals [3][4][5][6][7]. One of the virtual and numerical tools of dental biomechanics is finite element (FE) analysis which is commonplace in dental industry [8,9] and has increasingly been used in evolutionary biological studies [7,[10][11][12][13][14].
Nonetheless, it has been recently observed that several shortcomings in FE modeling exist [15], one of which are unrealistic loading conditions. Occlusal loads are usually simplified as point loads applied parallel or oblique with respect to the long axis of the tooth [9,[16][17][18], but do not take into account the actual topography of the occlusal surface.
Benazzi and colleagues [15] pointed out that during antagonistic occlusal contact wear facets should be considered and loaded in different directions and to different extents. The authors used the Occlusal Fingerprint Analyser (OFA) software [19] to determine the individual contact areas during the power stroke, i.e. the phase of the chewing cycle when upper and lower dentitions meet and antagonistic crowns usually develop wear facets due to attritional and abrasive loss of enamel.
The OFA software is developed to virtually test, record and verify the real physiological directions of the individual power stroke movements, also reflected in the patterns and spatial positions of physical wear facets on the tooth crowns [19]. After loading surface models of antagonistic crowns aligned in maximum intercuspation into the OFA software, the user determines an assumed pathway of the power stroke with beginning and end points for each direction of power stroke movements [19]. The software then moves the models of the lower jaw along the determined pathway and as soon as tooth-to-tooth contact occurs the models are automatically guided along their relief surfaces towards the next defined path-point of the power stroke, generating a real relief guided trajectory. The accuracy of the individual kinematics can be assessed afterwards by comparing the original contact patterns (wear facet position) on the crowns with the results of virtual contact areas, detected during the calculation of the pathway trajectory. The individuals`trajectory is recorded per user-defined time step. Based on OFA results, Benazzi and colleagues [15] selected for their FE analysis the power stroke contact areas (collision surfaces) recorded in three power stroke-representative time steps, i.e., during the incursive phase I, at maximum intercuspation contact (centric occlusion), and in the excursive phase II to evaluate the stress distribution in teeth during occlusion.
Even though this approach is by far more realistic than simply reducing the occlusal force to point loads, it explains just three static moments of a sequentially changing stress scenario in a chewing cycle.
Recent contributions tried to overcome this issue by combining FE analysis and information on the human chewing process acquired either by using a jaw motion-capture system from live patients [20], or by simulating digitally the closing phase of the masticatory cycle. The latter was done by moving vertically upwards and medially the lower molar toward the upper molar to achieve maximum intercuspation [21].
Despite the obvious merits (e.g., more realistic simulations), these approaches are based on patient teeth, which impose several crucial limitations: 1) often only the external surface of the crown and not the entire tooth morphology is considered (e.g. [20]); 2) the low resolution of the underlying image data to build a volumetric mesh or coarse resolution of the volumetric mesh itself (e.g. [21]); 3) the boundary conditions, which are applied on the tooth itself and thus can lead to unrealistic local stresses (e.g. [21]). Alternatively, the complete cranium or mandible is considered [22][23][24], but again to the detriment of volumetric mesh resolution and dental structure accuracy. In general, most contributions developed so far, which combine FE analysis and dynamic loading conditions, do not employ high resolution models to accurately identify the changes in occlusal contact areas during the simulation. Yet, these are needed to assess the force distribution on the tooth surface with cascading effects on stress/strain distributions in the whole tooth and dental supporting tissues.
In this contribution we used the full individual kinematic pathway reconstructed and recorded in a virtual environment during the power stroke (i.e., the trajectory of the mandibular molar relative to its antagonist) to run a non-linear dynamic FE crash colliding test in a high resolution digital model of the upper and lower first molars and their dental supporting tissues. To our knowledge this approach provides the most realistic simulation of how chewing forces are distributed over the tooth surface, which is fundamental for an accurate evaluation of the stresses and strains occurring in teeth and dental supporting tissues under physiological loading conditions.

Micro-Computed Tomography Scanning, Segmentation and 3D Reconstruction
In 2011, the skull of a young adult female (ID = S23) selected from the archaeological skeletal sample collected by Rudolf Poech in South Africa in 1907-1909 [25], currently stored at the Department of Anthropology of the University of Vienna, was micro-CT scanned at the Vienna Micro-CT Lab, Department of Anthropology, University of Vienna, with a Viscom X8060 μCT (scan parameters: 130kV, 100 mA). Volume data were reconstructed using isometric voxels of 60μm.
The micro-CT image data was cropped mesially and distally to the socket of the right mandibular and maxillary first molar (RM 1 and RM 1 , respectively) to reduce the size of the digital models. Segmentation of the dental and bone tissues (enamel, dentine, pulp, trabecular and cortical bone) was carried out in Avizo 7 (Visualization Sciences Group Inc.). The cementum was modeled as a 0.135 mm and 0.195 mm thick layer that envelops the external root surface of the lower and upper molars, respectively (average cementum thickness from [26] (Fig 1). The periodontal ligament (PDL) was created by filling the gap between the cementum and the alveolar socket, thus varying in thickness (Fig 1). Moreover, a layer of 0.1 mm was created between enamel and dentine [27] to reflect the cushioning zone between the two tissues (enamel-dentine junction, hereafter called EDJ), following recent findings that highlight the importance of this layer to accommodate mechanical stresses [14,28,29]. The final refinement of the digital model (i.e., triangles optimization) was carried out in Geomagic Studio 2012 (Geomagic, Inc) ( Fig 1A).
In order to define an orientation system for the digital models, the best-fit plane through the cervical outline at the enamel-cementum junction (i.e., cervical plane) of the RM 1 was computed in Geomagic Studio 2012. The RM 1 was aligned with the cervical plane parallel to the xy-plane of the Cartesian coordinate system and rotated around the z-axis, so that the mesial side was parallel to the y-axis (Fig 2) [15].

Occlusal Kinematics
The dental surface models of the RM 1 /RM 1 were imported into the OFA software and the RM 1 was adjusted in maximum intercuspation position on the oriented RM 1 . Following general  knowledge about functional anatomy and biomechanics of mastication [30], an assumed pathway for the power stroke was set for the lower molar coming from distobuccally upward and slightly anteriorly until maximum intercuspation, and then downward towards slightly mesiolingually. Thanks to collision detection, deflection and break free algorithms implemented in the OFA software, one model is allowed to move against the other, ultimately providing a kinematic analysis of the surface contacts between RM 1 and RM 1 during the power stroke [15]. In detail, after having set a starting point and an endpoint for the movements, as soon as the lower and upper teeth collide, the travelling model (RM 1 ) is deflected and guided by the antagonistic crown relief producing its individual power stroke trajectory. Contact areas are marked in different colors and the occlusal pathway trajectory is recorded, informing on the occlusal movements during the incursive (phase I) and excursive (phase II) of the power stroke ( Fig 3A  and 3B; S1 Video). The calculated collisions were compared to the wear facets (macrowear pattern) present on the occlusal surface of the real tooth in order to verify the accuracy of the generated trajectory. The starting and endpoints of the trajectory calculation were repeated three times until the resulting collisions matched with the functional macrowear pattern on the original specimen. The verified individuals`pathway trajectory was exported as a list of Cartesian coordinates (x,y,z) and used in the FE software to guide the RM 1 in the crash colliding test.
The OFA software allows for recording the occlusal pathway trajectory (or power stroke, i.e. phase I and phase II), but the duration of the chewing cycle (i.e., closing phase + power stroke + opening phase) is automatically, yet unphysiologically determined by the software [19]. Therefore, we used information gathered from the literature to calibrate the duration of each of the phases of the chewing cycle. Specifically, modern human jaw movements recorded with a magnet-sensing jaw-tracking system suggest that the power stroke may last about 289 ms (measured from [30]; his Fig 1-15). Since in the OFA simulation the closing and opening phases were computed to be 0.5479 and 0.3493 times the power stroke, respectively (Fig 4), we estimated a total movement time of 548.311ms.

Volumetric Mesh Generation and FE Analysis
The surface models were imported into HyperWorks Software (Altair Engineering, Inc.), where volumetric meshes were created (Fig 1B). Two simulations with 4-noded 2,846,131 tetrahedral elements (hereafter called "coarse" mesh) and 8,832,739 tetrahedral elements (hereafter called "fine" mesh), respectively, were carried out to verify the convergence of the results. For the hard tissues (enamel, dentine, EDJ, cementum, bone), tetrahedral elements with 6 rotational and translational degrees of freedom (DOF) at each node were used [32], hence improving accuracy of the simulation. For the soft tissues (pulp, PDL) each node had 3 DOFs and tetrahedral elements with constant pressure formulation for incompressible material were used. The number of elements arranged across the 0.15-0.85 mm PDL thickness varies between 2 and 8, accordingly. Only a single layer of elements was arranged across the 0.1 mm thick EDJ. Adaptive meshes were not required to capture the thin geometric features of PDL and EDJ, because in both cases the mesh size satisfied the mesh quality requirements (e.g., warpage, aspect ratio).
Material property values (elastic modulus E, Poisson's ratio and density ν) were collected from the literature (Tables 1 and 2). The same mechanical properties and density were used for cortical and trabecular bone (hereafter called "bone"), following the suggestion by Hall [43].
For the EDJ, E and ν were computed as the average values between enamel and dentine [14]. All materials were considered homogeneous, linearly elastic and isotropic, which is an evident simplification but it is comparable with current studies in dental biomechanics [44].
Regarding the boundary constraints, the superior surface of the maxilla was totally fixed in X, Y and Z directions, as during mastication the maxillary bone does not move. In contrast, since the mandible moves during mastication and is subjected to torsional loading and sagittal bending [45,46], for all the surface nodes of the mesial and distal surfaces of the mandible section the velocities calculated in X, Y and Z directions were imposed simultaneously (see below for details about the velocity).
The explicit scheme was used to solve the governing equations [47]. This method is efficient for highly non-linear problems, especially when dealing with complex contacts and large deformations. Moreover, it is particularly reasonable for FE models with several million elements, for which the implicit scheme has problems with computational cost, convergence or contact failure [48].
In FE non-linear simulations, the movement time (see above) is subdivided in multiple time steps (Δt), which should be small to guarantee 1) a linearly approximated behavior of the structure and 2) a stability of algorithm. In general, Δt is automatically calculated by the software. The initial default value is Δt = 2.1E-09s, but to save computational time Δt was scaled to 5.0E-08s for the open and close phases (mass scaling, i.e. the software increases the density of the components of the model; see [49]). The default Δt assigned by HyperWorks was used for the power stroke.
While the time scaling (i.e., increasing velocity) does not affect the linear elastic material used for the FE models (i.e., there is no strain rate dependency), high scaled velocity may introduce nonrealistic dynamic effects that can alter the accuracy of the solution. However, limiting scaled velocity to less than 1% of the wave speed in the simulated materials does not affect significantly the accuracy [49]. After several trials with different scaling factors (145, 100, 80, 60 and 50), the time scaling factor was set to Fsc = 50, because it assured the smallest acceptable kinetic energy in the closing phase. The maximum velocity after scaling was still smaller than 1% speed of sound (Vsound) of pulp tissue (Table 2), the latter was computed following [42]: where E = Young's modulus (MPa), v = Poisson's ratio and ρ = density (kg/mm 3 ). Moreover, we introduced a diagonal damping matrix C, proportional to mass matrix M, in the dynamic equation to damp spurious oscillations: where: K is the stiffness matrix; the time-dependent vectors u(t), _ u(t) and ü(t) are nodal displacements, velocities and accelerations, respectively; F(t) is time-dependent external load vector; B is the relaxation value (default value 1); T is the period to be damped (imposing T = 0.0035s, close to the highest frequency of the model from PDL).
The direct integration method used in the HyperWorks solver package is derived from Newmark time integration scheme. Velocity and displacement at time t n+1 = t n + Δt are calculated: The integration scheme used by HyperWorks' solver RADIOSS is based on the central difference integration algorithm with coefficients β = 0 and γ = 0.5. It is conditionally stable when the time step Δt is small enough to integrate accurately the response in the highest frequency component.
A penalty type formulation [32], which is accounted for by the type 7 interface of Hyper-Works software, allows sliding between a master surface (occlusal surface of the M 1 enamel) and a set of slave nodes (occlusal surface of the M 1 enamel). This interface has spring stiffness when a slave node penetrates the interface gap. At a contacted node the elastic contact force has two components: normal and tangent. The normal component of the contact force is a function of contact-related parameters as: with -g is the gap value (i.e., to determine when contact between slave nodes and master surface occurs), which by default was defined by the software HyperWorks as one tenth of the smallest side of all master element segments (g = 6.636E-03 mm); -p is the penetration of a slave node into a master segment (g ! p ! 0). Penetration is checked in each time step of the power stroke. When a penetration is registered, the force proportional to the penetration depth (6) is applied to resist the penetration; -K is the interface spring stiffness: K = min(K m , K S , K min ). The master stiffness is computed as The slave stiffness is an equivalent nodal stiffness computed as ; Minimum stiffness is K min = 84100, equal to Young's modulus of enamel. In the above formulae St FAC = 1 (stiffness factor); Bu is the bulk modulus of enamel material; S is the master element area; V is the master element volume. Overall, the interface spring stiffness is related to material parameters, ultimately favoring the efficacy of the penalty method, because the stiffness of contacting enamels has the same order of magnitude.
We introduced a Coulomb friction law, using a friction coefficient μ = 0.2, a value found for wet conditions [50]. The tangent component of the contact force is then calculated as: A nodal contact force is the resultant force of a normal component (6) and a tangent component (7). Summation of contact forces at each node produced by the enamel contact interface provides the contact force between the lower and upper enamel caps, which is equivalent to the chewing force. During the power stroke values of the chewing force are changing according to the changes of contact area size and spatial position between enamel caps. Therefore, the power stroke chewing force is time dependent, whereas in both the closing and opening phase the chewing force is equal to zero.
Gravity force was constantly applied on all parts of the model, even though we observed that gravity effect is not important for this simulation.
Results for were qualitatively and quantitatively evaluated according to the first maximum principal stresses criterion for brittle materials (as stress wave propagation is fundamental to investigate dynamic effects of chewing load on materials) [51,52], but for the PDL (i.e., soft tissue) the first principal strain was evaluated.

Sensitivity of FE Methodology
When experimental data are not available, fundamental checks should be undertaken to reduce the error of the explicit simulation (S1- S5 Figs). First, the mesh size of the FE model must be fine (namely small) to prevent unphysical stress even when the velocity is high. To address this issue, the average mesh size value is h (height) = 0.13mm, with a total number of 8,832,739 tetrahedral elements. Second, the quality of the simulation's results was estimated through the evaluation of the 1) energy balance (S1 Fig); 2) ratio of kinetic energy to internal energy (S2 Fig); and 3) time history of kinetic energy, internal energy and contact energy [48] (S1 Text). Moreover, differences between the two simulations (using the "coarse" and "fine" meshes, respectively) were <5% for displacement (S3 Fig) and internal energy (S4 Fig), confirming the convergence of the results. Only results obtained using the "fine" mesh are reported below.

Results
The maximum contact force achieved by solving the system is 923N (Fig 5). The force increased during phase I, reaching its greatest magnitude during maximum intercuspation contact. Even though the force reduced during phase II, values remained relatively high (Fig 5).
The displacements of the RM 1 and RM 1 along x-, y-and z-axes are summarized in Figs 6 and 7, respectively. For the RM 1 , the larger contact areas on the lingual cusps displaced the tooth mainly in the negative direction of the x-and y-axes, i.e. mesiolingually (Fig 6; S2 Video).
Due to dynamic effects the displacements of the RM 1 during the closing and opening phases are different from zero. However, during the power stroke the dynamic effects were constrained by contacts between antagonists, and the RM 1 displaced in the positive direction of the x-and mainly y-axes, i.e. distobuccally (Fig 7).
The distribution of maximum principal stresses on the enamel followed the same pattern for both RM 1 and RM 1 . Contact pressure causes tensile stresses concentrated in the grooves of the occlusal surface, and their magnitude and area of interest reached its peak during maximum intercuspation (Fig 8; S3 Video). At mid-time of the power stroke (t % 0.375s) the RM 1 is mostly moving in the negative direction of the y-axis, causing high contact pressure and a sudden increase of tensile stresses on the distal cusps (mainly the hypoconulid), affecting also the grooves that surround these cusps (Fig 8B). As a consequence, tensile stresses increased also in the mesiolingual wall of the root (Fig 9) and, to a lower extent, in the superior alveolar margin of the mandible (S6 Fig; S3 Video). In the PDL tensile strains were highest in the distobuccal region (Fig 9).
Towards the end of the closing cycle, the RM 1 continues moving in the negative direction of the y-axis, leading to maximum contact pressure on the hypoconid (at t % 0.425s). Consequently, tensile stresses increased in the occlusal basin of the crown (Fig 8), in the lingual aspect of the root (Fig 9) and, more generally, in the upper half of the cortical bone of the mandible  As far as the PDL is concerned, tensile strains remained higher in the buccal side than in the lingual one (Fig 9).

Discussion and Conclusions
Currently, most of the FE simulations of teeth are limited to static loading conditions [7,12,13,53], which reflect only a partial view of the dental hard and soft tissues' functional response. Few studies combine FE methods and dynamic loading conditions [19,20], but the analyses were restricted to coarse volumetric meshes that simplify the complex dental architecture and, importantly, did not envision the physiological sequential changes of occlusal contacts occurring during the power stroke. Yet, this information is pivotal to simulate realistic loading conditions.
The dynamic approach presented here is the logical advance following recent contributions that combined FE methods and occlusal macrowear analysis using high resolution digital models of teeth and dental supporting tissues [14,15,53,54]. Using a non-linear dynamic FE crash colliding test based on the trajectory computed in the OFA software, we were able to compute tensile stresses and displacements of the FE models during the complete power stroke. The maximum contact force registered solving the system (923N) is in agreement with maximum values obtained during in vivo bite force measurements [55][56][57][58]. Moreover, our findings support previous FE studies that suggest both the importance of the PDL for load transfer to the alveolar wall [59] and the increase of tensile stresses in the upper half of the mandible, mainly in the lingual side [53].
This approach will be very useful in research pertaining to clinical dentistry, because the individual direction and position of occlusal forces are important factors for, e.g., estimating the long-term success of dental and endodontic treatments, dental implant designs, dental prostheses. A precise and individual occlusal surface reconstruction is crucial for a functional and balanced stress distribution in the dental arches with an optimum chewing efficiency, during the maintenance phase for patients with periodontitis, as well as for an understanding of the etiology of some very common dental pathologies, such as non-carious cervical lesions and abnormal tooth wear [60][61][62]. Moreover, our dynamic approach will enhance interpretations of investigative results related to evolutionary biology of teeth, e.g. to unravel the adaptive processes and significance of dental and resulting masticatory variations in respect to the development of various dental features, e.g. such as protostylid and mid-trigonid crests as well as enamel thickness and root shape.
It is important to note that the large size of the FE models (which is pivotal to maintain a detailed representation of the various tissues) and the complex computation to solve the nonlinear equations require powerful computer facilities and is very time consuming. We are confident, however, that future technological advancements will address this practical issue, ultimately making the entire process less demanding and allowing us to implement food items between the occlusal surfaces in the simulation as well. Ultimately, experimentally derived in vitro or in vivo data on tooth function are required to substantiate our findings.