
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13555-4
10.1016/j.heliyon.2024.e37524
e37524
Research Article
Finite element model reveals the involvement of cartilage endplate in quasi-static biomechanics of intervertebral disc degeneration
Zhang Yujun a1
Pan Yanli a1
Mao Xinning a1
He Du a
Zhang Liangping a
Cheng Wei ab
Zhu Chengyue abc
Zhu Hang ab
Zhang Wei ab
Jin HongTing hongtingjin@163.com
ad⁎⁎
Pan Hao harper1966@alu.zcmu.edu.cn
ac⁎
Wang Dong hz-wdong@alu.zcmu.edu.cn
abc⁎⁎⁎
a Department of Orthopaedics, Hangzhou TCM Hospital Affiliated to Zhejiang Chinese Medical University (Hangzhou Hospital of Traditional Chinese Medicine), Hangzhou 310000, Zhejiang Province, China
b Department of Orthopaedics, Hangzhou Dingqiao Hospital, Huanding Road NO 1630, Hangzhou 310021, Zhejiang Province, China
c Institute of Orthopaedics and Traumatology, Hangzhou Traditional Chinese Medicine Hospital Affiliated to Zhejiang Chinese Medical University, Tiyuchang Road NO 453, Hangzhou 310007, Zhejiang Province, China
d Institute of Orthopaedics and Traumatology, The First Affiliated Hospital of Zhejiang Chinese Medical University, Hangzhou, China
⁎ Corresponding author. Department of Orthopaedics, Hangzhou TCM Hospital Affiliated to Zhejiang Chinese Medical University Hangzhou 310000, Zhejiang Province, China. harper1966@alu.zcmu.edu.cn
⁎⁎ Corresponding author. Institute of Orthopaedics and Traumatology, The First Affiliated Hospital of Zhejiang Chinese Medical University, Hangzhou, China. hongtingjin@163.com
⁎⁎⁎ Corresponding author. Department of Orthopaedics, Hangzhou TCM Hospital Affiliated to Zhejiang Chinese Medical University Hangzhou 310000, Zhejiang Province, China. hz-wdong@alu.zcmu.edu.cn
1 These authors contributed equally to this work.

05 9 2024
30 9 2024
05 9 2024
10 18 e3752425 4 2024
27 8 2024
4 9 2024
© 2024 The Authors. Published by Elsevier Ltd.
2024

https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Background and objective

The intrinsic link between the compositional and structural attributes and the biomechanical functionality is evident in intervertebral discs. However, it remains unclear from a biomechanical perspective whether cartilage endplate (CEP) degeneration exacerbates intervertebral disc degeneration.

Methods

This study developed and quantitatively validated four biphasic swelling-based finite element models. We then applied four quasi-static tests and simulated daily loading scenarios to examine the effects of CEP degradation.

Results

Under free-swelling conditions, short-term responses were prevalent, with CEP performance changes not significantly impacting response proportionality. The creep test results showed the more than 50 % of the strain was attributed to long-term responses. Stress-relaxation testing indicated that all responses increased with disc degeneration, yet CEP degeneration's impact was minimal. Daily load analyses revealed that disc degeneration significantly reduces nucleus pulposus pressure and disc height, whereas CEP degeneration marginally increases nucleus pressure and slightly decreases disc height.

Conclusions

Glycosaminoglycan content and CEP permeability are critical to the fluid-dependent viscoelastic response of intervertebral discs. Our findings suggest that CEP contributes to disc degeneration under daily loading conditions.

Keywords

Finite element
Cartilage endplate
Degeneration of the intervertebral disc
Quasi-static biomechanics
Daily loading
==== Body
pmc1 Introduction

As one of the most common musculoskeletal conditions, low back pain (LBP) is primarily caused by intervertebral disc degeneration (IVDD). Induced by IVDD, LBP has a global incidence rate of 9.4 %. Over the past decade, related healthcare costs have increased by 43 %, significantly impairing individuals' quality of life and imposing substantial economic burdens on society [1].Intervertebral discs (IVDs) principally utilizes endplate diffusion for the transport of nutrients and metabolites, thereby maintaining a stable internal environment [2]. Endplate degeneration leads to reduced permeability, obstructing this essential exchange and aggravating IVDD [3]. This phenomenon underscores the interplay between the disc's compositional changes and its biomechanical functions, highlighting the importance of exploring the biomechanical properties of cartilage endplates under varying conditions.

The nucleus pulposus (NP) is centrally located within the intervertebral discs (IVDs), enveloped by the cartilaginous endplate (CEP) and annulus fibrosus (AF), crucial for maintaining the disc's biomechanical function [4].NPs and AFs are mainly composed of water, glycosaminoglycan (GAG), and collagen fibers. Collagen fibers form a mesh in the IVD, providing the tissue with high tensile and stiffness strength. GAG is integrated into the collagen fiber mesh, maintaining water content of the IVD, providing osmotic and swelling pressure, and maintaining the height and swelling of the IVD under heavy pressure or impact [5]. The proportions of composition between the two tissues were different. NPs are similar to gelatinous tissue with a higher GAG ratio. The AF's collagen fibers are meticulously organized in cross-layered concentric sheets, approximately 20 layers thick [6]. Under typical conditions, the AF and NP collaborate to ensure the IVD possesses the mechanical properties needed to withstand stress on the vertebral body [7]. When degeneration occurs, the significant loss of GAG in these structures reduces the IVD's hydration capacity. In the meantime, the decrease in type II collagen fibers further reduces the mechanical integrity and elasticity, leading to a decreased structural support ability [8]. The CEP, a transparent cartilage layer about 0.6 mm thick, connects the IVD to the neighboring vertebral bodies [9]. It serves as a transition layer that alleviates stress at the disc-bone interface, effectively balancing the flexibility of the disc with the rigidity of the vertebral bone. Its low permeability maintains pressure within the disc's interstitial fluid, forming a barrier that limits fluid movement and enables the disc to withstand compressive forces [10]. Additionally, the CEP prevents disc bulging into adjacent vertebral bodies and provides structural anchorage for the fibers of the inner AF and NP [11]. CEP contains about 60 % water, with the remaining components primarily consisting of type II collagen and proteoglycans [12]. The tensile modulus of the CEP is primarily influenced by collagen content, while GAG content significantly influences osmotic properties. Multiphysics models indicate that GAG plays a crucial role in maintaining disc stability by effectively regulating fluid flux between the CEP and NP through the CEP [13]. CEP degeneration is thought to play a role in the early stages of IVDD.When the CEP undergoes degeneration, such as calcification [14] and dehydration [15], a decrease in GAG content can lead to a decrease in permeability, which hinders the material exchange of IVDs and affects their normal physiological function [16,17]. However, CEP is often overlooked in physical simulations and simplified as boundary conditions for IVD [18]. Moreover, it is currently unknown whether CEP degeneration will exacerbate IVDD by altering its biomechanical properties.

