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

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. , , , Our lab previously developed a disc model with such properties and validated it against uniaxial, quasi‐static compression tests. 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 ; 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. , , , , 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. , , , , 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). , , , , 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. , Bovine discs removed from the endplates and then cut radially, opened up as the internal residual strain is released. 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. , , 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). , While many models have included residual strain contributions due to swelling, , , , , , , , , , , 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).

FIGURE 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 pre‐strain. 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)

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.

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

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

TABLE 1

Material properties for outer annulus fibrosus (OAF), inner annulus fibrosus (IAF), nucleus pulposus (NP), and cartilage endplates (CEP)

	Matrix properties (Cortes+ 2012, Cortes+ 2014, DeLucca+ 2016)						Fiber properties (Elliott+ 2001, O'Connell+ 2009, Jacobs+ 2013, Cassidy+ 1989)
	Modulus E m (MPa)	Poisson ratio ν	Exponential stiffening coefficient β m	Hydraulic permeability k 0 (mm4/Ns)	Exponential strain dependence M	Fixed charge density (mM)	Modulus E f (MPa)	Transition stretch ratio λ o	Fiber angle θ (± degrees)	Toe region power law exponent β f
OAF	0.018	0.24	3.4	0.0047	5.75	44	15.6	1.028	31.0	4
AFtrans	0.023	0.20	2.8	0.0036	4.60	50	10.3	1.025	38.5	4
IAF	0.026	0.16	2.1	0.0025	3.50	55	6.9	1.023	41.5	4
NPtrans	0.045	0.20	1.5	0.0016	2.71	217	3.0	1.020	44.5	4
NP	0.065	0.24	0.95	0.00056	3.79	379	N/A	N/A	N/A	N/A
CEP	0.305	0.18	0.29	0.00056	3.79	248	N/A	N/A	N/A	N/A

Note: AFtrans is a 1‐element wide transitional layer between the OAF and IAF; similarly, NPtrans is a 1‐element wide transitional layer between the IAF and the NP. The material properties are from tissue‐level experimental work and were either calculated or fit via inverse finite element modeling.

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 degeneration grade 3. , 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).

FIGURE 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

Constitutive models and material properties

Matrix

The AF, NP, and CEP all included extrafibrillar matrix with a Holmes‐Mow constitutive equation. Holmes‐Mow matrix is an isotropic, hyperelastic material model which is widely used to model disc. , , , , 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, β is the matrix stiffening coefficient, and λ and μ are the Lamé parameters which can be related to the Young's modulus E and Poisson's ratio ν as follows:

E , ν, and β were specified for each constituent based on tissue‐level compression tests (Table 1). ,

Permeability and Donnan swelling

Holmes‐Mow permeability was utilized to achieve isotropic, strain‐dependent permeability according to the following permeability tensor: , 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). , ,

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 ) with reference to the initial fixed‐charge density (), and the osmolarity of the external bath (): 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 is specified for each constiuent based on tissue‐level experiments (Table 1).

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: , where the initial stretch parameter was defined as ; similarly, the instantaneous stretch was , and the initial strain energy density was .

The fiber modulus (E ), toe‐region power law exponential (β ), and transition stretch ratio (λ ) for each AF region, were specified based on tensile tissue testing of the AF (Table 1) , , and the fiber angles were based on optical microscope image analysis. 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. , , 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).

TABLE 2

Test case specific protocols that followed after the pre‐test conditions shown in Figure 3

Test case	Test type	Preload	Test control	Test condition	Experimental data source
Slow ramp	Quasi‐Static Uniaxial	Load at ΔH	Axial load (z)	1 N/s to 2000 N	O'Connell+ 2011
Creep	Quasi‐static uniaxial	Load at ΔH	Axial load (z)	200 N/s to 1000 N, then hold	O'Connell+ 2011
Stress relaxation	Quasi‐static uniaxial	50 N preload	Axial displacement (z)	5% compression in 5 seconds, then hold	Yoder+ 2014 Showalter+ 2016
Axial compression	Dynamic multiaxial	270 N preload	Axial load (z)	18 N/s to 900 N	DeLucca+ 2019
Torsion	Dynamic multiaxial	270 N preload	Axial Rotation (Rz)	0.08°/s to 3°	DeLucca+ 2019
Bending	Dynamic multiaxial	270 N preload	Sagittal Rotation (Ry)	0.08°/s to 3°	DeLucca+ 2019
Flexion	Dynamic multiaxial	270 N preload	Coronal Rotation (Rx)	0.08°/s to 3°	DeLucca+ 2019

