Finite element models can reproduce the effect of nucleotomy on the multi-axial compliance of human intervertebral discs

Abstract

Finite element (FE) models can unravel the link between intervertebral disc (IVD) degeneration and its mechanical behaviour. Nucleotomy may provide the data required for model verification. Three human IVDs were scanned with MRI and tested in multiple loading scenarios, prior and post nucleotomy. The resulting data was used to generate, calibrate, and verify the FE models. Nucleotomy increased the experimental range of motion by 26%, a result reproduced by the FE simulation within a 5% error. This work demonstrates the ability of FE models to reproduce the mechanical compliance of human IVDs prior and post nucleotomy.

Keywords:

• Intervertebral disc
• in vitro testing
• finite element modeling
• nucleotomy
• magnetic resonance imaging

1. Introduction

The intervertebral disc (IVD) is an essential part of the spine, providing the required flexibility and distributing the loads on the adjacent vertebrae. The kinematic and mechanical functions of the IVD are determined by its histologically and functionally distinct nucleus pulposus (NP), annulus _brosus (AF), and endplate (EP) tissues. In a healthy disc, the highly hydrated NP supports compressive loads while the surrounding AF supports tensile stresses caused by torsion or bending and constrains the lateral expansion of the NP by resisting tensile hoop stresses (White and Panjabi 1990).

IVD degeneration is an aberrant, cell-mediated response to progressive structural failure (Adams and Roughley 2006). NP and AF show aggrecan degradation, change in collagen type I to type II ratio, decrease in water and glycosaminoglycans (GAG) concentration, reduction in the concentration of viable cells (Johnstone and Bayliss 1995; Bibby et al. 2001) and increased synthesis of matrix-degrading enzymes (Roughley 2004). Eventually, these alterations lead to height loss (Brinckmann 1986; Goel et al. 1986), load redistribution within the disc (Adams et al. 1996) and on the underlying vertebra (Maquer et al. 2015). At advanced stages of degeneration, fissures and cracks appear in the disorganized AF (Vernon-Roberts 1980), leading to posterior disc bulging and ultimately spinal stenosis (Adams and Roughley 2006), which all have been clinically associated with disc herniations or discogenic pain (Brisby 2006; Yang et al. 2015). From a strictly mechanical point of view, disc degeneration generally increases joint laxity and stiffness (Ellingson et al. 2013). Hence, the ability to estimate the mechanical behaviour of an IVD could potentially be linked to its state of degeneration and thus be used as a diagnostic tool for the detection of early degeneration.

In silico models, especially finite element analysis (FEA), are increasingly employed to investigate the effect of such degenerative alterations in both composition and structure of the IVD (Gilbertson et al. 1995; Schmidt et al. 2013). The latest models were recently used to show the impact of degeneration on the NPs swelling capacity (Yang and O’Connell 2019), to determine the bone remodelling mechanisms of the vertebral bodies following a nucleotomy (Calvo-Echenique et al. 2019), or to detect early degenerative changes via magnetic resonance imaging (MRI)-based models (Chetoui et al. 2019).

Numerical models require experimental data for calibration and validation purposes. Unfortunately, inducing reproducible IVD degeneration on human discs in vitro to generate such data is difficult. On the other hand, nucleotomy, a procedure consisting in the partial or complete removal of the NP, is simple and easy to reproduce. Nucleotomy is a radical procedure that drastically changes the internal structure of the IVD. It was shown to cause disc de-pressurization, EP strain reduction (Broc et al. 1997; Frei et al. 2001; Heuer et al. 2008), increase in the disc’s neutral zone under cyclic compressive loads (Showalter et al. 2014), and to affect the deformation patterns of the laminae within the AF (Seroussi et al. 1989; Meakin and Hukins 2000; Meakin et al. 2001). Therefore, nucleotomy constitutes an ideal model to study how the disc’s mechanical behaviour is influenced by structural changes.

To allow the development and implementation of precision, patient-specific, diagnostic tools for detecting early degenerative changes in the IVD, require careful validation of the model’s ability to predict the mechanical response of both the intact and the damaged IVD under a full range of daily spinal loading. This study aims to demonstrate the potential of MRI-based FE models to estimate the behaviour of the intact IVDs under a wide range of applied loading conditions, and its ability to successfully predict the effect of severe damage to its constitute functional tissue, for example, the significant loss of nucleus function introduced via a nucleotomy procedure (Figure 1).

