Porohyperelastic modeling of intervertebral disc tissues

Disc tissues are often studied as deformable continuous media. Provided a reasonably large scale is chosen, differential calculus and material properties averaging become possible for the mathematical treatment of a material point as a homogenized mix of lower scale components and structures. A simplification of IVD material considers: (1) A solid phase made up of structural macromolecules (collagen, elastin, and proteoglycans) as well as cells, and (2) a fluid phase that consists of water and solutes. Following this approach, mixture or biphasic theory^{Reference Mow, Kuei, Lai and Armstrong10} considers that two particles (i.e., solid and fluid) occupy each material point, whereas the Biot poroelastic theory^{Reference Biot11} treats each material point as a homogenized combination of the two phases. However, both theories are equivalent when incompressible fluid and solid constituents are assumed.^{Reference Simon, Wu, Carlton, Evans and Kazarian12,Reference Huyghe and Janssen13}

In poroelastic IVD models, the apparent solid stiffness of the tissues is that of a compressible, drained porous skeleton, and it largely influences disc mechanics predictions.^{Reference Malandrino, Planell and Lacroix14} An improvement of the poroeleastic theory needs to account for geometrical nonlinearity effects under large deformations.^{Reference Wagner and Lotz15–Reference Malandrino, Noailly and Lacroix17} Known as porohyperelasticity,^{Reference Simon, Kaufmann, McAfee, Baldwin and Wilson18} the resulting theory derives this stiffness through a strain energy density function, the simplest form of which is the Neo-Hookean model (i.e., the Helmholtz free energy of a molecular network with Gaussian chain-length distribution):^{Reference Treloar19}

(1)$${W_{{\rm{ISO}}}} = {G \over 2}\left( {{I_1} - 3} \right) + {K \over 2}{\left( {J - 1} \right)^2}.$$

W _ISO is the isotropic strain energy density; G and K are the shear and bulk modulus of the porous solid; I ₁ is the first invariant of the Left Cauchy-Green deformation tensor of the continuum; and J is the determinant of the deformation gradient tensor. This strain energy density has been used to model the behavior of NP, CEP, and AF.^{Reference Malandrino, Noailly and Lacroix17,Reference Jacobs, Cortes, Peloquin, Vresilovic and Elliott20}

In AF, the concentric lamellae present specific fiber orientation patterns,^{Reference Zhu, Gu, Wu, Gong, Zhu, Jiang and Cao21,Reference Holzapfel, Schulze-Bauer, Feigl and Regitnig22} the mechanical importance of which has been highlighted in both experimental studies^{Reference Ebara, Iatridis, Setton, Foster, Mow and Weidenbaum23} and simulations.^{Reference Malandrino, Noailly and Lacroix17,Reference Noailly, Planell and Lacroix24} Hence, an anisotropic strain energy density term, W _ANISO, has often been added to W _ISO (Equation 1), W _hyper = W _ANISO+W _ISO. Though W _ANISO can account for fiber-matrix interactions,^{Reference Caner, Guo, Moran, Bazant and Carol25,Reference Reutlinger, Bürki, Brandejsky, Ebert and Büchler26} one of its simplest forms applied to the AF^{Reference Eberlein, Holzapfel and Schulze-Bauer27,Reference Eberlein, Holzapfel and Fröhlich28} considers that the hyperelastic fibers are perfectly bound to the matrix, as initially proposed by Holzapfel, Gasser, and Ogden.^{Reference Gasser, Ogden and Holzapfel29,Reference Holzapfel, Gasser and Ogden30} This approximation allowed the local AF mechanics to be reproduced with fiber angle spatial distributions in line with reported anatomical descriptions.^{Reference Malandrino, Noailly and Lacroix17} The simulation of viscoelastic fibers (e.g., as done by Schroeder and colleagues)^{Reference Schroeder, Wilson, Huyghe and Baaijens31,Reference Schroeder, Huyghe, van Donkelaar and Ito32} has been deemed necessary in order to properly simulate the short-term response of IVD specimens under compression.^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33}

In the porohyperelastic continuum, the fluid phase contribution to the total stress tensor affects the volumetric part and is quantified by the pore pressure, p, so that the total stress reads as:

(2)$${\rm{\sigma }} = {2 \over J}{{\partial W_{{\rm{hyper}}} } \over {\partial b}}b - pI,$$

