Data-driven Classification of the 3D Spinal Curve in Adolescent Idiopathic Scoliosis with an Applications in Surgical Outcome Prediction

Adolescent idiopathic scoliosis (AIS) is a three-dimensional (3D) deformity of the spinal column. For progressive deformities in AIS, the spinal fusion surgery aims to correct and stabilize the deformity; however, common surgical planning approaches based on the 2D X-rays and subjective surgical decision-making have been challenged by poor clinical outcomes. As the suboptimal surgical outcomes can significantly impact the cost, risk of revision surgery, and long-term rehabilitation of adolescent patients, objective patient-specific models that predict the outcome of different treatment scenarios are in high demand. 3D classification of the spinal curvature and identifying the key surgical parameters influencing the outcomes are required for such models. Here, we show that K-means clustering of the isotropically scaled 3D spinal curves provides an effective, data-driven method for classification of patients. We further propose, and evaluate in 67 right thoracic AIS patients, that by knowing the patients’ pre-operative and early post-operation clusters and the vertebral levels which were instrumented during the surgery, the two-year outcome cluster can be determined. This framework, once applied to a larger heterogeneous patient dataset, can further isolate the key surgeon-modifiable parameters and eventually lead to a patient-specific predictive model based on a limited number of factors determinable prior to surgery.

the surgical intervention has not been fully developed 9 . The successful prediction of the outcome can improve the patient satisfaction and reduces the costs and risks associated with a need for revision surgery 24 .
As the 3D curvature of the spine can explain a large amount of variability in the spinal alignment among the patients, a 3D classification of the overall spinal curve, as opposed to the segmental uni-planar measurements of the curve, can specify the true curvature of the spinal deformity in AIS patients. Similarly, the 3D surgical outcome evaluation of spine can assure that the harmony of the spinal configuration in the three anatomical planes is not overlooked 9 . Patient classification based on the 3D spinal curve, both pre-and post-operatively, can reduce the errors associated with the subjective 2D classification by an objective use of the 3D overall shape of the spine.
The current study aims to develop an objective, data-driven framework for classification-based outcome prediction of the spinal surgery in AIS. More specifically, the study herein focuses on determining pathways, comprised of the variables that can be identified prior to the surgery or modified during the surgery, and statistically examine the association of these pathways to the 3D shape of the spine at two-year follow-ups. The underlying assumption of this approach is that the biomechanical changes induced during the surgery can be described using a combination of the 3D shape of the spine at early post-operative and the location of the vertebral levels that are immobilized via a metal rod and are fused within two-year post-operative thus changing the mechanical behaviour of the spine in the instrumented section.
To identify these pathways and their relationships with the outcome of surgery, we hypothesize that the combination of the 3D curvature of the spine before surgery and the biomechanical impact of the surgery on the spinal alignment (the early post-operation 3D shape of the spine plus the fusion levels) can significantly determine the 3D shape of the spine at two-year follow-ups in the AIS population (Fig. 1).

Patient population.
A total number of 67 AIS with a main right thoracic curve (Lenke 1 and 2) age between 10-18 years were selected consecutively and retrospectively. All patients were selected from one hospital center, the Children's Hospital of Philadelphia, who underwent a posterior selective thoracic fusion with segmental derotation of the vertebrae. All patients had biplanar spinal stereoradiography at three visits: the X-ray images were registered within one-week pre-operation (PO), early (within a month) post-operation (EP), and at two-year after surgery (2Y). The exclusion criteria were previous spinal surgery, vertebral supernumerary, neuromuscular conditions, and musculoskeletal conditions other than scoliosis. Twenty non-scoliotic adolescents verified by The framework for 3D classification, path assignment, and outcome determination: (A) the 3D reconstruction of the spine from bi-planar X-rays is generated. (B) The vertebral centroids are calculated. The spline is generated by linearly interpolating the T1 to L5 spinal vertebral centroids 3D positions and is scaled in all directions (i.e., isotropically) such that the unit height is achieved. (C) K-means clusters are used to group the preoperative (PO) 3D scaled spinal curves. Three clusters are identified in our cohort (see Methods). (D) I. Fusion levels (F) are specified based on the position of the upper and lower fused vertebrae (UIV and LIV). Six fusion levels are identified in our cohort. II. K-means clusters are used to group the early post-operative (EP) 3D scaled spinal curvatures. Two clusters are identified in our cohort. (E) K-means clusters are used to group the 3D scaled spinal curves at two-year follow-up (2Y). Three clusters are identified in our cohort. The assigned treatment path for one patient has been shown with solid lines, indicating that PO 3 -EP 1 -F 6 path leads to 2Y 1 cluster at two-year follow-up (see Methods). The sagittal, frontal, and axial views of the PO, EP, and 2Y clusters are shown in Fig. 3. The 3D shapes of these clusters are shown in Fig. 4 bi-planar spinal X-rays and clinical examinations retrospectively were included as the control group. The institutional review board (IRB) at the Children's Hospital of Philadelphia approved the study procedures and the research was performed in accordance with the relevant guidelines and regulations. A waiver of consent/parental agreement was granted by the IRB.
Data collection and compilation. The 3D reconstruction of the spine was generated in a commercially available software, SterEOS 2D/3D (EOS imaging, Paris, France) for pre-and post-operative X-rays 25,26 . The 3D reconstruction of the vertebral bodies was used to calculate the 3D coordinates of the vertebral centroids (X, Y, Z) in the global coordinate system of the spine 27 . An isotropic scale factor, in three dimensions, was used to normalize the spinal height [0-1] of each patient in the AIS cohort ( Fig. 1). The same process, i.e., 3D reconstruction, center extraction, and scaling was performed for the control cohort. The Z levels of the average scaled spines in controls were used to interpolate the (X, Y) coordinate of the consecutive vertebral centroids in the AIS cohort, resulting in obtaining (X, Y) coordinates at equal Z levels for all the patients. This two-step normalization process was performed to scale all the spines in the Z-direction, eliminating to need to process the variabilities in the Z-direction and make comparison between the curves at the same Z-levels possible.