The study of IVD includes cell experiments, animal models, organ cultures, and cadaver sample experiments. However, given the intricate structure of the IVD and the limited availability of clinical samples, current methodologies have limitations. Exploring novel research pathways to understand the etiology of IVDD is imperative. Finite element analysis (FEA) offers a mean to simulate or analyze the behavior of structures or components. FEA can combine complex material properties with geometric shapes to model the multiaxial mechanics of IVDs and tackle research problems that cannot be experimentally tested [19]. Within FEA, the biphasic swelling theory offers significant advantages by describing the behavior of both fluid and solid phases and accounting for material swelling induced by electrochemical potential differences [20]. Additionally, the solid mixture theory connects the tissue's mechanical behavior to its individual components. These features make the biphasic swelling mixture theory an optimal framework for characterizing the biomechanical behavior of intervertebral soft tissues, thereby enhancing the accuracy and predictive power of FEA in IVD research [21]. In recent years, the finite element analysis of IVD has mostly focused on the impact of NP and AF degeneration on IVD [22].However, from a biomechanical perspective, it remains uncertain whether the occurrence of cartilage endplate degenerations (EPDs) will lead to intervertebral disc degeneration.

This study aims to employing swelling-based biphasic FEA models to analyze the biomechanical behavior of different IVD conditions under quasi-static loading conditions. Subsequent daily load analyses will approximate the biomechanical behavior under realistic conditions, offering new insights into the etiology for IVDD. This study not only reveals the critical role of CEP degeneration in IVDD, providing new insights and theoretical foundations that advance our understanding of the mechanisms of disc degeneration, but also offers technical and data support for the improvement of artificial discs and the design of patient-specific disc prostheses.

2 Models and methods

2.1 Computational models

In our study, we created four unique intervertebral disc models to represent different stages of endplate and disc degeneration, as detailed in Table 1(The detailed information on model assumptions, parameter selection, and verification process can be found in the supplementary documents S1).Table 1 Experimental grouping.

Table 1	Group1	Group2	Group3	Group4	
cartilaginous endplate	No degeneration	degeneration	No degeneration	degeneration	
intervertebral disc	Mild degeneration or health	Moderate to severe degeneration	

Utilizing Hypermesh software (Altair Engineering, Inc., Troy, Michigan, USA) alongside a tailor-made Python script, we constructed the initial volumetric model of IVD. The process began with generating a two-dimensional quad mesh for the transverse cross-section, which served as a foundation for creating three-dimensional 8-node hexahedral elements. The geometry of the IVD model (as illustrated in Fig. 1) was inspired by the average dimensions of seven human L4-L5 IVDs [23]. The histological and MR data [9] were consistent with the specifications for the CEP and bony endplates (BEPs), with CEP comprising four layers at a height of approximately 600 μm, and the BEP consisting of three layers at a thickness of 800 μm. BEP fibers predominantly align parallel to the disc. Subsequently, the order of nodes in each hexahedral element belonging to the AF was adjusted, enabling the definition of collagen fiber orientation through local coordinate frames. The AF was segmented into 20 concentric lamellae across 12 vertical layers, featuring a typical cross-ply arrangement of collagen fibers. The completed model boasts 69,368 elements, including 8976 in the NP and 23,040 in the AF, with all neighboring components unified at their interfaces. As depicted in Fig. 1A, the model's mean height is 11.1 mm, with a total volume of 20,682 mm3, and the NP contributes to 36 % of the volume. The superior surface area measured 1919.5 mm2, with the NP accounting for 40 % [12] of this volume."The centroid of the NP is offset toward the rear of the IVD model by 10 % [24], as shown in Fig. 1B and C.Fig. 1 Schematic illustrations of the IVD include: an isometric view (A) depicting the disc with the height of 11.1 mm; an axial section (B), providing a profile view; and a sagittal section (C), detailing the AF, NP, CEPs, and VBs.

Fig. 1

2.2 Constitutive models

The vertebral bodies (VBs), CEP, AF, and NP all utilize structure-based solid mixture constitutive laws. These laws effectively link behavior with its constituents. The constitutive behavior of VBs was modeled using a compressible neo-Hookean framework. This approach can express the hyperelastic strain-energy density as follows:(1) W=μ2(I1−3)−μlnJ+λ2(lnJ)2,

Here, the determinant of the deformation gradient tensor(F), denoted by J, the first invariant of the right Cauchy–Green deformation tensor (C=FTF) is represented by I1, Young's modulus (E) and Poisson's ratio (ν) are related to the Lamé constants λ and μ.(2) λ=Ev(1+v)(1−2v),μ=E2(1+v).

In this work, the VBs were modeled as a biphasic neo-Hookean materials with the Young's modulus E = 12 GPa and Poisson's ratio ν = 0.3 [25]. In addition, a constant permeability of k0 = 5 mm4/N was chosen [26], in order to allow fluids to flow through.

The SCEP, ICEP, AF, and NP, each comprise an isotropic, compressible, and hyperelastic solid matrix, characterized by the Holmes-Mow model. The expression for strain-energy density within the Holmes-Mow matrix framework is as follows:(3) Ψm(I1,I2,J)=λ+2μ4βm(eβmλ+2μ[−(λ−2μ)(I1−3)+λ(I2−3)−(λ+2μ)ln(J2)]−1)

The second invariant of the right Cauchy-Green deformation tensor is I2, Characterizing the rate of the exponential stiffening is βm. These parameters, along with Young's modulus (Em), Poisson's ratio (νm) and other material constants, are specified in Table 2 [27,28].

The solid matrix of the IVD also incorporates an isotropic and deformation-dependent permeability behavior, as described by the Holmes-Mow permeability model [29].(4) k=k0(J+φ0w−1φ0w)2eM2(J2−1)I,

where k is the permeability tensor, k0 is the initial isotropic hydraulic permeability in the reference configuration, φ0w is the water volume fraction, and M is the exponential coefficient of permeability. Both M and k0 were determined by fitting to confined compressive experimental data [[27], [28], [30]]. (Table 2, Table 3).Table 2 Material properties used in the constitutive models were obtained from experimental data in existing studies [23,[27], [28], [30], [31], [32], [33]]. SCEP (superior cartilaginous endplate) and ICEP (inferior cartilaginous endplate) are defined accordingly.