where b is the Left Cauchy-Green deformation tensor. W _hyper can be represented by, for example, Equation 1 in the isotropic case. In the AF, it would be complemented by the previously mentioned additional term W _ANISO, considering anisotropic energy contributions stemming from the oriented fibers.

In order to take into account the osmotic effects generated by the proteoglycan fixed charge densities in disc tissues (osmo-porohyperelasticity), p has been proposed to be the sum of an osmotic pressure, Δπ, and a water chemical potential,^{Reference Schroeder, Wilson, Huyghe and Baaijens31} μ^f:

(3)$$p = \Delta {\rm{\pi }} + {{\rm{\mu }}^{\rm{f}}}.$$

By convention and according to Equation 2, p is positive in compression. The techniques proposed to simulate Δπ are described in the following sections of the present article. As for μ^f, its action is weighted by the continuum porosity ϕ, according to Darcy’s law, which relates the spatial gradients of μ^f to the fluid velocity relative to the solid phase, u_f:

(4)$$u_f {\rm{\phi }} = k \cdot \nabla {\rm{\mu }}^{\rm{f}} ,$$

where k is the hydraulic permeability tensor, often defined via a porosity-dependent function based on in vitro measurements.^{Reference Holmes and Mow34}

If both the fluid and solid phases are assumed to be incompressible, the porosity varies with J from an initial (undeformed) porosity value ϕ:

(5)$${\rm{\phi }} = 1 - J^{ - 1} \left( {1 - {\rm{\phi }}_{\rm{0}} } \right).$$

Under mechanical loads, the local porosity of the continuum depends on the ability of the porous solid skeleton to deform volumetrically (e.g., on K in Equation 1, and on the capacity of water to flow out of the tissue, as described by Equation 4).

Electrochemomechanical models including osmotic effects

Donnan theory

The swelling of the proteoglycan gel and the stiffness of the collagen structure ensure that the structure of the disc is prestressed.^{Reference Galbusera, van Rijsbergen, Ito, Huyghe, Brayda-Bruno and Wilke35} While the swollen gel is under compression, the collagen fibers are under tension, resulting in an IVD compressive strength of more than 9000 N. The mechanism of swelling of the proteoglycan gel is achieved by electrostatic charges along the macromolecules. Indeed, the negative fixed charge of the proteoglycans attracts positive counter ions, typically sodium, which in turn attracts water. This attraction is quantified in the form of an osmotic pressure, which according to the classical van ’t Hoff equation is proportional to the difference in ionic concentration in the disc and the peripheral vasculature:

(6)$$\Delta {\rm{\pi }} = RT\left( {{c_ + } + {c_ - }} \right),$$

where R and T are the universal gas constant and absolute temperature, respectively. Since soluble salts diffuse rapidly, the cationic concentration c ₊ and the anionic concentration c _– are often assumed to be in equilibrium with the salt concentration c within the peripheral vasculature, in terms of osmolarity:

(7)$$c_ + + c_ - = \sqrt {c_{{\rm{fc}}}^2 + 4c^2 } .$$

In the above equation, c _fc represents the fixed charge density of the tissue per unit fluid volume, which typically depends on deformation and porosity (Equation 5):

(8)$$c_{{\rm{fc}}} = {{{\rm{\phi }}_0 c_{{\rm{fc}}0} } \over {J{\rm{\phi }}}},$$

where c _fc0 is the initial fixed charge density per unit of fluid volume. The above equation demonstrates that the fixed charge density increases as more fluid is expelled from the disc, and the increasing osmotic pressure (Equation 6) adds to the effective disc stress (Equations 2 and 3). The immunity of the disc against tearing is believed to result from this osmotic prestressing, while disc herniation may be triggered by loss of Donnan osmotic stress.^{Reference Wognum, Huyghe and Baaijens36}

Triphasic mixture theory

As an extension of the biphasic model, the triphasic mixture theory has been used to describe the electrochemomechanical behavior of the IVD.^{Reference Lai, Hou and Mow42–Reference Yao and Gu44} In addition to the solid and fluid mixture, a solute phase (nutrients, ions) was considered.^{Reference Lai, Hou and Mow42} This theory no longer assumes electrochemical equilibrium of the whole disc with the vascular osmolarity, as expressed in Equation 6. This equilibrium is restricted to the outer boundary of the disc, which requires explicit modeling of the ionic diffusion-convection within the disc. Hence, governing equations were based on the balance of linear momentum for the mixture and on the conservation of mass for each of the phase or solute species.

