A novel approach for tetrahedral-element-based finite element simulations of anisotropic hyperelastic intervertebral disc behavior

Intervertebral discs are microstructurally complex spinal tissues that add greatly to the flexibility and mechanical strength of the human spine. Attempting to provide an adjustable basis for capturing a wide range of mechanical characteristics and to better address known challenges of numerical modeling of the disc, we present a robust finite-element-based model formulation for spinal segments in a hyperelastic framework using tetrahedral elements. We evaluate the model stability and accuracy using numerical simulations, with particular attention to the degenerated intervertebral discs and their likely skewed and narrowed geometry. To this end, 1) annulus fibrosus is modeled as a fiber-reinforced Mooney-Rivlin type solid for numerical analysis. 2) An adaptive state-variable dependent explicit time step is proposed and utilized here as a computationally efficient alternative to theoretical estimates. 3) Tetrahedral-element-based FE models for spinal segments under various loading conditions are evaluated for their use in robust numerical simulations. For flexion, extension, lateral bending, and axial rotation load cases, numerical simulations reveal that a suitable framework based on tetrahedral elements can provide greater stability and flexibility concerning geometrical meshing over commonly employed hexahedral-element-based ones for representation and study of spinal segments in various stages of degeneration.

1 Introduction

The human spine is a complex arrangement of passive and active tissues (hard and soft) that provides postural control, flexibility of motion, and protects the spinal cord (Kibler et al., 2006; Agur and Dalley, 2009; Brittanica, 2014). Among these tissues, the hydrated soft collagenous intervertebral discs (IVD) separating adjacent vertebrae of the spine form major load-bearing elements (Widmer et al., 2020) that provide cushioning, regulate force distribution, and facilitate motion between spinal vertebrae (Humzah and Soames, 1988; Adams and Roughley, 2006; Roberts et al., 2006; Widmer et al., 2019). This multi-component complex structure in conjunction with the mechanical loads that it experiences during various physical activities (Schultz and Andersson, 1981; Schultz et al., 1982) results in complex internal load transfer mechanisms, which are expected to influence spinal pathologies such as back disorders and pain (Pearcy et al., 1984; Kumar, 1990; Dvořák et al., 1991; Manchikanti, 2000; Thiese et al., 2014) as well as secondary complications after surgical interventions such as adjacent segment disease (Bertagnoli, 2011), pseudoarthrosis (Steinmann and Herkowitz, 1992), and screw loosening (Bredow et al., 2016). In this context, finite element (FE) based models encompassing various spinal components have gained greater attention in recent decades to study spine biomechanics (Noailly et al., 2005; Schmidt et al., 2006; Campbell et al., 2016; Dreischarf et al., 2014; Zander et al., 2009; Schmidt et al., 2013; Jaramillo et al., 2015; Maquer et al., 2015; del Palomar et al., 2008; Ayturk et al., 2010; Zander et al., 2017) with increasing applications towards pre-clinical/surgical studies (Baroud et al., 2003; Rohlmann et al., 2007; Boccaccio et al., 2008; Talukdar et al., 2021), evaluating the influence of intervertebral disc degeneration (Schmidt et al., 2007b; Ayturk et al., 2012; Cegoñino et al., 2014), and towards subject-specific investigations (Widmer, 2020; Pickering et al., 2021). These computationally powerful tools are particularly effective in combining hierarchic intricacies of complex spinal systems with material and geometrical non-linearities and a wide range of loading scenarios (Schmidt et al., 2013). However, the performance of these FE models (i.e., accuracy, computational efficiency, and robustness) is greatly influenced by 1) component-specific material models and the corresponding model parameters (in particular for soft tissues like the IVD), 2) accuracy and discretization of the three-dimensional geometry, and 3) numerical solution techniques employed.

Due to complex inner microstructure and internal hierarchy (see Section 2.1), the IVD exhibits highly non-linear behavior (Markolf and Morris, 1974; Goel et al., 1995; Ebara et al., 1996; Eberlein et al., 2001; Wagner and Lotz, 2004). Various mathematical models were proposed in the literature, with different mathematical treatments for the nucleus pulposus (NP) and annulus fibrosus (AF) components of the IVD. Example modeling approaches are linearised-elasticity-based models (Harkness, 1961; Haut and Little, 1972; Ueno and Liu, 1987; Smit et al., 1997; Baroud et al., 2003; Polikeit et al., 2003), non-linear composite models with a one-dimensional description of collagen fibers (Schmidt et al., 2006; 2007b; Rohlmann et al., 2007; Zander et al., 2009; Dreischarf et al., 2014; Pickering et al., 2021), microstructure-based continuum non-linear models (Eberlein et al., 2001, 2004; Ayturk et al., 2012, 2010; Jaramillo et al., 2015), and micromechanical models (Ghezelbash et al., 2021). Continuum material models are often employed to model the IVD as they are particularly advantageous to represent its material non-linearity while implicitly incorporating its anisotropy to the desired degree of accuracy (Holzapfel et al., 2015, 2005) besides offering relative ease for FE implementation. However, ambiguity with the corresponding material parameters prevails because various parameters suggested in the literature (Moramarco et al., 2010; del Palomar et al., 2008; Jaramillo et al., 2015; Wagner and Lotz, 2004; Eberlein et al., 2001, 2004) were each attuned to a specific set of experiments that were largely uni-directional in nature (Markolf and Morris, 1974; Goel et al., 1995; Ebara et al., 1996; Eberlein et al., 2001; Wagner and Lotz, 2004). This suggests a need for robust calibration using multi-directional experimental data to increase the model’s predictive abilities.

Various image-processing-based methods are in use to generate FE models for spinal segments. These approaches use computed tomography (CT) images of vertebrae (Moramarco et al., 2010; del Palomar et al., 2008; Rohlmann et al., 2007; Eberlein et al., 2004; Schmidt et al., 2006; Pickering et al., 2021; Jaramillo et al., 2015) and magnetic resonance imaging (MRI) scans of IVDs (Maquer et al., 2015, 2014), including those based on automatic segmentation (Caprara et al., 2021). Hexahedral elements (HE) often feature in the subsequent discretization of the resulting geometry, in particular of healthy IVDs (Eberlein et al., 2001, 2004; Jaramillo et al., 2015; Moramarco et al., 2010; Rohlmann et al., 2007; Pickering et al., 2021; Baroud et al., 2003; Schmidt et al., 2006, 2007b,a; Zander et al., 2009, 2017; Cegoñino et al., 2014). This is because, 1) HE mesh regular geometries easily and with fewer elements, while offering high solution accuracy (Benzley et al., 1995), 2) they can be arranged parallel to the IVD circumference with one of the local coordinate axes being tangent to it and therefore simplifying the identification of local collagen fiber directions, and 3) HE in the form of hourglass controlled reduced integration overcome the volumetric locking problem in incompressible soft materials (Bathe, 2006; Hughes, 2012).

However, degeneration-induced changes significantly affect the geometry of IVDs. Specifically, there can be a considerable drop in height, an accumulation of tears in the annulus region, and endplate effects (Adams and Roughley, 2006; Widmer, 2020). Furthermore, olisthesis or dorsal disc narrowing can result in strongly skewed or wedged-shaped IVDs. Such complex geometries are extremely difficult to reproduce with (homogeneously sized) HE and if meshed inaptly this can cause numerical instabilities. This mandates using computationally expensive and laborious high-quality HE meshes. In this regard, tetrahedral elements (TE) are commonly employed in literature to discretize complex geometries with relative ease due to their superior flexibility and adaptability (Bathe, 2006; Hughes, 2012; Schneider et al., 2019). Furthermore, volumetric locking can be addressed with a suitable choice of integration schemes (e.g., reduced, selective reduced) (Bathe, 2006; Hughes, 2012). Refined approaches in the context of volumetric locking issues with TE elements have also been explored (Joldes et al., 2009; Pagani et al., 2014). Finally, while implicit finite element analysis is generally faster for linear problems, explicit numerical solution techniques are often selected over implicit methods in addressing quasi-static problems (with negligible inertial effects) in FE methods (FEM) because 1) no iterations are required to evaluate solution variables, 2) evaluation of computationally expensive inverse stiffness is not required, 3) material and geometric non-linearities, as well as contacts, are handled better, 4) high levels of efficiency are possible with parallelization for analyses solved using multiple processors. Yet, an optimal choice of the time step is paramount to ensure the desired accuracy of numerical solutions while maintaining computational efficiency. While this time step is deformation-dependent (Ogden, 1997), it is traditionally prescribed as a suitably small constant (dependent on the material parameters) in typical explicit FE solvers such as Ansys LS-DYNA and Radioss for ease of implementation in a range of problems in mechanics (Freed et al., 2005; Wu et al., 2020). However, in highly non-linear anisotropic hyperelastic materials like IVD tissue, this time step can be noticeably influenced by the state of deformation and local material symmetry, suggesting a re-evaluation of the traditional approach.

