Multiaxial validation of a finite element model of the intervertebral disc with multigenerational fibers to establish residual strain

Abstract Finite element models of the intervertebral disc are used to address research questions that cannot be tested through typical experimentation. A disc model requires complex geometry and tissue properties to be accurately defined to mimic the physiological disc. The physiological disc possesses residual strain in the annulus fibrosus (AF) due to osmotic swelling and due to inherently pre‐strained fibers. We developed a disc model with residual contributions due to swelling‐only, and a multigeneration model with residual contributions due to both swelling and AF fiber pre‐strain and validated it against organ‐scale uniaxial, quasi‐static and multiaxial, dynamic mechanical tests. In addition, we demonstrated the models' ability to mimic the opening angle observed following radial incision of bovine discs. Both models were validated against organ‐scale experimental data. While the swelling only model responses were within the experimental 95% confidence interval, the multigeneration model offered outcomes closer to the experimental mean and had a bovine model opening angle within one SD of the experimental mean. The better outcomes for the multigeneration model, which allowed for the inclusion of inherently pre‐strained fibers in AF, is likely due to its uniform fiber contribution throughout the AF. We conclude that the residual contribution of pre‐strained fibers in the AF should be included to best simulate the physiological disc and its behaviors.


| INTRODUCTION
Finite element models of the intervertebral disc aim to incorporate its complex geometry and material properties to predict the disc's multiaxial mechanics and address research questions that cannot be tested through experimentation. The constituent relations used in a model of the disc should describe the nonlinear, anisotropic, osmotic and biphasic properties in order to effectively simulate the mechanical behavior of the physiological disc. [1][2][3][4] Our lab previously developed a disc model with such properties and validated it against uniaxial, quasi-static compression tests. 4 In our model, disc constituent material properties were obtained from individually testing the disc's tissue components. The model as a whole was then validated with organ-scale testing. The model replicated mechanical outcomes in uniaxial quasi-static slow ramp, creep, and stress relaxation tests very well 4 ; however, when subsequently applied to dynamic, multiaxial tests, the model was unable to match the nonlinear experimental responses ( Figure S1).
The disc is comprised of an inner nucleus pulposus (NP) and outer layers of annulus fibrosus (AF). The NP is rich in proteoglycans that enable it to maintain a high water content; this fluid is essential to the NP pressurization which enables the disc to withstand high compressive loads. [5][6][7][8][9] The AF also has high proteoglycan content to support compression, but importantly, it consists of fibrous, concentric layers with an angle-ply fiber structure which enable the AF to withstand circumferential tension. 5,6,[10][11][12] To the superior and inferior of the NP and AF are porous cartilage endplates (CEP) which control fluid and nutrient exchange between the disc and vertebral bodies (VB). 8,[13][14][15][16] The combination of these components enables the disc to withstand large, multi-axial, dynamic loading during regular daily activities.
Significant residual strain is present in the AF, but its inclusion in finite element models is inconsistent and potentially insufficient. 17,18 Bovine discs removed from the endplates and then cut radially, opened up as the internal residual strain is released. 18 These outcomes demonstrate that the disc's residual strain is partially released in the absence of swelling pressure and further by the relaxation and recoil of the AF fibers. Therefore, the disc possesses two types of residual strain, 1) due to swelling, and 2) due to inherently pre-strained fibers in the AF. The high proteoglycan content in the disc sets up a fixed charge density gradient that draws and holds fluid in the disc space enabling a buildup of swelling pressure. The fibers present in the AF develop an inherent pre-strain as a result of growth and development. 17,19,20 The experiments conducted on physiological discs have mechanical contribution from both of these residual strains whereas model simulations do not. The material properties used in models are acquired from tissue tests which require the extraction of individual disc components. Tissue extraction releases residual strains and leaves their contribution unaccounted for in constituent testing and in the resultant material properties used in finite element models ( Figure 1). 4,21 While many models have included residual strain contributions due to swelling, 1,2,4,[22][23][24][25][26][27][28][29] none to date have explicitly included residual contribution from pre-strained AF fibers. We hypothesize our previous model lacked sufficient residual strain which inhibited its performance in dynamic axial compression and torsion testing scenarios ( Figure S1).
We sought to update our previous model and develop a process to establish pre-strained AF fibers in the model. In order to incorporate pre-strained fibers we utilized the multigeneration feature in FEBio, which allows for constitutive relationships of the model to be established at different reference configurations. 30 The first objective of this study was to incorporate residual strain due to swelling and multigeneration fibers in a finite element model of human disc and validate it against uniaxial quasi-static and multi-axial dynamic tests. We compared a swelling-only model, with AF residual strain caused by swelling, and a multigeneration model, with AF residual strain caused by both swelling and inherent AF fiber pre-strain. A parametric analysis was conducted for the multigeneration model with several parameter combinations to investigate the impact of fiber-induced residual strain, and to determine the values best suited for simulating experimental outcomes and by extension, the physiological disc.
The second objective of this study was to simulate the bovine disc incision experiment by Michalek and coauthors to demonstrate the model's residual strain. 18 The AF opening gap following radial incision is a measure of residual strain. We simulated the radial incision experiment using both the swelling-only and multigeneration models, calculated the AF opening gap, and compared it against the previously published experimental outcomes. 18