The total stress of the mixture can be expressed as given by Equation 2. As for the fluid (J^f) and solute (J^α) fluxes relative to the solid phase, they are coupled and expressed in terms of the water chemical potential, μ^f, and the solute (electro)chemical potential, ${\widetilde {\rm{\mu }}^{\rm{\alpha }}}$:

(9)$$J^{\rm{f}} = - k\left( {{\rm{\rho }}_{\rm{T}}^{\rm{w}} \nabla {\rm{\mu }}^{\rm{f}} + \sum\limits_{\rm{\alpha }} {H_{\rm{\alpha }} c_{\rm{\alpha }} M_{\rm{\alpha }} \nabla \widetilde{\rm{\mu }}^{\rm{\alpha }} } } \right)$$

(10)$$J^{\rm{\alpha }} = H_{\rm{\alpha }} c_{\rm{\alpha }} J^{\rm{f}} - {{D_{\rm{\alpha }} {\rm{\rho }}_{\rm{\alpha }} } \over {RT}}\nabla \widetilde{\rm{\mu }}^{\rm{\alpha }} ,$$

where ${\rm{\rho }}_{\rm{T}}^{\rm{w}}$ is the true mass density of water; H _α , c _α , M _α , D _α, and ρ_α are the convection coefficient, the concentration (per unit water volume), the molar weight, the diffusivity, and the apparent mass density of solute α, respectively. Equation 9 is a more sophisticated version of Equation 4. As for the electrically charged solutes, Na⁺ (c ₊) and Cl^– (c _–) concentrations are related through the electroneutrality condition, c ₊= c _–+ c _fc , where c _fc is the fixed charge density of the tissue, as defined in Equation 8.

Simulating disc tissue degeneration and damage

The production of macromolecules changes as discs degenerate. The use of the Donnan theory in a poromechanical IVD model that simulated the behavior of cracks has suggested that tissue damage and proteoglycan depletion are tightly linked together.^{Reference Wognum, Huyghe and Baaijens36} In general, both triphasic theory and Donnan theory-based models can explicitly simulate proteoglycan depletion through the fixed charge density variable (Equations 8, 9, 10).^{Reference Zhu, Jackson and Gu47,Reference Ruiz, Malandrino, Lacroix, Ito and Noailly48} The collagen balance is also altered as degeneration progresses,^{Reference Coventry, Ghormley and Kernohan49–Reference Iatridis, Setton, Foster, Rawlins, Weidenbaum and Mow51} which was simulated by expressing ϕ₀ (Equation 5) as a function of the total collagen content, in order to simulate disc tissue dehydration as a function of fibrosis.^{Reference Schroeder, Wilson, Huyghe and Baaijens31}

Finite element (FE) disc models, based on the osmo-porohyperelastic constitutive framework presented in the previous sections of this article, reproduced the time-dependent behavior of healthy specimens or tissue samples.^{Reference Malandrino, Noailly and Lacroix17,Reference Jacobs, Cortes, Peloquin, Vresilovic and Elliott20,Reference Schroeder, Wilson, Huyghe and Baaijens31,Reference Malandrino52} However, the values of porohyperelastic parameters such as K, G (Equation 1), and k (Equations 4 and 9) depend on multiple interactions between proteoglycan, collagen, water, and tissue damage. Even using electrochemomechanical models, it proved difficult to correlate these values to pathophysiological tissue variations.^{Reference Schroeder, Huyghe, van Donkelaar and Ito32} A wide range of clinical scenarios of degenerative changes (e.g., calcification and water and disc height loss) was simulated through randomized nonlinear poroelastic FE models.^{Reference Galbusera, Schmidt, Neidlinger-Wilke, Gottschalk and Wilke53} Realistic trends were predicted, such as early increase of segmental stiffness with nucleus, annulus and end plate degeneration, disc height loss, osteophytosis, and diffuse sclerosis. Moreover, the results allowed the validity of grading systems to be inferred based on those degenerative changes, which significantly affect the spinal biomechanics.

A later study explored the ability of osmo-porohyperelastic models to simulate disc tissue degenerative changes, from degeneration grade 2 to degeneration grade 4 on the Pfirrmann scale, through specimen-specific FE models. Model parameter values were numerically optimized against the results of in vitro measurements on the specimens modeled.^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} Remarkably, physically sound simulations of tissue degeneration were achieved only with a model inspired by the micromechanics of ellipsoidal crack inclusions.^{Reference Dormieux and Kondo54} While representing the crack density increase through a damage parameter, d, the degenerated disc model simultaneously accounted for a loss of mechanical competence via damaged shear (G _d) and bulk moduli (K _d) and via increased initial porosity of the damaged continuum ϕ_0d:

(11)$${G_{\rm{d}}} = {G \over {1 + {Q_{\rm{G}}}d}}$$

(12)$${K_{\rm{d}}} = {K \over {1 + {Q_{\rm{k}}}d}}$$

(13)$${\rm{\phi }}_{{\rm{0d}}} = {\rm{\phi }}_{\rm{0}} + {{4Xd{\rm{\pi }}} \over 3},$$

where Q _G and Q _K are derived from the undamaged Poisson’s ratio of the solid matrix, and X is the aspect ratio of the elliptical cracks. The general osmo-porohyperelastic framework (Equations 1–5) was modified with Equations 11–13. Competing with damage, optimized model parameters revealed an increase in the undamaged moduli G and K, and a decrease in the undamaged porosity with degeneration, according to the expected tissue fibrosis. Simultaneously, decreased swelling pressure in the NP accompanying degeneration was correlated with reported loss of proteoglycans. Finally, a high value of d in the AF for Grade 4 discs correlated with the presence of cracks revealed by magnetic resonance imaging (MRI)^{Reference Antoniou, Epure, Michalek, Grant, Iatridis and Mwale55–Reference Munter, Wasserman, Wu and Yousem58} and IVD softening from mild to severe degeneration.^{Reference Gardner-Morse59,Reference Panjabi, Krag and Chung60}

Cell nutrition and viability models

In order to simulate indirect mechanotransduction effects, models have incorporated a description of anaerobic glycolysis, the main source of energy for disc cells:^{Reference Maroudas, Stockwell, Nachemson and Urban61–Reference Holm, Maroudas, Urban, Selstam and Nachemson63} glucose is used to produce adenosine triphosphate, and lactic acid is the resulting metabolite. Known as the Pasteur effect, extremely low oxygen levels result in a fall in oxygen consumption and a rise in the rate of glycolysis.^{Reference Bibby, Jones, Lee, Yu and Urban64} Hence, oxygen, lactate, and glucose are necessarily coupled.

These couplings were measured by Bibby and co-workers^{Reference Bibby, Jones, Ripley and Urban65} who assessed the metabolism of bovine NP cell cultures through electrochemical and biochemical sensors. They found that (1) the ratio between lactic acid production, R _l, and glucose consumption, R _g, was very close to two, as expected for a glycolytic metabolism; (2) oxygen consumption, R _O2, and lactate production were interrelated through oxygen concentration and pH; and (3) starting from 7.4, pH decreased linearly with the lactate concentration.^{Reference Bibby, Jones, Ripley and Urban65,Reference Soukane, Shirazi-Adl and Urban66} The corresponding empirical equations have been integrated into several metabolic transport models of the IVD^{Reference Malandrino, Noailly and Lacroix41,Reference Zhu, Jackson and Gu47,Reference Huang and Gu67–Reference Malandrino, Noailly and Lacroix73} through standard diffusion-reaction models of the three solutes:

(14)$${\partial \over {\partial t}}\left( {\matrix{ {C_{{\rm{O}}_{\rm{2}} } } \cr {C_1 } \cr {C_{\rm{g}} } \cr } } \right) - \left( {\matrix{ {D_{{\rm{O}}_{\rm{2}} } } & {\rm{0}} & 0 \cr {\rm{0}} & {D_{\rm{t}} } & 0 \cr {\rm{0}} & 0 & {D_{\rm{g}} } \cr } } \right)\nabla ^2 \left( {\matrix{ {C_{{\rm{O}}_{\rm{2}} } } \cr {C_1 } \cr {C_{\rm{g}} } \cr } } \right) = \left( {\matrix{ {R_{{\rm{O}}_{\rm{2}} } } \cr {R_1 } \cr {R_{\rm{g}} } \cr } } \right),$$

where C _i, D _i, and R _i are respectively the concentrations, the tissue-homogenized (effective) diffusion coefficients, and the cellular reactions, for oxygen (i = O₂), glucose (i = g), and lactate (i = 1). While this equation is for diffusive transport, some authors have also incorporated advection, but found that this did not influence the solute concentrations throughout the IVD,^{Reference Ferguson, Ito and Nolte39,Reference Malandrino, Noailly and Lacroix41} in agreement with experimental evidence.^{Reference Holm, Maroudas, Urban, Selstam and Nachemson63,Reference Urban, Holm, Maroudas and Nachemson74,Reference Katz, Hargens and Garfin75} This outcome was correlated with the low hydraulic permeability of disc tissues leading to insignificant traveling distances of fluid-dragged particles.^{Reference Malandrino, Noailly and Lacroix41} As for couplings with mechanics, the volume changes updated through the osmo-porohyperelastic calculations affected both the diffusion distances and the diffusivity coefficients that depended on porosity (Equation 5).^{Reference Malandrino, Noailly and Lacroix41,Reference Zhu, Jackson and Gu47,Reference Huang and Gu67,Reference Jackson, Huang, Brown and Gu71,Reference Malandrino, Noailly and Lacroix73}

In Equation 14, the effective reactions, R _i, depend on a cell density parameter, ρ_cell, the value of which affects the mechano-transport calculations.^{Reference Malandrino, Noailly and Lacroix41} Diffusion chamber experiments have shown that IVD cells can survive at an extremely low oxygen level, but viability is compromised when low glucose levels and/or acidic pH occur.^{Reference Razaq, Wilkins and Urban76–Reference Horner and Urban78} More complex behaviors involving cell clusters were also observed but cell proliferation rates were slow.^{Reference Johnson, Eisenstein and Roberts79} Thus, models could reasonably consider ρ_cell as decreasing or constant, depending on environmental conditions, and cell viability was defined as the ratio between current and an initial cell density (i.e., ρ_cell/ρ_cell,0).

Based on experimental studies,^{Reference Bibby, Fairbank, Urban and Urban77,Reference Horner and Urban78,Reference Bibby and Urban80,Reference Guehring, Wilde, Sumner, Grünhagen, Karney, Tirlapur and Urban81} Shirazi-Adl et al.^{Reference Shirazi-Adl, Taheri and Urban82} and Jackson et al.^{Reference Jackson, Huang, Brown and Gu71} first incorporated cell viability into numerical analyses by assuming that ρ_cell varies instantaneously with glucose concentrations in the IVD. The lack of time-dependence of cell death was an important limitation, along with a possible cell resurrection if concentrations increased above the threshold. More recently, Zhu et al.^{Reference Zhu, Jackson and Gu47} incorporated a new constitutive model for disc cell viability into the triphasic theory (Equations 9 and 10). The rate of cell density change was based on the theoretical framework for microbial cell growth:^{Reference Lin, Lee, Lee and Koo83}

(15)$${{\rm{\rho }}_{{\rm{cell}}}} = {{\rm{\rho }}_{{\rm{cell,0}}}}{{\rm{e}}^{{\rm{ - \zeta }}t}},$$

where ζ is the rate of normalized cell density change, which depends on the solute concentrations as follows:

(16)$${\rm{\zeta }} = {\rm{\alpha }}\left( {{{C_{\rm{g}} - C_{{\rm{g,0}}} } \over {C_{\rm{g}} + k_1 }} - {{\left| {C_{\rm{g}} - C_{{\rm{g,0}}} } \right|} \over {C_{\rm{g}} + k_2 }}} \right) + {\rm{\alpha }}_{\rm{p}} {\rm{H}}\,{\rm{,}}$$

where C _g,0 is the minimum glucose concentration necessary for IVD cell survival. The positive parameters α, k ₁ and k ₂ were determined by the biology of IVD cells. This cell density model eliminated the problem of resurrection,^{Reference Lin, Lee, Lee and Koo83} but the effect of pH on cell viability was disregarded (i.e., α_pH = 0 in Equation 19). Hence, Malandrino and co-workers considered a constant α_pH that differed from zero when critical pH conditions were met.^{Reference Malandrino, Noailly and Lacroix73} When calibrated from experimental data, the critical levels, in computational models, were obtained to be 0.5 mM for glucose (C _g,0 in Equation 19) and 6.8 for pH,^{Reference Zhu, Jackson and Gu47,Reference Malandrino, Noailly and Lacroix73,Reference Shirazi-Adl, Taheri and Urban82} in agreement with independent experimental evidence.^{Reference Razaq, Wilkins and Urban76,Reference Bibby, Fairbank, Urban and Urban77}

