Abstract

The intervertebral disc is a complex structure that experiences multiaxial stresses regularly. Disc failure through herniation is a common cause of lower back pain, which causes reduced mobility and debilitating pain, resulting in heavy socioeconomic burdens. Unfortunately, herniation etiology is not well understood, partially due to challenges in replicating herniation in vitro. Previous studies suggest that flexion elevated risks of herniation. Thus, the objective of this study was to use a multiscale and multiphasic finite element model to evaluate the risk of failure under torque- or muscle-driven flexion. Models were developed to represent torque-driven flexion with the instantaneous center of rotation (ICR) located on the disc, and the more physiologically representative muscle-driven flexion with the ICR located anterior of the disc. Model predictions highlighted disparate disc mechanics regarding bulk deformation, stress-bearing mechanisms, and intradiscal stress–strain distributions. Specifically, failure was predicted to initiate at the bone-disc boundary under torque-driven flexion, which may explain why endplate junction failure, instead of herniation, has been the more common failure mode observed in vitro. By contrast, failure was predicted to initiate in the posterolateral annulus fibrosus under muscle-driven flexion, resulting in consistent herniation. Our findings also suggested that muscle-driven flexion combined with axial compression could be sufficient for provoking herniation in vitro and in silico. In conclusion, this study provided a computational framework for designing in vitro testing protocols that can advance the assessment of disc failure behavior and the performance of engineered disc implants.

1 Introduction

Mechanical failure of the intervertebral disc, such as lumbar disc herniation, is a common cause of lower back pain, which can affect 10% of the population annually [1]. Specifically, disc herniation can cause decreased mobility and debilitating pain and has been the principal cause for working-age individuals to undergo spinal surgeries, resulting in significant socioeconomic burdens [2,3]. Since the early 20th century, lumbar disc herniation has been the focus of spinal biomechanical and clinical research [4]; however, despite significant developments in joint-level testing techniques, challenges remain in repeatably inducing herniation in vitro, largely due to difficulties in replicating the multiaxial loads that the disc experience during physiological activities [5].

The range of viable lumbar disc in vitro mechanical tests and their clinical relevance are often limited by the capabilities of available testing equipment. For example, in vivo flexion and extension motions are mainly driven by active physiological structures (e.g., muscles), causing the instantaneous center of rotation (ICR) to be located at some distance away from the disc (i.e., muscle-driven) [6]. However, in vitro flexion or extension testing has been primarily conducted with the ICR located on the disc [5]. Nonphysiological torque-driven bending tests could contribute to the limited success in provoking in vitro herniation, making it more challenging for researchers to study the etiology and progression of disc herniation. For example, Adams and Hutton attempted to induce herniation by loading joint-level specimens under axial compression with a flexion angle. However, over 70% of the samples experienced non-herniation failure, with endplate junction failure being the most common failure mode observed [79]. More recently, six-degree-of-freedom loading devices highlighted the benefit and necessity of applying combined multiaxial loads to induce herniation. However, ∼50% of the motion segment specimens were still excluded due to endplate failure [5,10].

Finite element models have been an effective tool to complement experimental studies, providing predictions of mechanical behavior that are difficult or impossible to measure in the laboratory. Researchers have used various models of the intervertebral discs to investigate joint- and tissue-level stress and strain distributions. However, many of these models rely on single-phasic or poroelastic material descriptions that are not capable of describing Donnan equilibrium [1121]. The Donnan equilibrium is largely responsible for the swelling behavior observed in biological tissues and plays a pivotal role in tissue mechanics [22].

To address these limitations, we recently developed and validated a novel structure-based triphasic model for the bovine caudal motion segment. Model parameters were determined based on known physical or biochemical properties reported in the literature (e.g., collagen fiber stiffness and fixed charge density) [23]. The model also explicitly described individual annulus fibrosus collagen fiber bundles so subtissue-level mechanics could be directly investigated. The model accurately predicted disc mechanics across the joint, tissue, and subtissue scales and helped elucidate important load-bearing mechanisms between tissue phases (e.g., between fluid and solid phases) and tissue subcomponents (e.g., between fibers and extrafibrillar matrix), which are essential for understanding and assessing disc failure under multiaxial loading.

In this study, we employed our validated multiscale multiphasic structure-based model to address the current challenges in replicating physiological loading and disc failure in vitro. The objective of this study was twofold. The first objective was to investigate disc mechanics under torque- and muscle-driven flexion. The second objective was to relate those findings to clinical and experimental observations of disc failure. Although the current model was developed from bovine caudal disc geometry, the model findings corresponded well to clinical observations for bulging or herniated discs under muscle-driven flexion. Model predictions under torque-driven flexion also highlighted strain concentrations that corresponded well to disc failure behavior commonly observed in vitro for human and ovine discs. Thus, the findings from this study are considered translatable to relevant human disc biomechanics research.

2 Methods

2.1 Model Development.

Finite element models of the bovine caudal disc motion segment were created based on our previous work (Fig. 1(a)) [23]. To replicate motion segment samples prepared for most in vitro experiments, posterior structures, including facets joints and ligaments, were not included. The model geometry was created in Solidworks (Solidworks 2020). Finite element meshes were generated using abaqus and ansa preprocessor (abaqus 6.14; ansa 15.2.0). The appropriate mesh size was determined using results from our previous mesh convergence study [24]. Boundary and loading conditions were defined in FEBioStudio, and the fully developed models were solved by febio (febiostudio 1.5) [25]. Our prior work validated that proportional scaling did not significantly alter model predictions [23]. Thus, the disc joint was modeled at a 1:5 scale (∼2.1 × 106 tetrahedral elements) for computational efficiency due to limited accessible computing power (maximum of ∼200 × 106 nonzero entries in the stiffness matrix can be evaluated).

Fig. 1
(a) Schematics of the multiscale bovine caudal disc motion segment model. The top inset shows the cartilage endplate geometry. The bottom right inset details the angle-ply AF fiber structure (θ: fiber angle). This subfigure is adapted from Zhou et al. [23]. Radial variation of (b) AF fiber angle and solid volume fraction, and (c) tissue fixed charge density.
Fig. 1
(a) Schematics of the multiscale bovine caudal disc motion segment model. The top inset shows the cartilage endplate geometry. The bottom right inset details the angle-ply AF fiber structure (θ: fiber angle). This subfigure is adapted from Zhou et al. [23]. Radial variation of (b) AF fiber angle and solid volume fraction, and (c) tissue fixed charge density.
Close modal

Model geometry was determined based on data in the literature. A circular cross section was assumed in the transverse anatomical plane. Disc radius and height (not including bony endplates) were 2.85 and 1.40 mm, respectively (Fig. 1(a)). The nucleus pulposus (NP) was assumed to have the same circular cross section, with a ∼50% smaller radius (1.45 mm; Fig. 1(a)) [26]. The annulus fibrosus (AF) was created using the multiscale structure-based modeling approach validated in our previous work [23,27]. Particularly, seven concentric 0.2 mm-thick lamellae were modeled as fiber-reinforced angle-ply composites containing distinct materials for fiber bundles and the extrafibrillar matrix that occupy separate volumes (Fig. 1(a)) [28]. Native bovine AF structural features, including lamellar thickness, fiber bundle radius, and interfibrillar spacing, were maintained during the downscale to reduce the total number of elements required. This scaling approach has been widely applied and validated in human disc finite element models [11,18,19,29] and has been shown to improve computational efficiency while maintaining fiber volume fraction and preserving mesh quality for model convergence [24]. AF fiber bundles were modeled as uniformly distributed, full-length cylinders welded to the surrounding matrix [2931]. Although available data regarding bovine AF structures are limited in the literature, similarities between human and bovine AF structures have been reported [32]. Therefore, fiber bundle geometry from the human AF was applied, where fiber bundle radius was 0.06 mm and interfibrillar spacing within each lamella was 0.22 mm [28,33]. Fibers were oriented at ±45 deg to the transverse plane in the inner AF and decreased along the radial direction to ±30 deg in the outer AF (Fig. 1(a)—bottom right inset; Fig. 1(b)—turquoise circles) [34]. Cartilage endplates (CEP) covered the superior and inferior ends of the NP and the inner-to-mid-AF (Fig. 1(a)—cartilage endplate); spatial variations in CEP thickness were also incorporated (Fig. 1(a)—top inset) [35]. Bony endplates were modeled to cover both the superior and inferior ends of the disc (Fig. 1(a)—bony endplate). All interfaces were defined as welded interfaces [33]. To exclude the effect of mesh size on model-predicted mechanics, element size was held constant (Figure 1 available in the Supplemental Materials on the ASME Digital Collection).

Triphasic mixture theory was employed to account for Donnan equilibrium to properly describe the tissue water content and osmotic response [36]. The Holmes–Mow description was applied to model the strain-dependent tissue permeability (k) of the NP, AF, and CEP (Eq. (1)), where J was the determinant of the deformation gradient tensor (F), k0 represented hydraulic permeability in the reference configuration, φ0 represented tissue solid volume fraction, and M represented the exponential strain-dependence coefficient
k(J)=k0(Jφ01φ0)2e12M(J21)
(1)

Fixed charge density represented the proteoglycan content in the NP, CEP, and AF extrafibrillar matrix, driving the osmotic response. Radial variation in fixed charge density and AF solid volume fraction were determined based on previous literature and our recent work, where high-spatial-resolution measurements of bovine caudal disc biochemical composition were provided (Fig. 1(b)—grayscale circles; Fig. 1(c)) [37,38]. Collagen fiber bundles were assumed to have no active swelling capacity (i.e., zero fixed charge density). Free diffusivity (D0) and within-tissue diffusivity (D) terms for Na+ and Cl ions for the simulated phosphate-buffered saline solution were set based on data from Gu et al. [39] with a 100% ion solubility assumed (D0,Na+=0.00116 mm2/s; D0,Cl=0.00161 mm2/s; DNa+=0.00044 mm2/s; DCl=0.00069 mm2/s). The solution osmotic coefficient (0.927) was determined based on a linear interpolation of data reported in Partanen et al. [40].

