Sensitivity of Intervertebral Disc Finite Element Models to Internal Geometric and Non-geometric Parameters

Finite element models are useful for investigating internal intervertebral disc (IVD) behaviours without using disruptive experimental techniques. Simplified geometries are commonly used to reduce computational time or because internal geometries cannot be acquired from CT scans. This study aimed to (1) investigate the effect of altered geometries both at endplates and the nucleus-anulus boundary on model response, and (2) to investigate model sensitivity to material and geometric inputs, and different modelling approaches (graduated or consistent fibre bundle angles and glued or cohesive inter-lamellar contact). Six models were developed from 9.4 T MRIs of bovine IVDs. Models had two variations of endplate geometry (a simple curved profile from the centre of the disc to the periphery, and precise geometry segmented from MRIs), and three variations of NP-AF boundary (linear, curved, and segmented). Models were subjected to axial compressive loading (to 0.86 mm at a strain rate of 0.1/s) and the effect on stiffness and strain distributions, and the sensitivity to modelling approaches was investigated. The model with the most complex geometry (segmented endplates, curved NP-AF boundary) was 3.1 times stiffer than the model with the simplest geometry (curved endplates, linear NP-AF boundary), although this difference may be exaggerated since segmenting the endplates in the complex geometry models resulted in a shorter average disc height. Peak strains were close to the endplates at locations of high curvature in the segmented endplate models which were not captured in the curved endplate models. Differences were also seen in sensitivity to material properties, graduated fibre angles, cohesive rather than glued inter-lamellar contact, and NP:AF ratios. These results show that FE modellers must take care to ensure geometries are realistic so that load is distributed and passes through IVDs accurately.

Introduction

Intervertebral discs (IVDs) lie between vertebra in the spine and act to distribute loading while allowing the spine to bend and flex (Bogduk, 2005; Adams and Roughley, 2006). The IVD consists of the anulus fibrosus (AF), nucleus pulposus (NP) and the endplates that enclose the NP and AF above and below. As the IVD degenerates, its height reduces (Adams et al., 1987), the relative ratio of NP to AF cross-sectional area decreases (Adams et al., 1996), and the stiffness of individual components increases (Iatridis et al., 1998, 1999; O’Connell et al., 2009). Previously, these changes have been investigated using experimental techniques, however, it is challenging to measure stresses and strains within the IVD without disrupting it, and therefore, there has been an increasing trend towards the use of FE models for investigations of this kind.

Previous FE studies on human lumbar spines have explored the effect of altering geometric features of the IVD (Robin et al., 1994; Natarajan and Andersson, 1999; Noailly et al., 2007; Meijer et al., 2011; Niemeyer et al., 2012). All of these studies have found the disc height significantly affects the response of the IVD, and some have shown that NP position (Noailly et al., 2007), and endplate width and depth (Natarajan and Andersson, 1999; Meijer et al., 2011; Niemeyer et al., 2012) are important factors. Similarly, the nature of the contact between lamella in the AF in FE models has been shown to affect model response (Adam et al., 2015), and a method of using cohesive contact between lamella has been proposed by Mengoni et al. (2015). Recent FE studies on bovine (Mengoni et al., 2017), and human discs (Yang et al., 2019), have demonstrated the importance of the NP:AF ratio, few studies have investigated the effect of changes to the internal IVD geometry. This is likely due to the challenges associated with defining the boundary between the NP and AF. Recent advances in magnetic resonance imaging (MRI) allow accurate internal geometries to be identified (O’Connell et al., 2011; Tavana et al., 2020, 2021). In some cases, these geometries are substantially different to those that are used in current FE models, and the effect of these inaccuracies has not previously been investigated.

The aim of this study is to use high resolution MRI scans, and FE models to investigate the effect of altered geometries both at the boundary between the NP and the AF, and at the endplates on the response of the IVD. Specific aims include quantifying the effect of altered geometries on IVD;

(a) stiffness and strain distributions

(b) sensitivity to altered material properties

(c) sensitivity to graduated fibre bundle angles (increasing from outer to inner AF)

(d) sensitivity to modelling inter-lamellar behaviour with cohesive contact

(e) sensitivity to altered NP:AF ratio.

Materials and Methods

For this study, a vertebral body—disc—vertebral body specimen was dissected from the most caudal disc of a fresh-frozen bovine tail acquired from a local butcher. Soft tissue was removed before being scanned on a 9.4 T MRI scanner (Bruker BioSpec, Ettlingen, Germany) equipped with a volume RF resonator [T₂-weighted RARE sequence, coronal plane, resolution = (90 × 90) μm², slice thickness = 800 μm, 17 min scan time], such that internal geometries could be acquired.

Experimental Data

Experimental data, against which the FE models could be compared was obtained from previous literature (Newell et al., 2017). During these experiments ten bovine IVDs were axially compressed to 15% strain at a range of strain rates and the force-displacement response of each sample was recorded. For comparison purposes only data from mid-range strain rates (0.1/s) were used in this study.

FE Model Development

Non-linear, implicit, axisymmetric FE models using Marc (v2017, MSC Software, California, United States) were developed based on measurements from the mid-coronal slice of the high-resolution MRI images (Figure 1A).

FIGURE 1

Figure 1. (A) Middle coronal slice of the MRI of the bovine IVD. (B) Typical peripheral MRI slice from which fibre bundle angles were measured.

Geometry

Six FE models were developed, with two variations of endplate geometry (simple and segmented), and three variations of NP-AF boundary (linear, curved, and segmented) (Figure 2). Since the bovine tail disc is almost perfectly round all FE models were axisymmetric (Adam et al., 2015). The simple endplate geometries (Model 1.1, 1.2, and 1.3) had central and peripheral IVD heights measured from the mid-coronal MRI slice but a smooth curve between these two points. The segmented endplate geometries (Model 2.1, 2.2, and 2.3) were obtained from the mid-coronal MRI slice using Mimics (v16, Materialise, Leuven, Belgium) but also had central and peripheral IVD heights that matched those of the simple endplate geometry models. The linear NP-AF boundary models had a vertical, linear boundary between the NP and AF, while the curved NP-AF boundary models had a quadratically polynomial boundary. Both the linear, and curved NP-AF boundary models had a NP:AF ratio determined from measurements at the mid-height of the mid-coronal slice of the IVD. As with the endplates, the segmented NP-AF boundary was obtained from the mid-coronal MRI slice using Mimics.