Moreover, the nutrition-related diffusion-reaction and cell viability models (Equations 14–16) have been explicitly assessed: diffusion chamber experiments were compared to simulations in terms of viability profiles.^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33,Reference Shirazi-Adl, Taheri and Urban82} The decreased viability toward the center of the diffusion chamber was predicted satisfactorily for cell densities characteristic of the human IVD (i.e., 4–8 million cells ml^–1 at steady state^{Reference Malandrino, Noailly and Lacroix73,Reference Shirazi-Adl, Taheri and Urban82} and early time points, such as 0–3 days).^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} For these times and initial cell density ranges, the role of pH on cell death appeared negligible compared to glucose deprivation in IVD models.^{Reference Malandrino, Noailly and Lacroix73,Reference Shirazi-Adl, Taheri and Urban82} This outcome might be the consequence of a relatively low diffusivity of glucose compared to the other solutes considered in Equation 14. Yet, future studies are needed to understand the role of pH in amplifying death cascades expected for higher initial densities.

Model validation

In this review, we have highlighted couplings among disc condition, morphology, time-dependent mechanics, and cell nutrition. These interactions are local and difficult to apprehend through in vitro or in vivo studies only. Computational approaches clearly favor educated interpretations of complex phenomena, but the fundamental question of model validation is a difficult one. The use of osmo-porohyperelastic approaches has been shown to be valid for healthy IVD models, at the tissue^{Reference Schroeder, Huyghe, van Donkelaar and Ito32} and organ levels,^{Reference Simon, Wu, Carlton, Evans and Kazarian12,Reference Malandrino, Noailly and Lacroix17,Reference Jacobs, Cortes, Peloquin, Vresilovic and Elliott20,Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33,Reference Ferguson, Ito and Nolte39,Reference Peloquin, Yoder, Jacobs, Moon, Wright, Vresilovic and Elliott84–Reference Williams, Natarajan and Andersson86} and for mildly to moderately/severely degenerated disc models at the organ level.^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} The time-dependent and nonlinear mechanical response of the disc, including intradiscal pressure measurements could be reproduced. Importantly, the match between model predictions and experimental measurements can be increased (i.e., uncertainty within the 95% confidence interval calculated from the experimental results) through the use of patient-specific model geometries.^{Reference Jacobs, Cortes, Peloquin, Vresilovic and Elliott20,Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33}

The possibility to validate models against solute concentration measurements remains an open issue (see “limitations” section). However, IVD poromechanical and cell nutrition models can be used to simulate organ culture experiments that provide alternative source of data for validation.^{Reference Castro, Paul, Detiger, Smit, van Royen, Pimenta Claro, Mullender and Alves87,Reference Malandrino, Lacroix and Noailly88} For instance, in a pilot study, cell viability models predictions have been shown to match local experimental cell death data.^{Reference Malandrino, Lacroix and Noailly88} This validation has served to propose severe poroelastic tissue compaction as a relevant phenomenon in indirect disc mechanotransduction.

Major simulation results

Effect of poromechanical properties and simulated disc degeneration

Integrative studies can provide information on the respective effects of solid- and fluid-phase properties on the reactive transport of metabolites. Under compressive loads that simulated daily activity,^{Reference Malandrino, Noailly and Lacroix41} the load-dependency of solute transport throughout the disc was predicted to arise mainly from (1) the variation of the effective solute diffusivity in the continuum, and (2) the change of diffusion distances from and to the simulated vasculature at the disc periphery (Figure 2). A competitive interplay of (1) and (2) was identified: while tissue compaction limited diffusive transport in the pores of the disc matrix, distance reduction between cranial and caudal boundary vascular supplies enhanced diffusion. The latter beneficial effect was predominant in compliant models, in contrast to stiffer models.^{Reference Malandrino, Noailly and Lacroix41} Variations in porosity, initial cell density, and osmotic pressure in the NP mostly influenced the calculated solute distributions.

Figure 2.

