Surgical Outcome Prediction Using a Four-Dimensional Planning Simulation System With Finite Element Analysis Incorporating Pre-bent Rods in Adolescent Idiopathic Scoliosis: Simulation for Spatiotemporal Anatomical Correction Technique

An optimal surgical strategy for adolescent idiopathic scoliosis (AIS) is to provide maximal deformity correction while preserving spinal mobile segments as much as possible and obtaining a balanced posture. From a spatiotemporal deformity correction standpoint, we recently showed that anatomical four-dimensional (4D) spinal correction could be accomplished by curving the rod. In the surgical procedure, two rods are bent identically to confirm spinal anatomical alignment without referring to the intraoperative alignment of the deformity. Therefore, anatomically designed rods have been developed as notch-free, pre-bent rods for easier anatomical reconstruction. In addition to providing the best spinal instrumentation configurations as pre-bent rods, prediction of surgical outcome along with its biomechanical impact can be obtained by simulation of the surgical procedures with computer modeling. However, an objective model that can simulate the surgical outcome in patients with AIS has not been completely elucidated. The present study aimed to compare simulated deformity corrections based on our newly developed spatiotemporal morphological 4D planning simulation system incorporating pre-bent rods and actual deformity corrections in patients with AIS. A consecutive series of 47 patients who underwent anatomical posterior correction for AIS curves were prospectively evaluated. After multilevel facetectomy, except for the lowest instrumented segment, 11 types of pre-bent rods were used. Patient demographic data, radiographic measurements, and sagittal rod angles were analyzed within 1 week of surgery. Our simulation system incorporating pre-bent rods showed a significant correlation with the actual postoperative spinal alignment. The present study demonstrated the feasibility of our simulation system and the ability to simulate the surgical procedure using the pre-bent rods. The simulation system can be used to minimize the differences between the optimal and possible outcomes related to the instrumentation levels and rod shapes. Preoperative assumption of rod shape and length can contribute to a reduction in operative time which decreases blood loss and risk of infection. The results of the finite element analysis in the simulation system measured for each individual patient would also provide a more realistic representation of the surgical procedures.


INTRODUCTION
Adolescent idiopathic scoliosis (AIS) is the most commonly encountered pediatric musculoskeletal disease presenting a three-dimensional (3D) deformity of the spine. Standard measurement in scoliosis is the Cobb angle, which is the coronal plane angle measured between the vertebrae at the upper and lower bounds of the curve on a standing radiograph. Patients with severe (Cobb angle >40°-50°) or progressive curves may require surgery to correct the deformity.
An optimal surgical strategy is to provide maximal deformity correction while preserving spinal mobile segments to the best extent and obtaining a balanced posture. For instance, inadequate selection for instrumentation length may lead to a postoperative postural imbalance. In addition, although surgical technique as well as spinal instrumentation has been developed in which the 3D correction is achieved, there is still a possibility of implantrelated complications such as pedicle screw loosening, screw or spinal rod breakage, and pedicle fracture. Load levels of the screws and rods are important concerns in surgical outcomes. Furthermore, although a rod shape considerably affects postoperative spinal alignment (Salmingo et al., 2014;Kokabu et al., 2016;Sudo et al., 2016;Le Navéaux, et al., 2017), the rod-bending maneuver relies excessively on surgeons' experience. If the rod curvature does not match the patient's deformity and does not allow for deformity correction, such situations will lead to an inadequate correction or implant-related complications due to the overstress on the implant and spine . These issues require some innovative systems to assist surgery or predict the most probable outcome of surgery (Aubin et al., 2008;Sudo et al., 2021).
The typical thoracic AIS presents itself with thoracic hypokyphosis. Therefore, the surgical goal should be a correction of the thoracic kyphosis (TK) and achieves an anatomically correct thoracic curve. Post-surgery hypokyphosis can occur after using pedicle screw instrumentation. Several posterior surgical techniques have been developed to maintain and/or improve the TK (Clement et al., 2008;Sudo et al., 2014). However, next-generation surgical techniques are required in order to achieve true anatomical correction. In a healthy human population, the apex of the TK is typically located at T6-T8, when viewing standing sagittal films (Hasegawa et al., 2017). However, for some AIS, the postoperative apex of the TK is almost identical with the apex of the preoperative thoracic scoliosis , which is not anatomically correct. This insufficient correction resulting in a postoperative non-anatomical TK is thought to be, because the spinal rods are being bent to match the curvature of scoliosis. From the standpoint of spatiotemporal deformity correction, we recently showed that anatomical fourdimensional (4D) spinal correction could be accomplished by curving the rod (Figure 1) Sudo et al., 2021). In the surgical procedure, two rods are bent in a nearly identical fashion to confirm spinal anatomical alignment without reference to the intraoperative alignment of the deformity Sudo et al., 2021). Consequently, pre-bent rod geometries were obtained from intraoperative tracings of the rod shapes, and optimized rod shapes were derived using iterative closest point method followed by hierarchical cluster analysis . Currently, 11 types of pre-bent cobalt-chrome (CoCr) alloy rods are available based on the deformity types and its lengths in Japan that can guide anatomical spinal correction regardless of the surgeons' experience ( Figure 1, Sudo et al., 2021).
In addition to providing the best spinal instrumentation configurations as pre-bent rods, prediction of surgical outcome along with its biomechanical impact can be obtained by simulation of the surgical procedures with computer modeling (Aubin et al., 2008). However, an objective model that can simulate the 3D outcome of the AIS surgery by considering the preoperative spinal alignment and the surgical intervention has not been completely elucidated (Pasha and Flynn, 2018). In addition, most planning tools in AIS surgery only simulate morphology-based changes of the spinal alignment, lacking the biomechanical analysis (Ferrero et al., 2008;Pasha and Flynn, 2018;Shao et al., 2018). A planning simulator based on spatiotemporal morphological postoperative 4D changes with a patient-specific finite element analysis (FEA) can allow surgeons to predict postoperative outcomes and effectively assist in performing AIS surgery (Galbusera et al., 2015;Wang et al., 2016;Le Navéaux et al., 2016;Cobetto et al., 2020;La Barbera et al., 2021;Galbusera et al., 2021).
The hypothesis of the present study was that our newly developed 4D planning simulation system incorporating prebent rods would significantly correlate with the actual postoperative spinal alignment after anatomical 4D spinal correction surgery. The current study aimed to compare simulated and actual deformity corrections in patients with AIS.

