Thirty-eight adolescent subjects [mean age 14 (SD 2); 27 females and 11 males] suffering from mild idiopathic scoliosis (CAs <24°) underwent radiographic examination in orthostatic position at IRCCS Istituto Ortopedico Galeazzi (Milan, Italy). The subject assent and the parental permission to use the anonymized radiological data were given by signing an informed consent approved by the local ethical committee. Digitized images of the thoracolumbar spine and pelvis were simultaneously acquired in coronal and sagittal planes with EOS imaging system (EOS Imaging, France). Since EOS system provides spatially calibrated images (Illés and Somoskeöy, 2012), no further calibration procedures were required. The images pair was manually processed through SterEOS software (EOS Imaging, France), which allowed to identify as “scoliotic,” the spine curves characterized by CAs >5°, and provided in addition the identification of the anatomical parameters. The following indexes were extracted for each subject: (i) number of scoliotic curves with related CA (in case that more than one curve had been identified, the most severe curve and related CA were chosen for the analyses); (ii) sacral slope (SS); (iii) pelvic incidence (PI); (iv) lumbar lordosis (LL); (v) thoracic kyphosis (TK) from T1 to T12; and (vi) maximum vertebral rotation in the axial plane (MAR). In addition, the Roussouly Type (RT) for the classification of the lumbar and pelvis sagittal alignment was manually determined (Roussouly et al., 2005).

The coronal and sagittal pair of biplanar images was manually processed with in-house script running with MATLAB (MathWorks Inc., Natick, MA, USA). Right-hand global reference system {x,y,z} was adopted considering axis “x” craniocaudally oriented, axis “y” posteroanteriorly, and axis “z” from right to left facing coronal plane (Figure 1). Since EOS system allows for simultaneous acquisition of true to size images in one-to-one scale, each point of the 3D space results jointly projected on the coronal and sagittal images (Figure 1). Accordingly, the identification of appropriate landmarks over the biplanar images allows for the calculation of geometrical parameters of spinal configuration of each vertebra in the 3D space using (i) location; (ii) dimensions; and (iii) rotation around the three axes. The Section “Appendix” explains in detail the reconstruction procedure. The 12 thoracic vertebrae (T1, …, T12), the 5 lumbar (L1, …, L5), and the first sacral (S1) were assessed (Figures 2A,B) obtaining the corresponding 3D geometric model (Figures 2C,D). It is worth noting that although visualized through a standardized set of reference anatomical vertebral meshes, the geometrical subject-specificity of the model is guaranteed by reconstructing vertebral location, dimensions, and rotations. The accuracy of the spine reconstruction method was previously checked by comparing the computed vertebral orientations with those obtained by SterEOS proprietary software, known suitable in providing accurate reconstruction of the spine (Glaser et al., 2012). The differences in reconstructing the orientation angles were found comparable with those presented by other authors when evaluating different reconstruction techniques (Pomero et al., 2004; Kadoury et al., 2009; Moura et al., 2011). Specifically, the difference in the thoracic section was 1.9° ± 0.4° (mean ± SD) for coronal rotation, 2.3° ± 0.8° for lateral rotation, and 3.7° ± 1.1° for axial rotation. Concerning the lumbar region, similar values were pointed out as 1.8° ± 0.6° for coronal rotation, 2.4° ± 0.9° for lateral rotation, and 3.6° ± 1.4° for axial rotation.

