A Finite Element Approach for Locating the Center of Resistance of Maxillary Teeth

The center of resistance (C_RES) of a tooth or segment of teeth is analogous to the center of mass of a free body. It is a term borrowed from the field of mechanics of rigid bodies. When a single force is applied at the C_RES, translation of the tooth in the direction of the line of action of the force occurs^1,2. The position of the C_RES depends not only on the tooth's anatomy and properties but also on its environment (e.g., periodontal ligament, surrounding bone, adjacent teeth). The tooth is a restrained body, making its C_RES similar to the center of mass of a free body. In the manipulation of appliances, most orthodontists consider the relationship of the force vector to the C_RES of a tooth or a group of teeth. Indeed, whether an object will display tipping or bodily movement when submitted to a single force is mainly determined by the location of the C_RES of the object and the distance between the force vector and the C_RES. If this can be accurately predicted, treatment results will be greatly improved. Thus, an accurate estimation of C_RES can greatly enhance the efficiency of orthodontic tooth movement.

For decades, the orthodontic field has been revisiting research regarding the location of the C_RES of a given tooth, segment, or arch^{1,2,3,4,5,6,7,8,9,10,11,12}. However, these studies have been limited in their approach in many ways. Most studies have determined the C_RES for only a few teeth, leaving out the majority. For example, the maxillary central incisor and the maxillary incisor segment have been evaluated quite extensively. On the other hand, there are only a few studies on the maxillary canine and first molar and none for the remaining teeth. Also, many of these studies have determined the location of the C_RES based on generic anatomical data for teeth, measurements from two-dimensional (2D) radiographs, and calculations on 2D drawings⁸. In addition, some of the current literature uses generic models or three-dimensional (3D) scans of dentiform models rather than human data^4,8. As orthodontics shifts into 3D technology for planning tooth movement, it is crucial to revisit this concept to develop a 3D, scientific understanding of tooth movement.

With technological advancements resulting in increased computational power and modeling capabilities, the ability to create and study more complex models has increased. The introduction of computed tomography scanning and cone-beam computed tomography (CBCT) scanning has thrust models and calculations from the 2D world into 3D. Simultaneous increases in computing power and software complexity have allowed researchers to use 3D radiographs to extract accurate anatomical models for use in advanced software to segment the teeth, bone, periodontal ligament (PDL), and various other structures^{7,8,9,10,13,14,15}. These segmented structures can be converted into a virtual mesh for use in engineering software to calculate the response of a system when a given force or displacement is applied to it.

This study proposes a specific, replicable methodology that can be utilized to examine hypothetical orthodontic force systems applied on models derived from CBCT images of live patients. In utilizing this methodology, investigators can then estimate the C_RES of various teeth and take into consideration the biological morphology of dental structures, such as tooth anatomy, number of roots and their orientation in 3D space, mass distribution, and structure of periodontal attachments. A general outline of this process is shown in Figure 1. This is to orient the reader to the logical process involved in generation of 3D tooth models for locating the C_RES.

An institutional review board exemption was obtained for evaluating CBCT volumes archived in the Division of Oral and Maxillofacial Radiology (IRB No. 17-071S-2).

1. Volume selection and criteria

1. Acquire a CBCT image of the head and face¹⁶.
2. Examine the image for tooth alignment, missing teeth, voxel size, field of view, and overall quality of the image.
3. Make sure the voxel size is not larger than 350 µm (0.35 mm).

2. Segmentation of the teeth and bone

1. Load the raw DICOM files of the CBCT image into Mimics software for segmentation (Figure 2). Click Image > Crop Project. Crop the image to include only the maxilla and maxillary teeth.
   NOTE: The field of view should be large enough to capture the maxilla and maxillary teeth. Make sure the image includes the tooth crowns, the hard palate up to the nasal floor, maxillary sinuses, facial surfaces of the maxillary teeth, and the posterior extent of the hard palate and maxillary tuberosity.
2. Right click on the tab for Mask and create a New Mask for the image. Rename the mask as UL1, UL2, …, UL7 for the left side and UR1, UR2, …, UR7 for the right side, based on the tooth of interest.
3. Identify the tooth of interest on the masked CBCT image (see views). Use the Clear Mask tool to erase the mask. The software might be unable to distinguish between the teeth and bone because the gray values of the two are similar.
   NOTE: The thresholding tool in Mimics is unable to segment the teeth and bone separately. Therefore, a different method for segmentation is required.