This project focuses on establishing a novel and robust FE-based model formulation for spinal segments using TE in a hyperelastic framework. To this end, 1) a microstructure-based continuum anisotropic material formulation is utilized for the simulation of AF behavior in an explicit-time-integration-based numerical framework. 2) An adaptive time-stepping approach is proposed as a computationally efficient approximation to a refined deformation-dependent alternative (Ogden, 1997). Finally, 3) the performance of linear HE and TE is evaluated in terms of their accuracy and stability during flexion, extension, lateral bending, and axial rotation, using spinal FE models with non-degenerated, moderately, and severely degenerated IVDs. Also, a material parameter set is estimated using experimental data of spinal segments during the above load cases (Widmer et al., 2020) and an inverse FE-based approach.

2 Methods

2.1 Continuum material formulations for the IVD

The internal microstructure of the IVD consists of an inner NP enclosed by an outer AF and cartilage endplates anchoring the IVD to the vertebrae. Both AF and NP are predominantly filled with water (65%–90% in AF (Marcolongo et al., 2017) and 70%–88% in NP (Humzah and Soames, 1988; Marcolongo et al., 2017)) and a proteoglycan matrix into which collagen fiber bundles are embedded (Hashizume, 1980; Inoue, 1981; Cassidy et al., 1989). Collagen fiber bundles in AF are arranged into several layers of concentric lamellae with alternating orientations ranging between 25° and 45° (Cassidy et al., 1989; Adams and Roughley, 2006; Ambard and Cherblanc, 2009; Malandrino et al., 2013) about the transverse plane, while in NP they are randomly oriented in a homogeneous matrix (Hashizume, 1980; Inoue, 1981). The arrangement of the collagen fibers in AF is theorized to resist the high tensile hoop loads resulting from the hydrostatic pressure transferred from the nucleus pulposus during spinal compression by helping to absorb and redistribute stresses. The complex internal hierarchy and microstructure can be linked to the experimentally determined non-linear mechanical behavior observed both at the component and the tissue level (Markolf and Morris, 1974; Goel et al., 1995; Ebara et al., 1996; Eberlein et al., 2001; Wagner and Lotz, 2004; Holzapfel et al., 2005). In this study, IVD is modeled as a multi-component system with individual constitutive material models for AF and NP, to incorporate its complex three-dimensional microstructure. To this end, both AF and NP are modeled as hyperelastic solid bodies utilizing an invariant-based formulation (Spencer, 1984; Holzapfel, 2002) wherein pressure and displacement are decoupled for numerical ease (Flory, 1961; Ogden, 1978; Simo and Hughes, 2006).

2.1.1 Kinematics and preliminaries

In line with the standard notation in continuum mechanics, let the configuration of a body $\mathcal{B}$ in the reference and current (at time t) states be denoted by ${\varkappa }_R(\mathcal{B})$ and ${\varkappa }_t(\mathcal{B})$, respectively. Each material point of $\mathcal{B}$ corresponds with the positions $\mathit{X}\in {\varkappa }_R(\mathcal{B})$ and $\mathit{x}\in {\varkappa }_t(\mathcal{B})$, which are linked through the mapping $\mathit{x}={\mathit{\chi }}_{{\varkappa }_R}(\mathit{X},t)$. The deformation gradient F(X, t), its determinant J(X, t), and the right Cauchy-Green tensor C(X, t) are defined through

$\mathbf{F}=\text{Grad}{\mathit{\chi }}_{{\varkappa }_R},J=\text{det}\mathbf{F}>0,\mathbf{C}={\mathbf{F}}^{\mathrm{T}}\mathbf{F},(1)$

where the dependence on location and time is understood. Densities in the current (ρ) and reference states (ρ₀) are related through ρ = Jρ₀. The strain energy density function (SEDF) Ψ of each hyperelastic material component depends on F through the right Cauchy-Green tensor C = F^TF which, in the decoupled form, becomes (Holzapfel et al., 2000)

$\mathrm{\Psi }=\bar{\mathrm{\Psi }}\left(\bar{\mathbf{C}}\right)+\tilde{\mathrm{\Psi }}\left(J\right),(2)$

where $\bar{\mathbf{C}}={\bar{\mathbf{F}}}^{\mathrm{T}}\bar{\mathbf{F}}$ and $\bar{\mathbf{F}}=J^{-1/3}\mathbf{F}$. $\bar{\mathrm{\Psi }}$ and $\tilde{\mathrm{\Psi }}$ represent the distortional and dilatational isotropic strain energies, respectively. Furthermore, for an anisotropic hyperelastic material reinforced with n families of fibers whose directions are specified by unit vectors M_i, i = 1, 2, … , n, the SEDF can be given as (Holzapfel and Weizsäcker, 1998; Holzapfel et al., 2000; Truesdell and Noll, 2004)

$\mathrm{\Psi }=\bar{\mathrm{\Psi }}\left(\bar{\mathbf{C}}\right)+\tilde{\mathrm{\Psi }}\left(J\right)+\widehat{\mathrm{\Psi }}\left(\bar{\mathbf{C}},{\mathbf{M}}_i\right),(3)$

where M_i≔M_i ⊗M_i and $\widehat{\mathrm{\Psi }}$ represents the distortional anisotropic strain energy. The second Piola-Kirchhoff stress (S) and the Cauchy stress (σ) are then obtained as

$\mathbf{S}=2\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}},\mathit{\sigma }=J^{-1}{\mathbf{F}\mathbf{S}\mathbf{F}}^{\mathrm{T}}(4)$

Finally, Cartesian components of the material $({\mathcal{C}}_{ijkl})$ and spatial (c_ijkl) elasticity tensors are obtained as (Ogden, 1997; Holzapfel, 2002)

${\mathcal{C}}_{ijkl}=2\frac{\partial S_{ij}}{\partial C_{kl}},c_{ijkl}=J^{-1}{F}_{iP}{F}_{jQ}{F}_{kR}{F}_{lS}{\mathcal{C}}_{PQRS}(5)$

2.1.2 Nucleus pulposus

The behavior of the NP is modeled using a compressible Mooney-Rivlin type formulation given as (Holzapfel, 2002)

${\mathrm{\Psi }}_{\mathrm{N}\mathrm{P}}=b_{10}\left({\bar{I}}_1-3\right)+b_{01}\left({\bar{I}}_2-3\right)+\frac{k_{\mathrm{N}\mathrm{P}}}{2}{\left(J-1\right)}^2,(6)$

where ${\bar{I}}_1≔\mathrm{t}\mathrm{r}\bar{\mathbf{C}}$ and ${\bar{I}}_2≔\mathrm{t}\mathrm{r}{\bar{\mathbf{C}}}^{-1}$ are kinematic invariants, b₁₀ and b₀₁ are positive material constants with the units of stress and k_NP is the bulk modulus. Herein, it is noted that the NP is modeled as an isotropic material despite the presence of collagen fibers because of their random and homogeneous distribution in the matrix [see also Schmidt et al. (2006, 2007a,b); Ayturk et al. (2012)] The material model parameters for NP adopted from Schmidt et al. (2006, 2007a) and Ayturk et al. (2012) are listed in Table 1.

TABLE 1

TABLE 1. Material model parameters of NP and AF, where k_NP and k_AF are obtained from their respective Poisson’s ratios of 0.495 (Ayturk et al., 2012) and 0.45 (Goel et al., 1995).

2.1.3 Annulus fibrosus

To model the anisotropic mechanical response of the AF, a modified Mooney-Rivlin type formulation incorporating contributions from two collagen fiber families M₁ and M₂ is utilized based on (Eberlein et al., 2001; Eberlein et al., 2004; Moramarco et al., 2010). Therefore, $\mathrm{\Psi }(\bar{\mathbf{C}},J,{\mathbf{M}}_1,{\mathbf{M}}_2)$ from Eq. 3 becomes

${\mathrm{\Psi }}_{\mathrm{A}\mathrm{F}}=c_{10}\left({\bar{I}}_1-3\right)+c_{01}\left({\bar{I}}_2-3\right)+\frac{k_{\mathrm{A}\mathrm{F}}}{2}{\left(J-1\right)}^2+\sum _{i=1}^2\frac{a_1}{a_2}\left\{\text{exp}\left(a_2{⟨{\bar{I}}_{\mathrm{M}}^{\left(i\right)}-1⟩}^2\right)-1\right\},(7)$

where the invariant ${\bar{I}}_{\mathrm{M}}^{(i)}≔\mathrm{t}\mathrm{r}(\bar{\mathbf{C}}{\mathbf{M}}_i)$ measures the squared stretch of fibres along M_i. c₁₀, c₀₁ and a₁ are material parameters with units of stress, a₂ is a dimensionless constant, k_AF is the bulk modulus, and ⟨x⟩ = (|x| + x)/2. It is noted that only tensile stretch of collagen fibers is considered due to their crimped structure (Cassidy et al., 1989). In the current study, a homogeneous distribution of collagen fiber bundles is assumed, despite their alternating orientation in the lamellae, exploiting the periodic and concentric nature of the lamellae (Eberlein et al., 2001, 2004; Moramarco et al., 2010). Furthermore, M₁ and M₂ are assumed to be at an average orientation of ±30° about the transverse plane (Goel et al., 1995; Eberlein et al., 2001; Urban and Roberts, 2003; Moramarco et al., 2010). While the orthotropic material symmetry of the AF described by M₁ and M₂ is evident, the same can be described by the directions M₁ + M₂ and M₁ − M₂ with relative ease¹. The material model parameters of the AF used for the IVD FE model testing (Section 2.3.1) are listed in Table 1 and are in accordance with previous literature findings (Schmidt et al., 2006; 2007b; Ayturk et al., 2012).

2.2 Adaptive time step

Invoking the theory of infinitesimal waves and vibrations in unbounded materials in the context of finite deformations (Ogden, 1997), the acoustic tensor $\tilde{\mathbf{Q}}(\mathit{n})$ and the wave speed v of a plane given by v = mf (n⋅x − vt) are related through

$\rho v^2=\left[\tilde{\mathbf{Q}}\left(\mathit{n}\right)\mathit{m}\right]\cdot \mathit{m},{\tilde{Q}}_{jl}=A_{ijkl}{n}_i{n}_k,A_{ijkl}=c_{ijkl}+{\sigma }_{ik}{\delta }_{jl},(8)$

where unit vectors n and m denote the direction of wave propagation and polarisation of the wave, respectively. {n₁, n₂, n₃}, {m₁, m₂, m₃}, and δ_ij, respectively, represent the Cartesian components of m, n, and the identity tensor. A_ijkl is stiffness tensor. For longitudinal (P-) waves defined through n = m, Eq. 8 simplifies to

$v=\sqrt{\frac{J}{{\rho }_0}\left(c_{ijkl}+{\sigma }_{ik}{\delta }_{jl}\right)n_i{n}_j{n}_k{n}_l}(9)$

Noteworthy, the local wave speed in anisotropic materials is strongly influenced by the local fiber directions (Eqs 5, 7, 8). To this end, let n₁ = cos ϕ sin θ, n₂ = sin ϕ sin θ, and n₃ = cos θ without loss of generality and with ϕ ∈ [0, π) and θ ∈ [0, π). Then, the maximum local wave speed (v_max) and the corresponding direction of propagation for a given state of deformation (F, σ) can be deduced by maximizing a suitable objective function $\widehat{v}(\varphi ,\theta )$² i.e.,

$v_{\max }=\max \widehat{v}\left(\varphi ,\theta \right)(10)$

While the estimation of wave speed using Eq. 10 is essential for numerical simulations using the explicit FEM (see Section 2.3), it also increases the overall computational effort. Therefore, in this study v_max is approximated as

$v_{\max }=\max \tilde{v}\left(\tilde{\mathit{M}}\right),\tilde{\mathit{M}}\in \mathcal{M}≔\left\{{\mathit{M}}_1,{\mathit{M}}_2,{\mathit{M}}_1+{\mathit{M}}_2,{\mathit{M}}_1-{\mathit{M}}_2,{\mathit{M}}_1\times {\mathit{M}}_2\right\}(11)$

i.e, as the maximum of the wave speeds along the set of directions $\mathcal{M}$, by exploiting (i) the local orthotropy of the AF and (ii) substantially higher collagen fiber stiffness over the matrix as a rectification over the commonly used constant dilatational wave speed approach (Freed et al., 2005; Wu et al., 2020). We discuss this in Section 3, wherein the choice of $\mathcal{M}$ is compared against the theoretically estimated propagation direction (Eq. 10).

2.3 FE modeling

Numerical simulations were performed on three different FE models of lumbar spine segments representing various stages of IVD degeneration (Supplementary Figure S1).

The 3D geometrical mesh information of individual bony structures (i.e., cranial and caudal vertebra of each segment) was obtained from manual segmentation of CT images (Philips Brilliance 64, Philips Healthcare, Cleveland, OH, United States) using the 3D Slicer (V4.8.1) (3D Slicer, 2021; Fedorov et al., 2012) software (Figure 1). Statistical shape models were transformed onto this outcome by utilizing specific landmarks and invoking a non-rigid registration approach (Caprara et al., 2021). This information was then utilized in conjunction with custom-made scripts in MATLAB^® (The MathWorks Inc., Natick, MA, United States) to generate the 3D geometry of the enclosed IVD. For this purpose, the surface geometry of the cranial vertebra’s lower endplate and the caudal vertebra’s upper endplate were used. The circumferential profile of the IVD was shaped by a modest radial translation of the corresponding nodes to shape a gentle outward curvature. The location, shape, and size of the NP near the center of the IVD were defined based on anatomical studies (Pooni et al., 1986; O’Connell et al., 2007).

FIGURE 1

FIGURE 1. Finite element model generation of lumbar spine segments. (A) CT images of human cadaveric spines were used to obtain (B) vertebral 3D models through segmentation. (C) The endplates of the resulting vertebral surfaces were used to create the intervertebral disc geometry in between. Local collagen fiber directions M₁ and M₂ of the AF are displayed. {e₁, e₂, e₃} and $\{{\widehat{\mathit{e}}}_1,{\widehat{\mathit{e}}}_2,{\widehat{\mathit{e}}}_3\}$ are respectively the global and local orthonormal bases. (D) Finally, the components of the spinal segment were assembled into one FE model composed of bony and soft tissue (translucent upper vertebra for illustration purposes only).

In this study, the load cases of 1) flexion, 2) extension, 3) lateral bending, and 4) axial rotation in the three FE models of spinal segments were considered, as illustrated in Figure 2. To this end, the upper cranial vertebra was subjected to a moment of 5 Nm about an axis through the centroid of the IVD and normal to the sagittal plane for (1) and (2), frontal plane for (3), and transverse plane for (4). The IVD was rigidly connected to both vertebrae at their respective interfacing surfaces through a nodal tie constraint and the lower caudal vertebra was allowed only to rigidly translate in the moment plane. The purpose of this last constraint was to closely replicate the conditions of previously conducted biomechanical experiments (Widmer et al., 2020; Cornaz et al., 2021).

FIGURE 2

FIGURE 2. Load cases of (A) flexion/extension in the sagittal plane, (B) lateral bending in the frontal plane, and (C) axial rotation in the transverse plane for an applied moment of 5 Nm about the illustrated axis of rotation through the geometric center of the IVD.

The critical time step for Radioss explicit FE solver is estimated as

$\mathrm{\Delta }{t}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\mathcal{F}\frac{l_c}{v_{\mathrm{m}\mathrm{a}\mathrm{x}}},(12)$

Invoking Eq. 11 where l_c and $\mathcal{F}=0.9$ are the characteristic element length deduced from the FE mesh and a multiplicative adjustment factor, respectively. The choice of 0.9 for this multiplication factor was based on a rational compromise between maintaining a near-optimal computational effort and providing a reasonable buffer below the estimated minimum time step.

2.3.1 IVD degeneration and solid element type

We tested the described approach for IVD modeling by generating three different spinal segment geometries corresponding to L4-L5, L1-L2, and L2-L3 encompassing non-degenerated, moderately, and severely degenerated IVDs, respectively. Numerical simulations were performed on these three different FE models of lumbar spine segments with the IVD geometry being volumetrically discretized once with linear HE and once with linear TE. While the volumetric discretization of the IVD using hexahedral elements was performed in MATLAB, Hypermesh [HyperMesh, version 2017.2, Altair Engineering Inc., Troy, United States (HyperMesh, 2017)] was utilized for the tetrahedral-element-based discretization of the same.

