The Simulation of Muscles Forces Increases the Stresses in Lumbar Fixation Implants with Respect to Pure Moment Loading

Simplified loading conditions such as pure moments are frequently used to compare different instrumentation techniques to treat spine disorders. The purpose of this study was to determine if the use of realistic loading conditions such as muscle forces can alter the stresses in the implants with respect to pure moment loading. A musculoskeletal model and a finite element model sharing the same anatomy were built and validated against in vitro data, and coupled in order to drive the finite element model with muscle forces calculated by the musculoskeletal one for a prescribed motion. Intact conditions as well as a L1-L5 posterior fixation with pedicle screws and rods were simulated in flexion-extension and lateral bending. The hardware stresses calculated with the finite element model with instrumentation under simplified and realistic loading conditions were compared. The ROM under simplified loading conditions showed good agreement with in vitro data. As expected, the ROMs between the two types of loading conditions showed relatively small differences. Realistic loading conditions increased the stresses in the pedicle screws and in the posterior rods with respect to simplified loading conditions; an increase of hardware stresses up to 40 MPa in extension for the posterior rods and 57 MPa in flexion for the pedicle screws were observed with respect to simplified loading conditions. This conclusion can be critical for the literature since it means that previous models which used pure moments may have underestimated the stresses in the implants in flexion-extension and in lateral bending.


INTRODUCTION
Spinal fixation has become a consolidated treatment for severe degenerative spinal disorders such as adult scoliosis, fixed sagittal imbalance, and high-grade spondylolisthesis (Ha et al., 2008;Casaroli et al., 2019;Galbusera et al., 2020). Despite the generally high success rates of spine surgeries nowadays, biomechanical complications such as hardware failure and loosening are relatively frequent (Kuklo et al., 2001;Tsuchiya et al., 2006;Kebaish, 2010). The literature indeed includes several studies in which different spinal fixation techniques have been investigated in terms of stresses and strains in the instrumentation (Fleischer et al., 2012;Burns et al., 2016;Sutterlin et al., 2016;Casaroli et al., 2019;Cunningham et al., 2019;Zhang and Zhu, 2019), which can be considered as indicators of the risk of biomechanical complications.
Most of the in vitro and finite element (FE) studies have been conducted using simplified loading conditions, usually consisting of pure moments, in some cases in combination with compressive forces, which are easier to implement than more realistic conditions involving muscle forces (Stokes and Gardner-Morse, 1995;Wilke et al., 2001;Kim et al., 2014;Zhang and Zhu, 2019). Although several studies confirmed that applying simplified or realistic loading conditions provides the same motion in intact spines (Rohlmann et al., 2009;Han et al., 2011;Zhang and Zhu, 2019), the effect on the instrumentation stresses has never been documented.
As regards the identification of realistic loads, software for musculoskeletal (MSK) modelling such as AnyBody (AnyBody Technology, Aalborg, Denmark) and OpenSim (Stanford University, Stanford, US) provide pre-built models able to predict the muscle forces for any imposed motion of the body segments using algorithms based on inverse dynamics (Bruno et al., 2015;Bassani et al., 2017;Benditz et al., 2018;Liu et al., 2018;Liu et al., 2019). Such models are based on the equations of motion of rigid bodies and cannot therefore be used to estimate the stresses in the implants or in biological structures (Zhu et al., 2013;Arshad et al., 2016;Liu et al., 2018;Bassani et al., 2019). However, the computed muscle forces can be applied as loading condition from the MSK model to the FE model, from which detailed information regarding the hardware stresses can be extracted (Bassani et al., 2017;Liu et al., 2018;Bassani et al., 2019;Liu et al., 2019). Such strategy, coupling FE and MSK models, has never been used to investigate the instrumentation stresses after spinal fixation, and can represent an advantageous approach to determine if the simplified loading conditions consisting of pure moments used in the majority of the available studies are good enough to accurately estimate the hardware stresses in physiological conditions. The aim of this study is therefore to develop and validate coupled FE-MSK models with and without instrumentation and to explore the differences, in terms of hardware stresses, between simplified (pure moment without follower load) and realistic loading conditions.

Intact Model
The three-dimensional (3D) geometry of T10-T12 thoracic vertebrae, lumbar vertebrae, and pelvis of the body model from the AnyBody Managed Model Repository (AMMR, version 2.0.0) in the standing posture was used to construct a FE model (Figure 1) made of linear tetrahedral elements after the surfaces were cleaned using MeshLab software (http://www. meshlab.net). Then, the triangular elements of the surfaces of the vertebral endplates were extruded in order to obtain a volume made by tetrahedral elements representing the discs, which were divided into nucleus pulposus and annulus fibrosus based on anatomical data from the literature (Zhong et al., 2014). The annulus included collagen fibers modeled with nonlinear springs; ligaments were also modelled with the same type of elements. The material properties of the bones, intervertebral discs and ligaments were obtained with a calibration procedure based on data reported in the literature (Schmidt et al., 2007). Pure moments of 7.5 Nm in extension, flexion, lateral bending, and axial rotation were applied as simplified loading condition to the upper endplate of T10 through a set of rigid beam elements ( Figure 1). The acetabula were completely fixed during all the simulation. The range of motion (ROM) calculated at all levels were then compared to in vitro data in order to validate the FE model under simplified loading conditions (Cook et al., 2015;Lindsey et al., 2018).
Then, the obtained ROMs were used as inputs for the MSK model of the thoracolumbar spine with articulated ribcage (Figure 1) developed and validated by Ignasiak et al. (Ignasiak et al., 2016). The muscles simulated in this MSK model were transversus, spinalis, semispinalis, erector spinae, obliquus internus, obliquus externus, psoas major, multifidi, and quadratus lumborum. In this model, the rotational stiffness of the intervertebral joints was calibrated in order to match the linear moment-rotation behavior of the FE model, including the effects of all the joint structures (including facet joint and ligaments). This stiffness does not account for compressive FIGURE 1 | The combination of the intact FE model with the intact MSK model. Simplified loading conditions are applied to the FE model in order to obtain the motion. This motion is imposed as input to the MSK model and muscles forces that contribute to give that motion are predicted by this model. Muscles forces are then applied to the FE model to realistically simulate the loading conditions. The motion under realistic loading conditions is then obtained.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 745703 loading or coupling. The pelvis was constrained to the ground. Extension, flexion, and left and right lateral bending movements were simulated by imposing intervertebral rotations matching those calculated with the validated FE model under simplified loading conditions. The muscle forces were calculated for each simulated motion by inverse static analysis ( Figure 1). After that, the obtained muscles forces were modelled at each level (from T10 to the pelvis) in the FE model as concentrated loads to simulate the realistic loading condition (Figure 1), removing the pure moment which implemented the simplified loading conditions; the upper endplate of the T10 vertebra was subjected to a 3D translation equal to that predicted by the MSK model, using a set of rigid beam elements. Moreover, reaction loads from the MSK model were applied to the upper endplate of the T10. Such loading conditions were applied to the FE model for all the investigated motions: extension, flexion, left and right lateral bending. Finally, a validation of the FE model under realistic loading conditions was performed by comparing the ROMs obtained under realistic loading conditions ( Figure 1) and the one imposed as input for the MSK model. Ideally, for the validation the ROMs values had to be equal. Moreover, the reaction moment at the upper endplate of T10 was compared between MSK model and FE model in order to provide an extra validation.

Instrumented Model
From the intact FE model, an instrumented model was derived. This model included a posterior lumbar fixation in which pedicle screws and rods were inserted in the lumbar region between L1 and L5 vertebrae. 3D models of rods and screws were created with commercial software Solidworks (Dassault Systèmes, Waltham, MA, USA). Rods had a circular section with a diameter of 5.5 mm. Pedicle screws had a length of 45 mm and a diameter of 6 mm. The instrumentation was modelled in titanium with an elastic modulus of 110 GPa and a Poisson coefficient of 0.3. Simplified loading conditions and boundary conditions were the same as in the intact model ( Figure 1).
Similarly, an instrumented MSK model was created from the intact one. Spinal fusion was modeled by introducing rigid kinematic and kinetic constraints from L1 to L5 vertebrae, guaranteeing rigid connection (no relative motion between vertebrae) and full force and moment transmission (Ignasiak et al., 2018;Ignasiak, 2020). The validation of this model was performed against an in vivo study by Rohlmann (Rohlmann et al., 1997), in which bending moments in the posterior rods were measured by means of instrumented implants. Extension, flexion, and left and right lateral bending movements were simulated imposing the intervertebral rotations calculated with the instrumented FE model under simplified loading conditions. From this model, muscle forces were obtained for the four different motions, as for the intact model ( Figure 1).
Realistic loading conditions were then simulated in the instrumented model as for the intact model ( Figure 1). The FE model under realistic loading conditions was then validated by comparing the intervertebral motion calculated with it with the one imposed as input for the MSK model. As for the intact model, the reaction moment at the upper endplate of T10 was compared between those predicted by the MSK and the FE models in order to provide an extra validation.

Validation and Output
The outputs of the intact models that were calculated for two types of loading conditions were: 1) the ROM of L1-S1 vertebrae and SIJ; 2) the reaction moment at the upper endplate of T10 joint. The 1) values were used to validate the intact FE model under simplified loading conditions (Cook et al., 2015;Lindsey et al., 2018). 1) and 2) values were used to validate the intact FE model under realistic loading conditions.
The outputs of the instrumented models that were calculated for two types of loading conditions were, in addition to those of the intact model, also: 3) the maximal von Mises stresses in L1-L5 pedicle screws; 4) the maximal von Mises stresses in the posterior rods between the pedicle screws in L1 and L5. The 1) values were used to validate the instrumented FE model under simplified loading conditions. The 1) and 2) values were used to validate the instrumented FE model under realistic loading conditions. The 3) and 4) values were used to compare the hardware stresses between the finite element model under simplified loading conditions and those under realistic loading conditions.

Validation of the Intact Model
The ROMs calculated with the intact FE model under simplified loading conditions showed a tendency toward a higher rigidity with respect to the literature (Cook et al., 2015;Lindsey et al., 2018) (Figure 2); despite this, the predicted values were inside the standard deviations of the in vitro data in flexion-extension and axial rotation, except for the L1-L2 ROM in axial rotation which was higher than the corresponding experimental finding. In lateral bending, the ROM was approximately equal to the lower limit of the standard deviations of the in vitro studies, except for the SIJ ROM which was in agreement with the value found in the literature.
The ROM calculated with the intact model under realistic loading conditions revealed values similar to those calculated with the intact model under simplified loading conditions in flexion-extension and lateral bending, as expected ( Figure 3A); nevertheless, some relatively small differences were observed. For instance, negligible differences up to 0.7°were found in lateral bending. The reaction moment at the upper endplate of T10 showed a maximal relative difference of 3.5% between the FE model and the MSK model.

Validation of the Instrumented Model
The instrumented FE model showed that the ROM of instrumented levels was negligible with respect to the case without instrumentation, as seen in computational and in vitro studies (Rohlmann et al., 2007;Dmitriev et al., 2008). The reaction moments at the instrumented levels obtained with the MSK model were inside the range of fixator load measurements assessed in vivo (Rohlmann et al., 1997), demonstrating the plausibility of the results.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 745703 The comparison of the ROMs under simplified and realistic loading conditions revealed very similar values in the two cases, with small differences up to 0.3°in flexion-extension ( Figure 3B). The reaction moment at the upper endplate of T10 showed a negligible difference between the FE model and the MSK model.

Stresses in the Pedicle Screws
In extension, flexion, lateral bending on the left side, and lateral bending on the right side, the maximum stresses on the L1-L5 pedicle screws were higher when realistic loading conditions were applied to the model ( Figures 4A,B). In extension, the maximum stresses for the L1-L4 pedicle screws were higher than the values found with simplified loading conditions. With respect to simplified loading conditions, an increase of 44 MPa (5.6% with respect to the yield stress) for the right pedicle screw in L1 was found under realistic loading conditions. The maximum stresses for the L5 pedicle screws showed more comparable values, with changes up to up to 3 MPa. In flexion, results similar to the extension case were found, but with bigger differences; the highest difference was found for the left pedicle screw in L3 (7.2% with respect to the yield stress) ( Figure 4A). For lateral bending, higher stresses on all pedicle screws were predicted when realistic loading conditions were used; the maximal difference resulted to be 45 MPa (5.7% with respect to the yield stress) ( Figure 4B).

Stresses in the Posterior Rods
The maximal stresses on the left and right posterior rods had similar values among simplified and realistic loading conditions ( Figures 4A,B). Similar to the pedicle screws, the maximal stresses on the posterior rods were highest when realistic loading conditions were simulated. In flexion-extension, increases exceeding 99% were found with a maximal difference of 40 MPa in extension (5.1% with respect to the yield stress) ( Figure 4A). The same trend was observed in lateral bending, but with smaller differences up to 33 MPa (4.2% with respect to the yield stress) ( Figure 4B). FIGURE 2 | Validation of the intact FE model under simplified loading conditions. Predicted ranges of motion of L1-S1 and sacroiliac joints of the intact model under simplified loading conditions in flexion-extension (left), lateral bending (middle) and axial rotation (right), as compared with data from in vitro experiments, shown as mean and standard deviation (Cook et al., 2015;Lindsey et al., 2018).

DISCUSSION
This paper presents a preliminary biomechanical comparison between simplified and realistic loading conditions to determine the hardware stresses in a spinal fixation model, aimed at investigating if the implementation of a more realistic loading scenario has the potential to significantly affect the results. Simplified loading conditions consisting of pure bending moments, in some cases combined with compressive loads, are often preferred for in vitro and computational testing of spine specimens being easier to implement than more realistic conditions involving muscle forces, while ensuring better reproducibility. Nevertheless, this study demonstrated that applying simplified loads can result in an underestimation of the hardware stresses in the instrumented models, with potential implications about the safety of the implants (Wilke et al., 2001). In this study, the metrics used to evaluate the importance of the loading conditions were the stresses in the L1-L5 pedicle screws and those in the posterior rods. The results showed that realistic loading conditions increase the stresses in the hardware, up to 57 MPa in flexion (Figures 4A,B). For the pedicle screws in L1-L5, the stresses were higher for all four motions when realistic loading conditions were used. For the posterior rods, higher stresses were found for all conditions with realistic loading conditions ( Figures 4A,B). This is a very interesting result because the most common type of biomechanical failure of the fixation is indeed rod breakage (Yamanaka et al., 2015).
Since this approach (coupling FE and MSK modelling) is reported here for the first time for an instrumented spine model, no comparison of the current results with similar existing data can be performed. Despite this, previous studies investigated the validity of simplified loading conditions by comparing in vivo measurements with in vitro experimental tests. Wilke et al. compared the loads acting on an internal spinal fixator in 10 patients (in vivo) with an equivalent in vitro simulator under the application of pure bending moments (Wilke et al., 2001). They found good agreement for the loads acting in the internal fixator for axial rotation and lateral bending. For flexion ad extension, reasonable agreement was found only for the healthy spines instrumented with fixators, while for specimens in which a bone graft was implanted in the intervertebral space a lower agreement between in vivo and in vitro data was found. As regards the proposed FE-MSK approach, other studies in the literature have exploited such combination, but not to investigate instrumented scenarios. For instance, in a computational study by Liu and colleagues a coupled FE-MSK model was used to investigate the load-sharing in the lumbosacral spine, where muscle forces, as predicted by a MSK model, were used as loading conditions for the FE model (Liu et al., 2018;Liu et al., 2019).
This present study has some limitations. Axial rotation was not investigated; however, it is worth considering that the previous papers exploiting coupled FE-MSK models only considered flexion-extension motion, therefore the simulation of lateral bending motion still constitutes an advance with respect to the state-of-the-art (Rohlmann et al., 2006;Liu et al., 2018;Liu et al., 2019). Another limitation was the translation imposed to the most cranial vertebra, which was chosen as boundary condition after verifying that a pure load-controlled simulation driven by muscle forces and the reaction force calculated by the MSK model at the most cranial joint did not lead to convergence. This is however an improvement with respect to the method used by Liu et al., in which the L1 vertebra was subjected to a translation in the direction of the force equal to the one predicted by the MSK model in order to ensure quick convergence of the FE model, but the applied translation was adjusted if the difference of the reaction force in the MSK model and in the FE model was greater than predefined tolerances (Liu et al., 2018;Liu et al., 2019); in the present study, no adjustment was necessary. Besides, only one instrumented configuration was Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 745703 presented in this study while various instrumentation strategies are used in the clinical practice, depending on the condition of the spine of the patient; the simulation of other common configurations is indeed ongoing. Another limitation was that the simplified loading conditions included only pure bending moments without a compressive load mimicking the body weight, such as for example a follower load (Patwardhan et al., 2003). This simplification could justify the difference found in the hardware stresses between the two types of loading conditions, but it should be said that the use of pure moment without compressive loads is very common in vitro and computational studies investigating the stresses in the implants (e.g., Wilke et al., 2001;Galbusera et al., 2020), and pure moments have been recommended as the preferred method to test spinal implants in standardized laboratory tests (Wilke et al., 1998). However, we acknowledge that additional studies should be done using a follower load in combination with pure moments in order to simulate another commonly implemented set of simplified loading conditions. Moreover, it should be noted that the motion imposed as input for the MSK model was equal to the validated output of the corresponding FE model under simplified loading conditions, being therefore possibly different from the physiological motion of the spine. In this respect, gait analysis and fluoroscopy can be potential alternatives to determine the motion of the spine to be used as input for the MSK model (Haddas et al., 2018;Breen and Breen, 2020). Finally, it should be noted that the stiffness imposed to the MSK model does not account for compressive loading, but only for the pure moment applied to the FE model. In conclusion, hardware stresses resulted markedly higher when realistic loading conditions, consisting of muscles forces applied to several vertebrae from T10 to the pelvis, are used instead of simplified loading conditions. This conclusion has relevant biomechanical implications since it means that previous models which used pure moments may have underestimated the stresses in the implants in flexionextension and in lateral bending. Further studies, using different spinal fixation techniques and follower load, need to be done in order to understand if this combined method is more useful than a simplified one to predict implants failure.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MP, TB, and FG contributed to the conception and design of the study. MP, TB, and FG performed the analysis and interpretation of data; MP wrote the first draft of the manuscript; MP, TB, TV, and FG contributed to manuscript revision and approved the submitted version; MP, TB, TV, and FG contributed to an administrative, technical, or material support.

FUNDING
The study was supported by the Italian Ministry of Health ("Ricerca Corrente").