4. Click on the Multiple Slice Edit tool (Ctrl + M). Select the view (Axial, Coronal, or Sagittal). Manually highlight (i.e., draw) some of the slices as deemed necessary.
   NOTE: Highlighting more slices adds more detail to the structure.
5. Click the Interpolate tool to fill up the volume for the skipped slices and apply.
6. Generate the 3D volume for the tooth by right-clicking on the mask and selecting the option to calculate the 3D volume.
7. Repeat steps 2.2-2.6 for each tooth of the maxillary arch.
8. Select all the 3D maxillary teeth, UL7-UR7. Right-click to select Smoothing. Set the smoothing factor to 0.4 and iterations to 4.
9. To segment the maxillary bones right-click on the tab for Mask. Create a new mask for the image.
10. From the drop-down menu for predefined threshold sets, select Custom. Adjust the threshold value to include the complete maxillary bone. Be sure to check the Fill Holes box before applying the threshold.
    NOTE: Small holes of ≤1 mm in the cortical bone are acceptable, because they can be removed easily in later stages.
11. Click on the Dynamic Region Growing Tool to fill the large holes visible in the mask. Select the maxillary bone mask as the target for the tool in addition to selecting the Multiple Layer box. Use 50 for the Min and 150 for Max values. Hold down the Control key while clicking on the areas of cortical bone that were not highlighted in the mask.
12. Right-click on the maxillary bone mask for the Smooth Mask function. Repeat this step 3x for best results.
13. Generate the 3D volume for the maxilla by right-clicking on the mask and selecting the option to calculate the 3D volume.
14. Select the 3D maxillary bone. Right-click to select smoothing. Set the smoothing factor to ~0.4 and iterations to 4.
15. Select the 3D maxillary bone and right-click to select Wrap. Set 0.2 mm for the smallest detail and 1 mm for the gap closing distance. Check the Protect Thin Walls option. Press Ok.
16. Rename the 3D maxillary bone "Maxilla".

3. Cleaning and meshing

1. Select the 3D objects and copy (Ctrl + C).
2. Open the 3matic software, and paste (Ctrl + V) the selected 3D objects. They will appear in the object tree and work area of 3matic as a 3D structure (Figure 3).
3. Click the Fix tab from the tool bar and use the Smooth option. Under the Operations box select the desired 3D object(s) or entities and Apply the default parameters.
4. Click the Finish tab from the tool bar and use the Local Smoothing option. Under the Operations Box select the desired 3D object(s) or entities. Use the cursor to manually smoothen out the desired regions.
5. Duplicate the teeth. On the object tree select all teeth, right-click, and select Duplicate.
6. Select All Duplicated Teeth, group, and name the folder "group 1". The original set will serve as the final teeth for analysis.
7. For the duplicated teeth in group 1, click the Curve Module and the Create Curve option. Manually draw a curve around the cementoenamel junction (CEJ) for all duplicated teeth.
8. Select the Curve, Contour, and Border entities under the Smooth Curve option.
9. Separate the crown and root surfaces into their own parts by selecting the Split Surfaces by Curves option and left-clicking on the 3D object to select.
10. Generate the PDL from the root structure of the tooth by splitting the tooth into root and crown at the CEJ.
    1. Duplicate the 3D objects from group 1 (generated in step 3.6) as group 2. For group 2, in the object tree box, click on the Object. From the surface list delete the crown surface. Perform this step for all objects in group 2.
    2. For group 2, click on Design Module > Hollow. Apply the desired parameters (Table 1).
    3. Click on the Fix Module > Fix Wizard. Click on the individual parts, update, and follow the given directions.
    4. Repeat step 3.10.3 for all parts. Rename all parts in the group 2 as "UL1_PDL" to "UL7_PDL" and "UR1_PDL" to "UR7_PDL".
