On the relative relevance of subject-specific geometries and degeneration-specific mechanical properties for the study of cell death in human intervertebral disk models

Introduction

Early changes in the optimal biophysics of the intervertebral disk (IVD) are thought to be one of the major causes of degeneration, in contrast to the normal aging processes (Smith et al., 2011). As such, these changes are related to the long-term pathologies of the human spine and to the derived socio-economical burden. However, the complex IVD environment is difficult to study, and the explorations of cause and effect relationships in disk degenerative disease (DDD) still needs more efforts.

For instance, the nutritional supply for disk cells has been pointed out as a possibly important factor involved in DDD (Huang et al., 2014). The IVD is the largest avascular organ in humans and the metabolism of the relatively few disk cells has to ensure the maintenance of a large amount of extracellular matrix (ECM), while it mainly relies on proper diffusion of metabolites from and to the periphery of the disk. Such a situation can have deleterious effects in the innermost regions of the disk (Huang et al., 2014), which might be worsen by the fact that disk cells seem particularly prone to catabolic responses in presence of limited nutrition (Rinkler et al., 2010; Neidlinger-Wilke et al., 2012). Yet, disk cell biology should be considered along with the mechanical competence of the disk tissues.

For example, preliminary numerical analyses (Malandrino et al., 2013) suggested that mechanical compaction of the ECM could alter the metabolic transport throughout the disk, being this alteration a possible predominant cause of cell death under mechanical overloads. Reported as “indirect mechanotransduction,” this phenomenon has been pointed out as a possible important source of disturbance for ECM maintenance (Iatridis et al., 2006). However, the interplay between cell nutrition and IVD mechanics is difficult to explore in an integrated way.

Among all research techniques, numerical modeling stands for a promising tool to better understand the complex mechanobiology of the musculoskeletal system (Halloran et al., 2012; Henak et al., 2013). Existing finite element (FE) spine models revealed a growing capacity to integrate coupled effects of mechanical and biophysical phenomena, and of specific tissue and geometrical properties (Noailly and Lacroix, 2012). While the relation between spine biomechanics and implants has been an important target, the possibility of investigating the effect of non-mechanical factors in DDD and regenerative therapies is another major advantage of numerical models (Noailly et al., 2014).

Focusing on the important question whether nutrition could be a relevant actor in IVD pathogenesis, several FE models have been created in order to explore the implication of indirect mechanotransduction under different simulated loads and tissue conditions. Among the different results achieved, it was found that mechanical deformations might favor cell nutrition as long as disk tissue properties remain healthy (Malandrino et al., 2011). In contrast, FE analyses have suggested that tissue degenerative changes might induce nutrition-related cell death under physiological load magnitudes, when the boundary supply of nutrients is reduced, e.g., through hypothetical endplate sclerosis (Zhu et al., 2012; Malandrino et al., 2014b). As for treatment therapies, coupling cell nutrition-poro-mechanical models have been used to explore possible anabolic treatments, e.g., with growth factors, for regenerative strategies (Huang et al., 2012; Travascio et al., 2014).

Though model validation is an unresolved issue for such theoretical studies, the possibility to generate new educated guesses through FE simulations is a clear asset in order to define new experiments or alternative patient explorations/clinical result analyses. However, this asset can only be duly exploited if model interpretations are built on a sufficiently broad view of the range of results that may arise from the variability of the input hypotheses. Most of the disk poro-mechanical models coupled to cell nutrition models are based on different osmo-poro-hyperelastic formulations initially developed to simulate healthy tissues. The poro-mechanical response of degenerated tissues was then typically simulated by decreasing the porosity and the capacity to develop osmotic pressure, and by increasing the stiffness constants (Malandrino et al., 2011; Zhu et al., 2012). However, these models of degenerated tissues have never been validated to our knowledge. Disk degeneration involves proteoglycan depletion, a reduced production of collagen II that competes with tissue fibrosis, and cracks occurrence (Antoniou et al., 1996, 2013). Hence, the multiphysics changes that occur in degenerating tissues are complex, and might be difficult to capture only through simple parameter value adjustments. As such, we raise the question whether a poro-mechanical constitutive equation initially validated to simulate healthy tissues can simulate the degenerative changes cited above.

As for the geometrical effects, it has been shown that the prediction of solutes diffusion within IVD models depends on the specific geometry assumed (Motaghinasab et al., 2014). Disk height reduction has often been considered in the simulation of degenerated changes (Galbusera et al., 2011a; Malandrino et al., 2011, 2014b; Zhu et al., 2012), and evaluated as beneficial to cell nutrition (Galbusera et al., 2011a; Malandrino et al., 2014b). However, most of the models represented disk heights of about 10 mm in their healthy configuration, which is a typical lumbar disk dimension for the average population (Noailly et al., 2007), but might change patient-specifically. Actually, measuring disk height accurately is not straightforward, and results largely depend on the landmarks chosen from plane MRI or radiographs (Frobin et al., 1997; Tunset et al., 2013). Moreover, whether the respective progresses of disk height reduction and degeneration are correlated to each other has been questioned (Frobin et al., 2001). All these geometrical uncertainties about the healthy and degenerated IVD suggest that different interactions might be predicted between disk cell nutrition and disk degeneration, depending on the patient considered. Hence, quantifications of patient-specific 3D disk geometries (Neubert et al., 2014) along with multiphysics disk tissue changes appear necessary, in order to capture these interactions.

Accordingly, the present study stands for a preliminary exploration of patient- and degeneration-specific geometrical and mechanical factors in order to investigate possible onsets of disk degeneration related to cell nutrition issues. First, the ability of IVD poro-mechanical modeling to simulate disk degeneration effects was assessed. We focused on the physical significance of degeneration-dependent tissue model parameters values, achieved through an inverse optimization problem, i.e., by minimizing the differences between the predictions of specimen-specific geometrical models and in vitro measurements. In particular, we considered the need to simulate the effect of cracks on the calculated multiphysics response of the disk tissues. The best degeneration-specific material parameters found were used for an early examination of the relative importance of the degeneration-dependent mechanical competence of the ECM, in particular disk geometries.

Materials and Methods

Cadaveric Specimens

Four cadaveric human lumbosacral spines (L1-S1), i.e., Sp#1, Sp#2, Sp#4, and Sp#6 were provided by the Anatomy Department of the University Medical Center Utrecht, The Netherlands (Table 1). No additional information on the medical history was provided. Postmortem, specimens were stored at −25°C until the moment of use. L1-S1 specimens were harvested including all ligaments and para-spinal muscles. Anterior–posterior and lateral digital radiographs of the specimens were used to exclude severely degenerated spines. The specimens were scanned with both MRI and qCT for both grading and reconstruction protocols. CT protocol was performed at Polyclinique Saint Côme (Compiègne, France). Moreover, they were scanned in sagittal plane using the T2-weighted MRI protocol on a 1.5 T GE MRI machine. The slice thickness was 3 mm, the matrix was 512 × 512, and the field of view (FOV) was 360 mm × 360 mm. The repetition time (TR) and echo time (TE) were 3380 and 105.588 ms, respectively. The number of images was 19 and the acquisition time was about 3 min. All IVDs of all the L1-S1 specimens were graded by an experienced radiologist in terms of degeneration Pfirrmann score (Pfirrmann et al., 2001).

TABLE 1

Table 1. Specimen information.

Experimental Protocol

One L1-S1 specimen (Sp#4) was tested in vitro under flexion–extension, lateral bending, and axial rotation (see also Validation Study). Before testing, the specimens were allowed to thaw at room temperature overnight. Then, all muscles were carefully excised to prevent damaging the ligaments and cartilaginous tissues. The L1 and S1 vertebrae were embedded in Polymethyl Methacrylate. S1 was fully constrained in the testing machine (Wilke et al., 1994). Pure unconstraint moments of 5 Nm were applied on the upper plate with a constant loading rate of 1 deg/s in flexion–extension and lateral bending and 0.5 deg/s in axial rotation. These load directions were defined as principal directions, and the ranges of motion (RoM) were measured. Previous to measurements, a cycle of 5 Nm unconstrained flexion–extension was applied at 1 deg/s in order to assess the absolute neutral position. Then, two full cycles of preconditioning were applied in each principal direction, under loading conditions similar to those used for the mechanical analysis of the specimens. No follower load was applied, nor was the pressure in the IVDs measured.

Subsequently, creep compression experiments were performed on the L3-L4 segments of the four lumbar spines (Figure 1). The L3-L4 monosegments were excised from the L1-S1 specimens, and positioned on the testing machine in their neutral position. All degrees of freedom were fixed and locked in the gimbal system of the testing machine. Only vertical translations were allowed and measured, in response to the imposed axial compression. A compressive preload of 300 N was applied for 180 s, then released for 180 s. This cycle was repeated three times, in order to minimize the effect of IVD superhydration. Next, the load onto the spinal motion segment (SMS) was increased to 500 N in 10 s, and maintained for 3 h. The creep response was measured in terms of recorded displacements. All mechanical tests were performed in the spine-testing laboratory of the Institute of Orthopedic Research and Biomechanics (Ulm, Germany).

FIGURE 1

Figure 1. (A) Generic and example of a morphed FE models of the full L1-S1 osteo-ligamentous spine and (B) MRI-based FE reconstructions of the four L3-L4 IVDs used in the optimization study (note: the scale is not uniform from one model picture to another. Please refer to Table 5 for the dimensions).

FE Models Based on MR and qCT Images

Applying the protocol developed in the framework of the EU-funded project MySpine, the cadaveric specimen geometries of lumbar vertebrae and IVDs were reconstructed from MR and QCT images, and a generic FE model geometry (Figure 1A) was morphed to the subject-specific geometries. The IVD segmentation method has been presented in Castro-Mateos et al. (2014), and the vertebra segmentation in Castro-Mateos et al. (2015). The morphed FE meshes of the IVDs were representatives of the patient-specific proportions of nucleus pulposus (NP) and annulus fibrosus (AF) as detected from the MR images (Figure 1). The specific shape of the cartilage endplate (CEP) resulted from these adaptations.

The morphing of the generic model onto the 3D geometrical reconstructions of the specimens tested followed a method, the details of which are being reported independently (Pozo et al., in preparation). Briefly, the similarity energy was composed of two symmetric terms measuring the distances between the patient-specific surface geometry and the corresponding surface on the morphed model. The internal regularization energy was based on elastic energies relative to the element edge lengths and element solid angles. The transformation was parameterized with cubic B-splines, and the optimization algorithm was gradient-descent with direct line-search. The morphing method preserved the relative sizes of the elements as well as the mesh structure at material discontinuities, as previously determined through an extensive mesh convergence analysis of the IVD generic model (Ruiz et al., 2013). Such preservation aimed to maximize the convergence and avoid as much as possible poro-mechanical instabilities under physiologic load rates. It also guarantied the non-convexity of the elements.

The model included ligaments with a non-linear elastic tensile behavior, modeled through a power law for the toe region and through a linear stress–strain relationship for the linear region of the ligaments (Noailly et al., 2012). The model has three different groups of element-type: shell elements (vertebras and sacrum), hexahedral elements (disks and facets), and truss elements (ligaments). Facet cartilages and vertebral bodies were modeled as quasi-rigid bodies (i.e., ~30 GPa Young’s modulus), and facet contact was considered frictionless and resolved through a penalty algorithm with a penalty normal stiffness of 200 N/mm (Schmidt et al., 2009).

Constitutive Modeling of the IVD

The AF and the NP were both modeled as poro-hyperelastic materials, and the CEP as a poro-elastic material. In the poro-hyperelastic materials, the total stress tensor $\underset{=}{\sigma }$ generated by external loads was the superimposition of the porous solid stress and the fluid pore pressure, p, which were respectively derived from a strain energy density function W, and from Darcy’s law:

$W=\frac{G}{2}\left(I_1-3\right)+\frac{K}{2}{\left(J-1\right)}^2+W_{\mathrm{ANI}}$ (1)

$\underset{=}{\sigma }=\frac{1}{J}\frac{\partial W}{\partial \underset{=}{F}}{\underset{=}{F}}^{\mathrm{T}}-p\underset{=}{I}$ (2)

${\underset{¯}{u}}_{\text{f}}\mathrm{\varphi }=\underset{¯}{\underset{¯}{k}}\cdot \underset{¯}{\nabla p}$ (3)

In the above equations, G and K are respectively the shear and the bulk modulus of the drained porous solid skeleton. W_ANI is an anisotropic strain energy density term, different from zero only for the AF (see Eq. 6). $J=\det \underset{=}{F}$ is the volumetric strain with $\underset{=}{F}$ the deformation gradient tensor of the continuum, I₁ is the first strain invariant, $\underset{=}{I}$ is the second order unit tensor, ${\underset{-}{u}}_{\mathrm{f}}$ is the pore fluid velocity, and ϕ and $\underset{=}{k}$ are respectively the porosity and the hydraulic permeability tensor of the tissue. The two phases, i.e., fluid and solid, were assumed to be nearly incompressible, and the porosity varied with the deformation of the continuum, with respect to an initial value ϕ₀:

$\mathrm{\varphi }=1-J^{-1}(1-{\varphi }_{\text{0}})$ (4)

For all tissues (AF, NP, and CEP), an isotropic hydraulic permeability was considered and depended on the initial porosity, on the volumetric strain, and on an initial permeability, k₀ [5]:

$\underset{=}{k}=\left[k_0\left(\frac{(1-{\varphi }_{\text{0}})\mathrm{\varphi }}{{\mathrm{\varphi }}_0\left(1-\mathrm{\varphi }\right)}\right)\exp \left(\frac{M\left(J-1\right)}{2}\right)\right]\underset{=}{I}$ (5)

being M an empirical coefficient.

In the AF, the anisotropic and non-linear fiber-induced strengthening was taken into account through the W_ANI strain energy density term from Eq. 1 (Gasser et al., 2006):

$W_{\mathrm{ANI}}=\frac{K_1}{K_2}\underset{\mathrm{\delta }=1}{\stackrel{2}{\Sigma }}\left\{\exp \left[K_2{⟨{\overline{E}}_{\mathrm{\delta }}⟩}^2\right]-1\right\}$ (6a)

$\mathrm{with}{\overline{E}}_{\mathrm{\delta }}=\mathrm{\kappa }(I_1-3)+\left(\mathrm{1}-\mathrm{3\kappa }\right)[(a_{0{,}_{\text{δ}}}\otimes a_{0{,}_{\text{δ}}}:\underset{¯}{\underset{¯}{C}})-1]$ (6b)

$\mathrm{and}⟨{\overline{E}}_{\mathrm{\delta }}⟩=\frac{1}{2}\left(\left|{\overline{E}}_{\mathrm{\delta }}\right|+{\overline{E}}_{\mathrm{\delta }}\right)$ (6c)

This anisotropy term depended on fiber stiffness parameters, i.e., K₁ and K₂, and on the strain-like quantity ${\overline{E}}_{\mathrm{\delta }}$, with a_0,δ and $\underset{=}{C}$ being the unit vector in the direction δ and the right Cauchy–Green tensor. Hence, it represented the square of the stretch, active only in two fiber directions (δ = 1,2), respectively characterized by angle values opposite in sign, according to the criss-cross distribution of the fibers in the AF. The parameter κ describes the level of dispersion in the fiber directions, varying from 0 (fully aligned fibers) to 1/3 (fully random fibers). Equation 6c imposes the condition that only fibers with positive strains contribute to the AF stiffness.

Within the NP, proteoglycan-induced NP swelling was described by considering the fluid pressure as a sum of the water chemical potential, u_w, and a swelling pressure related term, ΔΠ, which was assumed constant during deformation (Wilson et al., 2005; Galbusera et al., 2011b):

$p=u_w+\mathrm{\Delta \Pi }$ (7)

Damage Modeling Approach

Undamaged drained shear (G) and bulk (K) moduli were both identical for the AF and for the NP (Schroeder et al., 2010). This equality assumed that the differences in deviatoric and volumetric stiffness between the two healthy tissues arise from the differences in oriented collagen and water contents, at low strains. According to the theoretical developments proposed by Dormieux and Kondo (2007), a dimensionless damage parameter, d, was used to set initially the effective damaged moduli G_d and K_d with both AF- and NP-specific values, based on a Mori–Tanaka homogenization estimate (Mori and Tanaka, 1973) for elastic constant in a medium with open cracks:

$G_d=\frac{G}{1+Q_{G}d}$ (8)

$K_d=\frac{K}{1+Q_{K}d}$ (9)

where $Q_G=\frac{32}{45}\frac{(1-\nu )(5-\nu )}{2-\nu }$ and $Q_K=\frac{16}{9}\frac{1-{\nu }^2}{1-2\nu }$ were only function of the undamaged Poisson’s ratio ν of the solid matrix, assumed to be 0.2 (Ferguson et al., 2004; Johannessen and Elliott, 2005). No damage was assumed for the CEP.

Exploiting the volume-based definition of d as a crack density factor multiplied by the cube of a typical crack length (Budiansky and O’Connell, 1976), we further introduced a modification of this micromechanics-based approach: the density of cracks altered the initial porosity as well, assuming thus that cracks were additional cavities saturated by fluid, and connected to the rest of the pores:

${\mathrm{\varphi }}_{0d}={\mathrm{\varphi }}_0+\frac{4\mathit{Xd}\mathrm{\pi }}{3}$ (10)

being ϕ_{0 d} the damaged initial porosity and X = 0.01 an aspect ratio that simulates ellipsoidal needle-shaped cracks, with a major semi axes two orders of magnitude larger that the minor semi axis.

Optimization of the Poro-Hyperelastic Material and Damage Constants

In order to obtain inversely a set of material properties as a function of the degeneration grade, the respective experimental creep responses of the three L3-L4 IVDs with degeneration grades from two to four were considered as gold standards. For this optimization study, vertebrae were not modeled, and the effect of the posterior elements on the vertical creep displacement was neglected. Thus, only the specimen-specific disk models were used, and for each IVD, an objective function was defined as the average least square difference between the experimental and the predicted creep displacements, u, as follows:

$S_{\mathrm{creep}}=\frac{{\displaystyle \sum _{j=1}^N}{\left(u_{\mathrm{EXP,}j}-u_{\mathrm{FEM,}j}\right)}^2}{N}$ (11)

where N is the number of points considered for the optimization along the experimental creep curve. FE-predicted displacements were obtained by imposing a load of 500 N, as used experimentally, which increased linearly from 0 to 10 s (short-term compressive response) and was subsequently maintained for 3 h (long-term creep response). No preconditioning cycles were simulated: the swelling pressure value in the NP at time zero was considered as a pre-stress in the initial condition, and was kept constant for the remaining simulated time.

Before undertaking the numerical optimizations, a design of experiment (DOE) approach was used to screen out the main factors affecting the displacement response in the short-term and in the long-term. A fractional factorial design was used (Malandrino et al., 2009), and two DOE studies were performed for one specimen-specific L3-L4 model (Sp#2, grade 3): while one study allowed examining the parameters that most affected the short-term compressive, the other tackled the long-term creep response. The influence of each parameter on each response was measured with the ANOVA standardized effect (Montgomery, 2001), and the level of significance α was set to 0.01. The parameter space was defined by: the initial CEP hydraulic permeability (k₀, see Eq. 5), the shear and bulk undamaged moduli (G, K, see Eqs 8 and 9), the AF and NP damage (d, see Eqs 8–10), the AF dispersion fiber (κ, see Eq. 6b), and the NP swelling pressure (ΔΠ, see Eq. 7). Extreme values were fixed for the exploration of the parameter space (Table 2). Values of parameters such as initial tissue porosities are already known to have an important effect on the poro-mechanical response of disk models and are unequivocally associated to disk degeneration. Hence, they were directly taken into account in the numerical optimizations without being previously explored through the DOE.

TABLE 2

Table 2. Range of values (high and low levels) for the parameters analyzed with the DOE sensitivity, and values used as components for the initial vector of parameters of the downhill simplex Python algorithm.

After having screened out the most important model parameters for both the short- and the long-term displacements, the specimen-specific optimizations were run by using the downhill simplex approach implemented in SciPy optimization toolbox (Oliphant, 2007). The initial vector to optimize was chosen based on literature data (Table 2). Since initial effective porosity ϕ_{0 d} was controlled through the damage parameter (Eq. 10), undamaged initial porosities ϕ₀ in Eq. 10 were fixed. According to biochemical measurements, values were respectively set to 0.75 and 0.8 for the AF and for the NP of grade 2 IVD, and to 0.7 and 0.75 for the AF and NP of disks with higher degeneration grades (Dao et al., 2014).

Optimizations were also run with no consideration of damage in order to assess the relevance of this parameter in the constitutive behavior. In this latter case, AF and NP initial porosities were optimized as independent parameters.

Validation Study

The reliability of the optimized grade-dependent material parameters was tested through an additional validation study for the L1-S1 specimen tested in vitro (Sp#4). The corresponding FE morphed model was loaded under unconstrained flexion, right lateral bending, and left axial rotation with moments of 5 Nm, simulating the experimental protocol detailed above. For each IVD, material properties were introduced in the FE model as a function of the Pfirrmann grade (see Table 4 in the Section “Result”). FE calculations were compared to the experiment in terms of both global RoM (rotation of L1 with respect to S1) and segmental RoM (e.g., rotation of L1 with respect to L2, etc.).

Mechano-Transport Simulations

Mechanical deformations determined changes in the amount of fluid/porosity and in the diffusion distances. A sequential mechanical-transport algorithm was used (Malandrino et al., 2011), in which the deformations from the mechanical simulation acted as template for all transport simulations. In the deformed geometry, diffusion of glucose, lactate, and oxygen were modeled through standard 3D diffusion-reaction:

$\frac{\partial }{\partial t}\left(\begin{array}{c}{C}_{{\text{O}}_2}\\ C_1\\ C_{\text{g}}\end{array}\right)-\left(\begin{array}{ccc}{D}_{{\text{O}}_2}& 0& 0\\ \text{0}& D_1& 0\\ \text{0}& 0& D_{\text{g}}\end{array}\right){\nabla }^2\left(\begin{array}{c}{C}_{{\text{O}}_2}\\ C_1\\ C_{\text{g}}\end{array}\right)=\left(\begin{array}{c}{R}_{{\text{O}}_2}\\ R_{\text{1}}\\ R_{\text{g}}\end{array}\right)$ (12)

where components C_i, D_i, and R_i are respectively the concentrations, the tissue-homogenized (effective) diffusion coefficients, and the reaction terms of oxygen (i = O₂), lactate (i = l), and glucose (i = g). Diffusion coefficients were strain-dependent and expressed as (Mackie and Meares, 1955; Gu et al., 2004):

$D_{\mathrm{i}}=-\frac{{\mathrm{\varphi }}^2}{{\left(2-\mathrm{\varphi }\right)}^2}{d}_i$ (13)

where d_i is the diffusivity of solute i in water, and ϕ is the current porosity, calculated from Eq. 4.

R_i terms in Eq. 12 stood for the interrelated cellular consumptions of oxygen and glucose and accumulation of lactate, proportional to the cell density ρ_cell at each point:

$R_{{\mathrm{O}}_2}=-\mathrm{\varphi }\frac{7.28{\mathrm{\rho }}_{\mathrm{cell}}}{S_{{\mathrm{O}}_2}}\left[\frac{C_{{\mathrm{O}}_{\mathrm{2}}}\left(\mathrm{pH}-\mathrm{4}\mathrm{.95}\right)}{1.46+C_{{\mathrm{O}}_{\mathrm{2}}}+4.03\left(\mathrm{pH}-\mathrm{4}\mathrm{.95}\right)}\right]$ (14)

$R_1={\text{ρ}}_{\text{cell}}{e}^{\left[-2.47+0.93\text{pH}+0.16C_{{\text{O}}_2}-0.0058C_{{\text{O}}_2}^2\right]}$ (15)

$R_{\mathrm{g}}=-\frac{1}{2}{R}_{\mathrm{l}}$ (16)

The pH was a continuous model variable, linearly related with the lactate concentration (Bibby et al., 2005; Soukane et al., 2007). Equations 14–16 were derived from in vitro tests on IVD cell cultures (Bibby et al., 2005). The three coupled diffusion-reaction models, defined through Eqs 12–16, were solved by the FE method through the implicit solver ABAQUS 6.12 (Simulia, Providence, RI, USA). Interrelations among Eqs 14–16 were taken into account through a finite time step approach (Malandrino et al., 2011, 2014b; Malandrino, 2012).

Exponential decays over time – starting from an initial cell density ρ_cell,0 characteristic of each tissue – simulated cell death, depending on the pH and on the glucose concentration within the tissue:

${\mathrm{\rho }}_{\mathrm{cell}}={\mathrm{\rho }}_{\mathrm{cell,0}}{e}^{-\left({\mathrm{\alpha }}_{\mathrm{g}}+{\mathrm{\alpha }}_{\mathrm{pH}}\right)t}$ (17)

where α_pH was −3.43 × 10⁻⁶ s⁻¹ if the pH was lower than 6.8 and 0 otherwise. Glucose-dependent cell death was implemented according to the cell death model proposed by Zhu et al. (2012):

${\mathrm{\alpha }}_{\mathrm{g}}=K\left(\frac{C_{\mathrm{g}}-0.5}{C_{\mathrm{g}}+k_1}-\frac{\left|C_{\mathrm{g}}-0.5\right|}{C_{\mathrm{g}}+k_2}\right)$ (18)

This formulation ensured that cell death always occurred below the critical concentration of glucose of 0.5 mM (Horner and Urban, 2001). Parameters K, k₁, and k₂ in Eq. 18 were taken directly from Zhu et al. (2012). The cell death model that further accounted for the contribution of α_pH (Eq. 17) was verified against its ability to reproduce temporal and spatial viability curves obtained in vitro, as a function of the cell densities characteristic of the NP and AF tissues, i.e., from 4 to 9 million cells/mL (Maroudas et al., 1975), and for the time periods used in the present study (0–3 days): a 26 mm width diffusion chamber – filled with cells embedded in 1% agarose gel – was simulated according to experimental tests on bovine nucleus cell viability (Horner and Urban, 2001). Oxygen, lactate, and glucose diffusivities accounted for gel porosity. Only half of a thin slice was reproduced due to chamber geometry. Initial conditions were applied throughout the chamber for the three metabolites to reproduce the experiment: an initial pH of 7.4 (i.e., an initial nil lactate concentration), an initial oxygen pressure of 21 kPa, and an initial glucose concentration of 5 mM. During simulations, these values were maintained only at the boundary since nutrient concentrations were maintained in the medium along the experiment (Horner and Urban, 2001).

For the mechano-transport simulations in the IVD models, a simplified daily load sequence was simulated, consisting in 8 h of 150 N compression (resting phase in a supine position) and 16 h of 750 N compression (physical activity phase). These levels of loads were chosen based on averages of activity and resting intradiscal pressures measured in vivo (Wilke et al., 1999). In order to avoid biases due to the transient initialization of the mechano-transport model, a preconditioning period of 3 days was simulated, both in terms of mechanical response and solute transport. Three additional days were subsequently simulated for cell viability calculations and for the analysis of the results. Solute boundary conditions were applied at the external CEP and AF edges with averaged values related to normal concentrations of metabolites in blood (Maroudas et al., 1975; Mokhbi Soukane et al., 2009).

The relevance of the particular IVD geometry combined to tissue mechanical competence was assessed by comparing the solute transport and cell viability results, respectively obtained with a grade 2 (Sp#4) and with a grade 4 (Sp#1) specimen-specific L3-L4 disk model. Also, a L5-S1 disk graded 4 from Sp#6 was analyzed. In all cases, grade-matching mechanical properties were used for both the AF and the NP, based on the constitutive parameter values extracted from the optimization previously described (see Results in Table 4). In order to assess further, the sole effect of mechanical competence, the grade 4 L3-L4 disk of Sp#1 was artificially modified into a “healthy” grade 2 disk in terms of material properties, and results were compared to those obtained with the grade-matching mechanical properties.

Results

Creep Experiments L3-L4

The behavior of the four specimens was different in absolute displacement, predominantly caused by the three preconditioning cycles (see Figure 2A). However, when looking at the height loss (0.87 mm) during the actual creep test, the outcome in displacement and time-dependent behavior was similar, except for Sp#6 (see Figure 2), which had a reduced initial height compared to the other disks (See dimensions reported in Table 5). All in all, the disk specimens reached about 20% of axial deformation with respect to their central height, except the large Sp#1 disk, the central part of which deformed by about 30%.

FIGURE 2

Figure 2. (A) Load history (left scale) and outcome of creep tests on the L3-L4 monosegments as measured in terms of axial displacements of the L3 vertebra (right scale). (B) Load history (left scale) and height loss when the effect of preconditioning was excluded (right scale).

Optimization of the Poro-Hyperelastic Material Properties with No Damage

Due its reduced height, Sp#6 disk model had extremely large element distortions and convergence issue impeded to achieve a successful optimization, independently of the constitutive material formulation. For the rest of disk specimen models, when damage was not considered (i.e., d = 0 in Eqs 8–10), the shear modulus increased with the degeneration grade. There were no other tendencies in terms of optimized parameters as a function of degeneration (Table 3).

TABLE 3

Table 3. No damage: optimized values as a function of the degeneration grades of the L3-L4 disks (see Table 1 and Figure 1 for specimen details) obtained through the minimization of the objective function (Eq. 11).

DOE Sensitivity Results

The DOE results indicated that among the seven parameters chosen, only CEP initial permeability could be excluded for further optimization, because its low significance on both the short- and the long-term displacement response (Figure 3). For the short-term simulations, the shear modulus, together with the NP and AF damage, and the dispersion parameter for AF fibers were the most important factors. For the long-term simulation, the AF dispersion parameter lost its significance, whereas the NP swelling pressure and the matrix undamaged bulk modulus had a substantially increased significance, though they were poorly significant in the short-term.

FIGURE 3

Figure 3. Pareto chart for the sensitivity of the displacement responses to the variation of the parameters chosen. (A) Short-term displacement and (B) long-term displacement responses. The threshold of significance (vertical line) for the ANOVA analyses is set to α = 0.01.

Optimization of the Poro-Hyperelastic Material Properties with Damage

Optimized undamaged shear and bulk moduli both increased with the degree of degeneration. Damage d increased with degeneration in the AF, and decreased in the NP. NP swelling decreased with degeneration (Table 4). NP damaged (or effective) shear modulus, calculated from Eq. 8 increased with degeneration while AF damaged shear modulus showed a decrease from grade 3 to grade 4 (Figure 4A), in agreement with the increase of AF damage. Both AF and NP, calculated with Eq. 9, increased with degeneration, although to a higher extent for the NP (Figure 4B). Initial damaged AF porosities slightly decreased with degeneration (Figure 4C) and NP swelling pressure switched from approximately 0.15 MPa to values reduced by about 40–50% (Table 4; Figure 4C).

FIGURE 4

Figure 4. Degeneration-dependent values of (A) the effective shear moduli, (B) the effective bulk moduli, and (C) the initial porosities and NP swelling pressures obtained after optimization.

TABLE 4

Table 4. Damage: optimized values as a function of the degeneration grades of the L3-L4 disks (see Table 1 and Figure 1 for specimen details) obtained through the minimization of the objective function (Eq. 11).

In Figure 5, the global performance of the set of material properties reported in Table 4 for the three disks optimized is shown. When compared to the experimental outcomes, the error was lower at the end of the long-term creep and higher at the end of the short-term creep (Figure 5). Figure 5 also shows that the importance of the osmotic pressure decreased from grade 2 to grade 4 (as also derivable from the optimized values of swelling pressure in Table 4), while the importance of the solid stress with respect to the fluid stress increased with the degeneration grade.

FIGURE 5

Figure 5. Comparisons between experiments and FE predictions in term of total vertical creep displacement for each of the optimized L3-L4 disks (top). Respective contributions of the osmotic, solid skeleton, and fluid stresses to the total stress in the central NP, predicted over the duration of the creep experiment for each of the optimized disks (down).

Validation Study

When incorporating specimen-specific distributions of the grade-dependent disk properties, the L1-S1 FE model morphed to the Sp#4 specimen geometry reproduced acceptably most of the in vitro flexion and right lateral bending RoM of the specimen, both globally and per segment (Figure 6). Under lateral bending, through the global RoM as well as the local RoM of three segments over five could be reproduced, the maximum rotations predicted for the L1-L2 grade 3 and for the L3-L4 grade 2 IVDs were about 40% larger than the experimental measures. As for the axial rotation load scenario, experimental results could not be reproduced by the FE predictions (data not shown, see Discussion).

FIGURE 6

Figure 6. Validation study by using the full L1-S1 FE model: (A) flexion results in terms of RoM (sagittal plane rotations in degrees) and comparison with the experimental response and (B) right lateral bending results in terms of RoM (frontal plane rotations in degrees) and comparison with the experimental response. All the curves shown are over time (transient response) and represent both the global (L1 vs. S1) and the segmental (e.g., L1 vs. L2, etc.) rotations.

Mechano-Transport Simulations

Figure 7 shows the viability profiles predicted by the diffusion chamber model in comparison to the related experimental data (Horner and Urban, 2001), for the two cell density levels relevant to the IVD, and at the specific time relevant to the present study (3 days of culture). In terms of distance to the solute source at which cells started to die, the computed cell viability profiles were in good agreement with the experimental data.

FIGURE 7

Figure 7. Cell viability profiles calculated in the simulated half-slice of the diffusion chamber and comparison with experimental results (Horner and Urban, 2001), for the two cell density levels relevant to the AF and NP tissues, i.e. (A) 4 million cells/mL and (B) 8 million cells/mL. Both numerical results and experimental data were considered at day 3 of culture.

When 6 days of mechano-transport (including the 3 days of preconditioning) were simulated, results in terms of glucose, lactate, oxygen, and cell viability depended strongly on both the patient-specific IVD geometry, and the degeneration-specific material properties. Only cell viability results are shown (Table 5). In the Sp#1 disk model, glucose levels dropped periodically below the critical concentration of 0.5 mM, and cells kept on dying over the 3 days of cell viability simulation, though progressively slower and slower along time (Figure 8). Surprisingly, when changing the material properties of Sp#1 disk model from grade 4 (original grade) to grade 2 (altered grade), the cell death area increased (Table 5).

TABLE 5

Table 5. Results of the mechano-transport simulations for different combinations of specimen-specific geometries and degeneration-specific tissue properties, and dimensions of the different disk geometries.

FIGURE 8

Figure 8. Glucose concentration and cell viability in the AF (boundary between the inner AF and the outer AF) for largest specimen studied, in terms of height (Table 5) (Sp#1, grade 4). Curves are shown over time for the last 3 days of mechano-transport simulation. Shaded areas correspond to the daily activity (16 h under 750 N compression) and non-shaded ones to the rest (8 h under 150 N compression).

Discussion

Validity of the FEM Subject-Specific Geometry

The correctness and accuracy of the complete subject-specific FEM models obtained after the morphing of the general model into the segmented geometries have been assessed by visual comparison and overlapping with the corresponding MRI images. In addition, the accuracy of the IVDs and the vertebrae segmentation was validated in Castro-Mateos et al. (2015), resulting in very good sub-voxel accuracies. Nevertheless, the craniocaudal concave shape of the anterior external AF of the Sp#6 disk model (Figure 1B) added to excessive element distortion affected the convergence and the capacity of the FE model to be successfully optimized. Such outcome might have arisen from a case-specific difficulty to define properly the AF external surface during the segmentation process.

Optimization of the Poro-Hyperelastic Material Properties with No Damage

Although the objective function was minimized and the inverse calculation successfully retrieved a set of parameters (Table 3), the behavior of these optimized parameters was not physically satisfactory. In particular, it has repeatedly been reported that porosity (tissue water content) of both AF and NP decreases with degeneration (Iatridis et al., 1997, 1998; Antoniou et al., 2013). This behavior was not observed for optimized parameters. Instead, very high and unrealistic porosities values were obtained for grade 3 and grade 4, probably indicative of numerical artifacts. Also, swelling pressure was at normal levels in the grade 4 IVD, in disagreement with the well known loss of proteoglycans with degeneration (Urban and Maroudas, 1981; Bibby et al., 2001). Other parameters, such as the CEP initial permeability and the AF fiber dispersion seemed to have a random effect on the optimization outcomes. These outcomes raised questions about the net effect and the physical significance of each constitutive parameter, and justified a DOE sensitivity analysis further performed with the additional inclusion of tissue damage.

DOE Sensitivity

Design of experiment sensitivity results pointed on the relative importance of the shear and bulk moduli in short- and long-term responses, respectively. Also, the significance of the damage parameter proved the phenomenological importance of this parameter, because of a direct relation to the effective elastic response of the tissue (see Eqs 8 and 9). Moreover, through the introduction of the possible presence of cracks, this approach paves the way for future studies targeting mechanistic-based description of degenerated and damaged disks. The DOE analyses also pointed out the dispersion parameter as an important factor in the anisotropic hyperelastic continuum model used (Gasser et al., 2006), under rapid loading. Though these sensitivity results correspond to a specific IVD, they were backed up by DOE analyses performed for other IVD specimen models and that included additional parameters, such as the AF permeability (data not shown).

Hydraulic permeability of the IVD tissues have been indicated as relevant parameters in other theoretical studies (Chagnon et al., 2010; Nikkhoo et al., 2013), while it has been shown that experimental AF permeability are not significantly affected by degeneration (Iatridis et al., 1998). Our sensitivity study prior to optimization confirmed these experimental outcomes and led us to exclude the CEP (and AF) permeability from the optimization studies. This choice was a compromise between the minimization of the parameter space and the accuracy of the results. However, not optimizing hydraulic permeability could be the cause of the lower optimization performance in short-term creep (Figure 5), and future optimization studies with an increased number of experimental constraints should overcome this limitation. Yet, having introduced a damage that relates with the porosity of a representative volume element (Eq. 10), also accounts indirectly for variations of the (effective) porosity-dependent hydraulic permeability (Eq. 5).

Optimization of the Poro-Hyperelastic Material Properties with Damage

Optimized poro-hyperelastic parameters combined with the concept of cracked continuum (damage) showed much more pronounced and meaningful trends than when damage was not considered (Table 3 vs. Table 4). Indeed, the increase of both the undamaged shear and bulk moduli correlates with reported fibrosis (Coventry et al., 1945; Iatridis et al., 1998; Pokharna and Phillips, 1998). Decrease of overall disk porosity and swelling pressure in the NP correlated with reported loss of fluid and drops in proteoglycans content. It could be argued that the number of specimens was limited and no statistical meaning can be attributed to the results. To this respect, further experimental–numerical explorations are planned for a next study. In the meantime, the high value of damage in the AF for severely degenerated disks (grade 4), shown in the present report, correlated very well with the presence of cracks often reported in MRI-based studies (Hilton et al., 1980; Munter et al., 2002; Sharma et al., 2009; Antoniou et al., 2013). These cracks in the AF would result in a net reduction of the stiffness of the disk. This reduction was captured by our description (Figure 4A) and has been reported in vitro comparing moderately (grade 3) and severely (grade 4) degenerated IVDs compressive responses (Panjabi et al., 1984; Gardner-Morse, 2004). In the NP, the nil damage prediction for grade 4 disks is obviously difficult to defend: it could indicate a numerical artifact related with the minor importance of the effective stiffness of the NP (compared to that of the AF) in controlling the global IVD vertical displacement (Malandrino et al., 2009) in severely degenerated IVDs. Also, fissures in the NP are less meaningful than those in the AF due to the relatively minor role of the NP shear resistance. As for the components that affect the effective volumetric stiffness of the NP, the model acceptably captured the reduced swelling pressure and hydration in the optimized grade 4 disks. Note that variations between grade 1 and grade 2 need to be interpreted carefully, since grade 1 properties were subjectively chosen among the data available from the literature.

All in all, optimized models reproduced nicely the creep behavior of the mono-segment specimen, in a general way. Analyzing the detail of the results, while models reproduced adequately the experimental axial displacements in the long-term, the constitutive description did not allow reproducing accurately the short-term creep experimental behavior. This could be ascribed to the absence of any explicit viscoelastic description for the collagen fibers. Nevertheless, this poro-hyper-viscoelastic description (Wilson et al., 2006) would have increased the total number of parameters, and it will be the subject of further evaluations in future studies. The solid phase contribution to total stress was shown to increase considerably when degeneration passed from grade 3 to grade 4. According to the known degeneration changes, our result illustrate increased fibrosis in the tissue, and importantly, a decreased capability of the fluid phase in the NP in bearing the load in grade 4 disks, as highlighted by other poro-elastic parametric studies (Chagnon et al., 2010).

When introduced in a whole L1-S1 FE model, the optimized degeneration-specific parameters of the disks, were able to reproduce most of the experimental responses of the specimen in terms of frontal and sagittal RoM. Notably, predictions were very accurate over time, at the global level (Figure 6). At the local level, two disk models resulted over-compliant compared to the measurements in lateral bending. Also, the FE predictions in terms of RoMs under axial rotation poorly compared with the experimental ones (results not shown). In the specific case of axial rotation, the following issues could have affected the outcomes of the comparisons: (i) the measured segmental ranges of motion were very low, hence indicating large probabilities of measurement errors; (ii) accordingly, the sum of the different intersegmental axial rotations did not correspond to the global amount of axial rotation (while it did for flexion and lateral bending loading scenarios). Thus, additional axial rotation experiments would be necessary to clarify this. But more generally, both axial rotation and lateral bending involve important coupled rotations that largely depend on the regional properties of the IVD, especially of the AF (Noailly, 2009; Noailly et al., 2011). Hence, local refinements of the disk properties might be planned based on a second step of optimizations that would use the rotation experimental data. According to the results of the present study, and to previous findings in spine model calibration (Schmidt et al., 2006), the optimized parameters should ensure a suitable model preconditioning for such a second optimization.

Mechano-Transport Simulations

Mechano-transport simulations showed that disk morphology was related to the onset of cell viability. In particular, disk height controlled the diffusion distances from the source of nutrients at the endplates. This distance-related cell death appeared as the most important factor in the present study, since it could explain by itself cell death in the IVD reconstructions tested with patient-specific geometry and different sets of material properties. Nevertheless, material properties influenced the mechano-transport results as well. For example, for the specific geometry of the Sp#1 model, initially grade 4, the introduction of grade 2 properties led to an increased consolidation, i.e., increased water loss and decreased solute diffusion (Eq. 13), which exacerbated cell death occurrence. Truly, associating grade 2 properties to a grade 4 disk that might already have reduced disk height could be expected to bias the interpretation of the calculations. Hence, we compared the consolidation of disks with large (Sp#1) and small (Sp#6, Sp#4) heights, for a similar set of material parameters and similar boundary loads. We found that the larger disk led to only slightly increased tissue consolidation compare to a small disk (small in terms of disk height), i.e., 2–3% of loss of volumetric fluid fraction in the critical zone of cell death during a daily cycle. According to Eq. 13, this loss of porosity would only slightly affect the diffusivity. Thus, the consolidation-dependent diffusion appeared as a minor effect compared to the large diffusion distance effect.

In previous mechano-transport models of the IVD, mechanical loads representative of daily activity induced fluctuations of glucose concentration and cell viability but did not cause cell death per se in disk models with a height of about 10 mm (Zhu et al., 2012; Malandrino et al., 2014b). These studies concluded on the paramount importance of endplate-related issues, such as calcifications or decrease in vascular according to other reported simulation outcomes (Galbusera et al., 2011a; Zhu et al., 2012; Malandrino et al., 2014a). Only in such cases, increased tissue dehydration/consolidation was reported to be an additional risk factor for cell death, related to tissue degeneration (Zhu et al., 2012; Malandrino et al., 2014b). In the present study, disk size-related diffusion distance emerges as an additional relevant factor to nutrition-induced cell viability, independently on tissue degeneration. This outcome is nicely supported by experimental measurements that negatively correlated NP cell density with disk height (Rodriguez et al., 2011), and naturally follows the conclusions of recent numerical explorations of the effect of disk size on the transport of externally administered drugs (Motaghinasab et al., 2014). These concerns can now be targeted patient-specifically from diagnostic images, and could highly impact the prognosis of DDD, including explanations of gender-related risk factors (males may have larger disks, in general).

Significance

On one hand, this study assessed the capacity of commonly used poro-mechanical disk tissue models to simulate both healthy and degenerated tissues. We found that the basic (osmo)poro-hyperelastic constitutive formulation was not able to capture simultaneously the respective effects of tissue dehydration, fibrosis, and crack occurrence. However, the introduction of a damage criterion that affected in a coordinated way the shear and the bulk stiffness, and the porosity, based on micromechanics rationales, allowed both successful simulation of known degenerated effects, and validation of the disk organ models. Coupling this new mechanical formulation to different specimen-specific IVD models confirmed the previously reported idea that AF and NP tissue degeneration alone might not cause nutrition-induced cell death, especially if cracks occur and can be filled with water. However, simulations suggested for the first time that large IVDs with healthy tissues might be particularly prone to degeneration because of cell nutrition issues mostly related to increased diffusion distances. This outcome strongly suggests that early degeneration might be related to disk morphology, i.e., a genetic characteristic unavoidably coupled to the spine mechanics.

Conflict of Interest Statement

The Associate Editor Fabio Galbusera declares that despite being affiliated with the same institution as the author, Hans-Joachim Wilke, the review process was handled objectively and no conflict of interest exists. 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.

Acknowledgments

European Commission (My SPINE 269909) is acknowledged for funding. The Authors are grateful to Dr. Fabrice Charleux from ACRIM-Polyclinique Saint Côme, Compiègne, France, for the MRI acquisitions and IVD grading. David Volkheimer, from the Institute of Orthopedic Research and Biomechanics, Ulm, Germany, is also acknowledged for his help with the mechanical tests.

References

Antoniou, J., Epure, L. M., Michalek, A. J., Grant, M. P., Iatridis, J. C., and Mwale, F. (2013). Analysis of quantitative magnetic resonance imaging and biomechanical parameters on human discs with different grades of degeneration. J. Magn. Reson. Imaging 38, 1402–1414. doi: 10.1002/jmri.24120

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Antoniou, J., Steffen, T., Nelson, F., Winterbottom, N., Hollander, A. P., Poole, R. A., et al. (1996). The human lumbar intervertebral disc: evidence for changes in the biosynthesis and denaturation of the extracellular matrix with growth, maturation, ageing, and degeneration. J. Clin. Invest. 98, 996–1003. doi:10.1172/JCI118884

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Bibby, S. R. S., Jones, D. A., Lee, R. B., Yu, J., and Urban, J. P. G. (2001). Pathophysiology of the intervertebral disc. J. Bone Spine 68, 537–542. doi:10.1002/jmri.20874

CrossRef Full Text | Google Scholar

Bibby, S. R. S., Jones, D. A., Ripley, R. M., and Urban, J. P. G. (2005). Metabolism of the intervertebral disc: effects of low levels of oxygen, glucose, and pH on rates of energy metabolism of bovine nucleus pulposus cells. Spine 30, 487–496. doi:10.1097/01.brs.0000154619.38122.47

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Budiansky, B., and O’Connell, R. (1976). Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81–97. doi:10.1016/j.jmbbm.2012.07.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Castro-Mateos, I., Jose, M. P., Peter, E. E., Del Rio, L., Lazary, A., and Frangi, A. F. (2014). 3D segmentation of annulus fibrosus and nucleus pulposus from T2-weighted magnetic resonance images. Phys. Med. Biol. 59, 7847–7864.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Castro-Mateos, I., Pozo, J. M., Lazary, A., and Frangi, A. F. (2015). 3D Vertebra segmentation by feature selection active shape model. Lecture Notes in Computational Vision and Biomechanics, Vol. 20.

Google Scholar

Chagnon, A., Aubin, C.-E., and Villemure, I. (2010). Biomechanical influence of disk properties on the load transfer of healthy and degenerated disks using a poroelastic finite element model. J. Biomech. Eng. 132, 111006. doi:10.1115/1.4002550

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Coventry, M. B., Ghormley, R. K., and Kernohan, J. W. (1945). The intervertebral disc: its microscopic anatomy and pathology part III. Pathological changes in the intervertebral disc. J. Bone Joint Surg. 27, 460–474.

Google Scholar

Dao, T. T., van Rijsbergen, M., Pouletaut, P., Charleux, F., Ito, K., and Ho Ba Tho, M. C. (2014). “In vitro assessment of micro-structural properties of intervertebral disc using 1.5T magnetic resonance T2 and ADC mappings,” in Proceedings of the 7th World Congress of Biomechanics, July 6-11, 2014, Boston, MA.

Google Scholar

Dormieux, L., and Kondo, D. (2007). Micromechanics of damage propagation in fluid-saturated cracked media. Revue Européenne de Génie Civil, 11, 37–41. doi:10.3166/REGC.11.945-962

CrossRef Full Text | Google Scholar

Eberlein, R., Holzapfel, G., and Schulze-Bauer, C. (2001). An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies. Comput. Methods Biomech. Biomed. Engin. 4, 209–229. doi:10.1080/10255840108908005

CrossRef Full Text | Google Scholar

Ferguson, S. J., Ito, K., and Nolte, L.-P. (2004). Fluid flow and convective transport of solutes within the intervertebral disc. J. Biomech. 37, 213–221. doi:10.1016/S0021-9290(03)00250-1

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Frobin, W., Brinckmann, P., Biggemann, M., Tillotson, M., and Burton, K. (1997). Precision measurement of disc height, vertebral height and sagittal plane displacement from lateral radiographic views of the lumbar spine. Clin. Biomech. 12, S1–S63. doi:10.1016/S0268-0033(96)00067-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Frobin, W., Brinckmann, P., Kramer, M., and Hartwig, E. (2001). Height of lumbar discs measured from radiographs compared with degeneration and height classified from MR images. Eur. Radiol. 11, 263–269. doi:10.1007/s003300000556

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Fujita, Y., Duncan, N., and Lotz, J. (1997). Radial tensile properties of the lumbar annulus fibrosus are site and degeneration dependent. J. Orthop. Res. 15, 814–819. doi:10.1002/jor.1100150605

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Galbusera, F., Mietsch, A., Schmidt, H., Wilke, H.-J., and Neidlinger-Wilke, C. (2011a). Effect of intervertebral disc degeneration on disc cell viability: a numerical investigation. Comput. Methods Biomech. Biomed. Engin. 16, 37–41. doi:10.1080/10255842.2011.619184

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Galbusera, F., Schmidt, H., Noailly, J., Malandrino, A., Lacroix, D., Wilke, H.-J., et al. (2011b). Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs. J. Mech. Behav. Biomed. Mater. 4, 1234–1241. doi:10.1016/j.jmbbm.2011.04.008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gardner-Morse, M. (2004). Structural behavior of human lumbar spinal motion segments. J. Biomech. 37, 205–212. doi:10.1016/j.jbiomech.2003.10.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gasser, T. C., Ogden, R. W., and Holzapfel, G. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface 3, 15–35. doi:10.1098/rsif.2005.0073

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Gu, W. Y., Yao, H., Vega, A. L., and Flagler, D. (2004). Diffusivity of ions in agarose gels and intervertebral disc: effect of porosity. Ann. Biomed. Eng. 32, 1710–1717. doi:10.1007/s10439-004-7823-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Halloran, J. P., Sibole, S., Van Donkelaar, C. C., Van Turnhout, M. C., Oomens, C. W. J., Weiss, J. A., et al. (2012). Multiscale mechanics of articular cartilage: potentials and challenges of coupling musculoskeletal, joint, and microscale computational models. Ann. Biomed. Eng. 40, 2456–2474. doi:10.1007/s10439-012-0598-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Henak, C. R., Anderson, A. E., and Weiss, J. A. (2013). Subject-specific analysis of joint contact mechanics: application to the study of osteoarthritis and surgical planning. J. Biomech. Eng. 135, 021003. doi:10.1115/1.4023386

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Hilton, R., Ball, J., and Benn, R. (1980). Annular tears in the dorsolumbar spine. Ann. Rheum. Dis. 39, 533–538. doi:10.1136/ard.39.6.533

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Horner, H. A., and Urban, J. P. (2001). 2001 Volvo Award Winner in Basic Science Studies: effect of nutrient supply on the viability of cells from the nucleus pulposus of the intervertebral disc. Spine 26, 2543–2549. doi:10.1097/00007632-200112010-00006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Huang, C.-Y., Travascio, F., and Gu, W. Y. (2012). Quantitative analysis of exogenous IGF-1 administration of intervertebral disc through intradiscal injection. J. Biomech. 45, 1149–1155. doi:10.1016/j.jbiomech.2012.02.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Huang, Y.-C., Urban, J. P. G., and Luk, K. D. K. (2014). Intervertebral disc regeneration: do nutrients lead the way? Nat. Rev. Rheumatol. 10, 561–566. doi:10.1038/nrrheum.2014.91

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Iatridis, J. C., MacLean, J. J., Roughley, P. J., and Alini, M. (2006). Effects of mechanical loading on intervertebral disc metabolism in vivo. J. Bone Joint Surg. Am. 88(Suppl. 2), 41–46. doi:10.2106/JBJS.E.01407

Pubmed Abstract | Pubmed Full Text | 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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Iatridis, J. C., Setton, L. A., Weidenbaum, M., and Mow, V. C. (1997). Alterations in the mechanical behavior of the human lumbar nucleus pulposus with degeneration and aging. J. Orthop. Res. 15, 318–322. doi:10.1002/jor.1100150224

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Johannessen, W., and Elliott, D. M. (2005). Effects of degeneration on the biphasic material properties of human nucleus pulposus in confined compression. Spine 30, E724–E729. doi:10.1097/01.brs.0000192236.92867.15

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Mackie, J. S., and Meares, P. (1955). The diffusion of electrolytes in a cation-exchange resin membrane. I. Theoretical. Proc. R. Soc. London. Ser. A Math. Phys. Sci. 232, 498–509. doi:10.1098/rspa.1955.0234

CrossRef Full Text | Google Scholar

Malandrino, A. (2012). Multi-Scale Biomechanical Study of Transport Phenomena in the Intervertebral Disc. PhD thesis, Universitat Politècnica de Catalunya, Barcelona.

Google Scholar

Malandrino, A., Lacroix, D., Hellmich, C., Ito, K., Ferguson, S. J., and Noailly, J. (2014a). The role of endplate poromechanical properties on the nutrient availability in the intervertebral disc. Osteoarthr. Cartil. 22, 1053–1060. doi:10.1016/j.joca.2014.05.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Malandrino, A., Noailly, J., and Lacroix, D. (2014b). Numerical exploration of the combined effect of nutrient supply, tissue condition and deformation in the intervertebral disc. J. Biomech. 47, 1520–1525. doi:10.1016/j.jbiomech.2014.02.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Malandrino, A., Lacroix, D., and Noailly, J. (2013). “Intervertebral disc cell death explained by metabolism-deformation couplings in a porohyperelastic finite element model,” in Poromechanics V: Proceedings of the Fifth Biot Conference on Poromechanics, (Reston, VA: American Society of Civil Engineers), 2193–2201.

Google Scholar

Malandrino, A., Noailly, J., and Lacroix, D. (2011). The effect of sustained compression on oxygen metabolic transport in the intervertebral disc decreases with degenerative changes. PLoS Comput. Biol. 7:e1002112. doi:10.1371/journal.pcbi.1002112

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Malandrino, A., Planell, J. A., and Lacroix, D. (2009). Statistical factorial analysis on the poroelastic material properties sensitivity of the lumbar intervertebral disc under compression, flexion and axial rotation. J. Biomech. 42, 2780–2788. doi:10.1016/j.jbiomech.2009.07.039

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Maroudas, A., Stockwell, R. A., Nachemson, A., and Urban, J. (1975). Factors involved in the nutrition of the human lumbar intervertebral disc: cellularity and diffusion of glucose in vitro. J. Anat. 120, 113–130.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Mokhbi Soukane, D., Shirazi-Adl, A., and Urban, J. P. G. (2009). Investigation of solute concentrations in a 3D model of intervertebral disc. Eur. Spine J. 18, 254–262. doi:10.1007/s00586-008-0822-7

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Montgomery, D. C. (2001). Design and analysis of experiments. Design 48, 158–158. doi:10.1198/tech.2006.s372

CrossRef Full Text | Google Scholar

Mori, T., and Tanaka, K. (1973). Average stress in the matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 21, 571–574. doi:10.1016/0001-6160(73)90064-3

CrossRef Full Text | Google Scholar

Motaghinasab, S., Shirazi-Adl, A., Parnianpour, M., and Urban, J. P. G. (2014). Disc size markedly influences concentration profiles of intravenously administered solutes in the intervertebral disc: a computational study on glucosamine as a model solute. Eur. Spine J. 23, 715–723. doi:10.1007/s00586-013-3142-5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Munter, F. M., Wasserman, B. A., Wu, H.-M., and Yousem, D. M. (2002). Serial MR imaging of annular tears in lumbar intervertebral disks. AJNR Am. J. Neuroradiol. 23, 1105–1109.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Neidlinger-Wilke, C., Mietsch, A., Rinkler, C., Wilke, H. J., Ignatius, A., and Urban, J. (2012). Interactions of environmental conditions and mechanical loads have influence on matrix turnover by nucleus pulposus cells. J. Orthop. Res. 30, 112–121. doi:10.1002/jor.21481

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Neubert, A., Fripp, J., Engstrom, C., Gal, Y., Crozier, S., and Kingsley, M. I. C. (2014). Validity and reliability of computerized measurement of lumbar intervertebral disc height and volume from magnetic resonance images. Spine J. 14, 2773–2781. doi:10.1016/j.spinee.2014.05.023

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Nikkhoo, M., Hsu, Y.-C., Haghpanahi, M., Parnianpour, M., and Wang, J.-L. (2013). A meta-model analysis of a finite element simulation for defining poroelastic properties of intervertebral discs. Proc. Inst. Mech. Eng. H 227, 672–682. doi:10.1177/0954411913480668

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Noailly, J. (2009). Model developments for in silico studies of the lumbar spine biomechanics. Ph.D. thesis, Technical University of Catalonia (UPC), Barcelona.

Google Scholar

Noailly, J., Ambrosio, L., Elizabeth Tanner, K., Planell, J. A., and Lacroix, D. (2012). In silico evaluation of a new composite disc substitute with a L3-L5 lumbar spine finite element model. Eur. Spine J. 21(Suppl. 5), S675–S687. doi:10.1007/s00586-011-1716-7

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Noailly, J., and Lacroix, D. (2012). “Finite element modelling of the spine,” in Biomaterials for Spinal Surgery – Part I: Fundamentals of Biomaterials for Spinal Surgery, eds L. Ambrosio and K. E. Tanner (Cambridge: Woodhead Publishing Ltd), 144–232.

Google Scholar

Noailly, J., Malandrino, A., and Galbusera, F. (2014). “Computational modelling of spinal implants,” in Computational Modelling of Biomechanics and Biotribology in the Musculoskeletal System, (Cambridge: Woodhead Publishing Ltd), 447–484.

Google Scholar

Noailly, J., Planell, J. A., and Lacroix, D. (2011). On the collagen criss-cross angles in the annuli fibrosi of lumbar spine finite element models. Biomech. Model. Mechanobiol. 10, 203–219. doi:10.1007/s10237-010-0227-5

Pubmed Abstract | Pubmed Full Text | 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 | Pubmed Full Text | CrossRef Full Text | Google Scholar

Oliphant, T. E. (2007). SciPy: open source scientific tools for Python. Comput. Sci. Eng. 9, 10–20. doi:10.1109/MCSE.2007.58

CrossRef Full Text | Google Scholar

Panjabi, M. M., Krag, M. H., and Chung, T. Q. (1984). Effects of disc injury on mechanical behavior of the human spine. Spine 9, 707–713. doi:10.1097/00007632-198410000-00010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Périé, D., Korda, D., and Iatridis, J. C. (2005). Confined compression experiments on bovine nucleus pulposus and annulus fibrosus: sensitivity of the experiment in the determination of compressive modulus and hydraulic permeability. J. Biomech. 38, 2164–2171. doi:10.1016/j.jbiomech.2004.10.002

CrossRef Full Text | Google Scholar

Pfirrmann, C. W., Metzdorf, A., Zanetti, M., Hodler, J., and Boos, N. (2001). Magnetic resonance classification of lumbar intervertebral disc degeneration. Spine 26, 1873–1878. doi:10.1097/00007632-200109010-00011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Pokharna, H. K., and Phillips, F. M. (1998). Collagen crosslinks in human lumbar intervertebral disc aging. Spine 23, 1645–1648. doi:10.1097/00007632-199808010-00005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Rinkler, C., Heuer, F., Pedro, M. T., Mauer, U. M., Ignatius, A., and Neidlinger-Wilke, C. (2010). Influence of low glucose supply on the regulation of gene expression by nucleus pulposus cells and their responsiveness to mechanical loading. J. Neurosurg. Spine 13, 535–542. doi:10.3171/2010.4.SPINE09713

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Rodriguez, A. G., Slichter, C. K., Acosta, F. L., Rodriguez-Soto, A. E., Burghardt, A. J., Majumdar, S., et al. (2011). Human disc nucleus properties and vertebral endplate permeability. Spine 36, 512–520. doi:10.1097/BRS.0b013e3181f72b94

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Ruiz, C., Noailly, J., and Lacroix, D. (2013). Material property discontinuities in intervertebral disc porohyperelastic finite element models generate numerical instabilities due to volumetric strain variations. J. Mech. Behav. Biomed. Mater. 26, 1–10. doi:10.1016/j.jmbbm.2013.05.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Schmidt, H., Heuer, F., Simon, U., Kettler, A., Rohlmann, A., Claes, L., et al. (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. Clin. Biomech. 21, 337–44. doi:10.1016/j.clinbiomech.2005.12.001

CrossRef Full Text | Google Scholar

Schmidt, H., Midderhoff, S., Adkins, K., and Wilke, H.-J. (2009). The effect of different design concepts in lumbar total disc arthroplasty on the range of motion, facet joint forces and instantaneous center of rotation of a L4-5 segment. Eur. Spine J. 18, 1695–1705. doi:10.1007/s00586-009-1146-y

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Schroeder, Y., Huyghe, J. M., van Donkelaar, C. C., and Ito, K. (2010). A biochemical/biophysical 3D FE intervertebral disc model. Biomech. Model Mechanobiol. 9, 641–650. doi:10.1007/s10237-010-0203-0

CrossRef Full Text | Google Scholar

Sharma, A., Pilgram, T., and Wippold, F. J. (2009). Association between annular tears and disk degeneration: a longitudinal study. AJNR Am. J. Neuroradiol. 30, 500–506. doi:10.3174/ajnr.A1411

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Smith, L. J., Nerurkar, N. L., Choi, K.-S., Harfe, B. D., and Elliott, D. M. (2011). Degeneration and regeneration of the intervertebral disc: lessons from development. Dis. Model. Mech. 4, 31–41. doi:10.1242/dmm.006403

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Soukane, D. M., Shirazi-Adl, A., Urban, J. P. G. (2007). Computation of coupled diffusion of oxygen, glucose and lactic acid in an intervertebral disc. J. Biomech. 40, 2645–54. doi:10.1016/j.jbiomech.2007.01.003

CrossRef Full Text | Google Scholar

Travascio, F., Elmasry, S., and Asfour, S. (2014). Modeling the role of IGF-1 on extracellular matrix biosynthesis and cellularity in intervertebral disc. J. Biomech. 47, 2269–2276. doi:10.1016/j.jbiomech.2014.04.046

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Tunset, A., Kjaer, P., Samir Chreiteh, S., and Secher Jensen, T. (2013). A method for quantitative measurement of lumbar intervertebral disc structures: an intra- and inter-rater agreement and reliability study. Chiropr. Man. Therap. 21, 26. doi:10.1186/2045-709X-21-26

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Urban, J. P. G., and Maroudas, A. (1981). Swelling of the intervertebral disc in vitro. Connect. Tissue Res. 9, 1–10. doi:10.3109/03008208109160234

CrossRef Full Text | Google Scholar

Wilke, H., Claes, L., Schmitt, H., and Wolf, S. (1994). A universal spine tester for in vitro experiments with muscle force simulation. Eur. Spine J. 3, 91–97. doi:10.1007/BF02221446

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wilke, H. J., Neef, P., Caimi, M., Hoogland, T., and Claes, L. E. (1999). New in vivo measurements of pressures in the intervertebral disc in daily life. Spine 24, 755–762. doi:10.1097/00007632-199904150-00005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wilson, W., Huyghe, J. M., and van Donkelaar, C. C. (2006). A composition-based cartilage model for the assessment of compositional changes during cartilage damage and adaptation. Osteoarthr. Cartil. 14, 554–560. doi:10.1016/j.joca.2005.12.006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Wilson, W., van Donkelaar, C. C., and Huyghe, J. M. (2005). A comparison between mechano-electrochemical and biphasic swelling theories for soft hydrated tissues. J. Biomech. Eng. 127, 158–165. doi:10.1115/1.1835361

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Zhu, Q., Jackson, A. R., and Gu, W. Y. (2012). Cell viability in intervertebral disc under various nutritional and dynamic loading conditions: 3d Finite element analysis. J. Biomech 45, 2769–2777. doi:10.1016/j.jbiomech.2012.08.044

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar