Computational model predicts risk of spinal screw loosening in patients

Pedicle screw loosening is a frequent complication in lumbar spine fixation, most commonly among patients with poor bone quality. Determining patients at high risk for insufficient implant stability would allow clinicians to adapt the treatment accordingly. The aim of this study was to develop a computational model for quantitative and reliable assessment of the risk of screw loosening. A cohort of patient vertebrae with diagnosed screw loosening was juxtaposed to a control group with stable fusion. Imaging data from the two cohorts were used to generate patient-specific biomechanical models of lumbar instrumented vertebral bodies. Single-level finite element models loading the screw in axial or caudo-cranial direction were generated. Further, multi-level models incorporating individualized joint loading were created. The simulation results indicate that there is no association between screw pull-out strength and the manifestation of implant loosening (p = 0.8). For patient models incorporating multiple instrumented vertebrae, CT-values and stress in the bone were significantly different between loose screws and non-loose screws (p = 0.017 and p = 0.029, for CT-values and stress, respectively). However, very high distinction (p = 0.001) and predictability (R2Pseudo = 0.358, AUC = 0.85) were achieved when considering the relationship between local bone strength and the predicted stress (loading factor). Screws surrounded by bone with a loading factor higher than 25% were likely to be loose, while the chances of screw loosening were close to 0 with a loading factor below 15%. The use of a biomechanics-based score for risk assessment of implant fixation failure might represent a paradigm shift in addressing screw loosening after spondylodesis surgery.


Introduction
Posterior spinal fixation of the lower spine is a gold standard surgical intervention and is employed for spine stabilization through the fusion of spinal segments. Pedicle screw loosening after spondylodesis surgery is a common complication that can lead to persistent pain, failed fusion, and the need for revision surgery [20,23]. Screw loosening has been reported to occur in about 10% of patients, and this percentage increases to above 60% for those who are affected by osteoporosis [5,10,21,28].

3
The ability to predict the risk of screw loosening through the examination of subject-specific biomechanical aspects could be groundbreaking for the improvement of surgical outcomes. It would allow to individualize the fusion strategy and to make treatment adaptations based on the predicted outcome of the intervention. For instance, patient care might be improved by using expandable or cement-augmented pedicle screws to increase holding power [6,12,13]. Also, the number of levels spanned by the fusion construct and the screw trajectory could be adapted [8]. Finally, anticipating long-term surgical outcomes may allow targeting at-risk patients with specific recommendations on cautious postsurgical behavior.
So far, there is no quantitative and reliable measure to identify level-specific risk of loosening besides a general evaluation of bone quality from image data [5,28]. However, such methods do not consider factors like patient-specific loading, which could be associated with screw fixation strength. Finite element (FE) and musculoskeletal modeling are well established methods in biomechanical research and can be complementary to each other [33]. However, until now, the clinical relevance and impact of such models have been hampered by the circumstance that the gained insight could usually not be directly translated into improved patient care.
We hypothesize that computational models allow for the prediction of patient-specific risk for screw loosening. To test this hypothesis, biomechanical models of varying complexity with respect to loading (simple uniaxial load or spinal loading derived from individualized musculoskeletal models) and spine representation (one or multiple vertebrae) were generated. These models were based on preoperative and post-operative imaging data of patients with and without diagnosed screw loosening.

Dataset
A retrospective database analysis was conducted to identify eligible patients who met the following inclusion criteria: • spondylodesis using pedicle screw instrumentation of the lumbar spine (including extension surgeries) between April 2014 and January 2019; • minimal follow-up time of half a year [32]; • no known bone disease (e.g. tumors of the spine); • available preoperative computed tomography (CT) scans and postoperative CT scans of the newly instrumented area of the spine (maximal slice thickness: 3 mm).
Further, vertebrae containing screws augmented with bone cement were not included and for multi-vertebrae models incorporating results from musculoskeletal analysis only patients with bi-planar x-ray images could be considered. For the case-cohort, screw loosening status was diagnosed by an experienced surgeon based on post-operative CT images and/or during revision surgery. A control group composed of vertebrae with non-diagnosed screw loosening in patients with comparable age and weight distribution was juxtaposed to the case group (Fig. 1).
In total, 26 patients (9 males, 17 females) were considered for single vertebra and 16 patients (6 males, 10 females) for multi-vertebrae modeling, respectively. Description of the patients and the included instrumentation can be found in Table 1.