Patient classification.
A K-means clustering method 28 was used to cluster the 3D scaled scoliotic spines.
Given that the Z coordinates of the T1-L5 spinal vertebrae for all patients are the same (as described above) the clustering is performed only on (X, Y) coordinates of the vertebrae of all patients. The number of clusters was determined by calculating the silhouette values using the Euclidian distances 29 . The K-means cluster analysis was performed on the scaled PO, EP, and 2Y spines resulting in determining clusters of patients at the three time-points. The fusion levels i.e., the upper and lower instrumented vertebrae (UIV and LIV) were recorded. All fusion levels were determined and a group number was assigned to each distinct fusion level, i.e., same UIV and LIV. The treatment paths were then determined using the PO cluster number (PO i ), EP cluster number (EP j ), and fusion level group (F m ), constructing a three-level path presented as PO i -EP j -F m (Fig. 1).
A multinomial regression 30 model was used to predict the outcome cluster at two-year (2Y k ) from the treatment paths for the cohort as follows: where i, j, k are the cluster number at each given time-point and m is the fusion group. The number of patients in each treatment path that had the same 2Y k was calculated to determine the occurrence of certain outcome for each of the identified treatment paths. Finally, to determine the differences between the clusters in terms of clinically measurable variables, the clinical measurements of the patients at PO, EP, and 2Y were measured and compared between clusters at each time-point. These variables are the parameters that either can be directly measured on the spinal X-rays or via a commercial software (SterEOS 2D/3D, EOS imaging, Paris, France). The visual descriptions of each of these measurements are presented in Fig. 2. These clinical parameters are proximal thoracic Cobb (PTC), proximal thoracic rotation (apical) (PTR), main thoracic Cobb (MTC), main thoracic rotation (apical) (MTR), lumbar Cobb (LC), lumbar rotation (apical) (LR), thoracic kyphosis (TK both between T1-T12 and T4-T12), lumbar lordosis (LL), pelvic incidence (PI), sacral slope (SS), pelvic tilt (PT), frontal balance (FB), and sagittal balance (SB). Data normality was tested using the Shapiro-Wilk test and the sample size was determined for a test of comparison of the means between the clusters for a type I error of 0.05 and type II error of 90%.

Results
Three clusters are identified at PO i , i = 1, 2, 3. The sagittal, frontal, and axial views of the spinal curvature in each cluster are shown in Fig. 3(A). The EP j spines are grouped in two clusters, j = 1, 2, Fig. 3(B). The spines at 2Y k follow-up are clustered into three groups, k = 1, 2, 3, shown in Fig. 3(C). Figure 4 shows an illustrative case in 3D for each of these curve types at different time-points. Table 1 summarizes the average and standard deviation (SD) of the clinical parameters of the PO, EP, and 2Y clusters. The statistical methods and the significant levels are shown in Table 1. Significant differences were observed between PO clusters in PTR, MTR, and TK, between EP clusters in MTC, LC, TK, and SB, and between 2Y clusters in LC, LR, TK, and SB ( Table 1).
The fusion levels and their assigned group number are listed in Table 2. A total number of 8 distinct fusion levels are identified. Two patients were fused between T3-T11 vertebral levels and one between T2-T10 levels. Due to the small number of patients with these two fusion levels, these patients are removed from the analysis (total of 3 patients). The remaining 64 patients and 6 fusion groups are included in the following analysis.
Classifying the 3D spinal curves into a limited number of PO and EP clusters, along with identifying a handful of distinct fusion levels provide the opportunity for identifying a limited number of treatment paths. Figure 1 shows the framework for determining these treatment paths. A patient example is shown in Figure 1; the patient belonged to cluster 3 at PO, was fused from T4 to L1 (Group 6, Table 2) and belonged to the cluster 1 at EP, thus was assigned to the treatment path PO 3 -FE 1 -F 6 . The two-year outcome cluster as determined by the 3D spinal curve belongs to cluster 2Y 1 . The number of patients in each 2Y outcome cluster for all the existing treatment paths are reported in Table 3.
In order to address whether the treatment paths can significantly predict the 2Y clusters, a multinomial logistic regression was performed. The details of the multinomial logistic regression for association between the assigned treatment path (predictors) and the 2Y outcomes (predicted) are described in the Online Supplementary Information. Likelihood ratio tests of the regression model show significant association between the predicted and predictor variables, i.e., treatment paths and 2Y outcomes, χ 2 = 132.52, p = 8.3825e-11.
The number of patients who belonged to the same treatment path and 2Y cluster are calculated in Fig. 5. In instances where different fusion groups resulted in the same 2Y outcomes, the fusion groups are merged. Out of the total number of 17 existing treatment paths, a total number of 10 paths had more than two patients (Fig. 5). These 10 treatment paths include 56 patients. At least 70% of the patients who belonged to each of these 10 treatment paths had the same two-year outcome (Fig. 5). As an example, a total number of 6 patients belonged to the P 3 -FE 1 -F 6 path (path 17), out of which 5 (83%) were clustered in 2Y 1 . The treatment path of these patients is also shown in Fig. 1.