Table 2Property	Symbol	SCEP	ICEP	AF	NP	
Initial water fraction (unitless)	φ0w	0.6	0.6	0.7	0.8	
Nonlinearity of permeability (unitless)	M	5.2	3.3	2.7 (5.7)a	1.92	
Initial permeability (10−4 mm4/N)	k0	8.24	2.87	16 (47)a	5.5	
Poisson's ratio (unitless)	νm	0.18	0.18	0.16 (0.24)a	0.24	
Osmotic coefficient (unitless)	Φ	1	1	1	1	
Modulus of fibers (MPa)	ξ1/Ef	7.01b	7.01b	3.0 (15.6)a	N.A.	
Rate of fiber stiffening (unitless)	βc/ βf	2.88c	2.88c	4	N.A.	
Rate of matrix stiffening (unitless)	βm	0.29	0.29	1.5 (3.4)a	0.95	
Initial fixed charge density (mmol/L)	cf0	−250	−450	−300 (−100)a	−400	
Critical stretch square (unitless)	I0	N.A.	N.A.	1.14 (1.06)a	N.A.	
Modulus of matrix (MPa)	Em	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 ξ1 represents the initial elastic modulus of fibers in the ICEP and SCEP.

c βc represents the nonlinearity in fiber strain-energy density within the ICEP and SCEP.

Table 3 Initial material property values were established based on experimental data (literature) [32,33,34,35].

Table 3Property	Symbol	Group1	Group2	Group3	Group4	
Nonlinearity of SCEP (unitless)	M	5.99	4.49	5.99	4.49	
Nonlinearity of ICEP (unitless)	M	3.79	3.08	3.79	3.08	
Permeability of SCEP (10−4 mm4/N)	k0	8.24	3.06	8.24	3.06	
Permeability of ICEP (10−4 mm4/N)	k0	2.87	1.62	2.87	1.62	
Initial fixed density of SCEP (mmol/L)	cf0	110	40	110	40	
Initial fixed density of ICEP (mmol/L)	cf0	110	40	110	40	
Fixed charge density of AF (mmol/L)	cf0	−300 (−100)a	−300 (−100)a	−150 (−100)a	−150 (−100)a	
Fixed charge density of NP(mmol/L)	cf0	−350	−350	−125	−125	
a The parameter varies linearly across the lamellae of the AF, with the inner lamella value in parentheses.

Within the framework of biphasic-swelling theory, the glycosaminoglycan (GAG) content is correlated with the fixed charge density. The fixed charge density, prior to deformation, is established in the reference configuration as follows:(5) cf0=cgZcMc,

where cg is the GAG content in water with the dimension mg/ml, and Mc and Zc are the molecular weight and number of charges of chondroitin-6-sulfate disaccharide. The latter two parameters were chosen as Mc=502.5 gram and Zc=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).

Upon applying external loads, the instantaneous fixed charge density changes according to the deformation Jacobian (J), initial water volume fraction (φsw) and initial fixed charge density (cf0)(6) cf=cf0φ0wJ−1+φ0w.

As deformation occurs, changes in the fixed charge density lead to variations in osmotic pressure.(7) po=RTΦ(cf2+cb2−cb),

The osmotic pressure equation incorporates the temperature (T), gas constant (R), osmotic coefficient(Φ), and the osmolarity of the external bath (cb).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 (pf) and solid stress (σs), the osmotic pressure (po) also contributes to the total Cauchy stress (σ)(8) σ=−(pf+po)I+σs,

where 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 expressed as follows:(9) Ψf(n,In)=ξ(n)(In−1)β(n),

The instantaneous square of the fiber stretch is represented by In , i.e., In=N⋅C⋅N, with N being the unit vector along the fiber axis in the reference configuration. n is the unit vector of fiber in the current reference, i.e., n=F⋅N/λn. In spherical coordinate system, the unit vector (n) measured by the longitudinal (φ) and latitudinal (θ) coordinates. With their help, the material parameters β(n) and λ(n) can be defined as(10) ξ(n)=(cos2θsin2φξ12+sin2θsin2φξ22+cos2φξ32)−12,

(11) ξ(n)=(cos2θsin2φβ12+sin2θsin2φβ22+cos2φβ32)−12,