FIGURE 2

Figure 2. Initial geometry of the six axisymmetric FE models overlaid on the MRI slice from which the geometries were obtained. Note the axis of symmetry in the models is along the mid-sagittal plane of the IVD. (A) Model 1.1—simple endplates, linear NP-AF boundary, (B) Model 1.2—simple endplates, curved NP-AF boundary, (C) Model 1.3—simple endplates, segmented NP-AF boundary, (D) Model 2.1—segmented endplates, linear NP-AF boundary, (E) Model 2.2—segmented endplates, curved NP-AF boundary, (F) Model 2.3—segmented endplates, segmented NP-AF boundary. Note the fibre bundle angles are relative to the transverse plane.

The AF was modelled with rebar elements to represent the collagen fibre bundles, surrounded by non-linear solid quadrilateral elements to represent the AF ground matrix. For the models with segmented NP-AF boundary, each lamella layer was segmented, while the number of rebar layers for the linear and curved models was set to the number of layers that could be identified on the mid-coronal slice of the MRIs. The NP was modelled with non-linear solid triangular elements. For this study endplates and vertebral bodies were assumed to be rigid as the effect of the endplates, and VBs on the behaviour of bovine IVD FE models has been found to be negligible (<2.5% difference in the peak force; Newell et al., 2017). A convergence study was performed on each specimen to ensure that the mesh density was sufficient.

Material Properties

The properties of each material used in the models are shown in Table 1. The AF matrix and the NP were assigned non-linear hyperelastic material properties (Mooney-Rivlin). The strain energy function for this material model is shown in Eq. 1, where W is the strain-energy density function, I₁ and I₂ are strain invariants, and C₁₀ and C₀₁ are material constants (Mooney, 1940):

TABLE 1

Table 1. Baseline material properties for components of each of the FE models.

$$W=C_{10}\left(I_1-3\right)+C_{01}(I_2-3)(1)$$

The AF fibre bundles were modelled using tension only rebar elements with a Young’s modulus (YM) of 415 MPa which was obtained from the results of an optimisation study (Newell et al., 2017). The fibre bundles were aligned at ±30.45° to the transverse plane which was an average of six measurements taken from different coronal MRI slices at regular intervals from the inner to outer AF (Figure 1B). The cross-sectional area of each fibre bundle was set to be 3.212 × 10^–2 mm², and spacing of the bundles was set to 4.35 bundles/mm (Marchand and Ahmed, 1990; Adam et al., 2015). The bulk modulus of the NP and AF ground matrix was set at 2,000 MPa to ensure near incompressibility.

Boundary Conditions

Replicating the experimental setup described by Newell et al. (2017), the inferior boundary of the IVD was fixed, and a displacement of 0.86 mm, which corresponded to a central disc axial strain of 15%, was applied to the superior boundary of the IVD. Since the models with segmented endplates had a shorter overall (or average) disc height (eventhough the central disc heights were kept the same), and the displacement of 0.86 mm was applied to all models, the models with segmented endplates were subjected to higher overall (or average) axial strains. Since the model was axisymmetric, nodes along the axis of symmetry were fixed in the radial direction.

Sensitivity Study

The influence of the AF C₁₀, AF C₀₁, fibre bundle YM, and fibre bundle angle values were investigated by varying their baseline value by ±20% and observing the effect on the peak force.

Graduated Fibre Bundle Angle

A number of previous FE studies have modelled fibre bundles with a constant orientation from outer to inner AF (Shirazi-Adl et al., 1986; Marchand and Ahmed, 1990; Natali and Meroi, 1990; Adam et al., 2015; Newell et al., 2017). However, in human IVDs fibre bundle angle has been reported to vary linearly from 28(outer) to 45(inner) in the radial direction relative to the transverse plane (Cassidy et al., 1989), and has been included in a number of other FE studies (Ayturk and Puttlitz, 2011; Schmidt et al., 2012). From the MRI scans, fibre bundle angles of 27.4° and 37.4° were measured at the outer most layer and inner most layer, respectively. In order to understand the effect of including a variation in fibre bundle angle, each layer of fibre bundles were assigned material properties individually in the model. The outer most layer was assigned 27.4° and the inner most layer was assigned 37.4°, while intermediate angles were varied linearly between these two angles.

Modelling Inter-lamellar Behaviour With Cohesive Contact

Interactions between lamellar in the AF were modelled using cohesive elements which allows the interfaces to be described by traction-separation laws. Normal cohesive stiffness (Knn = 1.18 MPa/mm), and tangential cohesive stiffness (Ktt and Kss = 1.31 MPa/mm) values were taken from Mengoni et al. (2015) who derived values using an optimisation algorithm and simulations of tension experiments where ovine AF samples were loaded radially. A stiffening factor in compression was assigned to ensure penetration between elements in adjacent lamellae was minimal. Preliminary investigations showed that a factor of 100,000 was sufficient to ensure that ±20% change resulted in a less than 2% change in peak force.

NP:AF Ratio

The NP:AF ratio was doubled (increasing the NP radius compared to the AF width but keeping the overall IVD width constant) in each of the models to investigate the effect of NP:AF ratio on the mechanics of the disc. Doubling the ratio ensured that there was a clear divergence from the baseline geometry while keeping within the physiological bounds of reported NP:AF ratios (O’Connell et al., 2007; Adam et al., 2015; Newell et al., 2017). In all baseline models the number of rebars/mm was defined based on measurements from Marchand and Ahmed (1990) (4.35 bundles/mm = 0.22 mm interbundle spacing) (Figure 3A). When increasing the NP:AF ratio care was taken to ensure that the total number of fibre bundles, was the same between the baseline, and adjusted NP:AF ratio models (Figure 3). This was achieved by calculating the difference between the circumference of each lamella of the baseline and the altered NP:AF ratio models using the horizontal distance from each lamella to the axisymmetric axis. The fibre bundle spacing was then adjusted to ensure each new model had the same fibre volume as the baseline models (Figure 3).