Discussion
Spinal fusion in AIS aims to realign the vertebrae in order to reduce the frontal deformity and curve rotation while imparting adequate sagittal profile 31 . New surgical maneuvers and instrumentation techniques allow for greater deformity correction of the curve in the three anatomical planes 32 . While much emphasis has been put on the pre-operative shape of the spine to guide the surgical decision-making, it has not been quantitatively determined whether a pathway which incorporates both the objectively identified pre-operative shape of the spine and the biomechanical changes in the spinal curve during the surgery can determine the long-term outcomes of the operation in AIS patients.
The 3D analysis of the spine has shown promising results in terms of determining subgroups of patients pre-operatively 33,34 . Given the importance of the 3D spinal curvature in postural assessment of the AIS patient 35,36 , a low-dimensional model that facilitates comparison between the 3D spinal shapes can improve the curve classification in AIS patients. The spinal evaluation in the three anatomical planes has underlined the parameters associated with the surgical outcome in AIS, which were not evident via the 2D evaluation of the spine 9,37 . Moreover, considering the interconnection of the pre-and post-operative 3D shape of the spine and the surgical outcomes, in the current study, we developed a quantitative pathway (patent pending) that can determine the outcome of the spinal surgery in AIS as a function of the 3D pre-operative spinal curve patterns as well as the modifiable factors during the surgery.
Here, to explore such pathways: (1) an objective classification of the 3D spinal curve was developed, and (2) the minimum key surgical parameters were identified. Regarding (1), the advantage of the clustering method based on the 3D shape of the isotropically scaled spines as opposed to defining the curves, based on a number of geometrical variables is reduced computational time and error. This allows recognizing the 3D patterns of the spinal curvature as opposed to the previous methods that extracted 3D variables of the curve for classification purposes 33 . Regarding (2), the 3D shape of the spine early after surgery in addition to the fusion level was used for characterizing the surgically induced changes in the spine. Using (1) and (2) in conjunction, a low-dimensional framework for assigning the most plausible outcome in terms of the 3D spinal curvature was introduced and evaluated in a group of 64 AIS patients with a thoracic curve. This approach determined the 3D shape of the spine in two-year follow-ups in 88% of the patients with an accuracy of at least 70%.
Selecting the fusion level is an ongoing debate in the surgical treatment of AIS [38][39][40] . The current clinical classification system and its modified versions attempt to address this issue 13,41,42 , yet complications and negative outcomes have challenged the surgical decision-making based on this guidelines 24,43,44 . Because the fusion level impacts the cost and time of the surgery 45,46 , patient's balance and posture after operation 39,47 , spinal flexibility, and disc degeneration 41,48 a quantifiable standardized tool for choosing the fusion levels is of critical need. Our proposed method lays out a pathway for outcome determination based on the shape of the spine before and after the surgery as the fusion levels changes (Fig. 1). As the clustering of the 2Y outcomes clearly distinguishes poor outcomes in term of the sagittal and frontal imbalances and residual curve at two-year post-operative ,   PO  PTC  PTR  MTC  MTR  LC  LR  TK T1-T12 TK T4-T12 LL L1-S1 PI  SS  PT  FB  SB 1 n = 24 6.6 ± 16.0 0.8 ± 0.5 57.6 ± 7.5 −8.4 ± 6.6 37.6 ± 6.5 8.1 ± 5.2 19.6 ± 10.6 § 7.6 ± 10.6 § 52.3 ± 11.4 54.7 ± 13.1 47.1 ± 9.1 7.6 ± 7.8 0.  ------EP  PTC  PTR  MTC  MTR  LC  LR  TK T1-T12 TK T4-T12 LL L1-S1 PI  SS  PT  FB   considering the treatment paths that results in a more favorable 2Y outcome can be achieved by selecting specific fusion levels (F). Comparing the clinical outcomes (2Y clusters) shows that achieving 2Y 2 at two-year follow-up is advantageous over clusters 1 and 3 due to the smaller change in the proximal junctional kyphosis, smaller residual lumbar curve, and both frontal and sagittal balances closer to the normative values 49 (Fig. 3), thus favoring the treatment paths ending to the 2Y 2 . As it is seen the fusion level can affect this outcome; for example, patients in clusters PO 3 and EP 1 can have 2Y 2 if they were to have fusion levels G3 (T4-T12) (Fig. 5), while other fusion levels (G6 (T4-L1) and G4(T2-L1)) resulted in a different 2Y outcome. The framework introduced here promises a surgical decision-making tool that can be personalized for physicians based on their preferred techniques. The rod length and curvature (using the EP clusters) can be determined prior to surgery to reduce the operation room time and number of staff which in turns reduces the surgical cost as well as blood loss and risk of infection 44 Table 3. in their satisfaction, self-image, and quality of life after surgery 52-54 , by exploring the possible outcomes, a more comprehensive expectation of the surgical outcomes can be achieved by the patient, patient's family, and their surgeon as they go through the process of surgery. This framework can be used as an assistive tool for the surgeons to modify several factors during the surgery in order to minimize the differences between the optimal and possible outcomes, similar to minimizing regret in decision theory 55,56 . Formulation of this decision-making problem with an embedded optimization approach will require additional factors that can possibly impact the outcomes.
The limitations of the current work include single-center and small patient dataset. This at some level limits incorporating various surgical technique, implant density, and variation in the implant including the material properties, rod diameter, and screw types. Additional surgical factors, which may vary between the surgeons and hospitals, remain to be evaluated via a larger multi-center study. While the impact of other surgical factors such as both magnitude and technique of the thoracic de-rotation, LIV alignment correction, number of osteotomies, and the reduction technique were not included, the early post-operative clustering of the spinal alignment in part accounts for these variables. A quantitative follow-up analysis of the spinal alignment as they relate to these surgical procedures is the subject of another study. Detailed analysis of the spinal curvature as considers the geometrical differences in the 3D spinal curve will be investigated in a larger database. Only patients with a right thoracic curve, the most common AIS deformity 49 , were included in this study. Application of a similar framework for surgical prediction outcomes in other types of spinal deformities and adult spinal surgical planning should be explored. As for any classification technique, the patients with borderline characteristics can adversely impact the outcome prediction. Yet, the K-means clustering classification in our cohort identified the outcomes with an accuracy of 70% or higher for 88% of the patients. More detailed clustering methods and its impact on the oucome prediction accuracy will be investigated. For a more heterogeneous cohort, Fuzzy clustering offers solution by assigning membership grades and further developing a probabilistic path for determining the outcome groups 57 . Despite the significant relationship between the intraoperative radiographs and the surgical outcome of AIS spinal fusion 58 , a study that associates the imparted intraoperative changes in the spine and the early post-operative alignment in a larger database can better guide surgical decision making.
The framework proposed here has an important implication in robotic surgery where a more quantifiable surgical goal can be established. Superimposition of both pre-operative and the aimed post-operative shape of the spine on the intraoperative images, not only can facilitate the screw placement but also guide imparting the required changes in the spinal profile in 3D. Furthermore, quantifying the trajectory of the changes in the spinal vertebrae position while quantifying the spinal flexibility through MRI 59 can closely calculate the required forces during the correction surgery. Quantifying these components can lead to a promising future for robotic surgery in pediatric spinal surgery.