Development and Validation of Two Intact Lumbar Spine Finite Element Models for In Silico Investigations: Comparison of the Bone Modelling Approaches

: Finite element (FE) analyses contribute to a better understanding of the human lumbar spine’s biomechanics and serve as an effective predictive tool. This study aims to present the development of two L1–L5 FE models using literature-based (LBM) and patient-speciﬁc (PSM) bone material assignment approaches. The geometry of the lumbar spine was developed based on quantitative computed tomography scans. The LBM and the PSM were compared under pure and combined loads. Various biomechanical parameters were investigated to validate the models. The total range of motion of the LBM in pure ﬂexion-extension, lateral bending, and axial rotation were 30.9 ◦ , 29 ◦ , and 13.7 ◦ , respectively, while for the PSM, it was 31.6 ◦ , 28.6 ◦ , and 14.1 ◦ . The required computational time of the PSM to complete against pure and combined loads were 12.1 and 16.6 times higher on average compared to the LBM. This study demonstrated that both models agree with experimental and in silico results, although the cumulative distribution of the stress and characterization of strain values showed a noteworthy difference between the two models. Based on these ﬁndings, the clinically-focused biomechanical FE studies must perceive the differences in internal mechanical parameters and computational demand between the different bone modelling approaches.


Introduction
In recent years, finite element analysis (FEA) has been recognised as an effective predictive tool for assessing spinal biomechanics [1]. The application of FEA can eliminate the difficulties, limitations and ethical concerns associated with in vitro experiments, furthermore, FEA is more cost-effective than in vivo and in vitro experimental studies [2,3]. The possibility of employing complex load cases and boundary conditions allows the investigation of biomechanical parameters that are difficult to measure by experimental methods [4]. Recently, in silico studies and trials using finite element methods are becoming more widespread [5][6][7], enabling patient-specific investigations and the development of new treatments and medical devices [8][9][10].
Numerous validated lumbar spine finite element (FE) models were published in the literature, which have varying degree of patient-specificity depending on the employed modelling strategy [2,11,12]. A commonly used technique is the "hybrid patient-specific" modelling approach, where the geometry is based on the patient's medical imaging data

Geometry of the FE Models
Two finite element models of the intact lumbar spine (L1-L5) were developed using quantitative computed tomography (QCT) images (Hitachi Presto, Hitachi Medical Corporation, Tokyo, Japan) of a 24-year-old male patient without any lumbar spine-related musculoskeletal pathology. The study involving the human participant was reviewed and approved by the National Ethics Committee of Hungary and the National Institute of Pharmacy and Nutrition (reference number: OGYÉI/163-4/2019). Written informed consent was obtained from the individual for the publication of any potentially identifiable images or data included in this article. All methods were carried out in accordance with relevant guidelines and regulations.
The QCT (voxel size of 0.6 × 0.6 × 0.6 mm 3 ) scans were used to obtain the geometry of both FE models, with previously reported imaging protocol [22]. The QCT dataset (in DICOM extension) was acquired from the hospital's PACS and subsequently anonymised with Clinical Trial Processor software (Radiological Society of North America, https:// www.rsna.org/ctp.aspx, accessed on 10 August 2022) [23]. The 2D images of the QCT scans were imported into Materialise Mimics software (Mimics Research, Mimics Innovation Suite v23.0, Materialise, Leuven, Belgium) to create a patient-specific geometry model by Hounsfield Unit-based segmentation, setting the lower and upper threshold to 226 HU and 1712 HU, respectively. The surface mesh was then exported to 3-Matic software (Mimics Research, Mimics Innovation Suite v21.0, Materialise, Leuven, Belgium) in STL (stereolithography) format. The surface mesh was smoothed (smooth factor: 0.7, iteration no.: 6, shrinkage compensation: on) and remeshed with a target edge length of 1 mm to improve mesh quality and remove spikes and cavities on the vertebral surfaces [24] ( Figure 1a).
The current study applied two bone material assignment approaches to create the literature-based model (LBM) (Figure 2a) and the patient-specific model (PSM) (Figure 2b). In the LBM, the vertebral body consisted of a 1 mm thick cortical shell, trabecular core, 0.5 mm thick bony endplates, and posterior elements [25,26]. The volumes were meshed adaptively with 1 mm linear tetrahedral elements in HyperWorks software (Altair Engineering, Inc., Troy, Michigan, United States) based on the surface mesh [27]. Vertebrae of the PSM were not further divided, but rather treated as a whole and were meshed with a uniform element size of 1 mm with linear tetrahedral elements.
Soft tissues, such as the facet joints, were identical in both FE models and were modelled as a 0.25 mm thick cartilage layer with an initial gap of 0.5 mm between adjacent surfaces [28] and then meshed with 6-node wedge elements [29]. The intervertebral discs (IVD) consisted of the nucleus pulposus (NP), the annulus fibrosus (AF), including both the ground substance (GS), the fibres, and the 0.5 mm thick cartilaginous endplates (CEP). The AF GS was constructed as six concentric layers [14] around the NP, which accounted for 45% of the IVD [30] and was moved to the posterior direction in accordance with the literature and general anatomy [31]. The AF GS and the NP were represented with 8-node hybrid hexahedral elements [14], while CEP was meshed with linear tetrahedral and pyramid elements (Figure 1b). The fibres were modelled with tension-only truss elements in a criss-cross pattern with an alternating angle of ±30 • to the axial plane of the IVD [32]. The cross-sectional area values of the fibre layers were calculated from the volume ratio of the fibres to the ground substance: 23% (outermost), 20%, 17%, 13%, 9%, and 5% (innermost), respectively [33] (Figure 1b). All seven major spinal ligaments were modelled: anterior longitudinal ligament, posterior longitudinal ligament, ligamentum flavum, interspinous ligament, supraspinous ligament, intertransverse ligament, and capsular ligaments ( Figure 1c).