Figure 1. Study overview: Three human lumbar IVDs were prepared (A) and scanned in a clinical 3 T MRI device (B). Using the acquired images, a finite element model was generated. The model is composed of isotropic NP and anisotropic AF elements including to distinct fibre families (C). The IVD joints were embedded in PMMA end-caps, and the mechanical compliance of the disc joints was measured in bending, flexion, extension, rotation (D) and tension (E). The measured load-displacement curves were used to calibrate the material properties of the FE model (F). Finally, a trans-endplate nucleotomy was performed and the testing procedure was repeated.

2. Materials and methods

2.1. Samples

Lumbar vertebral IVD’s are the most affected by disc degeneration (Saleem et al. 2013). Therefore, three lumbar, functional spine units were harvested from two fresh frozen cadaveric spines (A: L2-L3 and L4-L5 segments from a female, 32 years, and B: L2-L3 segments from a male, 60 years). With the sample still frozen, the posterior processes were removed and all soft tissues dissected with care not to injure the IVD. To isolate the disc joint, the cranial and caudal vertebral bodies were sectioned transversely at mid-body height. Each specimen was wrapped in saline-soaked gauze and stored in a double plastic bag at −20°C. Freezing can damage the disc tissue and potentially lead to a change in the its mechanical behaviour (Hongo et al. 2008). We assumed that such damage, if existing, would impact the IVD’s mechanical behaviour before and after nucleotomy similarly.

2.2. MR imaging

Each IVD was thawed, wrapped in saline soaked gauze for 4 hours at room temperature. The samples were then placed in a custom-made imaging chamber. To decrease susceptibility effects at the surface of the tissue (Benveniste et al. 2000) and prevent dehydration, the imaging chamber was filled with Manganese chloride (MnCl2, 10 mM). The chamber was located in a hand/wrist coil (MR2AGD330172, Siemens Healthcare, Germany) with the disc’s cranial-caudal axis aligned with the B0 field of the MRI scanner (Magnetom Prisma fit, 3 T, Siemens Healthcare, Germany). T1-weighted sequences (TR: 14.1 ms, TE: 5 ms) for IVD-background contrast and proton density (PD) weighted sequences (TR: 1700 ms, TE: 5.5 ms) for NP-AF contrast were used to acquire axial images (slice thickness: 0.4 mm, in-plane resolution: 0.285 mm) (Figure 3(A,D)). A spine surgeon evaluated the degeneration state of the discs using the Pfirrmann criterion (Pfirrmann et al. 2001).

2.3. Mechanical testing

Immediately after MRI imaging, the cranial and caudal vertebral bodies were embedded in polymethyl methacrylate (PMMA) cement (Technovit 3040, Heraeus Kulzer, Hanau, Germany). The bone-PMMA interface stability was further increased by inserting screws into the vertebral bodies before embedding (Wilke et al. 2002). A saline-soaked gauze was wrapped around the samples at all times to prevent drying of the sample.

2.3.1. Multi-axial testing

Each embedded disc joint was secured to a custom, computer-controlled, six degrees of freedom (DOF) spine testing system (Gédet et al. 2007). The disc was exposed to pure moments (up to 6 Nm, 0.5°/s, displacement controlled) in flexion, extension, and assuming symmetry, right lateral bending and right axial rotation. Thanks to a friction-free air-bearing system, the disc was allowed to freely deform about the remaining five DOF. Each loading scenario was repeated five times. An active optical tracking system (Optotrak Certus, Northern Digital, Canada, accuracy: 0.1 mm) was used to record the resulting six DOF motions of the disc joint while all forces and torques were recorded using a six DOF load cell (MC3A, AMTI, USA, accuracy: 2 N, 0.01 Nm). The error of the measured angles was shown to be below 0.2°.

2.3.2. Tensile testing

Although the application of tensile loading has limited clinical relevance, it imposes uniform tensile loading on the AF fibres, thus providing relevant data to calibrate or validate numerical models. For this purpose, the disc was secured to a hydraulic testing system (858 Mini Bionix II, MTS, Eden Prairie, MN) and five loading cycles (0 N to 270 N, loading rate = 10 N/s) were applied at the geometrical center of the IVD. The test frame was designed to allow the cranial vertebra to freely translate along the disc’s transversal plane with rotational motions constrained. The applied loads and resulting displacements were recorded using the MTS built-in sensors (load cell: model 662.20 D-04 and LVDT displacement sensor).

