Phase field modelling and simulation of damage occurring in human vertebra after screws fixation procedure

The present endeavor numerically exploits the use of a phase-field model to simulate and investigate fracture patterns, deformation mechanisms, damage, and mechanical responses in a human vertebra after the incision of pedicle screws under compressive regimes. Moreover, the proposed phase field framework can elucidate scenarios where different damage patterns, such as crack nucleation sites and crack trajectories, play a role after the spine fusion procedure, considering several simulated physiological movements of the vertebral body. A convergence analysis has been conducted for the vertebra-screws model, considering several mesh refinements, which has demonstrated good agreement with the existing literature on this topic. Consequently, by assuming different angles for the insertion of the pedicle screws and taking into account a few vertebral motion loading regimes, a plethora of numerical results characterizing the damage occurring within the vertebral model has been derived. Overall, the phase field results may shed more light on the medical community, which will be useful in enhancing clinical interventions and reducing post-surgery bone failure and screw loosening.


Introduction
Originated from mathematical techniques based on Γ−convergence [1,2] and tailored for the approximation of free discontinuity problems [3,4], the phase field regularization of brittle fracture proposed by Francfort and Marigo [5] has attracted a remarkable attention within the computational fracture mechanics community over the last decade.As compared to other computational methods for damage and fracture simulation in materials and components, such as for instance the Crack Band Model [6], the Smeared Crack Model [7], non-local and diffuse damage models [8,9], or gradient damage models [10], the phase field approach offers an elegant solution for problems involving linear elastic fracture mechanics.This solution is pursued through an energy minimization which, in the Γ−convergence limit, consistently reproduces the Griffith theory of fracture.The phase field approach to fracture, further analysed in [11,12], has been applied in a considerable series of works proposing comparisons with other non-local damage models and discussing several detailed aspects regarding the finite element implementation [13][14][15][16][17][18][19].In this context, it is worth recalling the fundamental contribution by Miehe and co-workers [20,21] that were the first to propose a robust finite element implementation of the phase field for brittle fracture, specialized to account for damage irreversibility and based on a suitable degradation mechanism able to simulate situations involving tensile stress states.
The state-of-the-art literature on phase field models clearly shows that the approach is mature for technical applications.In this regards, Tannè et al. [22] have recently assessed the capabilities of the phase field approach (see in particular the AT1 and the AT2 models) to predict crack nucleation from V-notches and from points with stress concentrations.Phase field modeling has numerous applications, especially in the field of biomedical engineering, particularly in the study of bone fractures.In this regard, fractures in the vertebrae can be challenging to analyze and treat due to their complex shapes and locations [23].However, phase field modeling can help by simulating the fracture process and shedding light on the underlying mechanisms involved.In particular, phase field modeling can be used to investigate the impact of vertebral fixation and screw placement on the healing of vertebral fractures.By incorporating the presence of screws and other fixation devices into the simulation, researchers can evaluate the mechanical stability of the fracture site and predict the potential for screw loosening or failure.Additionally, researchers can gain insights into the effects of different placement strategies on the fracture propagation process.
Quantitative and qualitative finite element models are being developed to evaluate damage patterns and predict crack propagation in biological tissues.So much so, recently, the capability of a plethora of phase field implementations for modelling and analysing crack growth in bone tissue has been successfully applied [24], thus phase field theory seems a promissory approach to assess bone fracture patterns.In this context, the biomechanical problem of characterizing damage and predicting crack trajectories in a human vertebra after the fixation procedure of pedicle screws (see Fig. 1(a)) can be tackled utilizing a phase field method.
Metallic pedicle screws are used in spine fusion procedures when the intervertebral discs were damaged by aging or any trauma, causing the vertebrae to rub against each other and compress the nerves that pass through them.Spinal fusion joins two or more diseased vertebrae together, preventing motion at the vertebral segment [25].The screws are inserted into the respective pedicle regions, connecting the screws through a vertical rod, providing a means for gripping onto a vertebral segment and limiting its motion, resulting in stable spine fixation.Bone damage and screw loosening may occur due to various post-operative events [26].Furthermore, the mechanical behavior of the pedicle screw insertion angle has been experimentally investigated [27,28].Nonetheless, finite element analysis allows for the assessment of numerous possible failure scenarios that may occur after the surgical procedure requiring screw insertion, as well as the character-ization of fracture damage and estimation of mechanical responses under varying screws fixation angles [29].Thereby, the impact of screw configuration angular parameters in fracture pattern and stress distributions were investigated via computational model considering a stress-based criterion [30,31].In the study by [23], the mechanical behaviour of instrumented metastatic vertebra in presence of pedicle screws was investigated.The findings showed that the prediction of mechanical properties was affected by the size and location of the metastasis.It was also concluded that a metastasis situated near the screws produces a higher fracture load response compared to a metastasis far from the screws.As a consequence, a finite element phase field model is hereof implemented in order to validate numerically the existing outcomes in the literature resulting from the variation of the screws insertion angle in a human vertebra, with the aim of further developing a coupled model with two or more fused vertebrae.
The paper is organized as follows: Section 2 provides the CAD models of the screws and the vertebra, which are discretized and have their material properties set.It also covers the virtual insertion of pedicle screws into the vertebra.Section 3 introduces the phase field finite element method strategy implemented in this study.The staggered and monolithic phase field models are detailed in Section 4. Section 5 presents the boundary conditions and parameter setup for the numerical phase field method.Additionally, it includes a mesh sensitivity analysis of the phase field framework applied to several discretized models of vertebra-screws.Finally, Section 6 presents the main results in terms of damage patterns.It considers various screws insertion trajectories for different vertebral movements and includes a discussion comparing the outcomes with other relevant works in this field.

