Fabio Galbusera1*
Tito Bassani1
Luigi La Barbera2
Claudia Ottardi2
Benedikt Schlager3
Marco Brayda-Bruno1
Tomaso Villa2Hans-Joachim Wilke3
- 1IRCCS Istituto Ortopedico Galeazzi, Milan, Italy
- 2LaBS, Department of Chemistry, Materials and Chemical Engineering, Politecnico di Milano, Milan, Italy
- 3Center of Musculoskeletal Research Ulm, Institute of Orthopedic Research and Biomechanics, Ulm University, Ulm, Germany
In decades of technical developments after the first surgical corrections of spinal deformities, the set of devices, techniques, and tools available to the surgeons has widened dramatically. Nevertheless, the rate of complications due to mechanical failure of the fixation or the instrumentation remains rather high. Indeed, basic and clinical research about the principles of deformity correction and the optimal surgical strategies (i.e., the choice of the fusion length, the most appropriate instrumentation, and the degree of tolerable correction) did not progress as much as the implantable devices and the surgical techniques. In this work, a software approach for the biomechanical simulation of the correction of patient-specific spinal deformities aimed to the identification of its biomechanical principles is presented. The method is based on three-dimensional reconstructions of the spinal anatomy obtained from biplanar radiographic images. A user-friendly graphical user interface allows for the planning of the desired deformity correction and to simulate the implantation of pedicle screws. Robust meshing of the instrumented spine is provided by using consolidated computational geometry and meshing libraries. Based on a finite element simulation, the program is able to predict the loads and stresses acting in the instrumentation as well as those in the biological tissues. A simple test case (reduction of a low-grade spondylolisthesis at L3–L4) was simulated as a proof of concept, and showed plausible results. Despite the numerous limitations of this approach which will be addressed in future implementations, the preliminary outcome is promising and encourages a wide effort toward its refinement.
Introduction
Spinal deformities are relatively frequent pathologies which may appear in different stages of the human life, from childhood and adolescence to maturity and old age. Deformities may occur in a single anatomical plane, e.g., hyperkyphosis (Zaina et al., 2009) and sagittal imbalance (Barrey et al., 2011), or may be fully three dimensional in nature, as in many cases of adolescent idiopathic scoliosis (Schlosser et al., 2014). Despite conservative treatment, such as physiotherapy or bracing, could be a valuable option to improve the quality of life and to stop the progression of the deformity (Zaina et al., 2009), corrective surgical treatment may be necessary in more severe cases.
In decades of technical developments after the first surgical corrections of spinal deformities (Harrington, 1962; Nachemson, 1971), the set of devices, techniques, and tools available to the surgeons has widened dramatically. Recent reports show cases in which severe scoliosis [e.g., Xie et al. (2012) and Cecchinato et al. (2015)], sagittal imbalance (Berjano and Aebi, 2015), or congenital deformities, such as hemivertebra (Zhuang et al., 2015), have been surgically treated with advanced techniques and success rates, which would have been inconceivable 10 or 20 years ago. Nevertheless, the rate of complications due to mechanical failure of the fixation or the instrumentation remains rather high, especially regarding the most severe cases [e.g., Bianco et al. (2014)]. Indeed, basic and clinical research about the principles of the deformity correction and the optimal surgical strategies (i.e., the choice of the fusion length, the most appropriate instrumentation, and the degree of tolerable correction) did not progress as much as the implantable devices and the surgical techniques.
There is general consensus that the key problem in this regard concerns the anatomical and biomechanical inter-patient variability of the spinal deformity (Aubin et al., 2007), which appears to be so significant that even well-established classification systems have a relatively low usefulness for the choice of the optimal surgical strategy for a specific patient. In a rather recent paper, Aubin et al. (2007) submitted radiographic and clinical data about five patients suffering from adolescent idiopathic scoliosis to six experienced spine surgeons who were asked to provide their preferred surgical planning, and obtained a large variability in the responses, which was attributed to the lack of clearly defined strategies or rational rules. Indeed, spinal deformity surgery is still often regarded as an “art” in which the principles cannot be formulated in a clear way, and the personal experience of the surgeon still plays a vital role.
A precious support for the identification of the optimal correction strategies may come in the future from biomechanical simulations. The research group lead by Carl-Éric Aubin at the Ecole Polytechnique de Montréal devoted a vast amount of resources for the development of a Spine Surgery Simulator (S3), which could simulate based on a biomechanical model the outcome of surgical correction of spinal deformities, thus supporting the surgeon in the identification of the optimal strategy for a specific patient with scientifically solid predictions (Aubin et al., 2008; Majdouline et al., 2009). The simulator was designed to achieve a high degree of user-friendliness and speed, so that it could be proficiently used in a clinical setting. Despite its limited availability to potentially interested surgeons, the research project has proved its technical feasibility and clinical relevance by means of a strict validation process against documented surgical results (Aubin et al., 2008).
In this work, an alternative software approach for the biomechanical simulation of surgical correction of patient-specific spinal deformities is presented. The method basically differs from S3 in its aims, since it was designed to be a tool for the identification of the biomechanical principles of correction rather than in the ability to provide fast bedside predictions to be used for a specific patient in a clinical setting. To these aims, the more general framework of the finite element method instead of multibody dynamics simulation was chosen in order to potentially provide better insight into detail aspects, such as local stresses and strains in bone, intervertebral disks, ligaments, and spinal instrumentation, at the cost of higher computational requirements.
Materials and Methods
Implementation and Software Tools
Similarly to a previous study (Vergari et al., 2015), the proposed method is based on three-dimensional (3D) reconstructions of the spinal anatomy carried out with sterEOS software on biplanar radiographic images, which is described and validated elsewhere (Somoskeoy et al., 2012; Ilharreborde et al., 2014). The method is implemented in C++ and makes use of the QT development platform1 to manage the graphical user interface. 3D rendering and user interactions are based on the libQGLViewer library.2 All basic operations on triangulated surfaces are carried out by the free GNU Triangulated Surface (GTS) library.3 A number of more sophisticated processing steps, specified in the following paragraphs, are performed with Meshlab software4 through its scripting interface. Tetrahedral meshing is carried out by using the free Tetgen library (Si, 2011). Finite element models are solved with the commercial ABAQUS package (Dassault Systèmes Simulia Corp., Johnston, RI, USA) and automatically postprocessed by means of its Python scripting interface.
Instrumented Spine Model
As mentioned above, the patient-specific finite element model of the instrumented spine is built based on the 3D reconstruction of biplanar radiographic images (Figure 1). Triangulated surfaces of all vertebrae and pelvis in their correct anatomical position and orientation are directly extracted from the DICOM file including the reconstruction generated by sterEOS. This piece of software uses algorithm-based parametric models of vertebrae and an inference method (Humbert et al., 2009) to generate the 3D reconstruction of the spine, which always starts from the same triangulated surfaces for each vertebra. Taking advantage of this choice, the localization of specific landmarks and areas can be easily performed by identifying specific nodes in the deformed triangulated surface, the numbering of which was determined in the original, undeformed surface by visual inspection (Figure 2). Since the node positions are morphed together with the vertebra, this method allows for an accurate and repeatable identification of landmark points and vertebral regions, such as endplates, insertions of ligaments and facet joints.
Figure 1. Biplanar radiographic images obtained with the EOS Imaging System of a patient suffering from adolescent idiopathic scoliosis (left). Patient-specific finite element mesh of the thoracolumbar spine built with the present method, and (right) detailed view of the L3–L4 spinal segment of the same patient.
Figure 2. Tracking of a specific node representing an anatomical landmark during the morphing procedure performed by sterEOS software, in which the surface mesh connectivity is not altered. The closest node is then identified after coarsening of the surface (right).
Pedicle screws are then positioned and oriented in the desired locations by using a graphical user interface (Figure 3). When the screw positioning is deemed optimal by the user, a new triangulated surface of the instrumented vertebra is generated by Boolean subtraction between the original vertebral surface and the surface of the screws. Tetrahedral finite element meshes of the vertebra and screws are then automatically generated. Prior to the generation of the 3D mesh, the user has the possibility to refine or coarsen the triangulated surface, in order to control the size of the final model.
Figure 3. Graphical user interface employed to properly position pedicle screws in each reconstructed vertebra.
When the volume meshes of two adjacent vertebrae have been generated, the intervertebral disk is automatically created. First, the surface nodes belonging to the endplate regions of the vertebrae are identified as described above. Then, the convex hull of the selected nodes is computed and resampled by means of the Poisson surface reconstruction algorithm (Hou et al., 2015) implemented in Meshlab. The elements belonging to a cylindrical region around the craniocaudal axis, the diameter of which is automatically calculated so that its volume is the 50% of that of the whole disk (Iatridis et al., 1996), are then identified as the nucleus pulposus (Figure 4), whereas the remaining elements constitute the annulus fibrosus. Annulus fibers are then created as non-linear springs arranged in order to have an average slope of ±30° with respect to the axial plane, by using a method described in detail elsewhere (Galbusera et al., 2006).
Figure 4. Finite element model of the intervertebral disk (left), showing the cylindrical nucleus pulposus and the collagen fibers in the annulus fibrosus. Assignment of the elements of the vertebra to distinct materials representing bone with different material properties (trabecular, cortical shell of the vertebral body and of the pedicles, posterior elements) (right).
Six ligament groups are modeled as non-linear springs: anterior longitudinal ligament, posterior longitudinal ligament, flaval ligament, interspinous ligament, supraspinous ligament, and capsular ligaments. Most ligaments were modeled by means of three spring elements, except for the capsular ligament which was modeled with 12 elements following the border of the facet capsule and the supraspinous ligament, which was represented by a single spring. The contact between the facet joints was modeled with GAPUNI elements (Vena et al., 2005), with a fixed initial clearance of 0.4 mm and no friction (Schmidt et al., 2012). Nodes representing the ligament insertions and facet joints were located by using the method described above and represented in Figure 2.
All bony tissues were modeled as linear isotropic continua (Figure 4). Similarly as done for endplates, ligaments, and facet joints, nodes on the vertebral surface belonging to distinct zones (vertebral bodies, pedicles and posterior elements) were identified. Appropriate materials were then associated to the volume elements of the vertebra based on a proximity criterion, i.e., each element was associated to the surface node closest to the element centroid, and thus to the corresponding material. For the vertebral bodies and pedicles, elements were classified as cortical if belonging to the outermost layer, or trabecular elsewhere. A single material was used to represent laminae and spinous, articular, and transverse processes.
In this first implementation, the same material properties taken from the literature were used in the whole model and for all patients (Table 1). A preliminary validation was conducted in order to assess the plausibility of the spinal flexibility obtained with a model based on EOS data of a subject without any spinal pathology. No attempts were done yet for a proper, more comprehensive validation of the modeling approach or for a patient-specific calibration of the material properties.
Table 1. Material properties used in the finite element models.
Planning of the Deformity Correction Surgery
After a finite element mesh of the spine instrumented with pedicle screws has been generated, the desired surgical correction of the deformity is planned manually with a dedicated graphical user interface (Figure 5A). The user has the possibility to translate and/or rotate a single vertebra or a group of vertebrae by using buttons and sliders, until a satisfying desired correction has been achieved. The task is facilitated by a multiwindow 3D viewer and the possibility to show or hide the original finite element mesh in the deformed configuration.
Figure 5. Graphical user interface employed to plan the deformity correction (A). The original, deformed spine is shown in dark gray, whereas the planned correction in green. Posterior rods (in blue) automatically designed in order to optimally fit the desired shape of the spine (B).
Posterior rods are created when the planning of the correction is concluded. The contour of the rods is automatically generated based on the corrected spinal configuration as a bicubic spline passing through the screw heads, which are assumed in this implementation to move rigidly with the screw shaft, thus modeling a monoaxial screw (Figure 5B). The splines are then meshed with linear beam elements with material and geometrical properties mimicking those of the posterior rods of interest.
The user has the possibility to select, among all the pedicle screws present in the model, which screws should be connected to the left and right posterior rods, on the two sides independently. Despite allowing for unrealistic strategies, i.e., in the surgical practice, all pedicle screws are connected to the rods, this choice allows for an easy and fast comparison of different correction strategies by using the same finite element mesh.
Furthermore, the user has the possibility to simulate discectomies at selected levels, as commonly done in scoliosis surgery to facilitate the correction of the curves (Waisman and Saute, 1997). In this case, the intervertebral disk is completely removed at the chosen levels, as well as the anterior and posterior longitudinal ligaments. No kinematic or contact interactions were defined between adjacent vertebral endplates.
The simulation of the correction maneuver is then performed, in three loading steps as described below (Figure 6): (i) in the first step, displacement boundary conditions are applied to the nodes on the screw heads to reach the corresponding positions on the posterior rods, thus implementing the correction. The spine is, therefore, strained, and reaction forces are acting on the screw heads. The rods are fixed in space and not loaded, since the engagement between screws and rods is not modeled in this first step. The lower extremities of the posterior rods are constrained in all degrees of freedom. (ii) Kinematic couplings are created between the screw heads and the corresponding nodes on the rods, by means of a Fortran user subroutine, MPC (ABAQUS 6.10 Documentation, 2010). Boundary conditions on the screw heads are deactivated. The strain energy of the spine is released to the rods, which become strained and loaded as well. Upper and lower extremities of the model are not constrained except for the relevant screw heads. (iii) In the third optional step, external loads (e.g., a compressive load or a flexion–extension moment) are applied as desired. Boundary conditions are not changed, and the kinematic couplings remain active.
Figure 6. Overview of the workflow to perform a biomechanical investigation of a spine deformity correction. First (“PLANNING”), the three-dimensional correction is planned as desired with the graphical user interface, and posterior rods with appropriate contours are generated. Then (“SIMULATION”), the finite element simulation is performed in three steps: (i) the correction maneuver; (ii) rods and pedicle screws are locked together; and (iii) external loads are applied, if desired. Finally (“DATA ANALYSIS”), relevant mechanical outputs are extracted from the simulation and analyzed.
Postprocessing and Data Analysis
When the finite element simulation is concluded, the ABAQUS output file is automatically processed by a Python script, which extracts the relevant results. First, nodal displacements are obtained and used to deform the mesh in the 3D viewer (Figure 5). Global orientations of the vertebrae before and after correction are geometrically calculated based on the original and deformed coordinates of landmark points. A direct comparison between the planned deformity correction and those actually achieved is, therefore, also possible.
Mechanical parameters describing the loads acting on the spinal instrumentation are also extracted. The forces acting between pedicle screws and vertebrae are calculated and plotted as vectors in the 3D viewer for an easy interpretation. Stresses and internal actions in the rods are also derived along the whole rod contours. In future implementations of the method, other relevant output values which may emerge as significant could be easily extracted via Python scripting.
Test Case
In order to test the plausibility of the results obtained with the approach, a simple case of spinal deformity has been selected as a benchmark. Biplanar radiographs of a patient suffering from grade I/II degenerative spondylolisthesis at L3–L4 (anterior slippage of 14 mm) have been collected. Apart from the listhesis, the patient had no other major spinal deformities. In the L2–L5 tract which was considered in the biomechanical model, minor coronal (6°) and axial (3°) asymmetries could be observed. Two possible surgical corrections were modeled: a partial reduction of the slippage of 7 mm and a complete reduction (Figure 7). In both cases, the correction maneuvers were purely translational, i.e., no rotational correction of the spinal shape was attempted. The spine was fixed in the corrected configuration by means of pedicle screws and posterior rods at L2–L5.
Figure 7. Patient-specific model used as a test case. In the pathological configuration (“ORIGINAL”), an anterior slippage of 14 mm is observable at L3–L4. Partial reduction to 7 mm is simulated by displacing L2 and L3 with the models “P–H,” “P–S,” and “P–D,” whereas complete correction is modeled with “C–H,” “C–S,” and “C–D.”
For each correction strategy, three cases were simulated: (i) no anterior surgery is performed, and the L3–L4 disk has the stiffness of that of a normal, healthy lumbar disk; (ii) no anterior surgery is performed, and the L3–L4 disk has half the stiffness of a normal disk (Table 1), thus modeling instability; and (iii) a complete discectomy is performed at L3–L4. Thus, six conditions were simulated: (i) partial reduction of the listhesis and healthy disk (“P–H”); (ii) partial reduction and “soft” disk (“P–S”); (iii) partial reduction and discectomy (“P–D”); (iv) complete reduction and healthy disk (“C–H”); (v) complete reduction and “soft” disk (“C–S”); and (vi) complete reduction and discectomy (“C–D”). For all models, the achieved correction was quantitatively compared to the desired correction in terms of vertebral orientations in the three anatomical planes. Furthermore, forces acting at the screw–bone interface as well as internal actions in the posterior rods were calculated.
Results
Usability of the Program
To test the robustness and usability of the newly developed computer program, biplanar radiographic images of eight patients suffering from mild adolescent scoliosis were processed and finite element meshes of the whole thoracolumbar spines were generated (Figure 8). After 3D reconstruction of the images with sterEOS software, the computer program allowed for a robust meshing of all spinal anatomies, requiring ~15–30 min of work by a skilled operator for each case. The planning of the desired correction could also be consistently performed in 5–10 min. The ABAQUS simulation of the correction maneuver required a variable time dependent on the number of nodes in the mesh, and ranged between 15 min for short spinal segments with a coarse mesh to several hours for the whole thoracolumbar spines. As expectable, convergence issues emerged in some cases of extreme corrections involving high local strains, especially in the intervertebral disks. Automated postprocessing of the ABAQUS solution by means of Python scripting required <1 min.
Figure 8. Finite element meshes of the thoracolumbar spines of eight patients suffering from adolescent idiopathic scoliosis, ranging from mild to severe grades.
Test Case – Correction Potential of the Surgical Treatment
All the various surgical strategies allowed for achieving a correction very similar to the planned one, with minor differences (Figure 9). In the sagittal plane, the maximal discrepancies were found for the orientation of L3, which ranged between 0° and 3° for the various configurations. As expectable, higher differences with respect to the planned alignment were found for the complete reduction scenarios, if compared to the partial reduction which required the application of a lower correction force. This effect was particularly pronounced for the C–H model, due to the higher stresses acting in the stiffer anterior structures.
Figure 9. Vertebral rotations predicted for the test cases in the sagittal (left), coronal (center) and axial (right) planes (in each panel around the axes highlighted in red) in the six configurations “C–H,” “C–S,” “C–D,” “P–H,” “P–S,” and “P–D”. “PLANNED” represents the vertebral orientations of the desired correction.
Similar trends were observed for rotations in the coronal and axial planes (Figure 9), with maximal deviations found in both cases at L3. Nevertheless, the magnitudes of these differences were in all cases below 0.5° and could be, therefore, considered negligible.
Loads on the Screws
Within the six configurations considered, the forces acting between each pedicle screw and the relevant vertebrae generally had the same directions and significantly differed only in magnitude (Figures 10 and 11). Coherently with the intuition, pull-out forces were found on the screws implanted in L3, whereas push-in forces with lower (on the left side of the instrumentation) or comparable magnitudes (on the right side) were predicted at the L4 level. Due to the proximity of the boundary conditions, negligible forces were calculated at L5, whereas the bone–screw interfaces at L2 were mainly loaded in the push-in, caudal direction.
Figure 10. Forces acting at the pedicle screw–bone interfaces, represented by the red arrows (sizes proportional to the magnitude of the forces), calculated for the six configurations. Pull-out forces are predicted at L2 (arrows directed outwards), whereas push-in forces at L3 (arrows directed inward).
Figure 11. Total forces (left) and pull-out forces (right) at the screw–bone interfaces for the screws on the left side (top) and those on the right side (bottom) predicted for the six configurations of the test case.
The highest forces were predicted with the C–H, whereas C–S and P–H provided very similar values. The magnitude of the highest calculated forces (>300 N) may suggest a risk of screw pull-out in L3 for the C–H model (Zdeblick et al., 1993; Pfeiffer et al., 1996), whereas cut-out could be safely excluded due to the direction of the force being mostly aligned with the screw shaft. Complete discectomy at L3–L4 (C–D and P–D) allowed for a major reduction of the loads on the screws. Minor differences were found between the left and right side of the instrumentation; these deviations could be attributed to the slight asymmetries in the model, especially regarding the coronal orientation of L4 and the axial orientation of L2.
Loads in the Rods
Internal effects in the posterior rods also showed similar patterns at the left and right sides of the instrumentation, despite with some non-negligible differences, which may be attributed to the asymmetries in the spinal shape (Figures 12 and 13). Similarly on the two sides, anteroposterior shear forces had a significantly higher magnitude with respect to axial and laterolateral shear forces at the level of the spondylolisthesis reduction (L3–L4). On the contrary, axial forces where dominant at L2–L3, whereas laterolateral forces were low or negligible at all levels. Among the six configurations, patterns similar to those predicted for the screw forces emerged, with maximal load values for C–H, followed by C–S, P–H, and P–S, whereas negligible rod loads were calculated for C–D and P–D.
Figure 12. Internal effects acting in the left rod of the test case. First row, from left to right: axial force, anteroposterior shear force, laterolateral shear force. Second row: moment in the sagittal plane (flexion–extension), in the coronal plane (lateral bending), in the axial plane (torsion). Sign conventions are indicated by the right arrows.
Figure 13. Internal effects acting in the right rod of the test case. First row, from left to right: axial force, anteroposterior shear force, laterolateral shear force. Second row: moment in the sagittal plane (flexion–extension), in the coronal plane (lateral bending), in the axial plane (torsion). Sign conventions are indicated by the right arrows.
Bending and twisting moments in the rods showed higher differences between the two sides (Figures 12 and 13). Flexion–extension bending moments reached significant values at L4–L5, especially on the left side (+53% with respect to the right side), but were marginal at the level of the reduction as well as lower, and with opposite sign, at L2–L3. This finding is coherent with the intuitive expectations, for which the rods should be loaded mostly in anteroposterior shear at L3–L4 and in flexion–extension at the adjacent levels. Non-negligible values of the torsion moment were predicted at L4–L5, especially on the left side, and at L3–L4 on the left side only.
In summary, all the six configurations allow for a satisfactory correction of the spondylolisthesis, but differ significantly regarding the loads on the instrumentation and at the bone–screw interface. The only critical scenario appears to be C–H, which may induce a non-negligible risk of screw pull-out. In general, internal effects on the posterior instrumentation have magnitudes lower or comparable to those observed in physiological loading (Rohlmann et al., 1995), and therefore the rods should not be in danger of mechanical failure.
Discussion
In this paper, a computational approach to simulate spinal deformity correction based on biplanar radiographic images of specific patients is presented. A user-friendly graphical interface allows the user to plan the desired correction strategy and to easily predict the achievable correction as well as mechanical variables, such as the stresses in the instrumentation and in the biological tissues. At this stage, the method is designed in order to be able to identify the biomechanical principles of deformity correction rather than as a bedside tool to be used to determine the optimal surgical strategy for a specific patient. In this respect, it differs from the previously mentioned S3 simulator, which is able to provide comparisons of various preoperative plannings including the various surgical maneuvers in a short time, being based on a relatively simple mechanical framework, such as multibody modeling. Nevertheless, the method also differentiates itself from other image-based programs, such as Surgimap (Akbar et al., 2013), which do not attempt to simulate the mechanics of the deformity correction but are limited to provide smart tools for measurements and comparisons.
As a proof of concept, we selected a simple exemplary case of low-grade spondylolisthesis which allows for an intuitive analysis of the results, without attempting either to provide a detailed analysis of a complex case or to perform an extensive investigation of spondylolisthesis reduction. Nevertheless, the predictions were plausible and respected the general principles of deformity correction. Since there were only minor coronal and axial deformities, screw loads, and internal effects on the rods were nearly symmetrical with respect to the sagittal plane. Pull-out screw forces acted on the screws implanted in the vertebra which was displaced in the posterior direction, whereas push-in forces were predicted in the lower vertebrae. At the level of spondylolisthesis reduction, the major internal effect in the rods was the anteroposterior shear force, while in the adjacent segments flexion–extension moments were dominant. In addition to this general, qualitative understanding of the loads acting in the level subjected to correction and in the adjacent segments, the model provided a quantification of the effects, which may be useful to estimate the risk of hardware failure and loosening, which is currently left to the experience of the surgeon. It should be noted that the finite element model could potentially be used for a deeper analysis of the mechanical response, such as the quantification of local stresses (Figure 14) as well as sensitivity analyses, within the limitations of the modeling approach described below.
Figure 14. Exemplary representation of the von Mises stresses in the L4 vertebra considered in the test case. External (A) and section (B) views highlighting the higher stresses in the pedicles. The section plane used to create the section view is shown in red.
Only a few published papers proposed methods for the estimation of the instrumentation loads during and after deformity correction. In the studies conducted with the S3 simulator, topics, such as the optimum screw patterns to minimize the forces required to achieve the desired corrections in adolescent scoliotic subjects (Wang et al., 2012b), the efficacy of some correction maneuvers (Wang et al., 2011; Martino et al., 2013) and of different instrumentations (Wang et al., 2012a) were investigated. Abe and coworkers (Abe et al., 2015) predicted the corrective forces in 20 adolescent patients based on the changes in rod geometry and finite element analysis, based on post-operative CT scans. More frequently, the biomechanics of spinal deformity has been investigated by means of finite element models of representative cases (e.g., Rohlmann et al. (2006) and Agarwal et al. (2015)], thus not considering the strong intersubject variability typical of most spinal deformities.
Convergence issues emerged in some cases involving highly strained areas, especially in the intervertebral disks, and were arguably related to the local mesh quality. However, the approach used for the geometrical generation of the intervertebral disks provides a smooth, high quality surface mesh, which constitutes an appropriate input for volume meshing (Hou et al., 2015). Nevertheless, local mesh refinement may improve convergence in selected cases, and should be taken into account as a possible future development.
Despite the potential of our approach in exceeding the limitations of the published works, many aspects relevant to the biomechanics of spinal deformity correction were not covered in the present implementation. First and foremost, only monoaxial pedicle screws were modeled. Polyaxial screws as well as hooks are fundamental tools which are used by most surgeons, even in hybrid constructs, and therefore support for them is needed. As a matter of fact, the revolutionary Cotrel–Dubousset approach to treat idiopathic scoliosis, which is still the base for many current instrumentation systems, included both polyaxial screws and hooks (Cotrel and Dubousset, 1984). Other surgical techniques which are not currently supported are spinal osteotomies, such as Ponte and pedicle subtraction osteotomies, which are widely used for the treatment of adult sagittal imbalance (Gill et al., 2008), and other specific instrumentations which may be used in some cases, such as interbody cages, crosslinks, double-rod systems with dominos and sublaminar wires.
Other limitations pertain to the modeling approach itself. Despite comprehensive biomechanical data are not available yet and preliminary studies appeared only very recently (Mannen et al., 2015), it is known that the rib cage has a considerable stiffness and may further reduce the flexibility of the trunk (Oda et al., 2002), with considerable effects on the forces necessary to achieve the correction and the loads acting on the instrumentation. Rib humps which may be present in severe scoliotic cases may also have an influence on the correction procedure (Clin et al., 2006), but the extent and relevance of this effect is nowadays unknown and cannot, therefore, be simulated.
No boundary conditions or other mechanical constraints were imposed to the upper and lower extremities of the spine model, which were, thus, free to move without any resistance. This assumption collides with the real corrections in which other anatomical structures, such as head and pelvis, are connected to the spine and limit its ability to be deformed during the correction maneuvers. Another simplification concerns the connection between screws and rods by means of the MPC user subroutine neglecting significant information about the shape of the instrumentation, and thus preventing to use rod stresses obtained in proximity of the screws to evaluate the risk of hardware failure in this region due to possible numerical artifacts and low accuracy.
Figure 14 highlights the approximate nature of the proximity criterion used to assign the material properties to the bone tissue. As a matter of fact, this approach does not allow for the precise definition of a cortical shell with a predefined thickness, and may cause irregular stress distributions. An approach including the explicit definition of separate cortical and trabecular volumes is currently being developed.
The material properties of the biological tissues were taken from the literature and were mostly derived from experiments conducted on healthy specimens. As a matter of fact, patients with deformities may exhibit a stiffening of the spinal soft tissues, to an extent which is strongly variable from patient to patient and from level to level (Lafon et al., 2010). At the same time, many elderly patients suffer from osteopenia or osteoporosis, which were shown to have a non-negligible prevalence also in adolescent scoliotic subjects (Ishida et al., 2015) and significantly reduce the stiffness and strength of the bony tissues. Therefore, it appears that an effort toward the incorporation of patient-specific material data as they could be derived from biomedical imaging (Lafon et al., 2010) is required in order to improve the accuracy and reliability of the numerical predictions, and is planned as a future development. Another possible strategy to integrate deformity-dependent material properties is given by intra-operative measurements of the segmental stiffness by means of instrumented forceps (Klockner et al., 2003; Reutlinger et al., 2012) or stepper motor-based equipment (Brown et al., 2002), which may be correlated with the radiological appearance of the spine and thus exploited to implement localized material properties. For the sake of simplicity, bone was modeled as a simple linear elastic continuum; a more sophisticated formulation including a failure criterion would allow for more accurate predictions of the risk of loosening of the instrumentation as well as of vertebral fracture.
The main limitation of the present work is, however, the lack of model validation. As a matter of fact, patient-specific biomechanical models of deformed spines are nearly impossible to validate, due to their inherent variability and the lack of data apart from radiological imaging or clinical assessment. Experimental testing with cadaveric spines has been extensively used for the validation of numerical models, but is largely unpractical for the investigation of spinal deformities. As a matter of fact, deformed spine specimens are very rare and exhibit a strong inter-specimen variability, which prevents a proper repeatability and statistical analysis of the results. Alternative methods [e.g., deforming physiological specimens, such as in Wilke et al. (2015)] are currently under development, but have significant limitations and cannot be considered as consolidated methods.
A proper verification of the finite element models is, however, possible and necessary before the predictions could be used to extract clinically relevant information (Viceconti et al., 2005). In particular, mesh sensitivity studies are required in order to ensure that the predictions are mesh independent. It should be noted that the method allows for an arbitrary refining or coarsening of the meshes, by using algorithms available in the Meshlab and TetGen computer programs.
In summary, taking into account our aim of determining the biomechanical principles of deformity correction rather than accurately modeling a surgical procedure for a specific patient, we believe the preliminary outcome of our approach is very promising, and encourages a wide effort toward its refinement.
Conflict of Interest Statement
The research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Funding
Funding by the Italian Ministry of Health (Project GR-2011-02351464) is gratefully acknowledged.
Footnotes
References
ABAQUS 6.10 Documentation. (2010). User Subroutines Reference Manual, Section 1.1.13, “User Subroutine to Define Multi-Point Constraints”. Johnston, RI: Sassault Systemes Simulia Corp.
Abe, Y., Ito, M., Abumi, K., Sudo, H., Salmingo, R., and Tadano, S. (2015). Scoliosis corrective force estimation from the implanted rod deformation using 3D-FEM analysis. Scoliosis 10, S2. doi:10.1186/1748-7161-10-S2-S2
Agarwal, A., Zakeri, A., Agarwal, A. K., Jayaswal, A., and Goel, V. K. (2015). Distraction magnitude and frequency affects the outcome in juvenile idiopathic patients with growth rods: finite element study using a representative scoliotic spine model. Spine J. 15, 1848–1855. doi:10.1016/j.spinee.2015.04.003
Akbar, M., Terran, J., Ames, C. P., Lafage, V., and Schwab, F. (2013). Use of Surgimap Spine in sagittal plane analysis, osteotomy planning, and correction calculation. Neurosurg. Clin. N. Am. 24, 163–172. doi:10.1016/j.nec.2012.12.007
Aubin, C. E., Labelle, H., Chevrefils, C., Desroches, G., Clin, J., and Eng, A. B. (2008). Preoperative planning simulator for spinal deformity surgeries. Spine (Phila Pa 1976) 33, 2143–2152. doi:10.1097/BRS.0b013e31817bd89f
Aubin, C. E., Labelle, H., and Ciolofan, O. C. (2007). Variability of spinal instrumentation configurations in adolescent idiopathic scoliosis. Eur. Spine J. 16, 57–64. doi:10.1007/s00586-006-0063-6
Barrey, C., Roussouly, P., Perrin, G., and Le Huec, J. C. (2011). Sagittal balance disorders in severe degenerative spine. Can we identify the compensatory mechanisms? Eur. Spine J. 20(Suppl. 5), 626–633. doi:10.1007/s00586-011-1930-3
Berjano, P., and Aebi, M. (2015). Pedicle subtraction osteotomies (PSO) in the lumbar spine for sagittal deformities. Eur. Spine J. 24(Suppl. 1), S49–S57. doi:10.1007/s00586-014-3618-y
Bianco, K., Norton, R., Schwab, F., Smith, J. S., Klineberg, E., Obeid, I., et al. (2014). Complications and intercenter variability of three-column osteotomies for spinal deformity surgery: a retrospective review of 423 patients. Neurosurg. Focus 36, E18. doi:10.3171/2014.2.FOCUS1422
Brown, M. D., Holmes, D. C., Heiner, A. D., and Wehman, K. F. (2002). Intraoperative measurement of lumbar spine motion segment stiffness. Spine (Phila Pa 1976) 27, 954–958. doi:10.1097/00007632-200205010-00006
Cecchinato, R., Berjano, P., Aguirre, M. F., and Lamartina, C. (2015). Asymmetrical pedicle subtraction osteotomy in the lumbar spine in combined coronal and sagittal imbalance. Eur. Spine J. 24(Suppl. 1), S66–S71. doi:10.1007/s00586-014-3669-0
Clin, J., Aubin, C. E., Parent, S., Ronsky, J., and Labelle, H. (2006). Biomechanical modeling of brace design. Stud. Health Technol. Inform. 123, 255–260.
Cotrel, Y., and Dubousset, J. (1984). A new technic for segmental spinal osteosynthesis using the posterior approach. Rev. Chir. Orthop. Reparatrice Appar. Mot. 70, 489–494.
Galbusera, F., Bellini, C. M., Anasetti, F., Ciavarro, C., Lovi, A., and Brayda-Bruno, M. (2011). Rigid and flexible spinal stabilization devices: a biomechanical comparison. Med. Eng. Phys. 33, 490–496. doi:10.1016/j.medengphy.2010.11.018
Galbusera, F., Fantigrossi, A., Raimondi, M. T., Sassi, M., Fornari, M., and Assietti, R. (2006). Biomechanics of the C5-C6 spinal unit before and after placement of a disc prosthesis. Biomech. Model. Mechanobiol. 5, 253–261. doi:10.1007/s10237-006-0015-4
Gill, J. B., Levin, A., Burd, T., and Longley, M. (2008). Corrective osteotomies in spine surgery. J. Bone Joint Surg. Am. 90, 2509–2520. doi:10.2106/JBJS.H.00081
Harrington, P. R. (1962). Treatment of scoliosis. Correction and internal fixation by spine instrumentation. J. Bone Joint Surg. Am. 44-A, 591–610.
Hou, W., Xu, Z., Qin, N., Xiong, D., and Ding, M. (2015). Surface reconstruction through poisson disk sampling. PLoS ONE 10:e0120151. doi:10.1371/journal.pone.0120151
Humbert, L., De Guise, J. A., Aubert, B., Godbout, B., and Skalli, W. (2009). 3D reconstruction of the spine from biplanar X-rays using parametric models based on transversal and longitudinal inferences. Med. Eng. Phys. 31, 681–687.
Iatridis, J. C., Weidenbaum, M., Setton, L. A., and Mow, V. C. (1996). Is the nucleus pulposus a solid or a fluid? Mechanical behaviors of the nucleus pulposus of the human intervertebral disc. Spine (Phila Pa 1976) 21, 1174–1184. doi:10.1097/00007632-199605150-00009
Ilharreborde, B., Dubousset, J., and Le Huec, J. C. (2014). Use of EOS imaging for the assessment of scoliosis deformities: application to postoperative 3D quantitative analysis of the trunk. Eur. Spine J. 23(Suppl. 4), S397–S405. doi:10.1007/s00586-014-3334-7
Ishida, K., Aota, Y., Mitsugi, N., Kono, M., Higashi, T., Kawai, T., et al. (2015). Relationship between bone density and bone metabolism in adolescent idiopathic scoliosis. Scoliosis 10, 19. doi:10.1186/s13013-015-0043-x
Klockner, C., Rohlmann, A., and Bergmann, G. (2003). Instrumented forceps for measuring tensile forces in the rod of the VDS implant during correction of scoliosis. Biomed. Tech. (Berl) 48, 362–364. doi:10.1515/bmte.2003.48.12.362
Lafon, Y., Lafage, V., Steib, J. P., Dubousset, J., and Skalli, W. (2010). In vivo distribution of spinal intervertebral stiffness based on clinical flexibility tests. Spine (Phila Pa 1976) 35, 186–193. doi:10.1097/BRS.0b013e3181b664b1
Lu, Y. M., Hutton, W. C., and Gharpuray, V. M. (1996). Can variations in intervertebral disc height affect the mechanical function of the disc? Spine (Phila Pa 1976) 21, 2208–2216; discussion 2217. doi:10.1097/00007632-199610010-00006
Majdouline, Y., Aubin, C. E., Sangole, A., and Labelle, H. (2009). Computer simulation for the optimization of instrumentation strategies in adolescent idiopathic scoliosis. Med. Biol. Eng. Comput. 47, 1143–1154. doi:10.1007/s11517-009-0509-1
Mannen, E. M., Anderson, J. T., Arnold, P. M., and Friis, E. A. (2015). Mechanical analysis of the human cadaveric thoracic spine with intact rib cage. J. Biomech. 48, 2060–2066. doi:10.1016/j.jbiomech.2015.03.021
Martino, J., Aubin, C. E., Labelle, H., Wang, X., and Parent, S. (2013). Biomechanical analysis of vertebral derotation techniques for the surgical correction of thoracic scoliosis. A numerical study through case simulations and a sensitivity analysis. Spine (Phila Pa 1976) 38, E73–E83. doi:10.1097/BRS.0b013e31827a641e
Oda, I., Abumi, K., Cunningham, B. W., Kaneda, K., and McAfee, P. C. (2002). An in vitro human cadaveric study investigating the biomechanical properties of the thoracic spine. Spine (Phila Pa 1976) 27, E64–E70. doi:10.1097/00007632-200202010-00007
Pfeiffer, M., Gilbertson, L. G., Goel, V. K., Griss, P., Keller, J. C., Ryken, T. C., et al. (1996). Effect of specimen fixation method on pullout tests of pedicle screws. Spine (Phila Pa 1976) 21, 1037–1044. doi:10.1097/00007632-199605010-00009
Pitzen, T., Geisler, F. H., Matthis, D., Muller-Storz, H., Pedersen, K., and Steudel, W. I. (2001). The influence of cancellous bone density on load sharing in human lumbar spine: a comparison between an intact and a surgically altered motion segment. Eur. Spine J. 10, 23–29. doi:10.1007/s005860000223
Reutlinger, C., Hasler, C., Scheffler, K., and Buchler, P. (2012). Intraoperative determination of the load-displacement behavior of scoliotic spinal motion segments: preliminary clinical results. Eur. Spine J. 21(Suppl. 6), S860–S867. doi:10.1007/s00586-012-2164-8
Rohlmann, A., Bergmann, G., Graichen, F., and Mayer, H. M. (1995). Telemeterized load measurement using instrumented spinal internal fixators in a patient with degenerative instability. Spine (Phila Pa 1976) 20, 2683–2689. doi:10.1097/00007632-199512150-00009
Rohlmann, A., Richter, M., Zander, T., Klockner, C., and Bergmann, G. (2006). Effect of different surgical strategies on screw forces after correction of scoliosis with a VDS implant. Eur. Spine J. 15, 457–464. doi:10.1007/s00586-005-0923-5
Schlosser, T. P., van Stralen, M., Brink, R. C., Chu, W. C., Lam, T. P., Vincken, K. L., et al. (2014). Three-dimensional characterization of torsion and asymmetry of the intervertebral discs versus vertebral bodies in adolescent idiopathic scoliosis. Spine (Phila Pa 1976) 39, E1159–E1166. doi:10.1097/BRS.0000000000000467
Schmidt, H., Galbusera, F., Rohlmann, A., Zander, T., and Wilke, H. J. (2012). Effect of multilevel lumbar disc arthroplasty on spine kinematics and facet joint loads in flexion and extension: a finite element analysis. Eur. Spine J. 21(Suppl. 5), S663–S674. doi:10.1007/s00586-010-1382-1
Si, H. (2011). TetGen: A Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator. Available at: http://wias-berlin.de/software/tetgen/
Somoskeoy, S., Tunyogi-Csapo, M., Bogyo, C., and Illes, T. (2012). Accuracy and reliability of coronal and sagittal spinal curvature data based on patient-specific three-dimensional models created by the EOS 2D/3D imaging system. Spine J. 12, 1052–1059. doi:10.1016/j.spinee.2012.10.002
Vena, P., Franzoso, G., Gastaldi, D., Contro, R., and Dallolio, V. (2005). A finite element model of the L4-L5 spinal motion segment: biomechanical compatibility of an interspinous device. Comput. Methods Biomech. Biomed. Engin. 8, 7–16. doi:10.1080/10255840500062914
Vergari, C., Ribes, G., Aubert, B., Adam, C., Miladi, L., Illharreborde, B., et al. (2015). Evaluation of a patient-specific finite-element model to simulate conservative treatment in adolescent scoliosis. Spine Deform. 3, 4–11. doi:10.1016/j.jspd.2014.06.014
Viceconti, M., Olsen, S., Nolte, L. P., and Burton, K. (2005). Extracting clinically relevant data from finite element simulations. Clin. Biomech. (Bristol, Avon) 20, 451–454. doi:10.1016/j.clinbiomech.2005.01.010
Waisman, M., and Saute, M. (1997). Thoracoscopic spine release before posterior instrumentation in scoliosis. Clin. Orthop. Relat. Res. 336, 130–136. doi:10.1097/00003086-199703000-00018
Wang, X., Aubin, C. E., Crandall, D., and Labelle, H. (2011). Biomechanical modeling and analysis of a direct incremental segmental translation system for the instrumentation of scoliotic deformities. Clin. Biomech. (Bristol, Avon) 26, 548–555. doi:10.1016/j.clinbiomech.2011.01.011
Wang, X., Aubin, C. E., Crandall, D., Parent, S., and Labelle, H. (2012a). Biomechanical analysis of 4 types of pedicle screws for scoliotic spine instrumentation. Spine (Phila Pa 1976) 37, E823–E835. doi:10.1097/BRS.0b013e31824b7154
Wang, X., Aubin, C. E., Labelle, H., Parent, S., and Crandall, D. (2012b). Biomechanical analysis of corrective forces in spinal instrumentation for scoliosis treatment. Spine (Phila Pa 1976) 37, E1479–E1487. doi:10.1097/BRS.0b013e3182706745
Wilke, H.-J., Mathes, B., Midderhoff, S., and Graf, N. (2015). Development of a scoliotic spine model for biomechanical in vitro studies. Clin. Biomech. 30, 182–187.
Xie, J., Wang, Y., Zhao, Z., Zhang, Y., Si, Y., Li, T., et al. (2012). Posterior vertebral column resection for correction of rigid spinal deformity curves greater than 100 degrees. J. Neurosurg. Spine 17, 540–551. doi:10.3171/2012.9.SPINE111026
Zaina, F., Atanasio, S., Ferraro, C., Fusco, C., Negrini, A., Romano, M., et al. (2009). Review of rehabilitation and orthopedic conservative approach to sagittal plane diseases during growth: hyperkyphosis, junctional kyphosis, and Scheuermann disease. Eur. J. Phys. Rehabil. Med. 45, 595–603.
Zdeblick, T. A., Kunz, D. N., Cooke, M. E., and McCabe, R. (1993). Pedicle screw pullout strength. Correlation with insertional torque. Spine (Phila Pa 1976) 18, 1673–1676. doi:10.1097/00007632-199309000-00016
Keywords: scoliosis, spinal deformity, finite element, deformity correction, spine biomechanics, patient specific
Citation: Galbusera F, Bassani T, La Barbera L, Ottardi C, Schlager B, Brayda-Bruno M, Villa T and Wilke H-J (2015) Planning the surgical correction of spinal deformities: toward the identification of the biomechanical principles by means of numerical simulation. Front. Bioeng. Biotechnol. 3:178. doi: 10.3389/fbioe.2015.00178
Received: 17 July 2015; Accepted: 19 October 2015;
Published: 03 November 2015
Edited by:
Lei Ren, University of Manchester, UKCopyright: © 2015 Galbusera, Bassani, La Barbera, Ottardi, Schlager, Brayda-Bruno, Villa and Wilke. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fabio Galbusera, ZmFiaW8uZ2FsYnVzZXJhJiN4MDAwNDA7Z3J1cHBvc2FuZG9uYXRvLml0