FIGURE 3

Figure 3. Schematic of axial cross-section of the IVD demonstrating the change in fibre bundle spacing required to allow a fair comparison between models when adjusting the NP:AF ratio. (A) Shows a baseline NP:AF ratio with a cut out showing the fibre bundle thickness (Tb), and interbundle spacing (Sb), (B) shows an increased NP:AF ratio with larger interbundle spacing. Bundles are only shown in the inner most lamella in both schematics, with a total of 12 bundles shown in both (A,B), although larger interbundle spacing in (B) resulting in the same fibre bundle volume as that shown in (A). Note these schematics are not to scale.

Results

Taking measurements from the MRIs, the sample had a central disc height of 8.41 mm, sagittal plane width of 25.56 mm, a coronal plane width of 24.36 mm, and an area of 489.30 mm². A convergence study was performed on each specimen to ensure that the mesh density was sufficient. This involved subdividing the number of elements and comparing the peak force obtained using the original mesh with that obtained with the subdivided mesh. If peak forces were within 1% the original mesh was considered converged, otherwise further subdivisions were performed until consistent (within 1%) peak forces were found. This resulted in models having an average of 1,163 ± 469 elements.

Stiffness and Stress Distributions

The responses obtained from all the FE models were stiffer than the average experimental response reported by Newell et al. (2017; Figure 4). A trend was seen for the stiffness to increase as the complexity of the model geometry increased. All models with segmented endplates had a greater stiffness than those with simple endplates.

FIGURE 4

Figure 4. Force-displacement response of each of the six models compared against average experimental response from Newell et al. (2017) The grey shaded region on the experimental data represents ± 1 standard deviation. Model 1.1—simple endplates, linear NP-AF boundary, Model 1.2—simple endplates, curved NP-AF boundary, Model 1.3—simple endplates, segmented NP-AF boundary, Model 2.1—segmented endplates, linear NP-AF boundary, Model 2.2—segmented endplates, curved NP-AF boundary, Model 2.3—segmented endplates, segmented NP-AF boundary.

In all models the predominant direction of the minimum principal strain at maximum displacement was axial (Figure 5). Axial strains were compressive throughout the IVD in all models, and particularly high along the endplates. Peak axial strains were lowest in the models with segmented endplates, in comparison to the simple endplates (-0.22, -0.34, -0.42, -0.51, -0.64, and -1.00, for Models 1.1, 1.2, 1.3, 2.1, 2.2, and 2.3, respectively). Peak axial and radial strains were seen at locations of high endplate curvature, particularly at the mid-AF-endplate boundary of the segmented endplate models. High compressive axial strains were generally seen at the NP-AF boundary at mid-height in all models. In all models a band of high circumferential strains was seen close to the mid-height of the disc but this band veered away from mid-height in the models with segmented internal geometry (Models 1.3 and 2.3).

FIGURE 5

Figure 5. Minimum principal, axial, radial, and circumferential strain distributions in each of the six models at maximum displacement. Vectors in the left column indicate minimum principal strain directions and the colours refer to magnitudes of the minimum principal strains.

Sensitivity to Altered Material Properties

The models were most sensitive to properties of the fibre bundles, particularly to fibre bundle angle (Figure 6D). The AF C₀₁ Mooney constant had a greater effect than the AF C₁₀ Mooney constant in all models (Figures 6A,B). Models with segmented endplates (2.1, 2.2, and 2.3), on average were less sensitive to changes in material properties in comparison to models with a simple endplate (Models 1.1, 1.2, and 1.3). Models 2.2 and 2.3 (segmented endplates with curved and segmented internal geometry, respectively) had similar sensitivities to all four material parameters, and both were less sensitive to the AF Mooney constants, and the fibre bundle angle, but more sensitive to the fibre bundle YM than Model 2.1 (segmented endplate, linear NP-AF boundary). Similarly, Models 1.2 and 1.3 (simple endplates with curved and segmented internal geometry, respectively) were less sensitive to the AF Mooney constants, and the fibre bundle angle, but more sensitive to the fibre bundle YM than Model 1.1 (simple endplate, linear NP-AF boundary).

FIGURE 6

Figure 6. Sensitivity study overview. The y-axis represents a percentage change in peak axial force from that obtained in the initial baseline run of each of the six models when a parameter was changed by ±20%. (A,B) are AF ground matrix material properties, (C,D) are AF fibre properties.

Effect of Graduated Fibre Bundle Angle

A comparison of the percentage change in peak force between the baseline and varied fibre bundle angle models is shown in Figure 7A. On average, varying the fibre bundle angle resulted in a 5.0 ± 2.8% decrease in peak force in all six models in comparison to the constant fibre bundle angle models. However, Models 2.2 and 2.3 (segmented endplates with curved and segmented internal geometry, respectively) were relatively insensitive to the variation in fibre bundle angle with a percentage change of 0.9 and 2.2%, respectively, compared to 5.9–11.6% for the other models.

FIGURE 7

Figure 7. Percentage changes in peak force compared to baseline models after (A) graduating the fibre bundle angles, (B) cohesive inter-lamellar contact, and (C) increasing the NP:AF ratio.

Sensitivity to Modelling Inter-lamellar Behaviour With Cohesive Contact

A comparison of percentage change in peak force between the glued, and cohesive inter-lamellar contact models is shown in Figure 6B. Allowing cohesive contact caused a decrease in stiffness in comparison to the glued models for all six geometries. There were only small differences in the percentage reduction in peak force compared to the baseline models between models, with an average reduction of 32.2 ± 2.2%.