Patient-specific modeling
The majority of the required steps for model generation were automatized with custom scripts in MATLAB (The MathWorks Inc., Natick, MA, USA). Lumbar vertebrae were manually segmented from preoperative and postoperative CT images (Fig. 2a). Additionally, screws and cages (whenever present) were segmented from the postoperative images (Fig. 2b). Image segmentation was performed by a professional segmentation and labeling company (Labelata GmbH, Zürich, Switzerland). The screw size (diameter and length) was determined from the surgery documentation if available, otherwise, the parameters were estimated from segmented screws. Because of image artifacts, the accuracy of the segmented screws was not sufficient to directly incorporate them as 3D surface models into the FE models. Instead, available template screw models from the MUST Pedicle Screw System (Medacta International SA, Castel San Pietro, Switzerland) were registered onto the segmentations with a rigid iterative closest point algorithm [3]. The Hounsfield unit (HU) distribution in the preoperative CT images was used to extract the heterogeneous distribution of bone mechanical properties within the intact vertebrae according to the method presented in Widmer et al. [31] ( Fig. 2c, Table 4). Variable tube voltage used for CT acquisition was accounted for by correcting the HU according to the results described in Afifi et al. [1]. For the 16 patients where upper body radiographs were available, musculoskeletal models were created and analyzed to determine the joint reaction force acting on the fusion construct (Fig. 2d).

Single-vertebra models
Pull-out finite element models were created for the uppermost vertebra of each patient following the experimentally validated procedure described in Widmer et al. [31] (Fig. 2e, Table 1). For every vertebra, there were two separate pull-out models: one with the left and one with the right screw. Screw extraction from the bone was simulated with constant velocity to obtain the force necessary to pull the implant out. In addition to axial screw pull-out, a second single-vertebra model type was analyzed. In this kind of model, both screws were concurrently present in one vertebra and while the implants were fixed, an imposed displacement was applied on the lower vertebral endplate in cranio-caudal direction (pulling the vertebra down with respect to the screws).

Multi-vertebrae models
Multi-vertebrae models were generated to address the question if additional insight is attainable with a more holistic finite element model of the spine. Specifically, model complexity was increased by including further implants (i.e. tulips, rods, and cages) and intervertebral discs. Patientspecific musculoskeletal models based on the individual anatomy, alignment, and mass distribution were generated to determine the physiological loading acting on the joints in the lumbar spine during standing. The resulting joint reaction forces were applied as loading conditions on the finite element models. For every patient with available bi-planar radiographs, one multi-vertebrae finite element model was simulated ( Fig. 2f; Table 1). The patient models included exclusively vertebrae of the lumbar spine. In-depth description of multi-vertebrae model generation can be found in the Appendix (Section A).
The explicit finite element solver RADIOSS (version 2017.2, Altair Engineering Inc., Troy, USA) was used for the analysis. For the single-vertebra simulations only the uppermost, newly instrumented vertebra was analyzed. For the multi-vertebrae simulations, the screw loosening status of all the lumbar vertebrae instrumented at index surgery was evaluated Table 1 Characteristics of patients included in the single-vertebra and the multi-vertebrae modeling groups Alignment parameters (lumbar lordosis, thoracic kyphosis, and sacral slope) were derived from annotated bi-planar radiography images [11]. The screw diameter was either 6 mm or 7 mm and the screw length was 40 mm, 45 mm, 50 mm, 55 mm, or 60 mm. p.p.: per patient. The indicated values correspond to: mean ± standard deviation