Note: Preload was enacted in axial force control with all other degrees of freedom fixed, the preload was held for 12 hours in accordance with the experimental protocols. For the test sequence, the test control was the only degree of freedom specified, all other degrees of freedom were fixed.

FIGURE 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

Human disc model protocol

The swelling‐only model did not include pre‐strained fibers, such that all residual strain was exclusively swelling‐induced. The swelling‐only protocol involved three stages: swelling at fixed height (H = 11 mm) (Methods 2.3.1, Figure 2B), axial preload (Methods 2.3.3, Table 2, Figure 7A‐D), and loading simulation for test case of interest (Methods 2.3.3, Table 2, Figures 8A‐D, 9A‐E).

The multigeneration model included residual strain contributions from both swelling and AF fiber pre‐strain. The fiber pre‐strain was dependent upon the disc geometry at the time of fiber placement. Two parameters: axial displacement (∆H) and torsional twist angle (Ω) were used to evaluate the impact of varying fiber pre‐strain on organ‐scale outcomes. We found a disc height of 11.5 mm (∆H=0.5 mm) was necessary to place fibers that produced reasonable uniaxial outcomes in slow ramp and creep (Figures S2 and S3) and this was therefore used in all subsequent multigeneration models herein. The twist angle was varied to Ω = 2° , 3° , 4° to investigate how torsional placement of fibers impacted model outcomes. To this end, the multigenerational protocol involved four stages: swelling with axial displacement (∆H = 0.5) (Methods 2.3.1, Figures 2C and 3B), multigenerational fiber deposition (Methods 2.3.2, Figure 3C‐G), axial preload (Methods 2.3.3, Table 2, Figure 7E‐H), and loading simulation for test case of interest (Methods 2.3.3, Table 2, Figures 8E‐H and 9F‐J).

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.

The multigeneration approach was used to add pre‐strained AF fibers to establish fiber residual strain independent of swelling residual strain. The extra‐fibrillar matrix was established at time zero in a non‐stressed state while the fibers were deposited in a stressed reference configuration that differed from the neutral state of the model. This provided for the disc model to possess stretched fibers creating fiber‐induced, circumferential residual stress in the AF. The stressed reference configuration was achieved by twisting the disc model about the spinal axis. The twist (Ω) was varied between model runs to 2°, 3°, and 4°.

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 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 was validated against uniaxial quasi‐static tests; slow ramp and creep tests mimicked the experimental protocol of Reference 47 and the stress relaxation test mimicked the protocol of Reference 48. These tests were repeated on the current models to ensure agreement and no loss in the current model's ability to simulate these outcomes. , , Furthermore, the present models were tested in multiaxial dynamic tests which mimicked the experimental protocol in Reference 49 for axial compression, torsion, bending, and flexion tests.

For the dynamic tests, the preload was calculated by the same process utilized in the experiment, where the cross‐sectional area of the disc was estimated by A = 0.84 * L right/left * L anterior/posterior and the preload necessary to achieve physiologically relevant nucleus pressurization was calculated as . , The disc length measurements (L) were taken from the images of the discs used to develop the model geometry (n = 7). The preload calculated for the dynamic tests was 270 N, this load was imposed and allowed to equilibrate in the simulation over a 12‐hour period, as was done in the experimental tests, prior to the test case. The test case loading/rotation imposed were those from the final experimental cycle of the slowest test frequency for each test condition (Table 2).

Human disc model outcomes

Nonlinear disc response

The average disc response and associated 95% confidence interval were calculated for all experimental test cases: slow ramp (n = 4), creep (n = 4), stress relaxation (n = 5), axial compression (n = 8), torsion (n = 8), bending (n = 8), and flexion (n = 8). , , , To quantify the model's ability to replicate experimental outcomes, a normalized mean square error (NMSE) was calculated as follows: 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). , 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 mid‐coronal 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 left‐right 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).

FIGURE 4

A, MRI of L4‐L5 disc following overnight 50 N load (Yoder 2014+, Showalter 2016+). Cyan dots on the superior/inferior edges were used to calculate the median disc height from a mid‐coronal MRI slice. Red reference lines connect the outermost points of the vertebral bodies. The lateral bulge was measured by the yellow line, from the end of the vertebral bodies reference (red line) to the outermost point of the disc. B, Mid‐coronal image of the multigeneration model after equilibration to 50 N load. The dashed cyan lines represent the superior/inferior disc surfaces which were used to calculate the median disc height. Disc bulge was calculated in the same manner as used for the MRI discs

Fiber stress and strain

The stress and strain of the AF fibers was quantified on an element‐by‐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) allowed for the calculation of fiber stress. For further details about the fiber stress and strain calculations please see supplemental section 6.5 Calculation of the Fiber Strain and Stress.

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. , , , 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, 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 model geometry was chosen such that the end‐of‐swelling diameter matched the experimentally measured mean bovine disc size. The swelling‐only model was initiated with fibers present in the AF, then allowed to swell to osmotic equilibrium as explained in Methods 2.3.1. The multigeneration model was swollen, then subjected to multigeneration fiber placement with a twist angle of 3°, according to the process previously explained in Methods 2.3.2. For both models, following swelling and fiber placement (multigeneration model only), the radial cut was simulated by release of planar boundary conditions along the radial interior surfaces on either side of the cut. The opening gap after the radial cut (Wi) was quantified and compared with the experimental data. Additional experimental work was conducted with partial and full nucleus detachment, but we were unable to simulate these experiment with our current model as they were conducted in open air while our model assumes immersion in a bath solution. ,

RESULTS

In general, for the human disc, both the swelling‐only and multigeneration models performed well in comparison to multiaxial experimental data. The multigeneration model with twist angle Ω = 3° produced the best set of outcomes across all seven test cases. For the bovine disc, both the swelling‐only and multigeneration models had a positive opening angle; the multigeneration model opening was closer to that seen experimentally. 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 quasi‐static and viscoelastic experimental data and to dynamic multiaxial experimental data and the error was calculated. In a quasi‐static axial slow ramp, the multigeneration models performed better (NMSE < 0.05) than the swelling‐only model (NMES = 0.25) (Figures 5A and 6A). Similarly, in quasi‐static creep, the multigeneration models (NMSE < 0.15) performed better than the swelling‐only model (NMSE = 0.54) (Figures 5B and 6A). In contrast, in stress relaxation, the swelling‐only model (NMSE = 0.06) performed better than the multigeneration models (NMSE > 0.11) (Figures 5C and 6A). The multiaxial dynamic response of all models matched the experimental response in axial compression (NMSE < 0.15), bending (NMSE = 0.01), and flexion (NMES < 0.1) (Figures 5D‐G and 6A). The dynamic torsion response varied across multigeneration models (0.15 < NMSE < 0.44) and the swelling‐only model performed adequately (NMSE = 0.39) (Figures 5E and 6A).

FIGURE 5

Experimental data ±95% confidence interval is shown with model outcomes for swelling‐only and multigeneration models with ΔH = 0.5 mm and variable twist angle (Ω). All models lie within the 95% confidence interval with for all tests. The multigeneration models perform better than swelling‐only in slow ramp (A), creep (B), and torsion (E). The swelling‐only model is better in stress relaxation (C) and axial compression (D). Increasing the twist angle in the multigeneration model (from Ω = 2° to 4°) improves the match to experimental data for torsion (E), but worsens the match for stress relaxation (C). All models sufficiently recapitulate bending (F) and flexion (G) outcomes

FIGURE 6

A, Normalized mean square error (NMSE) values are shown for the swelling‐only model and multigeneration model with varying twist angle (Ω). NMSE values closer to zero indicate a model response closer to the experimental mean. B and C, Disc height and bulge from the models equilibrated to a 50 N load compared to MRI data of discs under a 50 N load. B, All models were within expected height limits. C, All models are in the upper range of disc bulge

Three versions of the multigeneration model were tested with fiber twist angles of Ω = 2°, 3°, and 4°, to evaluate the impact of the multigeneration model input of fiber twist angle on the model's mechanical outcomes. The variation in twist angle had minimal impact on slow ramp, creep, axial compression, bending, and flexion responses, where the maximum difference across twist angles was only ΔNMSE = 0.05. However, variation in twist angle had a significant impact on the stress relaxation response (ΔNMSE = 0.21) and torsion response (ΔNMSE = 0.29). The stress relaxation response improved with decreasing the twist angle; conversely, the torsion response improved with increasing twist angle. The multigeneration model with twist angle Ω = 3° was selected as the optimal twist angle, this decision was primarily driven by the necessary compromise between the stress relaxation and torsion responses.

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 dynamic tests, for the swelling‐only model (Figure 7A‐D) and multigeneration model with Ω = 3° (Figure 7E‐H). For both models, the fiber strain profiles (Figure 7A,E) were similar between both fiber sets, therefore only one set is shown. The swelling‐only model had a concentration of fiber strain (Figure 7A) in the NP/AF transition and innermost AF which caused high circumferential Cauchy stress in the innermost AF layer (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).

FIGURE 7

Models after equilibration to 270 N (the preload for all dynamic multiaxial tests). The fiber strain profiles (A, E) were similar for both fiber sets, only one set shown. For the swelling‐only model there was a concentration of fiber strain (A) in the inner AF which caused axial (B) and circumferential (C) Cauchy stress in the same area. For the multigeneration model (Ω = 3°) there was uniform fiber strain (A) throughout the AF and a slight concentration of circumferential stress (G) in the disc posterior

Model state: Maximum axial compression

The fiber and disc stress and strain state was evaluated at maximum axial compression from the multiaxial dynamic compression test, for the swelling‐only model (Figure 8A‐D) and multigeneration model with Ω = 3° (Figure 8E‐H). The swelling‐only and multigeneration model profiles are similar to those seen at the end of the 270 N preload, but with greater magnitude as expected due to the increased axial load. The swelling‐only model had a high concentration of fiber strain (Figure 8A) in the inner most AF and moderate fiber strain through the remaining AF. The fiber stress (Figure S6E) followed a similar concentration pattern, with additional concentrations in the middle AF region, particularly at the superior and inferior surfaces. 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.

FIGURE 8

Models at maximum axial compression (from multiaxial dynamic axial compression test case); the fiber strain profiles (A, E) were similar for both fiber sets, only one set shown. For the swelling‐only model, there was a high concentration of fiber strain (A) in the inner most AF and moderate fiber strain through the middle and outer AF. The inner AF was mildly compressed in all directions (B‐D) and the fiber contribution is evident in the circumferential Cauchy stress (C) in the outer AF. The multigeneration model (Ω = 3°) had moderate fiber strain (E) throughout the AF except in the inner, posterior region. There was a moderate, uniform circumferential Cauchy stress (G) in the outer AF

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 9A‐E) and multigeneration model with Ω = 3° (Figure 9F‐J). The swelling‐only model at maximum axial torsion experienced highest fiber strain at the NP/AF transition and uniformly moderate fiber strain throughout the remaining AF for the fibers aligned in the direction of torsional rotation (Figure 9A). The fiber set aligned in the opposing direction (Figure 9B) experienced moderate fiber strain at the NP/AF transition; however, fibers in the remaining AF buckled under the torsional rotation and subsequently did not contribute fiber stress. Collectively, the total fiber stress (Figure S9F) was highest in the IAF, and there was a low circumferential Cauchy stress in IAF (Figure 9D). The multigeneration model (Figure 9F,G) once again exhibited uniform fiber contributions throughout the AF, compared to the swelling‐only model (Figure 9A,B). The fiber set aligned in the torsional direction (Figure 9F) experienced moderate fiber strain throughout the AF whereas the fiber set aligned in the opposing direction (Figure 9G) experienced minimal fiber strain in the anterior and posterior regions and the remaining AF fibers buckled. The total fiber stress (Figure S9F) was relatively uniform throughout the AF, with a slight concentration in the middle‐AF posterior region and a similar distribution was also seen in the circumferential Cauchy stress (Figure 9I).

FIGURE 9

Models at maximum torsional rotation (from multiaxial dynamic torsion test case). For the swelling‐only model, the fibers aligned in the torsion direction (A) experienced moderate fiber strain throughout the AF, and both fiber sets had a high concentration of strain at the NP/AF transition (A, B). There was slight circumferential Cauchy stress (D) in the outer AF. The multigeneration model (Ω = 3°) had moderate fiber strain throughout the AF for fibers aligned in the torsion direction (F) and minimal strain in fibers aligned in the opposing direction (G). The circumferential Cauchy stress (I) was mild in the outer AF and reduced outward

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 swelling‐only and multigeneration models.

FIGURE 10

The contribution of AF fibers in the swelling only model (A‐C) and multigeneration model (Ω = 3°) (D‐E) were evaluated across several loading scenarios. A single line of elements at the intersection of the mid‐coronal and mid‐transverse planes were chosen as a representative set. Fiber stretch less than 1 (gray area) indicates the fibers buckled and did not contribute any fiber stress

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

FIGURE 11

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

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

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. Later work in disc showed similar outcomes, where a disc cut radially opened as initial residual strains were released. The sources of residual stress exist in a hierarchy of micro‐, tissue‐, and organ‐scale contributions. The disc opening is driven by relaxation of tissue‐level residual stress which is initiated during early developmental stages. , 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. 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. , , , , 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 disc‐opening 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. 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 consistently exhibited greater Cauchy stress in the swelling‐only model compared to the multigeneration model (Figures 7, 8, 9). When subjected to the 270 N preload, the swelling‐only model exhibited a high circumferential stress in the inner AF, this hoop stress served to contain the NP pressurization (Figure 7A,C). Under the same conditions, the multigeneration model exhibited a uniform fiber contribution throughout the AF, the lack of inner AF concentration was because the fibers are deposited after the swelling simulation and their strain was therefore due largely to the multigeneration placement process, and not by the expansion of the NP during swelling (Figure 7E,G). The low magnitude of fiber contribution is expected under the 270 N load as this preload was chosen to establish a physiologically relevant intradiscal pressure, and under regular physiological conditions we do not expect high fiber engagement. The contribution of fibers was dependent on the discs' axial displacement and twist angle parameters, such that altering either would change the preloaded discs' stress and strain state.

The swelling‐only model consistently possessed high fiber contributions at the AF/NP boundary and had mixed positive and negative circumferential Cauchy stress in the remaining IAF. At maximum axial compression the swelling‐only model had concentrated fiber strain at the AF/NP boundary, though the IAF had both compressive and tensile circumferential stresses, while the OAF exhibited the greatest magnitude of tensile circumferential stress (Figure 8A,C). Under maximum torsional rotation the swelling‐only model exhibited moderate circumferential tension at the AF/NP transition and throughout the OAF. Collectively, the swelling‐only model appears to have an artificial concentration of stress and strain at the AF/NP transition which appears to shield the remaining AF from reasonably contributing to the model's overall mechanical response.

In the multigeneration model, the IAF was consistently in compression, while the outer AF was in moderate tension throughout with reducing magnitude toward the disc's outermost layers (Figure 8E,G). At maximum torsional rotation, the multigeneration model exhibited a moderate, uniform tensile stress distributed throughout the AF which reduced toward the outermost layers, as expected (Figure 9F,I). The multigeneration model does not exhibit inner concentrations and instead repeatedly features a smooth distribution of circumferential fiber contributions throughout the loading scenarios investigated. The lack of inner concentrations is consistent the prior work in artery, where inclusion of residual stress significantly reduced inner stress concentrations of the vessel following loading. The uniform fiber contribution enabled by the multigeneration model with an adequate twist angle (Ω ≥ 3°) seems necessary for recapitulating experimental outcomes in torsion; however, the increase in fiber contribution throughout the AF simultaneously worsens the model's ability to match stress relaxation outcomes.

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 pre‐strain 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). , The fixed charge density of all tissues, essential for simulating appropriate swelling, was determined from proteoglycan assessment in each tissue (Table 1). , The fiber properties were based on tensile tissue testing , , and fiber angles were from optical microscope assessment (Table 1). 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. There is also contradictory information in the literature as to whether there are regional differences in fiber angle. , Optical microscope analysis of AF tissue found no circumferential differences, only substantial radial variation in fiber angle ; however, later experimental work utilized surface marking and digital imaging which revealed regional variation in fiber angle. 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. 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 MRI‐based strain analysis. The underlying cause was suggested to be elongated rectangular elements in posterior unintentionally causing small discontinuities in the fiber distribution. 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 pre‐strained 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, 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 swelling and multigeneration fibers in a finite element model of intervertebral disc and validated human disc models against uniaxial quasi‐static and multiaxial dynamic tests. The swelling‐only model was within the confidence interval for all human disc outcomes, though it overestimated the opening of the incised bovine model. The multigeneration model provided better human disc model responses, closer to the experimental mean and had a bovine model opening angle within the experimental range. The use of multigeneration allowed for the inclusion of inherently pre‐strained fibers in AF, which produced uniform fiber contribution throughout the AF and parametric analysis found a twist angle of 3° and axial displacement of 0.5 mm were necessary for optimizing the multigeneration model outcomes. The inclusion of swelling and fiber‐induced residual strain in the multigeneration model was necessary for achieving a physiological residual strain state and for replicating disc mechanical behavior across uniaxial quasi‐static and multiaxial dynamic test cases.

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

Bovine disc model stress and strain states

The axial, circumferential and radial stress and strain were quantified for the bovine disc models preceding radial incision and immediately following the incision (Figure S10).

Calculation of the fiber strain and stress

AF fiber strain and stress were quantified to evaluate their contribution to the gross mechanical behavior of the disc. This was accomplished on an element‐by‐element basis within the layers which had residual stress from AF fibers. Fibers were aligned by specification of local [1, 2, 5] axis (Figure S11, see FEBio manual for further details), each element had node positions (x, y, z) which were used to determine the and directions at the start of the simulation (Figure S11):

These initial directions were used to calculate a set of orthonormal fiber vectors: Which were then compiled into a tensor, one for each element: The directional vector for the fibers in terms of global coordinates was then determined for time zero:

In order to calculate fiber stretch, it was necessary that the reference time be the multigeneration time, not time zero. To get the global position vector at the fiber deposition time: where F(x, u) was the deformation gradient at the time of deposition (u) with respect to time zero. With the global directional vector at deposition time known, it was necessary to determine the deformation gradient at the time of interest (t) with respect to the deposition time: Now with the appropriate deformation gradient, the Right Cauchy‐Green Strain was calculated: and the fiber stretch was determined: The fiber strain was defined as fiber stretch minus one such that: In addition to fiber stretch, the fiber stress was also quantified: where λ was the fiber stretch previously calculated, and all other properties including the fiber modulus (E ), toe‐region power law exponential (β ), and transition stretch (λ ) were specified for each constituent (Table 1).

CONFLICT OF INTEREST

The authors have no conflicts of interest to report.

AUTHOR CONTRIBUTIONS

All contributed to project concept. HRN, JFD, JMP performed modeling work. HRN, JMP, DME prepared manuscript. All contributed to manuscript revisions.

Figure S1 The Jacobs+2014 model was not able to recapitulate experimental outcomes in dynamic axial compression (A) and torsion (B). C, NMSE values closer to zero indicate a model response closer to the experimental mean and NMSE value greater than 1 indicate the model response was outside the experimental 95% confidence interval. D, The model equilibrated to a 50 N load had reasonable height compared to MRI data and had mild excess bulge (E).

Click here for additional data file.

Figure S2 Experimental data ±95% confidence interval is shown with model outcomes for the multigeneration models with Ω = 3° and variable axial displacement (ΔH). The model with ΔH = 0.5 mm was optimal for uniaxial quasi‐static slow ramp (A) and creep (B) outcomes. For uniaxial quasi‐static stress relaxation (C), the model with ΔH = 1 mm was best. There was minimal variation between the multigeneration model outcomes in multiaxial dynamic axial compression (D), bending (F), and flexion (G). The model with ΔH = 0 mm offered the best torsion response (E).

Click here for additional data file.

Figure S3 A, Normalized mean square error (NMSE) values are shown for the multigeneration model with varying axial displacement (ΔH) with twist angle Ω = 3°. NMSE values closer to zero indicate a model response closer to the experimental mean. Axial displacement ΔH = 0.5 mm was necessary to optimize slow ramp and creep outcomes. Stress relaxation was best at ΔH = 1.0 mm while torsion was best at ΔH = 0.0 mm. All model responses were reasonable in axial compression, bending, and torsion tests. B and C, Disc height and bulge from the models equilibrated to a 50 N load compared to MRI data of discs under a 50 N load. B, All models were within expected height limits. C, All models are in the upper range of disc bulge.

Click here for additional data file.

Figure S4 The stress and strain state of the swelling‐only model following equilibration to a 270 N preload, a relevant pre‐test state for the multiaxial dynamic tests. The fiber strain profile was similar for both fiber sets, only one set shown. Fiber strain and Cauchy stresses (A, F‐H) are identical to those shown in Figure 7A‐D, respectively.

Click here for additional data file.

Figure S5 The stress and strain state of the multigeneration model (Ω = 3°) following equilibration to a 270 N preload, a relevant pre‐test state for the multiaxial dynamic tests. The fiber strain profile was similar for both fiber sets, only one set shown. Fiber strain and Cauchy stresses (A, F‐H) are identical to those shown in Figure 7E‐H, respectively.

Click here for additional data file.

Figure S6 The stress and strain state of the swelling‐only model at maximum compression from the multiaxial dynamic axial compression test. The fiber strain profile was similar for both fiber sets, only one set shown. Fiber strain and Cauchy stresses (A, F‐H) are identical to those shown in Figure 8A‐D, respectively.

Click here for additional data file.

Figure S7 The stress and strain state of the multigeneration model (Ω = 3°) at maximum compression from the multiaxial dynamic axial compression test. The fiber strain profile was similar for both fiber sets, only one set shown. Fiber strain and Cauchy stresses (A, F‐H) are identical to those shown in Figure 8E‐H, respectively.

Click here for additional data file.

Figure S8 The stress and strain state of the swelling‐only model at maximum torsion from the multiaxial dynamic torsion test. Fiber strains and Cauchy stresses (A‐B, F‐H) are identical to those shown in Figure 9A‐E, respectively.

Click here for additional data file.

Figure S9 The stress and strain state of the multigeneration model (Ω = 3°) at maximum torsion from the multiaxial dynamic torsion test. Fiber strains and Cauchy stresses (A‐B, F‐H) are identical to those shown in Figure 9F‐J, respectively.

Click here for additional data file.

Figure S10 The bovine disc model immediately before (A‐C, G‐I) and following radial incision (D‐F, J‐L) for the swelling‐only model (A‐F) and multigeneration model (G‐L).

Click here for additional data file.

Figure S11 Disc model with anterior (A), posterior (P), superior (S), left (L) and right (R) directions labeled. A schematic of a single element with local node numbers and directional vectors 'a' and 'd' shown for a local material axis definition of [1, 2, 5] for more details see FEBio user manual.

Click here for additional data file.