| METHODS
We expanded upon our previous work to establish a swelling-only model and developed a novel multigeneration model which featured residual strain in the AF due to both swelling and pre-strained fiber contributions.
For both models, it is important to note that material properties were based on tissue-level experiments and were not tuned to fit organ-level experimental outcomes (Table 1). [4][5][6]13,21,[31][32][33][34] For the parametric study of the multigeneration model, we varied two parameters, axial displacement (ΔH) and twist angle (Ω), as the disc height and torsional rotation impacted the fiber pre-strain imposed and will be described in greater detail later in Section 2.3.2. Bovine and human discs possess similar water and proteoglycan content, and their mechanical responses scale with disc size; therefore, the same constitutive equations and material properties for the human AF and NP were also utilized for the bovine model and only geometric differences were included. [35][36][37][38] 2.1 | Initial geometry and mesh for human disc model The three-dimensional geometry of the model was created from the mean shape of seven human L4/L5 discs of Pfirrmann F I G U R E 1 AF tissue testing to determine material properties (red stress-strain response) is performed on excised tissue with residual strain released (blue circle); however, the intact disc is subject to significant prestrain. The disc finite element model is initiated in a stress-and strain-free reference configuration (blue circle) and the AF pre-strain can be increased by imposing residual strain due to swelling (green square) or due to both swelling and fiber pre-strain (black X); with the overarching goal of matching the disc state at the condition of interest (e.g., the start of motion segment testing, as in the tests utilized for validation herein) degeneration grade 3. 4,39 A custom Matlab script meshed the mean disc to 10,625 hexahedral elements and subsequent swelling conditions established the initial bulged shape. The disc model includes distinct constitutive models and material properties for the NP, AF, CEP, and VB ( Figure 2). The properties of the AF vary radially, therefore it was modeled as having an outer layer (OAF) and inner layer (IAF) with a transitional single element layer between them (AFtrans). Similarly, a single element transitional layer was created between the IAF and the NP (NPtrans) ( Figure 2; Table 1).

| Matrix
The AF, NP, and CEP all included extrafibrillar matrix with a Holmes-Mow constitutive equation. 40 Holmes-Mow matrix is an isotropic, hyperelastic material model which is widely used to model disc. 3,4,21,29,41 3,4,21,29,41 The strain-energy function for the Holmes Mow matrix is: where I 1 and I 2 are the right Cauchy-Green tensor invariants, J is the Jacobian of the deformation gradient, β m is the matrix stiffening coefficient, and λ and μ are the Lamé parameters which can be related to the Young's modulus E m and Poisson's ratio ν as follows: E m , ν, and β m were specified for each constituent based on tissuelevel compression tests (Table 1) F I G U R E 2 A, Sagittal view of the initial model geometry for with disc constituent labels for the outer annulus fibrosus (OAF), annulus fibrosus transitional layer (AFtrans), inner annulus fibrosus (IAF), nucleus pulposus transitional layer (NPtrans), nucleus pulposus (NP), cartilage endplates (CEP) and vertebral bodies (VB). B, Swelling-only model after swelling at fixed height (H = 11 mm), immediately prior to preload and test protocol. C, Multigeneration model after swelling and fiber deposition at prescribed height (H = 11.5 mm), immediately prior to preload and test protocol where k o is the isotropic hydraulic permeability in the reference state and M is the exponential strain-dependence coefficient, both of which are specified for each disc constituent based on tissue tests (Table 1). 6,13,21 The fluid content is established through Donnan equilibrium swelling pressure which established the pressure that would be produced if the matrix were to be populated with charged ions and surrounded by an external bath solution that contained counter-ions. The Donnan equilibrium response imposes the Cauchy stress for the material as: where π is the osmotic pressure, a function of the gas constant (R), temperature (T), osmotic coefficient (Φ), fixed-charge density (c F ) with reference to the initial fixed-charge density (c F 0 ), and the osmolarity of the external bath ( c * ): where the instantaneous fixed-charge density depends on the reference configuration as follows: The bath osmolarity was 300 mOsm/L and the reference state fixed charge density c F 0 is specified for each constiuent based on tissue-level experiments ( Table 1). 6