Musculoskeletal modeling characterizes bones as rigid segments connected by joints and muscles as tensile elements attached to segments. Through an inverse dynamic approach, muscles forces and intersegmental forces acting during the execution of specific imposed kinematics and under the action of known external loads are computed by minimizing muscles recruitment activation (Rasmussen et al., 2001; Damsgaard et al., 2006). In the present work, the 3D musculoskeletal model was processed with AnyBody software v.6 (AnyBody Technology, Denmark). Geometric properties describing vertebral location, dimension, and rotations were set in AnyBody through MATLAB-based script, importing the geometric parameters as reconstructed in Section “Geometric 3D Reconstruction.” Accordingly, the spine mechanical model consisted of the thoracic and lumbar region, as well as the sacrum (Figures 2E,F). Sacrum was constrained to ground. Since the model was aimed to an inverse static simulation, vertebral mass and moments of inertia were set to 0. Conversely, mass of the body related to the trunk weight of the subject was distributed along the whole spine according to the literature (Kiefer et al., 1997). Weight distribution proposed by Kiefer et al. was scaled to match the subject weight thus improving the subject-specificity of the model. Weight forces were applied anteriorly to vertebrae, whereas weights related to head and arms were applied to T1 upper end-plate and anteriorly to T4, respectively (Kiefer et al., 1997). A total of 89 muscles and their attachment sites were identified according to physiologically–anatomically appropriate values (Stokes and Gardner-Morse, 1999) and incorporated into the spine model [12 muscle elements for thoracic multifidus (MF), 20 lumbar MF, 5 longissimus pars lumborum, 4 iliocostalis (IC) pars lumborum, 12 longissimus pars thoracis, 8 IC pars thoracis, 11 psoas, 5 quadratus lumborum, 6 external oblique, and 6 internal oblique]. Muscles were modeled as single force components exerting only tensile forces. No force–length relationship and no force–velocity relationship of muscles were modeled, since the present work evaluated exclusively static postures. Muscle wrapping around vertebral surfaces was not allowed by the current AnyBody release. Intervertebral disks were modeled as spherical joints with three rotational degrees of freedom, located at the midpoint of the segment connecting vertebral centers. Two adjacent vertebrae result thus connected by the joint through rigid subsegments starting in vertebral centers and ending in the joint. AnyBody force-dependent kinematic (FDK) analysis method minimizes muscles recruitment to calculate the intersegmental forces acting in maintaining the assigned configuration. In addition, the computed forces are iteratively processed to calculate joint translational deformation (i.e., shifts of the subsegments endings) via strain-displacement linear relation with imposed joint stiffness as weighting coefficient. As summarized in Table 1, joints’ translational stiffness interpreting soft tissues, i.e., disks and ligaments, and apophyseal joints’ role was set according to Panjabi et al. (1976) describing thoracic segments and according to Gardner-Morse and Stokes (2004) evaluating lumbar segments in preload condition (i.e., 500 N axial compressive) to describe orthostatic position. Rotational joint stiffnesses were not set since not accountable for by the current AnyBody release. The spherical joints provide thus internal reaction forces but not reaction moments. No tendons or other passive element properties were modeled. Muscle recruitment was set in order to minimize muscle activation via a polynomial optimization criterion, according to optimal delays in the muscle fatigue while maximizing muscle synergism (Rasmussen et al., 2001). It is worth noting that in the proposed modeling approach, the subject-specificity was accounted from geometrical and weight-distrubution points of view. Indeed, the vertebral location, orientation, and dimension were reconstructed from the radiographic images, and the subject weight was appropriately distributed along spine.

For each joint, the intersegmental force F, acting on the caudal vertebra was evaluated and measured in global reference system (F_x, F_y, F_z), representing F_x axial compression, F_y posteroanterior shear, and F_z lateral shear. Before being compared with other relevant variables, the resulting F values were normalized by the corresponding subject weight force in order to avoid potential biases related to the differences in the intersubjects’ weight. Possible correlations between F computed at L4L5 (F^L4L5) and at L5S1 (F^L5S1) with SS, PI, LL, and RT were searched. The F of spinal level corresponding to the scoliotic curve apex (F^A) and the maximum of each component of F (F^M) found along the scoliotic curve were checked for correlations with CA, TK, and MAR. Pearson correlation coefficient or Spearman rank correlation coefficient in case that a normal distribution was not achieved was accounted for in the comparisons. Both coefficients range from −1 to 1, where 0 indicates null linear correlation and −1 and 1 indicate full negative and positive linear correlations, respectively. Statistical significance of the coefficients was tested according to two-tailed t-test or permutation distribution testing assessing Pearson and Spearman coefficients, respectively, considering 0.05 as significance level.