Material Properties
For the LBM, isotropic, homogeneous, and linear elastic material properties from the literature were assigned to the bony elements [16,34,35] (Figure 2a). The moduli of elasticity and the Poisson's ratios of the bony elements for the LBM are presented in Table 1.
In the PSM, locally defined material properties were applied to model bone heterogeneity, as the mechanical properties of the bone tissues were assigned to each element based on the QCT scans. Hydroxyapatite calibration phantoms (Hitachi Presto, Hitachi Medical Corporation, Tokyo, Japan) were used to determine the relationship between the Hounsfield Units (HU) and elastic moduli.

Material Properties
For the LBM, isotropic, homogeneous, and linear elastic material properties from the literature were assigned to the bony elements [16,34,35] (Figure 2a). The moduli of elas ticity and the Poisson's ratios of the bony elements for the LBM are presented in Table 1.
In the PSM, locally defined material properties were applied to model bone hetero geneity, as the mechanical properties of the bone tissues were assigned to each elemen based on the QCT scans. Hydroxyapatite calibration phantoms (Hitachi Presto, Hitach Medical Corporation, Tokyo, Japan) were used to determine the relationship between the Hounsfield Units (HU) and elastic moduli.
To assign the calculated elastic moduli to the elements, the volume meshes were im ported into Materialise Mimics software. Elements were placed in ten equal-sized groups based on the average HU value of the voxels contained by an element [41]. Elements with HU < 0 were placed in an eleventh group, representing fat-and marrow-like materials with lower density [42]. The median HU of each group was set as the representative value First, the HU was converted to equivalent bone mineral density (ρ QCT ) by linear calibration using calibration phantoms (1) with known densities of 50, 100, 150, and 200 mg/cm 3 [36]. Equivalent bone mineral density was considered to be equal to ash density (ρ ash ) (2) [11,37,38], while the apparent density (ρ app ) was calculated as 60% of ash density (3) [11,39]. The relationship between apparent density and HU (4) was calculated employing the Equations (1)-(3). The isotropic elastic modulus (E) was calculated from apparent density based on a relationship from the literature [40].
To assign the calculated elastic moduli to the elements, the volume meshes were imported into Materialise Mimics software. Elements were placed in ten equal-sized groups based on the average HU value of the voxels contained by an element [41]. Elements with HU < 0 were placed in an eleventh group, representing fat-and marrow-like materials with lower density [42]. The median HU of each group was set as the representative value and was assigned to all elements in the group. The Poisson's ratio of all bony tissues in the PSM was set to be 0.3 [37][38][39] (Figure 2b).
The material properties of the soft tissues were identical in the two models. The cartilaginous endplates were assigned with linear elastic material property taken from the literature [27]. NP and AF ground substance was modelled with an incompressible 2-parameter hyperelastic Mooney-Rivlin formulation. The mechanical behaviour of the annulus' collagenous fibres was represented by a nonlinear stress-strain curve from Shirazi-Adl [35]. Since the outer layers of the IVD are stiffer than the inner layers, the collagen fibres were weighted by scalar factors defined for each layer ( Figure S1) [43]. The facet cartilage was modelled with a hyperelastic Neo-Hooke formulation [27]. Tension-only, nonlinear uniaxial spring elements [44] were used to simulate the ligament behaviour (Appendix A). The material properties used in the current study are summarised in Table 1.

Loading and Boundary Conditions
To accurately validate the two finite element models, three different loading cases were applied in this study. First, a pure bending moment of 7.5 Nm was applied in the three anatomical planes at the cranial endplate of L1 to simulate flexion, extension, lateral bending, and axial rotation (Figure 2a,b) [48]. Second, to measure the IDP, a compressive follower load of 1000 N [49] was employed without any bending moments. The load was applied along a pre-optimised path through the centre of rotations of the vertebral bodies to minimise the presence of artefact bending moments (Figure 2a,b) [50]. The follower load technique considers both the effect of the upper body weight and the stabilising effect of the local muscle forces in compressive loading modes such as erect standing. Third, a combination of a compressive preload and a bending moment was used. The magnitude of the preload and moment for flexion, extension, lateral bending, and axial rotation were adopted from Dreischarf et al.: 1175 N with 7.5 Nm, 500 N with 7.5 Nm, 700 N with 7.8 Nm, and 720 N with 5.5 Nm, respectively [2]. During the simulations, the most caudal endplate of L5 was fixed in all degrees of freedom. The contact behaviour between the facet joint surfaces was assumed to be frictionless [51].

