Reducing the Brace Correction Stress on the Secondary Lumbar Curve Results in Excellent Muscle, Bone, and Disc Mechanical Performance: A Musculoskeletal Finite Element Simulation of AIS Patient With Rigo A3

2.1. Data Acquisition

A 13‐year‐old female patient with Rigo A3 scoliosis, a height of 152 cm, and a weight of 42 kg was included in the study. A Siemens dual‐source computed tomography (CT) scanner (Siemens Somatom Definition 2012B; Siemens Healthineers, Erlangen, Germany) was used for transverse scanning of the entire spine and pelvis, covering T1 through S1, with a layer thickness of 0.7 mm. All CT images were saved in DICOM format.

2.2. Building the Three‐Dimensional Musculoskeletal Model

The DICOM format CT images were imported into Mimics 21.0 software (Materialize, Leuven, Belgium) [18]. The model coordinates were defined as follows: coronal axis, X; sagittal axis, Y; and horizontal axis, Z. The cervical, thoracic, lumbar, sacral, and pelvic vertebral regions were identified based on the gray values of different tissues. The identified bones were imported into Geomagic Studio 2017 (Geomagic Inc. Research Triangle Park, NC, USA) in STL format to process the uneven surface of the medical image accumulation (Figure 1A). Each surface component was imported into SolidWorks 2020 (SolidWorks Corporation, MA, USA) to create a solid individual model fitting for the spine‐pelvis assembly [18]. The models included the articular cartilage, cortical bone, cancellous bone, fibrous annulus, nucleus pulposus and endplate, and muscle. The attachment positions of muscle bundles and bones in the FEM were added based on the classic skeletal muscle model. The muscle was idealized based on previous reports, and the physiologic cross‐sectional area of the skeletal muscle was determined to establish the trunk muscle model [21, 22, 23, 24, 25, 26, 27, 28, 29]. The curve of the whole extensor muscle [30] and the curve of the internal and external oblique muscles of the abdomen were simulated using the elliptical torso model [27]. Finally, the models were combined and saved in the SLDPRT part format.

Ansys Workbench 19.0 was used to perform the finite element analysis (ANSYS Inc. Cannon Sburg, USA). To ensure the accuracy and validity of the model and optimize the demand for computing resources, the solid model of each vertebra was divided into a high‐quality mesh using the adaptive dynamic mesh of the modeler. Bone, cartilage, and other tissues, including the anterior longitudinal, posterior longitudinal, Ligamenta flavum, supraspinous, and interspinous ligaments, were specified as the tetrahedral solid units of a network using Ansys Workbench 19.0. Spring elements were used to simulate ligaments between the vertebrae (Figure 1B). All components were defined as elastomers, and the bone structure and discs were modeled as isotropic elastic materials. The material properties of bone, IVD, and ligament were determined based on experimental literature [16, 17, 21, 31, 32, 33, 34, 35, 36]. The material parameters of bone and soft tissue are summarized in Table 1 [15]. Different materials can be distinguished according to their differences in elastic modulus and Poisson's ratio.

To observe the important effects of muscles and ligaments on the stability and mobility of the spine, a muscle ligament model was constructed. Muscle forces and reaction moments were estimated in the upright‐standing posture under load gravity simulation, which was constrained to equilibrium and stability conditions, and the scoliotic curvature and stress distribution were calculated (Figure 1B). According to Hueter–Volkmann's law, increased pressure on the concave side of the epiphysis and increased tension on the convex side change the growth rate of the vertebral body, manifested by concave inhibition and convex acceleration [37]. Loading on the spine includes gravitational forces acting on the vertebrae and muscle forces that maintain spinal stability [38]. Based on previous studies [39, 40], we assumed that the load on the vertebral body is closely related to the body's self‐gravity. The musculoskeletal FEM of the AIS‐simulated forces in the upright‐standing position of the patient is shown in Figure 2. The actual weight of the patient with AIS was 42 kg. The simulated forces act along the Z‐axis at the center of the vertebral body, constraining all degrees of freedom of the sacrum. The simulated forces allow T1 forward flexion, lateral flexion, and vertical translation to maintain balance. The equation used to calculate the loads on the patient at each vertebral level [20] is shown below.

where g = 9.81 m/s; Massi, mass corresponding to the ith vertebral body level; Ptoti, cumulative mass of the part above the level of vertebral i.

2.3. Mesh Sensitivity Analysis

Sensitivity analysis was performed after each parametric FEM was generated. The reason for the analysis is to determine the correct size of the part, as well as the type of part used and its numerical formula to ensure excellent results at the lowest computational cost [41, 42]. According to the previous study on the FEM of scoliosis, Ansys Workbench 19.0 was used to divide tetrahedral elements. We selected T7 as the research object to divide the number of four groups of mesh and conduct axial load simulation analysis [43]. It can be seen from Table 1 that the number of meshs has a corresponding impact on the results.