| Fibers
The AF constitutive model included fibers within the isotropic matrix to provide tensile stiffness. The nonlinear fibers were described by a strain energy density: 43,44 Ψ n I n ð Þ= where the initial stretch parameter was defined as I 0 = λ 2 0 ; similarly, the instantaneous stretch was I n = λ 2 n , and the initial strain energy den- The fiber modulus (E f ), toe-region power law exponential (β f ), and transition stretch ratio (λ o ) for each AF region, were specified based on tensile tissue testing of the AF (Table 1) 31,32,34 and the fiber angles were based on optical microscope image analysis. 5 Each AF element had two fiber sets, oriented at the positive and negative fiber angle (θ). The NP and CEP did not include fiber contributions.

| Rigid bodies and boundary conditions
The VB were assigned a Neo-Hookean constitutive relation with Young's modulus 10,000 MPa and Poisson's ratio of 0.3. 4,45,46 The superior and inferior surfaces of the model were defined by rigid bodies. The inferior rigid body was fixed in all degrees of freedom for the entirety of all testing scenarios. The superior rigid body was used to impose all loads, displacements and rotations as specified for each testing protocol ( Table 2).

| Human disc model protocol
The swelling-only model did not include pre-strained fibers, such that

| Swelling
All models were initiated with a bulged AF ( Figure 2A) and the fixed charge density in all tissues was zero. The fixed charge densities were then increased to full value (Table 1), altering the disc's stress state and resulting in a swollen disc geometry. The swelling-only model possessed both fiber sets at model initiation, time zero, while the multigeneration model did not have any fibers in the AF during initial swelling. For the swelling-only model, the disc height was fixed at H = 11 mm. For the multigeneration model, the disc height was increased while the tissue swelled to H = 11.5 mm (see Figure S2 for rationale). All models were allowed to reach the appropriate height and achieve Donnan equilibrium (Figures 2B, C and 3B).

| Fiber implementation (multigeneration model only)
The swelling-only model possessed both fiber sets at model initiation; the following fiber implementation process is relevant for only the multigeneration model. At the end of the swelling phase the disc was in a neutral position ( Figure 3B), it was then twisted −Ω degrees to a deformed state ( Figure 3C). While the disc was in the deformed state the first fiber set was positioned at +θ, according to angles specified in Table 1. The disc was then twisted back to the neutral position which stretched the first fiber set ( Figure 3D). From the neutral position the disc was twisted +Ω degrees to the opposing deformed state and the second fiber set was positioned at −θ ( Figure 3E). Then disc was returned to a neutral position, where both fiber sets were now stretched ( Figure 3F). Once both fiber sets were placed and the disc was in a neutral position, the material properties of the fibers were applied at T A B L E 2 Test case specific protocols that followed after the pre-test conditions shown in Figure 3 Test F I G U R E 3 Protocol schematic showing (A) the initial, unswollen state with superior reference marker (yellow). B, Axial displacement (ΔH = 0.5 mm) was applied and the disc was swollen to Donnan equilibrium. C, The disc was twisted −Ω from the neutral position and the first fiber set was deposited with negligible fiber modulus, then (D) it was returned to neutral position, which stretched the first fiber set. E, The disc was then twisted +Ω, and the second fiber set was deposited and then (F) it was then returned to the neutral position so that both fiber sets were stretched. G, Once both fiber sets were placed and stretched, the material properties of the fibers were increased to full value (see Table 1). This process established a disc model with residual stress contributions from both swelling and multigenerational AF fibers. H, The applied preload and test case protocols that followed are outlined in Table 2 their full values (Table 1) while boundary conditions were held fixed ( Figure 3G). This process ultimately resulted in the multigeneration model: a swollen disc model (H = 11.5 mm) with pre-strained crisscross fiber sets in the AF ( Figure 1B, Figure 3G). For comparison, the swelling-only model did not undergo this process and was, therefore, a swollen disc model (H = 11 mm) with crisscross fiber sets in the AF.

| Test cases
Seven test cases were utilized to validate the models and select the optimal multigeneration model parameters, axial displacement (ΔH), and twist angle (Ω). Our previous model 4 where Y(x) is the parameter of interest, for example in slow ramp Y (x) = Load(time). The mean square error between the model and experimental mean was normalized by the mean square error of the 95% confidence interval and experimental mean. This normalization enables evaluation of the model's ability to match experimental outcomes across different protocols. The closer the NMSE is to zero, the better the fit, and an NMSE greater than one indicates that the model outcome is outside of the 95% confidence interval for that particular test case.

| Model geometry
In addition to checking the model's test case outcomes, we also sought to validate the model's geometry to ensure it effectively mimicked physiological discs. The L4/L5 discs that underwent MRI and were used to construct the initial model geometry were also scanned following an overnight 50 N load (n = 7) ( Figure 4). 48,51 This load was imposed on the model following swelling and fiber placement (in the multigen model) such that the model's geometry could be compared to MRI data. The discs' height from MRI was calculated from a midcoronal slice, where points along the superior and inferior disc boundary were defined and the median distance between the superior and inferior markings was taken to be the disc height ( Figure 4A). The height of the disc model was calculated as the median difference between the superior and inferior disc surfaces ( Figure 4B). The leftright bulge for both the MRI discs and disc model was determined by measuring the disc's lateral extrusion beyond the endplates. The bulge was measured from an endplate reference line to the outermost disc boundary ( Figure 4A-B).

| Fiber stress and strain
The stress and strain of the AF fibers was quantified on an elementby-element basis in the fiber direction. The coordinate axes were rotated to align with the fiber direction and the instantaneous deformation gradient with respect to the deformation gradient at the time of fiber deposition was calculated. This ultimately enabled calculation of the right Cauchy-Green strain from which the strain along the fiber direction was determined. The fiber strain in conjunction with the fibers' defined material properties (Table 1)

| Cauchy stress and Lagrange strain in disc
The stress and strain state of the multigeneration model at key time points was quantified in terms of local anatomic axes. Cauchy stress and Lagrange strain were transformed from x, y, z directions to the axial, radial, and circumferential directions in accordance with our prior work. 4,39,48,51 The directional transformation enabled improved interpretation of the stresses and strains with respect to the disc structure.

| Bovine disc model protocol
In order to demonstrate our models' ability to mimic the observed bovine disc opening outcome, 18 we simulated the bovine disc as a cylinder. The constitutive equations and material properties were the same as those used for the human disc model ( Table 1) the multigeneration model opening was closer to that seen experimentally. 18 All human disc models had a run-time of less than an hour and the bovine disc model had a run-time of less than 10 minutes.

| Mechanical loading test cases for human disc model
The nonlinear load-displacement model response was compared to quasistatic and viscoelastic experimental data and to dynamic multiaxial

| Human disc model geometry
All models maintained a disc height within the experimental confidence interval after equilibration with a 50 N load ( Figure 6B). The models produced lateral bulge in the upper range of that seen in MRI; for the multigeneration model, bulge reduced as twist angle increased ( Figure 6C).

| Model state: End of 270 N preload
The fiber and disc stress and strain state was evaluated after equilibration to a 270 N load, representing the preload in the muliaxial  Figure 7C). Similarly, there was high fiber stress in the innermost AF ( Figure S4E). The multigeneration model did not have fiber strain concentrations, instead there was uniform fiber strain ( Figure 7E) throughout the AF and a uniform circumferential Cauchy stress gradient ( Figure 7G). The model also had a uniform fiber stress of~0.1 MPa ( Figure S5G).

| Model state: Maximum axial compression
The fiber and disc stress and strain state was evaluated at maximum The multigeneration model featured more uniform fiber strain ( Figure 8E) and fiber stress ( Figure S7E) profiles throughout the AF and the circumferential Cauchy stress ( Figure 8C) was highest in the OAF.

| Model state: Maximum torsion
The fiber and disc stress and strain state was evaluated at maximum torsion from the multiaxial dynamic torsion test, for the swelling-only model ( Figure

| AF fiber contributions across loading scenarios in human disc model
In order to compare fiber contribution across the disc, the fiber strain and stress in each element from the disc's coronal, mid-height was plotted across the disc for several loading conditions. Fiber stretch less than one indicates the fibers buckled and, therefore, do not contribute any fiber stress. In the swelling-only model fiber stretch generally decreases from the outermost OAF toward inner OAF ( Figure 10A) whereas the multigen model fiber stretch increases from outermost OAF toward the inner OAF ( Figure 10D). Both models tended to have peak fiber stretch at the NP transitional layer, with some cases possessing a secondary peak at the AF transition ( Figure 10A,D). The multigeneration model tended to have relatively uniform stretch in the OAF ( Figure 10D).
For the torsion test case, the models were rotated in the negative direction which results in both models exhibiting significant fiber stretch in the fiber set placed at −Ω ( Figure 10A,D) but minimal stretch and much buckling occur in the fiber set placed at +Ω (Figure 10B,E). In the swelling-only model, fiber stress was highest in the innermost IAF layer and fiber stress generally decreased from the outermost OAF inward ( Figure 10C). The multigeneration model also tended to have highest fiber stress in the innermost IAF layer, but fiber stress generally increased from the outermost OAF inward ( Figure 10F), unlike the swelling-only model ( Figure 10C). Collectively, the fiber stretch and fiber stress were found to exhibit substantially different profiles between the swellingonly and multigeneration models.

| Bovine disc opening angle
The initial opening of the multigeneration model (4.5 mm, Figure 11H) was within the experimental range (4.3 ± 1.8 mm) 18 ; the initial opening of the swelling-only model was greater than expected (11.5 mm, Figure 11D). As was seen in the human disc model, there is a high concentration of both fiber strain ( Figure 11A) and fiber stress ( Figure 11B) at the AF/NP transition and an additional fiber stress concentration at the inner/ outer AF transition. Following radial incision, the fiber strain ( Figure 11C) and stress ( Figure 11D) in the outer AF reduce dramatically; however the innermost AF layer maintains a high fiber stress concentration. The multigeneration model exhibited relatively uniform fiber strain ( Figure 11E) throughout the OAF, with a slight concentration in the innermost AF layer.
Similarly, there is a low, uniform fiber stress throughout the OAF with a slight concentration in the innermost AF layer ( Figure 11F). The multigeneration model after radial incision exhibits minimal fiber strain, predominantly at the cut plane ( Figure 11G) and fiber stress only at the AF/NP transition, but at a substantially reduced magnitude compared to the swelling-only model ( Figure 11D).

| DISCUSSION
This study successfully incorporated residual strain due to swelling and multigeneration fibers in a finite element model and validated the model against uniaxial quasi-static and multiaxial dynamic tests. We compared a swelling-only model, with AF residual strain due to swelling, and a multigeneration model, with AF residual strain due to both swelling and inherent AF fiber pre-strain. We conducted a parametric analysis for the multigeneration model with several parameter combinations to determine the optimal axial displacement (ΔH = 0.5 mm) and twist angle (Ω = 3 ) for simulating experimental outcomes and by F I G U R E 1 1 Bovine disc models before and after radial incision. Prior to incision, the swelling-only model had highly concentrated fiber strain at the AF/NP transition (A) and significant fiber stress at both the inner/outer AF and AF/NP boundaries (B). Following incision, the fiber strain reduced substantially, though a concentration at the AF/NP transition remained (C) and similarly, the fiber stress was still highly concentrated at the AF/NP boundary (D). The multigeneration model prior to incision exhibited a uniform, moderate fiber strain in the OAF and a minimal fiber strain in the IAF (E), there was a slight concentration of fiber stress at the AF/NP transition (F). After incision, the fiber strain (G) and stress (H) reduced with only minimal stress present at the AF/NP transition extension, the physiological disc. The models possessed material properties defined by tissue constituent tests and experimental work such that they were not tuned to fit any organ-scale outcomes investigated. The swelling-only model lacked the ability to tune fiber residual contributions. The multigeneration model provided a better overall response compared to experiments and had fiber residual strain contributions that could be adjusted independently of the fiber material properties. We were also able to match the residual strain opening gap from bovine experiments with the multigeneration model. 18

| Residual state
The residual state has been investigated in other fibrous soft tissues and multiple underlying mechanisms have been suggested. Fibrous soft tissue residual stress was initially investigated in artery, where it became evident that the unloaded state and stress-free state of a tissue were not the same. 54 Later work in disc showed similar outcomes, where a disc cut radially opened as initial residual strains were released. 18 The sources of residual stress exist in a hierarchy of micro-, tissue-, and organ-scale contributions. 55 The disc opening is driven by relaxation of tissue-level residual stress which is initiated during early developmental stages. 18,21 In the early development of the AF there are angled actin stress fibers present whose initial orientation has a significant impact on the subsequent alignment of the AF fibers. 19 As the disc grows, the highly organized and aligned fibers are stretched. The AF possesses radial gradients of fiber collagen content, fiber alignment, proteoglycans and water that ultimately result in fully developed discs possessing both inherent fiber pre-strain as well as swelling-induced residual strain. 5,8,18,19,21 In addition to growth, residual contributions in the disc are influenced by the loading state of the tissue, where less axial compression results in less fiber engagement ( Figure 7A,E) and greater axial compression results in more fiber engagement ( Figure 8A,E). Although the mechanisms may be similar, the use of multigeneration is intended to bring the fibers to a residual strain state similar to that observed in adult discs; our use of multigeneration was not intended to simulate the growth and development of fibers in the disc.
We found modeling the bovine disc incision including residual strain due to swelling alone caused the simulation to overestimate the discopening gap. The multigeneration model, with the inclusion of both swelling and AF fiber residual strain contributions, resulted in an opening gap that matched the experiments. 18 These outcomes support the necessity of residual strain considerations in disc finite element models and further support that two primary mechanisms of tissue-level residual strain in the disc are osmotic swelling and pre-strained AF fibers.

| Stress and strain state
We found that the NP was uniformly pressurized in all test cases and

| Parametric analysis of the multigeneration model
We conducted a parametric analysis of the multigeneration model parameters (axial displacement (ΔH) and twist angle (Ω)) because the model geometry at the time of fiber placement impacts the fiber prestrain imposed and subsequently, the model's performance in test cases. Ultimately, across all seven test cases, the multigeneration model that provided the best overall fit to experimental data was the model with ΔH = 0.5 mm and Ω = 3 . When a larger twist angle (Ω = 4 ) was used, the stress relaxation response was poor, due to the large initial load that propagated throughout the subsequent relaxation ( Figure 5C). When a smaller twist angle (Ω = 2 ) was imposed, the fit to torsion was poor because the model did not possess sufficient fiber contribution to support the disc in torsional rotation ( Figure 5E).
Therefore, we selected a twist angle of 3 as a compromise between these effects. In addition, through parametric analysis we found that an axial displacement of 0.5 mm was necessary to recapitulate experimental outcomes, particularly for quasi-static uniaxial slow ramp and creep tests (Figures S2 and S3).

| Material property selection
An advantage of our disc model is that the material properties are defined based on tissue testing and were not tuned to fit the organ-scale outcomes. Moreover, we included radial variation between the inner and outer AF and included transitional regions between the AF sections and the NP. The matrix properties and tissue permeability were determined from confined compression tissue testing (Table 1). 6,13 The fixed charge density of all tissues, essential for simulating appropriate swelling, was determined from proteoglycan assessment in each tissue (Table 1). 13,21 The fiber properties were based on tensile tissue testing 31,32,34 and fiber angles were from optical microscope assessment (Table 1). 5 Our material property definitions are rigorously rooted in experimental tissue tests, such that each tissue represented in our model accurately mimics the mechanical contribution of the physiological tissues and collectively enables our organ-scale model to optimally represent the whole, physiological disc.
We included radial variation in the material properties; however, we do not include regional circumferential differences between the anterior, posterior, and lateral regions, such as higher proteoglycan content in the posterior region. 8 There is also contradictory information in the literature as to whether there are regional differences in fiber angle. 5,33 Optical microscope analysis of AF tissue found no circumferential differences, only substantial radial variation in fiber angle 5 ; however, later experimental work utilized surface marking and digital imaging which revealed regional variation in fiber angle. 56 In addition, mathematical modeling suggested that the inner posterior region of the disc has substantially higher fiber angle than through the rest of the AF. 33 While our model does not possess regional variation, we account for radial differences in material properties which are greater in magnitude.

| Limitations
The organ-scale size and shape of the disc under load was generally reasonable with the exception of excess posterior protrusion (see Figure 8).
A similar model previously used by our lab found that the model showed the greatest difference from the MRI measurements in the posterior region where the model had excess radial strain compared to the MRIbased strain analysis. 48 The underlying cause was suggested to be elongated rectangular elements in posterior unintentionally causing small discontinuities in the fiber distribution. 48 The posterior protrusion might be corrected in future work by refining the mesh in that region, though in the present cases it does not impose unreasonable stress or strain concentrations nor does it impair the model's gross response.
We were unable to explicitly define fiber pre-strain in our biphasic model in FEBio, which limits the ability to control the pre-strain imposed on the AF fibers. We considered multiple placement methods to generate fiber pre-strain and ultimately found that the multigeneration feature with organ-scale axial displacement and rotation created an adequate deformation state for fiber placement to produce pre-strain. This technique to impose fiber pre-strain was not intended to represent physiological loading for either the human or bovine discs, it was purely to establish prestrained fibers in the AF. Fiber placement by this method was limited as all fibers were placed with the disc in the same organ-scale deformation state; however, explicit ability to impose particular pre-strain would greatly improve the control and precision of fiber pre-strain. FEBio was previously used to implement explicitly defined fiber pre-strain in a ligament model, 57 but we were unable to use their methods here due to the complexity of our biphasic model. Expansion of FEBio's pre-strain fiber method in future versions would allow for pre-strained fiber contributions to be more readily utilized in fibrous soft tissue models.

| CONCLUSION/RECOMMENDATION
In this study, we successfully incorporated residual strain due to swell- 6 | SUPPLEMENTAL

| Previous model
Our lab previously developed a disc finite element model, detailed in Reference 4. The model was validated against uniaxial quasi-static slow ramp, creep, and stress relaxation, but when expanded to multiaxial dynamic tests it was unable to recapitulate experimental outcomes ( Figure S1).

| Parametric analysis of multigeneration model: Axial displacement
The normalized mean square error (NMSE) was calculated as explained in the main manuscript Methods 2.4.1) Nonlinear Disc Response. Variations in axial displacement majorly impacted the uniaxial, quasi-static outcomes in creep and stress relaxation (ΔNMSE > 0.35) and mildly impacted slow ramp and dynamic torsion (0.25 < ΔNMES < 0.15) ( Figure S3). Dynamic multiaxial tests axial compression, bending, and flexion were minimally impacted by the axial displacement parameter (ΔNMSE < 0.1) ( Figure S3).
An axial displacement of 0.5 mm to achieve a fiber placement height of 11.5 mm was deemed necessary based largely on the slow ramp and creep outcomes. A fiber placement height of 11 mm over-estimated compressive displacement and a fiber placement height of 12 mm underestimated compressive displacement in both tests ( Figures S2 and S3). Furthermore, fiber placement height of 11.5 mm was necessary to compromise between the stress relaxation and torsion outcomes, where stress relaxation response improved with increasing fiber placement height but the torsion response improved with decreasing fiber placement height ( Figures S2 and S3).

| Human disc fiber and model stress and strain states
In addition to the fiber strain and Cauchy stresses in the main paper figures, the fiber stress as well as axial, circumferential and radial Lagrange strains were also evaluated for the swelling-only model at the end of preload ( Figure S4), maximum axial compression ( Figure S6), and maximum torsion ( Figure S8) and for the multigeneration model at the end of preload ( Figure S5), maximum axial compression ( Figure S7), and maximum torsion ( Figure S9).