2.4. Trans-endplate nucleotomy

After testing, a trans-endplate nucleotomy (Johannessen et al. 2006; Vresilovic et al. 2006) was employed to depressurize and mechanically remove parts of the NP. The centre of the NP with respect to the joint transverse plane was identified on the MRI images, and a 4.5 mm hole was drilled through the EP along the cranial-caudal axis. No further steps to remove NP material other than drilling were performed. This procedure preserved the structural integrity of the AF and ensured exposure of the NP cavity to environmental pressure. However, it did damage the vertebral EP. Once completed, the entire testing procedure was repeated.

2.5. Loading curve analysis

Our focus was on the role of NP and fibre arrangement on the multiaxial mechanical behaviour of the IVD. Therefore, we chose not the include any viscous effects in the present model. Accordingly, the acquired in vitro loading curves had to be post-processed to serve as ground truth for quasi-static models.

The measured data points (acquisition rate: 100 Hz) were cleaned with a moving average filter (kernel size: 5). The first three loading cycles of each loading scenario were discarded as preconditioning. Cycle four and five were averaged and the mean value of the hysteresis was computed at 20 locations, resulting in an average interval of 0.3° (0.03 mm). For the tensile loading case, the drift caused by the discs poro- and viscoelastic behaviour was subtracted before averaging the cycles. The resulting curves were used to calibrate the material constants. A third-order polynomial was fitted through the points (R² = 0.98) from which the minimum stiffness (S_min) at 0.1° (0.02 mm) and maximum stiffness (S_max) at the maximum displacement were computed (Figure 2). The ROM was defined as the displacement at the maximum moment (force).

Figure 2. Definition of the flexibility parameters used to compare the in vitro measurements with the FE predictions.

2.6. Image processing and model generation

The entire image processing, meshing and material mapping pipeline was developed using Matlab (V.15, Mathworks, USA). First, a multi-threshold Otsu algorithm (Otsu 1979) was applied to the T1-weighted MR images (Figure 3(A)) to extract the geometry of the IVD. For each image dataset, six thresholds were computed (Figure 3(B)). Regions corresponding to background or noise were removed, the others were merged into the IVD mask (Figure 3(C)). Morphological operations (erode and dilate) were applied to remove the islands, close the holes and smooth the mask. In a second step, the previously generated mask was used to generate quadratic, 10-node tetrahedral elements using the ISOMESH toolbox (Fang and Boas 2009). The element size was selected based on our previous study where the effect of mesh density on quality metrics, solution accuracy and CPU time were evaluated (Maquer et al. 2015). As all tested meshes were of good quality and within 5% difference for all loading cases (flexion, extension, torsion, bending and compression), the coarsest mesh (3 mm edge length) was used. Next, the segmentation segmentation procedure was repeated for the PD-weighted images (Figure 3(D,E)) to create the mask of the NP. This time, only the innermost region of the Otsu regions was kept (Figure 3(F)) before applying the morphological operations. Finally, all elements from the IVD mesh whose centroid was in the NP mask were assigned to the NP element set, the others to the AF element set.

Figure 3. The geometry of the IVD was defined from the T1-weighted MR image (A). An multi-region Otsu threshold was first applied (B). The non-connected voxels and holes were then removed with morphological operations to generate the IVD mask (C). The same procedure was applied to the proton density-weighted MR image (D, E) to generate the NP mask (F). The IVD mask was used to create a tetrahedral mesh (G). The elements belonging to the NP mask were assigned the NP material properties, the others were assigned the AF material properties, including a fibre orientation based on the outer shape of the AF (H).

The boundary conditions were defined to replicate experimental test conditions as follows; the caudal and cranial nodes were combined in two kinematic coupled sets corresponding to the vertebral end-plates. The caudal element set was constrained in all directions while the cranial set was assigned the moment/force of the individual loading cases. The cranial set was allowed to move freely along the remaining DOF.

2.7. Constitutive law and material mapping

Inspired by the Holzapfel-Gasser-Odgen model (Holzapfel et al. 2000; Eberlein et al. 2001), a time-independent hyper-elastic, fibre-reinforced material law was implemented as an Abaqus (6.13 Dassault Systemes, France) user material. The law was described as a free energy function (Equation (1)(1) $$\psi ={\psi }_{\mathit{GroundSubstance}}+{\psi }_{\mathit{fibres}}$$(1) ). It incorporates an isotropic term representing the ground substance (Equation (2)(2) $${\psi }_{\mathit{GroundSubstance}}(\mathbf{C})=\frac{\kappa }{2}{(J-1)}^2+\frac{\mu }{2}(J^{-\frac{2}{3}}tr(\mathbf{C})-3)$$(2) ), and an anisotropic term with two normalised vectors (a4 and a6 corresponding to the two main collagen fibre directions) representing the fibrous part of the IVD (Equation (3)(3) $${\psi }_{fibres}(C,a_4,a_6)=\frac{k_1}{2}\left[{(I_4-1)}^{2k_2}+{(I_6-1)}^{2k_2}\right]$$(3) ). The angles of these vectors in the transverse plane of the IVD were defined following the outer shape of the AF. The angle of the vectors with respect to the transverse plane was defined during the material constants optimisation procedure. Elements from the NP element set were assigned isotropic material properties (ψ_fibres=0). Those from the AF element set were assigned both the isotropic and anisotropic part of the material law. (1) $$\psi ={\psi }_{\mathit{GroundSubstance}}+{\psi }_{\mathit{fibres}}$$(1) (2) $${\psi }_{\mathit{GroundSubstance}}(\mathbf{C})=\frac{\kappa }{2}{(J-1)}^2+\frac{\mu }{2}(J^{-\frac{2}{3}}tr(\mathbf{C})-3)$$(2) (3) $${\psi }_{fibres}(C,a_4,a_6)=\frac{k_1}{2}\left[{(I_4-1)}^{2k_2}+{(I_6-1)}^{2k_2}\right]$$(3) (4) $$C=F^{T}F$$(4) (5) $$J=\mathit{det}{(\mathbf{C})}^{\frac{1}{2}}$$(5) (6) $$I_{4,6}(\mathbf{C},{\mathbf{a}}_{4,6})={\mathbf{a}}_{4,6}·\mathbf{C} {\mathbf{a}}_{4,6}$$(6) k1 and k2 correspond to the stiffness and non-linearity of the AF fibres with the invariants I4 and I6 describing the deformations and reorientations of the fibres derived from the deformation gradient F. At the material level, the fibres are only active under tension ${\psi }_{fibre{s}_{4,6}}$ = 0 if I_4,6 < 1). κ is the resistance against volumetric deformations (bulk modulus), μ the resistance against isochoric deformations (shear modulus). As highlighted by Nolan and colleagues (Nolan et al. 2014), the full anisotropic invariants (Equation (6)(6) $$I_{4,6}(\mathbf{C},{\mathbf{a}}_{4,6})={\mathbf{a}}_{4,6}·\mathbf{C} {\mathbf{a}}_{4,6}$$(6) ) were used.

2.8. Identification of the material parameters against the experimental tests

Based on the protocol of Maquer et al. (Maquer et al. 2015), a particle swarm optimisation (PSO) (w: 0.7, c1: 1.2, c2: 1.6, nb particles: 25) algorithm (Clerc and Kennedy 2002) was used to calibrate the material constants by minimizing the difference between the simulated and experimental data curves. The custom algorithm was implemented in Python (Python v2.7.13, Numpy v1.13).

The outcome of this optimisation provided the numerical values for the NP intact shear modulus (μ_np_i), the NP intact bulk modulus (κ_np_i), the AF shear modulus (μ_af), the AF bulk modulus (κ_af), the fibre stiffness (k1), the fibre non-linearity (k2) and the fiber angle (ϕ). With the assumption that the nucleotomy only affected the elements of the NP, the damaged NP shear modulus μ_np_d and the damaged NP bulk modulus κ_np_d were determined upon secondary optimisation procedure against the in vitro data after nucleotomy where all material constants belonging to the AF were kept from the intact calibration round. In total, 18 hours were required for the complete set of optimisation (Intel i7 quad-core, 3.2 GHz, 32 GB RAM). As previously described for the measurements, a third-order polynomial was fitted through the predicted load-displacement curves and S_min, S_max, and ROM were computed to allow for comparison between in vitro and in silico data.

3 Results

3.1. Comparison intact and denucleated load-displacement data

The in vitro load-displacement curves for both intact and denucleated IVDs are presented in Figure 4. In response to moments of 6 Nm, nucleotomy resulted in an average increase of 20% for flexion, extension and lateral bending, 28% for rotation, and 43% for tension. On average, minimum stiffness decreased by by 52%. Although no consistent direction of change was not observed for the maximal stiffness, it changed on average by 6% under all moment loading modes (Figure 5).

Figure 4. The FE model replicates the in vitro load-displacement curves. Black and gray represent the experimental data for the intact and damaged case, the dotted lines show the corresponding FE predictions.

Figure 5. Mean relative difference in ROM, minimal stiffness and maximal stiffness between intact and denucleated samples at 6 Nm (270 N) for all loading cases for the experimental data (black) and FE predictions (blue).

All IVDs were graded according to the Pfirrmann scale (Pfirrmann et al. 2001) as grade two. The older sample (B_L2L3) demonstrated decreased ROM and increased initial stiffness in all directions as well as reduced post nucleotomy alterations compared to the younger samples (A_L2L3 and A_L4L5).

3.2. Comparison experimental and simulated load-displacement data

Figure 4 summarizes all in silico load-displacement curves for intact and denucleated IVDs. The average discrepancies of the intact FE model were 4% for the ROM, 31% for the S_min and 14% for S_max compared to the experimental values. For the denucleated FE model, these were 5% for the ROM, 35% for S_min and 13% for S_max against the in vitro data. A detailed table of all computed flexibility parameters is given in Table A1 (appendix A). Table 1 summarizes all material constants for the intact and damaged case computed by the PSO calibration.

Table A1. S_max, S_min and ROM at 6 Nm (270 N) for intact and damaged discs, for all loading cases, experimental data and FE predictions.

		A_L2L3				A_L4L5				B_L2L3
		intact		damage		intact		damage		intact		damage
loading case		Exp	FE	EXP	FE	EXP	FE	Exp	FE	Exp	FE	Exp	FE
Flex.	ROM [°]	5.75	5.68	7.29	7.40	6.50	6.88	7.59	8.78	4.77	4.92	5.45	5.56
	S_max [Nm/°]	2.50	2.76	2.30	2.64	1.83	1.90	2.10	2.07	2.28	2.21	2.45	2.18
	S_min [Nm/°]	0.29	0.15	0.10	0.11	0.35	0.28	0.14	0.10	0.87	0.60	0.44	0.44
Ext.	ROM [°]	8.08	8.12	10.15	10.22	9.25	9.09	11.09	10.97	5.38	5.56	6.09	6.47
	S_max [Nm/°]	1.56	1.61	1.70	1.68	1.38	1.43	1.41	1.44	2.00	1.86	1.68	1.76
	S_min [Nm/°]	0.28	0.26	0.09	0.14	0.32	0.23	0.15	0.05	0.79	0.54	0.47	0.38
Bend.	ROM [°]	4.84	4.81	6.35	6.28	4.92	4.85	6.05	5.86	4.40	3.95	4.85	4.47
	S_max [Nm/°]	2.62	3.04	2.95	3.07	2.60	2.89	2.69	3.07	2.63	3.82	2.73	3.80
	S_min [Nm/°]	0.45	0.36	0.16	0.14	0.43	0.46	0.21	0.13	0.61	0.39	0.41	0.16
Rot.	ROM [°]	2.71	2.59	3.83	3.94	2.86	2.82	3.60	3.42	2.37	2.11	2.78	2.44
	S_max [Nm/°]	3.70	4.39	3.65	3.15	3.20	3.77	3.33	3.63	3.58	5.53	3.20	5.23
	S_min [Nm/°]	1.42	0.36	0.38	0.29	1.35	0.36	0.70	0.18	1.78	0.44	1.40	0.16
Tens.	ROM [mm]	0.43	0.44	0.67	0.67	0.52	0.48	0.78	0.71	0.42	0.40	0.53	0.53
	S_max [kN/mm]	930.8	842.2	815.6	783.9	844.8	800.7	706.4	771.6	681.3	679.4	595.7	552.6
	S_min [kN/mm]	362.6	332.5	97.9	102.6	292.5	279.9	138.2	72.8	594.7	572.7	386.3	409.9