11. In group 1, from the object tree box, click on Object. From the surface list delete the root surface.
12. Select Fill Hole Normal option and select the contour. Click on Bad Contour and Apply. The entire space will be filled.
13. Select the Design Module > Local Offset and select the entire crown surface. Check the following options: Direction (select external), Offset Distance (select 0.5), and Diminishing Distance (select 2.0). Apply.
14. Repeat step 3.13.
15. Repeat steps 3.11-3.14 for each tooth of the maxillary arch.
16. Remesh (Figure 3)
    1. Click the Remesh Module > Create Non-Manifold Assembly > Main Entity > Maxilla from the object tree. Select intersecting entity for all objects from 3.4 (original teeth) and select Apply.
    2. Click the Remesh Module. Split the non-manifold assembly.
    3. Repeat steps 3.16.1-3.16.2 using an intersecting entity as all objects from group 1 and Apply.
    4. As an optional step, only if required, select the Finish Module > Trim > Entity > Maxilla. Select the excess structure (i.e., noise) and Apply.
    5. Click the Fix Module > Fix Wizard > Maxilla > Update. Follow the directions given.
    6. Repeat step 3.16.1 using an intersecting entity as all objects from group 2 and Apply.
    7. Click the Remesh Module > Adaptive Remesh. Select all intersecting entities from 3.16.6 and Apply.
    8. Click the Remesh Module > Split Non-manifold Assembly.
    9. Click the Remesh Module > Create Non-manifold Assembly > Main Entity > Individual Object (PDL) from group 2 from the object tree. Select Intersecting Entity > Select Respective Object from step 3.4 (corresponding to that tooth type) and Apply.
    10. Click Remesh Module > Adaptive Remesh. Select the intersecting entity from 3.16.9 and Apply.
    11. Click Remesh Module > Split Non-manifold Assembly.
    12. Repeat steps 3.16.9-3.16.11 for each tooth.
17. Click the Remesh Module > Quality Preserving Reduce Triangles. In the object tree select all entities (i.e., teeth, PDLs, and Maxilla) and Apply.
18. Click Remesh Module > Create Volume Mesh > Select Entity. Choose Mesh Parameters.
19. Repeat step 3.18 for all entities (i.e., teeth, PDLs, and Maxilla).
20. Manually export the input(.inp) files from 3Matic to Abaqus (Figure 4).

4. Finite element analysis

NOTE: All custom Python scripts can be found in the supplemental attachments. They have been generated using the macro manager function in Abaqus.

1. Preprocessing setup
   1. Open Abaqus and select Standard Model. Click File > Set the Work Directory > Select Location for File Storage.
   2. Click File > Run Script and select Model_setup_Part1.py
   3. In the Model Directory specify the file path to load .inp files on Abaqus.
   4. Click on Models > Simulation > Parts > Maxilla > Surfaces.
   5. Name the surface in the dialog box "UL1 _socket".
   6. Under Select the Region of the Surface choose By Angle. Add "15" as the angle.
   7. Make sure all areas of the socket are selected. Press Done when completed.
   8. Repeat steps 4.1.4-4.1.7 for the individual sockets.
   9. Click on Models > Simulation > Parts. Then select UL1 > Surfaces. Name the surface "UL1".
   10. Under Select the Region of the Surface opt for "Individually". Select the tooth on the screen and press Done.
   11. Repeat steps 4.1.9-4.1.10 for all teeth.
   12. Click on Models > Simulation > Parts. Then select UL1_PDL > Surfaces. Name the surface "UL1_PDL_inner".
   13. Under Select the Region of the Surface choose By Angle. Add "15" as the angle.
       NOTE: If an error is found during the final simulation, reduce the angle and reselect the surface.
   14. Make sure the entire inner surface area of the PDL is selected. Press Done when completed.
   15. Select UL1_PDL > Surfaces. Name the surface "UL1_PDL_outer".
   16. Under Select the Region of the Surface choose By Angle. Add "15" as the angle.
       NOTE: If an error is found during the final simulation, reduce the angle and reselect the surface.
   17. Make sure the entire outer surface area of the PDL is selected. Press Done when completed.
   18. Repeat steps 4.1.13-4.1.19 for all PDLs.
   19. Click on File > Run Script and select Model_setup_Part2.py
   20. Click on Models > Simulation > BCs. Name BC_all, then select Step as Initial. Under category, select "Mechanical", and under "Types of Selected Step" select "Displacement/Rotation". Press Continue.
   21. Under Select Regions for the Boundary Condition select By Angle. Add "15" as the angle. Check Create Set. Select individual sockets for the 14 teeth. Press Done.
       NOTE: This helped simulate instantaneous tooth movement.
   22. Click on Models > Simulation > Assembly > Sets > Create Set. Name the set "U1_y_force".
   23. In Select the Nodes for the Set choose Individually.
       NOTE: A one Newton concentrated force was applied on a randomly selected tooth node in either the positive Y direction (simulating a distalization force) or the positive Z direction (simulating an intrusive force).
   24. Select a node at the center of the crown on the buccal surface of the upper central incisor (U1) and press Done.
   25. Click Sets > Create Set. Name the set "U1_z_force".
   26. Repeat steps 4.1.23-4.1.24.
   27. Repeat steps 4.1.22-4.1.26 for all teeth.
       NOTE: Before a set is generated for a particular tooth as in 4.1.25, go to Instance > Resume for that tooth.