Sensitivity to Altered NP:AF Ratio

Increasing the NP:AF ratio resulted in a reduction in peak force in all models in comparison to the baseline runs. Models 1.1 and 2.1 (simple EP, linear NP-AF boundary and segmented EP, linear NP-AF boundary) saw the greatest reduction in peak force in comparison to the curved, and segmented internal geometry models with the same endplate geometry (1.2 and 1.3, and 2.2 and 2.3, respectively). Reductions in peak force were higher in the segmented endplate models, compared to the simple endplate models (13.3 ± 4.7% vs. 7.2 ± 2.3%, respectively).

Discussion

FE modellers are required to find a balance between geometric accuracy and keeping computational complexity low. It is therefore common to simplify model geometry to increase the likelihood of convergence. In this study ultra-high field MRI was used to obtain accurate internal geometry of an intact IVD to evaluate its effect on the response of an IVD FE model.

The six models developed in this study all had geometry based on the same MR image of a bovine IVD, however, differences were seen in terms of the location of peak strain (Figure 5), and the overall stiffness of each model (Figure 4). Maximum strains were seen at locations of high curvature in the segmented endplate models, which could not be captured in the curved endplate models. An increase in model complexity resulted in increased model stiffness, with the average peak force doubling in the segment endplate models compared to the curved endplate models (919 ± 162 N and 1,845 ± 465 N, respectively), and the average peak force being 1.44, and 1.58 times larger than the linear NP-AF boundary model (Model 1.1—simple EP, linear NP-AF boundary and Model 2.1—segmented EP, linear NP-AF boundary) in the curved (Model 1.2—simple EP, curved NP-AF boundary and segmented EP, Model 2.3—curved NP-AF boundary) and segmented (Model 1.3—simple EP, segmented NP-AF boundary and Model 2.3—segmented EP, segmented NP-AF boundary) NP-AF boundary models, respectively. The high stiffnesses in the segmented endplate models was affected by the average height being 1.22 mm shorter than the curved endplate models (7.92 vs. 9.14 mm). Even though the central, and peripheral discs heights were the same in all models the change in average height meant that the overall applied strain was greater in the segmented endplate models compared to the simple endplate models. The differences in stiffness between the segmented and simple endplate models may have been smaller had the average, rather than just central and peripheral disc heights been kept consistent between the two approaches. Additionally, a linear or curved estimation of the NP-AF boundary created idealised strain distributions (Figure 5) that neglected the effects of non-uniform lamella geometries such as variations in width and curvatures.

The force-displacement response of Model 1.1 (simple endplates, linear NP-AF boundary) was closest to the experimental data obtained by Newell et al. (2017; Figure 4). This was expected since Model 1.1’s geometry was the most similar to the models in that study in that the endplates were curved, the boundary between the NP and AF was linear, and the AF fibre YM used in this study (415 MPa) was obtained through the optimisation process described in Newell et al. (2017) that ensured a close match between experimental and numerical results. The slightly stiffer response of Model 1.1 compared the experimental data is likely due to the relatively lower NP-AF ratios, obtained from the MRIs, being used in this study (0.66:1 for Model 1.1) compared to 3.72:1 used in Newell et al. (2017). As shown in Figure 7C, increasing the NP-AF ratio decreases model stiffness, therefore had a lower NP-AF ratio been used by Newell et al. (2017), their optimised fibre stiffnesses may have been lower and therefore a closer match between experimental and numerical response may have been seen in this study.

An increase in deviation from the experimental data was seen with increasing geometric complexity (Figure 4). This is likely due to all models using the same material properties, including AF fibre YM properties that were obtained in an optimisation study (Newell et al., 2017) where model geometry was most similar to Model 1.1 (simple EP, linear NP-AF boundary) in this study. Using the same optimisation algorithm described in (Newell et al., 2017) the material properties (AF fibre YM and AF C₁₀) of all the models used in this study can also be optimised to obtain a close match between numerical and experimental response. This results in the optimised values shown in Figure 8 where lower AF fibre YM and AF C₁₀ values were seen in the more complex models. For researchers who wish to model IVDs with accurate internal geometries it is likely that less stiff material properties are required to obtain an overall response that can be validated against experimental data. For example, the optimised AF fibre YM and AF C₁₀ values for Model 2.3 (segmented endplates, segmented NP-AF boundary) were 60.63 and 0.121 MPa, respectively, which is at the lower end of the range of values that have been used in previous FE studies [44–500 MPa for AF Fibre YM and 0.0146–0.7 MPa for the AF C₁₀ value (Newell et al., 2019b)].

FIGURE 8

Figure 8. Final (A) AF fibre YM (B) C₁₀ values following optimisation of these parameters in order to achieve a good fit between the computational and experimental data reported by Newell et al. (2017).

Irrespective of geometry, the models were most sensitive to the fibre YM, and the fibre angle (on average a 10.5 ± 1.8%, and 21.6 ± 12.0% change in peak force compared to the baseline runs when adjusted by 20% of the original fibre YM and fibre angle, respectively), and relatively insensitive to the AF ground matrix properties (C₁₀ and C₀₁ Mooney constants—Figure 6). Although care must be taken when comparing absolute values since the sensitivity study carried out here spanned a supra-physiologic range of force these findings are similar to those of Newell et al. (2017) in terms of Fibre YM who modelled bovine discs with a curved endplate, and linear NP-AF boundary (on average a 12.9% change in peak force). Interestingly, graduating the fibre bundle angles from inner to outer AF had a relatively small effect on the overall stiffness of the model with the average peak force of all models being 5.0 ± 2.8% lower compared to baseline runs, and the peak force being just 2.2% lower in the segmented endplate, segmented NP-AF boundary model (2.3). This suggests that future studies should focus on ensuring the accuracy of the average fibre bundle angle, rather than the accuracy of fibre bundle angles in individual lamella to ensure accurate load transfer through the IVD, particularly if a segmented NP-AF boundary and segmented endplates are used. However, previous FE studies have reported sensitivity to fibre angles in terms of IVD response to torsion (Yang and O’Connell, 2017), swelling (Yang and O’Connell, 2018, 2019a), and axial stiffness (Shirazi-Adl, 1989), with the general consensus being that anatomically relevant fibre angles are important for understanding the internal stress distributions through the disc. As demonstrated in this study, and in Stadelmann et al. (2018) it is possible to use high resolution MRI to obtain these angles non-invasively. There have been two recent FE studies investigating the effect of NP:AF ratio with Yang et al. (2019) finding a positive correlation between relative NP size and IVD stiffness, and Mengoni et al. (2017) finding a negative correlation. These differences could be due to modelling bovine rather than human discs, fibres at an orientation of ± 20 degrees compared to ± 43 degrees, or modelling 20 lamella layers compared to one in the Mengoni et al. (2017) and Yang et al. (2019) studies, respectively. In this study, doubling the NP:AF ratio (equivalent to on average, increasing the NP:Disc diameter from 0.38 ± 0.01 to 0.55 ± 0.01) reduced the peak force by 10.3 ± 4.7% supporting the findings of Mengoni et al. (2017) who saw a ∼32% decrease in peak force of bovine specimens at a similar axial displacement as the maximum in this study (0.86 mm—taken from Figure 6 in Mengoni et al., 2017) when increasing the NP:Disc diameter from 0.4 to 0.6. Conversely, Yang et al. (2019) found a 23.8% increase in normalised stiffness of human IVDs under axial compression when increasing the NP:Disc diameter from 0.35 to 0.6. Depending on the method of modelling AF fibres, adjusting the NP:AF ratio can affect the total fibre volume, since IVD FE models are sensitive to fibre properties (Figure 6C) it is important to ensure that a change in total fibre volume does not affect conclusions made about altering the NP:AF ratio. In this study this was accounted for by ensuring the fibre volume was consistent before and after adjusting the NP:AF ratio.

A network of elastic fibres and collagen cross-bridges exist between lamellae in the AF which provides resistance to shearing strains (Yu et al., 2002; Pezowicz et al., 2006; Schollum et al., 2008). Most current IVD FE studies omit this inter-lamellar behaviour and either glue the boundary between lamellae (Dreischarf et al., 2014; Yang and O’Connell, 2019a, b; Yang et al., 2019), or use continuum models (Jacobs et al., 2014; Showalter et al., 2016; Newman et al., 2021). Adam et al. (2015) investigated allowing adjacent lamellae to slide freely across each other, however, experimental studies have shown that inter-lamella shearing strain is due to skewing, rather than sliding (Michalek et al., 2009; Vergari et al., 2016). A number of studies have included inter-lamellar interactions in models of several lamellar layers through the incorporation of elements or fibres in a zone between lamellar (Labus et al., 2014; Derrouiche et al., 2019; Kandil et al., 2019, 2020; Ghezelbash et al., 2021; Tamoud et al., 2021), however, the computational complexity of some of these approaches can render them impractical for modelling the whole disc (Ghezelbash et al., 2021). Mengoni et al. (2015) derived normal and tangential cohesive stiffness values that have been assigned to cohesive elements in this study to represent inter-lamellar interactions. To our knowledge this is the first time that these elements have been applied to a full IVD model and the technique proved to have potential to model the inter-lamellar behaviour more physiologically than glued contact between lamellae. Introducing this cohesive behaviour reduced the IVD stiffness (on average 32 ± 2% reduction in peak force compared to the glued baseline models – Figure 7B), which, as expected is lower than the 40% reduction shown by Adam et al. when allowing total free sliding between lamellae.

Bovine samples were modelled in this study because they are almost perfectly round (O’Connell et al., 2007; Adam et al., 2015) providing the opportunity to model axisymmetrically and thus reducing computational cost compared to full 3-D models. Beckstein et al. (2008) found the normalised stiffness of bovine IVDs to be within 12% of human IVDs suggesting similarities in mechanical properties, however, future studies to investigate how the findings in this study relate to human IVDs is required, particularly since the internal geometry of just one bovine IVD has been modelled here. Modelling one IVD was sufficient to investigate the effect of the different modelling techniques deployed in this study, however, a study on a larger population with various internal geometries, for example investigating how the internal geometry changes with degeneration and how that affects FE model behaviour would be of interest. Modelling axisymmetrically significantly reduces computational time but does not allow modes of loading other than axial compression. This meant that some geometric intricacies were simplified, and the NP was assumed to be perfectly in the centre of the IVD where in fact it was offset by approximately 1.26 mm from the centre of the disc when measured on a mid-transverse slice of the MRIs. Yang et al. (2019) found little difference in disc joint stiffness under flexion, extension, and lateral bending when changing the NP:Disc area ratio from 0.21 to 0.6, however the study did not investigate the effect of endplate and internal geometry complexity under these modes of loading.

In this study IVDs were loaded at an intermediate strain rate of 0.1/s. Previous studies have demonstrated an increase in stiffness with strain rate (Smeathers and Joanes, 1988; Yingling et al., 1997; Race et al., 2000; Kemper et al., 2007; Costi et al., 2008; Newell et al., 2017), and a recent study has demonstrated that the NP has little effect on the response of IVDs at high loading rates (Newell et al., 2019a), which differs from its function at low strain rates (Seroussi et al., 1989; Meakin and Hukins, 2000). It is therefore possible that FE models are less sensitive to geometric simplifications at higher strain rates, although further work would be required to confirm this. In this study IVDs were only subjected to pure axial compression. In order to comprehensively understand the effect of internal geometry on the outcomes of FE models, future studies should extend this analysis to investigate IVD response in combined loading and in flexion-extension, axial rotation, and lateral bending.