Models and materials
The L4 lumbar vertebra model was obtained from a Computer Tomography (CT) scan images from the spine of 49-year-old female patient, as described in [25].The solid cylindrical pedicle screws have length of 40 mm.The major and minor diameters of screws are 6.5 mm and 4.3 mm, respectively (please see [25,29] for additional details).The CAD models of the pedicle screws virtually inserted in the L4 vertebra are depicted in Fig. 1.The phase field model has been evaluated for three different screws fixation paths, indicated by the vector ⃗  = ( 1 ,  2 ).Here,  1 represents the insertion angle in the cranio-caudal (CC) direction, and  2 indicates the insertion angle in the mediolateral (ML) direction [25], as shown in Fig. 1(a).In particular, the mesh sensitivity analysis is conducted using the screws insertion trajectory combination ⃗  = (−5, 0), which means that the screws are fixed with a negative angle of −5 ∘ in the CC direction and a neutral angle of 0 ∘ in the ML direction.
To import the combined vertebra-screws CAD model into the finite element environment FEniCS [32], uniform tetrahedral meshes are constructed.This is achieved by transferring the assembled STL model's triangular mesh into the pre-processing software HyperMesh [33], which generates the required tetrahedral solid.The solid is then assessed and converted into a MSH file using the Gmsh platform [34].Multiple refined meshes have been created, ranging from 60 thousand elements to 700 thousand elements.Please refer to Fig. 1(b) for visualization of these meshes.Bone density properties have been considered for the material of the L4 vertebra using CT data [35][36][37].In this regard, a constant Poisson's ratio of  = 0.3 was assumed, and isotropic and linear elastic material properties were adopted with a heterogeneous distribution of the Young's modulus .This distribution differentiates the cortical bone from the trabecular bone.For more details, please refer to [25,29].For trabecular bone, an elastic modulus ranging between 0 and 3 GPa was derived, while for cortical bone, the range was between 12 and 14 GPa, as shown in Fig. 2(a).In what regards to the pedicle screws, surgical procedures commonly use biomedical implants or replacements made of the titanium alloy Ti-6Al-4V [38][39][40].Therefore, the chosen mechanical properties for the pedicle screws are a Young's modulus of 110 GPa and a Poisson's ratio of 0.4.