All mechanical analyses were performed through a FE simulation of the corresponding boundary value problems with the explicit FE solver Radioss (2019). The load cases considered in this regard were implemented on a domain comprising the two vertebrae encompassing the IVD. Due to the considerable difference between the mechanical stiffness of the spinal vertebrae and the IVD (Lu et al., 1996; Baroud et al., 2003; Schmidt et al., 2006; Rohlmann et al., 2007), the former are modeled as rigid bodies. For both of the considered solid element types, a built-in Mooney-Rivlin-type material model in Radioss was used for NP, while the material model for AF followed the description of Sections 2.1 and 2.2 and was implemented through a user-defined material subroutine. Herein, the local collagen fiber directions (M₁ and M₂) in the local orthonormal basis $\{{\widehat{\mathit{e}}}_1,{\widehat{\mathit{e}}}_2,{\widehat{\mathit{e}}}_3\}$ are defined as

${\mathit{M}}_1=\left[\cos \mathrm{\Phi },\sin \mathrm{\Phi },0\right],{\mathit{M}}_2=\left[\cos \mathrm{\Phi },-\sin \mathrm{\Phi },0\right],(13)$

where Φ represents the orientation of the fibers about the local transverse plane spanned by ${\widehat{\mathit{e}}}_1$ and ${\widehat{\mathit{e}}}_3$. These vectors are then related to the element-specific coordinate system $\{{\bar{\mathit{e}}}_1,{\bar{\mathit{e}}}_2,{\bar{\mathit{e}}}_3\}$ through an orthogonal transformation ${\mathbf{Q}}_R:\{{\widehat{\mathit{e}}}_1,{\widehat{\mathit{e}}}_2,{\widehat{\mathit{e}}}_3\}\to \{{\bar{\mathit{e}}}_1,{\bar{\mathit{e}}}_2,{\bar{\mathit{e}}}_3\}$. The element-specific coordinate system is defined for the solver and based on the element edges and the sequence of the corresponding nodes (Radioss, 2019). The Euler angles of Q_R (Varshalovich et al., 1988) were determined in a pre-processing step performed in MATLAB with custom scripts using elemental edge (for HE) and orientation (for HE and TE) information and were provided as an initial input to the FE solver. Noteworthy, elemental collagen fiber orientations were determined by considering their tangency to the IVD outer circumference, a consequence of the concentric alignment of the lamellae (Holzapfel et al., 2005). Supplementary Table S1 summarizes the details of the FE discretization for all the models utilized in this study. Herein, it is noted that the vertebrae are discretized using shell-type elements with a thickness of 0.05 mm to associate non-zero inertia. Furthermore, a homogeneous FE mesh was considered to facilitate ease of meshing. Mesh dependency for the case of flexion is explored in Supplementary Figure S2 of the supplementary material wherein the mesh size was varied over 2.2 mm for coarse and 1.0 mm for refined meshes. Herein, it is noted that at 1.4 mm nominal size the maximum discretization error in relevant parameters was under 3% compared to the refined mesh. Therefore, 1.4 mm was chosen as a reasonable compromise between accuracy and computational efficiency (Supplementary Figure S2).

To compare the results of maximum longitudinal wave speeds determined by Eqs 10, 11 with each other, numerical simulations were first performed on the FE model of a spinal segment with non-degenerated IVD with TE and using a small value of Δt_crit³. Thereafter, the desired elemental state variables ({F,σ}^e) at the peak loading state were extracted. Finally, these output variables were used in conjunction with Eqs 10, 11 to estimate the corresponding elemental wave speeds. To this end, all four load cases were considered and material parameters based on the results of Eberlein et al. (2001); Ayturk et al. (2012); Goel et al. (1993, 1995) were assumed for AF.

Numerical simulations were performed on the FE models considering two different mesh element types, i.e., linear hexahedral and linear tetrahedral, and the corresponding results were compared with each other in terms of the distributions of pressure p, isotropic energy density $\bar{\mathrm{\Psi }}$, and anisotropic energy density $\widehat{\mathrm{\Psi }}$. p is defined as (Rajagopal, 2015)

$p=-\frac{tr\left(\mathbf{\sigma }\right)}{3}(14)$

Further, simulation results (i.e., stability and load-displacement) obtained using HE and TE types were compared for various degeneration states (Supplementary Figure S1A). Simulations providing results for the prescribed load without extreme deformation and buckling of single elements (causing the simulation to stop) were considered to be stable.

2.3.2 AF material parameters and inverse FEM

A set of material parameters of AF were estimated using inverse FEM and the calculations were performed on the spinal segment model with non-degenerated IVD meshed with TE elements (Supplementary Table S1). The experimental data for the four load cases were obtained from Widmer et al. (2020) and the average results for 31 non-degenerated spine segments were considered to be the experimental reference.

Briefly, the desired material parameter set (p_opt) of the AF is determined using an optimization algorithm implemented in MATLAB that iteratively minimizes the difference between numerically simulated responses from FE analyses performed in Radioss and the corresponding experimental data (see e.g., Ahn and Kim (2010); Kuravi et al. (2021); Böl et al. (2013)). To this end, the sequential quadratic programming algorithm (Wright and Nocedal, 2006) implemented in MATLAB’s fmincon function was invoked to minimize the objective function given as

$\mathcal{O}\left(\mathit{p}\right)=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(\frac{{\mathrm{\Theta }}_i^{\exp }-{\mathrm{\Theta }}_i^{\mathrm{s}\mathrm{i}\mathrm{m}}\left(\mathit{p}\right)}{{\mathrm{\Theta }}_n^{\exp }}\right)}^2},(15)$

where ${\mathrm{\Theta }}_i^{\exp }$ and ${\mathrm{\Theta }}_i^{\mathrm{s}\mathrm{i}\mathrm{m}}(\mathit{p})$, respectively, denote the experimental and simulated range of motion for a given parameter set p while n denoted the number of load increments. Figure 3 summarises the optimization algorithm.

FIGURE 3

FIGURE 3. Flowchart of the MATLAB-driven inverse-FEM-based optimization algorithm. The material properties characterizing the AF behavior are optimized towards the best agreement between experimental measurements and results of finite element simulations emulating the experiments. ROM: range of motion.

3 Results

In the first part of the results section, the outcome of time step estimation using a deformation-dependent and computationally efficient method are reported (Section 3.1). Thereafter, numerical simulation results using HE and TE are presented and compared with each other for the considered four load cases and for models with varying levels of degeneration (Section 3.2). Finally, material parameters based on experimental results obtained for non-degenerated IVDs through inverse FEM are reported (Section 3.3).

3.1 Time step estimation for explicit FEM

Figures 4A,B, respectively, depict probability density functions of the ratio of wave speeds and the corresponding angular separation estimated from Eqs 10, 11 for all elements of AF (Supplementary Table S1). The corresponding mean $(\bar{\mu })$ and standard deviation (SD) values are given in Table 2.

FIGURE 4

FIGURE 4. Comparing the theoretically predicted and approximated longitudinal wave speeds in non-degenerated AF. Plot (A) shows the ratio of wave speeds from Eq. 10 and the maximum from the directions in $\mathcal{M}$ and plot (B) shows the angular separation between the corresponding propagation directions, both expressed as probability density functions. $\bar{\mu }$ and SD, respectively, range in (1.033–1.038), (0.0281–0.0510) for (A), and (10.56–11.12), (4.75–6.25) for (B).

TABLE 2

TABLE 2. Comparing theoretical and approximated (P-) wave speed parameters.

3.2 Spinal model predictions for HE and TE

Figure 5–7, respectively, compare the state variables of mechanical pressure, isotropic (distortional), and anisotropic contributions to the SEDF of the numerical simulations obtained from a non-degenerated IVD meshed either with HE or with TE. Normalized histograms and empirical cumulative distribution functions (ECDF) were used to illustrate the distribution of the values in the elements. ECDF is a continuous function depicting the number of observations in percentile i.e., the percentage of observations that are less or equal to the value at a given point on X-axis.

FIGURE 5

FIGURE 5. Comparing the mechanical pressure expressed as probability (normalized histogram) and ECDF for flexion (A,E), extension (B,F), lateral bending (C,G), and axial rotation (D,H) for HE and TE types.

FIGURE 6

FIGURE 6. Comparing the isotropic contribution to the SEDF expressed as probability (normalized histogram) and ECDF for flexion (A,E), extension (B,F), lateral bending (C,G), and axial rotation (D,H) for HE and TE types.

FIGURE 7

FIGURE 7. Comparing the anisotropic contribution to the SEDF expressed as probability (normalized histogram) and ECDF for flexion (A,E), extension (B,F), lateral bending (C,G), and axial rotation (D,H) for HE and TE types.

Only a modest variation was inferred for pressure (Figures 5A–D), isotropic (Figures 6A–D) and anisotropic (Figures 7A–D) energy densities, illustrated as normalized histograms. The corresponding ECDF exhibited 3.8%–15.6% (Figures 5E–H), 1.6%–5.6% (Figures 6E–H), and 0.8%–7.9% (Figures 7E–H) variations, respectively, for the above load cases. For the applied torque of 5 Nm, the range of motion (ROM) differed moderately between the element types with a variation of 2.7%, 8.2%, 9.1%, and 9.9% for flexion, extension, lateral bending, and axial rotation, respectively.