Output parameters and statistical evaluation
Various parameters were assessed with respect to their ability to predict screw loosening. For pull-out simulations, the considered output was the maximal strength before screw fixation failure ( Fig. 2e; Widmer et al. [31]). For pull-down simulations, the force at peak velocity was considered. For multi-vertebrae simulations, the volumetrically weighted average of HU values, von Mises stress (Stress VonMises ), and loading factor in elements around the screw were evaluated. The loading factor is the relationship between the local predicted stress weighted by the local bone quality: In every vertebra, the left and the right screw had the same loosening condition: either they were both loose or both fixed. Therefore, the mean value between the left and the right pedicle was considered for the pull-out strength and the outcome from the multi-vertebrae models. This way, for every assessed parameter there was one data point per vertebra.
The Wilcoxon rank-sum test was used to evaluate if the outcome parameters were different between the vertebrae with loosening and those without loosening. Logistic regression models were fitted to the data and the respective goodness-of-fit was assessed with McFadden's Pseudo-R 2 (R 2 Pseudo ). Next, for every regression model, the receiver operating characteristic (ROC) curve was computed. Significance level α was set to 0.05.

Single-vertebra models
In the patient cohort considered for single-vertebra models, screw loosening was diagnosed for 22 screws in 11 vertebrae, while 30 screws in 15 vertebrae were assigned to the control group. Detailed description of the two cohorts can be found in Table 5.
The maximal pull-out force was not significantly different between the loosening and the control group (Z = 0.26, p = 0.7953), while this was the case with respect to the screws' resistance to loading in downwards direction (Z = 2.44, p = 0.0147). Results are listed and shown in Table 2 and Fig. 3, respectively.
Description of the predictive ability of the logistic regression fitted on the data can be found in Table 3 and these results are shown in Fig. 4 for the pull-out strength.

Multi-vertebrae models
Out of the 16 considered patients, 9 had at least 2 loose screws. Overall, 26 screws in 13 vertebrae were loose and therefore assigned to the case group. 36 screws in 18 vertebrae were classified as non-loose (control group). Details about the two cohorts are listed in Table 6.
The distributions of the case and control groups were significantly different for all three parameters (Hounsfield units: Z = 2.38, p = 0.0172; Stress VonMises : Z = 2.18, p = 0.0291; loading factor: Z = 3.30, p = 0.001). Results are listed and shown in Table 2 and Fig. 5, respectively. Table 3 contains quantitative information on the predictive ability of the three outcome parameters. Figure 6 shows the corresponding results for the loading factor.

Discussion
Health issues affecting the spine are ubiquitous in all age groups but mostly affect the elderly. With the global share of people being 65 years or older increasing from 6% in 1990 to a projected 16% in 2050, strategies improving longterm surgical outcomes become increasingly important and urgent [26]. A promising approach for achieving this is the improvement of surgical planning through a more accurate patient-specific risk assessment. To this end, we evaluated the predictive capacity of various outcome parameters of computational biomechanical models.
The simulation results from pull-out finite element models of 52 screws showed a poor capacity of distinguishing between loose and non-loose screws (Fig. 3a, Fig. 4, Table 3). Despite its frequent use for testing purposes, it Fig. 2 Schematic depiction of the steps involved in patient-specific model generations. a Vertebral bodies and relevant implants were segmented from CT images acquired before (Pre-OP CT) and after (Post-OP CT) the index surgery. b Screw size and position within the vertebra were detected from CT-segmentations and/or surgery documentation. This allowed determining an appropriate 3D surface model of the screw. c Material parameters (stiffness, density, strength) of the vertebral bone were extracted from pre-operative CT images and assigned to the volume elements of the model. d Annotated radiographs of the erect posture were used for the generation of individualized musculoskeletal models. Analyses of these models allowed computing the joint reaction force above the fused vertebral segments. e Pull-out and pull-down strength were computed for the uppermost instrumented vertebra of each patient. f Hounsfield units, Stress VonMises , and the loading factor in screw vicinity were determined from multi-vertebrae models loaded with the computed output from musculoskeletal models ◂ 1 3 is debated whether axial screw pull-out contributes significantly to the occurrence of screw loosening [9,14]. Craniocaudal cyclic loading is considered to be a more likely contributor to screw loosening development [4]. For this reason, in addition to pull-out models, we created finite element models of the same vertebrae and screws and loaded them along the body axis. The resistance to initial failure was significantly different between the case and the control screws (Fig. 3b). Overall, this data suggests that with respect to screw loosening, the screws' resistance to cranio-caudal load is more relevant than their fixation strength in the axial direction.
In this study, multi-vertebrae finite element models were successfully created for 16 patients with 62 screws added during index surgery. The models incorporated a variety of different patient-specific factors related to surgery. To simulate the load acting on the biological tissue and the implants in the lumbar spine region, the joint reaction forces determined with patient-specific musculoskeletal simulations were applied on the finite element models. The main output parameter of these simulations was the stress around the screws. Several studies have shown that there is an association between bone mechanical properties and the occurrence of screw loosening [5,28]. Also, it has been hypothesized that high stress in vertebral bone is associated with pedicle screw loosening [18]. The outcomes of the current study strengthen this view of things by providing further evidence for an important stress-dependency of screw loosening. Specifically, the results suggest that both the CT-values and the simulated stress within the bone are different around loose and fixed screws (Fig. 5a, b). However, the loading factor defined by the relationship between the predicted stress in the tissue and its resistance to failure was able to distinguish best between the case and control group (Fig. 5c). This factor is therefore a very good indicator of how close to structural failure the bone around the screw is and the results imply that both the local mechanical quality of the bone, as well as the predicted magnitude of its loading, are important parameters to assess the loosening risk. A Pseudo-R 2 of 0.358 suggests a good fit between the data and the regression model (values ranging between 0.2 and 0.4 signify a good fit [17]). Therefore, the regression model could provide indicative values about the chances of screw loosening occurring based on the loading factor (Fig. 6d, Table 3). According to our results, a loading factor higher than 25% is a likely indicator of screw loosening, while a value lower than 15% is associated with negligible chances of fixation failure. In clinical practice, a threshold established based on the optimal combination of sensitivity and specificity could be of interest. For this purpose, the chances of false-negative results (i.e. predict screw fixation when the screw becomes loose) need to be weighed against the risk of false-positive results (i.e. predict screw loosening even though the screw is fixed) (Fig. 6f).
Finite element model generation is usually a lengthy process. With exception of manual segmentation of CT images and annotation of radiographs, the here presented models were generated automatically. Through the automation of modeling methods, the inclusion of computer models as part of preoperative planning is becoming a realistic approach in future patient care.
Our study also had some limitations. First, the size of the dataset was relatively small. Yet, discernable effects were clear and promising. In the future, a longitudinal prospective study should be carried out to assess the reproducibility and the actual predictability of screw loosening based on a loading factor. Also, as is always the case when modeling biological tissues, some simplifications were present in the  (MPa) 0.628 (0.520-0.998) 0.367 (0.282-0 . 3 Simulated a pull-out strength and b"pull-down" strength for models incorporating loose (Case) and non-loose screws (Control) models. For example, ligamentous tissue was not considered. However, the stabilizing contribution of ligaments in a neutral position can be considered to be small [30]. For several patients, the sacrum was also included in the fusion but this was not captured in our models, where the most caudally modeled vertebra was L5. Among the patients, there was also considerable variability in features like the levels which were fused, the used type of cage, the length of the instrumentation, the screw size, and the properties of the acquired images. By introducing various spine-specific aspects in our biomechanical models, we could account for many of these variable elements. Also, screw surface models of only one manufacturing company were used, while screws from various providers were implanted in the considered patients. Still, the outer screw diameter and the screw length were properly rendered since size-related features have been shown to possess a major impact on the screw holding power [14]. Finally, several factors have been described to be associated with the etiology of failed pedicle screw fixation: among them, aspects linked to biomechanics (stress shielding, microfracture, strain at the implant-tissue interface), the presence of wear debris, and chronic infection [13,15]. Here we were solely focusing on biomechanical aspects. Yet, the results indicate that major aspects involved in the development of screw fixation failure were captured, suggesting that they are biomechanical in nature.

Conclusion and outlook
In conclusion, the immediate potential of the here presented methodology is the customized adaptation of therapy, such as the employment of expandable or cement-augmented screws. Further, the approach might allow assessing possibilities for optimizing screw positions and to evaluate newly designed implants with respect to their capacity of increasing the loading score. This ratio of local predicted stress to bone strength is also a candidate for the assessment of other implant-bone interactions. In the next step, a prospective study should be conducted for validation of the used pipeline and further investigation of loosening predictability. Finally, based on the results of previous studies and of the data resulting from this work, it seems advisable to reconsider the value of commonly conducted pull-out tests for the study of pedicle screw holding power within the spine.

Appendix A: Multi-vertebrae FE model generation
Vertebral level-specific statistical shape models (SSM) were non-rigidly fitted onto surface models of the vertebrae obtained by segmentation from CT scans [7]. The correspondence properties of the SSM meshes facilitated an automated and consistent positioning of soft tissue and the rendering of posterior structure resections. Between the vertebrae, either a cage or an intervertebral disc was positioned. When no cage was present in the postoperative CT scan, an intervertebral disc was placed between the vertebrae by extruding it from the endplates of the two adjacent vertebral bodies (with slight bulging in radial direction). Anatomical studies were used to define the size and the relative position of the nucleus pulposus within the annulus fibrosus [19,22].
Resections of bony parts in the posterior vertebral body due to surgical interventions (i.e. decompression or facetectomy) were modeled in a patient-specific and automatized approach as well.
Tulips and rods were automatically positioned based on the pedicle screw orientation and screw head positions. By means of a spline fitted through the tulips, rods with a 5.5 mm diameter were added to the model along the fused levels. Whenever present, the segmented interbody cages were incorporated into the models after simplifying their surface to a convex hull.
All implants were modeled as non-deformable components because of their considerably higher stiffness compared to the biological tissue. The vertebrae were meshed with tetrahedral elements surrounded by a cortical shell. Bone was modeled as an elasto-plastic material with heterogeneous properties replicating the distribution of the intensity values within the CT images [31]. For vertebrae that were already instrumented with screws at the time of preoperative image acquisition (i.e. before extension surgery), uniform material properties were ascribed to the trabecular bone (see Table 4). The material behaviors of the nucleus pulposus and the annulus fibrosus were characterized as hyperelastic material and fiber-reinforced hyperelastic material, respectively. Material property specifications are listed in Table 4.
The contact between the intervertebral discs and the bone was modeled as a tied interface and the contact surfaces between implants and bone were modeled with an impact interface based on the penalty method. Facet joints were introduced as zero-friction contacts.
In the multi-vertebrae finite element models, the lower endplate of the most caudally positioned vertebra was fixed along all degrees of freedom. Magnitude and direction of joint reaction forces obtained from semi-subject-specific Fig. 6 Assessment of loading factor as a variable for distinguishing between loose (Case) and non-loose (Control) screws. a-c Histogram and continuous probability distribution fitted on the data from the two considered cohorts. d Logistic regression. e Receiver operating characteristic curve. f Sensitivity and specificity of the prediction model with respect to the loading factor threshold. A threshold value of 24.9% is associated with a specificity of 1.00 and a sensitivity of 0.46. Alternatively, a threshold of 21.5% would lead to an almost equal specificity and sensitivity (i.e. 0.78 and 0.77, respectively) and a threshold of 15.3% would be linked to a sensitivity of 1.00 (and a specificity of 0.61) Table 4 Material parameters used in the finite element models Bone was modeled as a temperature and strain-rate-independent Johnson-Cook material. The nucleus pulposus was modeled as a Mooney-Rivlin material. Fibre-reinforced hyperelastic properties were assigned to the annulus fibrosus with the ground substance behavior being characterized as a Mooney-Rivlin material ν poisson's ratio; ρ density, E young's modulus, a yield stress, b hardening modulus, n the hardening exponent, IVD intervertebral disc, c 10 and c 01 constants of the Mooney-Rivlin hyperelastic material, a 1 and a 2 constants describing the fiber behaviour, α fiber orientation  multi-rigid body simulations were applied on the disc above the most cranially positioned vertebra. The standing patient position was considered for this purpose and the information required to render the individual patient anatomy, alignment, and mass distribution were derived from weight measurements and annotated bi-planar radiography images acquired with the EOS system (EOS imaging, Paris, France). A detailed description of the multi-rigid-body modeling approach can be found elsewhere [11,29]. Fusion of the segments was simulated by turning the segments affected by the spondylodesis surgery into rigid bodies, i.e. not allowing any relative motion between the fused vertebrae of a subject [16].