Mesh Convergence
To ensure that the mesh resolution did not affect the accuracy of model predictions, all elements of the bony parts of L4 vertebra was investigated with five different element lengths (0.607, 1, 1.65, 2.73, and 4.505 mm) for both modelling methods [11]. The edge length of the most refined mesh was determined based on the resolution of the CT. The mesh was considered convergent when the difference in the von Mises stress values was less than 10% compared to the most refined mesh [11,51]. As the results are most sensitive to axial rotation, a 7.5 Nm was applied on the upper vertebral endplate, while the lower endplate was fixed at all degrees of freedom [51,52].

Validation
The fully configured PSM and LBM with assigned material properties and applied loading and boundary conditions were exported from the pre-processor software Hyper-Works to Abaqus Standard (v2021, Dassault Systemes, Vélizy-Villacoublay, France) to solve the finite element simulations.
Both the literature-based and patient-specific FE models were validated. Outcomes were compared with in vivo [53][54][55][56] and in vitro experimental measurements taken from the literature [57][58][59][60]. Dreischarf et al. have shown that comparison with diverse finite element results can increase the quality and reliability of the FE model's prediction [2]. Therefore, the results of the current study were compared to multiple previously published well-accepted in silico data [2,29].
In the case of pure bending, the ROM, the intervertebral rotations (IVR), and the facet joint forces (FJF) were compared with in vitro [57,58,60] and in silico data [2,29] at each spinal level. In addition, the load-deflection curves of L1-L5 were plotted together with an in vitro measurement, and a range from numerical results. The ROM and IVR were calculated with an in-house developed MATLAB (MathWorks, Natick, Massachusetts, United States) script (Appendix B). Against pure compressive follower load, the intradiscal pressure (IDP) in the L4-L5 nucleus pulposus of both models was averaged over all elements and then compared with an in vitro [59] and a range of in silico results [2]. A combined compressive, and bending load was applied to the FE models to measure the IVR, FJF, and IDP values at each spinal level and compared with in silico [2] and in vivo measurements [53][54][55][56], as the combined loads are assumed to be the most representative for in vivo experiments [34].

Stress Distribution
In addition to the validation process, internal mechanical parameters were analysed to investigate the differences originating from the different modelling strategies applied. The Empirical Cumulative Distribution Functions (ECDF) of the von Mises stress values of all nodes in the vertebrae's bony parts were plotted as a function of the load and material models applied ( Figure S2). The ECDF is a step function that jumps up by 1/n at each n data point, where n is the number of von Mises stress values, practically equal to the number of nodes in the samples. At any specified von Mises stress value, ECDF gives the fraction of stress observations that are less or equal to that specified von Mises value [61]. Using MATLAB software, the ECDF estimated the cumulative distribution of the von Mises stress values in pure bending cases and both material models. The first and fifth lumbar vertebrae were excluded from this analysis to avoid data distortions due to load transfer and fixation.
ECDFs can describe the cumulative distribution of all nodes in the analysed sample but cannot characterise the maximum loaded state. To describe how the different bone material assignments affect the most compressed state of the vertebrae, 1% of the nodes with the smallest minimum principal (equivalent to highest compressive) strain from both models were analysed. Since PSM contained about tenfold more nodes than the LBM, a percentage-based selection was implemented instead of using an absolute number of nodes. In pure bending, all nodes of the bony regions were included in the sample, but L1 and L5 vertebrae were excluded to avoid compromised strain values due to load and boundary conditions applied to those regions. Boxplot figures with median strain values, whiskers, and outliers were plotted as a function of the loading cases for both LBM and PSM.

Mesh Convergence
The mesh convergence test analysed the percentage differences in the von Mises stress of the L4 vertebra between the finest and coarser meshes. The differences in stress induced by the coarser meshes compared to the finest mesh were 28.76%, 14.96%, 10.02%, and 9.07%, in order of decreasing edge length for the LBM and 46.98%, 30.57%, 16.99%, and 2.73% for the PSM, respectively ( Figure S3). For both the LBM and the PSM, the 1 mm mesh was the first to yield a difference of less than 10%, therefore found to converge for stress and was applied for both models.

Computational Times
The number of elements and the required computational time were compared. Due to the different FE meshing algorithms (LBM: adaptive meshing, PSM: uniform meshing), the PSM contained 12 million elements and the LBM 1.6 million elements. This difference resulted that, on average, the PSM took 12.1 times more to complete than the LBM against pure bending and compressive follower load. A similar trend was observed in the case of the combined load cases, namely that, on average, the PSM took 16.6 times more to finish the simulations ( Figure S4).  [58], the predictions of the current study were slightly out of range in extension and close to the mid-value in axial rotation. The agreement with the in silico literature data was excellent for all load cases [2] (Figure 3c).

Validation of the Lumbar Spine Models
IVR values occurred within the range of experimental measurements [60], except at the L2-3 and L4-5 levels in flexion, where the models slightly underestimated the in vitro measurement, although our results correlate well with the in silico data [29]. In extension, the IVR values of the PSM were 4.57%, 6.14%, 8.65%, and 6.32% greater at the L1-2, L2-3, L3-4, and L4-5 spinal levels compared to the LBM (Figure 4a-d).

Pure Compression Load
The IDP values under the compressive follower load were 0.99 and 1.02 MPa for the LBM and PSM, respectively. The results were in good agreement with both the experimental results [59] and the in silico range [2] (Figure 3d).

Combined Load
Against combined load, the IVR results of the LBM and PSM were consistent with the available literature data, although the results in flexion at L2-3, L3-4 and L4-5 were outside the in vivo range [53][54][55]. The in silico range [2] was not reached in flexion at L1-2, in extension at L2-3 and L4-5, and in axial rotation at L1-2, L2-3 and L3-4. Compared to the LBM, the PSM provided greater IVR values in extension at all spinal levels (Figure 4e-g).
The IDP values showed good agreement with the published mean FE mid-values, while the simulation results approximated the available in vivo measurements [56] adequately. The PSM induced higher intradiscal pressure values in all loading conditions at all spinal levels ( Figure 5).
The predicted facet joint forces were calculated at all spinal levels and compared to available in silico data [2]. The FJF values were not adaptable in flexion due to the lack of contact between the adjacent facet joint surfaces [2,51]. The outcomes of the current study were consistent with the existing literature data for all the applicable load cases. In agreement with other FE results, no contact occurred at L1-2, L2-3 and L3-4 in lateral bending [2]. The predicted FJF values were smaller at L1-2 and L2-3 in axial rotation than the FE ranges from the literature ( Figure 6).

Stress Distributions
The ECDF plots of the estimated cumulative distribution showed noteworthy differences between the models developed by two bone material assignment strategies (Figure 7a-d). Generally, the cumulative distribution of stress values in the low-stress range, from 0 MPa to 1 MPa, shows the same characteristics. It is also observed that the stress values in the LBM are more spread out and have a less steep curve than in the PSM. This difference in the curve steepness also means that proportionally, the LBM contains more high-stress elements.  [60]. The grey bars and their range correspond to the median and range of results of a calibrated FE simulation [29]. (e) IVR results of the LBM and PSM at different spinal levels against combined loads in flexion, (f) extension, (g) lateral bending, and (h) axial rotation. The white bars and their range correspond to the median, the minimum, and maximum values of multiple in vivo measurements [53][54][55]. The grey bars and their range correspond to the median and range of multiple validated FE models result [2].

Pure Compression Load
The IDP values under the compressive follower load were 0.99 and 1.02 MPa for the LBM and PSM, respectively. The results were in good agreement with both the experimental results [59] and the in silico range [2] (Figure 3d).  [60]. The grey bars and their range correspond to the median and range of results of a calibrated FE simulation [29]. (e) IVR results of the LBM and PSM at different spinal levels against combined loads in flexion, (f) extension, (g) lateral bending, and (h) axial rotation. The white bars and their range correspond to the median, the minimum, and maximum values of multiple in vivo measurements [53][54][55]. The grey bars and their range correspond to the median and range of multiple validated FE models result [2].
The strain values characterise the most compressed state of the bony tissues, as nodes with the smallest minimum principal (equivalent to highest compressive) strain were included in this analysis. In the LBM, the median compressive strain values were calculated as −0.0053, −0.0040, −0.0059, and −0.0023 for flexion, extension, lateral bending, and axial rotation, respectively. In contrast to the LBM, the PSM decreased median compressive strain values by 77.3%, 62.5%, 79.7%, and 39.1%, respectively (Figure 8). Although the outlier values reveals that the PSM spreads the strain values more evenly, as the box and whiskers are narrow, and contains more outliers than the LBM. While the median values are lower in the LBM, the outliers are even lower in the PSM, indicating that the PSM includes more nodes with higher compression.

Combined Load
Against combined load, the IVR results of the LBM and PSM were consistent with the available literature data, although the results in flexion at L2-3, L3-4 and L4-5 were outside the in vivo range [53][54][55]. The in silico range [2] was not reached in flexion at L1-2, in extension at L2-3 and L4-5, and in axial rotation at L1-2, L2-3 and L3-4. Compared to the LBM, the PSM provided greater IVR values in extension at all spinal levels ( Figure  4e-g).
The IDP values showed good agreement with the published mean FE mid-values, while the simulation results approximated the available in vivo measurements [56] adequately. The PSM induced higher intradiscal pressure values in all loading conditions at all spinal levels ( Figure 5). The predicted facet joint forces were calculated at all spinal levels and compared to available in silico data [2]. The FJF values were not adaptable in flexion due to the lack of contact between the adjacent facet joint surfaces [2,51]. The outcomes of the current study were consistent with the existing literature data for all the applicable load cases. In agreement with other FE results, no contact occurred at L1-2, L2-3 and L3-4 in lateral bending [2]. The predicted FJF values were smaller at L1-2 and L2-3 in axial rotation than the FE ranges from the literature ( Figure 6).