When the number of meshs is < 20,215, the difference of results is obvious. With the increase of the number of meshs, the stress–strain of the vertebrae divided by the number of adjacent meshs also decreases. The difference between the number of grids D and the number of grids C is the smallest, The difference between stress and displacement is 4.37% and 0.05%, respectively. Therefore, it can be considered that Group C satisfies the verification.

Therefore, it can be considered that in Group C, the meshes number of the vertebral body is 29,308, which meets the verification. The meshs quality check criteria are as follows: aspect ratio ≤ 5°, maximum Angle ≥ 30°, minimum angle ≤ 120°, Jacob ratio > 0.7, warpage ratio ≤ 0.6.At this time, the meshs of vertebrae and the whole torso is shown in Figure 3A. The overall torso meshes model is shown in Figure 3C.

2.4. Loading and Boundary Conditions

The intermodel contacts were defined based on the structure and motion characteristics of the spinal–pelvic region. The types of contact between models were defined in the connections. The contact type of the joint surfaces was defined as no separation, which does not allow the surfaces of the contact area to separate but allows small frictionless sliding along the contact surface. High pressure on the body causes discomfort and pressure ulcers. A reasonable stress range and area were selected to ensure acceptable local pressure generation by the body brace. Modifications due to human gravity and muscle forces were comprehensively considered, and correction forces within the range of the human pain threshold were guaranteed. Therefore, the corresponding ceiling thresholds for the area and correction force were set according to the model previously reported in this study (Table 2) [44, 45].

2.5. FEM Validation

Number of grids influences the results of the numerical calculation. Improving the grid quality gradually reduces the performance prediction error [46]. When grid increases have little influence on the calculated results, the numerical results are meaningful [17]. The model was meshed to ensure the accuracy of the calculations. The grid type was a tetrahedral mesh, and the grid size of the joint surfaces was set to 0.8 mm; the other grids were set at 1.5 mm. Therefore, the model was divided into 1,353,236 elements and 2,381,761 nodes.

The biomechanical properties of the FEM were validated to ensure the reliability of the mechanical response analysis during the brace correction simulation [47]. The geometry, range of motion (ROM) of the vertebral segments, and muscle axial force of the model were validated. To verify the geometric similarity of the model, the spinal imaging parameters of the model were measured after gravity loading (Figure 4). The model was adjusted to the coronal plane, and the Cobb angles of the main and near chest bends were measured using the x‐ray sheet Cobb measurement method. The model was adjusted to the sagittal plane to measure the thoracic kyphosis (TK) and lumbar lordosis (LL) angles. The model was adjusted to the axial position, the vertebra T2 and T7 were located, the respective axial vertebral rotations were measured, and the measured values were compared with the parameters from the patient's EOS imaging. The numerical difference between this model and the clinical imaging parameters was < 5°. Therefore, this model exhibited geometric similarity with the clinical data (Table 3).

The ROMs of the T1–T8 joints in this model were previously verified. All degrees of freedom were constrained on the lower surface of the T4 vertebra, and 4 Nm torques were applied to the upper surface of the T1 vertebra. Forward flexion and backward extension, left and right lateral flexion, and axial rotation were simulated to verify the T1–T4 segment ROM. Similarly, the ROM of T5–T8 in the model was verified [48, 49, 50]. The T1–T4 segment ROM simulation results were verified (Figure 5A); the stimulated ROM values in this model were consistent with anterior flexion, dorsal extension, left/right flexion, and left/right rotation ROM values from previous research models [48, 49, 50]. However, structural changes near the chest bend differed in this model; the left/right ROM values were slightly lower than the ROM values in the three preceding models but were within an acceptable range. Verification of the T5–T8 segment ROM (Figure 5B) revealed that the anterior flexion values and left/right rotation ROM simulated in the selected segments of this model were higher than those values obtained by Busscher et al. [48] and Xin et al. [50]. The location of the top of the main thoracic (MT) curvature, the different degrees and types of lateral curvatures, and the older age in the previous models may make the ROM values of this model slightly larger than those of previous models.

Herein, the musculoskeletal model incorporated an optimization algorithm based on balance and stability conditions to measure the axial force of muscles in the upright‐standing position and compared the overall concave/convex ratio of the axial force of muscles, including the concave–convex multifidus and longissimus muscles, with the musculoskeletal model of Kamal and Rouhi [51]. In our model, the ratios of the multifidus were 0.95 for the concave–convex muscle, 0.82 for the longissimus psoas, 0.89 for the longissimus pectoris, 0.75 for the quadratus psoas, and 0.86 for the psoas major. The overall ratio was 0.86, which was close to the concave/convex ratio of the paravertebral muscle group in the previous FEM (0.88). In summary, the new model is reliable and complies with the biomechanical simulations.