where ξi and βi (i=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 mainly oriented in the planes that are parallel to VBs and NP. Fibers along the lateral and anterior-posterior directions share the same modulus (ξ1=ξ2). Along the axial direction of the IVD, much smaller elastic modulus was assumed, e.g., ξ3=0.1ξ1. The fiber nonlinearity was assumed to be independent of coordinate directions, e.g., β(n)=βc, with βc being a constant. Following a regression analysis by fitting to tensile experimental data, the modulus ξ1=ξ2=10ξ3=7.01 MPa and the fiber nonlinearity βc=2.88 were obtained [30](Table 2).

The collagen fibers are modeled using a toe-linear stress-stretch constitutive law to capture their nonlinear mechanical properties in the AF. In terms of the square of the fiber stretch, the strain-energy density can be ued to expresse (In)(12) Ψf={0,In<1,ξβf(In−1)βf,1≤In≤I0,B(In−I0)−Ef(In−I0)+ξβf(I0−1)βf,I0≤In,

where I0 represent for the critical square of the fiber stretch (λ02) corresponding to the point that separates the linear and the toe regions, and ξ and B are both functions of the linear-region elastic modulus (Ef), I0 and the toe-region power-law coefficient (βf)(13) ξ=Ef4(βf−1)I0−32(I0−1)2−βf,B=ξ(I0−1)βf−1+Ef2In−12.

Parameters of the constitutive law such as Ef, I0 and βf were calibrated in terms of the experimental stress-stretch curve of single-lamella AF samples [31]. From the inner to the outer lamellae, Ef increases linearly, whereas I0 and β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 la, the average collagen fiber direction of the innermost lamella is θ = 40° with the transverse cross-section plane, and the element size of the inner boundary curve, is lin. These three parameters approximately satisfy(14) tanθ=lalin

2.3 Quasi-static tests, daily loads and boundary conditions

Before conducting quantitative analysis of the fluid dependent response of IVDs, free-swelling, stress-relaxation, creep, and slow-ramp are selected to verify the predictive ability of the established IVDs.

After testing four quasi-static testing, this paper simulated daily loads on the model. Daily loads were used to divide a normal day into two stages to further observe the IVD height and intramedullary pressure under different pressures [36]: there are 8 h of low level loading and 16 h of medium high level loading, respectively. Under low load cycles, the horizontal pressure of the load is 350N, while under medium load, the horizontal pressure is 1000N (representing moderate daily activities). The change in pressure load is generally completed within 10 s.

Two rigid points, named "bc-dot" and "loading-dot," were established to ensure precise application of compressive load, connecting to nodes on the superior and inferior VBs of the IVD.In the execution of all protocols, the bc-dot remains entirely stationary, while the loading-dot is restricted to axial movement along the IVD. Moreover, the outer surface IVDs is consistently subjected to a condition of zero fluidic pressure.

2.4 Calculation of intervertebral disc height and intramedullary pressure

2.4.1 Definition of the initial height

(15) Hinitial=Ssagtiallap

Where Ssagtial represents the sagittal cross-sectional area of the AF and NP complex, lap is the length of the anterior posterior diameter of the AF and NP. After measurement, the sagittal cross-sectional area of the AF and NP complex is 478.1 mm2, the length of the front and rear diameters is 43.1 mm, the height of the IVD is 11.1 mm.

2.4.2 Intervertebral disc height in hydrated state

(16) Hhydrated=Hinitial+dloading−dot

Where dloading−dot is the vertical displacement of the loading point.

2.4.3 Intranuclear pressure

During the hydration process, the intervertebral soft tissue will absorb water and swelling Simply averaging the intramedullary pressure of all units in the NP region is not the best practice. This article uses an average algorithm weighted by unit volume to calculate intramedullary pressure,(17) pnp=∑i=1Nwipi,

Where pnp is the average internal pressure of the NP, i is volume weight in units, wi=vi/∑i=1Nvi, pi is the volume of unit σ in a deformed configuration.

2.5 Validation

In order to validate the model, the data generated from quasi-static experiments are validated against the data generated from experiments [[37], [38], [39]]. Specifically, relative normalized mean square error (RNMSE) is used to quantify the difference between model predicted data and experimental data (S1.3). The experimental source literature on intervertebral disc nucleus pulposus pressure and intervertebral disc height are labeled in specific figure annotations.

2.6 Data analysis method

Simulations were performed using FEBio [40]. Additionally, a tailored Python script was created for some processing tasks.

The double Voigt model characterizes the nominal strain-time relationship derived from FEA (free-swelling and creep) [41]:(18) ε=σ0Ee+σ0E1c(1−e−tτ1c)+σ0E2c(1−e−tτ2c)

where σ0/Ee represents the completely elastic response, τc 1 and Ec 1 represent short-term response and τc 2 and Ec 2represent long-term response.

Appling the double Maxwell model to accurately fit the temporal stress behavior (stress-relaxation test) [42]:(19) σ=E∞+ε0E1se−tτ1s+ε0E2se−tτ2s

The equation is divided into the long-term (τs 2/Es 2), short-term (τs 1/Es 1) and residual (ε0/E∞) responses. Engineering equilibrium time (τ95 %) refers to the duration required to reach 95 % of the final response at the end of a simulation.

3 Result

3.1 Model validation

Comparing the model with existing literature data, it was found that the model has a high degree of consistency with the data (Fig. 2). Because the permeability and fixed charge density during endplate degeneration do not affect the result of hydration, the IVD height and intranuclear pressure of Group 1 and 2, and Group 3 and 4 are equal. From the comparison results (Fig. 2A and B), the simulation data are relatively close to the previous data, and the predicted results in this article are acceptable.Fig. 2 Comparison of simulation data with literature data to verify the model. A: The height of the IVD after free-swelling. The bar chart shows the disc height after 50N compression and hydration, compared with the literature value [43]. B: NP pressure after free-swelling (NP pressure is determined as a weighted average, using element volumes as weights, Refer to equation (17) for further details [39,[44], [45], [46]]).(C)The nominal strain of slow-ramp. (D) The nominal strain of creep, (E) Reaction Force in the stress-relaxation test. (C–E) Exp (mean) in the figure represents the literature data, while Exp (95 %) represents the confidence interval [37,38,47].

Fig. 2

3.2 Free-swelling

During free-swelling, the displacements of the IVD in each group were 3.34 mm, 3.31 mm, 1.59 mm, and 1.56 mm, respectively, with the compressive strains are of 30.2 %, 30.0 %, 14.5 %, and 14.2 %. Under these operating conditions, the contribution of the short-term response plays a dominant role (Fig. 4A), and the proportions of the short-term and long-term responses are not significantly related to changes in the endplate permeability. The nominal strain-time curve was numerically fitted using formula (18) (Fig. 3A), and the double Voigt model fitted the curve very well (R2>0.995, Table 4). Meanwhile, Table 4 shows that τ1candτ2c in the IVDD model were shorter than those in the normal model, while endplate degeneration led to the elongation of τ1candτ2c.Fig. 3 Strain-time curves (A, C) for free-swelling and creep tests; stress-strain curves (B, D) for slow-ramp and stress-relaxation tests.

Fig. 3

Fig. 4 Histograms depict each response in the (A) free-swelling, (B) creep, and (C) stress-relaxation tests.

Fig. 4

Table 4 Parameters obtained from free-swelling using the double Voigt model.

Table 4	σ0/Ee	σ0/E1	τ1c (hr)	σ0/E2	τ2c (hr)	R2	τ95% (hr)	
Group1	0	18.8	0.48	11.4	3.49	0.998	7.28	
Group2	0	18.1	0.55	11.9	4.34	0.998	9.26	
Group3	0	6.9	0.42	7.3	2.10	0.998	7.68	
Group4	0	7.6	0.53	6.3	3.13	0.998	10.33	
a: σ0/Ee:Complete elastic response. σ0/E1:Short-term response. τ1c (hr): Short-term response time constant. τ2c (hr): Long-term response. τ2c (hr): Long-term response time constant. R2: Goodness of fit. τ95% (hr): Time to reach 95 % of the final response.

3.3 Slow-ramp

In slow-ramp tests, as shown in Table 5 and Fig. 3B, the average nominal strains of the IVD in each group were approximately 24.2 %, 23.7 %, 27.4 %, and 26.3 %. Fig. 2C compares the predicted values and experimental results of various models under slow compression conditions. From the image, the predicted results of the model built in this paper are all within the confidence intervals of the experimental results, indicating that the model has high reliability for slow-ramp conditions.Table 5 The nominal compressive strain table at the time of Slow-ramp and Creep results (relative to the height of the hydrated equilibrium intervertebral disc, the creep test corresponding to this table corresponds to a compressive load of 1000N).

Table 5Slow-ramp	Group 1	Group 2	Group 3	Group 4	
Compressive Strain (%)	24.2	23.7	27.4	26.3	
Creep	
Compressive Strain (%)	40.8	39.4	52.7	50.0	

3.4 Creep

In creep test, the nominal compressive strains of each IVD group were approximately 40.8 %, 39.4 %, 52.7 %, and 50.0 % (Table 5). Using the dual Voigt model to fit the strain-time curve (Fig. 3C), it can be found that the model fits well with the data (R2>0.995, Table 6). Meanwhile, comparing the predicted values of the model with the experimental results, it is evident that the (Fig. 2D) model also has high reliability in creep testing. However, the specimens used in the experiment were degenerative IVD, so Group 1 and Group 2 are not included in the comparison. Finally, as shown in Table 6, the elastic response, long-term response, and short-term response all increase with intervertebral disc degeneration. Fig. 4B further provides the proportion of strain caused by each response, and it was found that there was no significant correlation with intervertebral disc and endplate degeneration.Table 6 Parameters obtained from creep using the double Voigt model.

Table 6	σ0/Ee	τ1c (hr)	σ0/E1	τ2c (hr)	σ0/E2	R2	τ95% (hr)	
Compressive load of 1000N	
Group1	12.0	0.17	6.7	5.84	21.4	0.997	17.95	
Group2	12.0	0.19	6.7	6.81	20.0	0.997	20.02	
Group3	14.0	0.24	9.8	6.75	28.0	0.997	19.86	
Group4	14.0	0.22	8.9	8.17	26.6	0.996	21.99	
a: σ0/Ee:Complete elastic response. σ0/E1:Short-term response. τ1c (hr): Short-term response time constant. σ0/E2: Long-term response. τ2c (hr): Long-term response time constant. R2: Goodness of fit. τ95% (hr): Time to reach 95 % of the final response.

3.5 Stress-relaxation

Fig. 2E shows the predicted results of the model and the experimental results. From the results, the model also has good reliability in the relaxation experiment. However, the specimen used in the experiment was a degenerative intervertebral disc, so Group 1 and Group 2 were not included in the comparison. The dual Maxwell model was used to numerically fit the nominal stress time curve shown (Fig. 3D), and all four curves were well fitted (R2>0.99, Table 7). From Table 7 showed that the residual component, long-term response, and short-term response all decreased with the increase of degradation, while the change in strain corresponding to endplate degeneration is relatively small. The further contribution data of Fig. 4C can indicate that there is no significant correlation between the contribution ratio and the degeneration of the endplate.Table 7 Parameters obtained from stress-relaxation using the double Maxwell model.

Table 7	E1sε0 (MPa)	τ1s (hr)	E2sε0 (MPa)	τ2s (hr)	R2	E∞sε0 (MPa)	τ95% (hr)	
5 % compressive norminal strain	
Group1	82.4	0.710	86.2	0.034	0.996	101.6	7.02	
Group2	81.7	0.826	88.0	0.040	0.995	100.8	8.63	
Group3	47.5	0.492	54.4	0.028	0.995	39.6	11.05	
Group4	43.5	0.629	56.2	0.036	0.996	39.4	12.06	
aE∞sε0 (MPa):Residual Component. E1sε0 (MPa):Short-term response. τ1s (hr): Short-term response time constant. E2sε0 (MPa): Long-term response. τ2s (hr): Long-term response time constant. R2: Goodness of fit. τ95% (hr): Time to reach 95 % of the final response.

3.6 Daily loads

According to a previous article [36], the study of the daily loads of the models are divided into two angles, compressive displacement and nominal strain (Fig. 5A and B). Both results are relative to the IVD height at the conclusion of the first daily load. As shown in the figure, both IVDD and EPD can cause the height of the IVD to be lower than the normal model. Finally, the entire process of NP pressure was fitted, and the IVDD leads to a significant decrease in pulposus pressure in the day night cycle, while EPD leads to higher pulposus pressure in the daily working state. (Fig. 5C).Fig. 5 Daily load: Curve graph of intervertebral disc height over time under the second daily load: (A) compressive displacement; (B) Nominal-strain. C: NP pressure.

Fig. 5

4 Discussion

This article used a FEA model to analyze the quasi-static mechanical behavior of IVDs under different endplate permeabilities and simulated the loading of IVD in daily human life. Our aim was to elucidate the biomechanical implications of EPD on IVD through finite element analysis.

In fluid-dependent viscoelastic responses, the content of GAG in the soft tissue of IVD is an important factor affecting fluid dependent viscoelastic response, followed by the permeability of CEP. In free-swelling scenarios, both long-term and short-term responses diminish with the IVDD. Conversely, EPD prolongs the time of each response. Variations in GAG content in degenerative tissues account for the differential responses observed between degenerative and healthy IVD. Furthermore, changes in CEP permeability increase response times, indicating that CEP degradation alters the fluid pathways, potentially facilitating fluid ingress into the IVD from the AF [39]. However, the changes caused by CEP degeneration may affect the long-term stability of IVD and increase the risk of degeneration. Specifically, reducing fluid flow through CEP can disrupt the uniform stress distribution within the IVD [48].During the slow-ramp test, the compressive load applied by the quasi-static protocol determines that the influence of CEP permeability can be overlooked. Consequently, variations in response are predominantly due to alterations in the GAG content within NP and AF. These changes affect the response by influencing the initial stiffness and deformation of IVD within the neutral zone [49].During the creep test, all responses exhibit a correlation with the IVDD, and individual contributions to total deformation are independent [50,51]. Within the four models, the AF primarily influences the long-term response, contributing to more than 50 % of it, marking a distinct contrast to the predominance of short-term response observed in free-swelling conditions. Earlier research has identified that the NP and CEPs chiefly govern short-term responses, whereas the AF is more influential in long-term responses [52]. Additionally, permeability's flow direction-dependency offers greater resistance to fluid outflow than inflow [53].This means that in creep test, the applied force overwhelms the osmotic pressure and quickly drives the fluid out of the system. Howvere, osmotic pressure is the sole force, directing fluid flow from the CEP into the IVD in free swelling. In fact, the hydraulic permeability of the fibrous ring is 10 times that of the cartilage endplate [27]. When the fluids of CEP and IVD remain stable, the remaining liquid will continue to exchange through the annulus fibrosis. This explains why AF accounts for over 50 % of the long-term response in creep. EPD results in diminished fluid discharge, subsequently lowering nominal compression. These changes affect the balance of pressure and the distribution of stress, leading to reduced mechanical stability of the IVD under long-term loading conditions and an increased risk of structural damage [54]. In stress-relaxation tests, stress is negatively correlated with the degeneration of intervertebral disc soft tissues. This suggests that the degeneration process weakens the aggregate modulus and swelling pressure of the disc soft tissues, leading to a more rapid decrease in stress under the same deformation conditions [55].Nevertheless, EPD does not significantly alter stress-relaxation outcomes. Although the primary fluid flow is anticipated to occur axially through the CEP, alterations in CEP properties could be offset by an enhanced radial flow through the AF [23].Furthermore, reactive force rapidly decreases to a stable value, leading to minor changes in the overall disc load response over time. In summary, high GAG content in IVD ensures sufficient fluid retention and swelling pressure, providing compressive strength and maintaining intervertebral disc height [56]. On the contrary, the increased permeability of degenerated CEP allows for excessive fluid outflow, reducing the ability of the intervertebral disc to withstand compressive loads and altering its viscoelastic properties [39].

To more accurately mimic the normal physiological activities of human body, we also examined the behavior of IVD under conditions that simulate daily loading.Because reference [36] only analyzes the second daily load, this article also analyzes the second daily load. In alignment with prior human studies, more than 50 % of the height loss is observed within the first hour of pressure increase, and approximately 80 % occurs within the first 3 h following compression [57]. In daily loads, the height of the IVD and the pressure inside the NP decrease with the degeneration of the IVD soft tissue. The degenerated NP loses its ability to bind water, resulting in a decrease in pressure inside the NP and intervertebral disc height [58]. The pressure inside the NP increases with the EPD, while the height of the IVD decreases. With EPD, the pressure within the NP escalates, accompanying a reduction in intervertebral disc height. Creep tests have previously shown that the EPD group exhibits a lower nominal compressive displacement compared to the control group. However, under daily loading conditions, the EPD group's intervertebral disc height is observed to be lower. This discrepancy might stem from the predominant role of disc hydration during a typical 24-h activity cycle. Reduced endplate permeability diminishes fluid outflow, leading to increased NP pressure. This decrease in disc height alongside elevated NP pressure are key indicators of intervertebral disc degeneration, suggesting that endplate degeneration could precipitate overall disc degeneration [59].

This study acknowledges certain limitations. Firstly, this study constructs a finite element model using image and assumes that it is in a state of osmotic dehydration.However, these images were pre-strained due to swelling, which could lead to discrepancies in the simulation results. Secondly, modeling fluid flow within an intervertebral disc under load encompasses a sophisticated three-dimensional challenge, characterized by anisotropy, non-uniformity, strain-dependent permeability, varying solid matrix modulus, along with heterogeneous fixed charge density and osmotic pressure [60]. This complexity necessitates the formulation of more refined and intricate constitutive models. Finally, we controlled the GAG content in CEP by setting the charge density. We only used one experimental data, but in fact, the content and structure of CEP are constantly changing during the degeneration process, which can affect biomechanical behavior.

Future research should consider further evaluating the biomechanical behavior of CEP and IVD under dynamic and cyclic loading conditions. In addition, the current intervertebral disc modeling process can further introduce ionic diffusion and convection constitutive relationships to explore the coupling mechanism between biochemical reactions and mechanical behavior of the IVD.

5 Conclusion

In summary, the following conclusions emerge from the current research. Firstly, the content of GAG in the soft tissue of IVD is an important factor affecting fluid dependent viscoelastic response, followed by the permeability of CEP. Secondly, degeneration of the CEP may alter the height and intramedullary pressure of IVD, potentially accelerating disc degeneration.

Data availability

Data included in article/supp. material/referenced in article.

Funding statement

This work was financially supported by a grant from 10.13039/501100014996 Health Commission of Zhejiang Province (2022KY997 ). This work was financially supported by a grant from The Medical and Health Science and Technology Program of Zhejiang Province (2020KY748 ). Project of Zhejiang provincial plan for TCM science and technology (2022ZB238 ).

CRediT authorship contribution statement

Yujun Zhang: Writing – review & editing, Writing – original draft. Yanli Pan: Software, Methodology. Xinning Mao: Methodology, Investigation. Du He: Writing – original draft, Methodology. Liangping Zhang: Investigation, Formal analysis. Wei Cheng: Investigation, Formal analysis, Data curation. Chengyue Zhu: Methodology, Investigation, Data curation. Hang Zhu: Supervision, Methodology. Wei Zhang: Methodology, Investigation. HongTing Jin: Formal analysis, Data curation, Conceptualization. Hao Pan: Funding acquisition, Formal analysis, Conceptualization. Dong Wang: Methodology, Formal analysis, Conceptualization.

Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:Hao Pan reports financial support was provided by Zhejiang Provincial Health Commission (2022KY997 ). Hao Pan reports financial support was provided by Project of Zhejiang provincial plan for TCM science and technology (2022ZB238 ). Hao Pan reports financial support was provided by The Medical and Health Science and Technology Program of Zhejiang Province (2020KY748 ). If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is the Supplementary data to this article.Multimedia component 1

Multimedia component 1

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e37524.
==== Refs
References

1 Cieza A. Global estimates of the need for rehabilitation based on the global burden of disease study 2019: a systematic analysis for the global burden of disease study 2019 Lancet (London, England) 396 2021 2006 2017 10.1016/s0140-6736(20)32340-0 33275908
2 Wong J. Nutrient supply and nucleus pulposus cell function: effects of the transport properties of the cartilage endplate and potential implications for intradiscal biologic therapy Osteoarthritis Cartilage 27 2019 956 964 10.1016/j.joca.2019.01.013 30721733
3 Fields A. Ballatori A. Liebenberg E. Lotz J. Contribution of the endplates to disc degeneration Current molecular biology reports 4 2018 151 160 10.1007/s40610-018-0105-y 30546999
4 Mohd Isa I. Intervertebral disc degeneration: biomaterials and tissue engineering strategies toward precision medicine Adv. Healthcare Mater. 11 2022 e2102530 10.1002/adhm.202102530
5 Buckwalter J. Aging and degeneration of the human intervertebral disc Spine 20 1995 1307 1314 10.1097/00007632-199506000-00022 7660243
6 Cassidy J. Hiltner A. Baer E. Hierarchical structure of the intervertebral disc Connect. Tissue Res. 23 1989 75 88 10.3109/03008208909103905 2632144
7 Lyu F. Painful intervertebral disc degeneration and inflammation: from laboratory evidence to clinical interventions Bone research 9 2021 7 10.1038/s41413-020-00125-x 33514693
8 Binch A. Fitzgerald J. Growney E. Barry F. Cell-based strategies for IVD repair: clinical progress and translational obstacles Nat. Rev. Rheumatol. 17 2021 158 175 10.1038/s41584-020-00568-w 33526926
9 Moon S. Evaluation of intervertebral disc cartilaginous endplate structure using magnetic resonance imaging Eur. Spine J. : official publication of the European Spine Society, the European Spinal Deformity Society, and the European Section of the Cervical Spine Research Society 22 2013 1820 1828 10.1007/s00586-013-2798-1
10 Jonathan K. Lixia Z. Thierry B. Hai Y. The region-dependent biphasic viscoelastic properties of human temporomandibular joint discs under confined compression J. Biomech. 43 2010 10.1016/j.jbiomech.2010.01.020
11 Katherine B.C. Cartilaginous endplates: a comprehensive review on a neglected structure in intervertebral disc research JOR Spine 6 2023 10.1002/jsp2.1294
12 Newell N. Biomechanics of the human intervertebral disc: a review of testing techniques and results J. Mech. Behav. Biomed. Mater. 69 2017 420 434 10.1016/j.jmbbm.2017.01.037 28262607
13 Ruiz Wills C. Foata B. González Ballester M. Karppinen J. Noailly J. Theoretical explorations generate new hypotheses about the role of the cartilage endplate in early intervertebral disk degeneration Front. Physiol. 9 2018 1210 10.3389/fphys.2018.01210 30283342
14 Grant M. Human cartilaginous endplate degeneration is induced by calcium and the extracellular calcium-sensing receptor in the intervertebral disc Eur. Cell. Mater. 32 2016 137 151 10.22203/ecm.v032a09 27452962
15 Antoniou J. The human lumbar endplate. Evidence of changes in biosynthesis and denaturation of the extracellular matrix with growth, maturation, aging, and degeneration Spine 21 1996 1153 1161 10.1097/00007632-199605150-00006 8727189
16 Maroudas A. Stockwell R. Nachemson A. Urban J. Factors involved in the nutrition of the human lumbar intervertebral disc: cellularity and diffusion of glucose in vitro J. Anat. 120 1975 113 130 1184452
17 Roberts S. Urban J. Evans H. Eisenstein S. Transport properties of the human cartilage endplate in relation to its composition and calcification Spine 21 1996 415 420 10.1097/00007632-199602150-00003 8658243
18 Hassan C. Lee W. Komatsu D. Qin Y. Evaluation of nucleus pulposus fluid velocity and pressure alteration induced by cartilage endplate sclerosis using a poro-elastic finite element analysis Biomech. Model. Mechanobiol. 20 2021 281 291 10.1007/s10237-020-01383-8 32949306
19 Hu B. Lv X. Chen S. Shao Z. Application of finite element analysis for investigation of intervertebral disc degeneration: from laboratory to clinic Current medical science 39 2019 7 15 10.1007/s11596-019-1993-7 30868485
20 A Pérez d.P. B C. M D. An accurate finite element model of the cervical spine under quasi-static loading J. Biomech. 41 2007 10.1016/j.jbiomech.2007.10.012
21 Y S. J M H. C C.v.D. K I. A biochemical/biophysical 3D FE intervertebral disc model Biomech. Model. Mechanobiol. 9 2010 10.1007/s10237-010-0203-0
22 Ayturk U. Gadomski B. Schuldt D. Patel V. Puttlitz C. Modeling degenerative disk disease in the lumbar spine: a combined experimental, constitutive, and computational approach J. Biomech. Eng. 134 2012 101003 10.1115/1.4007632
23 Jacobs N. Cortes D. Peloquin J. Vresilovic E. Elliott D. Validation and application of an intervertebral disc finite element model utilizing independently constructed tissue-level constitutive formulations that are nonlinear, anisotropic, and time-dependent J. Biomech. 47 2014 2540 2546 10.1016/j.jbiomech.2014.06.008 24998992
24 O'Connell G. Vresilovic E. Elliott D. Comparison of animals used in disc research to human lumbar disc geometry Spine 32 2007 328 333 10.1097/01.brs.0000253961.40910.c1 17268264
25 Williams J.R. Natarajan R.N. Andersson G.B. Inclusion of regional poroelastic material properties better predicts biomechanical behavior of lumbar discs subjected to dynamic loading J. Biomech. 40 2007 1981 1987 10.1016/j.jbiomech.2006.09.022 17156786
26 Ferguson S. Ito K. Nolte L. Fluid flow and convective transport of solutes within the intervertebral disc J. Biomech. 37 2004 213 221 10.1016/s0021-9290(03)00250-1 14706324
27 Cortes D. Jacobs N. DeLucca J. Elliott D. Elastic, permeability and swelling properties of human intervertebral disc tissues: a benchmark for tissue engineering J. Biomech. 47 2014 2088 2094 10.1016/j.jbiomech.2013.12.021 24438768
28 Cortes D. Elliott D. Extra-fibrillar matrix mechanics of annulus fibrosus in tension and compression Biomech. Model. Mechanobiol. 11 2012 781 790 10.1007/s10237-011-0351-x 21964839
29 Holmes M. Mow V. The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration J. Biomech. 23 1990 1145 1156 10.1016/0021-9290(90)90007-p 2277049
30 DeLucca J. Human cartilage endplate permeability varies with degeneration and intervertebral disc site J. Biomech. 49 2016 550 557 10.1016/j.jbiomech.2016.01.007 26874969
31 Holzapfel G. Schulze-Bauer C. Feigl G. Regitnig P. Single lamellar mechanics of the human lumbar anulus fibrosus Biomech. Model. Mechanobiol. 3 2005 125 140 10.1007/s10237-004-0053-8 15778871
32 Urban J.P.G. Maroudas A. The measurement of fixed charged density in the intervertebral disc BBA - General Subjects 586 1979 166 178
33 Urban J. McMullin J. Swelling pressure of the lumbar intervertebral discs: influence of age, spinal level, composition, and degeneration Spine 13 1988 179 187 10.1097/00007632-198802000-00009 3406838
34 Wu Y. Quantifying baseline fixed charge density in healthy human cartilage endplate: a two-point electrical conductivity method Spine 42 2017 E1002 E1009 10.1097/brs.0000000000002061 28699925
35 Roberts S. Menage J. Urban J. Biochemical and structural properties of the cartilage end-plate and its relation to the intervertebral disc Spine 14 1989 166 174 10.1097/00007632-198902000-00005 2922637
36 Schmidt H. Shirazi-Adl A. Galbusera F. Wilke H. Response analysis of the lumbar spine during regular daily activities--a finite element analysis J. Biomech. 43 2010 1849 1856 10.1016/j.jbiomech.2010.03.035 20394933
37 Grace D O.C. Nathan T J. Sounok S. Edward J.V. Dawn M.E. Axial creep loading and unloaded recovery of the human intervertebral disc and the effect of degeneration J. Mech. Behav. Biomed. Mater. 4 2011 10.1016/j.jmbbm.2011.02.002
38 Yoder J.H. Internal three-dimensional strains in human intervertebral discs under axial compression quantified noninvasively by magnetic resonance imaging and image registration J. Biomech. Eng. 136 2014 1110081 1110089 10.1115/1.4028250 25109533
39 Sun Z. Sun Y. Lu T. Li J. Mi C. A swelling-based biphasic analysis on the quasi-static biomechanical behaviors of healthy and degenerative intervertebral discs Comput. Methods Progr. Biomed. 235 2023 107513 10.1016/j.cmpb.2023.107513
40 Maas S. Ellis B. Ateshian G. Weiss J. FEBio: finite elements for biomechanics J. Biomech. Eng. 134 2012 011005 10.1115/1.4005694
41 Vergroesen P. Emanuel K. Peeters M. Kingma I. Smit T. Are axial intervertebral disc biomechanics determined by osmosis? J. Biomech. 70 2018 4 9 10.1016/j.jbiomech.2017.04.027 28579261
42 Abedi R. Shayegan S. Ataee G. Fatouraee N. Viscoelastic modeling of ovine cervical intervertebral disc through stress-relaxation Constant Strain Rate and Dynamic Loading Tests 2019
43 Newman H. DeLucca J. Peloquin J. Vresilovic E. Elliott D. Multiaxial validation of a finite element model of the intervertebral disc with multigenerational fibers to establish residual strain JOR spine 4 2021 e1145 10.1002/jsp2.1145
44 Sato K. Kikuchi S. Yonezawa T. In vivo intradiscal pressure measurement in healthy individuals and in patients with ongoing back problems Spine 24 1999 2468 2474 10.1097/00007632-199912010-00008 10626309
45 Nguyen A. Noninvasive quantification of human nucleus pulposus pressure with use of T1rho-weighted magnetic resonance imaging J. Bone Jt. Surg. Am. Vol. 90 2008 796 802 10.2106/jbjs.G.00667
46 Yang B. O'Connell G. Intervertebral disc swelling maintains strain homeostasis throughout the annulus fibrosus: a finite element analysis of healthy and degenerated discs Acta Biomater. 100 2019 61 74 10.1016/j.actbio.2019.09.035 31568880
47 Brent L.S. Novel human intervertebral disc strain template to quantify regional three-dimensional strains in a population and compare to internal strains predicted by a finite element model J. Orthop. Res. 34 2015 10.1002/jor.23137
48 Setton L. Zhu W. Weidenbaum M. Ratcliffe A. Mow V. Compressive properties of the cartilaginous end-plate of the baboon lumbar spine J. Orthop. Res. : official publication of the Orthopaedic Research Society 11 1993 228 239 10.1002/jor.1100110210
49 Shea M. Takeuchi T.Y. Wittenberg R.H. White A.A. 3rd Hayes W.C. A comparison of the effects of automated percutaneous diskectomy and conventional diskectomy on intradiscal pressure, disk geometry, and stiffness J. Spinal Disord. 7 1994 317 325 7949699
50 Johannessen W. Cloyd J. O'Connell G. Vresilovic E. Elliott D. Trans-endplate nucleotomy increases deformation and creep response in axial loading Ann. Biomed. Eng. 34 2006 687 696 10.1007/s10439-005-9070-8 16482409
51 Keller T.S. Spengler D.M. Hansson T.H. Mechanical behavior of the human lumbar spine. I. Creep analysis during static compressive loading J. Orthop. Res. 5 1987 467 478 10.1002/jor.1100050402 3681521
52 Vresilovic E. Johannessen W. Elliott D. Disc mechanics with trans-endplate partial nucleotomy are not fully restored following cyclic compressive loading and unloaded recovery J. Biomech. Eng. 128 2006 823 829 10.1115/1.2354210 17154681
53 D C.A. K I. S T. Direction-dependent resistance to flow in the endplate of the intervertebral disc: an ex vivo study J. Orthop. Res. 19 2002 10.1016/S0736-0266(01)00038-9
54 Yongren W. The region-dependent biomechanical and biochemical properties of bovine cartilaginous endplate J. Biomech. 48 2015 10.1016/j.jbiomech.2015.07.005
55 Antoniou J. Analysis of quantitative magnetic resonance imaging and biomechanical parameters on human discs with different grades of degeneration J. Magn. Reson. Imag. : JMRI 38 2013 1402 1414 10.1002/jmri.24120
56 Pieter-Paul A V. Albert J.v. d.V. Kaj S E. Jaap H.v.D. Theodoor H.S. The poro-elastic behaviour of the intervertebral disc: a new perspective on diurnal fluid flow J. Biomech. 49 2015 10.1016/j.jbiomech.2015.11.041
57 Reilly T. Tyrrell A. Troup J. Circadian variation in human stature Chronobiol. Int. 1 1984 121 126 10.3109/07420528409059129 6600017
58 Vergroesen P. Mechanics and biology in intervertebral disc degeneration: a vicious circle Osteoarthritis Cartilage 23 2015 1057 1070 10.1016/j.joca.2015.03.028 25827971
59 Eremina G. Smolin A. Xie J. Syrkashev V. Development of a computational model of the mechanical behavior of the L4-L5 lumbar spine: application to disc degeneration Materials 15 2022 10.3390/ma15196684
60 Jackson A.R. Yuan T.Y. Huang C.Y. Travascio F. Yong Gu W. Effect of compression and anisotropy on the diffusion of glucose in annulus fibrosus Spine 33 2008 1 7 10.1097/BRS.0b013e31815e4136 18165741
