A Robust Multiscale and Multiphasic Structure-Based Modeling Framework for the Intervertebral Disc

A comprehensive understanding of multiscale and multiphasic intervertebral disc mechanics is crucial for designing advanced tissue engineered structures aiming to recapitulate native tissue behavior. The bovine caudal disc is a commonly used human disc analog due to its availability, large disc height and area, and similarities in biochemical and mechanical properties to the human disc. Because of challenges in directly measuring subtissue-level mechanics, such as in situ fiber mechanics, finite element models have been widely employed in spinal biomechanics research. However, many previous models use homogenization theory and describe each model element as a homogenized combination of fibers and the extrafibrillar matrix while ignoring the role of water content or osmotic behavior. Thus, these models are limited in their ability in investigating subtissue-level mechanics and stress-bearing mechanisms through fluid pressure. The objective of this study was to develop and validate a structure-based bovine caudal disc model, and to evaluate multiscale and multiphasic intervertebral disc mechanics under different loading conditions and with degeneration. The structure-based model was developed based on native disc structure, where fibers and matrix in the annulus fibrosus were described as distinct materials occupying separate volumes. Model parameters were directly obtained from experimental studies without calibration. Under the multiscale validation framework, the model was validated across the joint-, tissue-, and subtissue-levels. Our model accurately predicted multiscale disc responses for 15 of 16 cases, emphasizing the accuracy of the model, as well as the effectiveness and robustness of the multiscale structure-based modeling-validation framework. The model also demonstrated the rim as a weak link for disc failure, highlighting the importance of keeping the cartilage endplate intact when evaluating disc failure mechanisms in vitro. Importantly, results from this study elucidated important fluid-based load-bearing mechanisms and fiber-matrix interactions that are important for understanding disease progression and regeneration in intervertebral discs. In conclusion, the methods presented in this study can be used in conjunction with experimental work to simultaneously investigate disc joint-, tissue-, and subtissue-level mechanics with degeneration, disease, and injury.

Introduction

Mechanical dysfunction of the intervertebral disc can lead to reduced mobility and debilitating pain (Adams and Roughley, 2006). Disc prolapse and herniation mostly occur in the posterolateral region, where stresses, strains, and intradiscal pressure in the annulus fibrosus (AF) are higher (Shah et al., 1978; Adams and Hutton, 1985; Steffen et al., 1998; O’Connell et al., 2007b; O’Connell et al., 2011; Wilke et al., 2016; Liu et al., 2017). The posterolateral region has also been linked to increasing bulging and protrusion of the nucleus pulposus under fatigue, with some discs experiencing full herniations (Wilke et al., 2016). Previous researchers have tracked progression of disc failure from bulging to herniation (Adams et al., 2000; Vernon-Roberts et al., 2007), but further investigation is limited due to experimental challenges in directly assessing in situ mechanics (e.g., fiber mechanics), which result in large variations in reported in situ fiber mechanics data. For example, earlier in vitro joint-level studies reported AF fiber strains that varied from ∼0.3 to 20% under axial compression, which may cause contradicting predictions regarding the likelihood of disc failure under physiological conditions (Shah et al., 1978; Stokes, 1987; Heuer et al., 2008a,b, 2012; Wang et al., 2009; Spera et al., 2011). Thus, despite recent advancements in experimental techniques, in situ fiber mechanics at the joint level remain poorly understood.

Human intervertebral disc cadaveric tissues are the benchmark for spine biomechanics research, but limited tissue availability and challenges in controlling for important variables, such as sex, age, and level of degeneration, can impact study designs (e.g., sample size) and confound results (Iatridis et al., 2005; Alini et al., 2008; Michalek and Iatridis, 2012; Costi et al., 2020). For these reasons, many researchers have resorted to large animal models, including ovine, porcine, and bovine, to investigate intervertebral disc biomechanics (Alini et al., 2008). Particularly, bovine caudal discs are more accessible than human discs, easier to handle than discs from smaller animals (e.g., rat and mouse discs), and have biochemical and mechanical properties similar to human discs (Demers et al., 2004; Beckstein et al., 2008; Showalter et al., 2012; Bezci et al., 2019). Furthermore, previous work demonstrated the effectiveness of using bovine discs to study the effect of injuries and degeneration by effectively inducing injuries (e.g., needle punctures) and degeneration (e.g., enzyme digestion) in the tissues in vitro (Korecki et al., 2008a; Roberts et al., 2008; Michalek and Iatridis, 2012). Despite improvements in availability, accessibility, consistency, and ease of manipulation, experimental limitations still prevent assessment of intradiscal deformations and stress distributions between disc components with injuries or degeneration. Instead, in vitro studies primarily assess joint-level bulk mechanics, compositional changes, or biological response (Oshima et al., 1993; Korecki et al., 2008a,b; Roberts et al., 2008; Walter et al., 2011; Michalek and Iatridis, 2012; Bezci et al., 2015, 2020a,b; Bezci and O’Connell, 2018). The growing wealth of data that can be obtained from the bovine caudal discs makes it an ideal animal model to develop a validated and comprehensive computational tool to assess in situ mechanics. Additionally, because of lower inter-specimen variability, bovine disc models can be more effectively and reliably validated with experimental data than human disc models.

Finite element models (FEM) have been used to complement experimental studies, providing a powerful tool for predicting hard-to-measure, three-dimensional mechanical and biochemical responses (Zhou et al., 2020c). Since the 1970s, FEMs have advanced the field of spinal biomechanics significantly by providing insights into disc joint-level mechanics and tissue-level stress and strain distributions (Shirazi-Adl et al., 1984; Shirazi-Adl, 1992; Galbusera et al., 2011a,b; Schmidt et al., 2013). However, many joint-level FEMs describe disc components as single-phasic elastic or hyperelastic materials and thus do not account for water content (Kurowski and Kubo, 1986; Kim et al., 1991; Rohlmann et al., 2006; Schmidt et al., 2007b), which is a primary constituent in all biological tissues and plays an important role in the tissue’s load-bearing capability (Ateshian et al., 1994). More recent models have accounted for tissue water content by describing disc components as poroelastic materials, which significantly advanced the field by enabling investigations into the stress-bearing role of the interstitial tissue water content, as well as tissue’s time-dependent behavior (Natarajan et al., 2006; Wilson et al., 2007; Galbusera et al., 2011a,b; Barthelemy et al., 2016; van Rijsbergen et al., 2018; Castro and Alves, 2020). However, these models have limited capability in describing the osmotic response, which has been shown to alter mechanical behavior and change with degeneration (Ishihara et al., 1996; Wognum et al., 2006; Wuertz et al., 2007).

In addition to the limitations in accounting for tissue’s fluid content and osmotic response, most FEMs are developed based on homogenization theory, where every model element includes a homogenized description of tissue subcomponents (e.g., fibers and extrafibrillar matrix) and, thus, does not accurately represent the heterogeneous AF native architecture, where fibers and extrafibrillar matrix are distinct materials that occupy separate volumes. As a result, these models are not capable of directly investigating subtissue-level mechanics (e.g., in situ fiber or interfibrillar stress and strain distributions; Yin and Elliott, 2005). To address some of these issues, we previously developed and validated a structure-based FEM of the AF that replicated its native tissue architecture, with fiber bundles modeled as a separate material from the extrafibrillar matrix (Zhou et al., 2020a). In this approach, model parameters directly represented tissue mechanical (e.g., modulus, Poisson’s ratio, etc.) or biochemical properties (e.g., proteoglycan content, referential hydraulic permeability, etc.). To account for tissue water content and osmotic behavior, triphasic mixture theory was employed to describe the swelling capacity of the extrafibrillar matrix (Lai et al., 1991; Ateshian et al., 2004). Our model was able to robustly and accurately predict multilamellar AF mechanics under various loading configurations and testing boundary conditions, including uniaxial tension, biaxial tension, and simple shear (Zhou et al., 2020a). More recently, by incorporating a structure-based fiber engagement analysis, we were also able to apply this model to explain the relationship between specimen geometry and AF tensile mechanics that was originally observed by Adams and Green (1993) and Zhou et al. (2020b).

The objective of this study was to expand our structure-based multiscale modeling-validation approach to study joint-level mechanics of the intervertebral disc under both healthy and degenerated conditions. Degeneration has been shown to alter subtissue-level fiber mechanics, which plays an important role in stress distributions, damage accumulation, and bulk tissue failure (Werbner et al., 2019). Understanding mechanisms of stress distribution within the disc and its subcomponents can help develop robust designs for tissue repair or replacement implants, such as tissue engineered discs. Therefore, we (1) developed and validated a joint-level FEM that was capable of investigating the multiscale and multiphasic structure-function relationship in bovine caudal discs, and (2) used the validated FEM to investigate the effect of loading condition and degeneration on multiscale disc mechanics at joint, tissue, and subtissue scales.

Materials and Methods

Model Development

FEMs were developed to represent a bone-disc-bone motion segment from the bovine tail (Figure 1A). Neighboring tissues (e.g., facet joints, ligaments, etc.) were not included in the model to minimize confounding effects and to more closely represent motion segment specimens prepared for experimental testing. Model geometry was created in Solidworks (2020) and finite element meshes were generated using ABAQUS and ANSA pre-processor (Abaqus 6.14; ANSA 15.2.0). Mesh size was determined based on results from our previous mesh convergence study (Zhou et al., 2020b). PreView was used to define boundary and loading conditions and the fully developed models were solved by FEBio (PreView 2.1; FEBio 2.8.5; Maas et al., 2012). Due to limited computational resources, the current available solver was only able to process a maximum of ∼200 million non-zero entries in the stiffness matrix. Thus, models created in this study were scaled down at 1:5 scale.

FIGURE 1

Figure 1. (A) Schematic of the multiscale, structure-based bovine caudal disc motion segment model. The extrafibrillar matrix and collagen fibers of the annulus fibrosus (AF) were modeled as distinct materials occupying separate volumes. Insets present the cartilage endplate geometry (top) and the angle-ply fiber structure (bottom right). (B) AF fiber angle and solid volume fraction from the inner AF (IAF) to the outer AF (OAF). (C) Fixed charge density distribution in healthy and degenerated (Degen) disc models.

To ensure that this scaling and the resulting changes in the number of AF lamellae modeled did not affect model predictions, preliminary work was performed to determine the effect of scaling ratio between 1:4 and 1:6 on model-predicted compressive and torsional mechanics. Compressive stress-strain behavior and normalized torsional stiffness-rotation response from the 1:4, 1:5, and 1:6 scale models were consistent (Supplementary Figure 1), suggesting that scaling and number of AF lamellae modeled did not affect model predictions when the model included enough AF lamellae. Thus, bovine caudal disc motion segment models were developed at 1:5 scale for computational efficiency (∼2.1 million elements). Finite element meshes of the model were shown in Supplementary Figure 2.

