We utilized a CT image-derived L4-L5 geometry that was manually reconstructed. This geometry was employed in a previously published study by Nicolini et al. (2022b). We scaled this geometry using average dimensions (Figure 1) from various experiments in the literature (Shirazi-Adl, 1994; Schmidt et al., 2006; Little et al., 2008; Kiapour et al., 2009; Zander et al., 2009; Busscher et al., 2010; Ayturk and Puttlitz, 2011; Liu et al., 2011; Park et al., 2013; Jaramillo et al., 2015; Bashkuev et al., 2020; Nicolini et al., 2022a). The geometry was divided into AF and NP, with the NP accounting for 44% of the total disc volume, as in established models (Schmidt et al., 2006; Ezquerro et al., 2011; Lavecchia et al., 2018; Wang et al., 2021a; Nicolini et al., 2022a). The center of the NP was positioned posterior to the geometric center of the disc (Noailly et al., 2007; Weisse et al., 2012). For detailed representation, the annulus was divided into five distinct subregions (anterior (A), anterior-lateral (B), lateral (C), posterior-lateral (D), and posterior (E)) (Figure 1), following the measurements of Holzapfel et al. (2005). Furthermore, the annulus geometry was partitioned into five layers, according to the methodology of Nicolini et al. (2022a).

A comprehensive description of the model implementation utilizing an anisotropic formulation for the AF is presented in a prior publication (Nicolini et al., 2022a). The NP was defined using a Mooney-Rivlin material model (Schmidt et al., 2006; Schmidt et al., 2007; Ezquerro et al., 2011; Tang and Rebholz, 2011; Yang and O’Connell, 2017; Cai et al., 2020; Wang et al., 2021a). This hyperelastic material model is derived from the strain energy density function as described in Eq. 1 (Freutel et al., 2014; Nicolini et al., 2022b): W = C 10 I 1 − 3 + C 01 I 2 − 3 + D J α − 1 2 . (1) The stiffness and compressibility of the matrix are determined by the material parameters C ₁₀, C ₀₁, and D. J _α defines the elastic volume strain, and I ₁ and I ₂ represent the first and second invariants of the right Cauchy-Green deformation tensor. To realize the anisotropic formulation for the AF, the HGO material definition was applied (Holzapfel et al., 2000; Eberlein et al., 2001). This model integrates an isotropic matrix described by the Neo-Hookean material model with anisotropic fiber contributions. The total strain energy density function is a sum of these elements, expressed in Eqs 2, 3: W = C 10 I 1 − 3 + 1 D J α 2 − 1 2 − ln J α + k 1 2 k 2 ∑ α = 1 N exp k 2 ⟨ E α ⟩ 2 − 1 , (2) with: E α = κ I 1 − 3 + 1 − 3 κ I 4 α − 1 . (3) In these expressions, C ₁₀ and D represent the material constants for the annular ground substance stiffness and compressibility. The parameters k ₁ and k ₂ characterize the exponential stress-strain relationship of collagen fibers, with κ adjusting their dispersion level. This modeling approach specifically accounts for fiber resistance in tension only. The term I _4α denotes pseudo-invariants associated with both the right Cauchy-Green deformation tensor and unit vectors in fiber directions. With N = 2, we considered two orientations of fibers within the AF.

The parameters k _1c , k _2c , k _1r and k _2r , as specified by Nicolini et al. (2022a), were utilized to represent the changes in fiber stiffness along both circumferential and radial directions (Figure 1). This approach is based on the findings of the studies by Ebara et al. (1996), Eberlein et al. (2001), Holzapfel et al. (2005), and Zhu et al. (2008). The introduction of the scaling factors α _r and α _c allowed the variation of the fiber angle in the radial and circumferential directions, respectively (Holzapfel et al., 2005; Zhu et al., 2008). The angle of the fibers is defined by their direction and the transverse plane (Nicolini et al., 2022a). A value of α _r = 0.1 indicates a 10% increase in fiber angle from layer 1 (most external layer) to layer 2, followed by 20%, 30%, and 40% increases in the transition to layers 3, 4, and 5, respectively. Similarly, α _c corresponds to the change in fiber angle in the circumferential direction within the same layer when moving from one subregion to the next, starting with subregion A. We conducted a mesh sensitivity study to choose appropriate mesh parameters that provide precise outcomes without significantly increasing the computational time. This analysis was based on the previous study from Nicolini et al. (2022b). We applied a pure moment load of 5 Nm in flexion and examined the outcomes for range of motion (RoM) and the simulation time across different mesh sizes and shape functions using hexahedral elements. Hybrid elements were implemented to account for the incompressible behavior of both NP and AF ground substance (Schmidt et al., 2006; Schmidt et al., 2007; Schmidt et al., 2012; Tang and Rebholz, 2011). The identified mesh settings were consistently assigned to all models to ensure an equal number of solid elements and nodes.