The complexity of the models used in this study were deliberately kept low to reduce the effects of confounding variables such as cartilage endplate properties, NP fluid phase (poroelasticity), NP swelling pressure, preload, vertebral bone properties, changes in the fibre volume fraction through the AF, and asymmetries in endplate shapes in planes other than the sagittal from which the geometry was segmented for Models 2.1—2.3 in this study. Additionally, the MRIs were obtained while the sample was well hydrated, and musculature had been removed. In vivo, internal geometries may change depending on posture and the time of the day and therefore the geometry used here represents the geometry at a single time point, and a single posture. Although outside the focus of this study, it is recommended that these factures, as well as those highlighted in this study are considered when developing more complex patient specific models of human discs that aim to replicate in vivo conditions.

Conclusion

Geometric simplifications in FE models of IVDs create idealised load transfer that may not be physiologic. These simplifications affect model response, particularly in terms of stiffness and strain distributions, sensitivity to average fibre angles, and sensitivity to modelling inter-lamellar contact. Therefore, defining more realistic internal and external geometry of the IVD can significantly affect IVD FE model response, and it is likely that a more realistic geometry leads to a more accurate strain distribution within the IVD. It is recommended that these geometric intricacies are incorporated into IVD FE models before material properties are optimised to develop a validated model. This is particularly important if the models are being used for clinical applications such as developing repair strategies that aim to replicate the mechanical behaviour of healthy discs.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

YD, UH, and NN designed the study. NB, ST, and NN designed the sequences and obtained the MRI images. YD and TR developed and ran the FE models. YD, ST, and TR analysed the data and drafted the manuscript which was edited by NN, NB, and UH. All authors approved the manuscript before submission.

Funding