Stress Distributions
The ECDF plots of the estimated cumulative distribution showed noteworthy diffe ences between the models developed by two bone material assignment strategies (Figur  7a-d). Generally, the cumulative distribution of stress values in the low-stress range, from 0 MPa to 1 MPa, shows the same characteristics. It is also observed that the stress value The ECDF plots of the estimated cumulative distribution showed noteworthy differences between the models developed by two bone material assignment strategies ( Figure  7a-d). Generally, the cumulative distribution of stress values in the low-stress range, from 0 MPa to 1 MPa, shows the same characteristics. It is also observed that the stress values in the LBM are more spread out and have a less steep curve than in the PSM. This difference in the curve steepness also means that proportionally, the LBM contains more highstress elements.

Discussion
The current study presents and compares two different bone material assign strategies commonly used in FE-based biomechanical investigations. Verification developed FE meshes was performed through a mesh convergence study to ensure p simulation predictions. Then, the developed models were validated using numerou mechanical parameters, such as ROM, IVR, IDP, and FJF. Furthermore, the required putational time, cumulative distribution of the stress values and the compressive values in the bony elements of the L2, L3, and L4 vertebra were assessed. The specifi challenges, and limitations of these techniques were highlighted.

Discussion
The current study presents and compares two different bone material assignment strategies commonly used in FE-based biomechanical investigations. Verification of the developed FE meshes was performed through a mesh convergence study to ensure proper simulation predictions. Then, the developed models were validated using numerous biomechanical parameters, such as ROM, IVR, IDP, and FJF. Furthermore, the required computational time, cumulative distribution of the stress values and the compressive strain values in the bony elements of the L2, L3, and L4 vertebra were assessed. The specificities, challenges, and limitations of these techniques were highlighted.