2. Model set up
   1. Click on Models > Simulation > Assembly > Instances. Select All Instances and click Resume.
   2. Click on Tools > Query > Point/Node. Select a node at the center of a randomly selected central incisor and press Done.
   3. Under the command center at the bottom of the page, copy the X, Y, and Z coordinates of the node selected in step 4.2.2.
   4. Under the vertical tool bar select Translate Instance and select the entire assembly (i.e., all instances) on the screen. Press Done.
   5. In the Select a Start Point for the Translation Vector box, paste the copied coordinates in step 4.2.3 or enter the X, Y, and Z values. Click Enter.
   6. Under Select an End Point for the Translation Vector or enter X,Y,Z: enter the coordinates "0.0", "0.0", and "0.0". Click Enter.
   7. For Position of Instance, press Ok.
   8. Click on Tools > Query > Point/Node and select a node directly above the midline of the center incisors. Enter Done.
   9. Under the command center at the bottom of the page, copy the X, Y, and Z coordinates of the node selected in step 4.2.8.
   10. Under the vertical tool bar select Translate Instance and select the entire assembly (i.e., all instances) on the screen. Enter Done.
   11. Paste the copied coordinates into the Select a Start Point for the Translation Vector - or Enter X,Y,Z box. Click Enter.
   12. Under Select an End Point for the Translation Vector - or enter X,Y,Z: insert the coordinates as copied in step 4.2.9. Change the X coordinate to 0.0. Click Enter.
   13. For Position of Instance, press Ok.
   14. Click on File > Run Script and select Model_setup_Part3.py. Insert or change material properties.
   15. Click on Models > Simulation > Materials and click Bone/PDL/Tooth. Insert tissue specific properties.
   16. Click on File > Run Script and select Functions.py.
3. Processing the model
   1. Click on File > Run Script and select Job_submission.py.
      NOTE: The job module is where the user sets up one or more actions on the model, and job manager is where model analysis is started, progress is shown, and completion is noted.
   2. In the dialog box titled Suppress All, enter the sides (L or R) of the teeth based on constraints (Under Models > Simulation > Constraints). Press Ok.
   3. In the dialog box titled Job Submission enter "Y" to run the analysis for the specified tooth/teeth. Press Ok.
   4. In the dialog box titled Directions for Analysis enter "Y" to specify force application. Press Ok.
4. Post processing for C _RES estimation
   1. Choose File > Run Script > Bulk_process.py.
   2. In the dialog box titled Analyze Multiple Jobs enter "Y" for the specified tooth/teeth. Press Ok.
   3. In the dialog box titled Directions for Analysis enter "Y" for specifying force application. Press Ok.
   4. In the dialog box titled Get Input enter Specific Tooth Number as outlined named Instances (e.g., UL1 or UL5, etc.). Press Ok.
   5. Check coordinates for the Force About Point and Estimated Location in the command box. If they are not similar, then repeat steps 4.3.1-4.4.4.
      NOTE: After the jobs for each step were run, a user-defined algorithm created in Python was run within the Abaqus interface to analyze the reaction force system and subsequent moments created as a result of load application. The algorithm automatically suggests a new node location to apply the load such that a moment of near-zero magnitude is created within the force system. This proceeds in an iterative process, until the node location that creates a moment closest to zero when a force is applied through it is found or estimated. The algorithm is described in detail in the Discussion section.

In order to verify segmentation and manual outlining as described in the Procedures section (step 2), a maxillary first molar was extracted from a dry skull, and a CBCT image was taken. The image processing and editing software Mimics was used for manually outlining the tooth as described in step 2. Subsequently, meshing was performed, the segmented models were cleaned with 3matic software, and they were imported to Abaqus for analysis. We did not find any significant difference in the linear and volumetric measurements made on the FE model of the tooth and the actual tooth measured in the lab (Supplemental Document 4).

To verify the validity of the user-defined algorithm in determining the C_RES of an object, a simplified model of a beam encased within a sheath was used in the initial stages of script creation (Figure 5A). The steel encasement was constrained to three degrees of freedom of displacement, and the nodes at the beam/sheath interface were tied together. Nodes for force application were selected at random, and the subroutine was applied in an iterative fashion until the solution converged. In the simplified model, a length of 30 units and a width of 10 units were encased in the sheath. By following the defined algorithm and its calculations, the C_RES of the model beam was predicted (Figure 5B). This agreed with the theoretical calculations (see Supplemental Document 3). Thus, the validity of the user-defined algorithm was developed and verified in this simplified model and was subsequently implemented for the determination of the C_RES of maxillary teeth.

Table 2 shows the material properties assigned to the structures. Differences in the modelling of the material properties of the PDL and bone could affect the final location of the C_RES of a tooth. PDL anisotropy related to fiber orientation, differences in Poisson's ratio, loading patterns, and magnitude can also make a difference. PDL was assigned nonlinear, hyperelastic properties according to the Ogden model (µ1 = 0.07277, α1 = 16.95703, D1 = 3 x 10^-7)^22,23. Specific densities were also assigned = 1.85 g/cm³ for bone; 2.02 g/cm³ for teeth; and 1 g/cm³ for PDL (i.e., the density of water, because PDL is mostly comprised of water)^24,25.

To standardize the force vectors and locate the position of the C_RES, a cartesian coordinate system was constructed (X-Y-Z) and defined by the following orientations: Y-axis (anteroposterior or labiolingual axis) oriented along the midpalatal suture with the posterior portion in the positive direction, Z-axis in the vertical direction (superio-inferior or occluso-gingival axis) with the superior or gingival portion of the model in the positive direction, and the X-axis in the transverse direction (buccolingual axis) with the buccal portion in the positive direction (Figure 6).

This coordinate system was applied in two ways: 1) A global coordinate system was established with its origin (O) located between the facial surfaces of the central incisors below the incisive papilla located on a line bisecting the inter-incisor and inter-molar widths in the X-Y plane; 2) Local coordinate systems were constructed with an origin 'R' for every tooth. The 'R' point specific for every tooth was defined as the geometric center on the buccal surface of the crown. This site was chosen to approximate the closest location where an operator might place a bracket to apply orthodontic forces. Representative results are shown in Figure 7.

The C_RES located relative to the global and local coordinate systems are shown in Table 3 and Table 4. The locations of the C_RES obtained along the X coordinate when a force system was applied along the Y and Z coordinates were different from each other (Table 5). However, the average difference was small (0.88 ± 0.54 mm).

Figure 1: Design flow chart. The three-step workflow for locating C_RES. Please click here to view a larger version of this figure.

Figure 2: Layout of the Mimics software displaying maxillary teeth in all the three views (X-Y-Z) and as a volumetric model. Please click here to view a larger version of this figure.

Figure 3: Steps involved for generating periodontal ligament (PDL) utilizing non-manifold assembly of the 3matic software. Remesh Module (A) Create non-manifold assembly, (B) Maxilla is set as the main entity, (C) PDL is set as the intersecting entity, (D) Adaptive Remesh, (E) Splitting the maxilla and the PDL, (F) Follow steps B-F for PDL as the main entity and the selected tooth as an intersecting entity, (G) Create volume mesh. Please click here to view a larger version of this figure.

Figure 4: The Abaqus software layout. Please click here to view a larger version of this figure.

Figure 5: Simplified model of the steel beam. (A) The beam encased in a steel sheath used to test the accuracy of the defined algorithm. (B) Location of C_RES of the encased beam as predicted by the defined algorithms. Please click here to view a larger version of this figure.

Figure 6: The coordinate system for C_RES estimation relative to a global point of origin (O) and local point of origin (R) for each tooth. This is an illustration for the maxillary second premolar. This method was utilized for every tooth in the arch. Please click here to view a larger version of this figure.

Figure 7: Three-dimensional representation of the C_RES of maxillary teeth. (A) Central incisor. (B) Lateral incisor. (C) Canine. (D) First premolar. (E) Second premolar. (F) First molar. (G) Second molar. Please click here to view a larger version of this figure.

Table 1: Hollow tool parameters.

Table 2: Material properties of the finite element model.

Table 3: Three-dimensional (X-Y-Z) location of the C_RES of maxillary teeth in relation to the global point O.