Numerical simulation results (i.e., load-displacement curves) from FE models for non-, moderately, and severely degenerated spinal segments (Supplementary Figure S1) using HE and TE types were compared in terms of stability. All simulations converged and yielded physically meaningful results except for severely degenerated instances, where the HE approach failed in flexion and generally yielded high computational times (in extension in particular).

3.3 Calibration of material constants for non-degenerated IVDs

The optimal constitutive model parameter set p_opt for the non-degenerated AF is listed in Table 3 and was deduced following the termination of the optimization algorithm (Figure 3). Figure 8 depicts the good agreement between the numerically simulated and the corresponding experimental data for all four load cases.

TABLE 3

TABLE 3. Material model parameters of the AF obtained through inverse FEM and with the experimental results published in Widmer et al. (2020).

FIGURE 8

FIGURE 8. Comparison of experimental results vs. numerically simulated responses for (A) flexion, (B) extension, (C) lateral bending, and (D) axial rotation load cases. The shaded grey area covers the range of mean experimental ROM (dotted black line) ± one standard deviation. Experimental data were obtained from 31 non-degenerated spine segments (Widmer et al., 2020).

4 Discussion

The goal of this work is to develop and validate state-of-the-art spinal FE models (developed from CT scan images) using a hyperelastic formulation for the IVDs and employing TE elements. To this end, continuum material models were used for the NP and the AF in FE models employing an adaptive refined time-stepping approach, in contrast to the traditional approaches (Freed et al., 2005; Graff, 2012; Wu et al., 2020). Furthermore, numerical simulations were performed and verified for their stability and accuracy for flexion, extension, lateral bending, and axial rotation load cases on FE models of pristine as well as moderately and severely degenerated spinal segments. Finally, a set of material constants meant to describe the average behavior of non-degenerated AF tissue was found with an optimization approach.

4.1 Time step estimation for explicit FEM

In a variety of problems in non-linear mechanics, Δt_cric is often prescribed as a small enough constant to ensure the stability and accuracy of the numerical solutions. However, in the presence of material, geometrical, and contact non-linearities, much smaller time steps are generally utilized which also increases the involved computational effort. Though in the current study the latter aspects are not involved, large deformations can be expected due to material non-linearities including influences from anisotropy. This can greatly impact Δt_cric owing to the state of deformation and stress. In this context, theoretically estimated Δt_cric (Eq. 10) can be employed as an alternative to choosing small but arbitrary time steps. However, such an approach involves a computationally expensive iterative optimization. To this end, the approximation method presented in this work (Eq. 11) provides an excellent alternative, in that the obtained results differ only modestly from the theoretical estimates, i.e., $<$10% and $\approx 10°$ in the magnitude of wave speed and the corresponding direction of propagation, respectively. Noteworthy, these differences were estimated considering all four load cases, thus involving over 2 × 10⁵ instances of {F,σ}^e (Supplementary Table S1), which is a reasonable sample size.

4.2 Spinal model predictions for HE and TE

Only moderate differences were observed between numerical simulation results obtained with HE and TE, although the former offers increased accuracy over the latter in general (Hughes, 2012). Herein, ROM, isotropic and anisotropic energy densities differed modestly between both element types for all load cases, while for mechanical pressure moderate differences were observed (in excess of 10% for axial rotation).

Numerical simulations on non-degenerated, moderately, and severely degenerated spinal segments reveal that TE-based FE models were stable and computationally efficient in all three scenarios for all four load cases. In contrast, HE-based FE models exhibited stability only for non-degenerated and moderately degenerated scenarios involving relatively smooth IVD geometries and failed in the case of a severely degenerated IVD with a highly skewed geometry. These results reinforce the superiority of TE elements for simulations involving complex geometrical shapes in line with Schneider et al. (2019) while maintaining a homogeneous element size. Although similar results can be expected from HE types with refined and inhomogeneous meshes, the generation of these can be labor-intensive and computationally expensive.

In light of the above results involving commonly encountered load cases and various degenerative states of IVD, it is proposed that the presented TE-based FE models with component-specific models and a revised explicit time step offer a robust and computationally efficient alternative to studying the mechanics of the spine in general.

4.3 Calibration of material constants for non-degenerated IVDs

The AF exhibits an inhomogeneous microstructure with varying collagen directions within the lamellae. Moreover, these directions along with the amount of water in the AF are also reported to be influenced by its degenerative states. For instance, collagen fiber bundles are less organized in severely degenerated IVDs. However, in this work, the AF was represented as a micro-structurally homogeneous hyperelastic material using a 3D continuum formulation with two constant preferential directions. While this is in line with many state-of-the-art approaches [e.g., Eberlein et al. (2001); Schmidt et al. (2006; 2007a); Jaramillo et al. (2015); Ayturk et al. (2010; 2012)], noteworthy, it is more appropriate to associate such a characterization of the IVD tissue with non-degenerated states. To this end, the corresponding material parameters of the IVD were estimated using the experimental data obtained from spinal segments in a non-degenerated state. Herein, mean curves were utilized for the optimization-algorithm-driven inverse FEM, for simplicity.

Furthermore, many previous studies have utilized uni-axial experimental data for calibrating material models. However, in view of the complex IVD microstructure, model calibration-based multi-axial experimental data is highly desirable for it not only enhances the predictive abilities for generic load cases but also towards highly valued applications in clinical and subject-specific studies (Widmer, 2020; Pickering et al., 2021; Fasser et al., 2022). Therefore, in this study mean experimental data from flexion, extension, lateral bending, and axial rotation loads were employed to appropriately estimate IVD material parameters utilizing a state-of-the-art inverse FEM approach employing explicit time stepping method and driven by an optimization algorithm. Noteworthy, in the case of degenerated IVDs, the anisotropy component of SEDF in Eq. 7 can be suitably altered by employing the generalized structural tensor approach (Gasser et al., 2006; Holzapfel et al., 2015) to represent various degrees of local anisotropy.

4.4 Limitations

Intervertebral disc degeneration is reported to cause irreversible morphological changes such as the appearance of inhomogeneous tears and delaminations, increased disorganization in the AF microstructure (Urban and Roberts, 2003; Adams and Roughley, 2006), biochemical changes such as a decrease in water content, stiffening of the AF (Ebara et al., 1996), and geometrical changes such as an irregular but substantial reduction in height (Frobin et al., 1997); see Urban and Roberts (2003); Adams and Roughley (2006) for more details. Addressing the morphological and biochemical aspects mandates a rigorous mathematical framework involving inelastic, time-dependent, and multi-phasic effects (e.g., Cegoñino et al. (2014); Yang and O’Connell (2019)) and is beyond the scope of the current work. In this study, readily available geometry details extracted from CT scans were considered and the performance of the corresponding spinal segment FE models was explored in the realm of hyperelasticity. Furthermore, it is noted that TE-based FE discretization employed here is known to exhibit overly stiff behavior in pure displacement formulation when Poisson’s ratio approaches 0.5. The use of higher-order elements with reduced integration methods and mixed element formulations can reduce volumetric locking and improve the accuracy of numerical simulations (Joldes et al., 2009), the exploration of which is beyond the scope of this contribution.

5 Summary and conclusion

In this work, a continuum hyperelastic anisotropic material model was utilized to represent the AF component of the IVD. To facilitate the computational efficiency of explicit time-stepping method for spinal segment FE models, a novel approach for adaptive time step approximation was used and its proximity to theoretical estimates was evaluated. Furthermore, the effectiveness of TE-based FE models over HE-based ones was verified for various load cases when dealing with complex shapes of degenerated IVDs. Finally, a material parameter set was determined using inverse FEM and experimental data from flexion, extension, lateral bending, and axial rotation of non-degenerated human cadaveric spinal segments. Integrating the proposed approach of time step estimation with appropriately formulated TEs enables time-efficient modeling of even highly degenerated human IVD anatomies. Such anatomies often present severe challenges to automated finite element meshing strategies and/or result in suboptimal computational efficiencies to severely limit the clinical translation of patient-specific FE modeling approaches. The proposed work offers a path forward to overcoming these obstacles, without compromising on numerical accuracy. The subroutine developed in this work and used for the description of the AF material model is incorporated into Altair Radioss and will be released in the next official release in 2023.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving human participants were reviewed and approved by Kantonale Ethikkommission, Kanton Zürich. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author contributions

