This study consists of three cohorts: individuals with non-specific low back pain (NS-LBP, n = 43), individuals who are considering spinal fusion due to adult spinal deformity (SD-LBP; n = 42), and controls (n = 26) (see Table 1 for details on demographics). Patients were recruited during scheduled appointments for low back pain or pre-surgical consultation due to degenerative spinal conditions that result in thoracolumbar malalignment. The control cohort was recruited from the clinicians, staff, and students at the recruitment site. This research study was conducted in accordance with the Declaration of the World Medical Association. All subjects signed the informed consent (UCSF IRB: 16-21015, 17-22291).

Patient movements were collected in the clinic before or after their regularly scheduled appointments. Subjects were provided an armless chair 17.75 inches high, with a Kinect 2 (Microsoft, Redmond, WA, United States) depth camera placed 84 inches in front of them, 45 inches off the ground. Subjects were told to place their feet so that their shanks were normal to the ground. Subjects were asked to perform three self-paced sit-to-stand actions keeping their arms by their side. Subjects were allowed to perform less than three actions if they moved their feet, used their arms, or if the investigator or subject wanted to stop the experiment.

The Kinect body tracking library (C++) was used to obtain estimates of body pose from the color and depth data (Zhang, 2012). Estimates of joint position were streamed into a log file using custom software for post-processing at 30 Hz. The complete modelling approach used in this study is covered in previous work (Matthew et al., 2019b; Matthew et al., 2019a). A quadruple-pendulum planar model was used to model the movements of the ankle, knee, hip, and L5S1 joints. Segment geometry and inertia were scaled based on the subject sex, height, and mass using the relationships from (Mcconville et al., 1980; Young et al., 1983; Chaffin et al., 2006; Dumas et al., 2007). Inverse kinematics was performed using an unscented Kalman filter to obtain smoothed estimates of joint position and the corresponding joint angle (Julier and Uhlmann, 1997). Angular derivatives were obtained using sequential second-order, low-pass Butterworth filtering at 5 Hz, and numerical differentiation.

The filtered kinematic time series data were used to estimate the forces, and torques at each joint were computed using the inverse dynamics formulation by Park et al. (1995). The world-frame body segment positions, velocities, and accelerations were computed in a forward recursion from the computed joint angles, velocities, and accelerations. The contact forces and torques at each joint were then computed in a backward recursion, using the inertial parameters from the allometrically scaled planar pendulum model. Estimated joint powers were calculated from these torques and the corresponding joint angular velocities. Numerical integration was used to compute the concentric and eccentric work for each joint (Zhang et al., 2000). A planar muscular model was used to find the total loading at the L5S1 joint using the approach from Chaffin et al. (2006). This model estimates the effective muscle force and abdominal pressure from the computed torques and angles at the L5S1 joint. All measures that are derived from this analysis are normalized by the height and mass of the subject using the convention proposed by Pinzone et al. (2016). From our prior validation studies, this approach has been shown to have an average joint position error of 1.80 cm compared to active motion capture, and 30N average error compared to floor-mounted force platforms (Matthew et al., 2019a; Matthew et al., 2019b). The angle and torque time series data from the Kinect system were found to be concordant (0.81–0.99, 0.89–0.97) compared to our baseline active motion capture system.

The biomechanical set of variables included trunk and lower limb positions and angles (L5S1, hip, knee, and ankle angles measured with two different reference points: joint and world), velocities and accelerations (anterior and vertical velocities and accelerations calculated for each joint and body segment (torso, pelvis, thigh, and shank) with three different reference points: joint, body and world), dynamics (torques and powers calculated around each joint and normalized to patient’s height and weight) as well as compressive, shear and maximum muscle forces calculated at the L5S1 joint. The maximum and minimum variants of joint angles, accelerations, velocities, torques and powers, and maximum muscle forces around L5S1 were compiled for the data analysis. There were a total of 118 biomechanical derived variables included in the analysis. (see Data Dictionary, in Supplementary Material S1 (pages 1–10), which has a full description of each variable).

Figure 1 shows the schematic of the overall analytical workflow. The statistical analysis was performed using R (R: The R Project for Statistical Computing, R Project, 2021) and syndromic R package (Torres-Espín et al., 2021) specifically developed to facilitate post hoc principal component analysis (PCA), interpretation, and visualization of the PCA solutions. For the analysis, the biomechanics data were averaged across trials and repetitions of sit-to-stand as PCA suffers from intra-subject correlation and assumes independence of observations. Several studies on the kinematics of human movement comparing linear and nonlinear analytical methodologies have demonstrated that nonlinear methods outperform linear in the percent variance accounted for in the data (Harbourne et al., 2009; Portnova-Fahreeva et al., 2020; Boe et al., 2021). Furthermore, if the variables are linearly related, then PCA and NLPCA both result in the same solution. Therefore, to consider the likelihood of nonlinear relationships between the large number of biomechanical variables in our data, we have chosen to use a nonlinear approach.