Being {x′,y′,z′} the local vertebral coordinates system, the F_x′ and F_y′ acting on T12 and on lumbar vertebrae from L1 to L5, and the muscle forces calculated with the proposed model for a subject with no scoliotic curves (CAs <5°) were compared with the corresponding values computed by El-Rich et al. (2004) with a symmetric spine model representing healthy subject. The compared subject (female, 16 years, 154 cm of height, 47 kg of weight) was chosen as the most appropriate non-scoliotic subject suitable to guarantee comparable total weight distributed force in the two models (i.e., 345 and 387 N, respectively). The height of the spine model (from T1 to the sacrum) in the evaluated subject and in the model processed by El-Rich et al. resulted 39 and 47 cm, respectively. Unlike El-Rich et al. approach, in the present model the spine symmetry was not guaranteed, and muscle forces to be compared were obtained as the average of the left and right specific muscles values. According to the setting of El-Rich et al., the comparison was performed evaluating two different conditions: (i) gravity alone and (ii) adding to T3 a 380 N gravity-oriented load, accurately placed according to El-Rich et al., simulating a weighted bar held in front with arms extended in gravity direction close to the body. It is worth considering that there is a substantial difference between the proposed approach and that of El-Rich et al. in modeling intervertebral disks since the former defines them as deformable spherical joints, while the latter as flexible beams. According to that, the comparison of disks’ reaction moments provided by El-Rich et al. was not feasible since differently from the flexible beams the spherical joints are not able to provide reaction moments.

When comparing disk pressure computed by the present model with in vivo measurements, for every subject the pressure at the specific level was calculated as the ratio between the intersegmental axial load, F_x′, acting on the caudal vertebra and the upper end-plate area of the caudal vertebra inferred from the geometrical reconstruction (see Appendix). Since recognized from in vitro studies that the pressure measured in the nucleus results higher than the average disk pressure (force divided by total cross-sectional area) by a factor of 1.54 (Brinckmann and Grootenboer, 1991; Nachemson, 1960), the computed pressure was accordingly corrected before being compared. Pressure at L4L5 was compared with in vivo measurements performed by Wilke et al. (2001) and by Sato et al. (1999). Thoracic pressure at T6T7, T7T8, T9T10, and T10T11 was compared with in vivo measurements obtained by Polga et al. (2004).

Intersegmental loads acting on T12 and in lumbar spine were compared with that computed by El-Rich et al. (2004) in two conditions: gravity alone and 380 N held in front load. As reported in Table 2, comparable values were pointed out assessing axial load, F_x′, and posteroanterior shear, F_y′. However, in the present model the maximum F_x′ during the two conditions were found at L1 level instead progressively increasing from T12 to L5 as reported by El-Rich et al. This finding is confirmed evaluating the axial load in global coordinates (see F_x in the lateral panel of Figure 3) describing the same subject in gravity alone condition. F_x shows a peak at the thoracolumbar junction between T12 and L1, while larger value exceeding that peak was pointed out at the L5–S1 lumbosacral junction only (Figure 3). We hypothesize that this result can be related to the distribution of the muscles’ forces necessary to achieve equilibrium, which depending on the muscles configuration can potentially promote an increase of the axial load on T12–L1 junction. Further investigations are necessary to better clarify this outcome.

Muscle forces were found generally lower for gravity alone condition and mainly larger during 380 N load with respect to El-Rich et al. The larger values found in 380 N load can be explained considering that El-Rich et al. evaluated the model in an optimized posture aiming at minimizing loads, different from that assessed during gravity alone. Conversely, postures were kept identical in the present model. It is thus expected that a non-optimized model would exhibit higher muscle activations. Furthermore, perfect matching in the muscle comparison was not viable since muscles’ insertion placements were not identical in the compared models. In addition to muscles connecting the pelvis to thoracic and lumbar vertebrae, the present model accounted for intervertebral muscles (i.e., additional IC and longissimus pars thoracic, thoracic MF and lumbar MF muscles) that were not modeled by El-Rich et al., which can induce differences in the muscles forces distribution. Moreover, El-Rich et al. accounted for the entire T1–T12 segment as a rigid body, whereas the present approach assessed the thoracic vertebrae as separated bodies connected by spherical joints.