2.6. Boundary and Loading Conditions for Brace Correction FEM Simulation

The Chêneau brace provides space for tissue displacement, growth, and respiratory movement by applying corrective force on the convex side of the spine and releasing pressure on the concave side. According to the correction mechanism of the Chêneau brace, the model's pelvis and sacral canals were restricted to X, Y, and Z axis of freedom, and T1 was allowed to do lateral bending, forward bending, and upward extension. Previous studies [52, 53] have shown that the corrective thrust was applied to the right chest to correct the right thoracic curve, the corrective reaction force was applied to the left axilla and the left lumbar side, and an upward extension force was applied to the axillae (400 N on the left and 200 N on the right) to simulate the “bending effect” and the “antigravity effect” of the brace (Figure 6).

The FEM was used to simulate the corrective force coupling of the brace in three regions of the torso, including the left axilla, left lumbar, and right lateral thoracic region (Table 4). The correction pressure of the left lumbar region (CPLL) was scaled down gradually, and the effects of the downward adjustment of the CPLL were simulated based on the three‐point correction force principle of the brace. According to our previous clinical study [20], lowering the correction pressure for the compensatory lumbar curvature induces a better correction rate for the MT curve. Therefore, in the simulation, the CPLL was reduced by 25%, 50%, 75%, and 100%, and the biological characteristics of the correction according to the maximum threshold setting were observed in other areas. The internal pressure threshold for the trunk tolerance was 0.035 MPa [45] and the axillary brace tolerance was 0.085 MPa [54]. Using these pressure threshold values, four different correct pressures were set (Table 3), including correct finite model element (CFME)1, CFME2, CFME3, and CFME4. The effective pressure area was determined based on previous studies, patient conditions, and clinical experience [38, 55]. The area of the correction force varied in different body parts (Table 4).

2.7. Biomechanical Analysis of the Brace Correction Simulation

The relevant parameters for the mechanical response of the corrected model were measured. First, the Cobb angle, TK, LL, and apical vertebral translation (AVT) were measured radiologically. The AVT in the proximal thoracic (PT) and MT regions were recorded as AVT_PT and AVT_MT, respectively. The correction rate for the Cobb angle was calculated using the following formula:

Second, the height ratios of the concave and convex sides of the upper and lower intervertebral spaces at the apex of vertebrae T7 and T2 (for T1–T2, T2–T3, T6–T7, and T7–T8) were measured and calculated. Third, the von Mises stress ratios for the concavities and convexities of the T1–T2, T2–T3, T6–T7, and T7–T8 discs were determined. The average value for the IVD stress was obtained by taking 10 points in the convex and concave, respectively (Figure 7) [56]. A larger concave/convex ratio of the height of the intervertebral space means that the concave side of the intervertebral space increased to create favorable conditions for disc remodeling. A smaller stress ratio for the concave side of the disc indicates that the originally overstressed concave side of the disc is released, breaking the vicious cycle of scoliosis progression by creating favorable remodeling through the disc. The ratio of axial forces in the concave and convex side muscles was measured. Based on the premise that the spinal skeleton returns to the midline, balanced axial force on the concave and convex sides of the muscle is beneficial to the construction of reasonable tensegrity.

3.1. Radiological Parameter

The radiological parameters when different combinations of corrective forces were applied in the model are given in Table 5. In the coronal plane, the Cobb angle correction rates were highest in the MT curve and lumbar curve when the CPLL was 0.035 MPa (100% threshold), followed by 75%, 50%, and 25% thresholds. In contrast, the correction rates of the PT curve were 25%, 50%, 75%, and 100% in the order of high to low. The distance reduction sequences for AVT_T7 were 100%, 75%, 50%, and 25%. The distance reduction sequences for AVT_T2 were 25%, 50%, 75%, and 100%. The TK angle for the sagittal plane was 7.4° when the CPLL was 0.09 MPa, which was the maximum value at 25%. A TK angle of < 10° is diagnosed as a flatback deformity. The remaining correction force pressure combination of TK was in the range of 10°–40°. In the simulation, LL decreased when the correction threshold of the left lumbar region reached 50% and 75%, which was still in the normal physiological range.

3.2. Intervals Height Concave to Convex Ratios

Variations in the concave/convex height ratios in the upper and lower IVDs at the apical vertebra when simulating four combinations of corrective pressures are shown in Figure 8. For the concave–convex ratio of the intervertebral space height of the vertebrae at the apex of the PT curvature, the most satisfactory performance was achieved at 25% CPLL, followed by 50%, 70%, and 100%. In the MT region, the best performance was achieved when CPLL was 100%, followed by 75%, 50%, and 25%. Lumbar curve measurements were most satisfactory when the CPLL was 75%, followed by 100%, 50%, and 25%.

3.3. The von Mises Stress Ratios From Concave to Convex in IVDs

The concave/convex stress ratios of the disc at the apex are shown in Figure 9. The most satisfactory performance in the concave/convex ratio of the upper and lower discs at the apex of the PT region was at 25% CPLL, followed by 75%, 50%, and 100% CPLL. At the apex of the MT, the most satisfactory performance was at 100% CPLL, followed by 75%, 50%, and 25% CPLL. The most satisfactory performance of the concave/convex ratio of the upper and lower discs at the apex of the lumbar curve was at 75% CPLL, followed by 50%, 100%, and 25%.

3.4. Muscles Axial Force in Concave to Convex Ratios

Changes in the axial stress ratio of muscles on the concave and convex sides after different corrective pressures loading are shown in Figure 10. After the correction group was applied, the ICPL/IO/EO/LGPL tended to be balanced; when the CPLL boundaries were 25%, 50%, 75%, and 100%, the balance changed from best to worst. The other muscles were more unbalanced than before the correction in each correction combination, especially when the LGPT reached 100%. Conversely, imbalances occurred in the axial force ratio of muscles, including an increase in MF and a decrease in RA, QL, and PM.

4.1. Main Findings

In this study, the FEM of “tensegrity” of AIS was constructed and effectively verified as the first innovative result of this study. Another unprecedented result from the simulation of brace intervention was that different combinations of corrective forces showed different characteristics in terms of radiological measurements, disc stress distribution and height, and muscle axial tension. Specifically, the combination with the highest correction pressure resulted in a high Cobb correction rate for the MT curve, whereas the correction rate for PT and the disc stress at the apex of PT and MT exhibited unsatisfactory characteristics. Interestingly, the biomechanical performance of Rigo A classification can be achieved by properly downregulating the correction force of the lumbar lateral curve. It was a combination of radiometric measurements, disc, and muscle axial force parameters that gave us evidence to support our initial clinical observations and hypotheses that were consistent. These results suggest that when brace intervention is used to break the vicious cycle, it should not only consider the correction rate of Cobb Angle, but also focus on the asymmetry of IVD and muscle axial force.

4.2. Verification for Tensegrity of FEM

The geometry of the EFM in this study is consistent with the radiological measurements, and the ROM verification of the main curve is close to the previous model, supporting the validity of the model. The ROM of the T1–T4 spinal segments was within a reasonable range and close to the ROM reported by previous authors [48, 49, 50]. The simulated T5–T8 spinal segment ROM in the left and right lateral flexion directions exhibited a 6° difference compared with the autopsy results of Busscher et al. [48] and the differences in the rest of the directions were in the range of 1°–2°. The experimental results were also similar but slightly lower than the results of Xin et al. [50] and Sheng et al. [49]. The differences may be related to the individual morphologic differences in age, gender, and spinal and other skeletal morphologic characteristics that affect joint mobility in AIS patients.

4.3. Simulation of Chêneau Brace Intervention

4.4. Limitations and Strengths

To our knowledge, this study is the first to developmentally construct a musculoskeletal FEM that embodies the tensegrity of the spine as a brace intervention simulation for AIS. It is possible to obtain excellent biomechanical performance by reducing the correcting force of the compensation curve. Musculoskeletal models fitted with biomechanical characteristics may enrich the multidimensional thinking of clinical intervention with braces. There are several limitations to this study. First, the effects of immediate correction of the brace were modeled, and the effects of long‐term correction of the brace were not analyzed. Immediate correction may correlate with long‐term treatment. However, to obtain a more detailed long‐term treatment effect, future studies should include dynamic simulations of individual growth. Second, no complete soft tissue model such as muscle was available to compare and validate the model. This experiment did not construct tissues such as internal organs and did not fully model the effects of human internal organs on the axial forces of spinal bones, IVDs, and muscles. Finally, the effects of individual patient differences cannot be excluded. In the future, a larger number of patients with more classifications will be included to enhance the generality of the findings.

4.5. Prospect of Clinical Application

Our simulation study has implications for brace correction in patients with challenging multiple curves. Theoretically, the design of the brace would be different depending on the type of scoliosis. As more scoliosis models are developed, the biomechanical effects of different types of scoliosis braces should be analyzed to improve the predictability of clinical outcomes.