Part of this work was funded by an Imperial College Research Fellowship for NN and an EPSRC DTP CASE Conversion Studentship for TR (EP/R513052/1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with one of the authors NN.

References

Adam, C., Rouch, P., and Skalli, W. (2015). Inter-lamellar shear resistance confers compressive stiffness in the intervertebral disc: an image-based modelling study on the bovine caudal disc. J. Biomech. 48, 4303–4308.

Google Scholar

Adams, M. A., and Roughley, P. J. (2006). What is intervertebral disc degeneration, and what causes it? Spine (Phila Pa 1976) 31, 2151–2161. doi: 10.1097/01.brs.0000231761.73859.2c

CrossRef Full Text | Google Scholar

Adams, M. A., Dolan, P., and Hutton, W. C. (1987). Diurnal variations in the stresses on the lumbar spine. Spine (Phila Pa 1976) 12, 130–137.

Google Scholar

Adams, M. A., McNally, D. S., and Dolan, P. (1996). “Stress” distributions inside intervertebral discs. The effects of age and degeneration. J. Bone Joint Surg. Br. 78, 965–972.

Google Scholar

Ayturk, U. M., and Puttlitz, C. M. (2011). Parametric convergence sensitivity and validation of a finite element model of the human lumbar spine. Comput. Methods Biomech. Biomed. Eng. 14, 695–705. doi: 10.1080/10255842.2010.493517

PubMed Abstract | CrossRef Full Text | Google Scholar

Beckstein, J. C., Sen, S., Schaer, T. P., Vresilovic, E. J., and Elliott, D. M. (2008). Comparison of animal discs used in disc research to human lumbar disc: axial compression mechanics and glycosaminoglycan content. Spine (Phila Pa 1976) 33, E166–E173. doi: 10.1097/BRS.0b013e31824d911c

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogduk, N. (2005). Clinical Anatomy of the Lumbar Spine and Sacrum, 4th Edn. London: Elsevier Health Sciences.

Google Scholar

Cassidy, J. J., Hiltner, A., and Baer, E. (1989). Hierarchical structure of the intervertebral disc. Connect. Tissue Res. 23, 75–88. doi: 10.3109/03008208909103905

PubMed Abstract | CrossRef Full Text | Google Scholar

Costi, J. J., Stokes, I. A., Gardner-Morse, M. G., and Iatridis, J. C. (2008). Frequency-dependent behaviour of the intervertebral disc in response to each of six degree of freedom dynamic loading: solid phase and fluid phase contributions. Spine (Phila Pa 1976) 33, 1731–1738. doi: 10.1097/BRS.0b013e31817bb116

PubMed Abstract | CrossRef Full Text | Google Scholar

Derrouiche, A., Zaïri, F., and Zaïri, F. (2019). A chemo-mechanical model for osmo-inelastic effects in the annulus fibrosus. Biomech. Model. Mechanobiol. 18, 1773–1790. doi: 10.1007/s10237-019-01176-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Dreischarf, M., Zander, T., Shirazi-Adl, A., Puttlitz, C., Adam, C., Chen, C., et al. (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, 1757–1766. doi: 10.1016/j.jbiomech.2014.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghezelbash, F., Eskandari, A. H., Shirazi-Adl, A., Kazempour, M., Tavakoli, J., Baghani, M., et al. (2021). Modeling of human intervertebral disc annulus fibrosus with complex multi-fiber networks. Acta Biomater. 123, 208–221. doi: 10.1016/j.actbio.2020.12.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Iatridis, J. C., Kumar, S., Foster, R. J., Weidenbaum, M., and Mow, V. (1999). Shear mechanical properties of human lumbar annulus fibrosus. J. Orthop. Res. 17, 732–737. doi: 10.1002/jor.1100170517

PubMed Abstract | CrossRef Full Text | Google Scholar

Iatridis, J. C., Setton, L. A., Foster, R. J., Rawlins, B. A., Weidenbaum, M., and Mow, V. C. (1998). Degeneration affects the anisotropic and nonlinear behaviors of human anulus fibrosus in compression. J. Biomech. 31, 535–544. doi: 10.1016/S0021-9290(98)00046-3

CrossRef Full Text | Google Scholar

Jacobs, N. T., Cortes, D. H., Peloquin, J. M., Vresilovic, E., and Elliott, D. (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, 2540–2546. doi: 10.1016/j.jbiomech.2014.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandil, K., Zaïri, F. F., Derrouiche, A., and Messager, T. (2019). Interlamellar-induced time-dependent response of intervertebral disc annulus: a microstructure-based chemo-viscoelastic model. Acta Biomater. 100, 75–91. doi: 10.1016/j.actbio.2019.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandil, K., Zaïri, F., Messager, T., and Zaïri, F. (2020). Interlamellar matrix governs human annulus fibrosus multiaxial behavior. Sci. Rep. 10:19292. doi: 10.1038/s41598-020-74107-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kemper, A. R., McNally, C., and Duma, S. M. (2007). The influence of strain rate on the compressive stiffness properties of human lumbar intervertebral discs. Biomed. Sci. Instrum. 43, 176–181.

Google Scholar

Labus, K. M., Han, S. K., Hsieh, A. H., and Puttlitz, C. M. (2014). A computational model to describe the regional interlamellar shear of the annulus fibrosus. J. Biomech. Eng. 136:51009. doi: 10.1115/1.4027061

CrossRef Full Text | Google Scholar

Marchand, F., and Ahmed, A. M. (1990). Investigation of the laminate structure of lumbar disc anulus fibrosus. Spine (Phila Pa 1976) 15, 402–410. doi: 10.1097/00007632-199005000-00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Meakin, J. R. R., and Hukins, D. W. L. W. (2000). Effect of removing the nucleus pulposus on the deformation of the annulus fibrosus during compression of the intervertebral disc. J. Biomech. 33, 575–580. doi: 10.1016/S0021-9290(99)00215-8

CrossRef Full Text | Google Scholar

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, 929–935. doi: 10.1097/BRS.0b013e3181fd7f7f

PubMed Abstract | CrossRef Full Text | Google Scholar

Mengoni, M., Kayode, O., Sikora, S. N. F., Zapata-Cornelio, F. Y., Gregory, D., and Wilcox, R. (2017). Annulus fibrosus functional extrafibrillar and fibrous mechanical behaviour: experimental and computational characterisation. R. Soc. Open Sci. 4:170807. doi: 10.1098/rsos.170807

PubMed Abstract | CrossRef Full Text | Google Scholar

Mengoni, M., Luxmoore, B. J., Wijayathunga, V. N., Jones, A., Broom, N., and Wilcox, R. (2015). Derivation of inter-lamellar behaviour of the intervertebral disc annulus. J. Mech. Behav. Biomed. Mater. 48, 164–172. doi: 10.1016/j.jmbbm.2015.03.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Michalek, A. J., Buckley, M. R., Bonassar, L. J., Cohen, I., and Iatridis, J. (2009). Measurement of local strains in intervertebral disc anulus fibrosus tissue under dynamic shear: contributions of matrix fiber orientation and elastin content. J. Biomech. 42, 2279–2285. doi: 10.1016/j.jbiomech.2009.06.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Mooney, M. (1940). A theory of large elastic deformation. J. Appl. Phys. 11, 582–592.

Google Scholar

Natali, A., and Meroi, E. (1990). Nonlinear analysis of intervertebral disk under dynamic load. J. Biomech. Eng. 112, 358–363. doi: 10.1115/1.2891196

CrossRef Full Text | Google Scholar

Natarajan, R. N., and Andersson, G. B. J. (1999). The influence of lumbar disc height and cross-sectional area on the mechanical rsponse of the disc to physiological loading. Spine (Phila Pa 1976) 24, 1873–1881.

Google Scholar

Newell, N., Carpanen, D., Evans, J. H., Pearcy, M. J., and Masouros, S. D. (2019a). Mechanical function of the nucleus pulposus of the intervertebral disc under high rates of loading. Spine (Phila Pa 1976) 44, 1035–1041.

Google Scholar

Newell, N., Carpanen, D., Grigoriadis, G., Little, J. P., and Masouros, S. D. (2019b). Material properties of human lumbar intervertebral discs across strain rates. Spine J. 19, 2013–2024.

Google Scholar

Newell, N., Grigoriadis, G., Christou, A., Carpanen, D., and Masouros, S. (2017). Material properties of bovine intervertebral discs across strain rates. J. Mech. Behav. Biomed. Mater. 65, 824–830.

Google Scholar

Newman, H. R., DeLucca, J. F., Peloquin, J. M., Vresilovic, E. J., and Elliott, D. M. (2021). Multiaxial validation of a finite element model of the intervertebral disc with multigenerational fibers to establish residual strain. JOR Spine e1145. doi: 10.1002/jsp2.1145

CrossRef Full Text | Google Scholar

Niemeyer, F., Wilke, H. J., and Schmidt, H. (2012). Geometry strongly influences the response of numerical models of the lumbar spine-A probabilistic finite element analysis. J. Biomech. 45, 1414–1423. doi: 10.1016/j.jbiomech.2012.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Noailly, J., Wilke, H.-J., Planell, J. A., and Lacroix, D. (2007). How does the geometry affect the internal biomechanics of a lumbar spine bi-segment finite element model? Consequences on the validation process. J. Biomech. 40, 2414–2425. doi: 10.1016/J.JBIOMECH.2006.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Connell, G. D. G., Guerin, H. L. H., and Elliott, D. M. (2009). Theoretical and uniaxial experimental evaluation of human annulus fibrosus degeneration. J. Biomech. Eng. 131:111007. doi: 10.1115/1.3212104

CrossRef Full Text | Google Scholar

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, 328–333. doi: 10.1097/01.brs.0000253961.40910.c1

CrossRef Full Text | Google Scholar

O’Connell, G., Vresilovic, E., and Elliott, D. (2011). Human intervertebral disc internal strain in compression: the effect of disc region, loading position, and degeneration. J. Orthop. Res. 29, 547–555.

Google Scholar

Pezowicz, C. A., Robertson, P. A., and Broom, N. D. (2006). The structural basis of interlamellar cohesion in the intervertebral disc wall. J. Anat. 208, 317–330. doi: 10.1111/j.1469-7580.2006.00536.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Race, A., Broom, N. D., and Robertson, P. (2000). Effect of loading rate and hydration on the mechanical properties of the disc. Spine (Phila Pa 1976) 25, 662–669. doi: 10.1097/00007632-200003150-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Robin, S., Skalli, W., and Lavaste, E. (1994). Influence of geometrical factors on the behaviour of lumbar spine segments: a finite element analysis. Eur. Spine J. 3, 84–90.

Google Scholar

Schmidt, H., Galbusera, F., Rohlmann, A., Zander, T., and Wilke, H. (2012). Effect of multilevel lumbar disc arthroplasty on spine kinematics and facet joint loads in flexion and extension: a finite element analysis. Eur. Spine J. 21, S663–S674. doi: 10.1007/s00586-010-1382-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Schollum, M. L., Robertson, P. A., and Broom, N. D. (2008). ISSLS prize winner: microstructure and mechanical disruption of the lumbar disc annulus. Spine (Phila Pa 1976) 33, 2702–2710. doi: 10.1097/BRS.0b013e31817bb92c

PubMed Abstract | CrossRef Full Text | Google Scholar

Seroussi, R. E., Krag, M. H., Muller, D. L., and Pope, M. H. (1989). Internal deformations of intact and denucleated human lumbar discs subjected to compression, flexion, and extension loads. J. Orthop. Res. 7, 122–131. doi: 10.1002/jor.1100070117

PubMed Abstract | CrossRef Full Text | Google Scholar

Shirazi-Adl, A. (1989). On the fibre composite material models of disc annulus – comparison of predicted stresses. J. Biomech. 22, 357–365.

Google Scholar

Shirazi-Adl, A., Ahmed, A. M., and Shrivastava, S. C. (1986). A finite element study of a lumbar motion segment subjected to pure sagittal plane moments. J. Biomech. 19, 331–350. doi: 10.1016/0021-9290(86)90009-6

CrossRef Full Text | Google Scholar

Showalter, B. L., DeLucca, J. F., Peloquin, J. M., Cortes, D. H., Yoder, J. H., Jacobs, N. T., et al. (2016). Novel human intervertebral disc strain template to quantify regional three-dimensional strains in a population and compare to internal strains predicted by a finite element model. J. Orthop. Res. 34, 1264–1273. doi: 10.1002/jor.23137

PubMed Abstract | CrossRef Full Text | Google Scholar

Smeathers, J. E., and Joanes, D. N. (1988). Dynamic compressive properties of human lumbar intervertebral joints: a comparison between fresh and thawed specimens. J. Biomech. 21, 425–433. doi: 10.1016/0021-9290(88)90148-0

CrossRef Full Text | Google Scholar

Stadelmann, M. A., Maquer, G., Voumard, B., Grant, A., Hackney, D., Vermathen, P., et al. (2018). Integrating MRI-based geometry, composition and fiber architecture in a finite element model of the human intervertebral disc. J. Mech. Behav. Biomed. Mater. 85, 37–42. doi: 10.1016/j.jmbbm.2018.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamoud, A., Zaïri, F., Mesbah, A., and Zaïri, F. (2021). A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus. Meccanica 56, 585–606. doi: 10.1007/s11012-020-01281-4

CrossRef Full Text | Google Scholar

Tavana, S., Baxan, N., Masouros, S. D. D., Freedman, B. A., Hansen, U. N., and Newell, N. (2021). The effect of degeneration on internal strains and the mechanism of failure in human intervertebral discs analysed using Digital Volume Correlation (DVC) and Ultra-High Field MRI. Front. Bioeng. Biotechnol. 8:610907. doi: 10.3389/fbioe.2020.610907

PubMed Abstract | CrossRef Full Text | Google Scholar

Tavana, S., Clark, J. N. N., Prior, J., Baxan, N., Masouros, S. D., Newell, N., et al. (2020). Quantifying deformations and strains in human intervertebral discs using Digital Volume Correlation combined with MRI (DVC-MRI). J. Biomech. 102:109604. doi: 10.1016/j.jbiomech.2020.109604

PubMed Abstract | CrossRef Full Text | Google Scholar

Vergari, C., Mansfield, J., Meakin, J. R., and Winlove, P. C. (2016). Lamellar and fibre bundle mechanics of the annulus fibrosus in bovine intervertebral disc. Acta Biomater. 37, 14–20. doi: 10.1016/j.actbio.2016.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., and O’Connell, G. D. (2017). Effect of collagen fibre orientation on intervertebral disc torsion mechanics. Biomech. Model. Mechanobiol. 16, 2005–2015. doi: 10.1007/s10237-017-0934-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., and O’Connell, G. D. (2018). Swelling of fiber-reinforced soft tissues is affected by fiber orientation, fiber stiffness, and lamella structure. J. Mech. Behav. Biomed. Mater. 82, 320–328. doi: 10.1016/j.jmbbm.2018.03.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., and O’Connell, G. D. (2019a). GAG content, fiber stiffness, and fiber angle affect swelling-based residual stress in the intact annulus fibrosus. Biomech. Model. Mechanobiol. 18, 617–630. doi: 10.1007/s10237-018-1105-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., and O’Connell, G. D. (2019b). Intervertebral disc swelling maintains strain homeostasis throughout the annulus fibrosus: a finite element analysis of healthy and degenerated discs. Acta Biomater. 100, 61–74. doi: 10.1016/j.actbio.2019.09.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., Lu, Y., Um, C., and O’Connell, G. (2019). Relative nucleus pulposus area and position alters disc joint mechanics. J. Biomech. Eng. 141, 1–11. doi: 10.1115/1.4043029

CrossRef Full Text | Google Scholar

Yingling, V. R., Callaghan, J. P., and McGill, S. M. (1997). Dynamic loading affects the mechanical properties and failure site of porcine spines. Clin. Biomech. 12, 301–305. doi: 10.1016/S0268-0033(97)00009-0

CrossRef Full Text | Google Scholar

Yu, J., Peter, C., Roberts, S., and Urban, J. (2002). Elastic fibre organization in the intervertebral discs of the bovine tail. J. Anat. 201, 465–475.

Google Scholar