Phase-field approaches to fracture with spectral decomposition
With reference to an arbitrary body occupying a domain Ω ∈ R   , with boundary Ω ∈ R   −1 , in the Euclidean space of dimension   , in which an evolving internal discontinuity Γ is postulated to exist, a material point is denoted by x and body forces by b : Ω → R   .Mixed conditions are prescribed along non-overlapping Neumann and Dirichlet regions of the boundary Ω  ∪ Ω  = Ω in the usual form where n denotes the outward unit normal to the boundary, u is the displacement field and  is the Cauchy stress tensor, while u and T are prescribed surface displacements and tractions.
The variational approach to brittle fracture, governing crack nucleation, propagation and branching, is set up through the definition of the free energy functional [20], embodying an additive decomposition between the elastic bulk energy Π Ω stored in the damaged body and the energy Π Γ necessary to nucleate and propagate a Griffith crack [41], defined as where   is the elastic strain energy density, function of strain , and   is the fracture energy, function of the displacement u and of the phase field variable .The latter parameter  ∈ [0, 1] is an internal state variable, ranging between 0 and 1 and representing isotropic damage, so that  = 0 is representative of the intact material, while  = 1 characterizes the fully damaged state.

The regularized variational formulation
Within the regularized framework of the phase field approach [12,20,42,43], the potential energy of the system can be decomposed into two terms as following where   (, ) is the energy density of the bulk, now function of the damaged parameter , and (, ∇) is the crack density functional, with ∇ denoting the spatial gradient operator.As a result, the total free energy density of the bulk ψ reads as The functional (, ∇) is assumed to be a convex function of  and its gradient ∇ and can be written, in agreement with the following expressions characterizing the AT2 model, respectively, as: where  0 stands for a regularisation characteristic length that can be related to the Young's modulus, the fracture toughness, and the tensile strength of the material, as specified in Section 5.
To avoid the development of damage in compression, so to allow fracture growth only under tensile stress states, the following 'tensile/compressive' decomposition is herein assumed for the energy density in the bulk   (, ) [20,[44][45][46][47] and included in the formulation: where () is a damage function that is assumed in the simple form where  is a residual stiffness (introduced to avoid ill-conditioning) and the positive and negative parts of the energy density are defined as where  and  are the Lamé constants, tr(•) denotes the trace operator and the positive and negative parts of the strain  ± are defined as follows.With reference to the spectral representation for the strain (with eigenvalues   and unit eigenvectors   ), denoted as the strain is additively decomposed as  =  + +  − , so that the tensile and compressive parts associated to  are and respectively, where the Macaulay bracket operator is defined for every scalar  as ⟨⟩ = ( + ||)/2.A standard derivation [48] leads Eq. ( 7) to the Cauchy stress tensor from the strain energy density: where I denotes the second-order identity tensor.The thermodynamic consistency of the above constitutive theory, in agreement with the Clausius-Duhem inequality, has been addressed in [20].

Weak form of the variational problem
The weak form corresponding to the phase field model for brittle fracture can be derived following a standard Galerkin procedure.In particular, the weak form of the coupled displacement and phase field damage problem according to Eq. ( 4) is: where v is the vector of the displacement test functions defined on H 1 0 (Ω),  stands for the phase field test function defined on H 1 0 (Ω).Eq. ( 11) holds for every test functions v and .The external contribution to the variation of the bulk functional in Eq. ( 11) is defined as follows:

Staggered solution scheme
The mechanical problem can be stated as: given the prescribed loading condition u  and T  at step , while the phase field problem is formulated as: find where }︁ is the strain history function, accounting for the irreversibility of crack formation [17,20].
To solve the quasi-static evolution problems for brittle fracture, isoparametric linear triangular finite elements are used for the spatial discretization, and a staggered solution scheme is considered.Staggered schemes based on alternate minimization exploit the convexity of the energy functional with respect to each individual variable u and  [49].
Here, an ad hoc developed solver has been implemented in the software FEniCS, see Alg. 1 for the staggered algorithm description.A series of benchmark tests taken from [20,50] has been carried out to validate the methodology.

Newton-Raphson procedure
Even if the mechanical problem has been split into Eqs.( 13) and ( 14), so that the phase field is reduced to a linear problem, nonlinearity still remains, because of the piece-wise linearity of the constitutive law, which includes a spectral decomposition of the strain.Therefore, a consistent linearization is required, so that the linear form defined by the residual can be written as: Given u  the current Newton-Raphson approximate solution at iteration , the correction u is therefore the solution of the following linear variational problem: find , and then iterate as u +1 = u  +u.The Jacobian entering the formulation is for details about the terms  − ,  + , please see [18].At this stage, it is fundamental to remark that, in order to predict crack trajectories in human vertebrae under tensile/compressive stress states, the phase field finite element method will be formulated.This formulation involves decomposing the strain energy density   (u, ) in Eq. ( 6) based on spectral diagonalization, as described in [20], into active and passive parts.This decomposition allows for the application of material response degradation only in tension.The variational formulation is then implemented Algorithm 1 Staggered iterative scheme for phase field fracture at a step  ≥ 1 1: Input: Displacements and phase field (u −1 ,  −1 ) and prescribed loads (u  , T  ): 2: Initialize (u 0 ,  0 ) := (u −1 ,  −1 ); 3: for  ≥ 1 staggered iteration do:

Monolithic solution scheme
A monolithic solver has also been also implemented in FEniCS.The tangent operator of the non-linear variational functional, is computed via the symbolic derivative derivative, the monolithic algorithm scheme can be visualized in Alg. 2. Algorithm 2 Monolithic iterative scheme for phase field fracture at a step  ≥ 1 1: Input: Displacements and phase field (u −1 ,  −1 ) and prescribed loads (u  , T  ): 2: Solve the coupled non-linear variational problem via Newton-Raphson iterative scheme:

Parameters calibration
The length scale  0 parameter is deeply inserted for modelling phase field, considering that for a sufficiently small length scale , the functional ( 11) converges to the total potential energy functional (2), in the sense that the global minimizers of Π  will also converge to that of Π.This entails that the length scale must be carefully chosen, rather than setting it arbitrarily.Having said that, on one hand, the characteristic length scale is related to the apparent material strength [51][52][53].Particularly, once the base material properties are attributed, namely Young's modulus , critical energy release rate   , then the characteristic length can be tuned as Consequently, the failure stress  max was considered according to [25,54,55], so the length scale can be evaluated through Eq. ( 18).On the other hand, the characteristic length scale can be interpreted, and therefore calibrated, as a structural factor.In particular, in a biological material, as cortical bone for instance, the length scale can be attributed to the lamellar microstructure, which may trap the crack tip within it, deflecting or stopping the leading edge of the crack [56].Moreover, structural changes in cortical bone due to ageing affect crack path and damage properties [35,57,58].
In essence, the characteristic length scale that smears the regularized crack can be obtained either from microscopic structural-related mechanisms [59,60], or from bone properties [61,62].As such, prior to determine the characteristic length scale through Eq. ( 18), a power-law equality is usually assumed by which the critical energy release rate   is derived from bone density properties, namely where the parameters have been set as  0 = 20.000MPa   0 = 0.01 N/mm, and  = 0.8 are the base elastic modulus, the base critical energy release rate, the bone ash density and the power-law exponent, respectively.For those reasons, when a mesh sensitivity analysis is being conducted in a biological material, the phase field parameters must also be carefully examined and validated.As such, it has been performed in [24,62], different values for the characteristic length scale have been considered, where the latter has calibrated phase field parameters after mechanical in-vitro experiments on human humeri bones.In addition, the critical energy release rate   plays a crucial role in phase field modelling for bone application, where throughout experiments, different values of   have been obtained at different aged cortical bones [57,61].Following these aspects, the critical energy release rate   and the length scale have been calibrated using Eqs.
(18) and 19, and validated from the experimental studies herein mentioned.

Boundary conditions, numerical implementation and mesh sensitivity analysis
Although spine models have been extensively exploited and validates in the literature [63], the present finite element phase field method focuses on a single vertebral body, as observed in other studies [29,64].In this study, the loadings are addressed on the inserted screws in the L4 vertebra, which will be simultaneously constrained.Throughout all the numerical implementations, the boundary conditions were set to replicate the principal movements permitted by the vertebral column.For this purpose, the pedicle screws head have been loaded to reproduce flexion (bending forward) and extension (bending backwards) (green arrows in Fig. 3), rotation (torsion/twisting) (red arrows in Fig. 3).The centre of the L4 vertebra has been assumed as the rotation axis centre for the twisting loading mode.Meanwhile, a compressive force of   = 5 N has been applied to the superior endplate of the vertebral body, while the inferior endplate has been constrained.Furthermore, the vertebral fracture patterns have been simulated by incrementally applying a quasi-static force to the screw heads for all the studied movements.Fundamentally, the total force applied at each loading step corresponds to 10% of the constant compressive load   .The numerical phase field procedures were implemented in a High-Performance Computing (HPC) system, utilizing parallel computing and dividing the model into three smaller tasks.The computational time ranged from 2 hours for the coarsest mesh to 30 hours for the finest mesh.Six cores were utilized, with an average of 200 GB of RAM per core.It is worth mentioning that the extension and flexion loading modes required more time for simulation, necessitating the use of additional cores and computational memory to improve the efficiency of these movement modes.
With the aim of achieving an appropriate mesh discretization for obtaining accurate numerical phase field responses for various physiological motions, a meticulous mesh sensitivity analysis has been conducted using the phase field finite element method.The vertebra-screws mesh model has been uniformly discretized at various refinements levels, ranging from 60.000 to 350.000.Fig. 4presents the numerical analysis comparing the outcomes in the flexion loading mode for different mesh refinements, with the screws insertion angle set as ⃗  = (−5, 0).Fig. 4(a) depicts the Load vs. Displacement curves juxtaposed, indicating that as the mesh becomes finer, the vertebra-screws model can bear more load.Additionally, Fig. 4(b) illustrates the relative error as the number of elements increases.According to [25], it can be noticed that mesh convergence is satisfactorily achieved at a relative error of 5%.Therefore, assuming a mesh size between 200.000 and 350.00 elements is acceptable for characterizing the fracture patterns in the simulated physiological movements.Nonetheless, in terms of characterizing the damage patterns, the present mesh sensitivity analysis revealed no significant distinctions among the different mesh refinements.Moreover, the computational cost aligns well with the numerical outcomes using a mesh with 200.000 elements.This is further supported by Fig 4(c) which displays the fractured volume as the loading step varies.These convergence analyses demonstrate good agreement with previous studies [25,29].In terms of fracture type, although the phase field model captured similar damage patterns for all mesh refinements, coarser meshes exhibited more rapid fracture spread within the cortical part.In accordance with the estimations made in [62], an examination of the role played by the characteristic length scale  0 is also developed here.Taking into consideration the previous mesh sensitivity analysis, with a mesh of 200.000 elements, a screws insertion angle of ⃗  = (−5, 0) in a flexion mode motion, Fig. 5 indicates that the apparent peak load from the phase field scheme increases as the characteristic length scale  0 is reduced, as expected based on previous results reported in the literature [65,66].

Characterization of loading modes
Continuing with the comparative investigation of the various screw insertion angles, namely ⃗  = (−5, 0), ⃗  = (+5, +5), and ⃗  = (−5, −5).Fig. 6 highlights how the screw configuration influences the mechanical responses of the vertebral body in different motion regimes.The Force vs. Displacement curves generated from the phase field finite element method are displayed in Figs.6(a), 6(c), and 6(e), representing the vertebral movements of flexion, extension, and torsion, respectively.In the three cases analyzed, the Load is computed as the vincular reaction on the bottom of the vertebra.In the cases of flexion and extension there is a linear elastic phase, followed by a softening due to the occurrence of fracture inside the vertebra, whereas in the case of rotation, the vincular reaction on the bottom is equal to zero, until a critical angle is reached and the fracture is activated.Regarding the damage patterns, the counter-clockwise torsion motions simulated for the insertion angles ⃗  = (−5, 0), ⃗  = (+5, +5), and ⃗  = (−5, −5) were characterized by the formation of asymmetric damage, where damage occurrence is more prominent on the side that experiences greater loading as the steps increase.Additionally, in terms of flexion and extension regimes, symmetrical patterns were observed for all three screw fixation angle combinations, as expected.These characterizations further validate the effectiveness of the present phase field modeling as a powerful tool for capturing damage in biological tissues.Moreover, they demonstrate good agreement with recent works found in the literature [25,29].

