Abstract

Intervertebral disc (IVD) degeneration and methods for repair and regeneration have commonly been studied in organ cultures with animal IVDs under compressive loading. With the recent establishment of a novel multi-axial organ culture system, accurate predictions of the global and local mechanical response of the IVD are needed for control system development and to aid in experiment planning. This study aimed to establish a finite element model of bovine IVD capable of predicting IVD behavior at physiological and detrimental load levels. A finite element model was created based on the dimensions and shape of a typical bovine IVD used in the organ culture. The nucleus pulposus (NP) was modeled as a neo-Hookean poroelastic material and the annulus fibrosus (AF) as a fiber-reinforced poroviscoelastic material. The AF consisted of 10 lamella layers and the material properties were distributed in the radial direction. The model outcome was compared to a bovine IVD in a compressive stress-relaxation experiment. A parametric study was conducted to investigate the effect of different material parameters on the overall IVD response. The model was able to capture the equilibrium response and the relaxation response at physiological and higher strain levels. Permeability and elastic stiffness of the AF fiber network affected the overall response most prominently. The established model can be used to evaluate the response of the bovine IVD at strain levels typical for organ culture experiments, to define relevant boundaries for such studies, and to aid in the development and use of new multi-axial organ culture systems.

Introduction

The intervertebral disc (IVD) is a soft connective tissue between vertebral bodies that enables flexible motion of the spine. The IVD consists of a proteoglycan-rich nucleus pulposus (NP) in the center and a collagen-type I-rich annulus fibrosus (AF) surrounding the NP. Cartilaginous endplates connect the IVD with the adjacent vertebrae. As a result of aging, trauma, or lesion, the IVD may start to degenerate, which manifests as proteoglycan, water and height loss, fiber dispersion, delamination, and fissures. Eventually, the degeneration may lead to NP herniating out of its cavity, which is commonly associated with severe back pain and often accompanying radiating pain to the limbs. IVD degeneration is a very prevalent condition, for example, in the population over the age of 50 years, over 90% have signs of degeneration at some discs of the spine [1].

Intervertebral disc degeneration and treatments to halt or reverse the degeneration are commonly studied in long-term whole-organ cultures under compressive loading [210]. Since whole human IVDs are scarcely available, bioreactor studies are using animal discs, such as caprine [8], ovine [2], rabbit [11,12], mouse [13], and most commonly bovine [47,10,1419]. Loading protocols for physiological and degenerative compressive loading have been established, such as cyclic loading between 0.02 and 0.2 MPa at 0.2 Hz for 2 h/day (physiological, able to maintain cell viability) and 0.32–0.5 MPa at 5 Hz for 2 h/day (degenerative) [5,6]. However, natural loading in the spine is not limited to compression but occurs in all the six degrees-of-freedom (translations in x, y, and z, rotations about x, y, and z axes), which likely has a large impact on the mode and extent of degeneration and the treatment success. Thus, bioreactors incorporating compression and torsion [19], or bending [15], have been developed, and our group is developing a new bioreactor system for long-term IVD organ culture, which can apply loads in all degrees-of-freedom [20]. At this stage, the system is intended to be used with bovine tail IVDs.

In the development of the bioreactor, predictions of the global and local responses of the tissue are needed to evaluate the relationship between applied global forces and moments and the resulting internal tissue strains and stresses. Global responses are needed in control system development (for force-controlled tests), whereas internal tissue strains and stresses can be used to evaluate physiological and degenerative levels of loading when planning experiments with complex multi-axial loads. Such predictions and evaluations can be efficiently made using finite element (FE) models at the organ and tissue scale (mm-scale). The load levels of interest for these evaluations are the physiological load (e.g., disc height change of ∼10% in compression [5,20]) and higher loads, which are likely to be degenerative. The loads applied in the bioreactor are dynamic, and therefore it is important to take the IVD viscoelasticity into account for these predictions. The viscoelasticity of the IVD arises from various factors, most importantly the fluid flow within and out of the tissue when it is under compression.

Previous computational models of the IVD at the organ and tissue level range in complexity from linear elastic approaches to complex nonlinear composition-based models. For example, the NP has been modeled as an incompressible fluid-filled cavity [21], linear elastic material [2124], hyperelastic material [2527], as well as fiber-reinforced poroviscoelastic with swelling [28,29]. The AF material models have, likewise, ranged from fiber-reinforced hyperelastic [2226,30] and poroelastic [31] to fiber-reinforced poroviscoelastic with swelling [2729,32]. The material property distribution within AF is commonly taken into account [2226], with outer layers having a stiffer material response than the inner layers, likely arising from the collagen type I distribution. The FE-model complexity should be adequate for the purpose, giving accurate enough results while avoiding unnecessary computational cost, which comes with increasing complexity and detail. For instance, in our preliminary analyses, a previous fiber-reinforced hyperelastic model of IVD [33], where the AF fiber network was viscoelastic, was unable to capture the time-dependent response of bovine IVD at physiological or higher strain levels in compression [34]. In this case, the viscoelasticity of the fiber network alone was not sufficient to model the viscoelastic behavior, thus porosity should be included to take into account the fluid flow, though it increases the computational cost.

The purpose of this study was to establish a finite element model of a representative healthy ex-vivo bovine tail IVD and validate it with experimentally measured compressive data. We created a model with a cylindrical geometry, taking vertebra shapes into account, and applied a novel fiber-reinforced poroviscoelastic material model, presented for the first time in this study. An extensive parametric study was conducted to investigate the effect of variation in material properties, as well as different modeling assumptions. The model created in this study is useful when evaluating the response of the bovine IVD to varying and complex loads.

Methods

Sample Preparation.

Seven bone-IVD-bone segments were cut from three fresh bovine tails obtained from a local slaughterhouse (Fig. 1(f)), using a previously established protocol [9,35]. One representative segment was later selected for finite element model creation. Muscles and ligaments were removed, and there was approximately 5 mm of bone on each side of the IVD. Samples were from levels C2–C6 and were 15.5–24.9 mm in diameter and 3.4–10.0 mm in height, and the age of the animals was 6–12 months. Dimensions were measured with a caliper (resolution 0.05 mm). The samples were immersed in phosphate-buffered saline solution (PBS) and put in fridge (+4 °C) until mechanical testing.

Fig. 1
Finite-element model creation and detailed features. (a) The annulus fibrosus (AF) was modeled as a fiber-reinforced poroviscoelastic material, and it incorporated 10 layers of lamellae with alternating fiber orientations (±30 deg with respect to the circumferential direction). (b) The model incorporated the bone shape typical of the bovine tail discs, as well as rigid boundary conditions. (c) The nucleus pulposus was modeled as a poroelastic material, located at the center of the (d) finite element model of the intervertebral disc (IVD). e) The AF had a radial distribution of material properties; the material property coefficient indicates relative value with respect to outer AF, for each lamella layer. (f) An IVD was isolated from a bovine tail and the outcome of the model was compared to (g) an experimental compressive stress-relaxation test.
Fig. 1
Finite-element model creation and detailed features. (a) The annulus fibrosus (AF) was modeled as a fiber-reinforced poroviscoelastic material, and it incorporated 10 layers of lamellae with alternating fiber orientations (±30 deg with respect to the circumferential direction). (b) The model incorporated the bone shape typical of the bovine tail discs, as well as rigid boundary conditions. (c) The nucleus pulposus was modeled as a poroelastic material, located at the center of the (d) finite element model of the intervertebral disc (IVD). e) The AF had a radial distribution of material properties; the material property coefficient indicates relative value with respect to outer AF, for each lamella layer. (f) An IVD was isolated from a bovine tail and the outcome of the model was compared to (g) an experimental compressive stress-relaxation test.
Close modal