In the second approach, the annular lamellae were defined by uniaxial rebar-reinforced membrane elements. This modification involved the integration of additional M3D4 elements into the solid elements of the discretized AF, resulting in two families of reinforced fibers per lamella with criss-cross directions. The geometric properties of the fibers were obtained from literature references (Cassidy et al., 1989; Marchand and Ahmed, 1990; Schmidt et al., 2006; Noailly et al., 2011). As in the first model, the same scaling factors α _r and α _c were used to determine the change in fiber angle in radial and circumferential directions. For this model, NP and AF ground substance were defined using a Mooney-Rivlin material model (Schmidt et al., 2006; Zhong et al., 2006; Little et al., 2008; Liu et al., 2011; Tang and Rebholz, 2011; Schmidt et al., 2012; Park et al., 2013; Shahraki et al., 2015; Wu et al., 2016; Yang and O’Connell, 2017; Cai et al., 2020). The fiber membrane was attributed with the same material properties as the AF ground substance (Wang et al., 2021b). The definition for the rebar elements included two different methods, resulting in two models that were analyzed. The linear rebar model used linear-elastic material properties, in accordance with the specification of Holzapfel et al. (2005) that these properties are only active in the tensile direction. The initial Young’s Modulus for the linear-elastic fiber definition was calculated by identifying a fitting linear stress-strain curve for the given data from Shirazi-Adl et al. (1986). The nonlinear rebar model utilized the hyperelastic Marlow material definition combined with the nonlinear stress-strain data from Shirazi-Adl et al. (1986). To account for the change in fiber stiffness in radial and circumferential directions, additional scaling factors were introduced for the rebar models. λ served as a scaling factor for the stress-strain curve (nonlinear rebar model) and the Young’s Modulus (linear rebar model) in the outermost layer of subregion A (Schmidt et al., 2006). The stiffness variation in the circumferential direction was defined by λ _c , reflecting the equidistant decrease from anterior to posterior. For example, λ _c = −0.1 means that the fiber stiffness in regions B, C, D, and E within the same layer was 10%, 20%, 30%, and 40% less than the fiber stiffness in region A, respectively. The change in fiber stiffness in the radial direction was similarly treated with λ _r .

Load introduction was achieved by defining a reference point located 10 mm above the disc surface, following the study of Nicolini et al. (2022a). This node was used to simulate the load application during the in vitro experiments (Heuer et al., 2007b; Nicolini et al., 2022a). A coupling constraint was employed to connect the upper surface of the IVD to the reference point (Wang et al., 2021a). The loads applied in the simulations increased from 0 to the final load value using a ramp function. To replicate the fixed lower surface of the caudal vertebra in the in vitro experiments, the FE model of the IVD was fixed by a coupling constraint with an anchored reference point below the lower surface of the IVD. The simulations were performed using a static solver configured for nonlinear material behavior.

Prior to calibration, a sensitivity analysis was performed to identify the model parameters to be calibrated. We implemented a design of experiments approach to evaluate the effect of different input parameters on RoM and intradiscal pressure (IDP) during four load cases: flexion, extension, lateral bending, and axial rotation. The parameter vectors from Eqs 4, 5 were selected to be investigated in the sensitivity analysis: x H G O = C 10 n , C 01 n , C 10 a , HGO , k 1 , k 2 , κ , α , k 1 c , k 2 c , k 1 r , k 2 r , α c , α r , (4) x r e b a r = C 10 n , C 01 n , C 10 a , MR , C 01 a , MR , λ , ν , α , λ c , λ r , α c , α r . (5) The index n is used for parameters related to the NP, while a is used for those of the AF. Additionally, HGO and MR are acronyms for the material models HGO and Mooney-Rivlin, respectively. As the NP and AF ground substance are defined as incompressible structures, we did not consider D, which specifies the compressibility of the material, in the sensitivity analysis. To determine the sensitivity of the models to the different input parameters, we adopted the one factor at a time (OFAT) method because of its ability to identify the gross effects of input parameters and its advantages in terms of speed and simplicity (Saltelli et al., 2007). Each parameter was stepwise adjusted four times within its range, while the other parameters remained constant. To establish the parameter ranges for our analysis, we conducted a review of the relevant literature. For each parameter, we obtained median values from the literature where they were available and applicable, as summarized in Table 1. These medians served as the central values for the parameter ranges. To ensure a reasonable distribution, the lower and upper limits of the ranges were set at 50% and 150% of the corresponding median values, respectively (Zander et al., 2017; Du et al., 2021). Given the limited literature describing scaling methods for adjusting fiber angle and stiffness in radial and circumferential directions, we selected the median within these ranges as the central value from the established definitions. Similarly, for κ, which defines the fiber dispersion in the HGO material model, we employed the median as the central value from the available parameter range. At the beginning of the sensitivity analysis, an initial simulation run was performed for each model using the median values for the input parameters across the four load cases, each with a pure moment of 5 Nm. The initial simulations provided reference results for RoM and IDP, which were used to evaluate the response variations of the models resulting from the adjustments made to the different parameters during the sensitivity analysis. Notably, since the two rebar models are defined identically except for their fiber stiffness, only the linear rebar model was included in the sensitivity analysis.

