2.2 Constitutive models
The vertebral bodies (VBs), CEP, AF, and NP all utilize structure-based solid mixture constitutive laws. These laws effectively link the behavior with its constituents. VBs' constitutive behavior was modeled using a compressible neo-Hookean framework. This approach allows the hyperelastic strain-energy density to be expressed as follows
$$\begin{array}{c}W=\frac{\mu }{2}\left({I}_{1}-3\right)-\mu \text{ln}J+\frac{\lambda }{2}{\left(\text{ln}J\right)}^{2},\#\left(1\right)\end{array}$$
where J represents the determinant of the deformation gradient tensor (F), \({I}_{1}\)signifies the first invariant of the right Cauchy–Green deformation tensor (\(C={F}^{T}F\)), λ and µ are Lamé constants that relate to Young’s modulus (E) and Poisson’s ratio (ν).
$$\begin{array}{c}\lambda =\frac{Ev}{\left(1+v\right)\left(1-2v\right)}, \mu =\frac{E}{2\left(1+v\right)}.\#\left(2\right)\end{array}$$
In this work, the VBs were modeled as a biphasic neo-Hookean material with the Young’s modulus E = 12 GPa and Poisson’s ratio ν = 0.321. In addition, a constant permeability of 𝑘0 = 5 mm4/Ns was chosen22, in order to allow fluids to flow through.
The SCEP, ICEP, AF, NP, each comprise an isotropic, compressible, and hyperelastic solid matrix, characterized by the Holmes-Mow model. The strain-energy density for the Holmes-Mow matrix is described as follows:
$${ {\Psi }}_{m}\left({I}_{1}, {I}_{2}, J\right)=\frac{\lambda +2\mu }{4{\beta }_{m}}\left({e}^{\frac{{\beta }_{m}}{\lambda +2\mu }\left[-\left(\lambda -2\mu \right)\left({I}_{1}-3\right)+\lambda \left({I}_{2}-3\right)-\left(\lambda +2\mu \right)\text{ln}\left({J}^{2}\right) \right]}-1\right) \left(3\right)$$
Denoting the second invariant of the right Cauchy-Green deformation tensor is \({I}_{2}\), Characterizing the rate of the exponential stiffening is \({\beta }_{m}\). These parameters, along with Young’s modulus (\({E}_{m}\)), Poisson’s ratio (\({\nu }_{m}\)) and other material constants,are specified in Table 2 23,24.
Table 2
Material properties for the constitutive models were determined based on experimental data from existing literature19,23,24,26,27,47,48.( SCEP and ICEP stand for superior and inferior cartilaginous endplates, respectively)
Property | Symbol | SCEP | ICEP | AF | NP |
Rate of matrix stiffening | 𝛽𝑚 (unitless) | 0.29 | 0.29 | 1.5(3.4)a | 0.95 |
Initial water fraction | \({\phi }_{0}^{w}\) (unitless) | 0.6 | 0.6 | 0.7 | 0.8 |
Nonlinearity of permeability | M (unitless) | 5.2 | 3.3 | 2.7(5.7)a | 1.92 |
Initial permeability | \({k}_{0}\) (10− 4 mm4/Ns) | 8.24 | 2.87 | 16(47)a | 5.5 |
Poisson’s ratio | 𝜈𝑚 (unitless) | 0.18 | 0.18 | 0.16(0.24)a | 0.24 |
Osmotic coefficient | \({\Phi }\) (unitless) | 1 | 1 | 1 | 1 |
Modulus of fibers | \(\xi\)1/𝐸f (MPa) | 7.01b | 7.01b | 3.0(15.6)a | N.A. |
Rate of fiber stiffening | \({\beta }_{c}/{\beta }_{f} \left(\text{u}\text{n}\text{i}\text{t}\text{l}\text{e}\text{s}\text{s}\right)\) | 2.88c | 2.88c | 4 | N.A. |
Initial fixed charge density | \({c}_{f0}\) (mmol/L) | −250 | −450 | −300(− 100)a | −400 |
Critical stretch square | \(I{}_{0}\) (unitless) | N.A. | N.A. | 1.14(1.06)a | N.A. |
Modulus of matrix | 𝐸𝑚 (MPa) | 0.305 | 0.305 | 0.045(0.018)a | 0.065 |
a the parameter exhibits a linear variation from the innermost to the outermost lamella of the AF, as indicated by the value in parentheses; |
b the initial elastic modulus of fibers within the SCEP and ICEP (\(\xi\)1); |
c Nonlinearity in fiber strain-energy density within the SCEP and ICEP. (\({\beta }_{c}\))。 |
The solid matrix of the IVD also incorporates an isotropic and deformation-dependent permeability behavior, as described by the Holmes-Mow permeability model. 25.
$$\begin{array}{c}k={k}_{0}{\left(\frac{J+{\phi }_{0}^{w}-1}{{\phi }_{0}^{w}}\right)}^{2}{e}^{\frac{M}{2}\left({J}^{2}-1\right)}I,\#\left(4\right)\end{array}$$
where\(\varvec{k}\) is the permeability tensor, \({k}_{0}\) is the initial isotropic hydraulic permeability in the reference configuration, \({\phi }_{0}^{w}\) is the water volume fraction, and \(M\) is the exponential coefficient of permeability. Both \(M\) and \({k}_{0}\) were determined by fitting to confined compressive experimental data 23,24,26. (Table 2和Table 3).
Within the framework of biphasic-swelling theory, glycosaminoglycan (GAG) content correlates with fixed charge density. The fixed charge density, prior to deformation, is established in the reference configuration as follows:
$$\begin{array}{c}{c}_{f0}=\frac{{c}_{g}{Z}_{c}}{{M}_{c}},\#\left(5\right)\end{array}$$
where \({c}_{g}\) is the GAG content in water with the dimension mg/ml, and \({M}_{c}\) and \({Z}_{c}\) are the molecular weight and number of charges of chondroitin-6-sulfate disaccharide. The latter two parameters were chosen as \({M}_{c}=502.5\) gram and \({Z}_{c}=2\) charges per repeating unit. Moreover, to reflect the loss of GAG content in degenerated IVD models, lower levels of initial fixed charge densities were assigned for the AF and NP (Table 3).
Table 3
Initial material property values were determined based on experimental data cited in the literature.47–50.
Property | Symbol | Group1 | Group2 | Group3 | Group4 |
Nonlinearity of SCEP | \(M\) (unitless) | 5.99 | 4.49 | 5.99 | 4.49 |
Nonlinearity of ICEP | \(M\) (unitless) | 3.79 | 3.08 | 3.79 | 3.08 |
Permeability of SCEP | \({k}_{0}\) (10− 4 mm4/Ns) | 8.24 | 3.06 | 8.24 | 3.06 |
Permeability of ICEP | \({k}_{0}\) (10− 4 mm4/Ns) | 2.87 | 1.62 | 2.87 | 1.62 |
Initial fixed density of SCEP | \({c}_{f0}\) (mmol/L) | 110 | 40 | 110 | 40 |
Initial fixed density of ICEP | \({c}_{f0}\) (mmol/L) | 110 | 40 | 110 | 40 |
Fixed charge density of AF | \({c}_{f0}\) (mmol/L) | -300 (-100)a | -300 (-100)a | -150 (-100)a | -150 (-100)a |
Fixed charge density of NP | \({c}_{f0}\) (mmol/L) | -350 | -350 | -125 | -125 |
a The parameter varies linearly across the lamellae of the AF, with the inner lamella value in parentheses. |
Upon applying external loads, the instantaneous fixed charge density changes according to the deformation Jacobian(\(J\)), initial water volume fraction (\({\phi }_{s}^{w}\)) and initial fixed charge density (\({c}_{f0}\))
$$\begin{array}{c}{c}_{f}=\frac{{c}_{f0}{\phi }_{0}^{w}}{J-1+{\phi }_{0}^{w}}.\#\left(6\right)\end{array}$$
As deformation occurs, changes in the fixed charge density lead to variations in osmotic pressure.
$$\begin{array}{c}{p}_{o}=RT\varPhi \left(\sqrt{{c}_{f}^{2}+{c}_{b}^{2}}-{c}_{b}\right),\#\left(7\right)\end{array}$$
The osmotic pressure equation incorporates the gas constant(\(R\)), temperature(\(T\)),osmotic coefficient(Φ),and the osmolarity of the external bath(\({c}_{b}\)).The entire intervertebral disc (IVD) model was submerged in a 0.15 M phosphate-buffered saline solution, reflecting ideal bath conditions as outlined in Table 2 for this study.
It is noted that, in addition to the fluid pressure (\({p}_{f}\)) and solid stress (\({\sigma }_{s}\)), the osmotic pressure(\({p}_{o}\)) also contributes to the total Cauchy stress (\(\varvec{\sigma }\))
$$\begin{array}{c}\sigma =-\left({p}_{f}+{p}_{o}\right)I+{\varvec{\sigma }}_{s},\#\left(8\right)\end{array}$$
where \(\varvec{I}\) is the identity tensor.
Besides the solid matrix, collagen fibers are embedded within the SCEP, ICEP, and AF, enabling them to withstand both compressive and tensile stresses. For the SCEP and ICEP, the behavior of these fibers is captured using an ellipsoidal fiber distribution model, whose strain-energy density function is defined as follows:
$$\begin{array}{c}{{\Psi }}_{f}\left(\varvec{n},{I}_{n}\right)=\xi \left(\varvec{n}\right){\left(I{}_{n}-1\right)}^{\beta \left(\varvec{n}\right)},\#\left(9\right)\end{array}$$
where\({I}_{n}\) represents the instantaneous square of the fiber stretch, i.e.,\({I}_{n}= \varvec{N}\cdot C\cdot \varvec{N}\), with \(\varvec{N}\) being the unit vector along the fiber axis in the reference configuration. \(\varvec{n}\) is the unit vector of fiber in the current reference, i.e., \(\varvec{n}=\varvec{F}\cdot \varvec{N}/{\lambda }_{n}\). In spherical coordinate system, the unit vector (\(\varvec{n}\)) is measured by the longitudinal (\(\phi\)) and latitudinal (\(\theta\)) coordinates. With their help, the material parameters\(\beta \left(\varvec{n}\right)\)and \(\lambda \left(\varvec{n}\right)\) can be defined as
$$\begin{array}{c}\xi \left(\varvec{n}\right)={\left(\frac{{\text{cos}}^{2}\theta {\text{sin}}^{2}\phi }{{\xi }_{1}^{2}}+\frac{{\text{sin}}^{2}\theta {\text{sin}}^{2}\phi }{{\xi }_{2}^{2}}+\frac{{\text{cos}}^{2}\phi }{{\xi }_{3}^{2}}\right)}^{-\frac{1}{2}},\#\left(10\right)\end{array}$$
$$\begin{array}{c}\xi \left(\varvec{n}\right)={\left(\frac{{\text{cos}}^{2}\theta {\text{sin}}^{2}\phi }{{\beta }_{1}^{2}}+\frac{{\text{sin}}^{2}\theta {\text{sin}}^{2}\phi }{{\beta }_{2}^{2}}+\frac{{\text{cos}}^{2}\phi }{{\beta }_{3}^{2}}\right)}^{-\frac{1}{2}},\#\left(11\right)\end{array}$$
where \({\xi }_{i}\)and \({\beta }_{i}\)(\(i=\text{1,2},3\)) stand for the initial moduli and fiber nonlinearity (power coefficient) along three principal axes of the ellipsoidal fibers. In the CEP, fibers are appear primarily oriented in the planes that are parallel to the VBs and the NP. Along the lateral and anterior-posterior directions, the fibers share the same modulus (\({\xi }_{1}={\xi }_{2}\)). Along the axial direction of the IVD, much smaller elastic modulus was assumed, e.g., \({\xi }_{3}=0.1{\xi }_{1}\). The fiber nonlinearity was assumed to be independent of coordinate directions, e.g., \(\beta \left(\varvec{n}\right)={\beta }_{c}\), with\({\beta }_{c}\) being a constant. Following a regression analysis by fitting to tensile experimental data, the modulus \({\xi }_{1}={\xi }_{2}=10 {\xi }_{3}=7.01\)MPa and the fiber nonlinearity \({\beta }_{c}=2.88\) were obtained26(Table 2).
In the AF, the collagen fibers embedded within are represented by a toe-linear stress-stretch constitutive law, aiming to capture their nonlinear mechanical properties. In terms of the square of the fiber stretch, the strain-energy density can be ued to expresse (\({I}_{n}\))
$$\begin{array}{c}{{\Psi }}_{f}=\left\{\begin{array}{c}0, {I}_{n}<1, \\ \frac{\xi }{{\beta }_{f}}{\left({I}_{n}-1\right)}^{{\beta }_{f}}, 1\le {I}_{n}\le {I}_{0}, \\ B\left({I}_{n}-{I}_{0}\right)-{E}_{f}\left(\sqrt{{I}_{n}}-\sqrt{{I}_{0}}\right)+\frac{\xi }{{\beta }_{f}}{\left({I}_{0}-1\right)}^{{\beta }_{f}}, {I}_{0}\le {I}_{n},\end{array}\right.\#\left(12\right)\end{array}$$
where \({I}_{0}\) represent for the critical square of the fiber stretch (\({\lambda }_{0}^{2}\)) corresponding to the point that separates the linear and the toe regions, and\(\xi\)and \(B\) are both functions of the linear-region elastic modulus (\({E}_{f}\)), \({I}_{0}\) and the toe-region power-law coefficient (\({\beta }_{f}\))
$$\begin{array}{c}\xi =\frac{{E}_{f}}{4\left({\beta }_{f}-1\right)}{I}_{0}^{-\frac{3}{2}}{\left({I}_{0}-1\right)}^{2-{\beta }_{f}}, B=\xi {\left({I}_{0}-1\right)}^{{\beta }_{f}-1}+\frac{{E}_{f}}{2}{I}_{n}^{-\frac{1}{2}}.\#\left(13\right)\end{array}$$
Parameters of the constitutive law such as \({E}_{f}\), \({I}_{0}\)and βf were calibrated in terms of the experimental stress-stretch curve of single-lamella AF samples 27. From the inner to the outer lamellae, \({E}_{f}\) increases linearly, whereas \({I}_{0}\) and \({\beta }_{f}\)decreases linearly (Table 2). In addition, the angle between the fiber axis and the transverse cross-section decreases from ± 40° in the innermost lamella to ± 28° in the outermost one.
The axial element size of the innermost lamella is determined according to the collagen fiber orientation;Assuming that the axial element size of the innermost lamella is 𝑙𝑎, the average collagen fiber direction of the innermost lamella is 𝜃 = 40∘ with the transvers cross-section plane, and the element size of the inner boundary curve, is 𝑙𝑖𝑛. These three parameters approximately satisfy
$$\begin{array}{c}tan\theta =\frac{{l}_{\text{a}}}{{l}_{in}}\#\left(14\right)\end{array}$$