Torque- and Muscle-Driven Flexion Induce Disparate Risks of In Vitro Herniation: A Multiscale and Multiphasic Structure-Based Finite Element Study

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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