M-RF: Formal analysis, computational investigation, software, methodology, writing–original draft. RK: Formal analysis, computational investigation, software, methodology, writing–original draft. MB: Writing–review and editing, software. JS: Writing–review and editing, supervision, project administration, funding acquisition. MF: Writing–review and editing, supervision, project administration, funding acquisition. JW: Methodology, writing–original draft, conceptualisation, supervision, project administration, funding acquisition.

Funding

This work is part of “SURGENT” under the umbrella of University Medicine Zurich/Hochschulmedizin Zürich. Open access funding was provided by ETH Zurich.

Acknowledgments

The authors thank Dr. Alexander E. Ehret (Empa, Dübendorf) for providing support with continuum mechanics related formulation.

Conflict of interest

MB was employed by Altair Engineering GmbH.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2022.1034441/full#supplementary-material

Footnotes

¹This property exploits the identical nature of collagen fibers along M₁ and M₂.

²For incompressible materials, m ⋅n = 0 and hence Eqs 9, 10 must be adjusted accordingly.

³One order of magnitude smaller than in Graff (2012), Freed et al. (2005), Wu et al. (2020).

References

3D Slicer (2021). Image computing platform. Available at: https://www.slicer.org/ (Accessed November 28, 2022).

Google Scholar

Adams, M. A., and Roughley, P. J. (2006). What is intervertebral disc degeneration, and what causes it? Spine 31, 2151–2161.

PubMed Abstract | CrossRef Full Text | Google Scholar

Agur, A. M., and Dalley, A. F. (2009). Grant’s atlas of anatomy. Philadelphia, Pennsylvania, USA: Lippincott Williams & Wilkins.

Google Scholar

Ahn, B., and Kim, J. (2010). Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations. Med. Image Anal. 14, 138–148. doi:10.1016/j.media.2009.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Ambard, D., and Cherblanc, F. (2009). Mechanical behavior of annulus fibrosus: A microstructural model of fibers reorientation. Ann. Biomed. Eng. 37, 2256–2265. doi:10.1007/s10439-009-9761-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayturk, U. M., Gadomski, B., Schuldt, D., Patel, V., and Puttlitz, C. M. (2012). Modeling degenerative disk disease in the lumbar spine: A combined experimental, constitutive, and computational approach. J. Biomech. Eng. 134, 101003. doi:10.1115/1.4007632

PubMed Abstract | CrossRef Full Text | Google Scholar

Ayturk, U. M., Garcia, J. J., and Puttlitz, C. M. (2010). The micromechanical role of the annulus fibrosus components under physiological loading of the lumbar spine. J. Biomech. Eng. 132, 061007. doi:10.1115/1.4001032

PubMed Abstract | CrossRef Full Text | Google Scholar

Baroud, G., Nemes, J., Heini, P., and Steffen, T. (2003). Load shift of the intervertebral disc after a vertebroplasty: A finite-element study. Eur. Spine J. 12, 421–426. doi:10.1007/s00586-002-0512-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bathe, K.-J. (2006). Finite element procedures. Berlin, Germany: Klaus-Jurgen Bathe.

Google Scholar

Benzley, S. E., Perry, E., Merkley, K., Clark, B., and Sjaardama, G. (1995). “A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis,” in Proceedings, 4th international meshing roundtable (Citeseer), Vol. 17, 179–191.

Google Scholar

Bertagnoli, R. (2011). Bewegungserhaltende wirbelsäulenchirurgie. Amsterdam, Netherlands: Elsevier Health Sciences.

Google Scholar

Boccaccio, A., Vena, P., Gastaldi, D., Franzoso, G., Pietrabissa, R., and Pappalettere, C. (2008). Finite element analysis of cancellous bone failure in the vertebral body of healthy and osteoporotic subjects. Proc. Inst. Mech. Eng. H. 222, 1023–1036. doi:10.1243/09544119jeim296

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M., Kruse, R., and Ehret, A. E. (2013). On a staggered iFEM approach to account for friction in compression testing of soft materials. J. Mech. Behav. Biomed. Mater. 27, 204–213. doi:10.1016/j.jmbbm.2013.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bredow, J., Boese, C., Werner, C., Siewe, J., Löhrer, L., Zarghooni, K., et al. (2016). Predictive validity of preoperative ct scans and the risk of pedicle screw loosening in spinal surgery. Arch. Orthop. Trauma Surg. 136, 1063–1067. doi:10.1007/s00402-016-2487-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Brittanica, T. E. (2014). “Vertebral column,” in Encyclopaedia britannica. Editor T. E. Britannica (Chicago, IL, USA: Encyclopaedia Britannica, Inc.).

Google Scholar

Campbell, J., Coombs, D., Rao, M., Rullkoetter, P. J., and Petrella, A. (2016). Automated finite element meshing of the lumbar spine: Verification and validation with 18 specimen-specific models. J. biomechanics 49, 2669–2676. doi:10.1016/j.jbiomech.2016.05.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Caprara, S., Carrillo, F., Snedeker, J. G., Farshad, M., and Senteler, M. (2021). Automated pipeline to generate anatomically accurate patient-specific biomechanical models of healthy and pathological fsus. Front. Bioeng. Biotechnol. 9, 636953. doi:10.3389/fbioe.2021.636953

PubMed Abstract | CrossRef Full Text | Google Scholar

Cassidy, J., Hiltner, A., and Baer, E. (1989). Hierarchical structure of the intervertebral disc. Connect. tissue Res. 23, 75–88. doi:10.3109/03008208909103905

PubMed Abstract | CrossRef Full Text | Google Scholar

Cegoñino, J., Moramarco, V., Calvo-Echenique, A., Pappalettere, C., and Pérez Del Palomar, A. (2014). A constitutive model for the annulus of human intervertebral disc: Implications for developing a degeneration model and its influence on lumbar spine functioning. J. Appl. Math. 2014, 1–15. doi:10.1155/2014/658719

CrossRef Full Text | Google Scholar

Cornaz, F., Widmer, J., Farshad-Amacker, N. A., Spirig, J. M., Snedeker, J. G., and Farshad, M. (2021). Biomechanical contributions of spinal structures with different degrees of disc degeneration. Spine 46, E869–E877. doi:10.1097/brs.0000000000003883

PubMed Abstract | CrossRef Full Text | Google Scholar

del Palomar, A. P., Calvo, B., and Doblaré, M. (2008). An accurate finite element model of the cervical spine under quasi-static loading. J. biomechanics 41, 523–531. doi:10.1016/j.jbiomech.2007.10.012

CrossRef Full Text | Google Scholar

Dreischarf, M., Zander, T., Shirazi-Adl, A., Puttlitz, C., Adam, C., Chen, C., et al. (2014). Comparison of eight published static finite element models of the intact lumbar spine: Predictive power of models improves when combined together. J. biomechanics 47, 1757–1766. doi:10.1016/j.jbiomech.2014.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Dvořák, J., Panjabi, M., Chang, D., Theiler, R., and Grob, D. (1991). Functional radiographic diagnosis of the lumbar spine: Flexion–extension and lateral bending. Spine 16, 562–571. doi:10.1097/00007632-199105000-00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebara, S., Iatridis, J. C., Setton, L. A., Foster, R. J., Mow, V. C., and Weidenbaum, M. (1996). Tensile properties of nondegenerate human lumbar anulus fibrosus. Spine 21, 452–461. doi:10.1097/00007632-199602150-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Eberlein, R., Holzapfel, G. A., and Fröhlich, M. (2004). Multi-segment fea of the human lumbar spine including the heterogeneity of the annulus fibrosus. Comput. Mech. 34, 147–163. doi:10.1007/s00466-004-0563-3

CrossRef Full Text | Google Scholar

Eberlein, R., Holzapfel, G. A., and Schulze-Bauer, C. A. (2001). An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies. Comput. methods biomechanics Biomed. Eng. 4, 209–229. doi:10.1080/10255840108908005

CrossRef Full Text | Google Scholar

Fasser, M.-R., Gerber, G., Passaplan, C., Cornaz, F., Snedeker, J. G., Farshad, M., et al. (2022). Computational model predicts risk of spinal screw loosening in patients. Eur. Spine J. 31, 2639–2649. doi:10.1007/s00586-022-07187-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedorov, A., Beichel, R., Kalpathy-Cramer, J., Finet, J., Fillion-Robin, J.-C., Pujol, S., et al. (2012). 3D Slicer as an image computing platform for the quantitative imaging network. Magn. Reson. imaging 30, 1323–1341. doi:10.1016/j.mri.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Flory, P. (1961). Thermodynamic relations for high elastic materials. Trans. Faraday Soc. 57, 829–838. doi:10.1039/tf9615700829

CrossRef Full Text | Google Scholar