Mechanical Testing of Bone-Intervertebral Disc-Bone Segments.

Compressive testing of the bone-IVD-bone segments was conducted with an Instron 5943 dynamic testing machine (Instron, Norwood, MA) equipped with a 1 kN load cell (Fig. 1(g)). Testing was performed at room temperature (approximately +21 °C) and samples were kept moist by spraying PBS onto the surface. The samples were preconditioned with 10 triangular loading-unloading cycles to 4% compressive strain with a 1%/s strain rate, after which they were allowed to recover for 1 min, and a relaxation experiment was started. The samples were subjected to a 4-step compressive stress-relaxation experiment using 4% strain increments (i.e., steps at 4, 8, 12, and 16% compressive strain), a ramp velocity of 1%/s, and 180 s of relaxation after each step. Strains were calculated based on the IVD height measured on the side with a caliper at four locations. The bones were considered rigid since they are much stiffer than the IVD. The compressive force-time graph was recorded and one of the segments, chosen for finite element model creation, was compared with the finite element model to validate the model predictions.

Material Model.

The AF lamellae were modeled as a fiber-reinforced poroviscoelastic homogenized continuum material, consisting of the additive contributions of a porous matrix material fully saturated with water and a reinforcing viscoelastic collagen fiber network (Fig. 1(a)). The total stress tensor σtot is given by
σtot=σm+SpI
(1)
where σm is the matrix material stress tensor, S is the collagen fiber network stress tensor and I is the identity tensor. Poroelastic theory is used to calculate the fluid pressure p [36]. The matrix material was modeled with a neo-Hookean model, where the stress is given by
σm=KlnJJI+GJ(FFTJ2/3I)
(2)
where bulk modulus K and shear modulus G are calculated based on Young's modulus E and Poisson's ratio ν as
K=E3(12ν)
(3)
G=E2(1+ν)
(4)
J is the determinant of the deformation gradient tensor F. The fluid flow within the porous material was modeled using Darcy's law [37]
q=kp
(5)
where the fluid flow rate q is calculated with permeability k and pressure gradient p. The collagen fiber network was modeled as an oriented tension-only viscoelastic material [33], where the stress tensor at time tn+1 is given by
S(tn+1)=μI4(I4γ/21+γ2(1I4))A0,ifI41,0otherwise
(6)
where I4 is the square of the fiber stretch, γ is a material parameter representing the nonlinearity of the fiber network, A0 is a structural tensor depending on fiber network orientation (A0=a0ao; fiber network orientation is ao=(100)T), and μ a time-dependent stiffness parameter given by
μ(tn+1)=μ+βh(tn+1)
(7)
where μ is the elastic modulus of the fibers, β is the elastic coefficient of a Maxwell element of the fiber network and h(tn+1) is given by
h(tn+1)=eΔt/τh(tn)+1eΔt/τΔt/τ(I4(tn+1)I4(tn))
(8)

where τ is a time constant. In summary, the behavior of the fiber network is thus described by four parameters: μ and γ representing the elastic modulus and nonlinearity of the fiber network, and β and τ representing the elastic coefficient and time constant of a Maxwell element.

The NP was modeled as a neo-Hookean poroelastic material, represented by Eqs. (1)(5) described above (Fig. 1(c)).

Finite Element Model of the Bovine Intervertebral Disc.

Out of the seven bovine IVDs extracted and experimentally measured, one representative sample was chosen for FE-model creation and validation. The chosen sample was representative in size to the IVDs used in the bioreactor application and exhibited typical relaxation behavior, in terms of amount and rate of relaxation. The IVDs in the bioreactor experiments are typically 17–25 mm in diameter, with a preference for larger discs to have more material in subsequent histological and biochemical analyses. This one-sample approach was chosen, instead of averaging all the samples as done in a similar study [38], as the disc size has a large effect on the force levels, and averaging such data could lead to inaccurate conclusions.

The IVD was modeled as a cylinder, with AF:NP diameter ratio of 12:5 based on a typical bovine IVD used in our laboratory (Figs. 1(b) and 1(d)). One FE-model was created with Abaqus 6.13-3 (Dassault Systèmes, Johnston, RI) with 3200 C3D8P elements in the AF and 1120 C3D8PH elements in the NP according to sample dimensions (diameter 23.7 mm, height 6.83 mm) and bone shape (Fig. 1(b)). The model had no distinct NP and AF parts; the IVD was a single continuum part, and the elements of NP and AF were assigned different materials. Displacements at the bottom surface nodes of the model were fixed, and the top surface nodes were rigidly tied to a reference point (Fig. 1(b)). The model was subjected to the four-step compressive stress-relaxation test described above, by applying the motion on the reference point, and the resulting force-time graph was recorded. The pore pressure at the outer surface of the AF was set to zero to allow for free fluid flow through the surface.

The AF was modeled with 10 layers of lamellae with alternating fiber orientations of ±30 deg with respect to the circumferential direction. Within the AF, the fiber network properties were assumed to be distributed radially according to the collagen type I content reported in the literature [39] (Fig. 1(e)). The reported contents were 290 mg/g (wet weight) in outer AF, 210 and 150 mg/g inside the AF, and 140 mg/g in inner AF, which were normalized by the content in outer lamella to give a coefficient ranging from 1.0 (outer lamella) to 0.48 (inner lamella). The coefficient at each of the 10 layers was linearly interpolated from the four points where collagen content was reported, corresponding to coefficients of 0.48, 0.52, 0.72, and 1.0 at layers 1, 4, 7, and 10. The material properties μ, γ, and τ at the 10 AF layers were then implemented by giving a value for the property at the outer lamella and multiplying by the coefficient. The material property β was assumed to inversely follow the collagen content distribution [40] so that the coefficient was 1 for inner lamella and 0.48 for outer lamella. The values for the fiber network properties μ, γ, τ, and β were determined by fitting the force-time output of the FE model to the experimentally measured force-time output, i.e., calibrated based on the experimental data. Changing of the parameter values was done systematically [4143] and the accurateness of fit was estimated based on the R2-value, and by inspecting the equilibrium, peak, and relaxation rate responses. The properties of the outer lamella layer were thus μ = 31.0 MPa, γ= 7.1, τ= 7.8 s, and β = 53.1 MPa, and at the inner lamella layer μ = 15.0 MPa, γ= 3.4, τ= 3.8 s, and β = 110 MPa.

The Young's modulus, Poisson's ratio and permeability of the AF matrix (solid volume fraction) were set to EAF = 0.035 MPa [44], νAF = 0.1 [31,45], and kAF = 11 × 10−15 m4/Ns [44], and for the NP they were fixed to ENP= 0.019 MPa [44], νNP= 0.49 [21] and kNP=32 × 10−15 m4/Ns [44].

The model was prestressed by giving the NP an isotropic thermal expansion (strain) of 0.08 and letting the disc equilibrate. This aimed to replicate the natural state of the disc without external loads, in which NP is pressurized and AF fibers are in tension [46], though the prestressed state of the disc (without external loads) is difficult to validate due to incomplete literature data.

Parametric Analyses.

A parametric study was conducted to investigate the effect of different material parameters and modeling choices on the resulting force-time output. These reveal the effect of variation of the material properties, AF fiber angle, model shape, and boundary conditions (rigid boundary conditions with bone shape versus plane-ended shape versus including 5 mm of bone into the model as in the compressive experiment or bone with a cross shape as in the bioreactor setup [20]), and the effect of disc diameter. The range of disc diameters (d = 17 to 25 mm) was chosen based on the IVDs typically used in bioreactor studies; the model created in this study (d = 23.7 mm) is thus among the large discs. The bone was modeled as a linear elastic material with Young's modulus of 220.5 MPa and Poisson's ratio of 0.399, as reported for bovine vertebrae [47]. The conducted analyses are summarized in Table 1.

Table 1

Summary of conducted parametric analyses

Parameter or propertySymbolInvestigated variation (original model value is in bold)Source for the original model valueRationale for the variation investigated
NP Young's modulusENP0.01; 0.019; 0.03 MPaExperimental [44]Variation comparable to the effect of digestion [44]
NP Poisson's ratioνNP0.47; 0.49; 0.495Computational studies [21]Behavior close to the assumption of incompressibility
AF Young's modulusEAF0.025; 0.035; 0.045 MPaExperimental [44]Variation comparable to the effect of digestion [44]
AF Poisson's ratioνAF0.05; 0.1; 0.4Computational studies [31,45]Covering the range measured experimentally for bovine tissue [60]
NP and AF permeability (investigated simultaneously)kNP, kAFkNP = 32 × 10−15 m4/Ns ± 20%Experimental [44]Variation similar to Ref. [61]
kAF = 11 × 10−15 m4/Ns ± 20 %
AF fiber network elastic modulusμ31.0 MPa ± 20 %Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network nonlinearityγ7.1 ±20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network elastic coefficientβ110 MPa ± 20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network time constantτ7.8 s ± 20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network fiber angle25, 30, 35, 40 degExperimental [52]Range in Ref. [52]
IVD model shape and boundary conditionsBone shape with rigid boundary conditions, plane-ended shape with rigid boundary conditions, 5 mm of bone, bone with crossExperimental, typical bovine IVD used in bioreactor experimentsPossible experimental and modeling setups
IVD diameterd17, 19, 21, 23, 23.7, 25 mmExperimental, typical bovine IVD used in bioreactor experimentsRange of diameters used in bioreactor experiments
Parameter or propertySymbolInvestigated variation (original model value is in bold)Source for the original model valueRationale for the variation investigated
NP Young's modulusENP0.01; 0.019; 0.03 MPaExperimental [44]Variation comparable to the effect of digestion [44]
NP Poisson's ratioνNP0.47; 0.49; 0.495Computational studies [21]Behavior close to the assumption of incompressibility
AF Young's modulusEAF0.025; 0.035; 0.045 MPaExperimental [44]Variation comparable to the effect of digestion [44]
AF Poisson's ratioνAF0.05; 0.1; 0.4Computational studies [31,45]Covering the range measured experimentally for bovine tissue [60]
NP and AF permeability (investigated simultaneously)kNP, kAFkNP = 32 × 10−15 m4/Ns ± 20%Experimental [44]Variation similar to Ref. [61]
kAF = 11 × 10−15 m4/Ns ± 20 %
AF fiber network elastic modulusμ31.0 MPa ± 20 %Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network nonlinearityγ7.1 ±20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network elastic coefficientβ110 MPa ± 20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network time constantτ7.8 s ± 20%Calibrated to match experimentVariation similar to Ref. [61]
AF fiber network fiber angle25, 30, 35, 40 degExperimental [52]Range in Ref. [52]
IVD model shape and boundary conditionsBone shape with rigid boundary conditions, plane-ended shape with rigid boundary conditions, 5 mm of bone, bone with crossExperimental, typical bovine IVD used in bioreactor experimentsPossible experimental and modeling setups
IVD diameterd17, 19, 21, 23, 23.7, 25 mmExperimental, typical bovine IVD used in bioreactor experimentsRange of diameters used in bioreactor experiments

NP = nucleus pulposus, AF = annulus fibrosus, IVD = intervertebral disc.

Results

The model captured well the experimentally observed equilibrium behavior at all strain levels, as well as the force peaks and relaxation phenomena at higher strain levels (12 and 16% strain, Fig. 2(a)). The equilibrium forces of the model and experiment were 5.8 and 9.8 N (1st step), 20.3 and 25.0 N (2nd step), 45.6 and 44.3 N (3rd step) and 69.4 and 68.2 N (4th step). The peak forces of the model and experiment were 93.7 and 119.2 N (3rd step) and 205.8 and 184.0 N (4th step). The initial rate of relaxation of the model at the 4th step was −7.5 N/s and of the experiment −7.7 N/s. The maximum principal strains during 3rd loading step were high in NP and in AF close to endplates (Fig. 2(b)), while the maximum principal stresses were highest in the AF (Fig. 2(c)). The pore pressure and fluid flow showed fluid outflow (Figs. 2(d) and 2(e)).

Fig. 2
(a) Comparison of the IVD finite element (FE) model and compressive relaxation experiment force-time graphs. (b)Maximum principal strain, (c) maximum principal stress, (d) pore pressure, and (e) fluid velocity within the IVD are shown during the 3rd step of the simulated experiment. Third step was chosen for visualization since the compressive strain level of 12% is close to the disc height change under physiological loading [6].
Fig. 2
(a) Comparison of the IVD finite element (FE) model and compressive relaxation experiment force-time graphs. (b)Maximum principal strain, (c) maximum principal stress, (d) pore pressure, and (e) fluid velocity within the IVD are shown during the 3rd step of the simulated experiment. Third step was chosen for visualization since the compressive strain level of 12% is close to the disc height change under physiological loading [6].
Close modal

The parametric study revealed that increasing the Young's modulus and the Poisson's ratio of NP or AF matrix material increases slightly the overall force levels (Fig. 3). Increasing permeability decreases the peak forces and makes relaxation considerably faster.

Fig. 3
Parametric study showing the effect of Young's modulus and Poisson's ratio of NP and AF matrix, and the effect of permeability of the tissue on the compressive force-time graph. The graph of the original model is shown thickened (blue), and arrows indicate the direction of change when the parameter value is increased. (Color version online.)
Fig. 3
Parametric study showing the effect of Young's modulus and Poisson's ratio of NP and AF matrix, and the effect of permeability of the tissue on the compressive force-time graph. The graph of the original model is shown thickened (blue), and arrows indicate the direction of change when the parameter value is increased. (Color version online.)
Close modal

Increasing the AF fiber network elastic stiffness parameters μ and γ highly increases the peak force levels and the following relaxation (Fig. 4). The AF fiber network elastic coefficient β minorly affects the peak force levels, while higher time constant τ makes initial relaxation slower. Altering fiber network parameters μ, γ, β, or τ did not considerably change the behavior at low strain levels. Fiber angle affects considerably the force levels; higher fiber angles result in lower forces (Fig. 5). Taking bone shape into account and using rigid boundary conditions in the model increases the force levels when compared with plane-ended shapes or models with elastic bone. Increasing model diameter leads to considerably increased force and stress levels, e.g., 15 N versus 54 N when comparing 17 and 25 mm IVDs at 3rd step equilibrium, corresponding to average compressive stresses of 0.066 and 0.11 MPa on the disc.

Fig. 4
Parametric study showing the effect of AF fiber network elastic stiffness parameters μ∞ and γ, and coefficient β and time constant τ on the compressive force-time graph. The graph of the original model is shown thickened (blue), and arrows indicate the direction of change when the parameter value is increased. (Color version online.)
Fig. 4
Parametric study showing the effect of AF fiber network elastic stiffness parameters μ∞ and γ, and coefficient β and time constant τ on the compressive force-time graph. The graph of the original model is shown thickened (blue), and arrows indicate the direction of change when the parameter value is increased. (Color version online.)
Close modal
Fig. 5
Analyses showing the effects of fiber angle, model shape and boundary conditions and IVD diameter on the resulting compressive force-time graph. The following shape and boundary conditions (BCs) were considered: bone shape with rigid BCs (original model, thickened, blue), plane-ended shape with rigid BCs (red), model including 5 mm of bone as in the compressive experiment (yellow) and model including 7 mm of bone with a cross shape as in the bioreactor setup (purple). (Color version online.)
Fig. 5
Analyses showing the effects of fiber angle, model shape and boundary conditions and IVD diameter on the resulting compressive force-time graph. The following shape and boundary conditions (BCs) were considered: bone shape with rigid BCs (original model, thickened, blue), plane-ended shape with rigid BCs (red), model including 5 mm of bone as in the compressive experiment (yellow) and model including 7 mm of bone with a cross shape as in the bioreactor setup (purple). (Color version online.)
Close modal

Discussion

The purpose of this study was to establish a finite element model of a bovine intervertebral disc to be able to evaluate the global and local responses of the tissue under various loads, that are physiologically relevant and will be applied in a new bioreactor. The proposed model was able to capture the equilibrium response at all strain levels, as well as the relaxation at higher strain levels, which is the intended application of this model. The model can be used to evaluate the response of the bovine IVD to aid in the design of a new bioreactor system and to help in future experiment planning.

The proposed model was able to capture the typical phenomenological behavior of an IVD in compression: the NP expands perpendicular to the load, pushing the AF, resulting in the lamellae being in tension in the fiber direction. Additionally, the model captures the fluid flow out of the tissue, which is an important determinant for the time-dependent compressive behavior. The model captures well the experimentally observed nonlinear elastic and viscoelastic behavior at the intended physiological and degenerative strain levels.

The high maximum principal strains during 3rd loading step close to endplates can be attributed to high strains resulting from NP expansion while the surfaces have rigid boundary conditions, allowing no deformation at the surfaces, leading to high shearing. The high maximum principal stresses in the AF, mostly due to extension in the stiff fiber direction, demonstrate its load-carrying function, and how the distribution of material properties results in a rather even distribution of the stresses within AF, as also demonstrated in human discs [25].

In the parametric analyses, it was found that the Poisson's ratios of NP or AF matrix did not have a major effect on the overall response. This can be attributed to the poroelasticity in the model: under fast compression, the model is initially incompressible (ν = 0.5, fluid has had no time to flow out), while without poroelasticity the Poisson's ratio would have a large effect. The Young's moduli of NP and AF matrix have likewise a small, almost negligible effect as they are small relative to the material stiffness in the lamella fiber direction. The permeability, on the other hand, has a large effect on the relaxation response with higher permeability allowing faster fluid flow and thus faster relaxation. Any changes in permeability, for example, due to degeneration, thus affect largely the dynamic behavior of the disc.

The elastic stiffness of the fiber network, represented by μ and γ, was affecting most the overall response of the disc. Higher values lead to higher peak forces and higher magnitude and rate of relaxation. When the disc is compressed, the fibers are strained relatively fast to high magnitudes due to effectively incompressible instantaneous response, resulting in high peak forces, but as the fluid flows out, the tissue contracts and fiber strains become lower. This explains the high effect of elastic fiber parameters on the overall elastic and viscoelastic response.

The fiber network viscoelastic parameters β and τ had an effect only in the initial relaxation response lasting for a few seconds. One could argue that viscoelasticity of the collagen fiber network could be omitted from the model to simplify the material, as has been done in IVD [2326] and, for example, in cartilage [48] and ligaments [49]. However, collagen fibers and the fiber network are viscoelastic [50,51], therefore viscoelasticity is considered in this model. Under different dynamic loading, for example, cyclic sinusoidal loading, the effect of β and τ could be higher.

The fiber angle, assumed to be 30 deg in this study, seems to affect the force levels for a given step with a relatively constant magnitude. Given the rather large effect of fiber angle, making this assumption more accurate with experimental measurement of the angles, for example, with polarized light microscopy [52], would increase the accuracy of the model and could capture the low-strain behavior better.

The cylindrical model without the bone shape but with rigid boundary conditions lowered considerably the force levels. This can be attributed to bone shape promoting the NP expansion toward the AF and resulting in a lower average height, effectively increasing the axial stiffness. Including either 5 mm of bone or the bone with a cross decreased the peak force levels to a small extent. This is attributed to the displacement-controlled test, in which the compression of the bone leads to lower disc strains compared with the rigid boundary conditions. The results suggest that the model's assumption of rigid boundary conditions is sufficient for the purpose of the model.

The disc diameter had a considerable effect on the force and average compressive stress levels. Although the influence on the forces was expected, having more than three times higher forces in a 25 mm disc than in a 17 mm disc might have significant implications on the control system design (in addition to natural variation). For example, in force-controlled tests, the control loop should be able to adapt sample-by-sample. These results also suggest that it is highly important to have the correct dimensions when doing sample-specific FE-modeling, also for the evaluation of the tissue stresses, as smaller diameters have been shown to increase AF fiber stresses [22,53], and smaller height stiffens the disc response [54,55].

The model presented in this study will aid in future bioreactor experiment planning. For example, one could compare different loading magnitudes to estimate whether changes would lead to large differences in terms of tissue strains or stresses, i.e., mechanical stimulus to the cells, if deviating from a certain loading protocol, such as the established physiological compressive loading of 0.02–0.2 MPa. Experimentally such investigations are laborious and time-consuming, while the model presented in this study could be used for a rapid assessment.

Limitations.

This study comes with certain limitations. There was only one experimentally measured IVD considered in this study, however, it was chosen as a representative sample out of seven measured IVDs. The parametric analyses address this limitation by giving a clear idea of how variation affects the IVD response. The FE model proposed in this study would capture the behavior of different discs, while the material property values could be different. If using the material property values defined in this study to model other discs, the response could differ due to natural variation. The material property values of the fiber network were assumed to be linearly dependent on the collagen content. However, in addition to the collagen content the microstructural organization and interactions with other constituents influence the material property values. Taking the distribution within AF into account has been shown to yield accurate results in human IVD models in comparison to constant values [25]. The permeabilities of NP and AF were assumed isotropic, though especially in the AF they may be different in fiber and transverse directions [44]. The model investigated in this study was rotationally symmetric and under compression, therefore flow in other directions than radial is minimal. The endplates in the model were assumed impermeable, though fluid flow through the endplates would occur to some extent during the relaxation times investigated in this study [56]. Taking this into account could make the relaxation under compressive loading slightly faster. The AF:NP diameter ratio was assumed to be 12:5, though the ratio may vary between IVDs of different size and the boundary between AF and NP may be difficult to determine accurately. The shapes of the IVD and NP were assumed to be cylindrical, an assumption that was accurate for our young IVDs and was also used in other studies with bovine IVDs [39,57]. However, with larger and more mature discs the NP shape may tend toward an ellipse [58], which results in different behavior in flexion–extension and lateral bending directions and could be considered in the future. The four parameters of the fiber network μ, γ, τ, and β that were calibrated using the compressive experimental data could be nonunique (another parameter combination could produce the same compressive force-time output). This uncertainty of nonuniqueness is mitigated to a large extent by the parametric study conducted. Nevertheless, we cannot exclude the possibility that the shown values, if nonunique, could be unable to capture the IVD response under other motions, such as axial torsion. The peak forces were under-predicted for the first two steps, likely due to nonlinear behavior of the AF fibers, as with the first two steps the AF strains and therefore material stiffness are lower than with the last two steps, influencing the IVD deformation and fluid outflow. The model should be validated against other motions and tissue strains; nonetheless, the model showed the ability to differentiate between physiological and degenerative compressive loading in a force-controlled setup, as well as other loads [59].

Conclusion.

The FE model presented in this study was able to capture the complex time-dependent phenomena of IVD under compression at physiological strain levels. The relaxation was mostly caused by the fluid flow out of the tissue, which was highly affected by the elastic properties of the AF fiber network and permeability. The model can be used to predict the response of bovine IVDs to aid in bioreactor system design and future experiment planning.

Funding Data

  • Swiss National Science Foundation (Grant No. 189915; Funder ID: 10.13039/501100001711).

  • AO Foundation (Funder ID: 10.13039/501100001702).

  • AO Spine (Funder ID: 10.13039/501100002732).

Data Availability Statement

The authors attest that all data for this study are included in the paper.

Nomenclature

a0 =

fiber network orientation

A0 =

structural tensor

d =

intervertebral disc diameter

E =

Young's modulus

F =

deformation gradient tensor

G =

shear modulus

I =

identity tensor

I4 =

square of the fiber stretch

J =

determinant of the deformation gradient tensor

K =

bulk modulus

k =

permeability

p =

pressure

q =

fluid flowrate

S =

stress tensor of the fiber network

β =

elastic coefficient of a Maxwell element of the fiber network

γ =

nonlinearity of the fiber network

μ =

time-dependent stiffness parameter

μ =

elastic modulus of the fibers

ν =

Poisson's ratio

σm =

matrix material stress tensor

σtot =

total stress tensor

τ =

time constant of the fiber network

References

1.
Teraguchi
,
M.
,
Yoshimura
,
N.
,
Hashizume
,
H.
,
Muraki
,
S.
,
Yamada
,
H.
,
Minamide
,
A.
,
Oka
,
H.
,
Ishimoto
,
Y.
,
Nagata
,
K.
,
Kagotani
,
R.
,
Takiguchi
,
N.
,
Akune
,
T.
,
Kawaguchi
,
H.
,
Nakamura
,
K.
, and
Yoshida
,
M.
,
2014
, “
Prevalence and Distribution of Intervertebral Disc Degeneration Over the Entire Spine in a Population-Based Cohort: The Wakayama Spine Study
,”
Osteoarthr. Cartil.
,
22
(
1
), pp.
104
110
.10.1016/j.joca.2013.10.019
2.
Gantenbein
,
B.
,
Grünhagen
,
T.
,
Lee
,
C. R.
,
Van Donkelaar
,
C. C.
,
Alini
,
M.
, and
Ito
,
K.
,
2006
, “
An In Vitro Organ Culturing System for Intervertebral Disc Explants With Vertebral Endplates: A Feasibility Study With Ovine Caudal Discs
,”
Spine (Phila. Pa. 1976)
,
31
(
23
), pp.
2665
2673
.10.1097/01.brs.0000244620.15386.df
3.
Gantenbein
,
B.
,
Illien-Jünger
,
S.
,
Chan
,
S.
,
Walser
,
J.
,
Haglund
,
L.
,
Ferguson
,
S.
,
Iatridis
,
J.
, and
Grad
,
S.
,
2015
, “
Organ Culture Bioreactors – Platforms to Study Human Intervertebral Disc Degeneration and Regenerative Therapy
,”
Curr. Stem Cell Res. Ther.
,
10
(
4
), pp.
339
352
.10.2174/1574888X10666150312102948
4.
Haglund
,
L.
,
Moir
,
J.
,
Beckman
,
L.
,
Mulligan
,
K. R.
,
Jim
,
B.
,
Ouellet
,
J. A.
,
Roughley
,
P.
, and
Steffen
,
T.
,
2011
, “
Development of a Bioreactor for Axially Loaded Intervertebral Disc Organ Culture
,”
Tissue Eng. Part C Methods
,
17
(
10
), pp.
1011
1019
.10.1089/ten.tec.2011.0025
5.
Lang
,
G.
,
Liu
,
Y.
,
Geries
,
J.
,
Zhou
,
Z.
,
Kubosch
,
D.
,
Südkamp
,
N.
,
Richards
,
R. G.
,
Alini
,
M.
,
Grad
,
S.
, and
Li
,
Z.
,
2018
, “
An Intervertebral Disc Whole Organ Culture System to Investigate Proinflammatory and Degenerative Disc Disease Condition
,”
J. Tissue Eng. Regen. Med.
,
12
(
4
), pp.
e2051
e2061
.10.1002/term.2636
6.
Li
,
Z.
,
Gehlen
,
Y.
,
Heizmann
,
F.
,
Grad
,
S.
,
Alini
,
M.
,
Richards
,
R. G.
,
Kubosch
,
D.
,
Südkamp
,
N.
,
Izadpanah
,
K.
,
Kubosch
,
E. J.
, and
Lang
,
G.
,
2020
, “
Preclinical Ex-Vivo Testing of Anti-Inflammatory Drugs in a Bovine Intervertebral Degenerative Disc Model
,”
Front. Bioeng. Biotechnol.
,
8
(
June
), pp.
1
23
.10.3389/fbioe.2020.00001
7.
Likhitpanichkul
,
M.
,
Dreischarf
,
M.
,
Illien-Junger
,
S.
,
Walter
,
B. A.
,
Nukaga
,
T.
,
Long
,
R. G.
,
Sakai
,
D.
,
Hecht
,
A. C.
, and
Iatridis
,
J. C.
,
2014
, “
Fibrin-Genipin Adhesive Hydrogel for Annulus Fibrosus Repair: Performance Evaluation With Large Animal Organ Culture, In Situ Biomechanics, and In Vivo Degradation Tests
,”
Eur. Cells Mater.
,
28
, pp.
25
38
.10.22203/eCM.v028a03
8.
Paul
,
C. P. L.
,
Schoorl
,
T.
,
Zuiderbaan
,
H. A.
,
Zandieh Doulabi
,
B.
,
van der Veen
,
A. J.
,
van de Ven
,
P. M.
,
Smit
,
T. H.
,
van Royen
,
B. J.
,
Helder
,
M. N.
, and
Mullender
,
M. G.
,
2013
, “
Dynamic and Static Overloading Induce Early Degenerative Processes in Caprine Lumbar Intervertebral Discs
,”
PLoS One
,
8
(
4
), p.
e62411
.10.1371/journal.pone.0062411
9.
Pfannkuche
,
J. J.
,
Guo
,
W.
,
Cui
,
S.
,
Ma
,
J.
,
Lang
,
G.
,
Peroglio
,
M.
,
Richards
,
R. G.
,
Alini
,
M.
,
Grad
,
S.
, and
Li
,
Z.
,
2020
, “
Intervertebral Disc Organ Culture for the Investigation of Disc Pathology and Regeneration–Benefits, Limitations, and Future Directions of Bioreactors
,”
Connect. Tissue Res.
,
61
(
3–4
), pp.
304
321
.10.1080/03008207.2019.1665652
10.
Korecki
,
C. L.
,
MacLean
,
J. J.
, and
Iatridis
,
J. C.
,
2008
, “
Dynamic Compression Effects on Intervertebral Disc Mechanics and Biology
,”
Spine (Phila. Pa. 1976)
,
33
(
13
), pp.
1403
1409
.10.1097/BRS.0b013e318175cae7
11.
Dudli
,
S.
,
Haschtmann
,
D.
, and
Ferguson
,
S. J.
,
2015
, “
Persistent Degenerative Changes in the Intervertebral Disc After Burst Fracture in an In Vitro Model Mimicking Physiological Post-Traumatic Conditions
,”
Eur. Spine J.
,
24
(
9
), pp.
1901
1908
.10.1007/s00586-014-3301-3
12.
Hartman
,
R. A.
,
Yurube
,
T.
,
Ngo
,
K.
,
Merzlak
,
N. E.
,
Debski
,
R. E.
,
Brown
,
B. N.
,
Kang
,
J. D.
, and
Sowa
,
G. A.
,
2015
, “
Biological Responses to Flexion/Extension in Spinal Segments Ex-Vivo
,”
J. Orthop. Res.
,
33
(
8
), pp.
1255
1264
.10.1002/jor.22900
13.
Ariga
,
K.
,
Yonenobu
,
K.
,
Nakase
,
T.
,
Hosono
,
N.
,
Okuda
,
S.
,
Meng
,
W.
,
Tamura
,
Y.
, and
Yoshikawa
,
H.
,
2003
, “
Mechanical Stress-Induced Apoptosis of Endplate Chondrocytes in Organ-Cultured Mouse Intervertebral Discs: An Ex Vivo Study
,”
Spine (Phila. Pa. 1976)
,
28
(
14
), pp.
1528
1533
.10.1097/01.BRS.0000076915.55939.E3
14.
Walter
,
B. A.
,
Korecki
,
C. L.
,
Purmessur
,
D.
,
Roughley
,
P. J.
,
Michalek
,
A. J.
, and
Iatridis
,
J. C.
,
2011
, “
Complex Loading Affects Intervertebral Disc Mechanics and Biology
,”
Osteoarthr. Cartil.
,
19
(
8
), pp.
1011
1018
.10.1016/j.joca.2011.04.005
15.
Beatty
,
A. M.
,
Bowden
,
A. E.
, and
Bridgewater
,
L. C.
,
2016
, “
Functional Validation of a Complex Loading Whole Spinal Segment Bioreactor Design
,”
ASME J. Biomech. Eng.
,
138
(
6
), p.
064501
.
16.
Chan
,
S. C. W.
,
Ferguson
,
S. J.
,
Wuertz
,
K.
, and
Gantenbein-Ritter
,
B.
,
2011
, “
Biological Response of the Intervertebral Disc to Repetitive Short-Term Cyclic Torsion
,”
Spine (Phila. Pa. 1976)
,
36
(
24
), pp.
2021
2030
.10.1097/BRS.0b013e318203aea5
17.
Chan
,
S. C. W.
,
Walser
,
J.
,
Ferguson
,
S. J.
, and
Gantenbein
,
B.
,
2015
, “
Duration-Dependent Influence of Dynamic Torsion on the Intervertebral Disc: An Intact Disc Organ Culture Study
,”
Eur. Spine J.
,
24
(
11
), pp.
2402
2410
.10.1007/s00586-015-4140-6
18.
Chan
,
S. C. W.
,
Walser
,
J.
,
Käppeli
,
P.
,
Shamsollahi
,
M. J.
,
Ferguson
,
S. J.
, and
Gantenbein-Ritter
,
B.
,
2013
, “
Region Specific Response of Intervertebral Disc Cells to Complex Dynamic Loading: An Organ Culture Study Using a Dynamic Torsion-Compression Bioreactor
,”
PLoS One
,
8
(
8
), p.
e72489
.10.1371/journal.pone.0072489
19.
Frauchiger
,
D. A.
,
Chan
,
S. C. W.
,
Benneker
,
L. M.
, and
Gantenbein
,
B.
,
2018
, “
Intervertebral Disc Damage Models in Organ Culture: A Comparison of Annulus Fibrosus Cross-Incision Versus Punch Model Under Complex Loading
,”
Eur. Spine J.
,
27
(
8
), pp.
1785
1797
.10.1007/s00586-018-5638-5
20.
Šećerović
,
A.
,
Ristaniemi
,
A.
,
Cui
,
S.
,
Li
,
Z.
,
Soubrier
,
A.
,
Alini
,
M.
,
Ferguson
,
S. J.
,
Weder
,
G.
,
Heub
,
S.
,
Ledroit
,
D.
, and
Grad
,
S.
,
2022
, “
Toward the Next Generation of Spine Bioreactors: Validation of an Ex Vivo Intervertebral Disc Organ Model and Customized Specimen Holder for Multiaxial Loading
,”
ACS Biomater. Sci. Eng.
,
8
(
9
), pp.
3969
3976
.10.1021/acsbiomaterials.2c00330
21.
Dreischarf
,
M.
,
Zander
,
T.
,
Shirazi-Adl
,
A.
,
Puttlitz
,
C. M.
,
Adam
,
C. J.
,
Chen
,
C. S.
,
Goel
,
V. K.
,
Kiapour
,
A.
,
Kim
,
Y. H.
,
Labus
,
K. M.
,
Little
,
J. P.
,
Park
,
W. M.
,
Wang
,
Y. H.
,
Wilke
,
H. J.
,
Rohlmann
,
A.
, and
Schmidt
,
H.
,
2014
, “
Comparison of Eight Published Static Finite Element Models of the Intact Lumbar Spine: Predictive Power of Models Improves When Combined Together
,”
J. Biomech.
,
47
(
8
), pp.
1757
1766
.10.1016/j.jbiomech.2014.04.002
22.
Natarajan
,
R. N.
, and
Andersson
,
G. B. J.
,
1999
, “
The Influence of Lumbar Disc Height and Cross-Sectional Area on the Mechanical Response of the Disc to Physiological Loading
,”
Spine (Phila. Pa. 1976)
,
24
(
18
), pp.
1873
1881
.10.1097/00007632-199909150-00003
23.
Sharabi
,
M.
,
Levi-Sasson
,
A.
,
Wolfson
,
R.
,
Wade
,
K. R.
,
Galbusera
,
F.
,
Benayahu
,
D.
,
Hans-Joachim
,
W.
, and
Haj-Ali
,
R.
,
2019
, “
The Mechanical Role of the Radial Fiber Network Within the Annulus Fibrosus of the Lumbar Intervertebral Disc: A Finite Elements Study
,”
ASME J. Biomech. Eng.
,
141
(
2
), p.
021006
.10.1115/1.4041769
24.
Sharabi
,
M.
,
Wertheimer
,
S.
,
Wade
,
K. R.
,
Galbusera
,
F.
,
Benayahu
,
D.
,
Wilke
,
H. J.
, and
Haj-Ali
,
R.
,
2019
, “
Towards Intervertebral Disc Engineering: Bio-Mimetics of Form and Function of the Annulus Fibrosus Lamellae
,”
J. Mech. Behav. Biomed. Mater.
,
94
(
February
), pp.
298
307
.10.1016/j.jmbbm.2019.03.023
25.
Marini
,
G.
, and
Ferguson
,
S. J.
,
2014
, “
Modelling the Influence of Heterogeneous Annulus Material Property Distribution on Intervertebral Disk Mechanics
,”
Ann. Biomed. Eng.
,
42
(
8
), pp.
1760
1772
.10.1007/s10439-014-1025-5
26.
Marini
,
G.
,
Studer
,
H.
,
Huber
,
G.
,
Püschel
,
K.
, and
Ferguson
,
S. J.
,
2016
, “
Geometrical Aspects of Patient-Specific Modelling of the Intervertebral Disc: Collagen Fibre Orientation and Residual Stress Distribution
,”
Biomech. Model. Mechanobiol.
,
15
(
3
), pp.
543
560
.10.1007/s10237-015-0709-6
27.
Yang
,
B.
, and
O'Connell
,
G. D.
,
2019
, “
Intervertebral Disc Swelling Maintains Strain Homeostasis Throughout the Annulus Fibrosus: A Finite Element Analysis of Healthy and Degenerated Discs
,”
Acta Biomater.
,
100
, pp.
61
74
.10.1016/j.actbio.2019.09.035
28.
Barthelemy
,
V. M. P.
,
van Rijsbergen
,
M. M.
,
Wilson
,
W.
,
Huyghe
,
J. M.
,
van Rietbergen
,
B.
, and
Ito
,
K.
,
2016
, “
A Computational Spinal Motion Segment Model Incorporating a Matrix Composition-Based Model of the Intervertebral Disc
,”
J. Mech. Behav. Biomed. Mater.
,
54
, pp.
194
204
.10.1016/j.jmbbm.2015.09.028
29.
Jacobs
,
N. T.
,
Cortes
,
D. H.
,
Peloquin
,
J. M.
,
Vresilovic
,
E. J.
, and
Elliott
,
D. M.
,
2014
, “
Validation and Application of an Intervertebral Disc Finite Element Model Utilizing Independently Constructed Tissue-Level Constitutive Formulations That Are Nonlinear, Anisotropic, and Time-Dependent
,”
J. Biomech.
,
47
(
11
), pp.
2540
2546
.10.1016/j.jbiomech.2014.06.008
30.
Kandil
,
K.
,
Zaïri
,
F.
,
Messager
,
T.
, and
Zaïri
,
F.
,
2020
, “
Interlamellar Matrix Governs Human Annulus Fibrosus Multiaxial Behavior
,”
Sci. Rep.
,
10
(
1
), pp.
1
14
.10.1038/s41598-020-74107-8
31.
Swider
,
P.
,
Pédrono
,
A.
,
Ambard
,
D.
,
Accadbled
,
F.
, and
Sales de Gauzy
,
J.
,
2010
, “
Substructuring and Poroelastic Modelling of the Intervertebral Disc
,”
J. Biomech.
,
43
(
7
), pp.
1287
1291
.10.1016/j.jbiomech.2010.01.006
32.
Schroeder
,
Y.
,
Wilson
,
W.
,
Huyghe
,
J. M.
, and
Baaijens
,
F. P. T.
,
2006
, “
Osmoviscoelastic Finite Element Model of the Intervertebral Disc
,”
Eur. Spine J.
,
15
(
Suppl. 3
), pp.
361
371
.https://www.academia.edu/27175485/Osmoviscoelastic_finite_element_model_of_the_intervertebral_disc
33.
Marini
,
G.
,
2015
,
Dynamic Characterization of the Intervertebral Disc
,
ETH Zürich
,
Zürich, Switzerland
.
34.
Ristaniemi
,
A.
,
Šećerović
,
A.
,
Grad
,
S.
, and
Ferguson
,
S. J.
,
2022
, “
A Novel Viscoelastic Bovine Intervertebral Disc Finite Element Model
,”
eCM Period
, Collection 1, p.
30
.
35.
Saravi
,
B.
,
Lang
,
G.
,
Grad
,
S.
,
Alini
,
M.
,
Richards
,
R. G.
,
Schmal
,
H.
,
Südkamp
,
N.
, and
Li
,
Z.
,
2021
, “
A Proinflammatory, Degenerative Organ Culture Model to Simulate Early-Stage Intervertebral Disc Disease
,”
JoVE
, (
168
), p.
e62100
.10.3791/62100-v
36.
Detournay
,
E.
, and
Cheng
,
A. H.
,
1994
, “
Fundamentals of Poroelasticity
,”
Int. J. Rock Mech. Min. Sci. Geomech. Abstr.
,
31
(
3
), pp.
138
139
.
37.
Holmes
,
M. H.
, and
Mow
,
V. C.
,
1990
, “
The Nonlinear Characteristics of Soft Gels and Hydrated Connective Tissues in Ultrafiltration
,”
J. Biomech.
,
23
(
11
), pp.
1145
1156
.10.1016/0021-9290(90)90007-P
38.
Ristaniemi
,
A.
,
Tanska
,
P.
,
Stenroth
,
L.
,
Finnilä
,
M. A. J.
, and
Korhonen
,
R. K.
,
2021
, “
Comparison of Material Models for Anterior Cruciate Ligament in Tension: From Poroelastic to a Novel Fibril-Reinforced Nonlinear Composite Model
,”
J. Biomech.
,
114
, p.
110141
.10.1016/j.jbiomech.2020.110141
39.
Bezci
,
S. E.
,
Werbner
,
B.
,
Zhou
,
M.
,
Malollari
,
K. G.
,
Dorlhiac
,
G.
,
Carraro
,
C.
,
Streets
,
A.
, and
O'Connell
,
G. D.
,
2019
, “
Radial Variation in Biochemical Composition of the Bovine Caudal Intervertebral Disc
,”
JOR Spine
,
2
(
3
), pp.
1
13
.10.1002/jsp2.1065
40.
Ristaniemi
,
A.
,
Šećerović
,
A.
,
Ferguson
,
S. J.
, and
Grad
,
S.
,
2022
, “
Viscoelastic Characterization of Bovine Annulus Fibrosus Lamellae
,”
9th World Congress of Biomechanics
, Taipei, Taiwan, July 10–14, p.
O-10067
.
41.
Korhonen
,
R. K.
,
Laasanen
,
M. S.
,
Töyräs
,
J.
,
Lappalainen
,
R.
,
Helminen
,
H. J.
, and
Jurvelin
,
J. S.
,
2003
, “
Fibril Reinforced Poroelastic Model Predicts Specifically Mechanical Behavior of Normal, Proteoglycan Depleted and Collagen Degraded Articular Cartilage
,”
J. Biomech.
,
36
(
9
), pp.
1373
1379
.10.1016/S0021-9290(03)00069-1
42.
Li
,
L. P.
, and
Herzog
,
W.
,
2004
, “
The Role of Viscoelasticity of Collagen Fibers in Articular Cartilage: Theory and Numerical Formulation
,”
Biorheology
,
41
(
3–4
), pp.
181
194
.https://content.iospress.com/articles/biorheology/bir297
43.
Li
,
L. P.
,
Herzog
,
W.
,
Korhonen
,
R. K.
, and
Jurvelin
,
J. S.
,
2005
, “
The Role of Viscoelasticity of Collagen Fibers in Articular Cartilage: Axial Tension Versus Compression
,”
Med. Eng. Phys.
,
27
(
1
), pp.
51
57
.10.1016/j.medengphy.2004.08.009
44.
Recuerda
,
M.
,
Périé
,
D.
,
Gilbert
,
G.
, and
Beaudoin
,
G.
,
2012
, “
Assessment of Mechanical Properties of Isolated Bovine Intervertebral Discs From Multi-Parametric Magnetic Resonance Imaging
,”
BMC Musculoskelet. Disord.
,
13
, pp.
1
14
.10.1186/1471-2474-13-195
45.
Argoubi
,
M.
, and
Shirazi-Adl
,
A.
,
1996
, “
Poroelastic Creep Response Analysis of a Lumbar Motion Segment in Compression
,”
J. Biomech.
,
29
(
10
), pp.
1331
1339
.10.1016/0021-9290(96)00035-8
46.
Duclos
,
S. E.
, and
Michalek
,
A. J.
,
2017
, “
Residual Strains in the Intervertebral Disc Annulus Fibrosus Suggest Complex Tissue Remodeling in Response to in-Vivo Loading
,”
J. Mech. Behav. Biomed. Mater.
,
68
(
February
), pp.
232
238
.10.1016/j.jmbbm.2017.02.010
47.
Lim
,
T. H.
, and
Hong
,
J. H.
,
2000
, “
Poroelastic Properties of Bovine Vertebral Trabecular Bone
,”
J. Orthop. Res.
,
18
(
4
), pp.
671
677
.10.1002/jor.1100180421
48.
Ebrahimi
,
M.
,
Ojanen
,
S.
,
Mohammadi
,
A.
,
Finnilä
,
M. A.
,
Joukainen
,
A.
,
Kröger
,
H.
,
Saarakkala
,
S.
,
Korhonen
,
R. K.
, and
Tanska
,
P.
,
2019
, “
Elastic, Viscoelastic and Fibril-Reinforced Poroelastic Material Properties of Healthy and Osteoarthritic Human Tibial Cartilage
,”
Ann. Biomed. Eng.
,
47
(
4
), pp.
953
966
.10.1007/s10439-019-02213-4
49.
Weiss
,
J. A.
,
Maker
,
B. N.
, and
Govindjee
,
S.
,
1996
, “
Finite Element Implementation of Incompressible, Transversely Isotropic Hyperelasticity
,”
Comput. Methods Appl. Mech. Eng.
,
135
(
1–2
), pp.
107
128
.10.1016/0045-7825(96)01035-3
50.
Shen
,
Z. L.
,
Kahn
,
H.
,
Ballarini
,
R.
, and
Eppell
,
S. J.
,
2011
, “
Viscoelastic Properties of Isolated Collagen Fibrils
,”
Biophys. J.
,
100
(
12
), pp.
3008
3015
.10.1016/j.bpj.2011.04.052
51.
Vergari
,
C.
,
Chan
,
D.
,
Clarke
,
A.
,
Mansfield
,
J. C.
,
Meakin
,
J. R.
, and
Winlove
,
P. C.
,
2017
, “
Bovine and Degenerated Human Annulus Fibrosus: A Microstructural and Micromechanical Comparison
,”
Biomech. Model. Mechanobiol.
,
16
(
4
), pp.
1475
1484
.10.1007/s10237-017-0900-z
52.
Michalek
,
A. J.
,
2019
, “
A Growth-Based Model for the Prediction of Fiber Angle Distribution in the Intervertebral Disc Annulus Fibrosus
,”
Biomech. Model. Mechanobiol.
,
18
(
5
), pp.
1363
1369
.10.1007/s10237-019-01150-4
53.
Lu
,
M. Y.
,
Hutton
,
W. C.
, and
Gharpuray
,
V. M.
,
1996
, “
Can Variations in Intervertebral Disc Height Affect the Mechanical Function of the Disc?
,”
Spine (Phila. Pa. 1976)
,
21
(
19
), pp.
2208
2216
.10.1097/00007632-199610010-00006
54.
Robin
,
S.
,
Skalli
,
W.
, and
Lavaste
,
E.
,
1994
, “
Influence of Geometrical Factors on the Behavior of Lumbar Spine Segments: A Finite Element Analysis
,”
Eur. Spine J.
,
3
(
2
), pp.
84
90
.10.1007/BF02221445
55.
Meijer
,
G. J. M.
,
Homminga
,
J.
,
Veldhuizen
,
A. G.
, and
Verkerke
,
G. J.
,
2011
, “
Influence of Interpersonal Geometrical Variation on Spinal Motion Segment Stiffness
,”
Spine (Phila. Pa. 1976)
,
36
(
14
), pp.
E929
E935
.10.1097/BRS.0b013e3181fd7f7f
56.
Ferguson
,
S. J.
,
Ito
,
K.
, and
Nolte
,
L. P.
,
2004
, “
Fluid Flow and Convective Transport of Solutes Within the Intervertebral Disc
,”
J. Biomech.
,
37
(
2
), pp.
213
221
.10.1016/S0021-9290(03)00250-1
57.
Newman
,
H. R.
,
Bowles
,
R. D.
, and
Buckley
,
M. R.
,
2018
, “
Viscoelastic Heating of Insulated Bovine Intervertebral Disc
,”
JOR Spine
,
1
(
1
), pp.
1
10
.10.1002/jsp2.1002
58.
O'Connell
,
G. D.
,
Vresilovic
,
E. J.
, and
Elliott
,
D. M.
,
2007
, “
Comparison of Animals Used in Disc Research to Human Lumbar Disc Geometry
,”
Spine (Phila. Pa. 1976)
,
32
(
3
), pp.
328
333
.10.1097/01.brs.0000253961.40910.c1
59.
Ristaniemi
,
A.
,
Šećerović
,
A.
,
Dischl
,
V.
,
Crivelli
,
F.
,
Heub
,
S.
,
Ledroit
,
D.
,
Weder
,
G.
,
Grad
,
S.
, and
Ferguson
,
S. J.
,
2023
, “
Physiological and Degenerative Loading of Bovine Intervertebral Disc in a Bioreactor: A Finite Element Study of Complex Motions
,”
J. Mech. Behav. Biomed. Mater
,
143
, p.
105900
.10.1016/j.jmbbm.2023.105900
60.
Farrell
,
M. D.
, and
Riches
,
P. E.
,
2013
, “
On the Poisson's Ratio of the Nucleus Pulposus
,”
ASME J. Biomech. Eng.
,
135
(
10
), p.
104501
.10.1115/1.4025180
61.
Du
,
Y.
,
Tavana
,
S.
,
Rahman
,
T.
,
Baxan
,
N.
,
Hansen
,
U. N.
, and
Newell
,
N.
,
2021
, “
Sensitivity of Intervertebral Disc Finite Element Models to Internal Geometric and Non-Geometric Parameters
,”
Front. Bioeng. Biotechnol.
,
9
(
June
), pp.
1
12
.10.3389/fbioe.2021.660013