Illustration of the two main interactions resulting from coupling disc poromechanics and diffusive transport: strain-dependent diffusivity (bottom, left) and diffusion distances (bottom, right).^{Reference Malandrino, Noailly and Lacroix41} The cartoon depicts an example of the competition between the two effects. Compressive loading causes matrix compaction in the disc center, which results in a reduced diffusivity of the representative volume element (RVE) for a given solute (bottom, left). Also, deformation causes a shortening of the distance between the end plates (i.e., the vasculature nutrient supply). This would cause enhanced diffusion in the RVE (bottom, right).

In general, mild to moderate disc degeneration causes an overall stiffening of the disc due to fibrosis and water loss.^{Reference Coventry, Ghormley and Kernohan49–Reference Iatridis, Setton, Foster, Rawlins, Weidenbaum and Mow51} Jackson et al.^{Reference Jackson, Huang, Brown and Gu71} and Zhu et al.^{Reference Zhu, Jackson and Gu47} investigated this stiffening and other changes (reduced disc height, water content, fixed charge density, and end plate permeability) on disc cell viability predictions. These simulated changes led to an overall decrease of glucose concentration in the disc, which in turn resulted in local cell death. In general, the NP was more strongly affected than the AF, with cell death being initiated at the AF–NP interface and spreading toward the center of the disc.

Decreased compliance impeded glucose concentrations from being enhanced by mechanical loads, despite simulated height reduction, in degenerated disc models.^{Reference Malandrino, Noailly and Lacroix41,Reference Malandrino, Noailly and Lacroix73} Hence, reduced strain-dependent diffusivities could easily exacerbate local cell death, as a result of simulated tissue dehydration and compaction.^{Reference Malandrino, Noailly and Lacroix73} However, simulating about 40% to 50% of decreased nutrient supply at the end plate was necessary to trigger cell death,^{Reference Zhu, Jackson and Gu47,Reference Malandrino, Noailly and Lacroix73,Reference Shirazi-Adl, Taheri and Urban82,Reference Jackson, Huang and Gu89} remaining a significant contribution to the overall effect of NP and AF condition^.

The end plate is the major route of diffusive exchanges between the IVD and the vascular system (Figure 1b), and plays a key role on disc nutrition.^{Reference Shirazi-Adl, Taheri and Urban82,Reference Jackson, Huang and Gu89–Reference Lee, Iatridis, Poveda and Alini91} It comprises the CEP (not vascularized) and a subchondral bone layer or bony end plate (BEP, vascularized). Both CEPs and BEPs have been integrated in osmo-porohyperelastic disc models coupled to diffusion-reaction and cell viability calculations.^{Reference Malandrino, Lacroix, Hellmich, Ito, Ferguson and Noailly92,Reference Noailly, Malandrino and Galbusera93} Under mechanical loads, microstructure-based variations of BEP porosity and hydraulic permeability had negligible effects on the distribution of metabolites within the disc, and limited effects on the fluid velocities across the juxtaposed CEP.^{Reference Malandrino, Lacroix, Hellmich, Ito, Ferguson and Noailly92} Also, variations in the CEP hydraulic permeability affected neither diffusion nor advection of metabolites.^{Reference Malandrino, Noailly and Lacroix41}

However, if cartilaginous calcification is simulated as a decrease in CEP porosity and/or diffusivity, the negative effect on cell viability became striking.^{Reference Shirazi-Adl, Taheri and Urban82,Reference Jackson, Huang and Gu89,Reference Noailly, Malandrino and Galbusera93} Models reasonably suggested that the net supply of glucose and oxygen at the BEP–CEP interface drives cell survival. Unfortunately, changes or damage to the end plate with degeneration remains controversial,^{Reference Rodriguez, Slichter, Acosta, Rodriguez-Soto, Burghardt, Majumdar and Lotz94} and simulations highlight the imperative need to better understand these changes as a complex result of microvasculature condition and matrix deposition.^{Reference Lories and Luyten95}

Effect of disc geometry

In a recent work that included mechano-transport and cell viability simulations, FE geometries of lumbar discs were derived from MRI of patients with different degeneration grades.^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} Grade-dependent material properties were additionally achieved for both the AF and the NP. While the study confirmed earlier simulation findings (see previous discussion about predicted cell death occurrence) for discs that were about 10 mm high, IVD models higher than 14 mm became prone to harsh nutrient conditions and cell death within the NP. For such large discs, the effect of diffusion distances became predominant. This striking impact of disc size was in line with both experimental^{Reference Rodriguez, Slichter, Acosta, Rodriguez-Soto, Burghardt, Majumdar and Lotz94} and other simulation findings.^{Reference Motaghinasab, Shirazi-Adl, Parnianpour and Urban101}Figure 3 summarizes the possible causes of IVD cell death detailed above according to the described mechano-transport phenomena.