Freed, A. D., Einstein, D., and Vesely, I. (2005). Invariant formulation for dispersed transverse isotropy in aortic heart valves: An efficient means for modeling fiber splay. Biomech. Model. Mechanobiol. 4, 100–117. doi:10.1007/s10237-005-0069-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Frobin, W., Brinckmann, P., Biggemann, M., Tillotson, M., and Burton, K. (1997). Precision measurement of disc height, vertebral height and sagittal plane displacement from lateral radiographic views of the lumbar spine. Clin. Biomech. 12, S1–S63. doi:10.1016/s0268-0033(96)00067-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gasser, T. C., Ogden, R. W., and Holzapfel, G. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface 3, 15–35. doi:10.1098/rsif.2005.0073

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghezelbash, F., Eskandari, A. H., Shirazi-Adl, A., Kazempour, M., Tavakoli, J., Baghani, M., et al. (2021). Modeling of human intervertebral disc annulus fibrosus with complex multi-fiber networks. Acta Biomater. 123, 208–221. doi:10.1016/j.actbio.2020.12.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Goel, V. K., Kong, W., Han, J. S., Weinstein, J. N., and Gilbertson, L. G. (1993). A combined finite element and optimization investigation of lumbar spine mechanics with and without muscles. Spine 18, 1531–1541. doi:10.1097/00007632-199318110-00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Goel, V., Monroe, B., Gilbertson, L., and Brinckmann, P. (1995). Interlaminar shear stresses and laminae separation in a disc: Finite element analysis of the l3-l4 motion segment subjected to axial compressive loads. Spine 20, 689–698. doi:10.1097/00007632-199503150-00010

PubMed Abstract | CrossRef Full Text | Google Scholar

Graff, K. F. (2012). Wave motion in elastic solids. Chelmsford, MA, USA: Courier Corporation.

Google Scholar

Harkness, R. (1961). Biological functions of collagen. Biol. Rev. 36, 399–455. doi:10.1111/j.1469-185x.1961.tb01596.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashizume, H. (1980). Three-dimensional architecture and development of lumber intervertebral discs. Acta Med. Okayama 34, 301–314. doi:10.18926/AMO/30545

PubMed Abstract | CrossRef Full Text | Google Scholar

Haut, R. C., and Little, R. W. (1972). A constitutive equation for collagen fibers. J. biomechanics 5, 423–430. doi:10.1016/0021-9290(72)90001-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A., Gasser, T. C., and Ogden, R. W. (2000). A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast. 61, 1–48. doi:10.1023/a:1010835316564

CrossRef Full Text | Google Scholar

Holzapfel, G. A., Niestrawska, J. A., Ogden, R. W., Reinisch, A. J., and Schriefl, A. J. (2015). Modelling non-symmetric collagen fibre dispersion in arterial walls. J. R. Soc. Interface 12, 20150188. doi:10.1098/rsif.2015.0188

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A. (2002). Nonlinear solid mechanics: A continuum approach for engineering science. Meccanica 37, 489–490. doi:10.1023/a:1020843529530

CrossRef Full Text | Google Scholar

Holzapfel, G. A., Schulze-Bauer, C. A., Feigl, G., and Regitnig, P. (2005). Single lamellar mechanics of the human lumbar anulus fibrosus. Biomech. Model. Mechanobiol. 3, 125–140. doi:10.1007/s10237-004-0053-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A., and Weizsäcker, H. W. (1998). Biomechanical behavior of the arterial wall and its numerical characterization. Comput. Biol. Med. 28, 377–392. doi:10.1016/s0010-4825(98)00022-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, T. J. (2012). The finite element method: Linear static and dynamic finite element analysis. Chelmsford, MA, USA: Courier Corporation.

Google Scholar

Humzah, M., and Soames, R. (1988). Human intervertebral disc: Structure and function. Anat. Rec. Hob. 220, 337–356. doi:10.1002/ar.1092200402

PubMed Abstract | CrossRef Full Text | Google Scholar

HyperMesh (2017). 2017.2. 1820 E. Big Beaver Rd., Troy, MI 48083, USA: Altair Inc.

Google Scholar

Inoue, H. (1981). Three-dimensional architecture of lumbar intervertebral discs. Spine 6, 139–146. doi:10.1097/00007632-198103000-00006

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaramillo, H., Gomez, L., and García, J. J. (2015). A finite element model of the l4-l5-s1 human spine segment including the heterogeneity and anisotropy of the discs. Acta Bioeng. Biomech. 17, 15–24.

PubMed Abstract | Google Scholar

Joldes, G. R., Wittek, A., and Miller, K. (2009). Non-locking tetrahedral finite element for surgical simulation. Commun. Numer. Methods Eng. 25, 827–836. doi:10.1002/cnm.1185

PubMed Abstract | CrossRef Full Text | Google Scholar

Kibler, W. B., Press, J., and Sciascia, A. (2006). The role of core stability in athletic function. Sports Med. 36, 189–198. doi:10.2165/00007256-200636030-00001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S. (1990). Cumulative load as a risk factor for back pain. Spine 15, 1311–1316. doi:10.1097/00007632-199012000-00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuravi, R., Leichsenring, K., Trostorf, R., Morales-Orcajo, E., Böl, M., and Ehret, A. E. (2021). Predicting muscle tissue response from calibrated component models and histology-based finite element models. J. Mech. Behav. Biomed. Mater. 117, 104375. doi:10.1016/j.jmbbm.2021.104375

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. M., Hutton, W. C., and Gharpuray, V. M. (1996). Do bending, twisting, and diurnal fluid changes in the disc affect the propensity to prolapse? A viscoelastic finite element model. Spine 21, 2570–2579. doi:10.1097/00007632-199611150-00006

PubMed Abstract | CrossRef Full Text | Google Scholar

Malandrino, A., Noailly, J., and Lacroix, D. (2013). Regional annulus fibre orientations used as a tool for the calibration of lumbar intervertebral disc finite element models. Comput. methods biomechanics Biomed. Eng. 16, 923–928. doi:10.1080/10255842.2011.644539

PubMed Abstract | CrossRef Full Text | Google Scholar

Manchikanti, L. (2000). Epidemiology of low back pain. Pain physician 3, 167–192. doi:10.36076/ppj.2000/3/167

PubMed Abstract | CrossRef Full Text | Google Scholar

Maquer, G., Laurent, M., Brandejsky, V., Pretterklieber, M. L., and Zysset, P. K. (2014). Finite element based nonlinear normalization of human lumbar intervertebral disc stiffness to account for its morphology. J. Biomech. Eng. 136, 061003. doi:10.1115/1.4027300

PubMed Abstract | CrossRef Full Text | Google Scholar

Maquer, G., Schwiedrzik, J., Huber, G., Morlock, M. M., and Zysset, P. K. (2015). Compressive strength of elderly vertebrae is reduced by disc degeneration and additional flexion. J. Mech. Behav. Biomed. Mater. 42, 54–66. doi:10.1016/j.jmbbm.2014.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcolongo, M., Sarkar, S., and Ganesh, N. (2017). “7.11 trends in materials for spine surgery,” in Comprehensive biomaterials II. Editor P. Ducheyne (Oxford: Elsevier), 175–198.

CrossRef Full Text | Google Scholar

Markolf, K. L., and Morris, J. M. (1974). The structural components of the intervertebral disc: A study of their contributions to the ability of the disc to withstand compressive forces. J. Bone Jt. Surg. 56, 675–687. doi:10.2106/00004623-197456040-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Moramarco, V., Del Palomar, A. P., Pappalettere, C., and Doblaré, M. (2010). An accurate validation of a computational model of a human lumbosacral segment. J. biomechanics 43, 334–342. doi:10.1016/j.jbiomech.2009.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Noailly, J., Lacroix, D., and Planell, J. A. (2005). Finite element study of a novel intervertebral disc substitute. Spine 30, 2257–2264. doi:10.1097/01.brs.0000182319.81795.72

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Connell, G. D., Vresilovic, E. J., and Elliott, D. M. (2007). Comparison of animals used in disc research to human lumbar disc geometry. Spine 32, 328–333. doi:10.1097/01.brs.0000253961.40910.c1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogden, R. (1978). Nearly isochoric elastic deformations: Application to rubberlike solids. J. Mech. Phys. Solids 26, 37–57. doi:10.1016/0022-5096(78)90012-1

CrossRef Full Text | Google Scholar

Ogden, R. W. (1997). Non-linear elastic deformations. Chelmsford, MA, USA: Courier Corporation.

Google Scholar

Pagani, M., Reese, S., and Perego, U. (2014). Computationally efficient explicit nonlinear analyses using reduced integration-based solid-shell finite elements. Comput. Methods Appl. Mech. Eng. 268, 141–159. doi:10.1016/j.cma.2013.09.005

CrossRef Full Text | Google Scholar

Pearcy, M., Portek, I., and Shepherd, J. (1984). Three-dimensional x-ray analysis of normal movement in the lumbar spine. Spine 9, 294–297. doi:10.1097/00007632-198404000-00013

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickering, E., Pivonka, P., and Little, J. P. (2021). Toward patient specific models of pediatric ivds: A parametric study of ivd mechanical properties. Front. Bioeng. Biotechnol. 9, 632408. doi:10.3389/fbioe.2021.632408

PubMed Abstract | CrossRef Full Text | Google Scholar

Polikeit, A., Nolte, L. P., and Ferguson, S. J. (2003). The effect of cement augmentation on the load transfer in an osteoporotic functional spinal unit: Finite-element analysis. Spine 28, 991–996. doi:10.1097/01.brs.0000061987.71624.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Pooni, J., Hukins, D., Harris, P., Hilton, R., and Davies, K. (1986). Comparison of the structure of human intervertebral discs in the cervical, thoracic and lumbar regions of the spine. Surg. Radiol. Anat. 8, 175–182. doi:10.1007/bf02427846

PubMed Abstract | CrossRef Full Text | Google Scholar

Radioss (2019). Reference guide. 1820 E. Big Beaver Rd., Troy, MI 48083, USA: Altair Inc.

Google Scholar

Rajagopal, K. (2015). Remarks on the notion of “pressure”. Int. J. non-linear Mech. 71, 165–172. doi:10.1016/j.ijnonlinmec.2014.11.031

CrossRef Full Text | Google Scholar

Roberts, S., Evans, H., Trivedi, J., and Menage, J. (2006). Histology and pathology of the human intervertebral disc. J. Bone Jt. Surg. 88, 10–14. doi:10.2106/jbjs.f.00019

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohlmann, A., Burra, N. K., Zander, T., and Bergmann, G. (2007). Comparison of the effects of bilateral posterior dynamic and rigid fixation devices on the loads in the lumbar spine: A finite element analysis. Eur. Spine J. 16, 1223–1231. doi:10.1007/s00586-006-0292-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Galbusera, F., Rohlmann, A., and Shirazi-Adl, A. (2013). What have we learned from finite element model studies of lumbar intervertebral discs in the past four decades? J. biomechanics 46, 2342–2355. doi:10.1016/j.jbiomech.2013.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Heuer, F., Simon, U., Kettler, A., Rohlmann, A., Claes, L., et al. (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. Clin. Biomech. 21, 337–344. doi:10.1016/j.clinbiomech.2005.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Kettler, A., Heuer, F., Simon, U., Claes, L., and Wilke, H.-J. (2007a). Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. Spine 32, 748–755. doi:10.1097/01.brs.0000259059.90430.c2

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Kettler, A., Rohlmann, A., Claes, L., and Wilke, H.-J. (2007b). The risk of disc prolapses with complex loading in different degrees of disc degeneration–a finite element analysis. Clin. Biomech. 22, 988–998. doi:10.1016/j.clinbiomech.2007.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, T., Hu, Y., Gao, X., Dumas, J., Zorin, D., and Panozzo, D. (2019). A large scale comparison of tetrahedral and hexahedral elements for finite element analysis. arXiv preprint arXiv:1903.09332.

Google Scholar

Schultz, A., Andersson, G., Ortengren, R., Haderspeck, K., and Nachemson, A. (1982). Loads on the lumbar spine. validation of a biomechanical analysis by measurements of intradiscal pressures and myoelectric signals. J. Bone Jt. Surg. 64, 713–720. doi:10.2106/00004623-198264050-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, A. B., and Andersson, G. B. (1981). Analysis of loads on the lumbar spine. Spine 6, 76–82. doi:10.1097/00007632-198101000-00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Simo, J. C., and Hughes, T. J. (2006). Computational inelasticity, vol. 7. Berlin, Germany: Springer Science & Business Media.

Google Scholar

Smit, T. H., Odgaard, A., and Schneider, E. (1997). Structure and function of vertebral trabecular bone. Spine 22, 2823–2833. doi:10.1097/00007632-199712150-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Spencer, A. J. M. (1984). “Constitutive theory for strongly anisotropic solids, continuum theory of the mechanics of fibre-reinforced composites,” in Continuum theory of the mechanics of fibre-reinforced composites. International centre for mechanical sciences (courses and lectures). Editor A. J. M. Spencer (Vienna: Springer), Vol. 282, 1–32.

CrossRef Full Text | Google Scholar

Steinmann, J. C., and Herkowitz, H. N. (1992). Pseudarthrosis of the spine. Clin. Orthop. Relat. Res. 284, 80–90. doi:10.1097/00003086-199211000-00011

CrossRef Full Text | Google Scholar

Talukdar, R. G., Mukhopadhyay, K. K., Dhara, S., and Gupta, S. (2021). Numerical analysis of the mechanical behaviour of intact and implanted lumbar functional spinal units: Effects of loading and boundary conditions. Proc. Inst. Mech. Eng. H. 235, 792–804. doi:10.1177/09544119211008343

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiese, M. S., Hegmann, K. T., Wood, E. M., Garg, A., Moore, J. S., Kapellusch, J., et al. (2014). Prevalence of low back pain by anatomic location and intensity in an occupational population. BMC Musculoskelet. Disord. 15, 283. doi:10.1186/1471-2474-15-283

PubMed Abstract | CrossRef Full Text | Google Scholar

Truesdell, C., and Noll, W. (2004). The non-linear field theories of mechanics. 3rd edn. Berlin, Germany: Springer. doi:10.1007/978-3-662-10388-3

CrossRef Full Text | Google Scholar

Ueno, K., and Liu, Y. K. (1987). A three-dimenstional nonlinear finite element model of lumbar intervertebral joint in torsion. J. Biomech. Eng. 109, 200. doi:10.1115/1.3138670

PubMed Abstract | CrossRef Full Text | Google Scholar

Urban, J. P., and Roberts, S. (2003). Degeneration of the intervertebral disc. Arthritis Res. Ther. 5, 120–130. doi:10.1186/ar629

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshalovich, D. A., Moskalev, A. N., and Khersonskii, V. K. (1988). “Description of rotation in terms of the euler angles,” in Quantum theory of angular momentum (Singapore: World Scientific), 21–23.

Google Scholar

Wagner, D. R., and Lotz, J. C. (2004). Theoretical model and experimental results for the nonlinear elastic behavior of human annulus fibrosus. J. Orthop. Res. 22, 901–909. doi:10.1016/j.orthres.2003.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Widmer, J., Cornaz, F., Scheibler, G., Spirig, J. M., Snedeker, J. G., and Farshad, M. (2020). Biomechanical contribution of spinal structures to stability of the lumbar spine–novel biomechanical insights. Spine J. 20, 1705–1716. doi:10.1016/j.spinee.2020.05.541

PubMed Abstract | CrossRef Full Text | Google Scholar

Widmer, J., Fornaciari, P., Senteler, M., Roth, T., Snedeker, J. G., and Farshad, M. (2019). Kinematics of the spine under healthy and degenerative conditions: A systematic review. Ann. Biomed. Eng. 47, 1491–1522. doi:10.1007/s10439-019-02252-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Widmer, J. (2020). Patient specific material mapping for functional outcome prediction and planning in spinal surgery. Ph.D. thesis (Zürich, Switzerland: ETH Zurich).

Wright, S., and Nocedal, J. (2006). Numerical optimization. New York, NY: Springer.

Google Scholar

Wu, S.-W., Jiang, C., Jiang, C., and Liu, G.-R. (2020). A selective smoothed finite element method with visco-hyperelastic constitutive model for analysis of biomechanical responses of brain tissues. Int. J. Numer. Methods Eng. 121, 5123–5149. doi:10.1002/nme.6515

CrossRef Full Text | Google Scholar

Yang, B., and O’Connell, G. D. (2019). Intervertebral disc swelling maintains strain homeostasis throughout the annulus fibrosus: A finite element analysis of healthy and degenerated discs. Acta biomater. 100, 61–74. doi:10.1016/j.actbio.2019.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Zander, T., Dreischarf, M., Timm, A.-K., Baumann, W. W., and Schmidt, H. (2017). Impact of material and morphological parameters on the mechanical response of the lumbar spine–a finite element sensitivity study. J. biomechanics 53, 185–190. doi:10.1016/j.jbiomech.2016.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zander, T., Rohlmann, A., and Bergmann, G. (2009). Influence of different artificial disc kinematics on spine biomechanics. Clin. Biomech. 24, 135–142. doi:10.1016/j.clinbiomech.2008.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar