ORIGINAL RESEARCH article

Front. Bioeng. Biotechnol., 07 June 2021

Sec. Biomechanics

Volume 9 - 2021 | https://doi.org/10.3389/fbioe.2021.685799

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

  • 1. Berkeley Biomechanics Laboratory, Department of Mechanical Engineering, University of California, Berkeley, Berkeley, CA, United States

  • 2. Department of Orthopaedic Surgery, University of California, San Francisco, San Francisco, CA, United States

Abstract

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

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), k0 represented hydraulic permeability in the reference configuration, φ0 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 (D0) 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 (D0,Na+ = 0.00116mm2/s; D0,Cl− = 0.00161mm2/s; DNa+ = 0.00044mm2/s; DCl− = 0.00069 mm2/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

NPAF
CEP
MatrixFibers
Fluid phaseφ00.2aSee Figure 1Ba0.4c,*
k0 × 10–16 (m4/Ns)5.5b64b64b5.6c,*
M1.92c,*4.8c,*4.8c,*3.79c,*
Solid phaseE (MPa)0.4b0.74b0.74b0.31g
ν0.24d0.16c,*0.16c,*0.18c,*
β0.95c,*3.3c,*3.3c,*0.29c,*
Elin: (MPa)N.A.N.A.600eN.A.
γN.A.N.A.5.95f,*N.A.
λ0N.A.N.A.1.05eN.A.

Triphasic material properties of the bovine caudal disc tissues.

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

* The parameter was determined based on experimental studies using matching human intervertebral disc tissues due to the lack of corresponding data obtained from bovine caudal disc tissues.

aBeckstein et al. (2008).

bPérié et al. (2005).

cCortes et al. (2014).

dFarrell and Riches (2013).

eFratzl et al. (1998), Gentleman et al. (2003), van der Rijt et al. (2006), and Shen et al. (2008).

fZhou et al. (2020a).

gWu et al. (2015).

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, I1 and I2 represented the first and second invariants of the right Cauchy-Green deformation tensor, C(C = FTF), 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, Elin represented the fiber modulus in the fiber linear region, and λ0 represented the transition stretch between the toe and linear regions (Holzapfel and Ogden, 2017). B was a function of γ, Elin., and λ0 (. 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: Elin, γ, and λ0). 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–“∗”).

Bony endplates were modeled as a compressible hyperelastic material using the Neo-Hookean description (Eq. 6). I1, I2, J were defined as above.Ebonyendplates 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).

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

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 CaseA, 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

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

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

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 7

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

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

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

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

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.

Statements

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

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

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.

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

  • 1

    AdamC.RouchP.SkalliW. (2015). Inter-lamellar shear resistance confers compressive stiffness in the intervertebral disc: an image-based modelling study on the bovine caudal disc.J. Biomech.484303–4308. 10.1016/j.jbiomech.2015.10.041

  • 2

    AdamsM. A.GreenT. P. (1993). Tensile properties of the annulus fibrosus.Eur. Spine J.2203–208. 10.1007/bf00299447

  • 3

    AdamsM. A.HuttonW. C. (1985). Gradual disc prolapse.Spine10524–531. 10.1097/00007632-198507000-00006

  • 4

    AdamsM. A.RoughleyP. J. (2006). What is intervertebral disc degeneration, and what causes it?Spine312151–2161. 10.1097/01.brs.0000231761.73859.2c

  • 5

    AdamsM. A.FreemanB. J.MorrisonH. P.NelsonI. W.DolanP. (2000). Mechanical initiation of intervertebral disc degeneration.Spine251625–1636. 10.1097/00007632-200007010-00005

  • 6

    AliniM.EisensteinS. M.ItoK.LittleC.KettlerA. A.MasudaK.et al (2008). Are animal models useful for studying human disc disorders/degeneration?Eur. Spine J.172–19. 10.1007/s00586-007-0414-y

  • 7

    ArmstrongC. G.MowV. C. (1982). Variations in the intrinsic mechanical properties of human articular cartilage with age, degeneration, and water content.J. Bone Joint Surg. Am.6488–94. 10.2106/00004623-198264010-00013

  • 8

    AteshianG. A.ChahineN. O.BasaloI. M.HungC. T. (2004). The correspondence between equilibrium biphasic and triphasic material properties in mixture models of articular cartilage.J. Biomech.37391–400. 10.1016/s0021-9290(03)00252-5

  • 9

    AteshianG. A.LaiW. M.ZhuW. B.MowV. C. (1994). An asymptotic solution for the contact of two biphasic cartilage layers.J. Biomech.271347–1360. 10.1016/0021-9290(94)90044-2

  • 10

    BarthelemyV. M. P.Van RijsbergenM. M.WilsonW.HuygheJ. M.Van RietbergenB.ItoK. (2016). A computational spinal motion segment model incorporating a matrix composition-based model of the intervertebral disc.J. Mech. Behav. Biomed. Mater.54194–204. 10.1016/j.jmbbm.2015.09.028

  • 11

    BecksteinJ. C.SenS.SchaerT. P.VresilovicE. J.ElliottD. M. (2008). Comparison of animal discs used in disc research to human lumbar disc: axial compression mechanics and glycosaminoglycan content.Spine33E166–E173.

  • 12

    Berger-RoscherN.CasaroliG.RascheV.VillaT.GalbuseraF.WilkeH. J. (2017). Influence of complex loading conditions on intervertebral disc failure.Spine42E78–E85.

  • 13

    Berg-JohansenB.HanM.FieldsA. J.LiebenbergE. C.LimB. J.LarsonP. E.et al (2018). Cartilage endplate thickness variation measured by ultrashort echo-time MRI is associated with adjacent disc degeneration.Spine43E592–E600.

  • 14

    BezciS. E.KlinebergE. O.O’ConnellG. 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. 10.1016/j.jmbbm.2017.09.022

  • 15

    BezciS. E.O’ConnellG. D. (2018). Osmotic pressure alters time-dependent recovery behavior of the intervertebral disc.Spine43E334–E340.

  • 16

    BezciS. E.LimS.O’ConnellG. D. (2020a). Nonlinear stress-dependent recovery behavior of the intervertebral disc.J. Mech. Behav. Biomed. Mater.110:103881. 10.1016/j.jmbbm.2020.103881

  • 17

    BezciS. E.NandyA.O’ConnellG. D. (2015). Effect of hydration on healthy intervertebral disk mechanical stiffness.J. Biomech. Eng.137:101007.

  • 18

    BezciS. E.TorresK.CarraroC.ChiavacciD.WerbnerB.LimS.et al (2020b). Transient swelling behavior of the bovine caudal disc.J. Mech. Behav. Biomed. Mater.112:104089. 10.1016/j.jmbbm.2020.104089

  • 19

    BezciS. E.WerbnerB.ZhouM.MalollariK. G.DorlhiacG.CarraroC.et al (2019). Radial variation in biochemical composition of the bovine caudal intervertebral disc.JOR Spine2:e1065.

  • 20

    CastroA. P. G.AlvesJ. L. (2020). Numerical implementation of an osmo-poro-visco-hyperelastic finite element solver: application to the intervertebral disc.Comput. Methods Biomech. Biomed. Eng.281–13. 10.1080/10255842.2020.1839059

  • 21

    CastroA. P. G.PaulC. P. L.DetigerS. E. L.SmitT. H.Van RoyenB. J.Pimenta ClaroJ. C.et al (2014). Long-term creep behavior of the intervertebral disk: comparison between bioreactor data and numerical results.Front. Bioeng. Biotechnol.2:56.

  • 22

    ChoiK.KuhnJ. L.CiarelliM. J.GoldsteinS. A. (1990). The elastic moduli of human subchondral, trabecular, and cortical bone tissue and the size-dependency of cortical bone modulus.J. Biomech.231103–1113. 10.1016/0021-9290(90)90003-l

  • 23

    CortesD. H.JacobsN. T.DeLuccaJ. F.ElliottD. M. (2014). Elastic, permeability and swelling properties of human intervertebral disc tissues: a benchmark for tissue engineering.J. Biomech.472088–2094. 10.1016/j.jbiomech.2013.12.021

  • 24

    CostiJ. J.LedetE. H.O’ConnellG. D. (2020). Spine biomechanical testing methodologies: the controversy of consensus vs scientific evidence.JOR Spine4:e1138.

  • 25

    DemersC. N.AntoniouJ.MwaleF. (2004). Value and limitations of using the bovine tail as a model for the human lumbar spine.Spine292793–2799. 10.1097/01.brs.0000147744.74215.b0

  • 26

    DreischarfM.ZanderT.Shirazi-AdlA.PuttlitzC. M.AdamC. J.ChenC. 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.471757–1766. 10.1016/j.jbiomech.2014.04.002

  • 27

    EskandariM.NordgrenT. M.O’ConnellG. D. (2019). Mechanics of pulmonary airways: linking structure to function through constitutive modeling, biochemistry, and histology.Acta Biomaterialia97513–523. 10.1016/j.actbio.2019.07.020

  • 28

    FarrellM. D.RichesP. E. (2013). On the poisson’s ratio of the nucleus pulposus.J. Biomech. Eng.135104501–104504.

  • 29

    FratzlP.MisofK.ZizakI.RappG.AmenitschH.BernstorffS. (1998). Fibrillar structure and mechanical properties of collagen.J. Struct. Biol.122119–122. 10.1006/jsbi.1998.3966

  • 30

    FrauchigerD. A.ChanS. C.BennekerL. M.GantenbeinB. (2018). Intervertebral disc damage models in organ culture: a comparison of annulus fibrosus cross-incision versus punch model under complex loading.Eur. Spine J.271785–1797. 10.1007/s00586-018-5638-5

  • 31

    GalbuseraF.SchmidtH.Neidlinger-WilkeC.GottschalkA.WilkeH. 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.20563–571. 10.1007/s00586-010-1586-4

  • 32

    GalbuseraF.SchmidtH.Neidlinger-WilkeC.WilkeH. 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.14729–739. 10.1080/10255842.2010.493522

  • 33

    GalbuseraF.Van RijsbergenM.ItoK.HuygheJ. M.Brayda-BrunoM.WilkeH. J. (2014). Ageing and degenerative changes of the intervertebral disc and their impact on spinal flexibility.Eur. Spine J.23324–332.

  • 34

    GentlemanE.LayA. N.DickersonD. A.NaumanE. A.LivesayG. A.DeeK. C. (2003). Mechanical characterization of collagen fibers and scaffolds for tissue engineering.Biomaterials243805–3813. 10.1016/s0142-9612(03)00206-0

  • 35

    GoelV. K.MonroeB. T.GilbertsonL. G.BrinckmannP. (1995a). Interlaminar shear stresses and laminae separation in a disc: finite element analysis of the L3-L4 motion segment subjected to axial compressive loads.Spine20689–698. 10.1097/00007632-199503150-00010

  • 36

    GoelV. K.RamirezS. A.KongW.GilbertsonL. 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.117266–271. 10.1115/1.2794180

  • 37

    GuW. Y.YaoH.VegaA. L.FlaglerD. (2004). Diffusivity of ions in agarose gels and intervertebral disc: effect of porosity.Ann. Biomed. Eng.32, 1710–1717. 10.1007/s10439-004-7823-4

  • 38

    HeuerF.SchmidtH.WilkeH. J. (2008a). Stepwise reduction of functional spinal structures increase disc bulge and surface strains.J. Biomech.411953–1960. 10.1016/j.jbiomech.2008.03.023

  • 39

    HeuerF.SchmidtH.WilkeH. J. (2008b). The relation between intervertebral disc bulging and annular fiber associated strains for simple and complex loading.J. Biomech.411086–1094. 10.1016/j.jbiomech.2007.11.019

  • 40

    HeuerF.SchmidtH.KäferW.GrafN.WilkeH. J. (2012). Posterior motion preserving implants evaluated by means of intervertebral disc bulging and annular fiber strains.Clin. Biomech.27218–225. 10.1016/j.clinbiomech.2011.09.004

  • 41

    HolzapfelG. A.OgdenR. W. (2017). Biomechanics: Trends in Modeling and Simulation, Vol. 316. Berlin: Springer.

  • 42

    HolzapfelG. A.Schulze-BauerC. A. J.FeiglG.RegitnigP. (2005). Single lamellar mechanics of the human lumbar anulus fibrosus.Biomech. Model. Mechanobiol.3125–140. 10.1007/s10237-004-0053-8

  • 43

    IatridisJ. C.MacLeanJ. J.RyanD. A. (2005). Mechanical damage to the intervertebral disc annulus fibrosus subjected to tensile loading.J. Biomech.38557–565. 10.1016/j.jbiomech.2004.03.038

  • 44

    IsaacsJ. L.VresilovicE.SarkarS.MarcolongoM. (2014). Role of biomolecules on annulus fibrosus micromechanics: effect of enzymatic digestion on elastic and failure properties.J. Mech. Behav. Biomed. Mater.4075–84. 10.1016/j.jmbbm.2014.08.012

  • 45

    IshiharaH.McNallyD. S.UrbanJ. P.HallA. C. (1996). Effects of hydrostatic pressure on matrix synthesis in different regions of the intervertebral disk.J. Appl. Physiol.80839–846. 10.1152/jappl.1996.80.3.839

  • 46

    KimY. E.GoelV. K.WeinsteinJ. N.LimT. H. (1991). Effect of disc degeneration at one level on the adjacent level in axial mode.Spine16331–335. 10.1097/00007632-199103000-00013

  • 47

    KoreckiC. L.CostiJ. J.IatridisJ. C. (2008a). Needle puncture injury affects intervertebral disc mechanics and biology in an organ culture model.Spine33235–241. 10.1097/brs.0b013e3181624504

  • 48

    KoreckiC. L.MacLeanJ. J.IatridisJ. C. (2008b). Dynamic compression effects on intervertebral disc mechanics and biology.Spine331403–1409. 10.1097/brs.0b013e318175cae7

  • 49

    KurowskiP.KuboA. I. Z. O. H. (1986). The relationship of degeneration of the intervertebral disc to mechanical loading conditions on lumbar vertebrae.Spine11726–731. 10.1097/00007632-198609000-00012

  • 50

    LaiW. M.HouJ. S.MowV. C. (1991). A triphasic theory for the swelling and deformation behaviors of articular cartilage.J. Biomech. Eng.113245–258. 10.1115/1.2894880

  • 51

    LiY.FesselG.GeorgiadisM.SnedekerJ. G. (2013). Advanced glycation end-products diminish tendon collagen fiber sliding.Matrix Biol.32169–177. 10.1016/j.matbio.2013.01.003

  • 52

    LiuQ.YangX. P.LiK.YangT.YeJ. D.ZhangC. Q. (2017). Internal strains of anulus fibrosus in the intervertebral disc under axial compression load.Biomed. Res.83483–3486.

  • 53

    LüsseaS.ClaassenH.GehrkeT.HassenpflugJ.SchünkeM.HellerM.et al (2000). Evaluation of water content by spatially resolved transverse relaxation times of human articular cartilage.Magn. Reson. Imaging18423–430. 10.1016/s0730-725x(99)00144-7

  • 54

    MaasS. A.EllisB. J.AteshianG. A.WeissJ. A. (2012). FEBio: finite elements for biomechanics.J. Biomech. Eng.134:011005.

  • 55

    MalandrinoA.NoaillyJ.LacroixD. (2013). Regional annulus fibre orientations used as a tool for the calibration of lumbar intervertebral disc finite element models.Comput. Methods Biomech. Biomed. Eng.16923–928. 10.1080/10255842.2011.644539

  • 56

    MarchandF. R.AhmedA. M. (1990). Investigation of the laminate structure of lumbar disc anulus fibrosus.Spine15402–410. 10.1097/00007632-199005000-00011

  • 57

    MaroudasA.MuirH.WinghamJ. (1969). The correlation of fixed negative charge with glycosaminoglycan content of human articular cartilage.Biochimica et Biophysica Acta177492–500. 10.1016/0304-4165(69)90311-0

  • 58

    MartinJ. T.MilbyA. H.ChiaroJ. A.KimD. H.HebelaN. M.SmithL. 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.102473–2481. 10.1016/j.actbio.2014.02.024

  • 59

    MatcherS. J.WinloveC. P.GangnusS. V. (2004). The collagen structure of bovine intervertebral disc studied using polarization-sensitive optical coherence tomography.Phys. Med. Biol.491295–1306. 10.1088/0031-9155/49/7/016

  • 60

    MichalekA. J.IatridisJ. C. (2012). Height and torsional stiffness are most sensitive to annular injury in large animal intervertebral discs.Spine J.12425–432. 10.1016/j.spinee.2012.04.001

  • 61

    MichalekA. J.BuckleyM. R.BonassarL. J.CohenI.IatridisJ. 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.422279–2285. 10.1016/j.jbiomech.2009.06.047

  • 62

    MimuraM.PanjabiM. M.OxlandT. R.CriscoJ. J.YamamotoI.VasavadaA. (1994). Disc degeneration affects the multidirectional flexibility of the lumbar spine.Spine191371–1380. 10.1097/00007632-199406000-00011

  • 63

    MonacoL. A.DeWitte-OrrS. J.GregoryD. E. (2016). A comparison between porcine, ovine, and bovine intervertebral disc anatomy and single lamella annulus fibrosus tensile properties.J. Morphol.277244–251. 10.1002/jmor.20492

  • 64

    NatarajanR. N.WilliamsJ. R.AnderssonG. B. (2006). Modeling changes in intervertebral disc mechanics with degeneration.J. Bone Joint Surg. Am.88(Suppl._2), 36–40. 10.2106/jbjs.f.00002

  • 65

    NewellN.Rivera TapiaD.RahmanT.LimS.O’ConnellG. D.HolsgroveT. P. (2020). Influence of testing environment and loading rate on intervertebral disc compressive mechanics: an assessment of repeatability at three different laboratories.JOR Spine3:e21110.

  • 66

    NguyenA. M.JohannessenW.YoderJ. H.WheatonA. J.VresilovicE. J.BorthakurA.et al (2008). Noninvasive quantification of human nucleus pulposus pressure with use of T1ρ-weighted magnetic resonance imaging.J. Bone Joint Surg. Am.90796–802. 10.2106/jbjs.g.00667

  • 67

    O’ConnellG. D.JohannessenW.VresilovicE. J.ElliottD. M. (2007a). Human internal disc strains in axial compression measured noninvasively using magnetic resonance imaging.Spine322860–2868. 10.1097/brs.0b013e31815b75fb

  • 68

    O’ConnellG. D.VresilovicE. J.ElliottD. M. (2007b). Comparison of animals used in disc research to human lumbar disc geometry.Spine32328–333. 10.1097/01.brs.0000253961.40910.c1

  • 69

    O’ConnellG. D.VresilovicE. J.ElliottD. M. (2011). Human intervertebral disc internal strain in compression: the effect of disc region, loading position, and degeneration.J. Orthop. Res.29547–555. 10.1002/jor.21232

  • 70

    OshimaH.IshiharaH.UrbanJ. P. G.TsujiH. (1993). The use of coccygeal discs to study intervertebral disc metabolism.J. Orthop. Res.11332–338. 10.1002/jor.1100110304

  • 71

    PartanenJ. I.PartanenL. J.VahteristoK. 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. Data622617–2632. 10.1021/acs.jced.7b00091

  • 72

    PaulC. P.EmanuelK. S.KingmaI.Van Der VeenA. J.HolewijnR. M.VergroesenP. P. A.et al (2018). Changes in intervertebral disk mechanical behavior during early degeneration.J. Biomech. Eng.140:091008.

  • 73

    PériéD.KordaD.IatridisJ. 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.382164–2171. 10.1016/j.jbiomech.2004.10.002

  • 74

    PfannkucheJ. J.GuoW.CuiS.MaJ.LangG.PeroglioM.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.61304–321. 10.1080/03008207.2019.1665652

  • 75

    RobertsS.MenageJ.SivanS.UrbanJ. P. (2008). Bovine explant model of degeneration of the intervertebral disc.BMC Musculoskelet. Disord.9:24.

  • 76

    RobinsonR. A.StokesR. H. (1949). Tables of osmotic and activity coefficients of electrolytes in aqueous solution at 25 C.Trans. Faraday Soc.45612–624. 10.1039/tf9494500612

  • 77

    RohlmannA.ZanderT.SchmidtH.WilkeH. J.BergmannG. (2006). Analysis of the influence of disc degeneration on the mechanical behaviour of a lumbar motion segment using the finite element method.J. Biomech.392484–2490. 10.1016/j.jbiomech.2005.07.026

  • 78

    SatoK.KikuchiS.YonezawaT. (1999). In vivo intradiscal pressure measurement in healthy individuals and in patients with ongoing back problems.Spine242468–2474. 10.1097/00007632-199912010-00008

  • 79

    SchmidtH.GalbuseraF.RohlmannA.Shirazi-AdlS. A. (2013). What have we learned from finite element model studies of lumbar intervertebral discs in the past four decades?J. Biomech.462342–2355. 10.1016/j.jbiomech.2013.07.014

  • 80

    SchmidtH.HeuerF.DrummJ.KlezlZ.ClaesL.WilkeH. J. (2007a). Application of a calibration method provides more realistic results for a finite element model of a lumbar spinal segment.Clin. Biomech.22377–384. 10.1016/j.clinbiomech.2006.11.008

  • 81

    SchmidtH.HeuerF.SimonU.KettlerA.RohlmannA.ClaesL.et al (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus.Clin. Biomech.21337–344. 10.1016/j.clinbiomech.2005.12.001

  • 82

    SchmidtH.KettlerA.HeuerF.SimonU.ClaesL.WilkeH. J. (2007b). Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading.Spine32748–755. 10.1097/01.brs.0000259059.90430.c2

  • 83

    SchollumM. L.RobertsonP. A.BroomN. D. (2010). How age influences unravelling morphology of annular lamellae–a study of interfibre cohesivity in the lumbar disc.J. Anat.216310–319. 10.1111/j.1469-7580.2009.01197.x

  • 84

    ShahJ. S.HampsonW. G.JaysonM. I. (1978). The distribution of surface strain in the cadaveric lumbar spine.J. Bone Joint Surg. Br.60246–251. 10.1302/0301-620x.60b2.659474

  • 85

    ShapiroE. M.BorthakurA.GougoutasA.ReddyR. (2002). 23Na MRI accurately measures fixed charge density in articular cartilage.Magn. Reson. Med.47284–291. 10.1002/mrm.10054

  • 86

    ShenZ. L.DodgeM. R.KahnH.BallariniR.EppellS. J. (2008). Stress-strain experiments on individual collagen fibrils.Biophys. J.953956–3963. 10.1529/biophysj.107.124602

  • 87

    Shirazi-AdlA. (1992). Finite-element simulation of changes in the fluid content of human lumbar discs. mechanical and clinical implications.Spine17206–212. 10.1097/00007632-199202000-00015

  • 88

    Shirazi-AdlS. A.ShrivastavaS. C.AhmedA. M. (1984). Stress analysis of the lumbar disc-body unit in compression. a three-dimensional nonlinear finite element study.Spine9120–134. 10.1097/00007632-198403000-00003

  • 89

    ShowalterB. L.BecksteinJ. C.MartinJ. T.BeattieE. E.OríasA. A. E.SchaerT. P.et al (2012). Comparison of animal discs used in disc research to human lumbar disc: torsion mechanics and collagen content.Spine37E900–E907.

  • 90

    SkaggsD. L.WeidenbaumM.IatridisJ. C.RatcliffeA.MowV. C. (1994). Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus.Spine191310–1319. 10.1097/00007632-199406000-00002

  • 91

    SperaD.GenoveseK.VoloshinA. (2011). Application of stereo-digital image correlation to full-field 3-D deformation measurement of intervertebral disc.Strain47e572–e587.

  • 92

    SteffenT.BaramkiH. G.RubinR.AntoniouJ.AebiM. (1998). Lumbar intradiscal pressure measured in the anterior and posterolateral annular regions during asymmetrical loading.Clin. Biomech.13495–505. 10.1016/s0268-0033(98)00039-4

  • 93

    StokesI. A. (1987). Surface strain on human intervertebral discs.J. Orthop. Res.5348–355. 10.1002/jor.1100050306

  • 94

    SvenssonR. B.SmithS. T.MoyerP. J.MagnussonS. P. (2018). Effects of maturation and advanced glycation on tensile mechanics of collagen fibrils from rat tail and Achilles tendons.Acta Biomaterialia70270–280. 10.1016/j.actbio.2018.02.005

  • 95

    UrbanJ. P. G.MaroudasA. (1979). The measurement of fixed charged density in the intervertebral disc.Biochimica et Biophysica Acta586166–178. 10.1016/0304-4165(79)90415-x

  • 96

    UrbanJ. P.McMullinJ. F. (1988). Swelling pressure of the lumbar intervertebral discs: influence of age, spinal level, composition, and degeneration.Spine13179–187. 10.1097/00007632-198802000-00009

  • 97

    van der RijtJ. A.Van Der WerfK. O.BenninkM. L.DijkstraP. J.FeijenJ. (2006). Micromechanical testing of individual collagen fibrils.Macromol. Biosci.6697–702. 10.1002/mabi.200600063

  • 98

    van RijsbergenM.van RietbergenB.BarthelemyV.EltesP.LazáryÁLacroixD.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 One13:e0200899. 10.1371/journal.pone.0200899

  • 99

    VergariC.ChanD.ClarkeA.MansfieldJ. C.MeakinJ. R.WinloveP. C. (2017). Bovine and degenerated human annulus fibrosus: a microstructural and micromechanical comparison.Biomech. Model. Mechanobiol.161475–1484. 10.1007/s10237-017-0900-z

  • 100

    Vernon-RobertsB.MooreR. J.FraserR. D. (2007). The natural history of age-related disc degeneration: the pathology and sequelae of tears.Spine322797–2804. 10.1097/brs.0b013e31815b64d2

  • 101

    WalterB. A.KoreckiC. L.PurmessurD.RoughleyP. J.MichalekA. J.IatridisJ. C. (2011). Complex loading affects intervertebral disc mechanics and biology.Osteoarthr. Cartil.191011–1018. 10.1016/j.joca.2011.04.005

  • 102

    WangS.XiaQ.PassiasP.WoodK.LiG. (2009). Measurement of geometric deformation of lumbar intervertebral discs under in-vivo weightbearing condition.J. Biomech.42705–711. 10.1016/j.jbiomech.2009.01.004

  • 103

    WerbnerB.SpackK.O’ConnellG. D. (2019). Bovine annulus fibrosus hydration affects rate-dependent failure mechanics in tension.J. Biomech.8934–39. 10.1016/j.jbiomech.2019.04.008

  • 104

    WerbnerB.ZhouM.O’ConnellG. D. (2017). A novel method for repeatable failure testing of annulus fibrosus.J. Biomech. Eng.139:111001.

  • 105

    WilkeH. J.KienleA.MaileS.RascheV.Berger-RoscherN. (2016). A new dynamic six degrees of freedom disc-loading simulator allows to provoke disc damage and herniation.Eur. Spine J.251363–1372. 10.1007/s00586-016-4416-5

  • 106

    WilkeH. J.NeefP.CaimiM.HooglandT.ClaesL. E. (1999). New in vivo measurements of pressures in the intervertebral disc in daily life.Spine24755–762. 10.1097/00007632-199904150-00005

  • 107

    WilsonW.HuygheJ. M.Van DonkelaarC. C. (2007). Depth-dependent compressive equilibrium properties of articular cartilage explained by its composition.Biomech. Model. Mechanobiol.643–53. 10.1007/s10237-006-0044-z

  • 108

    WognumS.HuygheJ. M.BaaijensF. P. (2006). Influence of osmotic pressure changes on the opening of existing cracks in 2 intervertebral disc models.Spine311783–1788. 10.1097/01.brs.0000227267.42924.bb

  • 109

    WuY.CisewskiS. E.SachsB. L.PellegriniV. D.KernM. J.SlateE. H.et al (2015). The region-dependent biomechanical and biochemical properties of bovine cartilaginous endplate.J. Biomech.483185–3191. 10.1016/j.jbiomech.2015.07.005

  • 110

    WuertzK.UrbanJ. P. G.KlasenJ.IgnatiusA.WilkeH. J.ClaesL.et al (2007). Influence of extracellular osmolarity and mechanical stimulation on gene expression of intervertebral disc cells.J. Orthop. Res.251513–1522. 10.1002/jor.20436

  • 111

    YangB.O’ConnellG. D. (2019). Intervertebral disc swelling maintains strain homeostasis throughout the annulus fibrosus: a finite element analysis of healthy and degenerated discs.Acta Biomaterialia10061–74. 10.1016/j.actbio.2019.09.035

  • 112

    YinL.ElliottD. M. (2005). A homogenization model of the annulus fibrosus.J. Biomech.381674–1684. 10.1016/j.jbiomech.2004.07.017

  • 113

    YuJ.PeterC.RobertsS.UrbanJ. P. (2002). Elastic fibre organization in the intervertebral discs of the bovine tail.J. Anat.201465–475. 10.1046/j.1469-7580.2002.00111.x

  • 114

    YuJ.TirlapurU.FairbankJ.HandfordP.RobertsS.WinloveC. P.et al (2007). Microfibrils, elastin fibres and collagen fibres in the human intervertebral disc and bovine tail disc.J. Anat.210460–471. 10.1111/j.1469-7580.2007.00707.x

  • 115

    ZhouM.BezciS. E.O’ConnellG. D. (2020a). Multiscale composite model of fiber-reinforced tissues with direct representation of sub-tissue properties.Biomech. Model Mechanobiol.19745–759. 10.1007/s10237-019-01246-x

  • 116

    ZhouM.WerbnerB.O’ConnellG. 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. 10.1016/j.jmbbm.2020.104292

  • 117

    ZhouM.WerbnerB.O’ConnellG. D. (2020c). Historical review on combined experimental and computational approaches for investigating annulus fibrosus mechanics.J. Biomech. Eng.142:030802.

  • 118

    ZouJ.YangH.MiyazakiM.MorishitaY.WeiF.McGovernS.et al (2009). Dynamic bulging of intervertebral discs in the degenerative lumbar spine.Spine342545–2550. 10.1097/brs.0b013e3181b32998

Summary

Keywords

finite element modeling, multiscale modeling, multiphasic modeling, structure-based modeling, structure-function relationship, bovine caudal disc, intervertebral disc degeneration

Citation

Zhou M, Lim S and O’Connell GD (2021) A Robust Multiscale and Multiphasic Structure-Based Modeling Framework for the Intervertebral Disc. Front. Bioeng. Biotechnol. 9:685799. doi: 10.3389/fbioe.2021.685799

Received

25 March 2021

Accepted

10 May 2021

Published

07 June 2021

Volume

9 - 2021

Edited by

Luca Cristofolini, University of Bologna, Italy

Reviewed by

André P. G. Castro, Universidade de Lisboa, Portugal; Diane Wagner, Indiana University, Purdue University Indianapolis, United States

Updates

Copyright

© 2021 Zhou, Lim and O’Connell.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Grace D. O’Connell,

This article was submitted to Biomechanics, a section of the journal Frontiers in Bioengineering and Biotechnology

Disclaimer

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

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Share on WeChat

Scan with WeChat to share this article

Article metrics