Nonlinear PCA (NLPCA) was conducted by solving for the Gifi loss function through optimal scaling and alternating least squares as implemented in the princal () function of the Gifi package in R (Mair et al., 2019). All variables were ordinally restricted limiting to monotonic transformations. Nonlinear transformations were performed by a b-spline of degree two and three knots placed into the quartiles of the data. To determine the number of significant principal components to retain for further analysis, 1,000 permutations of the resulting NLPCA solution were performed using a permutation test (permD method with 1,000 permutations) with the overall goal to reduce the number of components while maximizing the variance accounted for (Torres-Espín et al., 2021). Based on the results of the permutation test the first 7 PCs had significant signal above random noise (Supplementary Figure S1 shows the results of the permutation test). The final number of 3 PCs were chosen for construct validation established by applying the following criteria: Kaiser rule (Kaiser, 1960), Scree plot (Cattell, 1966), and factor specification based on top PC standardized loadings (correlations between the vector defined by a variable and the PC with an absolute correlation coefficient ≥0.5). The stability of the NLPCA solution was assessed by performing bootstrapping on the approximations of the NLPCA solution generated through linear PCA on nonlinear transformations of the data. Specifically, 300 balanced bootstrap samples were randomly drawn and 300 PCA solutions generated. Pattern matching statistics (root mean square difference in PC loading patterns, the coefficient of congruence, the Pearson product moment correlation coefficient) were calculated for the first 3 PCs between the PCA solution and each one of the bootstrapped PCA, generating an empirical distribution for each metric for each PC. The average and non-parametric 95% confident interval for each metric was calculated (Supplementary Table S1 shows the bootstrapping results). Following construct validation which is an expertise-guided process of determining whether the pattern of variable association determined by the loadings conform to an explainable interpretation based on the current knowledge and understanding of these variables, the first 2 PCs are discussed as containing pertinent information (Figure 2). The full list of loadings is provided in the barmap plot in Supplementary Figure S2.

ANOVA was used to determine the between-group difference of the stable principal component scores as well as patient demographics (age, BMI, VAS scores, and Oswestry Disability Index (ODI)). Pairwise comparisons were conducted by Tukey’s post hoc test. Chi-square statistic was used to asses between-group differences in participant sex distribution. Linear regression models (lm () function in R) were used to assess the relationship between principal components (PC1, PC2, and PC3) and patient-reported outcomes of pain (Visual Analogue Scale, VAS) and ODI with a covariate adjustment for age and body mass index (BMI). The absence of collinearity between the PCs, age, and BMI was confirmed with the variance inflation factor (vif () function in the car R package (Fox and Weisberg, 2019). Statistical significance for all analysis was set at α = 0.05.

The analysis was performed on a convenience sample of patients with non-specific LBP, patients with adult spinal deformity, and healthy participants. There was no significant difference in male vs. female sex distribution across groups. There were significant between-group differences in age and BMI (control and NS-LBP groups only). Patients with adult spinal deformity had significantly higher VAS scores than patients with NS-LBP, however, there were no significant differences in ODI between the two patient groups. Neither VAS nor ODI measures were assessed in healthy controls, therefore, the control group was not included in the patient-reported outcome analysis (Table 1).

Unsupervised machine learning on full body biomechanics (Figures 1A,B) identified divergent movement strategies in LBP patients compared to healthy controls (Figures 2). The NLPCA workflow consists of multiple steps starting from raw participant-level data curation (Figure 1C), which then is subjected to optimal scaling transformation (Figure 1D), unsupervised pattern detection (Figures 1E,F), and expert-augmented construct validation (Figure 1G). This helps identify biomechanical parameters related to pain, enabling us to cluster participants based on their movement phenotypes (Figures 1H,I). Two global movement strategies represented by the principal components were deemed to capture a significant proportion of the variance (see methods for details).

The first principal component (PC1) accounted for 36.8% of variability within the dataset (Figure 2A). PC1 captured three different components of STS associated with either the control participant strategy (positive loadings) or patient movement (negative loadings). Maximum torso, pelvis, and thigh kinematics (accelerations and velocities) in the world and body frame, as well as peak L5S1, hip, and knee joint kinetic variables had significant positive correlations with PC1. World thigh maximum flexion angle, L5S1, hip minimum flexion angles, and knee minimum flexion torque had significant negative correlations with PC1 (Figure 2A). Based on this loading pattern, PC1 was regarded as a “vertical rise” STS strategy (positive loadings = full body vertical kinetics) and final standing posture (negative loadings = joint/segment kinematics).

There was a significant group difference in PC1 scores (F (2, 108) = 7.55, p =0.0008) (Figure 1B). Specifically, control participants had significantly greater PC1 scores as compared to SD-LBP group (p = 0.0005). There were no significant differences in PC1 scores between NS-LBP and controls (p = 0.066), and the two patient groups (p =0.148). Taken together these results suggest that lower full body kinetics and forward flexed standing posture are the main differentiating/defining features of SD-LBP patient STS movement when compared to healthy control participants.

The second principal component (PC2) accounted for the next 19.3% of variability in the data (Figure 1C). PC2 loadings featured several kinematic, kinetic, dynamic, and muscle force variables and for interpretability purposes, the loadings were grouped based on their biomechanical identity as following. Positive loadings included: world torso and pelvis maximum flexion angle, L5S1 and hip maximum flexion and maximum sagittal vertical axis (SVA) characterize body flexion or forward-leaning; thigh and shank maximum world anterior acceleration, shank maximum world anterior velocity, shank maximum body anterior and vertical velocities altogether encompassing lower limb forward kinetics; L5S1 lumbar anterior shear and compressive maximum forces, L5S1 sacral compressive maximum force and L5S1 extensor muscle maximum force which characterize spine loading; L5S1 and hip maximum flexion powers and L5S1 and hip maximum extension torques representing total effort; pelvis and torso maximum body posterior accelerations and velocities, thigh maximum world posterior velocities capturing posterior body swing (Figure 2C). The variety of biomechanical modules captured by PC2 comprises a leaning stand-up strategy and its biomechanical consequences (e.g., overall effort and spinal loading). There was a significant group difference in PC2 scores (F (2,108) = 30.091, p <0.0001) (Figure 2D). The control group had significantly lower PC2 scores as compared to NS-LBP (p =0.0001) and SD-LBP group (p <0.0001). NS-LBP has significantly lower PC2 scores as compared to SD-LBP (p =0.0004). The positive PC2 scores in the patient groups suggest that the positive PC2 loadings are more descriptive of the patient participants’ STS strategy. Thus, patients utilize a leaning forward–swing back STS strategy that is more laborious (requires greater extension torques) and more strenuous (results in higher spine loadings) (Figure 2E). Overall, NLPCA has separated the control and 2 patient groups into three overlapping but distinct clusters across PC1 and PC2 space (Figure 2F).

The third principal component (PC3) captured 10.4% of variability in the data. Kinematic, dynamic as well as kinetic variables have loaded into PC3 in a less cohesive fashion, capturing the remaining STS movement components. Positive loadings captured downward or caudal velocity of the pelvis and thigh body segments, the resulting flexion torque around the L5S1, knee and hip joints as well as world torso and pelvis maximum extension angles (leaning back). Negative loadings captured lower body maximum kinematics (acceleration and deceleration). Overall, these movement components do not link or group into a strategy as do PC1 and PC2 limiting the interpretability and establishment of its construct validity (Supplementary Figure S3A). Furthermore, there were no significant group differences in PC3 scores (F (2,108) = 1.78, p = 0.17) (Supplementary Figure S3B); therefore, they were excluded from further analysis.

The construct PC4–PC7 were excluded from the main analysis following Kaiser and Scree plot rules. However, based on the permutation test PC4–PC7 have significant loadings/information. P7 did not contain any loadings above 0.5 and was excluded from further investigation. The loadings of PC4–PC6, which together account for additional 16.1% variance in the data and capture the information related to participants ankle biomechanics, are reported in Supplementary Figure S4.

Linear regressions were used to determine whether biomechanical parameters and/or compensatory movement strategies captured across PC1 and PC2 can predict patient-reported outcomes including pain (VAS) and disability (ODI). We found a significant association between PC1 and VAS, as well as PC1 and ODI (Figures 3A–C) when controlling for BMI and age by including them as co-factors. The higher PC1 scores corresponding to greater speed of movement had a negative correlation with VAS and ODI scores, suggesting that patients with worse pain and greater disability move slower. In contrast, there were no significant associations between PC2 and VAS or PC2 and ODI (Figures 3D–F).

There are several important limitations in this study. The analysis was performed on a convenience sample, in which age and BMI were not controlled. The two patient groups were significantly older and had higher BMI as compared to control subjects. Aging is associated with increased BMI, as well as overall reduced fitness, including muscle capacity. To control for the potential confounders, age and BMI were included as co-factors in the linear regressions of PC scores with VAS and ODI outcomes, and no significant contribution of either age or BMI was found on the significant correlation between PC1 scores with VAS and ODI, suggesting that pain and disability are more strongly related to the observed movement phenotypes than either age or BMI. In addition, LBP is a multidimensional heterogeneous disease and, therefore, patient phenotyping using even the most comprehensive biomechanics measures does not fully address the multifactorial nature of LBP. Deep patient phenotyping should be ideally performed on many outcome measures, including comprehensive patient-reported outcomes collected through validated NIH questionnaires (e.g., PROMIS), muscle electromyography, neuroimaging, quantitative sensory testing, etc., that capture the different aspects of the biopsychosocial model of LBP.