Locally Defined Material Properties
Clinical application of bone imaging for the in vivo assessment of structural properties involves, among others, dual X-ray absorptiometry (DXA), quantitative computed tomography (QCT), quantitative ultrasonometry (QUS), and high-resolution magnetic resonance (HR-MR) [62]. These non-invasive and non-destructive techniques are coupled with different advantages and disadvantages. DXA and QCT need ionising radiation, while HR-MR and QUS apply no radiation load on the patient [63]. Using the QUS technique in vivo, the number of measurement sites is limited; there is also little opportunity for calibration standardisation [63,64]. In contrast to QUS, the QCT technique can provide accurate data on the 3D BMD of the trabecular and cortical compartments via calibration. HR-MR acquisition takes much more time and is technically demanding, although, in contrast to DXA, the apparent microstructure of the bone is assessable, likewise with QCT [62,65].
In addition to these techniques, finite element models employing locally defined material properties have been used to evaluate fracture risk and ultimate strength and have shown superior results in comparison with densitometry measurements alone [62,65]. Due to the relative availability of QCT scans and the possibility of accurately modelling locally defined material properties in FE models, the current study employed QCT-derived materials in the patient-specific model.
Multiple studies have defined empirical relationships for assigning local material properties based on QCT for different anatomic sites by comparing the QCT-based bone density with the results of mechanical testing [40,[66][67][68]. According to Carpenter et al. [69], the QCT-based bone density can be directly converted into local mechanical properties when using calibration phantoms, which approach has been used in several finite element studies in various anatomical areas [19,37,[70][71][72]. The finite element models typically used voxel-based hexahedral elements [37,[72][73][74] or tetrahedral elements that follow the smoothed geometry [19,36,70,71,75]. Density assignment is more directly implemented in voxel-based models, where voxels are directly converted into hexahedral elements, and thus the density of an element is equal to the density of the corresponding voxel of the QCT scan [37]. On the contrary, tetrahedral-based models better approximate the complex structure and geometry of vertebrae than hexahedral ones [76]. However, as Imai et al. [77] stated, in these models, the density of an element is the average density of the voxels contained in one element, which clearly shows why QCT resolution and mesh resolution is important [11]. Furthermore, Morgan et al. [40] concluded that modulusdensity relationships depend on the anatomic site, which means there is no universal modulus-density relationship, so it is essential that the used relationship comes from the same anatomical site as the one being investigated. Accordingly, in our research, besides applying a thin-sliced QCT imaging protocol [22] and accurate tetrahedral mesh resolution [11], we implemented a commonly used, well-established correlation that has been determined based on in vitro measurements of vertebrae [40].

Mesh Convergence
Although in the literature, the displacement [78] and strain energy [52,76] are also used as parameters of mesh convergence, the von Mises stress was chosen since mesh resolution is more sensitive to internal mechanical parameters than displacement [52]. Based on the results, the PSM is more sensitive to the mesh resolution in terms of stress due to its element-based material assignment strategy ( Figure S3). In contrast, in the LBM, relatively minor stress changes occurred due to the homogeneous modelling of different anatomical sections. Both the LBM and PSM gave results less than 10% at 1 mm element size, agreeing with the findings of Costa et al. [11], who applied an identical imaging protocol and similar bone material assignment strategy. The required computational time of the models with the finest mesh increased significantly compared to the 1 mm results, which also reinforced the choice of the 1 mm element since the elevated computational demand does not couple with higher accuracy ( Figure S3).

Computational Times
In general, applying the literature-based modelling approach reduced the required CPU time since the PSM notably needed more time to complete for both the pure and combined loads. In the case of the patient-specific bone material assignment, the computational cost dramatically increased due to the different meshing algorithms applied for the two models. However, it is worth mentioning that the patient-specific model required less time to develop than the literature-based model considering the soft tissues are the same for both models. Since segmentation is required for both models, the difference is due to the design of the bony parts, as the calibration between Hounsfield Units and densities is based on segmented calibration phantoms that take less time than the separation of the cortical part, the trabecular part, the posterior elements, and the bony endplates for each vertebra.

Model Validation
Results of the current study were compared with available in vitro and in silico data from the literature. Dreischarf et al. described that the validation via diverse results increase the accuracy and reliability of the developed models' predictions [2]. Therefore, validation was performed in the current study utilising several published experimental and in silico results. The LBM and PSM yielded similar ROM results against pure bending moment of 7.5 Nm indicating that the material assignment of the bony structures slightly influences the spinal kinematics of the spine. Similar findings were described by Vena et al., stating that for FE simulations, the material properties of the bony structures have a minor influence on the spinal kinematics; thus, vertebral bodies could be modelled as rigid bodies [79]. Moreover, Remus et al. presented, calibrated, and validated a hybrid finite element-multibody model with rigid vertebral bodies, indicating that the hybrid simulation model is powerful and efficient for providing valid mechanical responses. Their study confirmed that the soft tissues such as the spinal ligaments and the IVD have primarily influence the kinematic properties of the spine [29]. Thus, nonlinear modelling of such structures is crucial, and therefore in the current study, the spinal ligaments and components of the intervertebral discs were modelled with nonlinear materials. The two models were in a good agreement with the published in silico and in vitro studies in L1-L5 ROM of flexion-extension, lateral bending, and axial rotation. However, the nonlinear load-displacement curve (Figure 3b) demonstrates that in flexion, both the patient-specific and literature-based models are slightly outside the in vitro range, while in extension, the LBM is within, the PSM is slightly outside the range. The same trend is reported by Dreischarf et al. [2], that majority of their assessed FE models overestimate the motion in extension and underestimate it in flexion, moreover, good agreement was found with the reported in silico results for lateral bending and axial rotation as well [2].
The FJF depends on several patient-specific factors, such as the underlying geometry of the individual cartilage surfaces, the gap size between the adjacent surfaces, or the thickness of the cartilage layer [80]. There is limited data available in the literature on the in vitro measurement of FJF, and these studies have some limitations and uncertainties [51]. Nevertheless, in pure bending, the results of the current study show good agreement with the available in vitro [58] and in silico [2] data, although a slight overestimation of the in vitro result appeared in extension (Figure 3c). The current study results indicate that the FJF in lateral bending is notably lower than in extension and axial rotation, which may originate from the patient-specific anatomy of the current study's FE model agreeing with the findings of Woldtvedt et al. [80]. In flexion, no force was transmitted due to the distancing cartilage surfaces during the load [34].
The IDP predictions of the FE models in pure compressive follower load showed excellent agreement with the median in vitro results at 300 N and 1000 N (Figure 3d). Among the published FE models in the literature, the NP is often modelled as incompressible fluid-filled cavity [25,49], which allows simulating the effect of the initial hydrostatic pressure in the IVD [81]. Nonetheless, as Dreischarf et al. pointed out, most NPs modelled with fluid-filled cavities neglect this effect since IDP values measured in the unloaded state are less significant [2].
The IVR values were obtained and evaluated for both pure and combined load cases (Figure 4a-h). Wilke et al. previously reported that 7.5 Nm pure bending moment is the most appropriate for simulating the loads on the lumbar section, including effects of muscle forces [48]. Panjabi et al. performed a human cadaveric test [60], while Remus et al. presented a hybrid FE-multibody analysis to investigate the biomechanical properties of the human lumbar spine [29]. Accordingly, the IVR results of these in vitro and in silico studies have been used as a basis for comparison in the case of pure bending load. The employed combined loads were determined by Rohlmann et al. and Dreischarf et al., using follower load and bending moment combination that most closely approximates the average of in vivo measurements reported in the literature [34,82,83]. In the case of the combined load, the average of three in vivo measurements [53][54][55] was used as a basis for comparison when evaluating the IVR results. In both load cases, the FE predictions of the current study agree with the experimental results in extension, lateral bending, and axial rotation, while the simulations underestimated the experimental results in flexion.
The predictions of the current study are comparable to the in silico findings reported by Remus et al. and Dreischarf et al. [2,29].
The predicted IDP values for combined load were compared to the mean of eight well-accepted in silico results [2]. Moreover, Wilke et al. performed an in vivo analysis in a non-degenerative L4-L5 IVD [56]. Their study was designed to provide a validation dataset for in silico or in vitro investigations. Excellent agreement with the in vivo and in silico results indicates that both models can predict the intradiscal pressure values accurately.
In the case of combined load, the FJF values were compared only to those mean in silico results reported by Dreischarf et al. due to the lack of relevant experimental results [2]. In general, the validated in silico studies show a significant variance regarding the FJF's magnitude, indicating that this force is sensitive to the FE modelling techniques and the patient-specific anatomy. This finding is also supported by Schmidt et al., who found that the FJF is highly dependent on the individual geometry of facet joint surfaces [84]. In flexion, no force occurred due to the lack of contact. In extension, the magnitude of the FJF increases caudally, agreeing to the results reported by Remus et al. [29]. In lateral bending, FJF occurred only at the L4-L5 spinal level, as it was reported in several models investigated by Dreischarf et al. [2].

Stress Distributions
The plot of the cumulative distribution of stress values highlights the fundamental difference between the LBM and PSM. The differences in the high-stress (above 1 MPa) range's distribution indicate that the application of patient-specific bone materials yields substantially different internal mechanical properties. The difference in the slope between the distribution curves is the largest for axial rotation, which agrees with the findings of Ayturk. et al. [52], i.e., the stress distribution of finite element models is most sensitive to axial rotation. The analysis of the minimum principal strain values also concludes that the two models show significant differences in internal mechanical parameters under the same loading conditions. In the LBM model, the strain values cluster around a peak with only a few outliers, while in the PSM, no such strain peak is observed. The clustering effect explains why the LBM gave lower median values than the PSM; since investigating the outliers, one can see that the PSM contains bony tissues with higher compression. Based on these findings, using medical imaging-based bone material properties could be beneficial in biomechanical investigations where it is essential to analyze internal mechanical parameters, i.e., vertebral fracture risk or vertebral strength assessment [11,17,18]. These conclusions confirm the findings of Mosleh et al., who found a good correlation between the fracture loads estimated by QCT-based FE models and in-vitro fracture loads [85]. Moreover, a biomechanical study by Rezaei has shown how QCT combined FE simulations can estimate fracture outcomes in single vertebral models [86]. Analogously, literature-based linearly elastic homogeneous material models can only be used carefully in cases where bone material quality significantly influences the parameters to be investigated.

Limitations, Significance
Like most finite element analyses, the current study also has limitations and simplifications worth highlighting. By its very nature, the FE method cannot address biological or physiological phenomena, as it solely works on a mechanical principle. The viscoelastic and poroelastic characteristic of the IVD was not considered since a quasi-static analysis was performed, and different time-dependent behaviours, such as creep, were not the interest of the current study. Additionally, the presented models have the potential for further development, with considering the effect of muscle forces, increasing the number of material groups in patient-specific bone material assignment, or modelling the NP as an incompressible fluid cavity and include initial hydrostatic pressure field. Although the models could be improved by addressing these simplifications, we believe that the relationship between the models and their results would not change; thus, the presented and validated models are suitable for drawing conclusions and for further biomechanical investigations.

Conclusions
The current study aimed to present two commonly used bone material assigning techniques to develop a finite element model of the healthy human lumbar spine. Mesh resolution was investigated to ensure accurate model predictions. The results of the FE models were compared with in vivo, in vitro, and well-accepted in silico data from the literature to validate them. Validation was performed considering ROM, IVR, IDP, and FJF variables against pure, and combined loads. Based on the parameters and loads investigated, both the LBM and PSM can be used to model the biomechanical properties, such as ROM, IVR, FJF, and IDP of the lumbar spine, as they are in good agreement with the previously published results in most investigated variables. However, substantial differences in computational time and cumulative distribution of von Mises stress values were observed between the literature-based and the patient-specific models.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app122010256/s1, Figure S1: The nonlinear stress-strain relationships of the six different annulus fibre layers in the intervertebral disc. Layer 1 is the outermost and layer 6 is the innermost layer; Figure  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The extracted data and the applied code of this investigation are appended as Supplementary Material to this article. In addition to the code, its detailed description, figures and file structure are shared for the sake of clarity and reproducibility. The authors believe that the "Materials and Methods" section, the appendices, and the developed in-house code cover the investigated topics satisfactorily.

Acknowledgments:
The support of CAD-Terv Kft. in using Abaqus, especially of István Nadj and Gábor Brezvai, is highly acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The attachment points and orientation of the elements modelling the ligaments were adopted from a previously published study, MySPINE [22]. All seven major spinal ligaments were modelled: anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), ligamentum flavum (LF), interspinous ligament (ISL), supraspinous ligament (SSL), intertransverse ligament (ITL), and capsular ligaments (CL). In the current study the ligaments were modelled with multiple spring elements in either parallel-series connection (ALL, PLL) or in parallel connection only (ISL, CL, ITL, LF, PLL, SSL) (Table A1) [44].  [47] (Table A2). As presented, Rohlmann et al. [47] defined the spring-like behaviour of the ligaments with strain dependent stiffness value of a single spring, describing the progressive characteristic. Nonlinear spring characteristic should be defined in Abaqus with force and elongation value pairs assigned to the corresponding SPRINGA element sets. To assign the material properties correctly, the following calculations were applied: For the sake of simplicity, 1 was chosen as an upper limit of the ligament strains, since 100% strain is not reached during the simulation, thus the calculation remains stable. From this, the elementary force can be calculated (the force which is assigned to each element): Beside the force values Abaqus requires the corresponding elongation values (strain dependent and elementary) as well: where: F-force (N) p-number of parallel chains (-) s-number of elements in a chain (-) K n -stiffness (N/mm) ε n -strain (-) D-elongation (mm) L 0 -average length of one element (mm) The applied material properties defined by Rohlmann [47] were plotted in Figure A1, where the progressive behaviour of the ligaments is visualised. For the sake of simplicity, 1 was chosen as an upper limit of the ligament strains, since 100% strain is not reached during the simulation, thus the calculation remains stable. From this, the elementary force can be calculated (the force which is assigned to each element): Beside the force values Abaqus requires the corresponding elongation values (strain dependent and elementary) as well: where: -force (N) -number of parallel chains (-) -number of elements in a chain (-) -stiffness (N/mm) -strain (-) -elongation (mm) -average length of one element (mm) The applied material properties defined by Rohlmann [47] were plotted in Figure A1, where the progressive behaviour of the ligaments is visualised. Figure A1. The nonlinear force-strain relationships of the seven spinal ligaments applied in the current study. Figure A1. The nonlinear force-strain relationships of the seven spinal ligaments applied in the current study.

Appendix B
Calculation of the investigated variables from the finite element analysis is challenging. The current study considered two force, two pressure, and four rotation variables. The force (FJF, FJF_SEG) and pressure (IDP, IDP_SEG) variables can be calculated using average and maximum calculations after summary of the exported data. The calculation of the rotation variables (IVR, IVR_COMB, ROM, ROM_MAX) was based on virtual markers. The markers were nodes coupled to the vertebrae with rigid coupling. Three markers were used for each vertebra: middle points of the upper and lower endplate and a point in the middle sagittal plane of the spinous process ( Figure A2). Appl. Sci. 2022, 12, x FOR PEER REVIEW 20 of 26

Appendix B
Calculation of the investigated variables from the finite element analysis is challenging. The current study considered two force, two pressure, and four rotation variables. The force (FJF, FJF_SEG) and pressure (IDP, IDP_SEG) variables can be calculated using average and maximum calculations after summary of the exported data. The calculation of the rotation variables (IVR, IVR_COMB, ROM, ROM_MAX) was based on virtual markers. The markers were nodes coupled to the vertebrae with rigid coupling. Three markers were used for each vertebra: middle points of the upper and lower endplate and a point in the middle sagittal plane of the spinous process ( Figure A2). For each vertebra, the middle points of the endplates define an axis and the three markers define a plane, so a local coordinate system can be created according to Panjabi and White [30] ( Figure A3). The definition of the axes means the sign of the rotations is as follows: positive: flexion, left rotation, right bending; negative: extension, right rotation, left bending.  For each vertebra, the middle points of the endplates define an axis and the three markers define a plane, so a local coordinate system can be created according to Panjabi and White [30] ( Figure A3). The definition of the axes means the sign of the rotations is as follows: positive: flexion, left rotation, right bending; negative: extension, right rotation, left bending.

Appendix B
Calculation of the investigated variables from the finite element analysis is challenging. The current study considered two force, two pressure, and four rotation variables. The force (FJF, FJF_SEG) and pressure (IDP, IDP_SEG) variables can be calculated using average and maximum calculations after summary of the exported data. The calculation of the rotation variables (IVR, IVR_COMB, ROM, ROM_MAX) was based on virtual markers. The markers were nodes coupled to the vertebrae with rigid coupling. Three markers were used for each vertebra: middle points of the upper and lower endplate and a point in the middle sagittal plane of the spinous process ( Figure A2). For each vertebra, the middle points of the endplates define an axis and the three markers define a plane, so a local coordinate system can be created according to Panjabi and White [30] ( Figure A3). The definition of the axes means the sign of the rotations is as follows: positive: flexion, left rotation, right bending; negative: extension, right rotation, left bending.  The rotation variables were calculated based on these local coordinate systems (LCSs) of the vertebrae: the rotation of a segment (IVR, IVR_COMB) was calculated as the rotation of the upper vertebra's LCS expressed in the lower vertebra's LCS. The rotation of the whole lumbar spine (ROM, ROM_MAX) was calculated likewise: the rotation of the top vertebra's LCS expressed in the bottom vertebra's LCS. To calculate these rotations, the calculation of a rotation matrix for each LCS is necessary since the coordinates of the markers are given in the global coordinate system.
Code S1 contains a folder structure ( Figure A4) and a MATLAB code structure ( Figure A5) which help with the calculation and visualisation of the above-mentioned variables. The user should load data to the right places in the RAW_DATA folder based on its type and run INTACT_LUMBAR.m in the MATLAB folder. When the code finished, the figures can be found in the PICTURES folder. Appl. Sci. 2022, 12, x FOR PEER REVIEW 21 The rotation variables were calculated based on these local coordinate systems (L of the vertebrae: the rotation of a segment (IVR, IVR_COMB) was calculated as the tion of the upper vertebra's LCS expressed in the lower vertebra's LCS. The rotation o whole lumbar spine (ROM, ROM_MAX) was calculated likewise: the rotation of th vertebra's LCS expressed in the bottom vertebra's LCS. To calculate these rotation calculation of a rotation matrix for each LCS is necessary since the coordinates of the m ers are given in the global coordinate system.
Code S1 contains a folder structure ( Figure A4) and a MATLAB code structure ure A5) which help with the calculation and visualisation of the above-mentioned v bles. The user should load data to the right places in the RAW_DATA folder based o type and run INTACT_LUMBAR.m in the MATLAB folder. When the code finished figures can be found in the PICTURES folder. The purpose and content of the folders: