A Novel Fiber-Reinforced Poroviscoelastic Bovine Intervertebral Disc Finite Element Model for Organ Culture Experiment Simulations

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 [2–10]. 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 [4–7,10,14–19]. 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 [21–24], hyperelastic material [25–27], as well as fiber-reinforced poroviscoelastic with swelling [28,29]. The AF material models have, likewise, ranged from fiber-reinforced hyperelastic [22–26,30] and poroelastic [31] to fiber-reinforced poroviscoelastic with swelling [27–29,32]. The material property distribution within AF is commonly taken into account [22–26], 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.

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.

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.

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

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 [41–43] and the accurateness of fit was estimated based on the R²-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⁻¹⁵ m⁴/Ns [44], and for the NP they were fixed to $ENP$= 0.019 MPa [44], $νNP$= 0.49 [21] and $kNP$=32 × 10⁻¹⁵ m⁴/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.

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.

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

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.

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.

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 [23–26] 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.

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

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.

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

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

$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