To evaluate the simulation results, we applied the sensitivity score as formulated by Karajan and Ehlers (2007) and shown in Eq. 6: S R , P = percentage of the variation of the response percentage of the variation of the parameter . (6) For instance, a sensitivity score of S _R,P = 0.4 indicates that a 10% increase in the given parameter results in a 4% higher response. Each parameter received eight sensitivity scores, the result of considering four different load cases and two result metrics. These eight scores were calculated individually by averaging the scores from the four variations per parameter, as presented in Eq. 7 S P , L , M = 1 n ∑ i = 1 n Δ R L , M , i Δ P i . (7) S _P,L,M represents the sensitivity score for parameter P, load case L (flexion, extension, lateral bending or axial rotation), and result metric M (RoM or IDP). ΔP _i denotes the i-th variation of parameter P, and ΔR _L,M,i stands for the i-th variation of the response for load case L and result metric M arising from ΔP _i . Additionally, n specifies the number of variations per parameter, here 4. To assess the need for calibration of each specific parameter, we applied the critical score S _P,L,M = 0.1, following the recommendation by Karajan and Ehlers (2007), to guide our evaluation process. If the absolute value of at least one sensitivity score for a parameter was equal to or greater than 0.1, the parameter was included in the calibration procedure. Parameters without a sensitivity score exceeding the critical threshold were excluded from the calibration process.

The calibration process was formulated as an optimization problem to minimize the difference between FE simulation predictions and experimental data (Schmidt et al., 2006). We used in vitro experimental data obtained from the study conducted by Heuer et al. (2007b) in our calibration. Their study involved a stepwise reduction analysis at the L4-L5 segment using eight human specimens. The RoM was examined under four different load cases: flexion, extension, lateral bending, and axial rotation. Moment loading ranging from 1 Nm to 10 Nm was applied during the experiments. Since the results for the last reduction stages were based on only three specimens at 10 Nm, we limited our calibration process to data up to 7.5 Nm. As our study was focused on the IVD, we used only the experimental results from the reduction stage labeled “w/o ALL” which represented the entire disc and vertebral bodies. The calibration process was designed for error minimization, using the R-squared function to quantify the difference between experimental and numerical RoM data. The RoM curves for each loading direction were used to calculate the R-squared values.

The R-squared function is defined as: R 2 = 1 − ∑ i = 1 n y i − y ̂ i 2 ∑ i = 1 n y i − y ̄ i 2 , (8) where y _i represents the experimental RoM, y ̂ i the numerical RoM, y ̄ i the mean of the experimental RoMs, and n the number of observations (Mittelhammer, 2013). To determine the overall agreement for each model, we took the mean of the R² values for the different load cases.

We modified the previously published genetic algorithm developed by Nicolini et al. (2022a) to accommodate different types of model implementations for the fiber reinforcement, distinguishing between the use of the HGO fiber model and the rebar models, as visualized in Figure 2. The calibration of each model was divided into two parts to reduce the number of parameters calibrated per step and increase the efficiency of the process (Arora, 2017). The calibration strategy for the HGO fiber model involved an initial step, denoted as Cal. 1a in Figure 2, in which we calibrated the material properties of the NP, AF ground substance, and fiber stiffness using constant values for the fiber angle, which were obtained as median values from the given ranges. In this step, we used an R² threshold of 0.85. Subsequently, we further refined the model by optimizing the scaling and variation of the fiber angles in the second step (Cal. 1b), increasing the R² threshold to 0.9 (Nicolini et al., 2022a). To ensure consistency in defining the NP across different approaches, we utilized calibrated values for the NP from the HGO fiber model in the rebar models as well. The linear rebar model was subjected to a two-step calibration process with the specified R² thresholds (0.85 and 0.9). The stiffness of the ground substance and fibers were calibrated initially (Cal. 2a), followed by optimizing the fiber angle and its variation (Cal. 2b). Subsequently, the nonlinear rebar model was calibrated, excluding the AF ground substance parameters already established in the linear rebar model calibration, since the differences between these two models are limited to the fiber definition. This resulted in calibrating only the fiber stiffness (Cal. 3a) and fiber angles (Cal. 3b) of the last model, again using the two-step approach with the mentioned thresholds. We established calibration bounds for most parameters based on the sensitivity analysis (Table 2). For some parameters, as outlined in the subsequent text, we defined the bounds differently. The potential values for κ were between 0 and 0.33, following the range defined by the material model. The parameters varying the fiber stiffness in both radial and circumferential directions were set to a physiological range between −0.2 and 0, with a maximum change of 80% from subregions A to E and layers 1 to 5, respectively, following the findings of Holzapfel et al. (2005). The scaling factor for fiber stiffness in the rebar models was selected in the range of 0.3–2, as suggested by Schmidt et al. (2006). In accordance with the research of Holzapfel et al. (2005), the values of α _r could vary from 0 to 0.2, while α _c could range from 0 to 0.3.

The genetic algorithm was configured with a population size of 20. In each generation, six individuals were selected based on their R² values. These chosen individuals were then utilized to create the next generation. Through a crossover process that combined pairs of the selected individuals with their parameter values, four new individuals were created. Another four individuals were defined by mutation to introduce variability. In the mutation process, one of the six selected individuals was chosen to undergo a random alteration in which a parameter value was changed within its defined range. To avoid the potential stagnation of solutions and to ensure a diverse population, six new individuals were introduced per generation by immigration. The optimization procedure started by generating a random initial population within the established parameter ranges. The algorithm iterated through selection, crossover, mutation, and immigration until either the convergence criteria (R² threshold) was met or the maximum number of generations was reached. The maximum number of generations for the algorithm varied with the calibration steps: 20 generations for steps with more parameters (Cal. 1a and Cal. 2a), and 10 generations for the remaining steps with fewer parameters. A more detailed description of the calibration algorithm can be found elsewhere (Nicolini et al., 2022a).

The calibrated models were validated against two different published in vitro experimental datasets to verify their accuracy and reliability. First, we used the IDP measurements from the study of Heuer et al. (2007a). For this purpose, the models were simulated in their final parameter set resulting from the calibration with pure moment loading up to 7.5 Nm in flexion, extension, lateral bending, and axial rotation. Furthermore, in a secondary validation step, our models were evaluated using the experimental RoM results from Jaramillo et al. (2016). This additional dataset was included to increase the variety of specimens and enhance the robustness of the validation process. They conducted a stepwise reduction analysis on the L4-S1 section using a pure moment loading with up to 8 Nm on five human specimens. Specifically, we considered the results obtained at the reduction stage labeled “Wout_All”, where only the IVD remained. To determine the agreement between the numerical predictions and the experimental data, we employed the R² function, and calculated values for each model by applying Eq. 8, thereby quantifying their agreement with the validation data. We also evaluated the efficiency of different models by examining the computational time associated with each model. To achieve this, we considered the simulations of the different load cases with an applied moment of 8 Nm from the validation using the dataset of Jaramillo et al. (2016). The simulations were performed on an AMD Ryzen 7 7700X (8-core processor) with a clock speed of 4.5 GHz. GPU-acceleration, provided by an NVIDIA GeForce RTX 4090, was utilized to enhance computational efficiency.

To improve the robustness of our calibration approach, we conducted an additional calibration sequence by altering the order in which the models were calibrated. The second calibration started by modifying the material properties of NP, ground substance of AF, and the fiber stiffness parameters in the linear rebar model. Subsequently, we made further adjustments by calibrating the fiber orientation parameters. Using the established NP and AF ground substance parameter values from the linear rebar model, we calibrated the fiber stiffness and orientation in two distinct phases for the nonlinear rebar model. Then, we proceeded to calibrate the HGO fiber model using the NP parameter values from the linear rebar model. Applying a similar two-step approach, we separately optimized the AF properties and fiber orientation parameters. Following this recalibration process, we compared the R² values and the material parameter configurations obtained with those from our initial calibration.