Patient Selection
After institutional review board approval, data from a consecutive series of 47 patients who underwent 4D anatomical correction surgery for AIS curves between 2019 and 2021 were prospectively evaluated; all patients had a Cobb angle of ≤90°. We did not define the lower limit of the Cobb angle. However, patients with severe (Cobb angle >40°-50°) and/or progressive curves were included. The Ethics Committee of Hokkaido University Hospital approved this research including any relevant details. All methods were performed in accordance with the relevant guidelines and regulations. Written consents were obtained from all the subjects, and when applicable from their guardians. The exclusion criteria were neuromuscular, congenital, and other syndromic scoliosis.
Standing posteroanterior radiographs were recorded preoperatively and within 1 week after surgery. Regarding Cobb measurements, the end vertebrae levels were determined on preoperative radiographs and measured on subsequent radiographs to maintain consistency for statistical comparisons (Cidambi et al., 2012;Sudo et al., 2013). The angle of rotation of the main thoracic (MT) and/or thoracolumbar/lumbar (TL/L) apical vertebra was determined on computed tomography (CT) images (Cidambi et al., 2012;Silvestre et al., 2013;Sudo et al., 2014). Internal studies of the present interrater and intrarater reliability have demonstrated high kappa statistics for all continuous measures (0.90-0.98).