Figure 3.

Hierarchical schematization of the main causes leading to nutrient deprivation and cell death in the intervertebral disc, mediated through the main poromechanics-transport phenomena discussed in this article: possible degenerative onsets include tissue degenerative processes,^{Reference Coventry, Ghormley and Kernohan49–Reference Iatridis, Setton, Foster, Rawlins, Weidenbaum and Mow51} adverse mechanical loads,^{Reference Zhu, Jackson and Gu47} disc geometrical changes,^{Reference Malandrino, Noailly and Lacroix41,Reference Malandrino, Noailly and Lacroix73} and end plate issues^{Reference Shirazi-Adl, Taheri and Urban82,Reference Jackson, Huang and Gu89–Reference Noailly, Malandrino and Galbusera93} (blue region). These causes interact with a landscape of several disc properties (magenta region), which in turn activates poromechanics-dependent diffusive-reactive transport phenomena (green region). The competition between these two basic poromechanics-transport interactions (see also Figure 2) determines the net amount of metabolites across any given disc representative volume element and over time. Among metabolites, the glucose concentration has been identified as the major cause responsible for cell death in the discs.^{Reference Razaq, Wilkins and Urban76–Reference Horner and Urban78}

Limitations

For the verification of IVD poromechanical models, studies have shown that care should be taken that fluid pressures and velocities are calculated under load velocities similar to those in vivo.^{Reference Stokes, Chegini, Ferguson, Gardner-Morse, Iatridis and Laible102,Reference Ruiz, Noailly and Lacroix103} Also, mathematical artifacts must be avoided, as a model with a given set of parameters can match the experiment even when this set is devoid of physical meaning.^{Reference Malandrino, Planell and Lacroix14,Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} As for metabolic transport and cell viability models, comparisons to in vivo measurements of metabolite concentrations^{Reference Bartels, Fairbank, Winlove and Urban104} have been performed.^{Reference Malandrino, Noailly and Lacroix41,Reference Magnier, Boiron, Wendling-Mansuy, Chabrand and Deplano70} However, because these measurements are challenging and produce noisy data,^{Reference Bartels, Fairbank, Winlove and Urban104} model validation remains largely unresolved. Furthermore, experiments that determined the reaction terms in Equation 14 and those that were used to validate the cell viability model (Equations 15 and 16) were performed on NP cells only,^{Reference Bibby, Jones, Ripley and Urban65,Reference Horner and Urban78} and the energy metabolisms of AF and NP cells might differ.^{Reference Salvatierra, Yuan, Fernando, Castillo, Gu, Cheung and Huant105,Reference Fernando, Czamanski, Yuan, Gu, Salahadin and Huang106} Including this difference in modeling could be particularly relevant, since several studies have pointed out reduced glucose availability at the inner AF–NP interface,^{Reference Malandrino, Noailly and Lacroix41,Reference Zhu, Jackson and Gu47,Reference Malandrino, Noailly and Lacroix73} due to local tissue consolidation under external compressive loads. Also, whether in vitro measurements can be extrapolated to in vivo situations is rarely addressed, though critical to model validation.

Despite these limitations, IVD poromechanics-transport simulations have allowed for understanding complex interactions, as discussed in the previous sections. Most interactions arise from fundamental physical laws, independent of the particular tissue/organ modeled (e.g., tissue consolidation involves water loss and diffusivity decrease, and tissue deformation involves changes in concentration gradients and diffusion distances). However, the relative relevance of each phenomenon in the disc models is tied to the parameter values. Accordingly, a detailed picture of the sensitivity of the model to parameter changes is mandatory, and might even be the core of complex model exploitation, so as to explore mechanisms that cannot be measured directly.

In this article, robust arguments of discussion and alternative validation paths have been derived. Importantly, the ability of the diffusion-reaction and cell viability models to predict reliable cell distributions over time and across space^{Reference Malandrino, Pozo, Castro-Mateos, Frangi, van Rijsbergen, Ito and Noailly33} adds to the credibility of disc nutrition calculations. Also, some simulation-based findings such as the possible relevance of both porosity and disc size in predisposition to degeneration are clinically assessable through cohort studies and image-based measurements.