Table 4: Three-dimensional (X-Y-Z) location of the C_RES of maxillary teeth in relation to a local point R for each tooth whose C_RES is being evaluated. Here, R is the geometric center of the buccal surface of the crown.

Table 5: Variation in the center of resistance position along the X-axis when force is applied along the Y-(Fy) and Z (Fz)-axes.

Supplemental Document 1: Python scripts of the algorithms utilized for the FEA. Please click here to view this file (Right click to download).

Supplemental Document 2: An overview of the analysis of the force system. Please click here to view this file (Right click to download).

Supplemental Document 3: Theoretical estimation of the center of mass of a simple beam encased in a sheath. Please click here to view this file (Right click to download).

Supplemental Document 4: A finite element model of an extracted maxillary first molar. Please click here to view this file (Right click to download).

This study shows a set of tools to establish a consistent workflow for finite element analysis (FEA) of models of maxillary teeth derived from CBCT images of patients to determine their C_RES. For the clinician, a clear and straightforward map of the C_RES of the maxillary teeth would be an invaluable clinical tool to plan tooth movements and predict side effects. The finite element method (FEM) was introduced in dental biomechanical research in 1973¹⁷, and since then has been applied to analyze the stress and strain fields in the alveolar support structures^{6,7,8,9,10,11,12}. As evidenced by the number of steps outlined in the workflow (Figure 1), creating finite element models is a complex task. Therefore, certain aspects of the methodology had to be simplified.

First, tooth movement only in the alveolar socket was considered by assuming that resorption and apposition of the alveolar bone did not occur. This type of displacement is called primary⁴ or instantaneous tooth movement¹⁸. It has been observed that the PDL is a critical entity in instantaneous tooth displacement. Bone and teeth could be reasonably assumed to be rigid to define PDL stresses for tooth movement¹⁵. Therefore, for this study the stress distribution was constrained within the tooth socket. The Create Boundary Condition tool allows the user to set boundary conditions for the model or apply constraints. Selected points are assigned zero degrees of freedom to ensure that the model stays rigid in that area. Consequently, the analysis time for calculating bone deformation and remeshing the solid elements of the deformed alveolar bone done in previous studies, was eliminated^19,20.

Second, an attempt to keep image resolution at moderate levels was made. CBCT image voxel size was 0.27 mm. This not only kept the radiation dosage at a minimum but also reduced the computational burden for assembling the global stiffness matrix for tetrahedral elements. However, the downside was that the CBCT resolution was insufficient to accurately and distinctly capture the PDL on the scans. This was largely because the average PDL thickness is around 0.15 mm-0.38 mm (average: 0.2 mm)²¹ and the image voxel size was 0.27 mm. This shortcoming with the CBCT scans created two issues: 1) The PDL could not be segmented on its own; and 2) Segmenting the bone and teeth using thresholding was not possible due to the lack of a distinct gray value change between the two. As a result, the software was unable to distinguish between the teeth and bone because the gray values were similar. In other words, Mimics was unable to segment the teeth and bone separately. Therefore, a different method of segmentation was developed. After attempting numerous tools, such as the region growing or split tool in Mimics, it was determined that the best way to segment the teeth was by manually highlighting the tooth structure on each slice of the CBCT. Here the Multiple Slice Edit Tool offered an efficiency advantage. Instead of having to manually highlight every slice, the user only has to highlight some of the slices. For this reason, it was the best method for segmenting the teeth, as it provided the greatest accuracy in getting good images of the anatomy of the teeth in a consistent manner.

Because Mimics was unable to segment the PDL in due to the low resolution of the CBCT images, it was necessary to grow the PDL from the root structure of the tooth. This required splitting the tooth into root and crown at the CEJ. Once grown, the constructed PDL was essentially two surfaces parallel to each other spaced 0.2 mm apart, where one surface was in intimate contact with the bone and the other with the root. It was critical that the surfaces were tied together in the finite element analysis so that a load added to a tooth was propagated through the PDL to the bone. The engineering software rejected models whose surfaces were too far apart or intersected too much, as this made connecting the surfaces impossible and invalidated the FEA model.

Third, all model surfaces were kept relatively smooth and free of small surface topography that is insignificant for the overall model analysis, such as a projection of extra bone off the buccal cortical surface. Fine elements on projections of anatomy add unnecessary complication to the mesh of the final model by decreasing the size of the elements in complicated areas of fine anatomy, thereby increasing the number of elements in the model. Smaller and more numerous elements increase the computing effort in the final finite element analysis.