Surgical Procedures
While the end vertebrae were to be considered part of the instrumentation levels, the selection of the upper or lower instrumented vertebrae was dependent on several preoperative anatomical conditions. Shoulder balance and anatomical TK determine the vertebra that was selected for the upper instrumented vertebra (UIV); T2 was selected if the radiographic shoulder height (RSH) was positive, T3 if RSH was between −5 and 0 mm, and T4 if RSH was < −5 mm . However, in case with TK < 20°and T5 or T6 upper-end vertebra, the UIV selected was T4 to create anatomical TK . The lowest instrumented vertebra (LIV) depends on the lumbar modifiers. For a lumbar modifier A or B, the last vertebra touching the center sacral vertebral line was the LIV (Matsumoto et al., 2013). In the case of lumbar modifier C, LIV was determined at L3 (Sudo et al., 2021).
Side-loading polyaxial pedicle screws (CVS spinal system; Robert Reid, Tokyo, Japan) were inserted. Our previous studies have shown that multilevel facetectomy and screw density on the concave side rather than the convex side significantly impact scoliosis correction and TK restoration . This means it is important to place as many screws as possible on the concave side. On the convex side, screws should be placed at 1) the UIV, 2) the upper-end vertebra, 3) the lower-end vertebra, 4) the LIV, and 5) at the apex of scoliosis and its periapical lesions. All-level facetectomy was performed in all patients except for the lowest instrumented segment to avoid pseudoarthrosis at this site . For pre-bent rods, CoCr alloy rods (φ 5.5 mm) were bent identically to duplicate the postoperative anatomical TK Sudo et al., 2021). The apex was anticipated to be at T6-T8 for the postoperative TK . The rod configurations were split into two types of shapes: single curve and double curves. In the case that LIV was L1 or above, the single-curve rods were applied Sudo et al., 2021), and the TL/L region remained straight. When the LIV was L2 or L3, the double-curve rods were applied Sudo et al., 2021). Each shape was provided by increments of 3 cm. After connecting to the screw heads, the rods were simultaneously rotated. During the rod derotation maneuver, the present technique helped prevent the hypokyphotic deformation of the rod compared with the simple single-rod derotation maneuver or direct vertebral rotation technique (Sudo et al., 2014). The simultaneous rod rotation maneuver does not intend to manipulate vertebral rotation at each level separately and works to correct the rotational deformity not at each segment separately but in the entire instrumentation area simultaneously (Sudo et al., 2014). After 90°rod rotation, several screw heads were tightened to lock the rods. The presence of a mark on the rod helped confirm 90°rotation. Distraction force was first applied on each screw head on the concave side of the thoracic curve, so that not only scoliosis but also TK could be corrected more effectively by lengthening the posterior column. Subsequently, compression force was applied segmentally on the convex curve. In situ rodbending procedure was not performed (Sudo et al., 2014;Sudo et al., 2016;Sudo et al., 2018).

Biomechanical Model of the Spine
For each patient, a custom spinal finite element model (FEM) was constructed based on the preoperative CT Digital Imaging and Communications in Medicine (DICOM) data ( Figure 1). The software ANSYS 19.2 (ANSYS JAPAN, Tokyo, Japan) was used to model the spine and perform surgical simulations. The collected raw data in the DICOM format were imported into Mimics research 19.0 (Materialize, Leuven, Belgium) to generate 3D vertebral models in a standard triangle language (STL). Subsequently, the STL data generated were imported into ANSYS 19.2 in the form of solid 3D structure (Peng et al., 2020;Zhou et al., 2020), which was built of 10-node tetrahedral element meshes (Ulrich et al., 1998). Screws were positioned and oriented in the desired locations. A new triangulated surface of the instrumented vertebra was generated by Boolean subtraction between the original vertebral surface and the surface of the screws to represent the insertion of screws into the vertebrae. Tetrahedral finite element meshes of the vertebrae and screws were then automatically generated (Galbusera et al., 2015). Eleven types of beam element rods were selected based on the deformity types and its lengths, and positioned for the screws.
Vertebrae were considered as rigid elements in the model to avoid penetration of bone structures. Spinal rods were modeled with a cross-sectional diameter of 5.5 mm, Young's modulus of 420 GPa, and Poisson's ratio of 0.3 according to CoCr alloy properties (Yamada et al., 2020). The screws were modeled as cylinders with a length of 30 mm, cross-sectional diameter of 5.5 mm, Young's modulus of 5,000 MPa, and Poisson's ratio of 0.3 (Zhou et al., 2020). Rods and screws were modeled with materials of isotropic elastic linear material properties. The number of elements in the implants was 790 in the rod and 88 in each screw (Shin et al., 2018).
Connections between the geometries were defined to simulate spatiotemporal morphological postoperative 4D changes. For Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2021 | Volume 9 | Article 746902 stability of simulation, intervertebral discs were completely removed, and each vertebra was connected with joint element having intervertebral stiffness (Galbusera et al., 2015). The stiffness matrix components used in this study (Table 1) were using the ANSYS software. We set the stiffness matrix components based on the literature (Argoubi and Shirazi-Adi, 1996). The intervertebral stiffness was calculated from reaction forces and moments according to the relative amount of translational and rotational displacements between two vertebrae (Senteler et al., 2016). The values were same through T1/2 to L4/5. However, the value of stiffness matrix component was reduced to be representative of facetectomy (Oda et al., 2002) where facetectomy was performed. The connection between vertebra and screw was set as a joint where boundary condition was defined with all translational and rotational degrees constrained. A spring contact model was defined between screw head and rod to simulate the rod being captured in the screws. Finally, a spring contact model was defined between the concave and convex rods to simulate the rotation of rods simultaneously. In this study, distractioncompression force was not applied on each screw head.

Surgical Simulations
During the whole simulation process, boundary conditions were imposed on the spinal FEM to mimic conditions observed in a surgical setting and ensure simulation convergence. Boundary conditions of the spinal FEM were defined with the sacrum fixed and T1 free to rotate and to translate in the caudocranial direction, allowing possible lengthening of the spine during the simulation of the correction process. For multilevel facetectomy simulation, contacts between posterior facets were neglected at instrumented levels except for the lowest instrumented segment to mimic their surgical removal. Regarding rod rotation maneuver, connecting concave and convex rods to screw head was simulated by setting the spring length between the rod and the screw to zero. Furthermore, the power delivered to the concave rod was gradually increased towards concave (66 ± 106 N; range, 0-300 N) and dorsal sides of scoliosis (922 ± 169 N; range, 400-1,000 N) to perform a 90°rotation of the rod. Then, final locking of the screws was simulated. For each screw, null relative translations and rotations between the appropriate rod node and the screw head were imposed. Therefore, the rods could not slide or rotate anymore into the screw heads. Moreover, all external constraints (displacement or forces) were released at this step, and the new equilibrium state was computed. The model was left free to reach equilibrium at the end of the simulation. In this surgical simulation, distractioncompression force was not applied on each screw head. Hence, intraoperative surgical steps described in the simulation included precisely the same as actual surgery except for distractioncompression procedure and screw length and diameter. However, the stiffness matrix does not correctly describe the real condition because this model did not consider preoperative curve flexibility.

Simulation Data Analysis
The spinal profile, quantified in terms of coronal MT Cobb angle, TK and TL/L lordosis, and apical vertebral rotation angle, was monitored over the course of the surgery simulation process. Von Mises stress, which is an equivalent stress, was shown as the reaction forces on the concave and convex rods at the end of the correction (Shin et al., 2018;Peng et al., 2020). Axial forces at the bonescrew interface were analyzed for all the pedicle screws with respect to their simulated final positions and the postoperative ones. The direction of these forces was along the direction of the FIGURE 2 | Rod angle before and after implantation (A) Prior to implantation, the angle between the proximal and distal tangential line was measured (θ1). Postoperative implant rod geometry (θ2) was obtained after the surgery using computed tomography (B) and the simulation model (C).
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2021 | Volume 9 | Article 746902 screw's shaft. They were supposed to be responsible for the pullout phenomenon; therefore, the expression "pull-out forces" was used to describe them. This first phase study did not analyze other forces and moments components such as medio-lateral forces and screw bending moments (responsible for pedicle wall breach and screw breakage, respectively) because the pull-out phenomenon is likely observed compared to the other phenomena in actual AIS surgery (Abul-Kasim and Ohlin, 2014; Oda et al., 2021).

Analysis of Rod Configuration
The angle between the cranial and caudal tangential lines was obtained before implantation (θ1) (Figure 2). Similarly, the postoperative rod angle was obtained (θ2) using reconstructed sagittal CT images and simulation models (Cidambi et al., 2012;Kokabu et al., 2016;Sudo et al., 2021). The angle of rod deformation was defined as the difference between θ1 and θ2 (θ1-θ2) Sudo et al., 2021).

Statistical Analysis
All data were presented as means ± standard deviation and range. The entire cohort was first analyzed and then further assessed based on selective thoracic fusion to L1 in ten Lenke 1 A patients to confirm the feasibility of our simulation system and the ability to simulate the uninstrumented lumbar segments. Repeated-measures analysis of variance was used to compare differences among the standing radiographs, CT images, and simulation data. Data were checked for normality and equality of variances, and Bonferroni post-hoc analysis was used to set the significance level at 0.05. Comparisons of radiographic quantitative variables and rod angles were performed using MannWhitney U test or paired t-test. Spearman's correlation coefficient analysis was used to assess relationships between CT images and simulation models. Data analyses were performed using JMP statistical software for Windows (version 14; SAS, Inc, Cary, NC, United States). p < 0.05 was considered statistically significant.

Patient Demographic Data
Demographic data are summarized in Tables 2, 3. The cephalad-instrumented vertebrae ranged from T2 to T6, and the caudal-instrumented vertebrae ranged from T12 to L3. Preoperative standing radiographic MT and TL/L curves averaged 52°and 37°, respectively, and TK angle was 17°, whereas the lumbar lordosis angle was of 47°. Postoperatively, MT and TL/L curves averaged 11°and 9°, respectively, and TK angle was 30°, whereas the lumbar lordosis angle was 48°. The average MT and TL/L curve correction rate was 81% (range, 64-98%) and 75% (range, 36-98%), respectively. The average preoperative MT and TL/ L vertebral rotation angles measured on CT images were 18°e ach, which decreased after surgery to an average of 11°each.

Comparison Among Standing Radiographs, CT Images, and Simulation Models
Preoperatively, there were significant differences between standing radiographs and CT images or simulation models in both coronal and sagittal plane data (p < 0.05). However, there were no significant differences among postoperative standing radiographs, CT images, and simulation models in both coronal and sagittal plane data (p > 0.05). Regarding postoperative vertebral rotation angle, there was no significant difference between CT images and simulation model (p > 0.05).

Correlation Analysis and Accuracy Evaluation
Spearman's correlation coefficient analysis showed that simulated coronal and sagittal plane data as well as vertebral rotation angle were significantly correlated with those of data on postoperative CT images (p < 0.001, Figure 3). Mean absolute error and root mean squared error between CT images and simulation model are summarized in Table 4. Simulated coronal and sagittal plane data as well as vertebral rotation angle were predicted within 5°compared to actual postoperative measurements.

Subgroup Analysis
There were no significant differences between postoperative CT images and simulation models in all coronal, sagittal, and axial plane data in Lenke 1 A patients ( Table 5, p > 0.05). Regarding mean absolute error and root mean squared error between CT images and simulation model, the simulated coronal and sagittal plane data as well as vertebral rotation angle were predicted within 5°compared to actual postoperative measurements ( Table 6).

Analysis of Rod Stress and Screw Forces
We estimated rod stress and screw force in simulation models. The models showed that peak stress was located near the apex of the curve in both single and double curve type rods (Figure 4). In addition, another peak was located near the extremities of the instrumented segments. Table 7 shows estimated rod stress and screw force in one simulation model. The rod stress was significantly higher at the concave side compared to the convex side (p < 0.05). Regarding pullout forces on screws, peak forces were located near the apex of the MT curve at the concave side in both single and double curves. In addition, another peak was located near the LIV at both concave and convex sides in double curve. The pull-out forces on screws were significantly higher at the concave side compared to the convex side (p < 0.05).

Implant-Rod Angles of Curvature
The rod deformation angle was significantly higher on the concave side than on the convex side (p < 0.001; Table 8).
There were no significant differences between postoperative CT images and simulation model in the rod deformation angle at both concave and convex sides (p 0.129 and p 0.237, respectively).

DISCUSSION
Our personalized finite element spinal biomechanical model and its simulated response to surgical instrumentation allowed to evaluate the effect of the 4D anatomical correction technique on AIS deformity correction and on loads in the instrumentation. Not only geometric aspects of the deformity correction but also biomechanical results were simulated using existing pre-and postoperative information of 47 patients and 11 types of pre-bent rods. Despite the heterogeneity in the cohort of patients with various Lenke types, this study analyzed a cohort of patient-specific surgery models and found that our newly developed 4D planning simulation system incorporating pre-bent rods showed a significant correlation with the actual postoperative spinal alignment after anatomical 4D spinal correction surgery. The simulated measurements were all within 5°agreement with the clinical values, equivalent to the generally accepted clinical error of 5° (Majdouline et al., 2009). The present study demonstrated the feasibility of our simulation system and the ability to simulate the surgical procedure using the pre-bent rods. The preoperative assumption for rod shape and length will help reduce operative time, thereby decreasing blood loss and risk of infection. However, in this model, the values of intervertebral stiffness matrix components were same through thoracic to lumbar spines and not personalized. We consider this as a main factor that may justify imperfect correlations of the slope and intercept. It may be important to consider the personal stiffness matrix for improving the capability of the model in predicting the perfect values. Other factors affecting the 11 ± 6 (1-21) 11 ± 6 (0-23) 0.719 Thoracolumbar/lumbar apical vertebra (°) NA 11 ± 7 (1-30) 10 ± 6 (1-26) 1.000 All data expressed as means ± SD and range.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2021 | Volume 9 | Article 746902 correlations could be that boundary conditions applied during simulation were not the same as during the actual surgery for each patient, and there was an issue that clinical measurements of the parameters of interest were not perfect. In addition, we cannot provide a quantitative justification regarding the prediction accuracy because no steps were used to ensure the credibility of our model in predicting the right values. Although computational models are increasingly used to support surgical planning, varying levels of model verification and validation limit the level of confidence in their predictive potential (Poncelas et al., 2021). Recently, Poncelas et al. performed a credibility assessment of their model to investigate proximal junctional failure in clinical cases with adult spine deformity using ASMEV&V40 standard (Poncelas et al., 2021). We should also assess the credibility of our model for AIS surgery using the recommended strategies in the future.
In the present study, the surgical simulations were conducted using the DICOM CT scans in a supine position and approximated with the surgical procedures performed on patients lying prone. In addition, the postoperative corrected angles were measured clinically using the DICOM CT data in a supine position. Consequently, supine CT scans were obtained after surgery, while the patient was still recovering; therefore, it was not yet load-bearing, and provided a better comparison
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2021 | Volume 9 | Article 746902 between the clinical and simulated measurements (Little et al., 2013). The difference between the present simulation model and reality was small in the instrumented segments. Conversely, there is a possibility that the results for the uninstrumented regions are not accurate because the gravity and postural control were not simulated in the standing posture (Robitaille et al., 2009). Due to the aforementioned reasons, we performed comparisons using standing radiographs. Consequently, while there were significant differences in preoperative measurement values between standing radiographs and simulation model, there were no significant differences in postoperative measurement values between the standing radiographs and simulation model, indicating that the simulation model can predict postoperative spinal alignment in standing position including uninstrumented lumbar segments in the case with selective thoracic instrumentation for Lenke 1 A curves.
There have been computational studies that similarly predicted 3D correction and implant loads Le Navéaux et al., 2016;La Barbera et al., 2021;Galbusera et al., 2021). These studies determined how several instrumentation parameters such as screw density and rod contouring angle affected correction and stress in the instrumentation. Because our previous studies have shown that multilevel facetectomy and screw density on the concave side significantly impact the amount of scoliosis correction and also TK restoration, especially in preoperative hypokyphotic (TK <15°) thoracic spine Sudo et al., 2016), we currently attempt on inserting the screw as much as possible on the concave side. Also, rod curvatures were limited to up to 11 types. We have analyzed the correlation between preoperative rod angle and rod stress in the apex of the MT curve or TL/L curve. In addition, the correlation between postoperative TK in the simulation model and screw density on the concave side or preoperative rod angle in patients with hypothoracic (TK <15°) Lenke one curves was analyzed. The results showed that there were no significant correlations between the simulated correction and rod stress and instrumentation parameters (Supplementary figures S1, S2). Therefore, we could not confirm that clinical observation reported by Sudo et al. (2016) and Kokabu et al. (2016) was confirmed by the present biomechanical models. Because there was no range of numbers for the screw density and rod curvature, as well as sufficient quantity of sample numbers, there may be limitations in statistical analysis.
However, because both pre-and postoperative rod measurements were available in the current study, it was possible to utilize the initial rod shape and simulate its elastic deformation precisely. In this simulation model, the maximum power delivered to the rod was 1,000 N. We previously documented that the notch-free, pre-bent CoCr alloy rod (φ5.5 mm) showed an approximated linear loaddisplacement curve under 1000 N of load (Yamada et al., 2020). Due to these reasons, we opine that almost only elastic deformation occurred. However, elastoplastic phenomena involved in rod contouring may be considered in future studies to better elucidate whether both elastic and plastic deformation may be involved in vivo to explain any change in rod shape in vivo. Although it has been demonstrated that the amount of curvature incorporated into the rods before their insertion impacts TK restoration (Salmingo et al., 2014;Kokabu et al., 2016;Sudo et al., 2016;Le Navéaux, et al., 2017), there are a few studies that have estimated the loads in the rod during and after deformity collection (Galbusera et al., 2015). To our knowledge, this is the first study to show that the simulation model can predict the deformation of the implanted rod. The simulation model determined significant changes in the rod contours, especially on the concave side which has been clinically reported (Sudo et al., 2021). In addition, based on the changes in rod geometry and FEA, the highest stress was found at the apex of the rod curvature and the extremities of the instrumented levels, which is in agreement with previous results (Belmont et al., 2001;Aubin et al., 2008;Abe et al., 2015). Qualitative understanding of the stress in the rods is useful to estimate the risk of implant failure and loosening intraoperatively and/or postoperatively, which currently depends on the surgeon's experience (Galbusera et al., 2015). Further understanding of bonescrew forces in AIS instrumentation is essential as high-stresses at the bonescrew interface can cause screw loosening or breakage. Shear forces on the screws are more relevant to be reported rather than only pullout values. In addition, because we did not observe intraoperative complications such as screw pull-out and bending, we could not analyze the predicted forces in cases with complications. However, in this study, the forces generated at the bonescrew interface (peak 497 N) were lower than the thoracic pedicle screws pull-out forces of approximately 800 N reported in experimental studies (Liljenqvist et al., 2001). The present study validated only the geometrical aspects, and more investigations are needed to validate the model in terms of forces at the bonescrew interfaces. However, the implants tested by Liljenquist et al. were monoaxial pedicle screws which were different from the screws in the current study (poly-axial screws). Bonescrew forces were higher for monoaxial screws than polyaxial screws, indicating that in patients with large and stiff spinal deformities or in patients with compromised bone quality, screws with more degrees of freedom would offer better perspective to reduce bonescrew connection failure . Although the thresholds may serve as a comparison in the present study, therefore, those simulated cases exceeding the threshold may not be considered necessarily unsafe, it is likely that the anatomical correction technique may be used safely.
Our study has some limitations. First, we applied boundary conditions only to the pelvis and did not include the cervical spine, ribs, and scapulae. Although this simplification of the real spine represents a condition wherein the vertebral levels are not entirely fixed, including the cervical vertebrae could demonstrate a more naturalistic behavior of the uninstrumented spinal segment (Majdouline et al., 2012). Second, we only simulated one diameter although the screw diameter is known to have the highest effect on the force to failure compared to screw length (Cho et al., 2010;Bianco et al., 2019). Since this was a first phase study to simulate the surgical procedure, the simulation was maintained as simple as possible. Additionally, model validation was purely based on the final rod geometry and the main spinal curves. However, to ensure that the simulation model can predict postoperative alignment, a more detailed validation on single vertebra position and rotation, together with individualized screw's models and trajectories, would be needed. Third, this study did not analyze other forces and moments such as medio-lateral forces and screw bending moments. The other loading components may play a role in other clinically relevant failure modes and may be addressed in the future. Finally, in this surgical simulation, preoperative curve flexibility was not considered, and distractioncompression force was not applied on each screw head. Nonetheless, the simulation model can predict postoperative surgical alignments. However, All data expressed as means ± SD and range. Positive value in the screw force means pull-out force. UIV, upper instrumented vertebra; MT, main thoracic; LIV, lowest instrumented vertebra, TL/L, thoracolumbar/lumbar. <0.001 Simulation model 6.5 ± 3.2 (0.2-12.6) 2.5 ± 2.9 (−3.9-7.6) <0.001 All data expressed as means ± SD and range. Δθ was defined as the difference between θ1 and θ2 (θ1-θ2).
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2021 | Volume 9 | Article 746902 we are currently improving the model to incorporate both preoperative curve flexibility based on bending radiographs and distractioncompression procedures, as well as actual screw length and diameter. There may be a limit to improve the comparison between model prediction and actual postoperative correction; however, the incorporation may further improve the estimation of rod stress and screw force because this information will contribute to set stiffness matrix components and connections between the geometries.

CONCLUSION
Our newly developed 4D planning simulation system incorporating pre-bent rods showed a significant correlation with the actual postoperative spinal alignment after anatomical 4D spinal correction surgery. The present study demonstrated the feasibility of our simulation system and the ability to simulate the surgical procedure using pre-bent rods. The FEA results in the simulation system measured for each individual patient would also provide a more realistic representation of the surgical procedures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Hokkaido University Hospital. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.