Staggered vs monolithic numerical schemes
A numerical comparison has been conducted to assess the performances of both the staggered and monolithic numerical schemes, as described in Secs.4.1 and 4.3, respectively, in terms of their accuracy in reproducing the fracture patterns and the Load vs. Displacement curves.
For this comparative analysis, the discretized vertebra-screws model consisting of 300.000 elements has considered, with the screws inserted at an angle combination of ⃗  = (−5, 0), in a counter-clockwise torsion motion.It is notable that the fracture pattern obtained from the monolithic method bears resemblance to the damage patterns observed in the staggered simulation, see Fig. 9. Furthermore, the fractioned volume outcomes are similar for both approaches.It can also be observed that, in terms of accuracy, the Load vs Displacement curves are more accurate in the monolithic scheme due to the fact that the binding reaction on the vertebra endplate is zero and deviates very little.Nevertheless, considering the computational time consumed, the advantages of using the staggered scheme outweigh the potential benefits offered by the monolithic technique.While the staggered simulation took approximately 3.5 hours to complete, the monolithic running time exceeded 15 hours.

Conclusions
A careful assessment of the fracture patterns and associated mechanisms in a human vertebra after the insertion of pedicle screws was scrutinized by implementing a finite element phase field method.Throughout the study, several aspects were pondered to accurately investigate the fractures and crack trajectories, such as maximum yield stress, critical energy release rate, and the characteristic length scale.
The analysis of various fracture types was ensured by conducting a mesh sensitivity study, which provided an optimal mesh size balancing simulation running time and the model's capability to reproduce outcomes within a small relative error.These findings were supported by comparing Force vs. Displacement curves for all considered mesh refinements at the extension vertebral motion mode and a screws insertion angle configuration of ⃗  = (−5, 0).It is worth mentioning that the damage responses in the vertebra from the developed phase field model are affected when different configurations of the pedicle screws fixation angle are simulated.
In addition, a comparison between phase field finite element approaches was also performed.Essentially, the staggered phase field model makes headway in terms of optimizing computational time consumption in relation to characterizing the damage within the vertebra-screws model, when compared to the results obtained from the monolithic phase field scheme.
(a) Identification of CC and ML directions of the pedicle screws virtually inserted.(b) Uniform computational discretized vertebra-screws meshes.From left to right: 100 Dof, 350 Dof, 700 Dof.
(a) Young's modulus (GPa) distribution of the vertebra.The red color represents the high stiffness layer given by the cortical bone, while the blue color indicates low stiffness layer given by the trabecular bone.(b)Young's modulus (GPa) distribution of the pedicle screws.The red color illustrates the high stiffness assumed for the pedicle screws.

Fig. 3 :
Fig. 3: Boundary conditions and loading regimes employed on the vertebra-screws models for implementation.
(a) Force vs. Displacement curves.(b) Relative error as function of the number of degrees of freedom.(c) Vertebral fractured volume as function of the loading steps.
(a) Back view of the evolution of the fracture pattern under vertical extension loading.(b) Back view of the evolution of the fracture pattern under vertical flexion loading regime.(c) Back view of the evolution of the fracture pattern under a counter-clockwise axial rotation loading regime applied on the head of the screws.

Fig. 7 :
Fig. 7: Evolution of the damage pattern caused by loading regimes applied on the pedicle screws head: (a) extension motion; (b) flexion motion; and (c) counter-clockwise axial rotation motion.
(a) Applied loading on a counter-clockwise axial rotation loading regime applied on the head of the screws as function of the loading step.(b) Fractured volume of the vertebra on a counterclockwise axial rotation loading regime applied on the head of the screws as the loading step increases.

Fig. 8 :
Fig. 8: Staggered and monolithic phase field models contrasted on a counter-clockwise axial rotation regime.

Fig. 9 :
Fig. 9: Damage pattern comparisons between staggered and monolithic phase field models on a counterclockwise axial rotation regime.