Optimization of Loading Condition for Maxillary Molar Intrusion with Midpalatal Miniscrews by Using Finite Element Analysis

An anterior open bite is one of the most difficult malocclusions in orthodontic treatment. For such malocclusion, orthodontic miniscrew insertion into both buccal and palatal alveolar regions has been indicated for molar intrusion, but it involves a risk of tooth root injury. To solve the problem, a midpalatal miniscrew-attached extension arm (MMEA) is adopted. However, this method causes palatal tipping of the molar because intrusive loads were applied only from the palatal side. Currently, a transpalatal arch is added to avoi0d tipping movement, but it induces the patient’s discomfort. Hence, the objective of this study was to evaluate the loading conditions for maxillary molar intrusion without tipping movement, only by MMEA through finite element (FE) analysis. FE models of maxillary right first molar and surrounding tissues were created. Three hook positions of MMEA were set at 6.0 mm perpendicular intervals in the occluso-apical direction along the mucosal contour. An intrusive unit load was applied from the palatal side of the molar, and various counter loads were applied from the buccal side. An optimal counter load for molar intrusion without palatal tipping was observed in each hook position. In conclusion, an ideal maxillary molar intrusion can be achieved only by MMEA with an optimal counter load.


Introduction
Orthodontic treatment is performed to obtain optimal occlusion in the functional harmony as well as well-balanced dental and facial esthetics in each individual [1]. When an orthodontic force is applied, alveolar bone resorption occurs at the compression side of tooth, whereas bone formation is induced in the tension side. Thus, tooth movement is observed in specific directions [2]. During orthodontic treatment, undesirable side effects are also observed, such as pain caused by chewing, loss of tooth vitality, and root resorption [3]. These side effects are induced by molecular mechanisms. For example, the first and second ones are reported to be associated with early phase reduction in alkaline phosphatase activity in dental pulp tissue [4].
In orthodontic treatment, malocclusion with an anterior open bite is considered one of the most difficult problems to treat because it occurs due to interaction between numerous etiologic factors such as genetic, dental, skeletal, functional, soft tissue, and habit [5]. Generally, clinical features such as increased lower facial height, increased gonial and mandibular plane angles, and increased molar dentoalveolar height have been associated with anterior open bites [6]. Among various treatment approaches, molar intrusion is recommended for anterior open-bite patients [7]. Recently, treatment with temporary anchorage devices (TADs), including miniplates and miniscrews, has been established [8], and an application of TADs to molar intrusion has improved the prognosis of treatment for anterior open bites by reducing the uncertainty of the outcome owing to patient cooperation [9][10][11][12]. In addition, molar intrusion with TADs reduces a need for orthognathic surgery in borderline open-bite cases [13].
Some studies recommended the application of miniscrews in both palatal and buccal alveolar regions of the molars, because the intrusive load can be alternated independently in each region and molar movement was controlled precisely [14,15]. However, this method causes root proximity to the miniscrews as the intrusion progresses, and it induces failure of miniscrews [16]. Moreover, this method inevitably involves serious risks of sinus perforation and injury to tooth roots, nerves, blood vessels, etc., since the devices require drilling into the bone tissue for fixation [17,18]; therefore, placement sites of miniscrews are restricted.
As an effective solution to these problems, the midpalatal area is adopted as a placement site of miniscrews to perform maxillary molar intrusion. The use of midpalatal miniscrews is considered to be a useful and less invasive procedure that can ensure good stability [19]. One approach with a midpalatal miniscrew is to apply an intrusive load from a midpalatal miniscrew-attached extension arm (MMEA) ( Figure 1) [19,20]. Currently, in molar intrusion with MMEA, a transpalatal arch (TPA) is applied together because using only MMEA will cause a palatal tipping movement of the molar. However, TPA induces patient's discomfort and jeopardizes oral hygiene, thus alternative methods to TPA have to be considered. Hence, the objective of this study was to optimize the loading condition for maxillary molar intrusion avoiding tipping movement by using only MMEA through finite element (FE) analysis.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 2 of 13 In orthodontic treatment, malocclusion with an anterior open bite is considered one of the most difficult problems to treat because it occurs due to interaction between numerous etiologic factors such as genetic, dental, skeletal, functional, soft tissue, and habit [5]. Generally, clinical features such as increased lower facial height, increased gonial and mandibular plane angles, and increased molar dentoalveolar height have been associated with anterior open bites [6]. Among various treatment approaches, molar intrusion is recommended for anterior open-bite patients [7]. Recently, treatment with temporary anchorage devices (TADs), including miniplates and miniscrews, has been established [8], and an application of TADs to molar intrusion has improved the prognosis of treatment for anterior open bites by reducing the uncertainty of the outcome owing to patient cooperation [9][10][11][12]. In addition, molar intrusion with TADs reduces a need for orthognathic surgery in borderline open-bite cases [13].
Some studies recommended the application of miniscrews in both palatal and buccal alveolar regions of the molars, because the intrusive load can be alternated independently in each region and molar movement was controlled precisely [14,15]. However, this method causes root proximity to the miniscrews as the intrusion progresses, and it induces failure of miniscrews [16]. Moreover, this method inevitably involves serious risks of sinus perforation and injury to tooth roots, nerves, blood vessels, etc., since the devices require drilling into the bone tissue for fixation [17,18]; therefore, placement sites of miniscrews are restricted.
As an effective solution to these problems, the midpalatal area is adopted as a placement site of miniscrews to perform maxillary molar intrusion. The use of midpalatal miniscrews is considered to be a useful and less invasive procedure that can ensure good stability [19]. One approach with a midpalatal miniscrew is to apply an intrusive load from a midpalatal miniscrew-attached extension arm (MMEA) ( Figure 1) [19,20]. Currently, in molar intrusion with MMEA, a transpalatal arch (TPA) is applied together because using only MMEA will cause a palatal tipping movement of the molar. However, TPA induces patient's discomfort and jeopardizes oral hygiene, thus alternative methods to TPA have to be considered. Hence, the objective of this study was to optimize the loading condition for maxillary molar intrusion avoiding tipping movement by using only MMEA through finite element (FE) analysis. The geometry of the maxillary first molar was created based on published data [21]. The size of the crown of the molar was fixed as 10.0 mm mesio-distally, 12.0 mm buccopalatally, and 7.0 mm occluso-apically. The buccal and palatal roots were represented as conical structures, 12.0 mm in height. The maximum radius of the palatal root was designated as 4.0 mm and that of the buccal roots was 3.0 mm. These roots partially overlapped in the trifurcation area ( Figure 2A). A layer structure was formed in the order of PDL, cortical bone, and trabecular bone from the root surface. The PDL models had a thickness of 0.3 mm [22]. The cortical and trabecular bone models had a rectangular parallelepiped shape, with dimensions of 18.0 mm mesio-distally, 16.0 mm bucco-palatally, and 15.0 mm occluso-apically. The surface layer of 1.0 mm in thickness was defined as the cortical bone [23], and the others were defined as the trabecular bone ( Figure 2B). The palatal mucosal contour was drawn on the plane that is perpendicular to palatal cortical bone surface and through the mesio-distal midline of the surface. Three landmarks were plotted at 6.0 mm intervals in the occluso-apical direction according to a mucosal thickness that measured with computed tomography (CT) images [24], and these values were 2.3 mm, 3.0 mm, and 4.7 mm in the order from cervical region. These landmarks were connected with curves and the palatal mucosal contour was created ( Figure 2C). The three-dimensional XYZ coordinates were constructed as follows: the center of the occlusal surface of the molar was defined as the origin, the mesio-distal direction as the x-axis, the bucco-palatal direction as the y-axis, and the occluso-apical direction as the z-axis ( Figure 2).
No ethics approval was required because no animal experiments or human studies were involved in this research. To simulate maxillary molar intrusion with MMEA, FE models were created using an integrated FE analysis software (Femap with NX Nastran: Siemens PLM, USA and CANADA) on a personal computer (ThinkCentre M92p; Lenovo, Hong Kong). The models consisted of the maxillary right first molar and its surrounding tissues, including periodontal ligament (PDL), cortical bone, trabecular bone, and palatal mucosa. These models were created as three-dimensional structures, but only the palatal mucosa was created as a two-dimensional contour.
The geometry of the maxillary first molar was created based on published data [21]. The size of the crown of the molar was fixed as 10.0 mm mesio-distally, 12.0 mm buccopalatally, and 7.0 mm occluso-apically. The buccal and palatal roots were represented as conical structures, 12.0 mm in height. The maximum radius of the palatal root was designated as 4.0 mm and that of the buccal roots was 3.0 mm. These roots partially overlapped in the trifurcation area ( Figure 2A). A layer structure was formed in the order of PDL, cortical bone, and trabecular bone from the root surface. The PDL models had a thickness of 0.3 mm [22]. The cortical and trabecular bone models had a rectangular parallelepiped shape, with dimensions of 18.0 mm mesio-distally, 16.0 mm bucco-palatally, and 15.0 mm occluso-apically. The surface layer of 1.0 mm in thickness was defined as the cortical bone [23], and the others were defined as the trabecular bone ( Figure 2B). The palatal mucosal contour was drawn on the plane that is perpendicular to palatal cortical bone surface and through the mesio-distal midline of the surface. Three landmarks were plotted at 6.0 mm intervals in the occluso-apical direction according to a mucosal thickness that measured with computed tomography (CT) images [24], and these values were 2.3 mm, 3.0 mm, and 4.7 mm in the order from cervical region. These landmarks were connected with curves and the palatal mucosal contour was created ( Figure 2C). The three-dimensional XYZ coordinates were constructed as follows: the center of the occlusal surface of the molar was defined as the origin, the mesio-distal direction as the x-axis, the bucco-palatal direction as the y-axis, and the occluso-apical direction as the z-axis ( Figure 2).

Finite Element Models
The models consisted of tetrahedral solid elements with a side of 1.0 mm. The numbers of nodes and elements for each model are listed in Table 1. The mechanical properties utilized in the models, including Young's modulus and Poisson's ratio (Table 2), were applied according to previous studies [25]. The interface between each material was fully

Finite Element Models
The models consisted of tetrahedral solid elements with a side of 1.0 mm. The numbers of nodes and elements for each model are listed in Table 1. The mechanical properties utilized in the models, including Young's modulus and Poisson's ratio (Table 2), were applied according to previous studies [25]. The interface between each material was fully bonded. The elements on the lateral and bottom surfaces of the cortical and trabecular bones were fully constrained.  When performing molar intrusion with MMEA, a position of the hook is considered to affect the behavior of the molar. To evaluate such influence, we set three vertical positions of the hook (coronal; Co, middle; Md, and apical; Ap) at 6.0 mm perpendicular intervals in the occluso-apical direction along the palatal mucosal contour ( Figure 2C). Regarding an intrusive load, a point of application for the load was 3.5 mm above the cervix along the medio-distal midline on the palatal surface of the crown. A unit load of 1.0 N was applied from the point to each hook position. In addition, a load to the buccal direction was also applied to simulate expansion by the buccal arch wire. The load was defined as a "counter load" to avoid palatal tipping of the molar. A point of application for the counter load was on the buccal surface of the molar, and it was symmetrical to the application point of the intrusive load with the XZ plane. The counter load was applied perpendicularly to the buccal surface, and its magnitude was varied from 0 N to 1.5 N in 0.1 N steps ( Figure 3A).