A compressible hyperelastic Holmes–Mow material description was used to describe NP, CEP, and AF extrafibrillar matrix mechanics (Eqs. (2)(4)) [41]. In the equations, I1 and I2 represented the first and second invariants of the right Cauchy–Green deformation tensor, C(C=FTF), E represented Young's modulus, ν represented Poisson's ratio, and β represented the exponential stiffening coefficient. AF collagen fibers were modeled using the same compressible hyperelastic Holmes-Mow as the ground matrix but were reinforced with a tension-only power-linear fiber description to account for AF nonlinearity and anisotropy (Eq. (5)). The AF collagen fibers were individually modeled instead of being embedded into the fiber-reinforced constitutive relationship, such that AF subtissue-level mechanics, such as interfibrillar stress and strain distributions, can be explicitly investigated. In Eq. (5), γ represented the power-law exponent in the toe region, Elin represented the fiber modulus in the linear region, and λ0 represented the transitional stretch between the toe and linear regions. Additionally, B was a function of γ, Elin, and λ0 (B=Elin2((λ021)2(γ1)+λ02). Collagen fiber properties were determined based on uniaxial tensile data for type I collagen [42,43]
W(I1,I2,J)=12c(eQ1)
(2)
Q=β(1+ν)(12ν)E(1ν)[(E1+νEν(1+ν)(12ν))(I13)+Eν(1+ν)(12ν)(I23)(E1+ν+Eν(1+ν)(12ν))lnJ2]
(3)
c=E(1ν)2β(1+ν)(12ν)
(4)
ψn(λn)=0λn<1Elin4γ(γ1)(λ021)2γ(λn1)γ1λnλ0Elin(λnλ0)+B(λn2λ02)+Elin4γ(γ1)(λ021)2γ(λn1)γλn>λ0
(5)
Bony endplates were modeled using a compressible hyperelastic material with the Neo-Hookean description (Eq. (6)). I1, I2, J were defined as above.Ebonyendplates and νbonyendplates represented Young's modulus (12,000 MPa) and Poisson's ratio (0.3), based on reported data in the literature [44,45]
Wbonyendplates(I1,I2,J)=Ebonyendplates4(1+νbonyendplates)(I13)Ebonyendplates2(1+νbonyendplates)lnJ+Ebonyendplatesνbonyendplates(1+νbonyendplates)(12νbonyendplates)(lnJ)2
(6)
All model parameters were directly obtained from our previous work that developed and validated the model for the bovine caudal disc motion segment (Table 1) [23]. Bovine tissue properties were used when corresponding data were available. When bovine data were not available, matching human disc properties were used, as previous studies have shown similarities between healthy human and bovine disc mechanical and biochemical properties (Table 1—“*”) [38,49,50].
Table 1

Triphasic material properties used in the model

AF
PhaseNPMatrixFibersCEP
Fluidφ00.2aFigure 1(b) a0.4c,h
k0 × 10−16 (m4/Ns)5.5b64b64b5.6c,h
M1.92c,h4.8c,h4.8c,h3.79c,h
SolidE (MPa)0.4b0.74b0.74b0.31g
ν0.24d0.16c,h0.16c,h0.18c,h
β0.95c,h3.3c,h3.3c,h0.29c,h
Elin (MPa)N.A.N.A.600eN.A.
γN.A.N.A.5.95f,hN.A.
λ0N.A.N.A.1.05eN.A.
AF
PhaseNPMatrixFibersCEP
Fluidφ00.2aFigure 1(b) a0.4c,h
k0 × 10−16 (m4/Ns)5.5b64b64b5.6c,h
M1.92c,h4.8c,h4.8c,h3.79c,h
SolidE (MPa)0.4b0.74b0.74b0.31g
ν0.24d0.16c,h0.16c,h0.18c,h
β0.95c,h3.3c,h3.3c,h0.29c,h
Elin (MPa)N.A.N.A.600eN.A.
γN.A.N.A.5.95f,hN.A.
λ0N.A.N.A.1.05eN.A.

NP: nucleus pulposus; AF: annulus fibrosus; CEP: cartilage endplate; φ0: solid volume fraction; k0: referential hydraulic permeability; M: exponential strain-dependence coefficient for permeability; E: Young's modulus; ν: Poisson's ratio; β: exponential stiffening coefficient of the Holmes–Mow model; Elin: collagen fiber bundle linear-region modulus; γ: collagen fiber bundle toe-region power-law exponent; λ0: collagen fiber bundle toe- to linear-region transitional stretch.

a

Beckstein et al. [37].

b

Périé et al. [46].

c

Cortes et al. [41].

d

Farrell and Riches [47].

e

Van der Rijt et al. [43], Shen et al. [42].

f

Zhou et al. [27].

g

Wu et al. [48].

h

The parameter was determined based on experimental studies using matching human intervertebral disc tissues due to the lack of corresponding data obtained from bovine caudal disc tissues; N.A.: not applicable.

2.2 Loading and Boundary Conditions.

All models were loaded in three steps (Fig. 2(a)). First, free swelling in the 0.15 M phosphate-buffered solution was applied until equilibrium. Then, a 0.5 MPa of axial compression was applied, which was immediately followed by a 5-deg flexion. During axial compression and flexion, all degrees-of-freedom were fixed for the bottom bony endplate (Fig. 2(a)—fixed boundary condition). The flexion angle was determined based on human lumbar spine range of motion data [6]. Due to the symmetry in bovine caudal disc geometry, only flexion was simulated.

Fig. 2
(a) Loading schematics demonstrate the model orientation, boundary condition, and loading conditions defined in the three loading steps. (b) Schematic of model mid-frontal plane demonstrating the orientation and the inner and outer AF location. (c) Schematic of model mid-transverse plane demonstrating the orientation, the posterolateral region, the rim, and the locations where the AF bulging was evaluated.
Fig. 2
(a) Loading schematics demonstrate the model orientation, boundary condition, and loading conditions defined in the three loading steps. (b) Schematic of model mid-frontal plane demonstrating the orientation and the inner and outer AF location. (c) Schematic of model mid-transverse plane demonstrating the orientation, the posterolateral region, the rim, and the locations where the AF bulging was evaluated.
Close modal

To model torque-driven flexion, the ICR was located on the line of symmetry on the top bony endplate (not including the edge; Fig. 3(a)). For muscle-driven flexion, ICRs were located on the same line of symmetry but at some distance anterior of the disc edge (Fig. 3(b)). A total of 10 cases were investigated, where cases A to C were considered torque-driven, and cases D to J were considered muscle-driven (Fig. 3(c)). The distance between the center of the top bony endplate and the ICR was defined as ICR distance (Fig. 3—ICR distance). To examine the effect of flexion on initiating disc herniation, axial rotation was not included to avoid potential confounding effects. All models were simulated using steady-state analyses and model outputs were evaluated at equilibrium.

Fig. 3
Schematics of (a) torque-driven flexion, where the instantaneous center of rotation (ICR) is located on the top bony endplate along the line of symmetry, and (b) muscle-driven flexion, where the ICR is located on the same anatomical transverse plane along the line of symmetry, but away from the disc. The distance between the center of the top bony endplate and the ICR was defined as the ICR distance. (c) The ICR location for the 10 cases investigated.
Fig. 3
Schematics of (a) torque-driven flexion, where the instantaneous center of rotation (ICR) is located on the top bony endplate along the line of symmetry, and (b) muscle-driven flexion, where the ICR is located on the same anatomical transverse plane along the line of symmetry, but away from the disc. The distance between the center of the top bony endplate and the ICR was defined as the ICR distance. (c) The ICR location for the 10 cases investigated.
Close modal

2.3 Data Analysis: Disc Mechanics Under Torque- and Muscle-Driven Flexion.

The magnitude of the torque and corresponding force required to achieve the 5-deg flexion were calculated. Force magnitudes were calculated as the torque divided by the corresponding ICR distance, which served as the lever arm. Intradiscal deformation was assessed by evaluating strains in the z-direction (Fig. 2(b)), which represented changes in disc height, and AF bulging at mid-disc height (Fig. 2(c)—solid red circles). The average disc height was calculated as the average of anterior and posterior disc height after flexion. Absolute AF bulging (i.e., AF radial displacement after flexion, [mm]) and relative AF bulging (i.e., [AF radial displacement postflexion]/[inner or outer radius of the AF ring]×100%) were evaluated using the post-swelling configuration as the reference configuration to better mimic the reference configuration used for previous experimental studies [51]. AF buckling was noted when the AF radius at mid-disc height became smaller than it was near the endplates.

Average disc solid stress (i.e., stress taken by the tissue solid phase) and fluid pressure (i.e., stress taken by the tissue fluid phase) were evaluated before and after flexion. The relative contribution of solid stress and fluid pressure was evaluated by normalizing each term by the total stress, which was defined as the sum of the two terms [36]. Effective Lagrangian strain, effective solid Lagrangian stress, fluid pressure, and maximum shear Lagrangian strain distributions were evaluated at the midfrontal plane. Within the febio environment, effective stresses and strains are comparable to Von Mises stresses and strains. Average NP fluid pressure was also evaluated.

2.4 Data Analysis: Predicting Risk of Herniation.

Model predictions of in vitro disc herniation were determined using two failure criteria based on both the AF effective strain and AF fiber stretch [16,17,52]. We evaluated the risk of in vitro herniation mainly based on AF failure mechanics due to the clinical prevalence of AF failure in the posterolateral disc region, the availability of AF failure mechanics data in the literature, and the lack of failure mechanics data characterized at the disc-bone interface. The average effective strain values were calculated in the posterolateral inner and outer AF before and after the applied flexion (Fig. 2(b)—inner and outer AF; Fig. 2(c)—posterolateral region) and were compared to the effective failure strain threshold reported in the literature [52]. Particularly, the range of effective strain that initiated failure in the AF was defined as 0.4–0.6, and the percentage of failed elements after flexion was calculated as the number of AF elements with an effective strain value above the threshold (i.e., 0.5, calculated as the average of the upper and lower bound for the failure initiation range previously defined) divided by the total number of AF elements in the respective region. Due to the consistent mesh size applied, the percentage of failed elements was considered equivalent to the failed tissue volume. The average AF fiber stretch in the posterolateral inner and outer AF with the corresponding percentage of failed elements was similarly calculated. The AF fiber stretch failure threshold (failure initiation range: 1.15–1.25; the failed element percentage was calculated using the fiber stretch value of 1.20) was determined based on values reported in previous joint-level studies [16,17,53].

AF maximum shear strain has been a commonly used metric to characterize disc mechanical response under combined loading in joint-level models [16,17,54]. To examine maximum shear strain as a candidate for a failure criterion, the average maximum shear strain in the posterolateral inner and outer AF with the corresponding failed element percentage was calculated. The threshold (failure initiation range: 0.3–0.5; the failed element percentage was calculated using the maximum shear strain value of 0.4) was determined based on values reported in previous joint-level modeling studies [16,17,54].

Although the risk of herniation in vitro was mainly evaluated using AF-based failure criteria, the average effective strain was also evaluated at the rim (i.e., outer cartilage endplate locating at the bone-AF interface) in the posterolateral region to help investigate the causes for the commonly observed endplate junction failure in vitro (Fig. 2(c)—posterolateral region; rim). The range of effective strain that initiated failure in the rim was defined as 0.5–0.7, and the percentage of failed elements was calculated using the effective strain value of 0.6 [55].

3 Results

3.1 Disc Mechanics Under Torque- and Muscle-Driven Flexion.

The torque magnitude required to achieve the 5-deg flexion increased nonlinearly as the ICR distance increased, except for case A, whose ICR was located at the center of the top bony endplate (Fig. 4). However, the corresponding force magnitude required followed a parabolic trend, reaching the minimum with an ICR distance between 3 and 6 mm, which was in the range of muscle-driven flexion (Fig. 4—muscle-driven).

Fig. 4
The torque and force magnitudes required to achieve 5 deg flexion
Fig. 4
The torque and force magnitudes required to achieve 5 deg flexion
Close modal

The average disc height was maintained under torque-driven flexion but was increased under muscle-driven flexion (Fig. 5(a)). Specifically, under torque-driven flexion, the posterior AF experienced tensile strains while the anterior AF experienced compressive strains in the z-direction (Fig. 5(a)—diagonal circles). However, under muscle-driven flexion, both the posterior and anterior sides of the disc experienced tensile z-strains (Fig. 5(a)—solid circles). Overall, both the anterior and posterior disc height increased linearly with ICR distance, resulting in large disc height differences between torque- and muscle-driven cases (Fig. 5(a)). For example, the average disc height for the case I (height: 2.2 mm) was ∼50% greater than that for case A (1.4 mm; Fig. 5(b)).

Fig. 5
(a) The anterior, posterior, and average disc height at 5 deg flexion. The horizontal dashed line highlights the average pre-flexion disc height. (b) Disc mid-frontal plane z-strain distributions for five representative cases. Anterior and posterior disc heights were labeled.
Fig. 5
(a) The anterior, posterior, and average disc height at 5 deg flexion. The horizontal dashed line highlights the average pre-flexion disc height. (b) Disc mid-frontal plane z-strain distributions for five representative cases. Anterior and posterior disc heights were labeled.
Close modal

The inner and outer radius of the AF ring were 1.56 and 2.94 mm in the reference configuration (Fig. 6(a)). Assessment of AF radial displacement at the mid-disc height suggested outward bulging for both the inner and outer AF after compression (Fig. 6(b)). Under torque-driven flexion, both the inner and outer AF bulged outward on the posterior side, while the anterior AF experienced inward bulging (Fig. 6(c)—diagonal triangles and circles; Fig. 6(d)—cases A and B). Under muscle-driven flexion, the relative inward bulging for the posterior AF and anterior outer AF increased with ICR distance (Fig. 6(c)—solid orange circles and triangles; solid blue circles); however, bulging in the anterior inner AF was relatively consistent across all muscle-driven cases, ranging from −2.3% to −3.5% (Fig. 6(c)—solid blue triangles). Buckling in the posterior outer AF was first observed in case I, when the relative inward bulging exceeded 4% (Fig. 6(d)—cases I). Buckling in the anterior AF was not observed for any cases investigated.

Fig. 6
(a) AF bulging was evaluated using the post-swelling configuration as a reference. Data in red and blue suggest outward and inward bulging compared to the reference configuration. (b) Model-predicted AF bulging after compression. Relative bulging is only shown for one side due to symmetry. (c) Relative bulging in the inner and outer AF evaluated in the posterior and anterior regions. Positive and negative relative bulging suggest outward and inward AF bulging compared to the reference configuration. The horizontal dashed line represents the relative disc bulging threshold (0%), below which the AF was predicted to bulge inward. (d) Disc midfrontal cross sections demonstrating post-flexion AF bulging for five representative cases.
Fig. 6
(a) AF bulging was evaluated using the post-swelling configuration as a reference. Data in red and blue suggest outward and inward bulging compared to the reference configuration. (b) Model-predicted AF bulging after compression. Relative bulging is only shown for one side due to symmetry. (c) Relative bulging in the inner and outer AF evaluated in the posterior and anterior regions. Positive and negative relative bulging suggest outward and inward AF bulging compared to the reference configuration. The horizontal dashed line represents the relative disc bulging threshold (0%), below which the AF was predicted to bulge inward. (d) Disc midfrontal cross sections demonstrating post-flexion AF bulging for five representative cases.
Close modal

Fluid pressure contributed significantly to the disc's overall stress-bearing capability. Before flexion was applied, the average fluid pressure was 0.24 MPa, which corresponded to 47% of the total stress (Fig. 7). As the ICR distance increased, the average fluid pressure decreased nonlinearly from 0.25 MPa in case A to 0.08 MPa in case J (Fig. 7(a)—blue bars), while the average solid stress followed a parabolic trend, reaching its minimum in case F (0.11 MPa) and then increasing with ICR distance, reaching its maximum in case J (0.43 MPa; Fig. 7(a)—black bars). As a result, under torque-driven flexion, the relative solid stress and fluid pressure contributions were comparable and not altered by the applied flexion (Fig. 7(b)—diagonal bars). However, under muscle-driven flexion, the relative fluid pressure contribution decreased pseudo-linearly with increasing ICR distance, from 55% in case D to 15% in case J (Fig. 7(b)—solid bars).

Fig. 7
Model-predicted (a) solid stress and fluid pressure, as well as (b) their relative contribution to the total stress taken by the disc pre- and post-flexion. The horizontal dashed line in (b) highlights the relative contribution before flexion.
Fig. 7
Model-predicted (a) solid stress and fluid pressure, as well as (b) their relative contribution to the total stress taken by the disc pre- and post-flexion. The horizontal dashed line in (b) highlights the relative contribution before flexion.
Close modal

Before flexion, high maximum shear strains and effective strains were observed near the rim (Figs. 8(a) and 8(b)—“*”). The solid stress was mainly absorbed by the AF (Fig. 8(c)), while the fluid pressure was concentrated in the NP, with an average NP fluid pressure of 0.47 MPa (Fig. 8(d)). Torque-driven flexion had a minimal impact on the maximum shear strain, effective strain, and fluid pressure distributions (Figs. 8(a), 8(b), and 8(d)) but resulted in greater stresses concentrated in the anterior outer AF (Fig. 8(c)). By contrast, muscle-driven flexion resulted in higher maximum shear strains, effective strains, and effective solid stresses in the posterior AF and at the disc-bone boundary (Figs. 8(a)8(c)). Additionally, the NP fluid pressure decreased pseudo-linearly with increasing ICR distance (Fig. 8(d); Fig. 8(e)—solid black circles). For example, the average NP fluid pressure for the case I was 0.19 MPa, representing a 60% decrease from case A (0.47 MPa) and the preflexion configuration (0.47 MPa). Overall, under muscle-driven flexion, changes in strain, stress, and fluid pressure increased with increasing ICR distance, with changes in effective strains being more apparent than changes in maximum shear strains.

Fig. 8
Pre- and post-flexion disc midfrontal plane (a) maximum shear strain, (b) effective strain, (c) effective solid stress, and (d) fluid pressure distributions for five representative cases. In the pre-flexion configuration, “*” in (a) and (b) highlight strain concentrations. Average nucleus pulposus (NP) fluid pressure was labeled in (d). (e) Average NP fluid pressure values at 5-deg flexion. The horizontal line highlights the average NP fluid pressure value before flexion.
Fig. 8
Pre- and post-flexion disc midfrontal plane (a) maximum shear strain, (b) effective strain, (c) effective solid stress, and (d) fluid pressure distributions for five representative cases. In the pre-flexion configuration, “*” in (a) and (b) highlight strain concentrations. Average nucleus pulposus (NP) fluid pressure was labeled in (d). (e) Average NP fluid pressure values at 5-deg flexion. The horizontal line highlights the average NP fluid pressure value before flexion.
Close modal

3.2 Predicting Risk of Herniation.

The average effective strain in the posterolateral inner and outer AF was 0.25 and 0.26 before flexion, and no elements were predicted to fail (Fig. 9(a)). Torque-driven flexion had a negligible effect on the AF effective strain, regardless of the AF location (Fig. 9(a)). By contrast, muscle-driven flexion increased effective strain in both the posterolateral inner and outer AF, with the strain magnitude increasing pseudo-linearly with ICR distance (Fig. 9(a)). Based on the effective strain criterion, the posterolateral outer AF was predicted to fail before the inner AF (Fig. 9(a)—solid versus hollow bars). Noticeably, the effective strain in the posterolateral outer AF reached 0.54 in case I, resulting in 65% of elements exceeding the failure threshold (Fig. 9(a)—red “+” symbol; Fig. 9(a) – case I, solid bar and circle). However, the effective strain in the posterolateral inner AF and the percentage of failed elements never exceeded 0.45 and 5% for all cases (Fig. 9(a)—hollow bars and circles).

Fig. 9
The average (a) effective strain and (b) fiber stretch with the corresponding failed element percentage evaluated in the posterolateral (post-lat) inner and outer AF. The boxes represent the range where tissue failure was expected to initiate. The failed element percentage was calculated using the failure threshold highlighted by the horizontal dashed lines, above which tissue failure was highly expected. The “*,” “̂,” and “+” represent failure initiation in the post-lat IAF fibers, OAF fibers, and bulk OAF.
Fig. 9
The average (a) effective strain and (b) fiber stretch with the corresponding failed element percentage evaluated in the posterolateral (post-lat) inner and outer AF. The boxes represent the range where tissue failure was expected to initiate. The failed element percentage was calculated using the failure threshold highlighted by the horizontal dashed lines, above which tissue failure was highly expected. The “*,” “̂,” and “+” represent failure initiation in the post-lat IAF fibers, OAF fibers, and bulk OAF.
Close modal

The average fiber stretch in the posterolateral inner and outer AF was 1.11 and 1.07 before flexion, and no elements were predicted to fail (Fig. 9(b)). The average fiber stretch in the inner and outer AF increased pseudo-linearly with ICR distance, and failure was predicted to occur earlier in the inner AF fibers than the outer AF fibers, regardless of the type of flexion applied (Fig. 9(b)—hollow versus solid bars). For the inner AF, the average fiber stretch was 1.21 in case D, resulting in over 60% of the elements exceeding the AF fiber stretch failure threshold (Fig. 9(b)—red “*” symbol; Fig. 9(b)—case D, hollow bar and circle). For the outer AF, the average fiber stretch did not reach the failure threshold until case H, where ∼50% of the elements were predicted to fail (Fig. 9(b)—red “̂” symbol; Fig. 9(b)—case H, solid bar and circle).

The average effective strain in the rim was 0.55 before flexion, with only 5% of elements predicted to fail (Fig. 10(a)). The average effective strain was consistent for cases A to F and then increased pseudo-linearly with increasing ICR distance, from 0.57 in case F to 1.63 in case J (Fig. 10(a)). Interestingly, the percentage of failed elements followed a parabolic trend, where more than 50% of rim elements failed in cases A and B (torque-driven), as well as in cases G to J (muscle-driven).

Fig. 10
(a) The average effective strain evaluated in the rim, and (b) the average maximum shear strain evaluated in the inner AF and outer AF in the posterolateral (postlat) region with the corresponding failed element percentage. The boxes represent the range where failure was expected to initiate. The failed element percentages were calculated using the failure threshold highlighted by the horizontal dashed lines, above which tissue failure was highly expected. The “+'s” in (a) represent cases with >50% failed element percentage.
Fig. 10
(a) The average effective strain evaluated in the rim, and (b) the average maximum shear strain evaluated in the inner AF and outer AF in the posterolateral (postlat) region with the corresponding failed element percentage. The boxes represent the range where failure was expected to initiate. The failed element percentages were calculated using the failure threshold highlighted by the horizontal dashed lines, above which tissue failure was highly expected. The “+'s” in (a) represent cases with >50% failed element percentage.
Close modal

The average maximum shear strain in the posterolateral inner and outer AF were both 0.14 before flexion, and no elements were predicted to fail (Fig. 10(b)). AF average maximum shear strain never exceeded the failure threshold, regardless of the AF anatomical region and ICR distance (Fig. 10(b)).

4 Discussion

This study used a finite element modeling approach to investigate the risk of tissue failure leading to herniation under flexion. We employed our structure-based model that was previously validated under single and combined loading conditions to evaluate multiscale disc mechanics under torque- and muscle-driven flexion. The torque-driven models are intended to replicate the commonly used in vitro flexion testing setup with the ICRs located on the disc. By contrast, physiologically representative flexion motions driven by muscle contractions result in ICRs located anterior of the disc and were simulated by the muscle-driven models [6]. The risk of herniation was assessed based on posterolateral AF failure, which was considered a major precursor for herniation. Model simulations demonstrated vastly different disc mechanics under the two flexion setups. Our findings illustrated that by shifting the instantaneous center of rotation to the anterior of the disc, the more physiologically representative muscle-driven flexion placed the disc at a higher risk for herniation through posterolateral AF failure, which is representative of clinical observations [3]. Under torque-driven flexion, strains were more concentrated in the rim. This finding helped explain the more commonly observed endplate junction failure in vitro, which contributed to the limited success researchers have had in provoking in vitro herniation in the past several decades [5,710].

Replicating herniation in vitro is essential for researchers to study herniation etiology and in vivo failure mechanisms related to mechanical overloading. Investigating disc failure is important for understanding and assessing workplace risks (e.g., factory workers) and for evaluating the performance of engineered implants [56]. Our failure criterion was defined in vitro herniation based on both bulk AF strain and AF fiber stretch in the posterolateral region. Under torque-driven flexion, neither the bulk AF strain nor the AF fiber stretch in the posterolateral region exceeded their respective failure threshold (failed elements <20%; Fig. 9—torque-driven). Thus, herniation through the posterolateral AF was not predicted for any torque-driven flexion cases; however, for cases A and B, the effective strain in the rim exceeded the failure threshold while the failed element percentage exceeded 50% (Fig. 10(a)), suggesting that torque-driven flexion most likely initiated failure from the endplate instead of the AF. These observations agree well with in vitro herniation studies, where endplate junction failure instead of herniation has been the main provoked failure mode under combined loading [3,5,710].

Applying more physiologically relevant muscle-driven flexion increased the likelihood of herniation through posterolateral AF (Figs. 8 and 9). The risk of in vitro herniation increased greatly with ICR distance. In these cases, failure was predicted to first occur in the inner AF before the outer AF fibers and bulk AF (Fig. 9(b)). These predicted failure locations were consistent with clinical observations for herniated discs [3]. Together with the predicted failure mode (i.e., endplate junction failure) under torque-driven flexion, the predictive power of our model and the failure criterion applied was demonstrated. However, caution is still needed when interpreting these results for experimental study design such that an ICR distance that may result in unwanted endplate failure under muscle-driven flexion can be avoided. For example, case H had similar effective strains, fiber stretches, and corresponding failed element percentages in the posterolateral AF as case J, but the average effective strain in the rim for case H was ∼40% smaller (Fig. 10(a)), which reduced the risk of premature rim failure and made it a more preferred option than case J.

Interstitial fluid plays a pivotal role as a stress-bearing mechanism in hydrated soft tissues, including articular cartilage, meniscus, and the intervertebral disc [23,5759]. In healthy cartilage, fluid pressurization can contribute to more than 80% of the total stress. A loss of fluid pressurization, which occurs with aging and degeneration, can lead to excessive stress on the tissue solid phase, making it more susceptible to damage [58]. In this study, we observed that the mode of flexion had a large impact on the overall fluid contribution. Particularly, our model predicted a ∼50% fluid pressurization contribution before flexion, and torque-driven flexion did not alter the relative fluid pressure contribution (∼50%; Fig. 7(b)—diagonal bars). However, muscle-driven flexion decreased the relative fluid contribution with ICR distance, suggesting a reduced protective role from interstitial fluid. For example, the relative fluid contribution was 15% for case J, representing a ∼70% decrease compared to the preflexion configuration (Fig. 7(b)—solid bars). The decrease in fluid contribution under muscle-driven flexion was paired with increases in joint-level solid stresses (Fig. 7(a)) and AF strains (Fig. 9), which increased the overall possibility of tissue- and joint-level failure, making the disc more susceptible to herniation.

Despite extensive research on disc joint-level failure mechanics, a consensus has not been reached regarding a failure criterion, due to challenges in observing and accurately characterizing tissue damage in vivo or in situ. In addition to the effective strain- and fiber stretch-based failure criterion applied in this study, we evaluated maximum shear strain as a potential failure criterion to predict tissue failure. AF maximum shear strains and the failed element percentage never exceeded the failure threshold (Fig. 10(b)). Thus, failure was not predicted for any loading condition, making maximum shear strain an ineffective failure criterion. By contrast, an agreement between model-predicted failure location (i.e., endplate for the torque-driven models and posterolateral AF for the muscle-driven models) with in vitro and clinical observations demonstrated the predictive power of our current failure criterion based on AF local effective strain and AF fiber stretch. Additionally, in our previous work, the local effective strain was shown to be an effective and accurate predictor for bulk AF tissue failure, with a 90% agreement between model predictions and experimental observations [52]. Regardless, in vitro experiments replicating the torque- and muscle-driven flexion models are required to validate the model predictions obtained in this study to fully evaluate the accuracy and robustness of the failure criterion applied.

Axial rotation combined with flexion and axial compression has been shown to increase the risk of herniation in vitro; thus, it was recommended that a combination of at least these three loading modalities was applied for repeatable herniation [5,10,60,61]. We chose to not include axial rotation in this study, as the ICR location for axial rotation is a variable that also requires parametric evaluation. Interestingly, though axial rotation was not included in the loading protocol, our model predicted herniation failure through the posterolateral AF for at least three out of seven muscle-driven cases (i.e., cases H to J). This could potentially make in vitro assessment of herniation more accessible to a wider range of researchers due to less demanding testing equipment (i.e., the tester does not need to support simultaneous rotation around the transverse and sagittal axis). Furthermore, joint-level failure mechanical tests could potentially benefit from a simpler testing protocol, as differences in mechanical test setups and protocols can introduce hard-to-identify variations in measured mechanics, making it difficult to compare data across groups [59,62,63].

Model predictions highlighted disparate multiscale disc mechanics, not only with different flexion setups, but with respect to different ICR locations. For example, cases D and I were both considered muscle-driven but may represent different activities or postures (Fig. 11). Particularly, model simulations showed large differences in bulk deformation, stress-bearing mechanisms, and intradiscal stress and strain distributions between cases D and I (Figs. 58), resulting in differences in failure risk and failure behavior predicted (Figs. 9 and 10). Thus, it is reasonable to assume that disc mechanics would vary considerably with ICR location under other physiologically relevant degrees-of-freedom, including axial rotation and lateral bending. Similar to common flexion and extension tests, most mechanical testing protocols defined the ICR on the disc for axial rotation and lateral bending [5,64]. Thus, by employing a similar study design framework, the current model can be further applied to investigate variations in disc mechanics under rotation and lateral bending.

Fig. 11
Disc deformation and stress–strain distribution under two physiological flexion postures. The ICR is located (a) close to the body, and (b) away from the body.
Fig. 11
Disc deformation and stress–strain distribution under two physiological flexion postures. The ICR is located (a) close to the body, and (b) away from the body.
Close modal

One limitation to this study was that the risk of in vitro herniation was only evaluated based on AF strain, AF fiber stretch, and maximum shear strain; however, previous studies have suggested that AF failure might be driven by stress [65], strain energy density [66], or a combination of stress and strain (i.e., Tsai–Wu damage criterion) [67]. The inclusion of a wider range of failure criteria could further improve the robustness of model predictions, generating more accurate and precise conclusions regarding failure initiation and progression. Secondly, a welded contact was assumed between the fibers and matrix as well as between interlamellar interfaces, and other contact mechanisms were not assessed. Although the contact mechanism is not well understood, it is likely that they could change with failure initiation and progression, thus altering tissue stress and strain distributions [6870]. Additionally, generalized disc material and geometric properties based on average measurements reported in the literature were assumed and employed in the models developed in this study. However, disc mechanical and nutrient transport behaviors have been shown to be sensitive to variations in disc geometries and morphologies [71,72]. Thus, future studies that intend to investigate multiscale mechanics of discs of specimen-specific geometries or morphologies should consider conducting corresponding sensitivity analyses a priori to help evaluate the effect of geometry and material property on model predictions.

Finite element modeling provides a powerful and effective tool for assessing multiscale and multiphasic disc mechanics under loading conditions that are difficult to setup experimentally. The multiscale and multiphasic structure-based model used in this study demonstrated significant differences in mechanical behavior and risk of failure from torque- and muscle-driven flexion. Specifically, model results highlighted the effectiveness of muscle-driven flexion in provoking herniation in vitro. In conclusion, this study provided a potential computational framework for designing improved in vitro mechanical testing protocols for the intervertebral disc, which can advance the assessment of disc failure both in vitro and in vivo.

Funding Data

  • National Science Foundation (NSF) (Grant Nos. 1760467 and 1751212; Funder ID: 10.13039/100000001).

References

1.
Yao
,
M.
,
Xu
,
B. P.
,
Li
,
Z. J.
,
Zhu
,
S.
,
Tian
,
Z. R.
,
Li
,
D. H.
,
Cen
,
J.
,
Cheng
,
S. D.
,
Wang
,
Y. J.
,
Guo
,
Y. M.
, and
Cui
,
X. J.
,
2020
, “
A Comparison Between the Low Back Pain Scales for Patients With Lumbar Disc Herniation: Validity, Reliability, and Responsiveness
,”
Health Qual. Life Outcomes
,
18
(
1
), pp.
1
12
.10.1186/s12955-020-01403-2
2.
Katz
,
J. N.
,
2006
, “
Lumbar Disc Disorders and Low-Back Pain: Socioeconomic Factors and Consequences
,”
JBJS
,
88
(
suppl_2
), pp.
21
24
.10.2106/JBJS.E.01273
3.
Schroeder
,
G. D.
,
Guyre
,
C. A.
, and
Vaccaro
,
A. R.
,
2016
, “
The Epidemiology and Pathophysiology of Lumbar Disc Herniations
,”
Semin. Spine Surg.
,
28
(
1
), pp.
2
7
.10.1053/j.semss.2015.08.003
4.
Truumees
,
E.
,
2015
, “
A History of Lumbar Disc Herniation From Hippocrates to the 1990s
,”
Clin. Orthop. Relat. Res.
,
473
(
6
), pp.
1885
1895
.10.1007/s11999-014-3633-7
5.
Wilke
,
H. J.
,
Kienle
,
A.
,
Maile
,
S.
,
Rasche
,
V.
, and
Berger-Roscher
,
N.
,
2016
, “
A New Dynamic Six Degrees of Freedom Disc-Loading Simulator Allows to Provoke Disc Damage and Herniation
,”
Eur. Spine J.
,
25
(
5
), pp.
1363
1372
.10.1007/s00586-016-4416-5
6.
White
,
A. A.
III
, and
Panjabi
,
M. M.
,
1990
, “
Clinical Biomechanics of the Spine
,” Lippincott, Philadelphia, PA.
7.
Adams
,
M. A.
, and
Hutton
,
W. C.
,
1983
, “
The Effect of Fatigue on the Lumbar Intervertebral Disc
,”
J. Bone Joint Surg. Brit. Vol.
,
65-B
(
2
), pp.
199
203
.10.1302/0301-620X.65B2.6826631
8.
Adams
,
M. A.
, and
Hutton
,
W. C.
,
1983
, “
The Mechanics of Prolapsed Intervertebral Disc
,”
Int. Orthop.
,
6
(
4
), pp.
249
253
.10.1007/BF00267146
9.
Adams
,
M. A.
, and
Hutton
,
W. C.
,
1985
, “
Gradual Disc Prolapse
,”
Spine
,
10
(
6
), pp.
524
531
.10.1097/00007632-198507000-00006
10.
Berger-Roscher
,
N.
,
Casaroli
,
G.
,
Rasche
,
V.
,
Villa
,
T.
,
Galbusera
,
F.
, and
Wilke
,
H. J.
,
2017
, “
Influence of Complex Loading Conditions on Intervertebral Disc Failure
,”
Spine
,
42
(
2
), pp.
E78
E85
.10.1097/BRS.0000000000001699
11.
Shirazi-Adl
,
S. A.
,
Shrivastava
,
S. C.
, and
Ahmed
,
A. M.
,
1984
, “
Stress Analysis of the Lumbar Disc-Body Unit in Compression. A Three-Dimensional Nonlinear Finite Element Study
,”
Spine
,
9
(
2
), pp.
120
134
.10.1097/00007632-198403000-00003
12.
Kurowski
,
P.
, and
Kubo
,
A. I. Z. O. H.
,
1986
, “
The Relationship of Degeneration of the Intervertebral Disc to Mechanical Loading Conditions on Lumbar Vertebrae
,”
Spine
,
11
(
7
), pp.
726
731
.10.1097/00007632-198609000-00012
13.
Kim
,
Y. E.
,
Goel
,
V. K.
,
Weinstein
,
J. N.
, and
Lim
,
T. H.
,
1991
, “
Effect of Disc Degeneration at One Level on the Adjacent Level in Axial Mode
,”
Spine
,
16
(
3
), pp.
331
335
.10.1097/00007632-199103000-00013
14.
Shirazi-Adl
,
A.
,
1992
, “
Finite-Element Simulation of Changes in the Fluid Content of Human Lumbar Discs. Mechanical and Clinical Implications
,”
Spine
,
17
(
2
), pp.
206
212
.10.1097/00007632-199202000-00015
15.
Rohlmann
,
A.
,
Zander
,
T.
,
Schmidt
,
H.
,
Wilke
,
H. J.
, and
Bergmann
,
G.
,
2006
, “
Analysis of the Influence of Disc Degeneration on the Mechanical Behaviour of a Lumbar Motion Segment Using the Finite Element Method
,”
J. Biomech.
,
39
(
13
), pp.
2484
2490
.10.1016/j.jbiomech.2005.07.026
16.
Schmidt
,
H.
,
Kettler
,
A.
,
Rohlmann
,
A.
,
Claes
,
L.
, and
Wilke
,
H. J.
,
2007
, “
The Risk of Disc Prolapses With Complex Loading in Different Degrees of Disc Degeneration–A Finite Element Analysis
,”
Clin. Biomech.
,
22
(
9
), pp.
988
998
.10.1016/j.clinbiomech.2007.07.008
17.
Schmidt
,
H.
,
Kettler
,
A.
,
Heuer
,
F.
,
Simon
,
U.
,
Claes
,
L.
, and
Wilke
,
H. J.
,
2007
, “
Intradiscal Pressure, Shear Strain, and Fiber Strain in the Intervertebral Disc Under Combined Loading
,”
Spine
,
32
(
7
), pp.
748
755
.10.1097/01.brs.0000259059.90430.c2
18.
Galbusera
,
F.
,
Schmidt
,
H.
,
Neidlinger-Wilke
,
C.
,
Gottschalk
,
A.
, and
Wilke
,
H. J.
,
2011
, “
The Mechanical Response of the Lumbar Spine to Different Combinations of Disc Degenerative Changes Investigated Using Randomized Poroelastic Finite Element Models
,”
Eur. Spine J.
,
20
(
4
), pp.
563
571
.10.1007/s00586-010-1586-4
19.
Galbusera
,
F.
,
Schmidt
,
H.
,
Neidlinger-Wilke
,
C.
, and
Wilke
,
H. J.
,
2011
, “
The Effect of Degenerative Morphological Changes of the Intervertebral Disc on the Lumbar Spine Biomechanics: A Poroelastic Finite Element Investigation
,”
Comput. Methods Biomech. Biomed. Eng.
,
14
(
8
), pp.
729
739
.10.1080/10255842.2010.493522
20.
Barthelemy
,
V. M. P.
,
Van Rijsbergen
,
M. M.
,
Wilson
,
W.
,
Huyghe
,
J. M.
,
Van Rietbergen
,
B.
, and
Ito
,
K.
,
2016
, “
A Computational Spinal Motion Segment Model Incorporating a Matrix Composition-Based Model of the Intervertebral Disc
,”
J. Mech. Behav. Biomed. Mater.
,
54
, pp.
194
204
.10.1016/j.jmbbm.2015.09.028
21.
Castro
,
A. P. G.
, and
Alves
,
J. L.
,
2021
, “
Numerical Implementation of an Osmo-Poro-Visco-Hyperelastic Finite Element Solver: Application to the Intervertebral Disc
,”
Comput. Methods Biomech. Biomed. Eng.
,
24
(
5
), pp.
538
550
.10.1080/10255842.2020.1839059
22.
Ehlers
,
W.
,
Karajan
,
N.
, and
Markert
,
B.
,
2009
, “
An Extended Biphasic Model for Charged Hydrated Tissues With Application to the Intervertebral Disc
,”
Biomech. Model. Mechanobiol.
,
8
(
3
), pp.
233
251
.10.1007/s10237-008-0129-y
23.
Zhou
,
M.
,
Lim
,
S.
, and
O'Connell
,
G. D.
,
2021
, “
A Robust Multiscale and Multiphasic Structure-Based Modeling Framework for the Intervertebral Disc
,”
Front. Bioeng. Biotechnol.
,
9
, p.
452
.10.3389/fbioe.2021.685799
24.
Zhou
,
M.
,
Werbner
,
B.
, and
O'Connell
,
G. D.
,
2021
, “
Fiber Engagement Accounts for Geometry-Dependent Annulus Fibrosus Mechanics: A Multiscale, Structure-Based Finite Element Study
,”
J. Mech. Behav. Biomed. Mater.
,
115
, p.
104292
.10.1016/j.jmbbm.2020.104292
25.
Maas
,
S. A.
,
Ellis
,
B. J.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2012
, “
FEBio: Finite Elements for Biomechanics
,”
ASME J. Biomech. Eng.
,
134
(
1
), p. 011005.10.1115/1.4005694
26.
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
(
3
), pp.
328
333
.10.1097/01.brs.0000253961.40910.c1
27.
Zhou
,
M.
,
Bezci
,
S. E.
, and
O'Connell
,
G. D.
,
2020
, “
Multiscale Composite Model of Fiber-Reinforced Tissues With Direct Representation of Sub-Tissue Properties
,”
Biomech. Model. Mechanobiol.
,
19
(
2
), pp.
745
759
.10.1007/s10237-019-01246-x
28.
Marchand
,
F.
, and
Ahmed
,
A. M.
,
1990
, “
Investigation of the Laminate Structure of Lumbar Disc Anulus Fibrosus
,”
Spine
,
15
(
5
), pp.
402
410
.10.1097/00007632-199005000-00011
29.
Goel
,
V. K.
,
Monroe
,
B. T.
,
Gilbertson
,
L. G.
, 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
(
6
), pp.
689
698
.10.1097/00007632-199503150-00010
30.
Michalek
,
A. J.
,
Buckley
,
M. R.
,
Bonassar
,
L. J.
,
Cohen
,
I.
, and
Iatridis
,
J. C.
,
2009
, “
Measurement of Local Strains in Intervertebral Disc Anulus Fibrosus Tissue Under Dynamic Shear: Contributions of Matrix Fiber Orientation and Elastin Content
,”
J. Biomech.
,
42
(
14
), pp.
2279
2285
.10.1016/j.jbiomech.2009.06.047
31.
Schollum
,
M. L.
,
Robertson
,
P. A.
, and
Broom
,
N. D.
,
2010
, “
How Age Influences Unravelling Morphology of Annular Lamellae–a Study of Interfibre Cohesivity in the Lumbar Disc
,”
J. Anat.
,
216
(
3
), pp.
310
319
.10.1111/j.1469-7580.2009.01197.x
32.
Yu
,
J.
,
Tirlapur
,
U.
,
Fairbank
,
J.
,
Handford
,
P.
,
Roberts
,
S.
,
Winlove
,
C. P.
,
Cui
,
Z.
, and
Urban
,
J.
,
2007
, “
Microfibrils, Elastin Fibres and Collagen Fibres in the Human Intervertebral Disc and Bovine Tail Disc
,”
J. Anat.
,
210
(
4
), pp.
460
471
.10.1111/j.1469-7580.2007.00707.x
33.
Adam
,
C.
,
Rouch
,
P.
, and
Skalli
,
W.
,
2015
, “
Inter-Lamellar Shear Resistance Confers Compressive Stiffness in the Intervertebral Disc: An Image-Based Modelling Study on the Bovine Caudal Disc
,”
J. Biomech.
,
48
(
16
), pp.
4303
4308
.10.1016/j.jbiomech.2015.10.041
34.
Matcher
,
S. J.
,
Winlove
,
C. P.
, and
Gangnus
,
S. V.
,
2004
, “
The Collagen Structure of Bovine Intervertebral Disc Studied Using Polarization-Sensitive Optical Coherence Tomography
,”
Phys. Med. Biol.
,
49
(
7
), pp.
1295
1306
.10.1088/0031-9155/49/7/016
35.
Berg-Johansen
,
B.
,
Han
,
M.
,
Fields
,
A. J.
,
Liebenberg
,
E. C.
,
Lim
,
B. J.
,
Larson
,
P. E.
,
Gunduz-Demir
,
C.
,
Kazakia
,
G. J.
,
Krug
,
R.
, and
Lotz
,
J. C.
,
2018
, “
Cartilage Endplate Thickness Variation Measured by Ultrashort Echo-Time MRI is Associated With Adjacent Disc Degeneration
,”
Spine
,
43
(
10
), pp.
E592
E600
.10.1097/BRS.0000000000002432
36.
Lai
,
W. M.
,
Hou
,
J. S.
, and
Mow
,
V. C.
,
1991
, “
A Triphasic Theory for the Swelling and Deformation Behaviors of Articular Cartilage
,”
ASME J. Biomech. Eng.
, 113(3), pp.
245
258
.10.1115/1.2894880
37.
Beckstein
,
J. C.
,
Sen
,
S.
,
Schaer
,
T. P.
,
Vresilovic
,
E. J.
, and
Elliott
,
D. M.
,
2008
, “
Comparison of Animal Discs Used in Disc Research to Human Lumbar Disc: Axial Compression Mechanics and Glycosaminoglycan Content
,”
Spine
,
33
(
6
), pp.
E166
E173
.10.1097/BRS.0b013e318166e001
38.
Bezci
,
S. E.
,
Werbner
,
B.
,
Zhou
,
M.
,
Malollari
,
K. G.
,
Dorlhiac
,
G.
,
Carraro
,
C.
,
Streets
,
A.
, and
O'Connell
,
G. D.
,
2019
, “
Radial Variation in Biochemical Composition of the Bovine Caudal Intervertebral Disc
,”
JOR Spine
,
2
(
3
), p.
e1065
.10.1002/jsp2.1065
39.
Gu
,
W. Y.
,
Yao
,
H.
,
Vega
,
A. L.
, and
Flagler
,
D.
,
2004
, “
Diffusivity of Ions in Agarose Gels and Intervertebral Disc: Effect of Porosity
,”
Ann. Biomed. Eng.
,
32
(
12
), pp.
1710
1717
.10.1007/s10439-004-7823-4
40.
Partanen
,
J. I.
,
Partanen
,
L. J.
, and
Vahteristo
,
K. P.
,
2017
, “
Traceable Thermodynamic Quantities for Dilute Aqueous Sodium Chloride Solutions at Temperatures From (0 to 80) C. Part 1. Activity Coefficient, Osmotic Coefficient, and the Quantities Associated With the Partial Molar Enthalpy
,”
J. Chem. Eng. Data
,
62
(
9
), pp.
2617
2632
.10.1021/acs.jced.7b00091
41.
Cortes
,
D. H.
,
Jacobs
,
N. T.
,
DeLucca
,
J. F.
, and
Elliott
,
D. M.
,
2014
, “
Elastic, Permeability and Swelling Properties of Human Intervertebral Disc Tissues: A Benchmark for Tissue Engineering
,”
J. Biomech.
,
47
(
9
), pp.
2088
2094
.10.1016/j.jbiomech.2013.12.021
42.
Van Der Rijt
,
J. A.
,
Van Der Werf
,
K. O.
,
Bennink
,
M. L.
,
Dijkstra
,
P. J.
, and
Feijen
,
J.
,
2006
, “
Micromechanical Testing of Individual Collagen Fibrils
,”
Macromol. Biosci.
,
6
(
9
), pp.
697
702
.10.1002/mabi.200600063
43.
Shen
,
Z. L.
,
Dodge
,
M. R.
,
Kahn
,
H.
,
Ballarini
,
R.
, and
Eppell
,
S. J.
,
2008
, “
Stress-Strain Experiments on Individual Collagen Fibrils
,”
Biophys. J.
,
95
(
8
), pp.
3956
3963
.10.1529/biophysj.107.124602
44.
Choi
,
K.
,
Kuhn
,
J. L.
,
Ciarelli
,
M. J.
, and
Goldstein
,
S. A.
,
1990
, “
The Elastic Moduli of Human Subchondral, Trabecular, and Cortical Bone Tissue and the Size-Dependency of Cortical Bone Modulus
,”
J. Biomech.
,
23
(
11
), pp.
1103
1113
.10.1016/0021-9290(90)90003-L
45.
Goel
,
V. K.
,
Ramirez
,
S. A.
,
Kong
,
W.
, and
Gilbertson
,
L. G.
,
1995
, “
Cancellous Bone Young's Modulus Variation Within the Vertebral Body of a Ligamentous Lumbar Spine—Application of Bone Adaptive Remodeling Concepts
,”
ASME J. Biomech. Eng.
, 117(3), pp.
266
271
.10.1115/1.2794180
46.
Périé
,
D.
,
Korda
,
D.
, and
Iatridis
,
J. C.
,
2005
, “
Confined Compression Experiments on Bovine Nucleus Pulposus and Annulus Fibrosus: Sensitivity of the Experiment in the Determination of Compressive Modulus and Hydraulic Permeability
,”
J. Biomech.
,
38
(
11
), pp.
2164
2171
.10.1016/j.jbiomech.2004.10.002
47.
Farrell
,
M. D.
, and
Riches
,
P. E.
,
2013
, “
On the Poisson's Ratio of the Nucleus Pulposus
,”
ASME J. Biomech. Eng.
,
135
(
10
), p.
104501
.10.1115/1.4025180
48.
Wu
,
Y.
,
Cisewski
,
S. E.
,
Sachs
,
B. L.
,
Pellegrini
,
V. D.
, Jr.
,
Kern
,
M. J.
,
Slate
,
E. H.
, and
Yao
,
H.
,
2015
, “
The Region-Dependent Biomechanical and Biochemical Properties of Bovine Cartilaginous Endplate
,”
J. Biomech.
,
48
(
12
), pp.
3185
3191
.10.1016/j.jbiomech.2015.07.005
49.
Demers
,
C. N.
,
Antoniou
,
J.
, and
Mwale
,
F.
,
2004
, “
Value and Limitations of Using the Bovine Tail as a Model for the Human Lumbar Spine
,”
Spine
,
29
(
24
), pp.
2793
2799
.10.1097/01.brs.0000147744.74215.b0
50.
Alini
,
M.
,
Eisenstein
,
S. M.
,
Ito
,
K.
,
Little
,
C.
,
Kettler
,
A. A.
,
Masuda
,
K.
,
Melrose
,
J.
,
Ralphs
,
J.
,
Stokes
,
I.
, and
Wilke
,
H. J.
,
2008
, “
Are Animal Models Useful for Studying Human Disc Disorders/Degeneration?
,”
Eur. Spine J.
,
17
(
1
), pp.
2
19
.10.1007/s00586-007-0414-y
51.
O'Connell
,
G. D.
,
Johannessen
,
W.
,
Vresilovic
,
E. J.
, and
Elliott
,
D. M.
,
2007
, “
Human Internal Disc Strains in Axial Compression Measured Noninvasively Using Magnetic Resonance Imaging
,”
Spine
,
32
(
25
), pp.
2860
2868
.10.1097/BRS.0b013e31815b75fb
52.
Werbner
,
B.
,
Zhou
,
M.
, and
O'Connell
,
G.
,
2017
, “
A Novel Method for Repeatable Failure Testing of Annulus Fibrosus
,”
ASME J. Biomech. Eng.
,
139
(
11
), p. 111001.10.1115/1.4037855
53.
Heuer
,
F.
,
Schmidt
,
H.
, and
Wilke
,
H. J.
,
2008
, “
The Relation Between Intervertebral Disc Bulging and Annular Fiber Associated Strains for Simple and Complex Loading
,”
J. Biomech.
,
41
(
5
), pp.
1086
1094
.10.1016/j.jbiomech.2007.11.019
54.
Amin
,
D. B.
,
Moawad
,
C. M.
, and
Costi
,
J. J.
,
2019
, “
New Findings Confirm Regional Internal Disc Strain Changes During Simulation of Repetitive Lifting Motions
,”
Ann. Biomed. Eng.
,
47
(
6
), pp.
1378
1390
.10.1007/s10439-019-02250-z
55.
Danso
,
E. K.
,
Honkanen
,
J. T. J.
,
Saarakkala
,
S.
, and
Korhonen
,
R. K.
,
2014
, “
Comparison of Nonlinear Mechanical Properties of Bovine Articular Cartilage and Meniscus
,”
J. Biomech.
,
47
(
1
), pp.
200
206
.10.1016/j.jbiomech.2013.09.015
56.
Yan
,
Y.
,
Fan
,
H.
,
Li
,
Y.
,
Hoeglinger
,
E.
,
Wiesinger
,
A.
,
Barr
,
A.
,
O'Connell
,
G. D.
, and
Harris-Adamson
,
C.
,
2021
, “
Applying Wearable Technology and a Deep Learning Model to Predict Occupational Physical Activities
,”
Appl. Sci.
,
11
(
20
), p.
9636
.10.3390/app11209636
57.
Proctor
,
C. S.
,
Schmidt
,
M. B.
,
Whipple
,
R. R.
,
Kelly
,
M. A.
, and
Mow
,
V. C.
,
1989
, “
Material Properties of the Normal Medial Bovine Meniscus
,”
J. Orthop. Res.
,
7
(
6
), pp.
771
782
.10.1002/jor.1100070602
58.
Ateshian
,
G. A.
,
Lai
,
W. M.
,
Zhu
,
W. B.
, and
Mow
,
V. C.
,
1994
, “
An Asymptotic Solution for the Contact of Two Biphasic Cartilage Layers
,”
J. Biomech.
,
27
(
11
), pp.
1347
1360
.10.1016/0021-9290(94)90044-2
59.
Werbner
,
B.
,
Zhou
,
M.
,
McMindes
,
N.
,
Lee
,
A.
,
Lee
,
M.
, and
O'Connell
,
G. D.
,
2022
, “
Saline-Polyethylene Glycol Blends Preserve In Vitro Annulus Fibrosus Hydration and Mechanics: An Experimental and Finite-Element Analysis
,”
J. Mech. Behav. Biomed. Mater.
,
125
, p.
104951
.10.1016/j.jmbbm.2021.104951
60.
Veres
,
S. P.
,
Robertson
,
P. A.
, and
Broom
,
N. D.
,
2009
, “
The Morphology of Acute Disc Herniation: A Clinically Relevant Model Defining the Role of Flexion
,”
Spine
,
34
(
21
), pp.
2288
2296
.10.1097/BRS.0b013e3181a49d7e
61.
Veres
,
S. P.
,
Robertson
,
P. A.
, and
Broom
,
N. D.
,
2010
, “
The Influence of Torsion on Disc Herniation When Combined With Flexion
,”
Eur. Spine J.
,
19
(
9
), pp.
1468
1478
.10.1007/s00586-010-1383-0
62.
Newell
,
N.
,
Rivera Tapia
,
D.
,
Rahman
,
T.
,
Lim
,
S.
,
O'Connell
,
G. D.
, and
Holsgrove
,
T. P.
,
2020
, “
Influence of Testing Environment and Loading Rate on Intervertebral Disc Compressive Mechanics: An Assessment of Repeatability at Three Different Laboratories
,”
JOR Spine
,
3
(
3
), p.
e21110
.10.1002/jsp2.1110
63.
Costi
,
J. J.
,
Ledet
,
E. H.
, and
O'Connell
,
G. D.
,
2021
, “
Spine Biomechanical Testing Methodologies: The Controversy of Consensus Vs Scientific Evidence
,”
JOR Spine
,
4
(
1
), p.
e1138
.10.1002/jsp2.1138
64.
Bezci
,
S. E.
,
Klineberg
,
E. O.
, and
O'Connell
,
G. D.
,
2018
, “
Effects of Axial Compression and Rotation Angle on Torsional Mechanical Properties of Bovine Caudal Discs
,”
J. Mech. Behav. Biomed. Mater.
,
77
, pp.
353
359
.10.1016/j.jmbbm.2017.09.022
65.
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
(
3
), pp.
125
140
.10.1007/s10237-004-0053-8
66.
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
,”
ASME J. Biomech. Eng.
, 132(6), p.
061007
.10.1115/1.4001032
67.
Shahraki
,
N. M.
,
Fatemi
,
A.
,
Agarwal
,
A.
, and
Goel
,
V. K.
,
2017
, “
Prediction of Clinically Relevant Initiation and Progression of Tears Within Annulus Fibrosus
,”
J. Orthop. Res.
,
35
(
1
), pp.
113
122
.10.1002/jor.23346
68.
Bruehlmann
,
S. B.
,
Matyas
,
J. R.
, and
Duncan
,
N. A.
,
2004
, “
ISSLS Prize Winner: Collagen Fibril Sliding Governs Cell Mechanics in the Anulus Fibrosus: An in Situ Confocal Microscopy Study of Bovine Discs
,”
Spine
,
29
(
23
), pp.
2612
2620
.10.1097/01.brs.0000146465.05972.56
69.
Vergari
,
C.
,
Mansfield
,
J.
,
Meakin
,
J. R.
, and
Winlove
,
P. C.
,
2016
, “
Lamellar and Fibre Bundle Mechanics of the Annulus Fibrosus in Bovine Intervertebral Disc
,”
Acta Biomater.
,
37
, pp.
14
20
.10.1016/j.actbio.2016.04.002
70.
Szczesny
,
S. E.
,
Fetchko
,
K. L.
,
Dodge
,
G. R.
, and
Elliott
,
D. M.
,
2017
, “
Evidence That Interfibrillar Load Transfer in Tendon is Supported by Small Diameter Fibrils and Not Extrafibrillar Tissue Components
,”
J. Orthop. Res.
,
35
(
10
), pp.
2127
2134
.10.1002/jor.23517
71.
Sélard
,
É.
,
Shirazi-Adl
,
A.
, and
Urban
,
J. P.
,
2003
, “
Finite Element Study of Nutrient Diffusion in the Human Intervertebral Disc
,”
Spine
,
28
(
17
), pp.
1945
1953
.10.1097/01.BRS.0000087210.93541.23
72.
Schlager
,
B.
,
Niemeyer
,
F.
,
Galbusera
,
F.
,
Volkheimer
,
D.
,
Jonas
,
R.
, and
Wilke
,
H. J.
,
2018
, “
Uncertainty Analysis of Material Properties and Morphology Parameters in Numerical Models Regarding the Motion of Lumbar Vertebral Segments
,”
Comput. Methods Biomech. Biomed. Eng.
,
21
(
12
), pp.
673
83
.10.1080/10255842.2018.1508571

Supplementary data