Model geometry was determined based on data reported in the literature. At full scale, the radius and height of bovine caudal discs are 14.20 ± 0.85 mm and 6.90 ± 0.35 mm, respectively, assuming a circular cross section in the transverse plane (O’Connell et al., 2007a). Thus, the 1:5 scaled model radius and height (not including both bony endplates) were created at 2.85 and 1.40 mm, respectively (Figure 1A). The nucleus pulposus (NP) was assumed to have the same circular cross section in the transverse plane, but with a ∼50% smaller radius (1.45 mm; Figure 1A; O’Connell et al., 2007a). The AF was created using our previously reported structure-based modeling approach, where the tissue was described as a fiber-reinforced angle-ply composite containing distinct materials for fiber bundles and the extrafibrillar matrix (Figure 1A; Zhou et al., 2020a). Due to limited computational resources, the native bovine AF structural features, including lamellar thickness, fiber radius, and interfibrillar spacing, were preserved during scaling to reduce the total number of elements needed. This scaling approach, which has been widely applied and validated for human disc models (Shirazi-Adl et al., 1984; Goel et al., 1995a; Galbusera et al., 2011a,b), maintained fiber volume fraction and preserved mesh quality for model convergence and model predictions (Zhou et al., 2020b). As such, seven concentric AF layers were created (lamellar thickness = 0.2 mm; Adam et al., 2015). Fiber bundles were uniformly distributed, full-length cylinders welded to the surrounding matrix (Goel et al., 1995a; Michalek et al., 2009; Schollum et al., 2010). Due to the lack of bovine caudal disc anatomy data in the literature, fiber bundle geometry from the human AF was used, based on the similar collagen networks reported between human and bovine discs (Yu et al., 2002, 2007). Specifically, the fiber bundle radius was 0.06 mm, and interfibrillar spacing within each lamella was 0.22 mm (Marchand and Ahmed, 1990). Fiber angles were oriented at ± 45° to the transverse plane in the inner AF and decreased along the radial direction to ± 30° in the outer AF (Figure 1A–bottom inset; Figure 1B–turquoise circles; Matcher et al., 2004). Cartilage endplates (CEP) covered the superior and inferior ends of the NP and the inner-middle AF (Figure 1A–cartilage endplate); spatial variation in CEP thickness was included based on data in the literature (Figure 1A–top inset; Berg-Johansen et al., 2018). Bony endplates were modeled to cover the superior and inferior ends of the disc (Figure 1A–bony endplate). All interfaces were defined as welded interfaces (Adam et al., 2015).

Triphasic mixture theory was employed to account for tissue water content and osmotic response (Lai et al., 1991; Ateshian et al., 2004). The Holmes-Mow description was employed 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), k₀ represented hydraulic permeability in the reference configuration, φ₀ represented tissue solid volume fraction, and M represented the exponential strain-dependence coefficient. Tissue fluid phase model parameters were determined based on reported values for bovine tissues when available (Table 1–Fluid phase). AF solid volume fraction (i.e., 100% minus water content as a percentage) varied linearly along the radial direction, increasing from 0.2 in the inner AF to 0.3 in the outer AF (Table 1 and Figure 1B–grayscale circles). Fixed charge density represented proteoglycan content in the NP, CEP, and AF extrafibrillar matrix, allowing for osmotic swelling. Radial variation in fixed charge density was determined based on our recent work that provided high-spatial-resolution measurements of bovine caudal disc biochemical composition (Figure 1C–solid bars; Bezci et al., 2019). The collagen fiber bundles were assumed to have no swelling capability (i.e., zero fixed charge density). Free diffusivity (D₀) and within-tissue diffusivity (D) of Na⁺ and Cl^– were set based on data reported in Gu et al. (2004); 100% ion solubility was assumed (D_0,Na⁺ = 0.00116mm²/s; D_0,Cl⁻ = 0.00161mm²/s; D_Na⁺ = 0.00044mm²/s; D_Cl⁻ = 0.00069 mm²/s). The solution osmotic coefficient (0.927) was determined based on a linear interpolation of data reported in Robinson and Stokes (1949) and Partanen et al. (2017).

TABLE 1

Table 1. Triphasic material properties of the bovine caudal disc tissues.

$k\left(J\right)=k_0{\left(\frac{J-{\mathrm{\phi }}_0}{1-{\mathrm{\phi }}_0}\right)}^2{e}^{\frac{1}{2}M\left(J^2-1\right)}(1)$

To describe NP, CEP, and AF extrafibrillar matrix mechanics, a compressible hyperelastic Holmes-Mow material description was used (Eqs 2–4; Cortes et al., 2014). Particularly, I₁ and I₂ represented the first and second invariants of the right Cauchy-Green deformation tensor, C(C = F^TF), E represented Young’s modulus, v represented Poisson’s ratio, and β represented the exponential stiffening coefficient. AF collagen fibers were modeled using the same compressible hyperelastic Holmes-Mow ground matrix but reinforced with a power-linear fiber description to account for AF non-linearity and anisotropy (Eq. 5). γ represented the power-law exponent in the toe region, E_lin represented the fiber modulus in the fiber linear region, and λ₀ represented the transition stretch between the toe and linear regions (Holzapfel and Ogden, 2017). B was a function of γ, E_lin., and λ₀ ($B=\frac{E_{lin}}{2}(\frac{({\mathrm{\lambda }}_0^2-1)}{2(\mathrm{\gamma }-1)}+{\mathrm{\lambda }}_0^2)$. Solid phase parameters were determined based on bovine experimental studies when available (Table 1–solid phase), and collagen fiber properties were determined based on type I collagen uniaxial tensile test experimental data (Table 1–solid phase: E_lin, γ, and λ₀). For all material properties, data from healthy human discs was used when bovine properties were not available, due to similarities in tissue properties (Table 1–“^∗”).

$W(I_1,I_2,J)=\frac{1}{2}c\left(e^Q-1\right)(2)$

$$\begin{array}{l}Q=\frac{\text{β}(1+\text{ν})(1-2\text{ν})}{E(1-\text{ν})}[(\frac{E}{1+\text{ν}}-\frac{E\text{ν}}{(1+\text{ν})(1-2\text{ν})})(I_1-3)\\ +\frac{E\text{ν}}{(1+\text{ν})(1-2\text{ν})}(I_2-3)-(\frac{E}{1+\text{ν}}+\frac{E\text{ν}}{(1+\text{ν})(1-2\text{ν})})ln{J}^2](3)\end{array}$$

$c=\frac{E(1-\mathrm{\nu })}{2\mathrm{\beta }\left(1+\mathrm{\nu }\right)\left(1-2\mathrm{\nu }\right)}(4)$

${\mathrm{\lambda }}_n=\{\begin{array}{c}0\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} }{\lambda }_n<1\hfill \\ \frac{E_{lin}}{4\mathrm{\gamma }(\mathrm{\gamma }-1)}{({\mathrm{\lambda }}_0^2-1)}^{2-\mathrm{\gamma }}{({\mathrm{\lambda }}_n-1)}^{\mathrm{\gamma }}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}\le {\mathrm{\lambda }}_n\le {\mathrm{\lambda }}_0\hfill \\ E_{lin}\left({\mathrm{\lambda }}_n-{\mathrm{\lambda }}_0\right)+B({\mathrm{\lambda }}_n^2-{\mathrm{\lambda }}_0^2)\hfill \\ +\frac{E_{lin}}{4\mathrm{\gamma }(\mathrm{\gamma }-1)}{({\mathrm{\lambda }}_0^2-1)}^{2-\mathrm{\gamma }}{({\mathrm{\lambda }}_n-1)}^{\mathrm{\gamma }}\mathit{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}{\mathrm{\lambda }}_n>{\mathrm{\lambda }}_0\hfill \end{array}(5)$

Bony endplates were modeled as a compressible hyperelastic material using the Neo-Hookean description (Eq. 6). I₁, I₂, J were defined as above.E_{bonyendplates} and ν_{bony  endplates} represented the Young’s modulus (12,000 MPa) and Poisson’s ratio (0.3) of the bony endplates, which were determined based on reported data in the literature (Choi et al., 1990; Goel et al., 1995b; Dreischarf et al., 2014).

$$\begin{array}{l}{W}_{bonyendplates}(I_1,I_2,J)\\ =\frac{E_{bonyendplates}}{4(1+{\text{ν}}_{bonyendplates})}(I_1-3)-\frac{E_{bonyendplates}}{2(1+{\text{ν}}_{bonyendplates})}lnJ\\ +\frac{E_{bonyendplates}{\text{ν}}_{bonyendplates}}{(1+{\text{ν}}_{bonyendplates})(1-2{\text{ν}}_{bonyendplates})}(lnJ{)}^2(6)\end{array}$$

Multiscale Model Validation

Model robustness and accuracy (i.e., predictive power) were evaluated by simulating a range of loading modalities tested in experiments. All models were simulated using steady-state analyses and the model output were evaluated at equilibrium. Model-predicted properties were compared to experimental measurements at the joint, tissue, and subtissue levels.

Joint-Level Validation

At the joint level, resting intradiscal pressure, compressive mechanics, and torsional mechanics were evaluated for the motion segment model described in Section “Model Development.” Resting intradiscal pressure was defined as the average NP pressure after swelling and was compared to in vivo and in vitro intradiscal pressure data (Urban and McMullin, 1988; Ishihara et al., 1996; Sato et al., 1999; Wilke et al., 1999; Nguyen et al., 2008). Both human intervertebral disc and bovine caudal disc intradiscal pressure data were included for validation, because previous studies have shown similar results between the two species (Oshima et al., 1993; Ishihara et al., 1996; Alini et al., 2008).

Disc compressive and torsional mechanics were evaluated by applying loading protocols described in corresponding experimental studies (Beckstein et al., 2008; Showalter et al., 2012). After swelling (triphasic) in 0.15 M phosphate-buffered saline, compressive mechanics were evaluated by applying a 0.5 MPa axial compression. Boundary conditions at the top and bottom bony endplates were defined to represent boundary conditions reported in Beckstein et al. (2008). The normalized compressive stiffness was calculated as the slope of the model-predicted compressive load-displacement curve in the linear region, which was then normalized by the model geometry (i.e., cross-sectional area and height; Beckstein et al., 2008). Torsional mechanics were evaluated by applying a 0.5 MPa axial compressive preload immediately followed by a 10° axial rotation. Boundary conditions at the top and bottom bony endplates were defined to represent boundary conditions reported in Showalter et al. (2012). Normalized torsional stiffness was calculated by normalizing the slope of the torque-rotation curve between 7.5° and 10° by the model polar moment of inertia (Showalter et al., 2012; Bezci et al., 2018). The model was considered valid for predicting disc intradiscal pressure and stiffness when model-predicted values were within one standard deviation of reported mean values.

To assess the influence of including water content and osmotic response on predicted mechanical behavior, a 1:5 hyperelastic disc model, which is more commonly used in FEMs of the intervertebral disc, was created. In the model, all disc components were modeled using hyperelastic material descriptions, and its compressive stiffness was evaluated by applying a 0.5 MPa axial compression and calculating the slope of the linear region of the stress-strain curve.

Tissue-Level Validation

At the tissue level, both model-predicted AF mechanical properties and swelling properties were evaluated for model validation. A structure-based FEM was created for bovine multilamellar AF tissue specimens to simulate uniaxial tensile tests performed by Vergari and coworkers (Figure 2A; Vergari et al., 2017). After swelling (triphasic) in 0.15 M phosphate-buffered saline, a 1.1 uniaxial tensile stretch was applied along the circumferential direction (Figure 2A). Boundary conditions were defined to represent no slipping between the grips and the multilamellar tissue sample surface, as reported in Vergari et al. (2017). Tensile modulus was calculated as the slope of the stress-stretch curve at stretch ratios between 1.02 and 1.06 in 0.01 increments, as reported in the literature (Vergari et al., 2017). Tissue explant models of the NP and inner-middle AF were created to evaluate model-predicted swelling behavior in 0.15 M phosphate-buffered saline. Swelling ratios were calculated as the difference between post- and pre-swelling weight divided by the tissue pre-swelling weight and compared to data reported in Bezci et al. (2019). If model-predicted mechanical and swelling properties were within one standard deviation of reported mean values, the model was considered valid for predicting the respective behavior.

FIGURE 2

Figure 2. (A) Model validation schematic for multilamellar mechanics of bovine annulus fibrosus (AF). Model geometry and loading conditions were determined based on protocols reported in Vergari et al. (2017). (B) Model-predicted (Mod-pred) bovine AF multilamellar stress-stretch response compared to representative experimental (EXP) data from Vergari et al. (2017). (C) Model-predicted tensile modulus at five specified stretch ratios compared to experimental data from Vergari et al. (2017). (D) Model validation for single lamellar mechanics of bovine AF. Model geometry and loading conditions were determined based on protocols reported in Monaco et al. (2016). (E) Model-predicted bovine single lamellar stress-stretch response. (F) Model-predicted bovine AF single lamellar tensile mechanical properties compared to experimental data [mean (standard deviation)] from Monaco et al. (2016).

Subtissue-Level Validation

At the subtissue level, model-predicted AF mechanics were evaluated for model validation. A structure-based model was created for bovine single lamellar AF specimens to simulate uniaxial tensile tests performed by Monaco and coworkers (Figure 2D; Monaco et al., 2016). After swelling (triphasic) in 0.15 M phosphate-buffered saline, a 1.5 uniaxial tensile stretch was applied to the specimen transverse to the fiber direction (Figure 2D). Boundary conditions were defined to effectively replicate the flexible rake system applied in Monaco et al. (2016). Model-predicted uniaxial tensile mechanics were only assessed transverse to the fiber direction, because to the best of the authors’ knowledge, no studies have evaluated bovine single lamellar AF mechanics along the fiber direction analogous to Holzapfel and coworkers’ work using the human AF (Holzapfel et al., 2005). Tensile modulus was calculated as the slope of the stress-stretch curve in the linear region. The model-predicted mechanical properties, including modulus and the stress and strain at the end of the toe-region, were compared to experimental data (Monaco et al., 2016). The model was considered valid for predicting subtissue-level mechanics if the model-predicted mechanical properties were within one standard deviation of reported mean values.

Effect of Loading Condition on Multiscale Bovine Caudal Disc Mechanics

After validation, three loading conditions were applied to the motion segment model described in Section “Model Development” to evaluate the effect of loading condition on multiscale bovine caudal disc mechanics. All three cases were loaded in two steps. First, swelling in 0.15 M phosphate-buffered saline was simulated. Then, one of the three loading conditions was assessed, including Case A: 0.5 MPa axial compression, Case B: 10° axial rotation, and Case C: 0.5 MPa axial precompression followed by 10° axial rotation. For Case A, axial compression was simulated between 0–1.0 MPa, but only data from 0.5 MPa axial compression was presented, as it corresponded to experimental data reported in the papers that we compared and validated our model to (Beckstein et al., 2008; Showalter et al., 2012; Bezci et al., 2018). Additionally, the 0.5 MPa axial compression more closely mimicked the compressive stress observed in low-intensity daily activities (e.g., relaxed standing and sitting, walking, etc.; Wilke et al., 1999). For Cases B and C, disc height was not allowed to change during rotation. Model boundary conditions were defined as in Section “Multiscale Model Validation,” while Cases B and C shared identical boundary conditions. All models were simulated using steady-state analyses with the output evaluated at equilibrium. The effect of loading condition was evaluated at the joint, tissue and subtissue levels, as follows:

Joint-Level Mechanics

Average solid stress (i.e., stress absorbed by tissue solid matrix) and fluid pressure (i.e., stress absorbed by the tissue interstitial fluid) of the entire bovine caudal disc, including the NP, AF, and CEP, were evaluated for all three cases. The relative contribution of solid stress was evaluated as the solid stress divided by the total stress, which was calculated as the sum of solid stress and fluid pressure based on triphasic mixture theory (Lai et al., 1991). Similarly, the relative contribution of fluid pressure was calculated by normalizing the fluid pressure by the total stress.

Tissue-Level Mechanics

NP, AF, and CEP in situ swelling ratios were evaluated post-swelling. After the applied mechanical loading, average solid stress, strain, and fluid pressure in the NP, AF, and CEP were evaluated for all three cases. For each disc component, the relative contribution of the solid stress and fluid pressure to the total stress was evaluated. The total stress was calculated as the sum of the component’s solid stress and fluid pressure. Disc bulging of the inner and outer AF was assessed under 0.5 MPa axial compression (Case A) and was calculated by dividing the respective change in mid-disc-height radius with loading by the post-swelling disc radius (reported as a percentage value).

Subtissue-Level Mechanics

Average fiber stretch was evaluated within each AF lamellae after swelling and after loading. Swelling-induced fiber stretch was calculated as the post-swelling fiber length divided by the initial fiber length. Post-loading fiber stretch was calculated as the post-mechanical loading fiber length divided by the post-swelling fiber length. Average solid stress in the fibers and extrafibrillar matrix was evaluated post-loading. The relative solid stress contribution of collagen fibers and extrafibrillar matrix to the overall AF solid stress, which was calculated as the sum of fiber and matrix solid stress, was also assessed. Additionally, post-loading fiber solid stress profiles along the fiber length from the inferior to the superior end of the disc were evaluated in both the inner- and outermost AF lamellae.

Effect of Degeneration on Multiscale Bovine Caudal Disc Mechanics

The effect of degeneration on multiscale disc mechanics was investigated under the three loading conditions evaluated in Section “Effect of Loading Condition on Multiscale Bovine Caudal Disc Mechanics.” Degeneration was achieved by reducing tissue proteoglycan content, which was simulated by reducing the fixed charge density in the NP, AF, and CEP (Adams and Roughley, 2006). Bovines are commonly slaughtered between 18 and 24 months and do not experience spontaneous degeneration within that timespan (Alini et al., 2008). Therefore, fixed charge density distribution for the degenerated disc was determined based on trends observed in degenerated human discs (Figure 1C–checkered bars; Urban and Maroudas, 1979; Beckstein et al., 2008; Bezci et al., 2019), as well as data reported from ex vivo degeneration models in relevant bioreactor studies (Castro et al., 2014; Paul et al., 2018). All model-predicted properties discussed and evaluated in Section “Effect of Loading Condition on Multiscale Bovine Caudal Disc Mechanics” were evaluated with degeneration. Additionally, model-predicted resting intradiscal pressure, normalized compressive stiffness, and normalized torsional stiffness were also calculated for the degenerated disc model and compared to available experimental data for a more rigorous model validation (Urban and McMullin, 1988; Sato et al., 1999; Showalter et al., 2012; Bezci et al., 2018). All models were simulated using steady-state analyses with the output evaluated at equilibrium.

Results

Multiscale Model Validation

Joint-Level Validation

Model-predicted intradiscal pressure value for the healthy disc was 0.17 MPa, which was within the range of reported experimental values (<0.90× standard deviation from reported mean values; Figure 3A–black diagonal bar vs. white bars enclosed by black lines).

FIGURE 3

Figure 3. (A) Model-predicted (Mod-pred) resting intradiscal pressure in healthy and degenerated (Degen) disc models compared to experimental (EXP) values. Data reported by Ishihara et al. (1996) (noted by *) were obtained from bovine caudal discs while data reported by the other listed studies were obtained from human intervertebral discs, which have shown to share comparable intradiscal pressure values. Variations were not reported in Wilke et al. (1999). (B) Representative model-predicted compressive (Comp) stress-strain response of hyperelastic (Hyper), healthy, and degenerated disc models under axial compression. (C) Model-predicted normalized (Norm) compressive stiffness (stiff) compared to EXP values. (D) Representative model-predicted torsional (tors) response of healthy and degenerated discs when evaluated for torsional mechanics. (E) Model-predicted normalized torsional stiffness compared to EXP values.

Model-predicted compressive stress-strain response was non-linear for healthy disc models developed using hyperelastic and triphasic mixture theory material descriptions, agreeing well with experimental observations (Figure 3B–solid lines). However, the hyperelastic disc model predicted a stiffer joint-level response than the triphasic model, which accounted for water content and osmotic behavior (Healthy). For the hyperelastic model, predicted normalized compressive stiffness was 12.52 MPa and did not agree with any available datasets (>1.2× standard deviations from reported means). Employing the triphasic material description resulted in a normalized compressive stiffness of 8.12 MPa, agreeing well with Beckstein et al. (2008) and two of three datasets collected, but not published, by Newell et al. (2020) (moduli calculated at a more relevant loading range than the previously published data, see Supplementary Figure 3). Model-predicted compressive stiffness was within 0.8 standard deviation of the reported mean for the three agreed datasets (Figure 3C–black diagonal bar vs. Beckstein et al., 2008 and Newell et al., 2020). However, our model was not able to accurately predict the compressive stiffness reported by the remaining dataset collected for Newell et al. (2020), which represents data from the authors’ own laboratory (18.74 ± 3.35 MPa, Supplementary Figure 3–Berkeley). The model-predicted compressive stiffness was >3.0 × standard deviations from the reported mean of this single dataset since the experimental data from our laboratory was higher than values reported by other institutes (Figure 3C–black diagonal bar vs. Newell et al., 2020).

A pseudo-linear torque-rotation response was observed for the healthy disc (Figure 3D–solid line). Model-predicted normalized torsional stiffness was 36 kPa/°, matching well with reported values (<0.75× standard deviation from the reported mean values; Figure 3E–black diagonal vs. white bars).

Tissue- and Subtissue-Level Validation

For multilamellar AF specimens, model-predicted stress-stretch response under uniaxial tension was non-linear, agreeing well with the literature (Figure 2B). Model-predicted tensile modulus agreed with the literature but tended to be on the higher end of reported values, particularly as stretch increased (Figure 2C). For single lamellar AF specimens, model-predicted stress-stretch response under uniaxial tension was also non-linear, agreeing well with the literature (Figure 2E). Model-predicted mechanical properties for the toe and linear regions were well within one standard deviation of the reported mean (<0.35× standard deviation from the reported mean; Figure 2F). Based on our model predictions, ex situ swelling ratio was 1.10 for the healthy NP tissue and 0.76 for the inner-middle AF, which were both within one standard deviation of the reported means (<0.88× standard deviation; Figure 4A).

FIGURE 4

Figure 4. (A) Model-predicted (Mod-pred) ex situ swelling ratios of the nucleus pulposus (NP) and the inner-middle annulus fibrosus (AF) compared to experimental (EXP) data reported by Bezci et al. (2019). (B,C) Model-predicted in situ swelling ratios of the NP, AF, and cartilage endplate (CEP) in healthy and degenerated (Degen) disc models. Relative changes in in situ swelling ratio with degeneration are labeled above corresponding neighboring bars.

Effect of Loading Condition on Multiscale Bovine Caudal Disc Mechanics

Joint-Level Mechanics

Fluid pressure contributed significantly to the disc’s overall load-bearing capacity, especially for loading conditions that incorporated axial compression. In healthy disc models, the average solid stress and average fluid pressure were both approximately 0.2 MPa under axial compression, resulting in relatively equal contribution to the total stress in the disc (Figure 5–Case A). Lower solid stress (0.11 MPa) and fluid pressure (0.13 MPa) were observed under axial rotation, but the relative contribution of solid stress and fluid pressure remained almost identical (Figure 5–Case B vs. A). Compared to Case A, the combined loading more than doubled the solid stress to 0.43 MPa but did not change the fluid pressure (0.24 MPa). Thus, the resulting relative contribution of the solid stress increased to 64% of the total stress (Figure 5–Case C vs. A).

FIGURE 5

Figure 5. Model-predicted (A) solid stress and fluid pressure, as well as (B) their relative contribution to the total stress taken by the disc in healthy and degenerated (Degen) models for Cases A, B, and C. Relative changes in solid stress or fluid pressure with degeneration are labeled above corresponding neighboring bars.

Tissue-Level Mechanics

Different applied boundary and loading conditions resulted in heterogeneous solid stress, fluid pressure, and strain distributions throughout the disc (Figure 6). Large solid stresses were observed in the outer AF, especially in Cases A and C (Figure 6A–“^∗”). Compared to Case A, the rotation-only loading condition resulted in lower solid stresses in all disc components (Figure 6A–Case B vs. A), where the solid stress in the NP, CEP, and AF decreased by more than 80, 67, and 42% (Figure 7A–Case B vs. A). Under combined loading, a two-fold increase in AF and CEP average solid stress was observed (Figure 7A–Case C vs. A: black and pink solid bars). However, the addition of rotation to axial compression did not change the NP solid stress (Figure 7A–Case C vs. A: green solid bar).

FIGURE 6

Figure 6. Representative post-loading disc mid-frontal (or coronal) plane (A) solid stress, (B) fluid pressure, and (C) strain distributions in healthy and degenerated (Degen) disc models. Black asterisks highlight stress concentrations. Black triangles point at strain concentrations.

FIGURE 7

Figure 7. Model-predicted post-loading average (A) solid stress, (B) fluid pressure, and (C) strain in the nucleus pulposus (NP), annulus fibrosus (AF), and cartilage endplate (CEP) in healthy and degenerated (Degen) disc models. Relative changes in NP, AF, and CEP solid stress, strain, or fluid pressure with degeneration are labeled above corresponding neighboring bars.

In situ swelling ratios for the NP, AF, and CEP were 0.25, 0.13, and 0.03, respectively (Figure 4B–Healthy; Figure 4C–black solid bars). Under axial compression, average fluid pressure was 0.14 MPa in the AF, which was ∼70% lower than that in the NP (0.47 MPa) and ∼60% lower than that in the CEP (0.36 MPa; Figure 7B–Case A: solid bars). Fluid pressure under the torsion-only loading was generally lower than that under the compression-only loading. Particularly, compared to Case A, NP and AF fluid pressure were both ∼40% lower while CEP fluid pressure was ∼60% lower (Figure 7B–Case B vs. A). Interestingly, compared to the compression-only loading condition, combining axial compression with rotation did not have a significant effect on the fluid pressure in any disc components (Figure 7B–Case C vs. A).

As expected, the relative fluid pressure to the total stress was significant and tissue-specific. Across all three loading conditions, fluid pressure accounted for more than 85% of the total stress in the NP and more than 70% of the total stress in the CEP (Figure 8–NP). The relative contribution of fluid pressure was smaller in the AF, but nevertheless accounted for 20–36% of the total AF stress (Figure 8–AF). Compared to the compression-only loading condition, the torsion-only loading resulted in a slight increase in the relative fluid pressure in the NP (Figure 8–Case B vs. A). However, the combined loading did not alter the relative solid stress or fluid pressure contribution in the NP but resulted in a ∼25% larger solid stress contribution in the AF (Figure 8–Case C vs. A). The relative solid and fluid contribution in the CEP was not affected by applied loading conditions (Figure 8–CEP).

FIGURE 8

Figure 8. Model-predicted relative contribution of solid stress and fluid pressure in the nucleus pulposus (NP), annulus fibrosus (AF), and cartilage endplate (CEP) in healthy and degenerated (Degen) disc models for Cases A, B, and C after the applied mechanical loading.

Large strains were observed at the AF-NP-CEP interface (i.e., the rim) and in the outer AF (Figure 6C–“^”). Under axial compression, NP and AF strains were comparable (0.16 and 0.13, respectively) and were approximately twofold greater than strains in the CEP (0.07; Figure 7C–Case A). Under axial rotation, strains in the NP decreased by ∼75%; however, AF and CEP strains increased by ∼20% (Figure 7C–Case B vs. A). Compression combined with rotation increased AF strains by 80% from 0.13 to 0.24 and increased CEP strains by more than 200% from 0.07 to 0.18. However, the combined loading did not greatly alter NP strains (∼5% change; Figure 7C–Case C vs. A).

Assessment of AF radial displacement at the mid-disc height under axial compression showed outward bulging for both the inner and outer AF after swelling (Figure 9A). In the outer AF, the relative outward bulging increased with applied load, reaching ∼1.8% under 0.5 MPa axial compression (Figure 9B–black solid circles). In the inner AF, the relative bulging reached a maximum of ∼0.4% under 0.2 MPa of compression but then decreased with additional applied compressive load (Figure 9B–red solid circles).

FIGURE 9

Figure 9. (A) Disc mid-frontal (or coronal) cross sections demonstrating the relative annulus fibrosus (AF) bulging in healthy and degenerated (Degen) disc models under axial compression. The relative AF bulging was calculated using the post-swelling 0 MPa configuration as the reference configuration (Ref config). (B) Relative bulging in the inner and outer AF in healthy and degenerated disc models. Positive and negative relative bulging suggest outward and inward AF bulging compared to the reference configuration, respectively. The gray horizontal dashed line represents the relative disc bulging threshold, below which the AF was predicted to bulge inward compared to the reference configuration.

Subtissue-Level Mechanics

The triphasic swelling step applied to all model cases prior to the applied mechanical loading resulted in an average swelling-induced fiber stretch of 1.05 in the inner AF and 1.02 in the outer AF. After applying 0.5 MPa of axial compression, the post-loading fiber stretch was ∼1.05 and was relatively consistent throughout the AF (Figure 10A–black solid circles). The magnitude of fiber stretch under the torsion-only loading was comparable, but there was a linear increase in fiber stretch from the innermost AF layer (1.04) to the outermost layer (1.07; Figure 10A–blue solid circles). Under the combined loading, the fiber stretch was nearly twofold greater than that under the single-axis loading conditions and was ∼1.10 through the AF (Figure 10A–red solid circles).

FIGURE 10

Figure 10. (A) Model-predicted average annulus fibrosus (AF) post-loading fiber stretch along the disc radial direction from the inner AF (IAF) to outer AF (OAF) in healthy and degenerated (Degen) disc models. The gray horizontal dashed line highlights the fiber stretch threshold, above which the fibers have a more significant chance of failure based on previous experimental observations. The threshold value was determined based on data reported by Skaggs et al. (1994) and Isaacs et al. (2014). The inset presents the relative (Rel) percentage change in average fiber stretch with degeneration along the disc radial direction. (B) Model-predicted post-loading average AF solid stress along the disc radial direction from the IAF to OAF. The inset presents the relative increase in fiber solid stress with degeneration for Case A.

Average fiber solid stress was relatively consistent throughout the AF under axial compression, ranging from 0.22 MPa in the inner AF to 0.29 MPa in the outer AF (Figure 10B–black solid circles). Under the rotation-only loading, fiber stress in the inner AF was 60% lower than the compression-only condition; however, large changes in fiber solid stress were not observed in the outer AF (Figure 10B–blue vs. black solid circles). Under the combined loading, fiber stress increased linearly from 0.37 MPa in the inner AF to 0.80 MPa in the outer AF. Compared to Case A, the fiber stress was increased by 70% in the inner AF and by 300% in the outer AF (Figure 10B–red vs. black solid circles). The solid stress of AF extrafibrillar matrix, as well as its observed trends with loading condition were both comparable to that of the fibers. Thus, across all three loading conditions, AF collagen fibers and extrafibrillar matrix contributed equally to the overall AF solid stress (Supplementary Figure 4).

Fiber solid stress profiles were tracked along the fiber length between the inferior and superior bony endplates. In all cases, fiber solid stress distributions were symmetric about the mid-transverse plane, due to disc symmetry (Figure 11). For Cases A and C, peak fiber solid stresses in the outer AF were observed right below the bony endplates, and peak fiber solid stresses in the inner AF were observed at the mid-disc height (Figure 11–Cases A and C: solid lines). By contrast, fiber stress was relatively consistent along the fiber length in both the inner and outer AF for Case B (Figure 11–Case B: solid lines). The combined loading amplified the fiber stress difference between the inner- and outermost lamellae, which shared comparable fiber stresses under the compression- or rotation-only loading conditions (Figure 11–solid black vs. gray lines).

FIGURE 11

Figure 11. Model-predicted post-loading annulus fibrosus (AF) fiber solid stress profiles along the fiber length from the inferior (Infer) to the superior (Super) bony endplates. The stress distributions were evaluated for the inner- and outermost AF layers in both healthy and degenerated (Degen) discs.

Effect of Degeneration on Multiscale Bovine Caudal Disc Mechanics

Joint-Level Mechanics

Resting intradiscal pressure decreased by ∼70% with degeneration (0.048 MPa) and was within the range of reported values (<0.10× standard deviation from the reported mean values; Figure 3A–red bars). Normalized compressive stiffness increased by ∼30% with degeneration (10.67 MPa; Figures 3B,C). Normalized torsional stiffness was approximately 37 kPa/°, which was not affected by degeneration (Figures 3D,E).

With degeneration, stresses were redistributed with the tissue solid component taking on more of the overall total stress (Figure 5–Degen vs. Healthy). Across the three loading conditions, degeneration increased solid stress by 18–66%, depending on the disc components, and the greatest relative increase with degeneration was observed in the compression-only loading condition (Figure 5A–checkered vs. solid bars). Fluid pressure decreased by ∼60% for all three loading conditions. Thus, the resulting relative contribution of solid stress increased from 45–65% in the healthy discs to 75–85% in the degenerated discs (Figure 5B–checkered vs. solid bars).

Tissue-Level Mechanics

As expected, degeneration reduced tissue swelling capability (Figures 4B,C–checkered vs. solid bars). The NP in situ swelling ratio reduced by >60%, decreasing from 0.25 to 0.09 with degeneration. Similarly, in situ AF swelling ratio decreased by ∼45% from 0.13 to 0.07 with degeneration. Interestingly, the CEP in situ swelling ratio became negative (−0.02) in the degenerated disc, indicating a loss of tissue volume after swelling (Figures 4B,C). The decrease in swelling capacity resulted in a 40–90% decrease in fluid pressure, depending on the tissue types and applied loading conditions. Particularly, large degeneration-induced fluid pressure decreases were mostly observed in the NP and CEP (Figure 7B–checkered vs. solids bars).

Similar to joint-level observations, degeneration redistributed stress in each disc component by decreasing the relative contribution of fluid pressure and increasing the relative contribution of solid stress (Figure 8–Degen vs. Healthy). The greatest stress redistribution was observed in the CEP, where the relative fluid pressure contribution decreased from ∼70–80% in the healthy discs to ∼20–50% in the degenerate discs. Noticeably, in Case B, the CEP relative fluid pressure contribution reduced by more than 75% from 83% in the healthy disc to 20% in the degenerate disc (Figure 8–CEP: checkered vs. solid bars). In the NP, the decrease in fluid contribution was relatively consistent for all three loading conditions. Particularly, degeneration reduced NP fluid contribution by ∼20–30%, decreasing from ∼85–95% in the healthy discs to ∼60–75% with degeneration (Figure 8–NP: checkered vs. solid bars). In the AF, the relative fluid pressure contribution decreased by ∼50% with degeneration, ranging from 11 to 17% in the degenerated discs compared to 20–36% in the healthy discs (Figure 8–AF: checkered vs. solid bars). Degeneration also increased the average strain in each disc components by ∼20–240%, with the largest increase observed in the CEP. Similar to the healthy disc, peak strains were observed at the AF-NP-CEP interface (i.e., the rim) and in the outer AF (Figure 6C–“^”).

The outer AF was still expected to bulge outward with the level of degeneration simulated in this study. Relative outward bulging for the outer AF at 0.5 MPa axial compression was ∼1%, which was ∼45% smaller than that in the degenerated disc (Figure 9–checkered vs. solid black circles). While the inner AF appeared to bulge outward slightly, calculating the relative change in radial displacement between the post-swelling and post-loading configuration showed that the inner AF moved inward toward the NP by 0.3% (Figure 9A–Degen; Figure 9B–checkered black circles). Although the inner AF moved toward the NP, collapse of the inner AF into the NP, which has been reported for more severely degenerated discs (Adams and Roughley, 2006), was not observed in our model.

Subtissue-Level Mechanics

Degeneration increased the average post-loading fiber stretch throughout the AF and had a greater impact on the inner AF than the outer AF (Figure 10A–checkered vs. solid black circles). For Case A, average fiber stretch decreased linearly from 1.10 in the inner AF to 1.07 in the outer AF (Figure 10A–checkered black circles), representing a 90% increase in fiber stretch in the inner AF and a 50% increase in the outer AF with degeneration (Figure 10A–inset: black circles). For Case B, the average fiber stretch was ∼1.08 and was relatively consistent throughout the AF (Figure 10A–checkered blue circles), where degeneration increased inner AF fiber stretch by more than 70% and increased outer AF fiber stretch by ∼20% (Figure 10A–inset: blue circles). Under the combined loading condition, average fiber stretch exceeded the 1.10 threshold in all AF lamellae, decreasing from 1.14 in the inner AF to 1.11 in the outer AF (Figure 10A–checkered red circles). However, although the inner AF fiber stretch increased by ∼50% with degeneration, the outer AF fiber stretch was not affected (Figure 10A–inset: red circles).

The overall increase in fiber stretch with degeneration did not result in a similar increase in fiber or extrafibrillar matrix solid stress. Under the compression-only loading, solid stress in the fibers increased by more than 40% in the inner AF and by ∼85% in the outer AF (Figure 10B–inset). However, the increases in both fiber and matrix solid stresses were smaller and not as consistent for Cases B and C (Figure 10B). Degeneration did not alter the AF fiber/matrix solid stress contribution (Supplementary Figure 4B), nor the pattern of stress distribution along the fiber length, but did increase the stress magnitude, with the largest increase observed for the compression-only loading (Figure 11–dashed vs. solid lines).

Discussion

This study developed and validated a multiscale and multiphasic structure-based finite element model of the bovine caudal disc motion segment. During development and validation, model parameters were determined based on tissue- or subtissue-level experimental data reported in the literature, as opposed to being calibrated to joint-level mechanics prior to validation. The model validation results highlight the model accuracy and robustness, as well as the advantages of employing the proposed multiscale, structure-based modeling-validation framework. After validation, the model was used to investigate the effect of loading condition and degeneration on solid stress, fluid pressure, and strain distributions at joint, tissue and subtissue scales. While only three loading conditions and one level of degeneration were assessed, results from this study demonstrate the model’s capability in investigating the shifts in disc load bearing or stress distribution mechanisms that can act to induce degenerative remodeling or damage accumulation.

Validation is critical for overall model performance, including accuracy and robustness. Most intervertebral disc models are only validated with respect to global disc measurements, such as axial displacement or intradiscal pressure. This limited validation approach can contribute to inaccurate model predictions, especially at tissue and subtissue scales, where model validation is not usually performed (Shirazi-Adl et al., 1984; Kim et al., 1991; Schmidt et al., 2007b; Galbusera et al., 2011a). Some studies calibrated model parameters, especially those associated with the AF, through optimization algorithms in order for the model predictions to fit experimental datasets measured in tests conducted under specific loading modalities (e.g., axial compression, flexion; Schmidt et al., 2006, 2007a; Malandrino et al., 2013); however, this framework requires models to be recalibrated for each new loading modality or disc geometry. The current study expanded upon our previously reported multiscale validation framework by performing model validation at joint, tissue, and subtissue levels (Zhou et al., 2020a; Figures 3, 4). A total of 16 validation cases were assessed and model-predicted properties agreed well with all but one dataset. Differences in joint stiffness between the outstanding dataset, which originate from our previous work, and our model predictions, are likely caused by the non-ideal machine compliance during experimental data collection (Newell et al., 2020). Importantly, model parameters were directly obtained from tissue- or subtissue-level experimental data and no adjustments were made to match tissue- or joint-level behavior. These results demonstrated the model’s predictive power and the effectiveness of the multiscale validation framework.

The structure-based modeling approach may improve clinical relevance and expand potential use for finite element models of the disc joint. At the tissue level, modeling discrete AF lamellae allowed for reproduction of radial variations in AF biochemical composition (i.e., proteoglycan content and water content). Describing variations in localized proteoglycan content is important for simulating and replicating morphological changes observed with degeneration, including the decrease in disc height, increased outward radial bulging, and inward bulging of the inner AF in severely degenerated discs (Yang and O’Connell, 2019). At the subtissue level, modeling collagen fiber bundles allowed us to explicitly evaluate fiber stress and strain distributions, rather than relying on indirect assessment, such as vector summation to evaluate fiber strain (Schmidt et al., 2007b). The separate fiber bundles generated more realistic predictions of in situ fiber mechanics and allowed for direct investigations into fiber-matrix interactions. For example, our findings demonstrate that a ∼50% decrease in proteoglycans caused a 40–90% increase in fiber stress when the disc was loaded under axial compression (Figure 10B–checkered vs. solid black circles). It should be noted that this study only assessed the moderate to severe degeneration level. Thus, additional work is needed to determine whether a decrease in only NP proteoglycan content, as observed in early degeneration, would result in similar increases in fiber stress.

Attributed to the structure-based modeling approach, the majority of our model parameters can be directly linked to tissue mechanical (e.g., modulus, Poisson’s ratio, etc.) or biochemical properties (e.g., water content, proteoglycan content, etc.; Table 1). Model parameters with physical significance help address concerns regarding overparameterization, which is a common issue associated with homogeneous finite element models, where model parameter calibration relies heavily on optimization algorithms (Yin and Elliott, 2005; Eskandari et al., 2019). Taken together, explicitly modeled disc structures with physically relevant model parameters benefit further investigations into disc joint behavior with degeneration, disease, or injury. For example, collagen fiber diameter and stiffness can be readily modified based on structural and mechanical changes noted with degeneration, or diseases such as diabetes (Adams and Roughley, 2006; Li et al., 2013; Svensson et al., 2018). Furthermore, the model can be easily modified to evaluate advanced tissue engineering designs (e.g., angle-ply disc replacements) before conducting costly and time-intensive in vivo studies in large animal models (Martin et al., 2014), or to help track time-dependent changes during bioreactor organ cultures (Frauchiger et al., 2018; Pfannkuche et al., 2020).

The importance of accounting for tissue water content and osmotic response was elucidated by assessing the relative stress contribution from tissue solid matrix and interstitial fluid (Figures 5–8). The contribution of fluid pressure plays a pivotal role in the disc’s load-bearing capacity (Adams and Roughley, 2006), but to the best of the authors’ knowledge, it has not been quantified. Inclusion of triphasic material properties allows for direct measurements of fluid pressure. Based on our model predictions for healthy discs, fluid pressure accounted for 35–55% of the total stress (Figure 5). More specifically, the fluid pressure contribution in the NP was greater than 85% (Figure 8), agreeing with previous findings for the healthy articular cartilage, which has a comparable fixed charge density and water content as healthy NP tissues (Maroudas et al., 1969; Armstrong and Mow, 1982; Lüssea et al., 2000; Shapiro et al., 2002). Degeneration reduced tissue swelling capacity, altering the disc’s load-bearing mechanism by shifting more stress to the tissue solid matrix (Figures 5, 8). This shift in stress-bearing was particularly noticeable under axial compression, where the decrease in fluid pressure (i.e., 0.13 MPa) was balanced by an equivalent increase in solid stress (Figure 5A–Case A). Despite the decrease in relative fluid pressure contribution with degeneration, fluid pressure still accounted for up to 25% of the total stress and contributed to more than 60% of NP stress (Figures 5, 8–checkered bars).

Models that do not incorporate tissue swelling describe stress as being entirely absorbed by the solid matrix (single-phasic hyperelastic material description), which likely contributed to overestimations in AF fiber stretch. For example, a previous model that employed single-phasic material descriptions for the disc predicted a fiber stretch of ∼1.12 under the rotation-only loading, even with the inclusion of posterior functional spinal structures (Schmidt et al., 2007b). However, experimental data on AF single lamellar tensile mechanics reported AF fiber bundle failure stretch as 1.14 ± 0.04 (Skaggs et al., 1994; Isaacs et al., 2014). Thus, such a model would suggest a relatively high likelihood of disc failure, contradicting to in vitro studies that showed low risk of disc failure under axial rotation (Berger-Roscher et al., 2017). The single-phasic material description may also help explain the overestimated compressive stiffness predicted by our hyperelastic model, as omission of water content and osmotic response led to higher AF solid matrix stress and larger fiber deformations that stiffened the disc joint (Figures 3B,C). Thus, our proposed model can potentially provide valuable insights into cell mechanobiology studies, as more accurate predictions of solid matrix stress and stretch data are required in order to apply physiological loading to cells or tissues in vitro (Martin et al., 2014).

The predictive power of our model was further demonstrated by evaluating the multiscale disc mechanics under different loading conditions and degeneration. Single-axis loading conditions (i.e., compression-only or rotation-only) resulted in a fiber solid stress <0.3 MPa and fiber stretch between 1.03 and 1.07 for the healthy disc model, which was comparable to in situ subfailure fiber stretch data obtained from photogrammetry-based studies (1.07–1.11; Heuer et al., 2008a,b; Heuer et al., 2012). Taken together, our model predictions for fiber stretch and stress suggest low risks of failure under the single-axis loading conditions, especially under axial rotation, as the average AF fiber stretch did not exceed 1.10 even with degeneration, which agrees well with recent six-degree of freedom testing results (Berger-Roscher et al., 2017). In contrast, multi-axis loading increased the likelihood of damage accumulation and disc failure as axial rotation combined with compression increased the average fiber stretch to 1.10 and almost tripled the average fiber stress in the outer AF from 0.3 to 0.9 MPa, which is much closer to the 1.0 MPa threshold reported in the literature (Skaggs et al., 1994; Holzapfel et al., 2005; Isaacs et al., 2014).

Degeneration increased the fiber stretch and fiber solid stress under all three simulated loading conditions, especially under the compression-only loading (Figure 10–Case A insets and Figure 11). Interestingly, under the combined loading, the average AF fiber stretch exceeded the 1.10 threshold for failure or significant damage accumulation (range: 1.11–1.14) but the average fiber solid stress still remained below 1.0 MPa. Taken together, these findings suggest that disc failures, especially those initiated in the AF (e.g., clefts, tears, etc.) may be strain-driven rather than stress-driven, agreeing with our previous tissue-level study (Werbner et al., 2017). Six degree of freedom testing machines provide the best approach for elucidating disc failure mechanisms in vitro (Costi et al., 2020). However, their high cost and complexity have limited their use. This model may provide a high-throughput approach to better understand the role of complex loading on damage accumulation and ultimate tissue failure (e.g., disc herniation).

Disc failure, especially those induced in vitro, have been commonly shown to occur through endplate fracture or annulus prolapse (Adams and Hutton, 1985; Wilke et al., 2016; Berger-Roscher et al., 2017). Across the three loading conditions evaluated, strain concentrations and peak fiber stresses were observed near the NP-AF-CEP interface and at the outer AF, especially in the degenerated disc (Figure 7C–“^”; Figure 11–gray solid lines). With degeneration, the CEP exhibited a volume loss post-swelling, likely caused by the compression from surrounding tissues due to differences in swelling capacities (Figure 4C). These results further highlight the NP-AF-CEP interface (i.e., the rim) as a weak link for disc failure. It should be noted that the flatter interface modeled between the CEP and the NP/AF was more representative of discs found in ovine, porcine, and human rather than bovine, which has a more concave CEP-NP-AF interface. Thus, it is within our expectations that our model-predicted peak stress and strain locations match well with in vitro failure locations observed in human and ovine discs (Adams and Hutton, 1985; Wilke et al., 2016; Berger-Roscher et al., 2017).

Although this study presents a strong validation and a robust modeling-validation framework, it is not without limitations. First, disc degeneration was simulated by only reducing tissue fixed charge density (i.e., proteoglycan content), without including any degeneration-related structural changes, such as AF lesions and decreased disc height. The omission of these structural or morphological changes might explain model predictions that contradicted previous experimental observations. For example, it has been widely accepted that degeneration results in higher disc flexibility in axial rotation, which was not predicted by our model within the simulated axial rotation range (Mimura et al., 1994; Galbusera et al., 2014). Additionally, previous experimental studies showed that annular bulging increases with degeneration and injury (Heuer et al., 2008b; Zou et al., 2009). While our model accurately predicted relative AF bulging in healthy discs (O’Connell et al., 2007b), it predicted that AF bulging decreased with degeneration (Figure 9–Degen vs. Healthy). Secondly, flexion/extension and lateral bending, which are important physiological loading modalities that have been shown to initiate disc failure at the CEP, were not assessed (Berger-Roscher et al., 2017). Ongoing and future work will include applying this multiscale, structure-based modeling-validation framework to human intervertebral discs to evaluate the risk of disc failure with early to moderate, or even more severe degenerative changes in tissue composition.

This study used a multiscale, structure-based modeling-validation framework to examine multiscale bovine caudal disc mechanics, including but not limited to fluid pressure, solid stress, and fiber stretch and strain. The model accurately predicted variations in disc mechanics under various loading conditions and with degeneration. Importantly, results from this study elucidated important load-bearing mechanisms and fiber-matrix interactions that are important for understanding disease progression and regeneration in intervertebral discs. In conclusion, the methods presented in this study can be used in conjunction with experimental work to simultaneously investigate disc joint-, tissue-, and subtissue-level mechanics with degeneration, disease, and injury.

Data Availability Statement

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

Author Contributions

MZ: conceptualization, methodology, software, validation, investigation, data collection and analysis, writing, review and editing, visualization, and project administration. SL: validation, investigation, writing, review and editing, and visualization. GDO: supervision, writing, review and editing, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

The work was supported by the National Science Foundation (CAREER BMMB: 1751212).

Conflict of Interest

The 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.

Acknowledgments

The authors thank Nicolas Newell from Imperial College London for his valuable advice. The authors also thank both Timothy P. Holsgrove from University of Exeter and Nicolas Newell for providing the permission to use the raw data collected for Newell et al. (2020).

Supplementary Material

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

References

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, 4303–4308. doi: 10.1016/j.jbiomech.2015.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, M. A., and Green, T. P. (1993). Tensile properties of the annulus fibrosus. Eur. Spine J. 2, 203–208. doi: 10.1007/bf00299447

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, M. A., and Hutton, W. C. (1985). Gradual disc prolapse. Spine 10, 524–531. doi: 10.1097/00007632-198507000-00006

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, M. A., and Roughley, P. J. (2006). What is intervertebral disc degeneration, and what causes it? Spine 31, 2151–2161. doi: 10.1097/01.brs.0000231761.73859.2c

CrossRef Full Text | Google Scholar

Adams, M. A., Freeman, B. J., Morrison, H. P., Nelson, I. W., and Dolan, P. (2000). Mechanical initiation of intervertebral disc degeneration. Spine 25, 1625–1636. doi: 10.1097/00007632-200007010-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Alini, M., Eisenstein, S. M., Ito, K., Little, C., Kettler, A. A., Masuda, K., et al. (2008). Are animal models useful for studying human disc disorders/degeneration? Eur. Spine J. 17, 2–19. doi: 10.1007/s00586-007-0414-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Armstrong, C. G., and Mow, V. C. (1982). Variations in the intrinsic mechanical properties of human articular cartilage with age, degeneration, and water content. J. Bone Joint Surg. Am. 64, 88–94. doi: 10.2106/00004623-198264010-00013

CrossRef Full Text | Google Scholar

Ateshian, G. A., Chahine, N. O., Basalo, I. M., and Hung, C. T. (2004). The correspondence between equilibrium biphasic and triphasic material properties in mixture models of articular cartilage. J. Biomech. 37, 391–400. doi: 10.1016/s0021-9290(03)00252-5

CrossRef Full Text | Google Scholar

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, 1347–1360. doi: 10.1016/0021-9290(94)90044-2

CrossRef Full Text | Google Scholar

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, 194–204. doi: 10.1016/j.jmbbm.2015.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

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, E166–E173.

Google Scholar

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, E78–E85.

Google Scholar

Berg-Johansen, B., Han, M., Fields, A. J., Liebenberg, E. C., Lim, B. J., Larson, P. E., et al. (2018). Cartilage endplate thickness variation measured by ultrashort echo-time MRI is associated with adjacent disc degeneration. Spine 43, E592–E600.

Google Scholar

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, 353–359. doi: 10.1016/j.jmbbm.2017.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezci, S. E., and O’Connell, G. D. (2018). Osmotic pressure alters time-dependent recovery behavior of the intervertebral disc. Spine 43, E334–E340.

Google Scholar

Bezci, S. E., Lim, S., and O’Connell, G. D. (2020a). Nonlinear stress-dependent recovery behavior of the intervertebral disc. J. Mech. Behav. Biomed. Mater. 110:103881. doi: 10.1016/j.jmbbm.2020.103881

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezci, S. E., Nandy, A., and O’Connell, G. D. (2015). Effect of hydration on healthy intervertebral disk mechanical stiffness. J. Biomech. Eng. 137:101007.

Google Scholar

Bezci, S. E., Torres, K., Carraro, C., Chiavacci, D., Werbner, B., Lim, S., et al. (2020b). Transient swelling behavior of the bovine caudal disc. J. Mech. Behav. Biomed. Mater. 112:104089. doi: 10.1016/j.jmbbm.2020.104089

PubMed Abstract | CrossRef Full Text | Google Scholar

Bezci, S. E., Werbner, B., Zhou, M., Malollari, K. G., Dorlhiac, G., Carraro, C., et al. (2019). Radial variation in biochemical composition of the bovine caudal intervertebral disc. JOR Spine 2:e1065.

Google Scholar

Castro, A. P. G., and Alves, J. L. (2020). Numerical implementation of an osmo-poro-visco-hyperelastic finite element solver: application to the intervertebral disc. Comput. Methods Biomech. Biomed. Eng. 28, 1–13. doi: 10.1080/10255842.2020.1839059

PubMed Abstract | CrossRef Full Text | Google Scholar

Castro, A. P. G., Paul, C. P. L., Detiger, S. E. L., Smit, T. H., Van Royen, B. J., Pimenta Claro, J. C., et al. (2014). Long-term creep behavior of the intervertebral disk: comparison between bioreactor data and numerical results. Front. Bioeng. Biotechnol. 2:56.

Google Scholar

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, 1103–1113. doi: 10.1016/0021-9290(90)90003-l

CrossRef Full Text | Google Scholar

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, 2088–2094. doi: 10.1016/j.jbiomech.2013.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Costi, J. J., Ledet, E. H., and O’Connell, G. D. (2020). Spine biomechanical testing methodologies: the controversy of consensus vs scientific evidence. JOR Spine 4:e1138.

Google Scholar

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, 2793–2799. doi: 10.1097/01.brs.0000147744.74215.b0

CrossRef Full Text | Google Scholar

Dreischarf, M., Zander, T., Shirazi-Adl, A., Puttlitz, C. M., Adam, C. J., Chen, C. S., 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. Biomech. 47, 1757–1766. doi: 10.1016/j.jbiomech.2014.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Eskandari, M., Nordgren, T. M., and O’Connell, G. D. (2019). Mechanics of pulmonary airways: linking structure to function through constitutive modeling, biochemistry, and histology. Acta Biomaterialia 97, 513–523. doi: 10.1016/j.actbio.2019.07.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Farrell, M. D., and Riches, P. E. (2013). On the poisson’s ratio of the nucleus pulposus. J. Biomech. Eng. 135, 104501–104504.

Google Scholar

Fratzl, P., Misof, K., Zizak, I., Rapp, G., Amenitsch, H., and Bernstorff, S. (1998). Fibrillar structure and mechanical properties of collagen. J. Struct. Biol. 122, 119–122. doi: 10.1006/jsbi.1998.3966

PubMed Abstract | CrossRef Full Text | Google Scholar

Frauchiger, D. A., Chan, S. C., Benneker, L. M., and Gantenbein, B. (2018). Intervertebral disc damage models in organ culture: a comparison of annulus fibrosus cross-incision versus punch model under complex loading. Eur. Spine J. 27, 1785–1797. doi: 10.1007/s00586-018-5638-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Galbusera, F., Schmidt, H., Neidlinger-Wilke, C., Gottschalk, A., and Wilke, H. J. (2011a). 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, 563–571. doi: 10.1007/s00586-010-1586-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Galbusera, F., Schmidt, H., Neidlinger-Wilke, C., and Wilke, H. J. (2011b). 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, 729–739. doi: 10.1080/10255842.2010.493522

PubMed Abstract | CrossRef Full Text | Google Scholar

Galbusera, F., Van Rijsbergen, M., Ito, K., Huyghe, J. M., Brayda-Bruno, M., and Wilke, H. J. (2014). Ageing and degenerative changes of the intervertebral disc and their impact on spinal flexibility. Eur. Spine J. 23, 324–332.

Google Scholar

Gentleman, E., Lay, A. N., Dickerson, D. A., Nauman, E. A., Livesay, G. A., and Dee, K. C. (2003). Mechanical characterization of collagen fibers and scaffolds for tissue engineering. Biomaterials 24, 3805–3813. doi: 10.1016/s0142-9612(03)00206-0

CrossRef Full Text | Google Scholar

Goel, V. K., Monroe, B. T., Gilbertson, L. G., and Brinckmann, P. (1995a). 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

CrossRef Full Text | Google Scholar

Goel, V. K., Ramirez, S. A., Kong, W., and Gilbertson, L. G. (1995b). Cancellous bone Young’s modulus variation within the vertebral body of a ligamentous lumbar spine—application of bone adaptive remodeling concepts. J. Biomed. Eng. 117, 266–271. doi: 10.1115/1.2794180

CrossRef Full Text | Google Scholar

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, 1710–1717. doi: 10.1007/s10439-004-7823-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Heuer, F., Schmidt, H., and Wilke, H. J. (2008a). Stepwise reduction of functional spinal structures increase disc bulge and surface strains. J. Biomech. 41, 1953–1960. doi: 10.1016/j.jbiomech.2008.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Heuer, F., Schmidt, H., and Wilke, H. J. (2008b). The relation between intervertebral disc bulging and annular fiber associated strains for simple and complex loading. J. Biomech. 41, 1086–1094. doi: 10.1016/j.jbiomech.2007.11.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Heuer, F., Schmidt, H., Käfer, W., Graf, N., and Wilke, H. J. (2012). Posterior motion preserving implants evaluated by means of intervertebral disc bulging and annular fiber strains. Clin. Biomech. 27, 218–225. doi: 10.1016/j.clinbiomech.2011.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A., and Ogden, R. W. (2017). Biomechanics: Trends in Modeling and Simulation, Vol. 316. Berlin: Springer.

Google Scholar

Holzapfel, G. A., Schulze-Bauer, C. A. J., 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

Iatridis, J. C., MacLean, J. J., and Ryan, D. A. (2005). Mechanical damage to the intervertebral disc annulus fibrosus subjected to tensile loading. J. Biomech. 38, 557–565. doi: 10.1016/j.jbiomech.2004.03.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Isaacs, J. L., Vresilovic, E., Sarkar, S., and Marcolongo, M. (2014). Role of biomolecules on annulus fibrosus micromechanics: effect of enzymatic digestion on elastic and failure properties. J. Mech. Behav. Biomed. Mater. 40, 75–84. doi: 10.1016/j.jmbbm.2014.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishihara, H., McNally, D. S., Urban, J. P., and Hall, A. C. (1996). Effects of hydrostatic pressure on matrix synthesis in different regions of the intervertebral disk. J. Appl. Physiol. 80, 839–846. doi: 10.1152/jappl.1996.80.3.839

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 331–335. doi: 10.1097/00007632-199103000-00013

PubMed Abstract | CrossRef Full Text | Google Scholar

Korecki, C. L., Costi, J. J., and Iatridis, J. C. (2008a). Needle puncture injury affects intervertebral disc mechanics and biology in an organ culture model. Spine 33, 235–241. doi: 10.1097/brs.0b013e3181624504

PubMed Abstract | CrossRef Full Text | Google Scholar

Korecki, C. L., MacLean, J. J., and Iatridis, J. C. (2008b). Dynamic compression effects on intervertebral disc mechanics and biology. Spine 33, 1403–1409. doi: 10.1097/brs.0b013e318175cae7

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 726–731. doi: 10.1097/00007632-198609000-00012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lai, W. M., Hou, J. S., and Mow, V. C. (1991). A triphasic theory for the swelling and deformation behaviors of articular cartilage. J. Biomech. Eng. 113, 245–258. doi: 10.1115/1.2894880

CrossRef Full Text | Google Scholar

Li, Y., Fessel, G., Georgiadis, M., and Snedeker, J. G. (2013). Advanced glycation end-products diminish tendon collagen fiber sliding. Matrix Biol. 32, 169–177. doi: 10.1016/j.matbio.2013.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., Yang, X. P., Li, K., Yang, T., Ye, J. D., and Zhang, C. Q. (2017). Internal strains of anulus fibrosus in the intervertebral disc under axial compression load. Biomed. Res. 8, 3483–3486.

Google Scholar

Lüssea, S., Claassen, H., Gehrke, T., Hassenpflug, J., Schünke, M., Heller, M., et al. (2000). Evaluation of water content by spatially resolved transverse relaxation times of human articular cartilage. Magn. Reson. Imaging 18, 423–430. doi: 10.1016/s0730-725x(99)00144-7

CrossRef Full Text | Google Scholar

Maas, S. A., Ellis, B. J., Ateshian, G. A., and Weiss, J. A. (2012). FEBio: finite elements for biomechanics. J. Biomech. Eng. 134:011005.

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 Biomech. Biomed. Eng. 16, 923–928. doi: 10.1080/10255842.2011.644539

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchand, F. R., and Ahmed, A. M. (1990). Investigation of the laminate structure of lumbar disc anulus fibrosus. Spine 15, 402–410. doi: 10.1097/00007632-199005000-00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Maroudas, A., Muir, H., and Wingham, J. (1969). The correlation of fixed negative charge with glycosaminoglycan content of human articular cartilage. Biochimica et Biophysica Acta 177, 492–500. doi: 10.1016/0304-4165(69)90311-0

CrossRef Full Text | Google Scholar

Martin, J. T., Milby, A. H., Chiaro, J. A., Kim, D. H., Hebela, N. M., Smith, L. J., et al. (2014). Translation of an engineered nanofibrous disc-like angle-ply structure for intervertebral disc replacement in a small animal model. Acta Biomater. 10, 2473–2481. doi: 10.1016/j.actbio.2014.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 1295–1306. doi: 10.1088/0031-9155/49/7/016

CrossRef Full Text | Google Scholar

Michalek, A. J., and Iatridis, J. C. (2012). Height and torsional stiffness are most sensitive to annular injury in large animal intervertebral discs. Spine J. 12, 425–432. doi: 10.1016/j.spinee.2012.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 2279–2285. doi: 10.1016/j.jbiomech.2009.06.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Mimura, M., Panjabi, M. M., Oxland, T. R., Crisco, J. J., Yamamoto, I., and Vasavada, A. (1994). Disc degeneration affects the multidirectional flexibility of the lumbar spine. Spine 19, 1371–1380. doi: 10.1097/00007632-199406000-00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Monaco, L. A., DeWitte-Orr, S. J., and Gregory, D. E. (2016). A comparison between porcine, ovine, and bovine intervertebral disc anatomy and single lamella annulus fibrosus tensile properties. J. Morphol. 277, 244–251. doi: 10.1002/jmor.20492

PubMed Abstract | CrossRef Full Text | Google Scholar

Natarajan, R. N., Williams, J. R., and Andersson, G. B. (2006). Modeling changes in intervertebral disc mechanics with degeneration. J. Bone Joint Surg. Am. 88(Suppl._2), 36–40. doi: 10.2106/jbjs.f.00002

PubMed Abstract | CrossRef Full Text | Google Scholar

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:e21110.

Google Scholar

Nguyen, A. M., Johannessen, W., Yoder, J. H., Wheaton, A. J., Vresilovic, E. J., Borthakur, A., et al. (2008). Noninvasive quantification of human nucleus pulposus pressure with use of T1ρ-weighted magnetic resonance imaging. J. Bone Joint Surg. Am. 90, 796–802. doi: 10.2106/jbjs.g.00667

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Connell, G. D., Johannessen, W., Vresilovic, E. J., and Elliott, D. M. (2007a). Human internal disc strains in axial compression measured noninvasively using magnetic resonance imaging. Spine 32, 2860–2868. doi: 10.1097/brs.0b013e31815b75fb

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

O’Connell, G. D., Vresilovic, E. J., and Elliott, D. M. (2011). Human intervertebral disc internal strain in compression: the effect of disc region, loading position, and degeneration. J. Orthop. Res. 29, 547–555. doi: 10.1002/jor.21232

PubMed Abstract | CrossRef Full Text | Google Scholar

Oshima, H., Ishihara, H., Urban, J. P. G., and Tsuji, H. (1993). The use of coccygeal discs to study intervertebral disc metabolism. J. Orthop. Res. 11, 332–338. doi: 10.1002/jor.1100110304

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 2617–2632. doi: 10.1021/acs.jced.7b00091

CrossRef Full Text | Google Scholar

Paul, C. P., Emanuel, K. S., Kingma, I., Van Der Veen, A. J., Holewijn, R. M., Vergroesen, P. P. A., et al. (2018). Changes in intervertebral disk mechanical behavior during early degeneration. J. Biomech. Eng. 140:091008.

Google Scholar

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, 2164–2171. doi: 10.1016/j.jbiomech.2004.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfannkuche, J. J., Guo, W., Cui, S., Ma, J., Lang, G., Peroglio, M., et al. (2020). Intervertebral disc organ culture for the investigation of disc pathology and regeneration–benefits, limitations, and future directions of bioreactors. Connect. Tissue Res. 61, 304–321. doi: 10.1080/03008207.2019.1665652

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, S., Menage, J., Sivan, S., and Urban, J. P. (2008). Bovine explant model of degeneration of the intervertebral disc. BMC Musculoskelet. Disord. 9:24.

Google Scholar

Robinson, R. A., and Stokes, R. H. (1949). Tables of osmotic and activity coefficients of electrolytes in aqueous solution at 25 C. Trans. Faraday Soc. 45, 612–624. doi: 10.1039/tf9494500612

CrossRef Full Text | Google Scholar

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, 2484–2490. doi: 10.1016/j.jbiomech.2005.07.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, K., Kikuchi, S., and Yonezawa, T. (1999). In vivo intradiscal pressure measurement in healthy individuals and in patients with ongoing back problems. Spine 24, 2468–2474. doi: 10.1097/00007632-199912010-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, H., Heuer, F., Drumm, J., Klezl, Z., Claes, L., and Wilke, H. J. (2007a). Application of a calibration method provides more realistic results for a finite element model of a lumbar spinal segment. Clin. Biomech. 22, 377–384. doi: 10.1016/j.clinbiomech.2006.11.008

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. (2007b). 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

CrossRef Full Text | Google Scholar

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, 310–319. doi: 10.1111/j.1469-7580.2009.01197.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, J. S., Hampson, W. G., and Jayson, M. I. (1978). The distribution of surface strain in the cadaveric lumbar spine. J. Bone Joint Surg. Br. 60, 246–251. doi: 10.1302/0301-620x.60b2.659474

PubMed Abstract | CrossRef Full Text | Google Scholar

Shapiro, E. M., Borthakur, A., Gougoutas, A., and Reddy, R. (2002). 23Na MRI accurately measures fixed charge density in articular cartilage. Magn. Reson. Med. 47, 284–291. doi: 10.1002/mrm.10054

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 3956–3963. doi: 10.1529/biophysj.107.124602

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirazi-Adl, A. (1992). Finite-element simulation of changes in the fluid content of human lumbar discs. mechanical and clinical implications. Spine 17, 206–212. doi: 10.1097/00007632-199202000-00015

CrossRef Full Text | Google Scholar

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, 120–134. doi: 10.1097/00007632-198403000-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Showalter, B. L., Beckstein, J. C., Martin, J. T., Beattie, E. E., Orías, A. A. E., Schaer, T. P., et al. (2012). Comparison of animal discs used in disc research to human lumbar disc: torsion mechanics and collagen content. Spine 37, E900–E907.

Google Scholar

Skaggs, D. L., Weidenbaum, M., Iatridis, J. C., Ratcliffe, A., and Mow, V. C. (1994). Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus. Spine 19, 1310–1319. doi: 10.1097/00007632-199406000-00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Spera, D., Genovese, K., and Voloshin, A. (2011). Application of stereo-digital image correlation to full-field 3-D deformation measurement of intervertebral disc. Strain 47, e572–e587.

Google Scholar

Steffen, T., Baramki, H. G., Rubin, R., Antoniou, J., and Aebi, M. (1998). Lumbar intradiscal pressure measured in the anterior and posterolateral annular regions during asymmetrical loading. Clin. Biomech. 13, 495–505. doi: 10.1016/s0268-0033(98)00039-4

CrossRef Full Text | Google Scholar

Stokes, I. A. (1987). Surface strain on human intervertebral discs. J. Orthop. Res. 5, 348–355. doi: 10.1002/jor.1100050306

PubMed Abstract | CrossRef Full Text | Google Scholar

Svensson, R. B., Smith, S. T., Moyer, P. J., and Magnusson, S. P. (2018). Effects of maturation and advanced glycation on tensile mechanics of collagen fibrils from rat tail and Achilles tendons. Acta Biomaterialia 70, 270–280. doi: 10.1016/j.actbio.2018.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Urban, J. P. G., and Maroudas, A. (1979). The measurement of fixed charged density in the intervertebral disc. Biochimica et Biophysica Acta 586, 166–178. doi: 10.1016/0304-4165(79)90415-x

CrossRef Full Text | Google Scholar

Urban, J. P., and McMullin, J. F. (1988). Swelling pressure of the lumbar intervertebral discs: influence of age, spinal level, composition, and degeneration. Spine 13, 179–187. doi: 10.1097/00007632-198802000-00009

PubMed Abstract | CrossRef Full Text | Google Scholar

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, 697–702. doi: 10.1002/mabi.200600063

PubMed Abstract | CrossRef Full Text | Google Scholar

van Rijsbergen, M., van Rietbergen, B., Barthelemy, V., Eltes, P., Lazáry, Á, Lacroix, D., et al. (2018). Comparison of patient-specific computational models vs. clinical follow-up, for adjacent segment disc degeneration and bone remodelling after spinal fusion. PloS One 13:e0200899. doi: 10.1371/journal.pone.0200899

PubMed Abstract | CrossRef Full Text | Google Scholar

Vergari, C., Chan, D., Clarke, A., Mansfield, J. C., Meakin, J. R., and Winlove, P. C. (2017). Bovine and degenerated human annulus fibrosus: a microstructural and micromechanical comparison. Biomech. Model. Mechanobiol. 16, 1475–1484. doi: 10.1007/s10237-017-0900-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Vernon-Roberts, B., Moore, R. J., and Fraser, R. D. (2007). The natural history of age-related disc degeneration: the pathology and sequelae of tears. Spine 32, 2797–2804. doi: 10.1097/brs.0b013e31815b64d2

PubMed Abstract | CrossRef Full Text | Google Scholar

Walter, B. A., Korecki, C. L., Purmessur, D., Roughley, P. J., Michalek, A. J., and Iatridis, J. C. (2011). Complex loading affects intervertebral disc mechanics and biology. Osteoarthr. Cartil. 19, 1011–1018. doi: 10.1016/j.joca.2011.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Xia, Q., Passias, P., Wood, K., and Li, G. (2009). Measurement of geometric deformation of lumbar intervertebral discs under in-vivo weightbearing condition. J. Biomech. 42, 705–711. doi: 10.1016/j.jbiomech.2009.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Werbner, B., Spack, K., and O’Connell, G. D. (2019). Bovine annulus fibrosus hydration affects rate-dependent failure mechanics in tension. J. Biomech. 89, 34–39. doi: 10.1016/j.jbiomech.2019.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Werbner, B., Zhou, M., and O’Connell, G. D. (2017). A novel method for repeatable failure testing of annulus fibrosus. J. Biomech. Eng. 139:111001.

Google Scholar

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, 1363–1372. doi: 10.1007/s00586-016-4416-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, H. J., Neef, P., Caimi, M., Hoogland, T., and Claes, L. E. (1999). New in vivo measurements of pressures in the intervertebral disc in daily life. Spine 24, 755–762. doi: 10.1097/00007632-199904150-00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, W., Huyghe, J. M., and Van Donkelaar, C. C. (2007). Depth-dependent compressive equilibrium properties of articular cartilage explained by its composition. Biomech. Model. Mechanobiol. 6, 43–53. doi: 10.1007/s10237-006-0044-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wognum, S., Huyghe, J. M., and Baaijens, F. P. (2006). Influence of osmotic pressure changes on the opening of existing cracks in 2 intervertebral disc models. Spine 31, 1783–1788. doi: 10.1097/01.brs.0000227267.42924.bb

CrossRef Full Text | Google Scholar

Wu, Y., Cisewski, S. E., Sachs, B. L., Pellegrini, V. D., Kern, M. J., Slate, E. H., et al. (2015). The region-dependent biomechanical and biochemical properties of bovine cartilaginous endplate. J. Biomech. 48, 3185–3191. doi: 10.1016/j.jbiomech.2015.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wuertz, K., Urban, J. P. G., Klasen, J., Ignatius, A., Wilke, H. J., Claes, L., et al. (2007). Influence of extracellular osmolarity and mechanical stimulation on gene expression of intervertebral disc cells. J. Orthop. Res. 25, 1513–1522. doi: 10.1002/jor.20436

PubMed Abstract | 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 Biomaterialia 100, 61–74. doi: 10.1016/j.actbio.2019.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, L., and Elliott, D. M. (2005). A homogenization model of the annulus fibrosus. J. Biomech. 38, 1674–1684. doi: 10.1016/j.jbiomech.2004.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Peter, C., Roberts, S., and Urban, J. P. (2002). Elastic fibre organization in the intervertebral discs of the bovine tail. J. Anat. 201, 465–475. doi: 10.1046/j.1469-7580.2002.00111.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Tirlapur, U., Fairbank, J., Handford, P., Roberts, S., Winlove, C. P., et al. (2007). Microfibrils, elastin fibres and collagen fibres in the human intervertebral disc and bovine tail disc. J. Anat. 210, 460–471. doi: 10.1111/j.1469-7580.2007.00707.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, M., Bezci, S. E., and O’Connell, G. D. (2020a). Multiscale composite model of fiber-reinforced tissues with direct representation of sub-tissue properties. Biomech. Model Mechanobiol. 19, 745–759. doi: 10.1007/s10237-019-01246-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, M., Werbner, B., and O’Connell, G. D. (2020b). Fiber engagement accounts for geometry-dependent annulus fibrosus mechanics: a multiscale, Structure-Based Finite Element Study. J. Mech. Behav. Biomed. Mater. 115:104292. doi: 10.1016/j.jmbbm.2020.104292

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, M., Werbner, B., and O’Connell, G. D. (2020c). Historical review on combined experimental and computational approaches for investigating annulus fibrosus mechanics. J. Biomech. Eng. 142:030802.

Google Scholar

Zou, J., Yang, H., Miyazaki, M., Morishita, Y., Wei, F., McGovern, S., et al. (2009). Dynamic bulging of intervertebral discs in the degenerative lumbar spine. Spine 34, 2545–2550. doi: 10.1097/brs.0b013e3181b32998

PubMed Abstract | CrossRef Full Text | Google Scholar