Materials
Number of the  When performing molar intrusion with MMEA, a position of the hook is considered to affect the behavior of the molar. To evaluate such influence, we set three vertical positions of the hook (coronal; Co, middle; Md, and apical; Ap) at 6.0 mm perpendicular intervals in the occluso-apical direction along the palatal mucosal contour ( Figure 2C). Regarding an intrusive load, a point of application for the load was 3.5 mm above the cervix along the medio-distal midline on the palatal surface of the crown. A unit load of 1.0 N was applied from the point to each hook position. In addition, a load to the buccal direction was also applied to simulate expansion by the buccal arch wire. The load was defined as a "counter load" to avoid palatal tipping of the molar. A point of application for the counter load was on the buccal surface of the molar, and it was symmetrical to the application point of the intrusive load with the XZ plane. The counter load was applied perpendicularly to the buccal surface, and its magnitude was varied from 0 N to 1.5 N in 0.1 N steps ( Figure 3A). Co. An intrusive unit load of 1.0 N was applied from the palatal side of the crown to the hook position. In addition, a counter load to the buccal direction was applied to simulate expansion by the buccal arch wire, and its magnitude was varied from 0 N to 1.5 N in 0.1 N steps. (B) Two measurement points were defined to evaluate the horizontal and vertical displacements of the molar. "Crown" was the center of the occlusal surface, and "Root" was the center of gravity among buccal and palatal root apical. The displacement of the molar was measured on each y-and z-axis according to the change in these points. Positive values in the y-and z-axis directions indicate displacement to the palatal side and intrusion, respectively. Co. An intrusive unit load of 1.0 N was applied from the palatal side of the crown to the hook position. In addition, a counter load to the buccal direction was applied to simulate expansion by the buccal arch wire, and its magnitude was varied from 0 N to 1.5 N in 0.1 N steps. (B) Two measurement points were defined to evaluate the horizontal and vertical displacements of the molar. "Crown" was the center of the occlusal surface, and "Root" was the center of gravity among buccal and palatal root apical. The displacement of the molar was measured on each y-and z-axis according to the change in these points. Positive values in the yand z-axis directions indicate displacement to the palatal side and intrusion, respectively.

Evaluation
Two measurement points were defined to evaluate the horizontal and vertical displacements of the molar. One was the "crown", which was the center of the occlusal surface; the other was the "root", which was the center of gravity among buccal and palatal root apical. The displacement of the molar was measured on each y-and z-axis according to the change in these points. Positive values in the y-axis and z-axis directions indicate displacement to the palatal side and intrusion, respectively ( Figure 3B).
Based on the result of the horizontal displacement of the crown, we calculated the optimal counter load without tipping movement of the crown during intrusion. The procedure of calculation is as follows: (i) A regression line was derived from the data of horizontal displacement of the crown for each hook position. The line can be written as: where y was the value of horizontal displacement, x was the value of the counter load, and a and b were coefficients of approximation.
(i) x value was calculated by substituting zero for y.
(ii) The x value was defined as the optimal counter load.
In addition, the von Mises stress distributions on the PDL around the palatal root apical (PDL-P) and that around the buccal root apical (PDL-B) were estimated to evaluate the biomechanical response.

Statistical Analysis
The data from horizontal and vertical displacements and stress distribution were analyzed using Pearson's correlation coefficient to evaluate the correlation between these data and counter load. Correlation coefficients of 0 < |r| ≤ 0.2 indicated no correlation, 0.2 < |r| ≤ 0.4 indicated a weak correlation, 0.4 < |r| ≤ 0.7 indicated a moderate correlation, and 0.7 < |r| < 1.0 indicated a high correlation. Further, for all results in horizontal and vertical displacements, linear regression was performed, and the regression coefficients were compared between the groups using an analysis of covariance (ANCOVA). Values of p < 0.05 were considered statistically significant. All statistical analyses were carried out using "R" software (version 4.0.2, http://www.r-project.org/ (accessed on 7 June 2021)).     Calculated regression lines and the results of ANCOVA were shown in Table 3. In horizontal dimension, the hook positions did not significantly affect the displacements of the crown and root (p > 0.05), whereas in the vertical dimension, Co showed significantly smaller displacements in both crown and root than the other hook positions (p < 0.05). In all analysis, displacements of x-axis direction were not observed.

Horizontal and Vertical Displacement
Hook positions indicate as follows: Co; coronal, Mi; middle, Ap; apical. The asterisk represents significant deference among groups (p < 0.05).

Optimal Counter Load
Based on the regression lines of horizontal displacement of the crown, the optimal counter loads were calculated as follows: 1.30 N for Co, 1.22 N for Mi, and 1.14 N for Ap. To demonstrate the comparison between no counter load and optimal counter load conditions, representative models of Mi are shown in Figure 5. Notable palatal tipping of the crown and excessive stress concentrations were observed on the root surface, the cortical bone around the palatal root, and the PDL around the root apical regions after intrusive load application without a counter load. On the other hand, application of the optimal counter load showed an intrusion in the z-axis direction with no palatal tipping of the crown ( Figure 5A). Furthermore, the stress concentrations mentioned above were dispersed. (Figure 5B-D).  With no counter load, notable palatal tipping of the crown (A(i)) and excessive stress concentrations were observed on the root surface (B(i)), the cortical bone around the palatal root (C(i)), and the PDL around the root apical regions (D(i)). With optimal counter load, intrusion avoiding tipping movement was observed (A(ii)), and stress concentrations were dispersed (B(ii),C(ii),D(ii)).

Von Mises Stress Distribution on PDL
Since the root geometry of our FE model was linearly symmetrical with the YZ plane and all loads that applied in this study were along to YZ plane, a similar stress distribution was observed in two buccal roots. In all hook positions, the calculated Pearson's correlation coefficient indicated a high negative relationship between the counter load and the stress value on PDL-P

Von Mises Stress Distribution on PDL
Since the root geometry of our FE model was linearly symmetrical with the YZ plane and all loads that applied in this study were along to YZ plane, a similar stress distribution was observed in two buccal roots. In all hook positions, the calculated Pearson's correlation coefficient indicated a high negative relationship between the counter load and the stress value on PDL-P and PDL-B (r = −0.9190 and p-value = 0 for Co, −0.7817 and p-value = 0.003 for Mi, and −0.7551 and p-value = 0.0049 in PDL-B; r = −0.9986 and p-value = 0 for Co, −0.9961 and p-value = 0 for Mi, and −0.9962 and p-value = 0 in for Ap in PDL-P), but the stress on PDL-P decreased drawing a curve, and that of PDL-B showed a parabolic-like change regardless of the hook positions ( Figure 6). In both PDL-P and PDL-B, Co showed lower stress value than the other hook position.

Von Mises Stress Distribution on PDL
Since the root geometry of our FE model was linearly symmetrical with the YZ plane and all loads that applied in this study were along to YZ plane, a similar stress distribution was observed in two buccal roots. In all hook positions, the calculated Pearson's correlation coefficient indicated a high negative relationship between the counter load and the stress value on PDL-P and PDL-B (r = −0.9190 and p-value = 0 for Co, −0.7817 and p-value = 0.003 for Mi, and −0.7551 and p-value = 0.0049 in PDL-B; r = −0.9986 and p-value = 0 for Co, −0.9961 and p-value = 0 for Mi, and −0.9962 and p-value = 0 in for Ap in PDL-P), but the stress on PDL-P decreased drawing a curve, and that of PDL-B showed a parabolic-like change regardless of the hook positions ( Figure 6). In both PDL-P and PDL-B, Co showed lower stress value than the other hook position.
In all analyses, no error messages were reported in output files.

Discussion
In orthodontic treatment, decision of loading condition may primarily depend on the clinical experience or subjective judgment of orthodontists [22]. Biomechanical evidence for various loading conditions would provide indispensable information to establish a treatment plan; however, it is not justified from an ethical aspect to repeat clinical trials for examining numerous loading conditions. Among various attempts to overcome this difficulty, FE analysis has been commonly utilized. FE analysis is a numerical method that In all analyses, no error messages were reported in output files.

Discussion
In orthodontic treatment, decision of loading condition may primarily depend on the clinical experience or subjective judgment of orthodontists [22]. Biomechanical evidence for various loading conditions would provide indispensable information to establish a treatment plan; however, it is not justified from an ethical aspect to repeat clinical trials for examining numerous loading conditions. Among various attempts to overcome this difficulty, FE analysis has been commonly utilized. FE analysis is a numerical method that calculates the stresses and determines the mechanical behavior of complex structures [26]. Recently, the predictability of FE analysis in orthodontic treatment has been validated by evaluation with clinical feedbacks [27]. Thus, we evaluated the loading condition for maxillary molar intrusion with MMEA through FE analysis.
In malocclusion with an anterior open bite, position of the first molar is strongly associated with vertical occlusal relationship [28]. Therefore, we focused on the maxillary first molar as a target to evaluate the loading condition during intrusion. Regarding the intrusive load, the magnitude of the load varies depending on clinical procedure of the orthodontist and the treatment objective. Even in previous FE studies performing molar intrusion, a large variety of the load was adopted [25,29,30]. Hence, the biomechanical evidence with the smallest unit of load was necessary to compare with other loads, and we set 1.0 N as an intrusive load.
Concerning the modeling procedure, we utilized simplified average model. Creating the simplified average model does not need much labor and time that are required in modeling based on CT images [31]. In addition, the simplified average model can generalize and evaluate the loading condition without considering numerous varieties of the teeth owing to tooth wear [28] or occlusal trauma [32] that were caused by anterior open bite. However, an estimation error due to simplified modeling should be considered. To evaluate the estimation error, the analysis results of molar intrusion without tipping movement were compared with a previous FE study using CT images. Cifter et al. constructed FE models of posterior dental segments based on CT images and an intrusive load of 3.0 N from buccal side was applied [25]. Regarding the ratio of the vertical displacement to the stress on the root surface of the first molar, the previous study demonstrated the ratio in a range from 1.3 to 2.1 µm/MPa. On the other hand, our analysis with a simplified FE model demonstrated the ratio in a range from 1.5 to 2.0 µm/MPa. Although these analysis Appl. Sci. 2021, 11, 11749 8 of 13 methods were different in terms of the number of teeth, the magnitude of the load, and modeling procedure, these ratios were almost in same ranges. It indicated that our analysis method using simplified model could produce the result almost in the same range as the FE method using CT images.
In displacement results of the crown and root, horizontal displacement had no significant difference in hook positions, on the other hand, vertical displacement showed significant difference between Co and the others (Table 3). Considering with the result of the line graph (Figure 4), in coronal hook position, the amount of intrusion was significantly smaller than the other hook positions. The reason for these results is explained as follows: owing to the sloped shape of the palatal mucosa, the increment from Co to Ap was 2.4 mm in the horizontal direction, whereas 12.0 mm was shown in vertical direction ( Figure 2C). Therefore, the vector of the intrusive load had much larger difference in vertical component among three hooks than that of horizontal component. In addition, the vector direction in Mi and Ap were similar compared with Co.
As shown in the result of displacements, application of the counter load avoids palatal tipping movement during intrusion with MMEA. However, since excessive counter load causes buccal tipping contrarily, the optimum counter load must be calculated. In the results of the optimal counter load, the required load was the largest in Co, followed by Mi and Ap because the vector of the intrusive load had the largest horizontal component in Co (Table 3). It indicates that heavy counter load should be applied to avoid tipping movements if the hook is placed near the cervical region owing to anatomical problems, such as palatal shape.
To obtain the optimal counter load, the selection of arch wire size and required amount of expansion should be investigated. To reveal them, we performed another analysis using cantilever wire models. For the cross-sectional shape of the wire, three types of rectangular arch wires utilized in clinical practice were prepared. The height and width of the cross sections were 0.016 × 0.022 (0.406 × 0.559 mm), 0.019 × 0.025 (0.483 × 0.635 mm) and 0.021 × 0.025 (0.533 × 0.635 mm), respectively. The length of the wire was set to 6.0 mm; that is, an inter-bracket span between second premolar bracket and first molar buccal tube in maxillary arch [33]. The material was stainless steel and its mechanical properties including Young's modulus and Poisson's ratio were 200.0 GPa and 0.3 [25]. The material was assumed to be homogenous, isotropic, and linearly elastic. The model consisted of tetrahedral solid elements with a side of 0.2 mm. The elements on one end of the wire were fully constraint, and optimal counter load of each hook was applied to the vertical edge in the other end, perpendicularly to the long axis of the model (Appendix A, Figure A1). The displacement in the loading part was measured and defined as the amount of required expansion to obtain each counter load. The results were shown in Appendix A, Table A1. Since the optimal counter load is the value for intrusive load of a unit load, the result increases linearly as the intrusive load increases. If the intrusive load reaches 4.0 N, which is the largest value in previous studies [34], the amount of required expansion is calculated to be only 0.276-0.312 mm even in the 0.016 × 0.022 wire that is the narrowest of the three. It indicates that the optimal counter load can be obtained with a small amount of expansion regardless of the wire size. However, this expansion load is not fully demonstrated in clinical practice when there is a large bracket/wire play, and it causes improper teeth movement [35]. Therefore, in molar intrusion with MMEA, arch wire that minimizes the bracket/wire play should be selected to achieve molar intrusion without tipping movement. In the case that 0.021 × 0.025 stainless steel wire is set on 0.022 × 0.028 slots, the wire should be expanded at first molar region in a ratio of 0.036-0.041 mm to intrusive load of 1.0 N.
In order to evaluate the effectiveness of intrusion with MMEA, we carried out additional analysis. Intrusion with miniscrews in the palatal and buccal region was analyzed using the same model. In addition to the intrusive load in Mi, symmetrical load with the XZ plane was applied in the buccal side. Appendix B, Figure A2 shows the result of the analysis. Comparing with the result of MMEA with optimal counter load in Mi, molar intrusion without tipping movement was also observed, but vertical displacement and stress on PDL were approximately one and a half times larger than those of MMEA. However, in clinical practice, bucco-palatal miniscrews cannot necessarily provide the intrusive loads symmetrically owing to the alveolar bone morphology and the location of the tooth root. In other words, molar intrusion with MMEA can be expected similar tooth movement to bucco-palatal miniscrews, and the risk of root resorption can be reduced although the amount of vertical displacement is smaller. Regarding the load direction, MMEA can determine the load direction freely by adjusting the shape of the arm. In addition, to perform maxillary molar intrusion in both sides, MMEA requires only two miniscrews, whereas intrusion with bucco-palatal miniscrews requires four miniscrews. It contributes to reduction in invasive and economic burdens on the patient. From these perspectives, we consider that MMEA would have some advantages for molar intrusion.
In orthodontic treatment, root resorption is a frequent consequence of tooth movement, especially intrusion and tipping movement [36,37]. As shown in Figure 5, the optimal counter load brought no tipping movement, and the stress value on the PDL became lower. Although it remains a matter of speculation because stress distribution would not be accurate owing to the simplified root shape, the risk of root resorption may be reduced under optimal load conditions [38]. Regarding the maximum von Mises stress values on the PDL, although ANCOVA cannot be performed since the results draw curves, the results appeared to be divided into two groups, i.e., Co and the others ( Figure 6). This characteristic was similar to the result of vertical displacement. Gupta et al. [39] described that intrusion caused larger stress on PDL than tipping movement when the load magnitude was same. Thus, the stress on the PDL might be affected by vertical displacement. A possible explanation of a parabolic-like change on PDL-B is that the palatal tipping movement disappeared, and buccal tipping movement occurred instead as the counter load increased. Since the stress changed along with the magnitude of tipping movement [39], it showed a parabolic-like change.
In our study, there were at least four limitations. The first is related to the mechanical property. A linear elastic analysis was performed for simplicity; however, the bone tissues and PDL were anisotropic and inhomogeneous. In fact, the mechanical properties of the PDL have not been fully clarified. PDL possesses anisotropic and non-linear characteristics owing to the orientation of collagen fibers [40]. However, tooth movement can be approximated to clinical results, even if the PDL was assumed to be a linear, isotropic material [41]. Furthermore, the difference between the linear and non-linear values can be ignored in the FE analysis for tooth intrusion within a load of 2.5 N [42]. Therefore, the PDL was assumed to be a linear, isotropic, and homogeneous material in the study. Secondly, due to the simplified geometries of the tooth, stress distribution of PDL could not be evaluated accurately. Stress distribution determines biological response of PDL that initiates tooth movement or PDL necrosis [43]. Hence, in this study, the biological response to stress was also uncertain. Thirdly, our FE model was simplified compared with clinical condition. Because we applied a single tooth model, transmission of intrusive load to adjacent teeth is not calculated. Therefore, the tipping movement and the required counter load in our model are larger than in a non-simplified model. However, if this model is interpreted as the condition with maximum consideration for possibility of tipping movement, our result can potentially provide the safest clinical index. The final imitation was that the counter load was considered in only one direction, but the load could be applied in various directions in clinical practice. Various directions of the counter load should be considered in further studies.
As an outlook for the FE study in orthodontics, we consider that it is necessary to clarify the relationship between the calculated stress and the biological response such as alveolar bone remodeling and tooth movement in detail. By introducing these mechanobiological evidence, the results of FE analysis would have biological support and become clinically indispensable information. Furthermore, they would support diagnosis and establish a treatment plan.

Conclusions
It was suggested that an ideal maxillary molar intrusion avoiding palatal tipping can be achieved by MMEA with an optimal counter load without an additional TPA. Co exhibited a significantly lesser molar intrusion in need of a greater counter load than Mi and Ap did, whereas similar palatal tipping was observed among the three hook positions. Size matching of the archwire and the bracket slot was more essential than archwire activation for expansion to obtain an optimal counter load. Nevertheless, further FE studies employing non-linear analysis with accurate geometry of the models and various counter load directions should be considered.  Conflicts of Interest: The authors declare no conflict of interest.
Appendix A Figure A1. In order to investigate the selection of arch wire size and required amount of expansion to obtain the optimal counter load, FE analysis using cantilever wire models was performed. (A) Cross sections and three-dimensional geometry of the model. (B) The result of analysis in ten times display (0.016" × 0.022" with optimal counter load of 1.22 N in Mi). The color bar indicates the amount of displacement. Table A1. Required expansion of arch wire to obtain optimal counter load for each hook position.  Appendix B Figure A1. In order to investigate the selection of arch wire size and required amount of expansion to obtain the optimal counter load, FE analysis using cantilever wire models was performed. (A) Cross sections and three-dimensional geometry of the model. (B) The result of analysis in ten times display (0.016" × 0.022" with optimal counter load of 1.22 N in Mi). The color bar indicates the amount of displacement. Table A1. Required expansion of arch wire to obtain optimal counter load for each hook position.