Table 1. Material constants resulting from the calibration of the models with the in vitro load-displacement data. The index _i stands for intact and _d for the denucleated variables with μ_np the NP shear modulus, κ_np the NP bulk modulus, μ_af the AF shear modulus, κ_af the AF bulk modulus, k1 the fibre stiffness, k2 the fibre non-linearity, and ϕ fibre angle.

	A_L2L3	A_L4L5	B_L2L3
μ_np_i [MPa]	1.0	0.9	1.69
μ_np_d [MPa]	0.04	0.04	1.4
μ_af [MPa]	0.07	0.01	0.01
κ_np_i [MPa]	684.4	699.9	880.8
κ_np_d [MPa]	0.01	267.2	31.9
κ_af [MPa]	620.3	519.0	358.6
k1 [MPa]	93.2	99.8	108.9
k2 [-]	1.45	1.43	1.57
ϕ [°]	25.2	31.9	32.0

3.3. Element size sensitivity

The sensitivity of the FE model to the element size was evaluated for the intact case only. Reducing the element size from 3 mm mean edge length to 2.5 mm resulted in an average increase of 3.7% in predicted ROM, a decrease of 3.5% in S_min, and a decrease of 3.8% in S_max. Increasing the element size to 3.5 mm resulted in a decrease of 5.9% in ROM, an increase of 3.0% in S_min and an increase of 5.2% in S_max (Table 2).

Table 2. Element size sensitivity. Native mean edge length is 3 mm, (−) corresponds to 2.5 mm mean edge length and (+) 3.5 mm mean edge length. All values are given in [%] change compared to the native element size.

		ROM	S_min	S_max
Flex	−	3.0	0.5	−7.4
	+	−5.3	5.2	4.5
Ext.	−	3.9	−6.9	−5.5
	+	−6.6	1.6	8.1
Bend.	−	7.6	−6.8	−2.9
	+	−4.7	−1.9	5.6
Rot.	−	3.4	−2.4	−2.8
	+	−1.7	2.0	1.6
Tens.	−	0.5	−0.7	−0.3
	+	−11	4.4	6.2

4. Discussion

A nucleotomy procedure was applied as a simple model to induce clinically relevant structural changes into human IVDs. The goal of the study was to develop a subject-specific IVD FE model with the ability mimic the alterations of the IVD’s overall mechanical behaviour. For this purpose, the multi-axial compliance of three human IVDs was measured prior and post nucleotomy in the laboratory. The material constants of the model were calibrated against the in vitro data and the model’s ability the reproduce the measured load-displacement curves was evaluated.

In vitro nucleotomy resulted in an ROM increase of average increase in ROM of 26% at 6 Nm (270 N). The minimal stiffness was reduced by 52% while the maximal stiffness did not follow a clear trend, but changed on average by 6%. These findings are following previous studies that observed a clear trend of increased ROM and a reduction in initial stiffness post nucleotomy. Applying a nucleotomy on a L4-L5 IVD, Wilke et al. reported at 7.5 Nm an average increase in ROM of 18% for flexion/extension, 12% for lateral bending and 21% for axial rotation (Wilke et al. 2002). The lower values by the Wilke study can be explained with the fact that, unlike our study, the facet joints were left in place. While successively reducing the anatomy of spinal segments, Heuer and colleagues found that the nucleotomy increased the ROM (at 5 Nm) by 37% for flexion, 49% for extension, 38% for bending and 47% for rotation (Heuer et al. 2008). Lafage et al. (2007) reported an increase in ROM (at 8 Nm) after nucleotomy of 15% for flexion and extension, 39% for lateral bending and 50% for axial rotation. While the nucletomy procedure was only performed on one disc, the testing procedure was applied to two complete motion segments, making direct comparisons more difficult. It’s also important to note that the nucleotomies in these studies were conducted via an incision through the AF, which might have further changed the mechanical response of the IVD (Johnstone and Bayliss 1995).

The fibre angles yielded by the optimisation procedure (25.2° to 32°) were in agreement with anatomical observations (Eberlein et al. 2004; Iatridis and Ap Gwynn 2004; Holzapfel et al. 2005; Pezowicz 2010). However, previous studies further found that the angle of the fibres varies from 23° in the anterior part to 45° in the posterior part and that the type of fibres and their associated stiffness change from the outer to the inner AF (Holzapfel et al. 2005). Although we did not account for these regional variations of the fibrous network, the FE model successfully reproduced the moment-deformation response curves measured experimentally for the intact and denucleated IVD. Our model estimated the measured ROM with an average error below 5%, S_min was estimated with an error of 33% and S_max with an error 14%. We suggest that the larger error in estimation for S_min was mainly caused by the low absolute values.

A sensitivity analysis revealed that altering the mean element edge length by ± 20% resulted on average in less than 5% change of the flexibility parameters (ROM, S_min and S_max). These findings are in agreement with our previous study (Maquer et al. 2015). These variations have two causes beyond the simple effect of mesh density for isotropic homogeneous materials. First, because of the homogenisation, the coarser the mesh, the “rougher” the in-plane fibre distribution. Second, the IVD was meshed as a single part and only later assigned NP or AF material constants. While this method allows for a simple, automatic and regular mesh generation, a change in element size can lead to a slightly different ratio between NP and AF volume (± 2%) and different geometry. We have employed this method for the sake of limiting the computing time associated with our optimisation strategy. A more robust method would be to apply a direct MRI-based material mapping as we previously suggested (Stadelmann et al. 2018). Still, based on the sensitivity tests, these variations were minimal.

Beyond the common constraints associated with the use of cadaver material in vitro, this study is limited by the fact that our FE model was not designed to simulate any time-dependent behaviour, which is present in the IVD (Humzah and Soames 1988; Keller et al. 1987). The focus of this feasibility study was to develop a model that can accurately estimate a multi-axial quasi-static response.

Our assumption of a fully reversible elastic behaviour is supported by a previous study demonstrating the sensitivity of angular motion to loading rate to be relatively low (Wilke et al. 1998). In our study, the hystereses of the acquired angular load-displacement curves accounted on average for 10.5% of the stored elastic energy and was found not to drift between the cycles. For the tensile loading scenario, we measured a small drift that was subtracted from the data to allow a comparison with a quasi-static model. Such a correction would not have been possible for the highly rate-dependent axial compression (Race et al. 2000). To avoid misinterpretations, no measurements and no simulations were performed in axial compression.

Even with the low number of samples, the experimental data revealed the large intra-specimen variability (up to 45% difference in ROM). These variations are most probably caused by differences in both the geometry (Maquer et al. 2014) and composition of the individual IVDs. While the geometry can be captured very well using clinical MRI scanners, the in vivo quantification of the IVD composition and fibre architecture remain a major challenge. For the present study, a specific set of material constants was calibrated for each IVD, limiting its use for potential clinical applications. Our findings call therefore for improved MR imaging to quantify collagen fibre density/orientation and proteoglycan content in order to directly assign IVD-specific material properties (Stadelmann et al. 2018). In the future, it might be possible to augment in vivo low resolution clinical MR images with a-priori information extracted from ex vivo high resolution images using machine learning techniques as it is already done with bone tissue and CT images (Chandran et al. 2017, 2019).

In conclusion, this study demonstrated that a rate-independent FE model of the IVD is able to reproduce up to five distinct load cases with a fair accuracy before and after nucleotomy. Nucleotomy is an extreme case of structural damage. FEA has therefore the potential to cope with any degenerative alterations of the disc structure, especially if combined with quantitative MR imaging.

Disclosure statement

The authors have no conflict of interest to declare.

Additional information

Funding

The research reported in this publication was supported by the Swiss National Science Foundation, grant number 147153 and the National Institutes of Health under award number R21AR066916.

Appendix A

Appendix B

Acronyms

AF	=	annulus fibrosus
DOF	=	degrees of freedom
EP	=	endplate
FEA	=	finite element analysis
ϕ	=	fibre angle
IVD	=	intervertebral disc
k1	=	fibre stiffness
k2	=	fibre non-linearity
κ_af	=	AF bulk modulus
κ_np_d	=	NP damaged bulk modulus
κ_np_i	=	NP intact bulk modulus
μ_af	=	AF shear modulus
μ_np_d	=	NP damaged shear modulus
μ_np_i	=	NP intact shear modulus
NP	=	nucleus pulposus
PD	=	proton density
PMMA	=	polymethyl methacrylate
PSO	=	particle swarm optimisation
S_max	=	maximum stiffness
S_min	=	minimum stiffness
TE	=	echo time
TR	=	repetition time