In order to determine appropriate mesh settings, we performed a mesh sensitivity study. The analyzed outcome metric, RoM, showed a converging trend as the mesh was refined (Figure 3). The use of quadratic shape functions instead of linear ones resulted in fewer nodes being needed to achieve stable results. A mesh consisting of 93,327 nodes with C3D20H elements was chosen, which showed a deviation of 0.43% for RoM compared to the results obtained with the most refined quadratic mesh tested (605,894 nodes). However, the simulation time for the selected mesh was only 4.5% of the time required for the most refined quadratic mesh.

The results of the sensitivity analysis are presented in Figure 4, where each plot shows the averaged sensitivity scores for one model and one result metric (RoM or IDP) across the four load cases. The presented results encompass the HGO fiber model and the linear rebar model. The nonlinear rebar model, which was calibrated only for variations in fiber stiffness and angle, was not included in this sensitivity analysis. In the HGO fiber model, almost all parameters, with the exception of C_01n , reached the critical sensitivity score (S _P,L,M = 0.1) for at least one of the outcome metrics. For this reason, we did not calibrate C_01n within the HGO fiber model. Instead, the median value obtained from the literature (Table 1) was used in the configuration of the model. It was observed that C_10n and the parameters affecting the fiber stiffness (k _1c , k _2c , k _1r , k _2r ) did not exceed the sensitivity threshold within the RoM for all loading scenarios. However, they had a significant impact on intradiscal pressure (IDP) for at least one load case. For the linear rebar model, sensitivity scores for C_01n and ν showed negligible influence on the model’s response across all load cases and outcome metrics. Following this, both parameters were excluded from the calibration of the rebar models, and median literature values were used instead. Similar to the HGO fiber model, C_10n had a negligible effect on RoM across all load cases, but it surpassed the sensitivity threshold in IDP during extension. In contrast, C_01a did not reach the sensitivity criteria for IDP. But it had a considerable effect on the RoM in the linear rebar model during extension.

Table 3 presents the calibrated material parameter configurations for each model. During calibration, the HGO fiber model achieved the envisioned initial agreement with the experimental data (R² ≥ 0.85). The second calibration step, which focused on refining the variation of the fiber angle in both circumferential and radial directions, resulted in further improvement (R² = 0.95). The linear rebar model, although exhibiting the specified agreement after the first calibration step, could not meet the defined final fit with the experimental RoM data. The nonlinear rebar model was able to predict the experimental data with an R² value of 0.87, but did not surpass the predefined R² threshold. Altogether, the HGO fiber model showed superior calibration performance compared to the rebar models. The detailed results, showing the R² values for different load cases, are presented in Table 4. Figure 5 compares the RoM curves over applied moments for the models with the experimental data.

To validate and compare the models, we used IDP measurements from the experiments of Heuer et al. (2007a). The corresponding R² values using these measurements are summarized in Table 5; Figure 6 illustrates the comparison of IDP curves over applied moments. In all load cases except axial rotation, the HGO fiber model and the linear rebar model accurately reproduced the experimental in vitro results, with the HGO model achieving a marginally higher agreement (R² = 0.87) compared to the linear rebar model (R² = 0.86). It was observed that the nonlinear rebar model showed reduced accuracy, especially in the context of lateral bending, resulting in a lower overall agreement with the experimental data. For additional validation and comparison, the experimental RoM data reported by Jaramillo et al. (2016) was included. The HGO fiber model had the highest agreement with this RoM reference data (R² = 0.85), underscoring its superiority in accuracy over the rebar models in this specific context. Table 6 contains the detailed R² values, while Figure 6 visually presents the RoM curves of the models for the four different load cases and compares them with the experimental data from Jaramillo et al. (2016). When evaluating the computational efficiency, as detailed in Table 6, the HGO fiber model was found to outperform the rebar models, achieving a considerably shorter mean computational time (191 s). In contrast, the simulations using the linear and nonlinear rebar models took considerably more time to complete, with mean times of 461 s and 478 s, respectively. These times were obtained from the simulations of the different load cases with an applied moment of 8 Nm.

To assess the reliability of our calibration method, we repeated the process with a modified order. This resulted in different material configurations compared to the first calibration procedure, as shown in Table 3. However, all three models achieved similar agreement with the experimental RoM data from Heuer et al. (2007b) for both calibrations. Table 4 reveals that the HGO fiber model again had the highest R² value (0.93), followed by the linear rebar model (R² = 0.88) and then the nonlinear rebar model (R² = 0.86).