The calculated disk pressures were found comparable with that obtained in vivo in lower thoracic section (Table 3) thus supporting model suitability. However, pressures calculated in lumbar and upper thoracic sections were found 60% larger and 50% lower, respectively, in comparison with the in vivo measurements (Table 3). Although these discrepancies could be interpreted as effects of spine deformation characterizing scoliosis, further investigations are necessary to guarantee comprehensive model validation. For example, assessing the present approach in unload configuration, comparing pressures with the in vivo values measured by Meir et al. evaluating scoliotic patients in supine position during reconstructive surgery (Meir et al., 2007).

The orientation of the intersegmental load, F, was found to be in principle aligned with spine curvature for every subject (Figure 3). This finding is in agreement with the follower load concept discussed in previous works (Patwardhan et al., 1999; Shirazi-Adl and Parnianpour, 2000) and supports the hypothesis that a load path tangent to the spine curve better sustains compression load (Patwardhan et al., 1999). However, this result can be mainly related to neglecting rotational joint stiffness because it is not allowed by AnyBody. Future developments are thus necessary to assess F orientation when adding intervertebral structures able to interpret passive rotational contributions (i.e., elastic elements).

For every subject, the incremental tendency of F_x in the craniocaudal direction, related to the progressive increment of body weight to be sustained, exhibited a peak between T12 and L1 (right panel of Figure 3). Since body weight forces act along “x” axis and represent the solely external forces to be balanced in maintaining posture, the F_x discontinuity can be related to the distribution of the muscles forces necessary to achieve equilibrium, which depending on the muscles’ configuration can potentially promote increased joint axial load on thoracolumbar junction. As stated in Section “Preliminary Model Validation” discussing model validation, further inquiries are necessary to better clarify this outcome.

The significant correlations found comparing SS, PI, and LL with the posteroanterior shear in L4L5 and L5S1 joints (FyL4L5 and FyL5S1, Figure 4A) support the suitability of the proposed model. The result is indeed in agreement with previous works (Galbusera et al., 2014) and suggests that the higher is the slope of sacrum and pelvis in sagittal plane, the more lumbosacral loads result anteroposteriorly oriented (i.e., showing negative values) in maintaining standing posture. That relation is verifiable observing the negative correlation between PI and FyL5S1 expanded in Figure 5A and can be elucidated through the example in Figure 3 where the F load acting on S1 in the sagittal projection appears anteroposteriorly oriented, and thus numerically negative. The correlation found between RT and FyL4L5 and FyL5S1 (Figure 4A) supports the relation identified between lumbosacral loads and anatomical sagittal parameters. By definition, Roussouly classification is indeed obtained evaluating SS together with lumbar curve characterization through the recognition of lumbar curve apex in sagittal plane (Roussouly et al., 2005) thus linking RT to both sacral and lumbar sagittal peculiarities.

The apical and the maximum loads inside the most severe scoliotic curve (F^A and F^M) were compared with the anatomical parameters in the three anatomical planes (i.e., CA related to coronal plane, TK to sagittal, and MAR to axial). In coronal plane the CA index, which measures scoliosis severity, was found not correlated either with F^A or F^M (Figures 4B,C). This unexpected lack of correlation can be related to mild scoliosis condition, not able to induce significant increasing in lateral loads. Further developments evaluating severe scoliosis are thus required to verify this hypothesis. In sagittal plane, TK index was found correlated with apical and maximum posteroanterior shears (FyA and FyM in Figures 4B,C). In that case, as expanded in Figure 5B, the more TK was pronounced (i.e., larger TK values) the more FyA resulted posteroanteriorly oriented. This finding confirms the relation between sagittal anatomical parameters and intersegmental posteroanterior shears, supporting the relation observed in Section “Correlating Lumbosacral Loads with Sagittal Parameters” between sagittal parameters and lumbosacral posteroanterior shears. In axial plane, the negative correlation found between MAR and maximum lateral shear FzM (Figure 4C) is expanded in Figure 5C. This relation is explicable by observing the spine curvature depicted in the example of Figure 3. Although not scoliotic (CA <5°), this spine curve is characterized by positive axial rotations of L1 and L2 vertebrae (i.e., right-hand rotations around craniocaudal “x” axis) providing negative F_z loads (right panel of Figure 3) oriented toward curve concavity. As inferable by Figure 5C, the larger was MAR the larger resulted F_z oriented concavely to the curve.

Considering that a more severe scoliotic deformation is usually associated with larger vertebral axial rotations (Stokes et al., 1987), the significant correlations pointed out in the three anatomical planes corroborate previous reports about the importance of the 3D aspects of scoliotic deformity (Graf et al., 1983).

The proposed modeling approach provides the following advances with respect to the state-of-the-art: (i) offers a non-invasive method to assess intersegmental loads, based on radiographic images evaluation; (ii) can be exploited to evaluate intersegmental loads in case of spine deformities, e.g., scoliosis, since characterizing each vertebral level of thoracolumbar spine; (iii) the reconstructed mechanical model is subject-specific from geometrical and weight-distrubution points of view, since vertebral geometrical location, orientation, and dimension are reconstructed evaluating the radiographic images, and subject weight is appropriately distributed along spine; (iv) the deformation of the intervertebral disks was accounted for through the AnyBody FDK method, which assesses force-dependent translational displacements by setting joint stiffnesses.

The present modeling approach has several limitations, which are as follows: (i) co-contraction of trunk muscles, known to be a physiological strategy to enhance spine stability (van Dieen et al., 2003), was neglected. Muscles were modeled as pure forces, not considering relevant aspects such as wrapping (Arjmand et al., 2006) and curved courses (Bazrgari and Shirazi-Adl, 2007); (ii) motion segments were modeled as kinematic joints with linear stiffness, thus neglecting real shape of the structure, the presence of posterior elements, and the non-linearity of biological soft tissues mechanics. However, the present model was designed to maintain standing posture thus avoiding the necessity of accounting for non-linear flexibility response of spinal segment (Abouhossein et al., 2011). The joint stiffness values chosen from the literature were measured in healthy young and middle-aged subjects without disk injuries or degenerative factors (Panjabi et al., 1976; Gardner-Morse and Stokes, 2004) since values from adolescents and scoliotic subjects are not available; (iii) rotational joint stiffness was neglected since not allowed by AnyBody. Further developments should include rotational parameters by purposely introducing intervertebral elastic elements; (iv) the description of spino-pelvic configuration was simplified to include only SS neglecting PI and hip axis; (v) the positions of muscle insertions were based on literature data and assumptions, which took into consideration only some aspects of patient-specific anatomy; (vi) rib cage modeling known affecting thoracic stiffness properties was neglected. Intra-abdominal pressure, recognized to play potential role on lumbar spine, was not accounted for since not expected to play significant role in relaxed standing.

In conclusion, the comparison with previous modeling works and in vivo studies partially fulfilled the preliminary validation purpose, but further investigations are necessary to clarify the observed incongruities and to guarantee comprehensive validation. When applying model approach to evaluate the relation between intersegmental loads and spine anatomical parameters in mild scoliotic subjects, encouraging results supporting model suitability were pointed out. Despite a number of limitations that suggest to be prudent and that will be overcome in future developments, the present method appears to be a promising tool. Once fully consolidated, it can allow the subject-specific non-invasive evaluation of a deformed spine, providing supplementary information to the routine clinical examination and supporting the surgical intervention planning.