The locations of the C_RES when the force was applied in the Y and Z directions were different, represented by the differences in their location along the X direction. However, the difference was small (Table 5) and was clinically as well as statistically insignificant. Therefore, the location of C_RES calculated in one direction may be used for the other. Previous work has also shown that when evaluated in 3D a single point for the C_RES is not observed^10,26,27. Therefore, it has been suggested that instead of having a definite C_RES a better terminology could be "radius of resistance". This difference can be attributed to a number of factors, such as root morphology, boundary conditions, material properties, and point of load application.

Analysis of Force Systems using Custom Algorithms
The mathematical concepts, theoretical derivations, and computer simulations for locating the C_RES of a tooth have been previously described in detail^27,28,29,30. In order to analyze the force systems created by the various loads applied and to predict the C_RES for the teeth, a custom algorithm was written and run within Abaqus (see supplemental coding files). This algorithm was written using Python, accepts data from the FEA software output database (.odb file) as input, processes the data, and provides values for the moments created in the system by the applied load. Additionally, it estimates node locations that result in the generation of a lower moment within the system. This allows the user to run the simulation in an iterative fashion until the estimations converge onto a single location.

The algorithm accesses the nodal coordinates, the total displacement of each node, and the reaction forces at each node as a result of the applied load in each step. Reaction forces in the same direction as the original load application and reaction forces in the opposite direction are summed at each of the nodes in the system to determine the aggregate force vectors acting on the tooth during the simulation. Resulting moments are calculated in relation to the point of force application for each reaction force at each node and are also summed in the same fashion as the reaction forces. Thus, an aggregate force vector in the same direction as the original load application and the resulting moment created by that force vector about the point of force application is calculated, as well as the force vector in the opposite direction and its resulting moment. Because the system is in static equilibrium, the sum of all forces and moments equals zero. However, the breakdown of the reaction forces and moments in this fashion allows for the calculation of the effective locations where these aggregate forces act as pivot points in the system, and the center point between these pivot points provides an approximation of a point of force application that is closer to the C_RES.

In order to perform these calculations, the magnitude of the resulting moments is divided by the magnitude of their respective forces to give the magnitude of the distance (R vector) from the pivot points to the point of force application. The direction of the R vector is determined via a cross product of the moment and force vectors, where all must be orthogonal to each other, and the unit vector is determined by dividing by the magnitude of the cross product. The unit vector R is multiplied by the magnitude of the R vector previously calculated to yield the overall estimation in 3D space of the coordinates of each pivot point relative to the original point of force application. The midpoint between these two vectors provides the estimation for the location of the next point of force application in the following iteration. Additional information is attached in Supplemental Document 2.

The estimation of the C_RES is determined when the resulting moments in the system add to approximately zero. For the current study, this determination is made by finding the lowest positive and negative X-components of the calculated moments and averaging the two. Due to the randomly generated location of the nodes, and the inherent distance between any two nodes (0.5 mm), it is difficult to find a location where a precise zero moment is generated (Table 5).

Limitations
Despite our best efforts, there are some limitations to this study. First, because the PDL could not be visualized on the CBCT, it could not be segmented on its own and was generated from the root surface of the tooth at a uniform thickness of 0.2 mm. Finite element studies have shown that uniform versus nonuniform modeling affects the outcome of the FEA, and that nonuniform modeling is superior^30,31. Second, the number of steps to create an accurate model was lengthy. This is a limitation in terms of how quickly models can be made, which limits the possibility of using these tools for personal treatment plans for patients on a case by case basis. Additionally, the software required to generate these models is expensive and limited to the resources available at an educational institution or a large business. Further, once the models were made, very powerful computing was necessary to run the FEA. Thus, this method cannot be a viable treatment planning tool until the technology necessary is widely available.

Future research should focus on using these models to perform finite element analyses on the maxillary teeth to determine the C_RES for the arch and groups of teeth, especially those groups of teeth typically manipulated in orthodontics, such as the anterior segment in an extraction case or a posterior segment for intrusion in open bite patients. Once the C_RES is determined for these models, additional models should be developed from additional CBCT images to add to the existing data. With a sufficient data pool of C_RES locations, heat maps could be generated to indicate a general position of the C_RES that could serve as an invaluable reference for clinicians.