Introduction
Although many attempts have been made to simulate orthodontic tooth movement using the finite element method, most were limited to analyses of the initial displacement in the periodontal ligament and were insufficient to evaluate the effect of orthodontic appliances on long-term tooth movement. Numeric simulation of long-term tooth movement was performed in some studies; however, neither the play between the brackets and archwire nor the interproximal contact forces were considered. The objectives of this study were to simulate long-term orthodontic tooth movement with the edgewise appliance by incorporating those contact conditions into the finite element model and to determine the force system when the space is closed with sliding mechanics.
Methods
We constructed a 3-dimensional model of maxillary dentition with 0.022-in brackets and 0.019 × 0.025-in archwire. Forces of 100 cN simulating sliding mechanics were applied. The simulation was accomplished on the assumption that bone remodeling correlates with the initial tooth displacement.
Results
This method could successfully represent the changes in the moment-to-force ratio: the tooth movement pattern during space closure.
Conclusions
We developed a novel method that could simulate the long-term orthodontic tooth movement and accurately determine the force system in the course of time by incorporating contact boundary conditions into finite element analysis. It was also suggested that friction is progressively increased during space closure in sliding mechanics.
Highlights
- •
We analyzed long-term tooth movement using contact conditions.
- •
A test bench to evaluate temporal change of M/F ratio was constructed in silico.
- •
We examined changes in the force system during space closure in sliding mechanics.
- •
Interaction between brackets and archwire must be modeled accurately for analysis.
- •
Friction at posterior teeth is progressively increased in sliding mechanics.
The prediction and planning of orthodontic tooth movement have largely depended on clinical experiences. If the long-term tooth movement can be accurately predicted, the treatment results will be greatly improved. With such a background, demand for simulating tooth movement under various loading conditions has been increasing in recent years. In orthodontic treatment, 2 phases of tooth movement are reiterated at each visit of a patient to the clinic. In the first phase, the application of orthodontic force induces the instantaneous deformation of the periodontal ligament (PDL), called “initial displacement.” Then, in the second phase, the stress distributed in the PDL initiates bone remodeling. Thus, long-term orthodontic tooth movement occurs as a result of biomechanical and histologic changes in the periodontal tissues produced by the reiterations of the 2 phases consisting of mechanical and biologic processes.
Many finite element (FE) studies have been performed to simulate tooth movement to improve therapeutic efficiency in orthodontic treatment. Most were, however, limited to the analysis of initial displacement. Although the initial displacement may be an indicator of the subsequent tooth movement after some bone remodeling, loading conditions change during tooth movement. In other words, as a tooth is tipped, the relationship between the line of action of the force and the location of the center of resistance (CR) of the tooth—the force system—is changing. It is therefore difficult to precisely predict overall tooth movement from the initial displacement. Incorporation of the bone remodeling process into the FE method is thus required to simulate the long-term tooth movement more accurately.
Only a few studies investigating bone remodeling with the FE method have been performed. Bourauel et al and Schneider et al simulated orthodontic tooth movement for a simplified single-root tooth using the FE method based on bone remodeling theories. In these simulations, alveolar bone geometries must be updated each time small amounts of bone remodeling are calculated; thus, such methods are time-consuming when trying to simulate the long-term movements of multiple teeth simultaneously. Kojima and Fukui and Kojima et al performed numeric analyses to simulate long-term tooth movements of multiple teeth by using special elements that represented teeth supported by PDLs in an FE model. In those studies, however, neither the play between the brackets and archwire nor the interproximal contact force was considered. A previous study showed that the play between the brackets and archwire has a large impact on the tooth movement pattern. Moreover, the contact forces between teeth are thought to substantially affect the force system acting on the tooth and its resultant movement. To our knowledge, no previous study has considered both the bracket-wire interface and the interproximal contact force in the numeric analysis of orthodontic tooth movement. An analysis incorporating these 2 mechanical conditions will determine the force system in sliding mechanics, which has been mechanically indeterminate when using conventional methods.
The purposes of this study were to develop a method to simulate long-term orthodontic tooth movement with contact boundary conditions based on the FE method and to determine the force system acting on each tooth in sliding mechanics using a model with realistic bracket slot, archwire, and tooth dimensions. See Supplemental Materials for a short video presentation about this study.
Material and methods
Three-dimensional (3D) images of a maxillary dentition were obtained from a dry skull using a multi-image microcomputed tomography scanner (3DX; J. Morita, Kyoto, Japan) with a voxel size of 80 μm. The data were exported as a DICOM file to 3D image processing and editing software (Mimics 10.02; Materialise, Leuven, Belgium). Then the DICOM data were reconstructed to 3D surface data and imported to FE analysis preprocessing and postprocessing software (Patran 2012.1; MSC Software, Los Angeles, Calif) as 3-node triangular elements. These triangular elements were converted to 4-node quadrilateral elements. After remeshing, the teeth were aligned to the ideal positions to make a normal dentition. Each interproximal space between adjacent teeth was set to be less than 0.01 mm.
The PDL was constructed on the root surface with 8-node hexahedral elements. The thickness was determined to be a uniform 0.2 mm, because a previous study showed that whether the PDL is uniform or not does not make a significant difference in the strain distribution in the PDL.
The coordinate system of the analysis model was defined so that the x-axis pointed in the transverse direction, the y-axis in the anteroposterior direction, and the z-axis in the vertical direction, with the positive direction to the left, posterior, and up, respectively ( Fig 1 ). All analyses were performed using an FE package (Marc 2014.1; MSC Software).
The PDL was modeled as a nonlinear (bilinear) isotropic material. An elastic modulus from 0.03 to 0.3 MPa and a Poisson’s ratio of 0.3 were assigned to the PDL according to previous studies. To develop a more simplified model for reducing the analysis time, the tooth material was assumed to be a rigid body, since an orthodontic force is small, thereby causing only a small amount of tooth deformation, which can be negligible. Therefore, we modeled the teeth using thin shell elements with a thickness of 3.0 mm and a Young’s modulus of 204 GPa, which is 10 times higher than the previously reported value. A Poisson’s ratio of 0.3 was also assigned. In a preliminary analysis, we verified that this tooth model had enough stiffness.
Figure 2 shows an algorithm for simulating long-term orthodontic tooth movement. As the first step, an orthodontic force is applied to the tooth model. Then the initial displacement is produced reflecting the deformed PDL. When we take measurements of the initial tooth displacement, the alveolar bone is assumed to be a rigid body, since its deformation is negligible compared with that of the PDL. According to this assumption, displacements are constrained at the outer surface of the PDL. This procedure omits the need to model the alveolar bone. Consequently, the analysis time for calculating bone deformation and remeshing the solid elements of the deformed alveolar bone, which had been done in previous studies, was eliminated.
As the second step, each node forming the outer surface of the PDL is displaced so that the PDL is restored to its original configuration and thickness of 0.2 mm. In other words, the geometry of the PDL is updated on the assumption that the orthodontic tooth movement is correlated with the initial displacement, since that assumption has been accepted by the widely recognized theory of orthodontic mechanics. Thus, the alveolar bone is remodeled so that the space between the tooth and the alveolus is kept constant during the long-term movement. In this manner, the alveolar bone remodeling is simulated after the first application of an orthodontic force.
Subsequently, the second orthodontic force is applied to the tooth model. At this time, displacements are constrained to the outer surface of the PDL whose geometry is updated through the process of the first remodeling, and initial displacement is again produced. Then, these 2 steps—the initial displacement and bone remodeling—are iterated to simulate orthodontic tooth movement after the long-term process of absorption and apposition of the alveolar bone. These procedures are executed automatically using custom-developed subroutines.
Contact boundary conditions, which used the node-to-segment function of the Marc software, were applied to each interproximal surface of the tooth. A contact surface was formed on the midsagittal plane for symmetric analysis ( Fig 3 ).
The method for evaluating 3D tooth movement is described as follows. Three-dimensional motion can be described as a combination of translation and rotational motion around an arbitrary point. The center of gravity is generally used as the center of rotation to describe motions in rigid body dynamics. Instead of the center of gravity, the CR is commonly applied in the mechanical theory of orthodontic tooth movement. Thus, the tooth movements were quantified by measuring the rotation of each tooth body and the translation of the CR. The rotation and translation were calculated from the entries of the transformation matrix between the initial position and the resulting position. However, previous studies have shown that the CR does not exist as a point in 3D space, and the direction of the applied force that translates a tooth does not necessarily correspond to that of a moving CR because of the geometric asymmetry of a tooth. Therefore, the CR is determined expediently using the following method ( Fig 4 ).
- 1.
Three arbitrary nodes on each tooth are chosen, and displacement boundary conditions are prescribed; ie, displacement U i is assigned to each node in the distal direction. Let U i be defined as U i = (Lsinθ i , Lcosθ i , 0), where L is 0.02 mm, and θ i is the angle between U i and the x-axis.
- 2.
Each reaction force acting on the 3 nodes that are subjected to displacement is calculated.
- 3.
A single counter force, which is balanced with a resultant force of the 3 force vectors, is determined. This force translates the tooth model in the U i direction. We confirmed that the tooth moves bodily on applying the force to the model.
- 4.
Another displacement U‘i=(Lsin(θi+90o),Lcos(θi+90o),0)
U i ‘ = ( Lsin ( θ i + 90 o ) , Lcos ( θ i + 90 o ) , 0 )
is assigned to the 3 nodes perpendicular to U i in the transverse direction. - 5.
Each reaction force acting on the 3 nodes is calculated. Then, a single counter force that translates the tooth model in the U‘i
U i ‘
direction is calculated in the same manner as 2 and 3. - 6.
The closest point between the lines of action of the 2 forces, which translate the tooth model in the U i and U‘i
U i ‘
directions, respectively, is calculated. Let CR i (center of resistance) be the point determined. - 7.
Each position of CR i is calculated on applying forces running in several different directions parallel to the occlusal plane, since it varies depending on the direction of U i . The ultimate position of the CR of the tooth model is determined as the point in which each coordinate value of CR i is averaged. Thus, the position of the CR is represented by the following equation: CR=∑360i=1CRi/360
CR = ∑ i = 1 360 C R i / 360
