https://orcid.org/0000-0001-8060-2471
Figures
Abstract
Studying different types of tooth movements can help us to better understand the force systems used for tooth position correction in orthodontic treatments. This study considers a more realistic force system in tooth movement modeling across different patients and investigates the effect of the couple force direction on the position of the center of rotation (CRot). The finite-element (FE) models of human mandibles from three patients are used to investigate the position of the CRots for different patients’ teeth in 3D space. The CRot is considered a single point in a 3D coordinate system and is obtained by choosing the closest point on the axis of rotation to the center of resistance (CRes). A force system, consisting of a constant load and a couple (pair of forces), is applied to each tooth, and the corresponding CRot trajectories are examined across different patients. To perform a consistent inter-patient analysis, different patients’ teeth are registered to the corresponding reference teeth using an affine transformation. The selected directions and applied points of force on the reference teeth are then transformed into the registered teeth domains. The effect of the direction of the couple on the location of the CRot is also studied by rotating the couples about the three principal axes of a patient’s premolar. Our results indicate that similar patterns can be obtained for the CRot positions of different patients and teeth if the same load conditions are used. Moreover, equally rotating the direction of the couple about the three principal axes results in different patterns for the CRot positions, especially in labiolingual direction. The CRot trajectories follow similar patterns in the corresponding teeth, but any changes in the direction of the force and couple cause misalignment of the CRot trajectories, seen as rotations about the long axis of the tooth.
Citation: Gholamalizadeh T, Darkner S, Søndergaard PL, Erleben K (2021) A multi-patient analysis of the center of rotation trajectories using finite element models of the human mandible. PLoS ONE 16(11): e0259794. https://doi.org/10.1371/journal.pone.0259794
Editor: Antonio Riveiro Rodríguez, University of Vigo, SPAIN
Received: April 26, 2021; Accepted: October 26, 2021; Published: November 15, 2021
Copyright: © 2021 Gholamalizadeh et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The dataset used in this study is made publicly available and can be obtained from Electronic Research Data Archive at the University of Copenhagen (https://doi.org/10.17894/ucph.04e91c97-5c5d-45d7-afd3-c5e0b5953f58). All other relevant information is provided in the manuscript.
Funding: This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 764644 awarded to TG (https://cordis.europa.eu/project/id/764644). This paper only contains the author’s views and the Research Executive Agency and the Commission are not responsible for any use that may be made of the information it contains. 3Shape A/S also provided financial support in the form of salaries for authors TG and PS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: Authors TG and PS are affiliated with 3Shape A/S, which provided support in the form of salaries and data for this study. The commercial affiliation of the aforementioned authors does not alter our adherence to PLOS ONE policies on sharing data and materials. Additionally, there are no patents, products in development or marketed products associated with this research to declare.
Introduction
Finite element (FE) modeling is a widely used computational method for the analysis of the reactions to real-world physical effects such as forces and biomechanical responses occurring during treatments in medicine and dentistry [1] that attempts to solve partial differential equations numerically, based on reconstructing the desired geometry and discretizing the domain into a finite mesh, with minimum need for clinical trials of patients [2]. One of the main goals in FE-based modeling of the tooth and its supporting complex, i.e., periodontal ligament and adjacent alveolar bone, is to improve tooth movement prediction performance in orthodontic treatments. Compared to clinical studies that aim, among other things, to assess the accuracy of digital planning in dentistry [3] or investigate anatomical characteristics for optimal occlusions [4] using different patient scans, computational modeling allows for more efficient plannings in orthodontic treatment for correcting dentofacial anomalies using more generic solutions.
Tooth movement modeling in an uncontrolled tipping scenario using a perpendicular loading system has been widely used for initial tooth movement simulations due to its simplicity [5]. However, in addition to the uncontrolled tipping movements, different tooth movements such as the pure translation, crown/root tipping, intrusion/extrusion, or a combination of them are typically required for tooth position correction in orthodontic treatments. Each of these movements can be described based on the position of the center of resistance (CRes) and center of rotation (CRot) with respect to the tooth geometry [6, 7]. For example, a pure translation, root tipping, and crown tipping can result in CRots located at infinity, crown, and root apex, respectively. The CRes of a tooth within its supporting complex can be seen as the center of the mass of a free rigid body [8]. To be more specific, the tooth CRes is a point at which applying any forces would always result in a pure translation of the tooth [6]. Although the location of the CRes has been investigated in the literature either experimentally in vivo [9–11], analytically [12, 13], or computationally using finite-element (FE) models [7, 8, 14], estimating the exact position of the CRes is challenging as it depends on the root shape, the anatomy and mechanical properties of the supporting complex [7, 15], the direction of the tooth movement and force [6, 7, 16–18]. More specifically, since the CRes may not exist as a single point in 3D space [8] due to the asymmetric geometries, it is proposed to use a volume of CRes instead [14].
The CRot is another important concept in tooth movement analysis and is a point at which the movement of the tooth can be defined as pure rotation. Different approaches can be used for determining the position of the CRes and CRot. One simple approach for determining the tooth CRes is exerting a couple in different directions. A couple consists of a pair of forces with the same magnitude and parallel line of action, but in opposite directions where they are not collinear [19]. Therefore, utilizing only a couple results in a pure rotation in the tooth where the position of the CRot coincides with the position of the CRes.
The CRot position depends on the load system applied to the tooth. The loading system can result in some moments and forces in different directions, and due to its complexity in a 3D coordinate system, it can be decomposed into three planes of the tooth. Hence, tooth movement can be analyzed by decomposing each moment-to-force ratio (M: F) on each principal plane of a tooth, where each M: F on a plane can be defined based on the force vectors located on the plane and moments perpendicular to the studied plane [20].
A classical theory on the relation between the applied load and the type of the tooth movement is Burstone’s formula [21]. In this theory, the effect of the various moment-to-force ratios (M: F) is studied on the position of the CRot for a canine with a parabolic root, while the force is applied perpendicular to the long axis of the tooth at the bracket level. The study introduces the following formula
(1)
where h is the distance from the alveolar crest to apex, and D is the perpendicular distance between the CRes and CRot. Note that M: F has a unit of millimeters where F corresponds to the magnitude of the applied force and M denotes the moment of the couple mc exerted to the tooth in order to counteract the moment caused by the applied force mf [19]. According to this theory, specific values of M: F ratios always correspond to specific types of tooth movement. That is, the location of the CRot for a specific tooth is only dependent on the moment-to-force ratio of that tooth.
The classical theory has been widely studied by using FE models [17, 20, 22]. Cattaneo et al. [22] study the influence of M: F on the same force magnitude and analyze the effect of applying various load magnitudes on a mandibular premolar and canine teeth with a constant M: F. The study reports that a constant M: F with different force magnitudes can result in different types of tooth movements. However, this finding does not follow Burstone’s formula, and this is due to considering a nonlinear material model for the PDL layer where the location of the CRot is dependent on both the M: F and the magnitude of the applied load to the teeth.
Different studies have utilized different M: F ratios for the pure translation of various teeth. For example, unlike the generally accepted M: F of 10 and 12 for pure translation of the premolar and canine, Cattaneo et al. [22] suggest M: F of 9 and 11, respectively, for a mandibular premolar and a canine. Also, it is recommended [18] to avoid using a universal M: F and CRes due to the patient-specific geometries, reporting M: F of approximately 8.8, 9.7, and 10.2, to simulate the pure translation of maxillary first premolar, lateral incisor, and canine.
Unlike the abovementioned studies that are limited to single-plane analysis, Savignano et al. [20] investigate the effect of the force system directions on the tooth movement of a maxillary first premolar by performing analysis on principal planes of the tooth. More specifically, the mesiodistal, buccolingual, and occlusal planes are used to study the location of the CRot using the projected axis of rotation on these planes for various M: F and directions. This study suggests a nonlinear relation between the direction of the force system and the projected CRots. However, the mentioned study is limited to a single tooth analysis with a limited load magnitude range applied to the CRes due to using a linear elastic material model for the PDL layer. Although two patients are considered in a recent work [23], an inter-patient analysis is missing in general.
The existing studies do not represent the effect of geometrical variations on the M: F ratios and the position of the CRot under the same analysis setup, boundary conditions, and loading system for different patients. To this end, in this work, we consider patient-specific full dentition computational models of three human mandibles to investigate the position of the CRots in different patients’ teeth. For a comprehensive and consistent inter-patient analysis, we register different patients’ teeth to a corresponding reference tooth of the first patient using an affine transformation and use this transformation to map the selected force directions and force application points on the reference teeth into the registered teeth domain. Moreover, we model the clinical forces with no specific limitation in the force magnitudes by using a hyperelastic material model for PDL tissue and investigate the position of the CRot in a 3D space to better represent the resulting rotation. Finally, we assess the influence of the couple directions on the positions of the CRots of a specific tooth and show that equally rotating the direction of the couple about the three principal axes results in different patterns for the CRot positions, especially in labiolingual direction. Modeling different types of tooth movements using a more realistic force system, i.e., tipping and couple forces applied to the tooth crown, and investigating the effect of these couple forces on the CRots can help us to better estimate the force systems needed to achieve the desired tooth movement.
Materials and methods
This section first reviews the details of the FE models used in this study including the geometry reconstruction, mesh generation, boundary conditions, contact definitions, and material models. Next, it focuses on describing an approach for setting up a consistent loading condition for different patients’ teeth, both for intra- and inter-patient analysis. Finally, it briefly specifies how the CRot is computed and illustrated in 3D.
Geometry reconstruction
Three patient-specific and anatomically accurate FE models of the human mandible composed of mandibular teeth, corresponding PDLs, and bone are considered in this work. To have enough geometrical variations in the dataset, the three patients’ scans are chosen of various crown heights, root lengths, and teeth sizes. The scans are used from the 3Shape A/S private dataset collected from different clinics by orthodontists as a part of treatment. More specifically, the anonymized Cone-Beam Computed Tomography (CBCT) scans of patients stored in Digital Imaging and Communications in Medicine (DICOM) format used in this study contain no sensitive personal information including gender and age or details of the scanner device utilized for image acquisition.
The utilized patient scans labeled as Patient 1, Patient 2, and Patient 3 have an isotropic voxel size of 0.3mm, 0.3mm, and 0.15mm, respectively. For a consistent and detailed geometry reconstruction, a cropped region of interest (ROI) of every patient scan including full dentition and mandible is upsampled to an isotropic voxel size of 0.15mm. Further details of the utilized scans can be found in Table 1. Note that the patient-specific geometries are reconstructed by importing DICOM files to 3DSlicer [24], resampling the ROIs, and segmenting the teeth and bone geometries in CBCT scans using a semi-automatic watershed algorithm.
The accuracy of the annotated scans is verified by a clinical expert and any required modifications are applied to the segmented regions accordingly. A general criterion for the verification is precise segmentation of roots, crowns, and bone cervical regions including the teeth sockets, and the miss-annotated regions indicated by the expert are revised until reaching the criterion. The segmented geometries are then exported as surface meshes in STL file format. Since the resolutions of the original scans are not high enough to properly represent the thin layer of the PDL, annotating the PDL layer using the available CBCT scans is not feasible. The width of the PDL layer is a tooth root dependent factor which can vary between 0.15mm and 0.38mm with an average suggested as 0.2mm [25]. This average thickness is used in the literature [25–27] to generate the PDL geometries by uniformly extruding the teeth roots in Meshmixer [28]. Therefore, we first re-mesh the teeth and bone geometries using a uniform and adaptive mesh, respectively. Next, the obtained bone surface is offset by 0.2mm in the reverse direction of the surface normals to create the required space for the PDL layers in between the teeth roots and the regenerated bone. Finally, the PDL geometries are produced by extruding the teeth roots, filling the generated gap in between the teeth and bone. It should be noted that the generated surface mesh for each PDL geometry includes a uniform mesh with two elements in the PDL thickness. The reconstructed geometries have been made publicly available and can be obtained from Electronic Research Data Archive at the University of Copenhagen under the OpenJaw Dataset (https://doi.org/10.17894/ucph.04e91c97-5c5d-45d7-afd3-c5e0b5953f58).
After preparing the surface meshes, high-quality tetrahedral meshes are generated by preserving the surface meshes using TetGen [29]. The mesh quality required for an FE analysis can vary depending on the application and utilized numerical methods [30]. In general, a regular tetrahedron has the highest mesh quality, and the main guideline is to avoid using low-quality tetrahedra with small or large dihedral angles, as they can affect the accuracy of the numerical methods [30]. Therefore, we use an upper limit constraint of 1.2 for the radius-edge ratio, as the ratio of the circumscribed sphere’s radius and the shortest edge of the linear tetrahedral meshes. Additionally, the quality of the generated volumetric meshes are assessed using four different quality measurements, i.e., the volume-edge ratio [31, 32], radius-edge ratio [33], radius ratio [34, 35], and volume-area ratio [34]. The volumetric meshes are then imported into the FEBio software package [36], which is open-source software for nonlinear FE analysis in biomechanics, to set up a reproducible FE model for each patient. A summary of the utilized mesh properties is provided in Table 2.
A mesh convergence study is performed to obtain an optimal mesh size for the FE models [5] by iteratively increasing the number of elements per step with a factor of 2. We continue the process until the relative stress error does not exceed 4% of maximum von Mises stress. This study is concerned with smooth bone meshes where there are no sharp elements on the meshes. Therefore, the model does not experience extreme local stresses in a single node or element. However, in nonsmooth geometries with sharp edges, more robust methods can be applied to avoid any extreme local stresses as outliers [37, 38]. Also, note that using more elements in the PDL thickness would exponentially increase the total mesh size of the full jaw models and the required computational time. Besides, based on our experience, using smaller elements would result in extreme element distortion, causing problems in FE model convergence especially for higher load values. Similar behavior is also seen in [39] for hexahedral elements.
Material properties
We assume teeth as rigid bodies with six degrees of freedom to simplify the proposed model. Also, we assume the center of mass as the CRes for each tooth which can automatically be calculated using FEBio user-defined options [40]. Since the deformation of the bone tissue is negligible under the orthodontic forces used in this study, no distinction is made for the cortical and trabecular bone [41, 42], and isotropic elastic material model is used for the homogeneous bone geometry. Due to the importance of the PDL tissue in transferring loads from the tooth to the alveolar bone [39, 43], a Mooney-Rivlin Hyperelastic (MRH) material model is used to simulate its nonlinear behavior [42, 44]. Furthermore, for the simplicity of the model and following the literature [14, 20, 23, 39], the gingival tissue is discarded from our computational model due to its extremely low elasticity modulus compared to that of the other tissues, i.e., PDL, bone, and tooth.
Boundary condition and contact definition
A rigid contact is defined between rigid teeth and PDLs, and a tied facet-on-facet contact with an augmented Lagrangian method is used in PDL-bone interfaces to model adhesion in these two interfaces, as their corresponding surfaces do not have any sliding or separation. Other parameters, such as boundary conditions, edge length of the elements, and material properties are chosen as mentioned in Table 2. Also, a Dirichlet boundary condition is applied to all nodes at the bottom surface of the bone in all directions to fix the jaw model.
Loading conditions for different scenarios
Three different scenarios with varying loading conditions are studied in this work, as shown in Fig 1. In the first scenario, various M: F values are examined for different patients’ teeth using a consistent loading system. In the second scenario, the influence of the direction of the couple is studied on a specific tooth. As our third scenario, the loading conditions are applied to the CRes, which makes us able to compare our simulation results with those obtained by Savignano et al. [20, 23].
A: An illustration of the principal axes and planes of premolar. B to D: Loading conditions used in three different scenarios.
Scenario 1: Consistent loading conditions for intra- and inter-patient analysis.
In the first scenario, to perform a systematic intra- and inter-patient analysis and to avoid potential errors caused by the selection of the load points and force directions, different teeth of each patient are registered to a reference tooth of the same type. To do so, the left teeth of one of the patients (here, Patient 1), with the Universal Numbering (UNN) ID of 18 through 24, are selected as the reference teeth. Next, the mandibular teeth of all other patients on both sides, as well as the right mandibular teeth of Patient 1, are all registered to the corresponding reference teeth. This allows us to analyze more tooth models of the same type and study the position of the CRots for different teeth of different patients.
Before any tooth registration, the jaw models of different patients are aligned based on rigid body transformation using the Iterative Closest Point (ICP) algorithm [45, 46]. Moreover, the Coherent Point Drift (CPD) algorithm [47] is used for teeth alignment based on an affine transformation of the corresponding left/right teeth of different patients. The registration processes are all performed in MATLAB by considering vertices of each surface mesh as a set of a point cloud. Compared to the ICP algorithm, CPD is computationally intensive, yet it is more robust to noise and outliers. Therefore, we prefer using CPD for teeth registration which includes fewer data points than jaw models. It should also be noted that all teeth on the right side of the mandibles, with the UNN of 25 through 31, are horizontally flipped with respect to the sagittal plane of each patient’s jaw before applying the affine transformation using the CPD method. This is to ensure that the labial surfaces of the registered teeth are properly aligned to the labial surfaces of the reference teeth.
After registering all patients’ teeth to the reference teeth, three points, at which the force and couple are applied, are selected on each reference tooth. One point is set for tipping force, and two points are used for the couple. The center of the middle third of the crown is set for the tipping load, and two points are selected with a 0.5mm distance below and top of the tipping load point for the labiolingual and linguolabial couple, respectively. The normal vector of the reference mesh at the tipping point is used for both directions of the force and the couple. More specifically, the direction of the computed normal vector is used for the linguolabial couple force on top, and the opposite direction is used for both the tipping force and the labiolingual couple. Nodal loads are utilized to simulate the mentioned forces in the FEBio framework.
The corresponding three load points are then obtained based on finding the closest mesh points in each registered tooth mesh. These points, as well as the direction of the forces, are finally transformed back to the original coordinates of each patient’s tooth. This would ensure unbiased loading conditions for the intra- and inter-patient analysis. Fig 2 displays how registration is employed for defining consistent loading conditions and calculating the CRots using FE models of mandibles from three patients.
An illustration of the utilized method for calculating the CRots using finite-element models of mandibles from three patients under consistent loading conditions.
Multi-patient analyses are conducted on different teeth to evaluate the effect of the couple’s magnitude on counteracting the moment produced by the tipping load. In this scenario, the simulations are performed with a constant tipping force of 1N and varying couple magnitudes of 1N to 13N with 1N intervals. Note that the applied force points are fixed for both the couple and the tipping load. Hence, increasing the load magnitude of the couple will increase the resulting moment on the tooth according to the force-moment relation M = F × d, where d is the distance between the applied force point and the CRes of the tooth. Finally, the corresponding CRots are computed for each patient’s teeth and the same transformations are used to align different teeth of the same type and their CRots.
Scenario 2: The effect of the couple’s direction on the position of the center of rotation.
In the second scenario, as an example, the left premolar of Patient 1 is selected to study the influence of the couple’s direction on the position of the CRot. The three load points, the initial load directions, and the direction and magnitude of the tipping load remain the same as in the first scenario. The couple’s directions are rotated about the mesiodistal, the labiolingual, and the long axis of the tooth with different degrees varying from 0 to 90 degrees with 10-degree increments. To investigate the trajectories of the CRots in 3D space, for each direction, we conduct experiments with incremental couple magnitudes of 1N to 13N with intervals of 0.0625N. This generates a dense representation of the CRot trajectories in 3D space.
Scenario 3: The effect of the force’s direction on tooth movement.
As our third scenario, we investigate the influence of the force’s direction on tooth movement and the position of the CRot, and compare our results with those obtained by Savignano et al. [20, 23]. More specifically, the effect of the force direction is investigated at the CRes of the left premolar of Patient 1, while the load direction changes in one of the three principal planes of the tooth, and a perpendicular moment is applied to the studied plane. Additionally, instead of using the couple for generating the moments on the tooth crown, we follow the same scenarios of the aforementioned study and use a load and a moment applied to the CRes, using prescribed rigid force constraint and prescribed rigid torque constraint, respectively, in FEBio framework. The direction of the load changes with 10-degree increments, and the M: F changes from -12 to 12 with 2mm intervals.
Computing the center of rotation
The axis of rotation is computed for each tooth by considering displacement vectors of two nodes arbitrarily selected on the tooth crown. To this end, the intersection line of the perpendicular bisector planes of the displacement vectors is obtained. Any point on the intersection axis can be assumed as the CRot in 3D space, and to find a unique CRot in 3D, the point with the closest distance to the CRes of the tooth is selected on the rotation axis. This allows for representing the axis of the rotation in 3D as a single point and for better analyzing the influence of the couple directions on the location of the CRot with respect to the CRes.
Results
We conduct our experiments on a 3.4 GHz processor with 64 GB of RAM, which takes about an hour to solve the FE model. First, we focus on the analysis of the position of the CRots for various couple magnitudes of all patients’ teeth. Second, the influence of the couple directions is studied on the position of the computed CRot of the left premolar of Patient 1. Third, we assess the influence of the force direction with fixed moment direction on the position of the CRots. In all the experiments, the position of the closest point to the CRes on the rotation axis is considered as a single point CRot in 3D coordinate system for each tooth. Fig 3 illustrates the position of the CRots of each tooth for varying magnitudes of the couple from 1N to 13N with 1N increments. As can be seen, the CRots follow parabolic shapes with varying slopes from one tooth type to another.
Calculated CRots for different couple magnitudes changing from 1N to 13N with 1N intervals shown in 3D space for all teeth with transparent teeth geometries, Left: Results of the three different patients. Right: A closeup view of the results for Patient 1.
To better represent the inter-patient results, the utilized patients’ teeth and their CRots were registered to the corresponding reference teeth. Fig 4 shows the registered teeth geometries, with transparent shapes located on either left or right side of the mandible, and their CRots for the utilized patients. Note that the third molars were missing in all cases. Hence, the different teeth types are represented from the second molars to the central incisors. As can be noted, the CRot trajectories experience similar patterns in each tooth type.
Trajectories of the obtained CRots for the different patients’ teeth registered to the reference teeth of the same types and shown as transparent teeth. A: The second molars. B: The first molars. C: The second premolars. D: The first premolars. E: The canines. F: The lateral incisors. G: The central incisors. A consistent loading condition is applied to all teeth of the same types, and the magnitude of the couple changes from 1N to 13N. Note that a few extreme outliers are discarded for a better illustration purpose.
In practice, more CRot data points are required for a better illustration of the trajectory curves in 3D space. This in turn requires sampling more data points by using smaller intervals of the couple magnitudes. As mentioned before, we conduct this experiment using the left premolar of Patient 1 and investigate the CRot trajectories influenced by the couple’s directions based on smaller increments in the magnitude of the couple. More specifically, the directions of the couple are rotated about the mesiodistal, labiolingual, and long axes of the tooth from 0 to 90 degrees with 10-degree increments. Fig 5 illustrates the CRot trajectories for different directions of the couple rotated about the principal axes of the tooth.
Directions of the couple rotated about the tooth principal axes with various angles and the corresponding CRots for different couple magnitudes. Left: An illustration of how the couple forces are rotated about mesiodistal, long axis, and labiolingual of the tooth. Right(A to C): The corresponding CRots for the rotated couple about mesiodistal, long axis, and labiolingual of the left premolar of Patient 1, respectively. Note that the initial direction of the couple and tipping load are set based on the surface normal direction.
As our third analysis, we examine the effect of the force direction on the position of the CRot by considering three principal planes of the tooth. We use the same scenarios followed by Savignano et al. [20], where the force and moment are applied to the CRes with a force of 1N magnitude. Fig 6, shows the effect of the load’s direction on the position of the CRot, where the rotation range of the direction of the force system can vary from one axis of the tooth to another with 10-degree increments. As it can be seen in both subfigures, the distance between the CRots of the same direction changes nonlinearly with linear increments of the M: F. The nonlinear variation increases when the force direction is rotated from the mesiodistal or labiolingual axis direction towards the long axis direction. Additionally, the distance of the CRots associated with the same M: F ratios changes nonlinearly for linear increments of the rotation degree in 3D space. This pattern is also shown by Savignano et al. [20] on 2D planes.
The force and moments are applied directly to the CRes of the left premolar of Patient 1. A: The CRots of the forces rotating from long axis towards the labiolingual direction of the tooth. B: The CRots of the forces rotating from the mesiodistal axis towards the long axis of the tooth. Note that each color corresponds to a different force direction, the M: F changes from -12N to 12N, and the distance between the CRots of the same color changes nonlinearly while the M: F values increment linearly.
Discussion
Scenario 1
As illustrated in Fig 4, the computed CRots for different patients follow similar trajectories for the teeth of the same types. It can also be seen that linear increments in the magnitude of the couples nonlinearly increases the distance between the neighboring CRots of each tooth. The same behavior was reported by Savignano et al. [20, 23]. Moreover, the CRot trajectories begin from the CRes of the tooth, and by increasing the M: F, the trajectories move towards the crowns of the teeth, and then, they change their directions to the back of the teeth. A schematic curve of the aforementioned path can be seen in part A of Fig 5.
The turning points can be found as outlying points on the top of each tooth with the highest distance from their neighboring CRots. These points are associated with specific M: F values where the direction of the tooth movement changes in the opposite direction. In other words, in the turning point, the generated moment by the couple counteracts the generated moment caused by the tipping load. In addition, it can be observed that, in some cases, the CRot trajectories are not well aligned, i.e., they seem to be rotated about the long axis of the tooth. This could be due to the registration errors in matching the teeth models of different patients to each other, most likely caused by the size and high geometry variations of teeth.
Scenario 2
The investigation of the effects of the direction of the loading conditions on the CRot trajectories in Fig 5 reveals that although rotating the couples about each axis results in similar CRot trajectories for different rotation angles, equally rotating the direction of the couple about the three principal axes results in different patterns for the CRot positions, especially in labiolingual direction. Besides, by increasing the magnitude of the couples in a specific direction, the distances between two neighboring CRots changes nonlinearly which is in line with the results of the state-of-the-art [20, 23], indicating that changing the M: F nonlinearly affects the position of the CRots.
As it can be observed in Fig 5, the CRot trajectories almost follow the same ∞-shaped curve by the rotation of the couple direction about the mesiodistal and long axes of the tooth with different angles (subfigure A and B). However, the sequence of the CRot points stops earlier when increasing the rotation angle. The situation in rotating the couple’s direction about the labiolingual axis of the tooth looks a bit different. That is to say, the CRot trajectories follow different ∞-shaped curves which are seemed to be rotated about the long axis.
Considering the obtained results in Fig 4 and subfigure C of Fig 5, one can deduce that any changes in the direction of force and couple can result in the misalignment of the CRot trajectories, and specifically, rotation of the CRots trajectories about the long axis of the tooth. In fact, the latter case shows how such systematic changes can affect the trajectory of the CRots, from which we can infer the reason for the misalignments that occurred in the former case. This behavior is somehow seen in a related work [20] where the CRot points scatter in the mesiodistal-long axes plane seems more incoherent than the other cases.
In this study, the direction of the load and couples are determined based on the normal of the tooth surface at the underlying region. As a result, the normal direction is not necessarily parallel to the labiolingual direction of the tooth. Therefore, in contrast to the literature where the force system is simplified in 2D planes by assuming force and couple directions parallel to the labiolingual direction, we investigate a more realistic clinical force system and study the influence of the couple’s direction in the tooth movement.
Scenario 3
This work benefits from using more realistic material models and scenarios, such as applying the forces on the tooth crown, for better modeling of the forces in clinical environments. A hyperelastic material model is used here to model the nonlinear behavior of the PDL tissue under orthodontic forces which allows us to model different scenarios without any prior assumptions for the range of the applied load. Furthermore, we propose using 3D space for representing the position of the CRot for the exerted force and moment to help with better analysis and understanding of the results in real-world applications. In contrast, Savignano et al. [20, 23] have represented the results in 2D planes and applied the loads and moments on the CRes, which is difficult to precisely be determined for different teeth. Besides, they have utilized a linear material model to the PDL tissue which confines the range of the applied load [48].
Future directions
Although an accurate geometry of the PDL can be obtained in vitro using micro-computed tomographies, these methods are known to be highly invasive and can only be applied to dead specimens. In contrast, reconstructing the PDL layer using CBCT scans obtained in vivo is a challenging process [49] due to the small width of the PDL layer compared to the commonly-used voxel size of the scans [1, 50]. Therefore, for a consistent analysis of different patients and following the literature, the PDL layers of this study are generated by extruding the teeth roots with a uniform thickness of 0.2 mm. Still, further investigation on the PDL is required to study, for example, the effect of nonuniform PDL geometry on the position of the computed CRots since different factors such as the patient’s age and the presence of periodontal disease can influence the shape and thickness of the PDL geometry [25], and consequently, the local stress maxima and tooth displacements [49]. In addition, the same material parameters are used in this study for the teeth, PDLs, and bones of different patients. It should however be noted that the initial and long-term teeth movements can vary from an individual to another based on differences in material properties, the density of the surrounding bone, and the rate of the bone remodeling process, apart from the geometrical differences.
In addition to an accurate geometry reconstruction, utilizing appropriate material models and parameters is important for accurately describing the mechanical behavior of the PDL. Although experimental studies investigating the biomechanical behavior of the PDL in vitro [51] and vivo [11, 52, 53] have indicated an anisotropic nonlinear viscoelastic behavior of the fibrous PDL [42], computational studies mostly consider simpler material models such as linear elastic, bilinear, and piecewise linear [43] models. The linear assumption confines the range of the applied load [48], preventing the simulation of realistic load systems in clinical environments. More complex material models such as hyperplastic and viscoelastic have been used in the literature [42, 54], and few studies [39] have considered collagen fibers in the PDL layer to simulate an anisotropic behavior of the PDL layer [55, 56]. However, more investigations are required to see the effectiveness of such complex material models applied to tooth movement modeling in terms of computational complexity and accuracy.
Conclusion
In this work, three patient-specific computational models of human mandibles, composed of full dentition, PDL layers, and jawbone, were utilized to assess the position of the CRots in different patients’ teeth against varying M: F ratios. The patients’ teeth were aligned for a consistent inter-patient analysis and FE simulations were performed under identical scenarios and boundary conditions for different teeth of each patient to obtain the CRots.
This work benefitted from using more realistic material models and scenarios for better modeling of the forces in clinical environments. A hyperelastic material model was used for the PDL tissue under orthodontic forces which allows for modeling different scenarios without any force constraints. The 3D space was used for representing the positions of the CRots for the exerted force and moment for better analysis and understanding of the results in real-world applications.
The influence of the couple’s magnitudes and directions on the positions of the CRots of the patients’ teeth were examined. The results show that the CRot trajectories could follow similar patterns in the corresponding teeth, but any changes in the direction of the force and couple could cause misalignments of the CRot trajectories that could be seen as rotations about the long axis of the tooth.
This work considered a more realistic force system in multi-patient tooth movement modeling by studying the CRot positions. The measured CRot position in relation to the tooth geometry can be used to infer the type of tooth movement, e.g., pure translation, uncontrolled tipping, and crown or root tipping. This, in turn, can be used in treatment planning software to assist clinicians to identify optimal forces required for achieving a desired patient-specific tooth movement.
Acknowledgments
The authors would like to thank Dr. Paolo Maria Cattaneo from Melbourne Dental School of the University of Melbourne for his valuable contribution in verifying the segmented CBCT scans used in this study.
References
- 1. Savignano R, Valentino R, Razionale A, Michelotti A, Barone S, D’Antò V. Biomechanical effects of different auxiliary-aligner designs for the extrusion of an upper central incisor: a finite element analysis. Journal of Healthcare Engineering. 2019;2019. pmid:31485303
- 2. Cervino G, Fiorillo L, Arzukanyan AV, Spagnuolo G, Campagna P, Cicciù M. Application of bioengineering devices for stress evaluation in dentistry: the last 10 years FEM parametric analysis of outcomes and current trends. Minerva Stomatologica. 2020;69(1):55–62. pmid:32214067
- 3. Lavorgna L, Cervino G, Fiorillo L, Di Leo G, Troiano G, Ortensi M, et al. Reliability of a virtual prosthodontic project realized through a 2D and 3D photographic acquisition: An experimental study on the accuracy of different digital systems. International Journal of Environmental Research and Public Health. 2019;16(24):5139.
- 4. Sambataro S, Cervino G, Fiorillo L, Cicciù M. Upper first premolar positioning evaluation for the stability of the dental occlusion: Anatomical considerations. Journal of Craniofacial Surgery. 2018;29(5):1366–1369. pmid:29521757
- 5.
Gholamalizadeh T, Darkner S, Cattaneo PM, Søndergaard P, Erleben K. Mandibular Teeth Movement Variations in Tipping Scenario: A Finite Element Study on Several Patients. arXiv preprint arXiv:201005258. 2020;.
- 6. Nyashin Y, Nyashin M, Osipenko M, Lokhov V, Dubinin A, Rammerstorfer F, et al. Centre of resistance and centre of rotation of a tooth: experimental determination, computer simulation and the effect of tissue nonlinearity. Computer Methods in Biomechanics and Biomedical Engineering. 2016;19(3):229–239. pmid:25761525
- 7. Meyer BN, Chen J, Katona TR. Does the center of resistance depend on the direction of tooth movement? American Journal of Orthodontics and Dentofacial Orthopedics. 2010;137(3):354–361. pmid:20197172
- 8. Viecilli RF, Budiman A, Burstone CJ. Axes of resistance for tooth movement: does the center of resistance exist in 3-dimensional space? American Journal of Orthodontics and Dentofacial Orthopedics. 2013;143(2):163–172. pmid:23374922
- 9. Sia SS, Koga Y, Yoshida N. Determining the center of resistance of maxillary anterior teeth subjected to retraction forces in sliding mechanics: an in vivo study. The Angle Orthodontist. 2007;77(6):999–1003. pmid:18004933
- 10. Sia S, Shibazaki T, Koga Y, Yoshida N. Experimental determination of optimal force system required for control of anterior tooth movement in sliding mechanics. American Journal of Orthodontics and Dentofacial Orthopedics. 2009;135(1):36–41. pmid:19121498
- 11. Yoshida N, Jost-Brinkmann PG, Koga Y, Mimaki N, Kobayashi K. Experimental evaluation of initial tooth displacement, center of resistance, and center of rotation under the influence of an orthodontic force. American Journal of Orthodontics and Dentofacial Orthopedics. 2001;120(2):190–197. pmid:11500662
- 12. Van Schepdael A, Geris L, Vander Sloten J. Analytical determination of stress patterns in the periodontal ligament during orthodontic tooth movement. Medical Engineering & Physics. 2013;35(3):403–410. pmid:23046973
- 13. Provatidis CG. An analytical model for stress analysis of a tooth in translation. International Journal of Engineering Science. 2001;39(12):1361–1381.
- 14. Oh MB, Mo SS, Hwang CJ, Chung C, Kang JM, Lee KJ. The 3-dimensional zone of the center of resistance of the mandibular posterior teeth segment. American Journal of Orthodontics and Dentofacial Orthopedics. 2019;156(3):365–374. pmid:31474266
- 15. Schneider J, Geiger M, Sander FG. Numerical experiments on long-time orthodontic tooth movement. American Journal of Orthodontics and Dentofacial Orthopedics. 2002;121(3):257–265. pmid:11941339
- 16. Schmidt F, Geiger ME, Jäger R, Lapatki BG. Comparison of methods to determine the centre of resistance of teeth. Computer Methods in Biomechanics and Biomedical Engineering. 2016;19(15):1673–1682. pmid:27113861
- 17.
Melsen B, Cattaneo PM, Dalstra M, Kraft DC. The importance of force levels in relation to tooth movement. In: Seminars in Orthodontics. Elsevier; 2007. p. 220–233.
- 18. Liao Z, Chen J, Li W, Darendeliler MA, Swain M, Li Q. Biomechanical investigation into the role of the periodontal ligament in optimising orthodontic force: a finite element case study. Archives of Oral Biology. 2016;66:98–107. pmid:26943815
- 19. Smith RJ, Burstone CJ. Mechanics of tooth movement. American Journal of Orthodontics. 1984;85(4):294–307. pmid:6585147
- 20. Savignano R, Viecilli RF, Paoli A, Razionale AV, Barone S. Nonlinear dependency of tooth movement on force system directions. American Journal of Orthodontics and Dentofacial Orthopedics. 2016;149(6):838–846. pmid:27241994
- 21. Burstone CJ. The biomechanics of tooth movement. Vistas in Orthodontics. 1962; p. 197–213.
- 22. Cattaneo PM, Dalstra M, Melsen B. Moment-to-force ratio, center of rotation, and force level: a finite element study predicting their interdependency for simulated orthodontic loading regimens. American Journal of Orthodontics and Dentofacial Orthopedics. 2008;133(5):681–689. pmid:18456141
- 23. Savignano R, Viecilli RF, Oyoyo U. Three-dimensional nonlinear prediction of tooth movement from the force system and root morphology. The Angle Orthodontist. 2020;90(6):811–822. pmid:33378522
- 24. Fedorov A, Beichel R, Kalpathy-Cramer J, Finet J, Fillion-Robin JC, Pujol S, et al. 3D Slicer as an image computing platform for the Quantitative Imaging Network. Magnetic Resonance Imaging. 2012;30(9):1323–1341. pmid:22770690
- 25. Li Y, Jacox LA, Little SH, Ko CC. Orthodontic tooth movement: The biology and clinical implications. The Kaohsiung Journal of Medical Sciences. 2018;34(4):207–214. pmid:29655409
- 26.
Barone S, Paoli A, Paolo N, Razionale AV, Giannese M. Mechanical and Geometrical properties assessment of thermoplastic materials for biomedical application. In: Advances on Mechanics, Design Engineering and Manufacturing. Springer; 2017. p. 437–446.
- 27. Seo JH, Eghan-Acquah E, Kim MS, Lee JH, Jeong YH, Jung TG, et al. Comparative Analysis of Stress in the Periodontal Ligament and Center of Rotation in the Tooth after Orthodontic Treatment Depending on Clear Aligner Thickness—Finite Element Analysis Study. Materials. 2021;14(2):324. pmid:33435457
- 28.
Schmidt R, Singh K. Meshmixer: an interface for rapid mesh composition. In: ACM SIGGRAPH 2010 Talks. ACM; 2010. p. 6.
- 29. Si H. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software (TOMS). 2015;41(2):1–36.
- 30.
Shewchuk J. What is a good linear finite element? interpolation, conditioning, anisotropy, and quality measures (preprint). University of California at Berkeley. 2002;73:137.
- 31. Liu A, Joe B. Relationship between tetrahedron shape measures. BIT Numerical Mathematics. 1994;34(2):268–287.
- 32. Misztal MK, Erleben K, Bargteil A, Fursund J, Christensen BB, Bærentzen JA, et al. Multiphase flow of immiscible fluids on unstructured moving meshes. IEEE Transactions on Visualization and Computer Graphics. 2013;20(1):4–16.
- 33.
Baker TJ. Element quality in tetrahedral meshes. In: 7th International Conference on Finite Element Methods in Flow Problems; 1989. p. 1018.
- 34. Freitag LA, Ollivier-Gooch C. Tetrahedral mesh improvement using swapping and smoothing. International Journal for Numerical Methods in Engineering. 1997;40(21):3979–4002.
- 35. Caendish JC, Field DA, Frey WH. An approach to automatic three-dimensional finite element mesh generation. International Journal for Numerical Methods in Engineering. 1985;21(2):329–347.
- 36. Maas SA, Ellis BJ, Ateshian GA, Weiss JA. FEBio: finite elements for biomechanics. Journal of Biomechanical Engineering. 2012;134(1). pmid:22482660
- 37. Żmudzki J, Chladek W. Stress present in bone surrounding dental implants in FEM model experiments. Journal of Achievements in Materials and Manufacturing Engineering. 2008;27(1):71–74.
- 38. Roda-Casanova V, Zubizarreta-Macho Á, Sanchez-Marin F, Alonso Ezpeleta Ó, Albaladejo Martínez A, Galparsoro Catalán A. Computerized Generation and Finite Element Stress Analysis of Endodontic Rotary Files. Applied Sciences. 2021;11(10):4329.
- 39. Ortún-Terrazas J, Cegoñino J, Santana-Penín U, Santana-Mora U, del Palomar AP. Approach towards the porous fibrous structure of the periodontal ligament using micro-computerized tomography and finite element analysis. Journal of the Mechanical Behavior of Biomedical Materials. 2018;79:135–149. pmid:29304428
- 40.
Maas S, Rawlins D, Weiss J, Ateshian G. FEBio 2.8 User Manual; 2018.
- 41. Ziegler A, Keilig L, Kawarizadeh A, Jäger A, Bourauel C. Numerical simulation of the biomechanical behaviour of multi-rooted teeth. The European Journal of Orthodontics. 2005;27(4):333–339. pmid:15961572
- 42. Qian L, Todo M, Morita Y, Matsushita Y, Koyano K. Deformation analysis of the periodontium considering the viscoelasticity of the periodontal ligament. Dental Materials. 2009;25(10):1285–1292. pmid:19560807
- 43. Cattaneo P, Dalstra M, Melsen B. The finite element method: A tool to study orthodontic tooth movement. Journal of Dental Research. 2005;84(5):428–433. pmid:15840778
- 44. Uhlir R, Mayo V, Lin PH, Chen S, Lee YT, Hershey G, et al. Biomechanical characterization of the periodontal ligament: Orthodontic tooth movement. The Angle Orthodontist. 2016;87(2):183–192. pmid:27542105
- 45. Chen Y, Medioni G. Object modelling by registration of multiple range images. Image and Vision Computing. 1992;10(3):145–155.
- 46. Besl PJ, McKay ND. A method for registration of 3-D shapes. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1992;14(2):239–256.
- 47. Myronenko A, Song X. Point set registration: Coherent point drift. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2010;32(12):2262–2275. pmid:20975122
- 48. Kawarizadeh A, Bourauel C, Zhang D, Götz W, Jäger A. Correlation of stress and strain profiles and the distribution of osteoclastic cells induced by orthodontic loading in rat. European Journal of Oral Sciences. 2004;112(2):140–147. pmid:15056111
- 49. Hohmann A, Kober C, Young P, Dorow C, Geiger M, Boryor A, et al. Influence of different modeling strategies for the periodontal ligament on finite element simulation results. American Journal of Orthodontics and Dentofacial Orthopedics. 2011;139(6):775–783. pmid:21640884
- 50. Dorow C, Schneider J, Sander FG. Finite element simulation of in vivo tooth mobility in comparison with experimental results. Journal of Mechanics in Medicine and Biology. 2003;3(01):79–94.
- 51. Sanctuary CS, Wiskott HA, Justiz J, Botsis J, Belser UC. In vitro time-dependent response of periodontal ligament to mechanical loading. Journal of Applied Physiology. 2005;99(6):2369–2378. pmid:16109833
- 52. Keilig L, Drolshagen M, Tran K, Hasan I, Reimann S, Deschner J, et al. In vivo measurements and numerical analysis of the biomechanical characteristics of the human periodontal ligament. Annals of Anatomy-Anatomischer Anzeiger. 2016;206:80–88. pmid:26395824
- 53. Yoshida N, Koga Y, Peng CL, Tanaka E, Kobayashi K. In vivo measurement of the elastic modulus of the human periodontal ligament. Medical Engineering & Physics. 2001;23(8):567–572. pmid:11719079
- 54. Natali AN, Carniel EL, Pavan PG, Sander FG, Dorow C, Geiger M. A visco-hyperelastic-damage constitutive model for the analysis of the biomechanical response of the periodontal ligament. Journal of Biomechanical Engineering. 2008;130(3). pmid:18532853
- 55. Fill TS, Toogood RW, Major PW, Carey JP. Analytically determined mechanical properties of, and models for the periodontal ligament: critical review of literature. Journal of Biomechanics. 2012;45(1):9–16. pmid:22014328
- 56. Komatsu K. Mechanical strength and viscoelastic response of the periodontal ligament in relation to structure. Journal of Dental Biomechanics. 2010;2010. pmid:20948569