# Biomechanical consequences of the intervertebral disc centre of rotation kinematics during lateral bending and axial rotation

Roman Allais, Antoine Capart, Anabela Da Silva, Olivier Boiron

► To cite this version:

Roman Allais, Antoine Capart, Anabela Da Silva, Olivier Boiron. Biomechanical consequences of the intervertebral disc centre of rotation kinematics during lateral bending and axial rotation. Scientific Reports, 2023, 13 (1), pp.3172. 10.1038/s41598-023-29551-7. hal-04041276

HAL Id: hal-04041276

https://hal.science/hal-04041276

Submitted on 20 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Biomechanical consequences of the intervertebral disc centre of rotation kinematics during lateral bending and axial rotation

Roman Allais<fnref n="1" />, Antoine Capart², Anabela Da Silva² & Olivier Boiron<fnref n="1" />

The location of the instantaneous centre of rotation (ICR) of a lumbar unit has a considerable clinical importance as a spinal health estimator. Consequently, many studies have been conducted to measure or estimate the ICR during rotations in the three anatomical planes; however the results reported are widely scattered. Even if some inter-subjects variability is to be expected, such inconsistencies are likely explained by the differences in methods and experiments. Therefore, in this paper we seek to model three behaviours of the ICR during lateral bending and axial rotation based on results published in the literature. In order to assess the metabolic and mechanical sensibility to the assumption made on the ICR kinematics, we used a previously validated three dimensional non-linear poroelastic model of a porcine intervertebral disc to simulate physiological lateral and axial rotations. The impact of the geometry was also briefly investigated by considering a 11° wedge angle. From our simulations, it appears that the hypothesis made on the ICR location does not significantly affect the critical nutrients concentrations but gives disparate predictions of the intradiscal pressure at the centre of the disc (variation up to 0.7 MPa) and of the displacement fields (variation up to 0.4 mm). On the contrary, the wedge angle does not influence the estimated intradiscal pressure but leads to minimal oxygen concentration decreased up to 33% and increased maximal lactate concentration up to 13%. While we can not settle on which definition of the ICR is more accurate, this work suggests that patient-specific modeling of the ICR is required and brings new insights that can be useful for the development of new tools or the design of surgical material such as total lumbar disc prostheses.

The position of the instantaneous centre of rotation (ICR) of a functional lumbar unit is a known metric of the spinal health and stability<fnref n="1" />,². In the past decades, extensive research has been conducted on the estimation of the ICR location during full range physiological motions but the results of these studies are surprisingly scattered. Ogston and colleagues, in their preliminary study³, described how to measure the centrode length (i.e. ICR trajectory) during flexion-extension in a clinical setting using radiographs. Meanwhile, Gertzbein et al. reported that intervertebral disc (IVD) degeneration at the L4/L5 level significantly impacts the length and pattern of the centrode during the same motion². More specifically, their results suggest that early stages of degeneration seem to increase the centrode length while higher stages decrease it. However, ICR location estimation through radiographs appears to be subject to considerable errors for rotation less than 5° in the sagittal plane, which raises the question of using an average centre of rotation instead⁴.

Regarding lateral bending (LB) of the lumbar spine, the average centre of rotation for right to left lateral bending appears to be located near the centre of the discs⁵. On the other hand, when studying right and left lateral bending separately two average centres of rotation are described in the study of Rolander⁶. Although the results are scattered, Rolander measured the average centre of rotation in the superior contralateral part of the disc (e.g. right side of the disc for left side bending), which was later corroborated by other researchers⁷. Nevertheless those studies only estimated an average centre of rotation. More recent works gave some insights on the motion of the ICR during lateral bending: for example, an upward shift towards the side of the bending for the L5/S1 joint

<footnote n="1">¹CNRS, Centrale Marseille, IRPHE, Aix Marseille Univ, 13013 Marseille, France. ²Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, 13013 Marseille, France. ℧email: roman.allais@centrale-marseille.fr</footnote>
was reported⁸. Finally, the numerical study of Schmidt and colleagues supported to some extent these findings as their finite elements model predicted the same trend for a L4/L5 joint⁹.

As for twisting (sometimes referred to as axial rotation (AR)), although not providing with numerical data, Cossette et al. measured the ICR in the posterior part of the disc with a trend of moving towards the side of the rotation with increasing rotation¹⁰. Quantitative results that refined the findings of Cossette and colleagues can be found in numerical studies: using a finite element model, a team of researchers reported the ICR trajectory as a function of increasing torque⁹. However, these results do not match the ones illustrated elsewhere where the ICR is displayed at the centre of the disc⁷.

While some of the variance in the ICR estimations found in the literature can be explained by inter-subjects variability, it should be borne in mind that the apparatus and experimental conditions used in the cited literature are diverse. One can indeed argue that results from in vivo, ex vivo and in silico studies can only be compared with care. The methods and key findings from the cited literature concerning the ICR locations during lateral bending and axial rotation are summed up in Table 1.

The estimation of the ICR of lumbar spine units recently received a renewed interest with the development of musculoskeletal models¹¹,¹². Musculoskeletal modeling can be used to help planning the surgery which will yield the best biomechanical outcomes; yet, the inverse dynamics resolution used in such models depends on the validity and modeling of the prescribed kinematics. A sensitivity study of the outcomes of a L4/L5 total disc replacement pointed out that a 5mm sagittal shift of the centre of rotation considerably impacts the predicted joint resultant, ligament and muscle forces during flexion-extension simulations¹¹. Senteler et al. refined the analysis with their parametric approach on the whole lumbar spine: they predicted that at 30° flexion the anterior shift of the ICR decreases the compression forces at a rate of 2.5%/mm¹². A proper understanding and an adequate estimation of the ICR location is therefore of prime interest for both clinical and modeling applications.

To get a right estimation of the ICR, the kinematics should be as close as possible to the physiological reality. There are two ways of modeling a rotation: through an applied moment⁹,¹³,¹⁴ or a prescribed displacement¹⁵. We deem the latter to better apprehend the physiology as it enforces the coupled rotations reported in the in vivo studies¹⁶⁻¹⁸ and it allows to specify the range of motion. On the other hand, with a moment applied, even though it would show the true kinematics of our disc, the range of motion can not be known beforehand and the ICR kinematics can not easily be constrained. As the goal of the study is to compare ICR behaviours reported in the literature, the kinematics of this specific disc are not of particular relevance.

Faced with the inconsistencies in the ICR locations and kinematics modeling, the aim of this paper is to provide with quantitative results on the ICR behaviour during lateral bending and axial rotation. Previous work has been published on the sensitivity of some biomechanical parameters to the position of an average centre of rotation during flexion/extension¹¹,¹². However, they made the convenient but non-physiological assumption of a fixed centre of rotation for their simulations. Others presented finite elements models, validated for LB and AR, but without considering their sensitivity to the ICR location¹⁹,²⁰. Finally, drawing meaningful comparisons between different numerical studies can be delicate as the anatomic completeness among the models varies. Also, we tried to tackle this issue by comparing the outputs of our previously validated porcine IVD model to its modified version which mimics the geometry of a L4/L5 human disc.

Hence, we seek here to assess the variation on both the mechanical and metabolic features of the IVD caused by the hypothesis made on the ICR kinematics. We considered 3 different locations of centre of rotation during lateral bending and axial rotation reported in the previous literature and we defined them in our finite element model of an IVD. For the sake of brevity, since less studies seem to have been conducted on the ICR kinematics during lateral bending and axial rotation, we chose to focus on those motions while flexion and extension could be the topic of a future study.

## Methods

**Finite elements model.** We used a previously published multiphysics three-dimensional nonlinear poroelastic model of a porcine lumbar intervertebral disc²¹,²². Succinctly, the geometry and porosity field in an unloaded state of the 4-month-old porcine disc were obtained using quantitative MRI using a 4.7-T horizontal superconducting magnet of a 47/30 Biospec Avance system (Bruker, Germany), equipped with a Bruker 120-mm BGA12SL (200 mT/m) gradient insert, using a proton RF coil (RF RES 200 1H Model 1PT5346; Bruker BioSpin MRI GmbH, Germany). The slice thickness was 0.75 mm in the transverse plane and 1 mm in the

<table><caption>**Table 1.** Summary of studies on the ICR location during lateral bending and axial rotation.</caption><thead><tr><td>Reference</td><td>Motion</td><td>Study type</td><td>Population</td><td>ICR location</td></tr></thead><tbody><tr><td>Dimnet et al.⁵</td><td>LB</td><td>In vivo</td><td>20 male subjects, 20–40 years old, levels T11–L5</td><td>IVD centre</td></tr><tr><td>Rolander⁶</td><td>LB</td><td>Ex vivo</td><td>38 cadavers (13 F/25 M), 4–76 years old, levels L1–L5</td><td>Superior, opposite side of the rotation</td></tr><tr><td>White and Panjabi⁷</td><td>LB</td><td>N/A</td><td>N/A</td><td>Opposite side of the rotation</td></tr><tr><td>Rousseau et al.⁸</td><td>LB</td><td>Ex vivo</td><td>12 cadavers (4 F/8 M), 50–64 years old, level L5/S1</td><td>Moving towards the side of the bending</td></tr><tr><td>Schmidt et al.⁹</td><td>LB</td><td>In silico</td><td>Finite-elements model of a L4/L5 level</td><td>IVD centre and moving towards the side of the bending</td></tr><tr><td>Cossette et al.¹⁰</td><td>AR</td><td>Ex vivo</td><td>12 cadavers, level L3/L4</td><td>Posterior side, moving towards the side of the rotation</td></tr><tr><td>White and Panjabi⁷</td><td>AR</td><td>N/A</td><td>N/A</td><td>IVD centre</td></tr><tr><td>Schmidt et al.⁹</td><td>AR</td><td>In silico</td><td>Finite-elements model of a L4/L5 level</td><td>IVD centre and moving posteriorly, towards the right side for an axial rotation to the left</td></tr></tbody></table>
sagittal plane²¹ and we used multi-slices multi-echo sequences²³. The main interest of this approach is to quantitatively and accurately retrieve the main internal biomechanical variables of the disc (i.e. the porosity field and the geometry). Then, the IVD was modeled by a biphasic fully saturated porous medium²⁴ with a swelling pressure²⁵ and a fiber-reinforced extracellular matrix to account for the anisotropy of the collagen network in the annulus fibrosus (AF)²⁶. The nucleus pulposus (NP), annulus fibrosus and cartilaginous endplates (CEP) were segmented using the porosity field as a proxy to identify the different anatomical compartments using the software Mimics 14.1 (Materialise NV, Leuven, Belgium). All the segmentations were performed by one operator to minimise the reconstruction errors (estimated to less than 3.5% on the disc volume). Each of these compartments had their own hyperelastic constitutive law that were fit to the measurements of relaxation tests made on the actual disc using an optimisation scheme and a transition area where all the physical parameters are smoothed²². This mechanical model is then coupled with a metabolic model²⁷ to predict the concentration fields of oxygen and lactate in response to a deformation; the glucose concentration can be derived from the energy metabolism equation assuming a fully glycolytic pathway. Pictures of the experimental study can be found in a previous publication²³ and the complete derivation of the model as well as the couplings between the metabolism and the deformation were already published²²,²⁸ but a summary of the main mechanical properties are recalled in Supplementary Table 1. However, human intervertebral discs display a wedge angle²⁹,³⁰, whereas our initial geometry does not. Therefore, we also performed simulations on a modified version of the original model to be closer to a human disc and to examine the sensitivity to the geometry. Using a CAD software, we modified the cartilaginous endplates to create a 11° wedge angle to mimic the geometry of a L4/L5 human disc³⁰. Removing materials would have further stray away our model from a human disc, that is why we chose to add materials instead. As a result, the posterior height was kept the same resulting in a 2.6mm increase in the anterior height and in a 7% increase of the disc volume. All the added volume was given the same properties as the original cartilaginous endplates; it should be noted that this assumption could lead to an oversized height of the endplate on the most anterior side but we do not think it would significantly alter our results as the cartilaginous end-plates are considerably stiffer than the nucleus and annulus and are therefore less prone to deformation. Both the original and modified geometry can be seen on Fig. 1 and more details on the meshes will be provided at the end of the section.

**Validation.** The validation of our initial model against experimental uniaxial relaxation tests was described in previous studies²¹,²². The validation is extended to lateral bending and axial rotation through a methodical comparison with previous models. Our results on the intradiscal pressure, displacement and concentration fields stay in the range of the values we observed during our experimental validation and were additionally compared with results from in vivo, in vitro and in silico studies.

**Initial state.** The imaging approach allows to recover the geometry and porosity field of the swollen disc at the physiological equilibrium. In order to retrieve the osmotic pressure and concentration fields corresponding to the imaged unloaded disc at this equilibrium, we ran a stationary boundary values problem. For the mechanical problem, no displacement is applied on the IVD and its lateral face is free to swell as a consequence of the osmotic pressure. For the nutrients diffusion problem, Dirichlet conditions are set on the three boundaries using the concentration values found in the plasma surrounding a human disc, as usually used in IVD finite-element models³¹⁻³³. The output of this simulation is then used as initial values for our dynamic simulations. All the simulations were performed using COMSOL Multiphysics 5.5 on a workstation with a Intel(R) Xeon(R) W-11955M CPU and 64GB of RAM.

**Motion description.** *Range of motion and displacement.* For both LB and AR, we chose to model the rotations by a prescribed displacement of the superior endplate while the inferior remains rigidly fixed. Thus, it makes possible to impose the coupled rotations and models a more physiological displacement. The displacement *u* underwent by the superior endplate is defined by Eq. (1). As for the amplitudes of rotations, we used

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two models of the intervertebral disc with their final mesh: initial geometry (left) and modified geometry (right). *x* axis points anteriorly, *y* points left and *z* points superiorly.</figcaption></figure>
<table><caption>**Table 2.** Prescribed rotations in the anatomical planes for LB and AR simulations. Source: Pearcy and Tibrewal¹⁶.</caption><thead><tr><th>Motion</th><th>Frontal (β)</th><th>Transverse (Ψ)</th><th>Sagittal (Θ)</th></tr></thead><tbody><tr><td>Lateral bending</td><td>2.5°</td><td>1°</td><td>-0.5°</td></tr><tr><td>Axial rotation</td><td>1.5°</td><td>1.5°</td><td>0</td></tr></tbody></table>

the ranges of motion measured by Pearcy and Tibrewal¹⁶. Briefly, they asked their subjects to perform full-range lateral and axial bending while they were taking biplanar radiographs to notably measure the three dimensional rotation between each lumbar level. The specific results for a L4/L5 disc we retained are reported in Table 2. Finally, in the vein of mimicking as close as possible a real motion, we smoothed the displacement over 4 s. At t = 0 s, the disc is in an unloaded state, at t = 2 s the maximal values presented in Table 2 are reached and at t = 4 s the disc is back in its initial state.

$$ \begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix} = \begin{pmatrix} \cos\theta \cos\Psi & \sin\beta\sin\theta\cos\Psi - \cos\beta\sin\Psi & \sin\beta\sin\Psi + \cos\beta\sin\theta\cos\Psi \\ \cos\theta\sin\Psi & \cos\beta\cos\Psi + \sin\beta\sin\theta\sin\Psi & \cos\beta\sin\theta\sin\Psi - \sin\beta\cos\Psi \\ -\sin\theta & \sin\beta\cos\theta & \cos\beta\cos\theta \end{pmatrix} \begin{pmatrix} x - c_x \\ y - c_y \\ z - c_z \end{pmatrix} - \begin{pmatrix} x - c_x \\ y - c_y \\ z - c_z \end{pmatrix} \quad (1) $$

With β, θ and Ψ as defined on Fig. 1.

**Lateral bending.** From there on, we define new variables: $bar_{x,y,z}$ and $c_{x,y,z}$ respectively stand for the spatial coordinates of the IVD barycentre and of the centre of rotation considered; $D_{sag}$, $D_{coro}$ and $H$ are the lengths of the IVD in respectively the antero-posterior, lateral and axial direction, passing through its barycentre. From the literature review displayed in Table 1, we kept three different definitions of the L4/L5 centre of rotation to implement them in our model:

LB-1: A centre of rotation set at the centre of the disc⁵:

$$ \begin{cases} c_x = bar_x \\ c_y = bar_y \\ c_z = bar_z \end{cases} \qquad (2) $$

LB-2: A centre of rotation fixed in the posterior, left side of the IVD for a right side bending⁶,¹⁵:

$$ \begin{cases} c_x = bar_x - \frac{D_{sag}}{4} \\ c_y = bar_y + \frac{D_{coro}}{5} \\ c_z = bar_z \end{cases} \qquad (3) $$

LB-3: A moving ICR, starting from the centre of the disc and moving towards the side of the bending⁹. However, Schmidt and colleagues reported the trajectory as a function of the increasing moment applied on the upper endplate of the L4 vertebra.

We therefore referred to the experimental data they used to calibrate their model¹³. These data, which are reported in Fig. 2A, give the relationship between the moment applied on the upper endplate of the L4 vertebra and the range of motion (RoM) of the L4/L5 segment. Finally, using this approach, we could model

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (A) Range of motion (RoM) over moment applied in the frontal plane for a L4/L5 human disc during lateral bending. Experimental values taken from Heuer et al.¹³. (B) Positions and trajectories of the 3 ICR definitions during LB. Green: hypothesis LB-1. Magenta: hypothesis LB-2. Red: hypothesis LB-3. (C) Positions and trajectories of the 3 ICR definitions during AR. Green: hypothesis AR-1. Magenta: hypothesis AR-2. Red: hypothesis LB-3.</figcaption></figure>
<table><caption>**Table 3.** Resolution of the considered meshes. D.o.F.: degrees of freedom.</caption><thead><tr><th colspan="4">Initial geometry</th><th colspan="4">Modified geometry</th></tr><tr><th>Mesh</th><th>Elements</th><th>D.o.F.</th><th>Min. element size</th><th>Mesh</th><th>Elements</th><th>D.o.F.</th><th>Min. element size</th></tr></thead><tbody><tr><td>Mesh 1</td><td>55,486</td><td>300,737</td><td>0.05 mm</td><td>Mesh 4</td><td>83,808</td><td>492,398</td><td>0.25 mm</td></tr><tr><td>Mesh 2</td><td>20,197</td><td>133,077</td><td>0.10 mm</td><td>Mesh 5</td><td>56,568</td><td>342,699</td><td>0.33 mm</td></tr><tr><td>Mesh 3</td><td>12,678</td><td>85,808</td><td>0.33 mm</td><td>Mesh 6</td><td>28,877</td><td>183,138</td><td>0.59 mm</td></tr></tbody></table>

the trajectory of the ICR as a function of the angle in the frontal plane. These three hypotheses are illustrated on Fig. 2B.

**Axial rotation.** For the axial rotation, we had the same approach as for LB: we kept three different definitions of the L4/L5 centre of rotation reported in previous works:

AR-1: A centre of rotation set at the disc barycentre⁷:

$$ \left\{ \begin{array}{l} c_x = bar_x \\ c_y = bar_y \\ c_z = bar_z \end{array} \right. \qquad (4) $$

AR-2: A centre of rotation in the superior ipsilateral side of the rotation¹⁵:

$$ \left\{ \begin{array}{l} c_x = bar_x - \frac{D_{sag}}{4} \\ c_y = bar_y + \frac{D_{coro}}{5} \\ c_z = H - h_{CEP} \end{array} \right. \qquad (5) $$

With $h_{CEP}$ the height of the superior cartilaginous endplate, which is 1.5 mm in our model.

AR-3: A moving ICR, starting from the centre of the disc and moving posteriorly in the right side of the disc⁹. Just as for LB, we used the calibration data¹³ to match the RoM to an equivalent moment in order to compute the trajectory in the transverse plane described by Schmidt and coworkers⁹. These three hypotheses are illustrated on Fig. 2C.

**Mesh convergence.** We performed a mesh convergence analysis to ensure mesh independence. By their off-centre locations, the centre of rotation kinematics as defined by hypotheses LB-2 and AR-2 are more likely to yield higher deformations and perturbations of both the mechanical and metabolic variables. We will therefore perform our mesh analysis on these specific kinematics for the initial and modified geometries.

Such a validation had already been done for our model for an uniaxial compression of the superior endplate²¹,²⁸. For the two geometries we built three meshes of tetrahedral elements whose resolutions are detailed in Table 3; Lagrange quadratic interpolation elements were used for the displacement field and linear interpolation elements were used for the pressure and concentrations. The choice of the optimal mesh for each geometry was based on a trade-off between computation time and relative errors on the displacement, porosity, concentrations and total pressure fields, the latter being defined as the sum of the osmotic and interstitial pressure fields.

**Anisotropic strain energy.** The consequences on the strain energy were also investigated. More specifically, we looked at the temporal variation of the anisotropic strain energy normalised by the disc volume and its proportion with respect to the total energy. Briefly, we define in the annulus two fibre families according to their mean orientation with respect to the transverse plane in order to model the anisotropy of the collagen fibres network. Then, two anisotropic strain energy densities, namely $W_{anis1}$ and $W_{anis2}$, are associated with each family and their sum represents the total anisotropic strain energy density; the details of the derivation was previously described²⁸ and the derivation is briefly recalled in Supplementary information 1. Finally, we computed the second Piola-Kirchhoff (PK) tensors associated with the anisotropic strain energy density ($PK_{aniso}$) and the one associated with the isotropic strain energy density ($PK_{iso}$). To investigate spatial variations as well, we looked at the ratio of the first eigenvalue of $PK_{aniso}$ over the first eigenvalue of $PK_{iso}$ at $t = 2$ s.

**Peclet number.** To draw a quantitative comparison between advection and diffusion, we computed the Peclet number associated with the oxygen and lactate. Half the posterior height of the disc was chosen as the characteristic diffusion length and the fluid velocity was derived using Darcy's law for a porous medium. Consequently, a global Peclet number was computed at each timestep using the spatial average of the fluid velocity.

## Results

**Mesh convergence.** To assess the mesh convergence, we took Mesh 1 and Mesh 4 as references and used two metrics to assess the numerical accuracy of the four other meshes. The relative error on the average value over the whole disc was used as a global metric and the relative error on the value at the centre of the disc as a
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Mesh analysis for LB (left) and AR (right) at the full range of motion. [.] is used to express the concentration.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (A–D) Empirical cumulative distribution functions of the displacement (A,C) and IDP fields (B,D). Solid line: initial geometry; dashed line: modified geometry. E: Critical nutrients concentration at *t* = 2 s for the initial (blue) and modified (orange) geometry. (F,G) Proportion of the anisotropic strain energy with respect to the total strain energy over time during LB (F) and AR (G). (H,I) Ratio of the first eigenvalue of *PK*<sub>aniso</sub> over the first eigenvalue of *PK*<sub>iso</sub> at *t* = 2 s during AR1 (H) and AR3 (I).</figcaption></figure>

local metric. The results are displayed on Fig. 3. The total pressure field is the variable which yields the highest average relative error, up to 4% computed when using Mesh 3. On the other hand, when looking at the local metric, the nutrients concentration shows the highest relative error, up to 10% for the oxygen on Mesh 6. Because of these large errors on local values, meshes 3 and 6 were considered not fine enough and discarded. Finally, considering their overall acceptable relative errors (< 4%) and computation time, Mesh 2 and Mesh 5 were kept.

**Initial state.** This conditioning step gives an average total pressure over the entire disc of 0.102 MPa with an increasing gradient as we move towards the centre of IVD. This is due to the higher concentration of fixed charges which consequently results in a higher osmotic pressure in the nucleus pulposus where the pressure reaches 0.203 MPa. The average oxygen concentration is 0.034 mol m⁻³ with a maximal value of 0.060 mol m⁻³ found at the boundaries of the AF and a minimal value of 0.014 mol m⁻³ inside the NP. The lactate concentration shows the opposite trend: values range from 0.9 mol m⁻³ at the external borders of the AF to 4.12 mol m⁻³ at the centre of the disc with an average of 2.566 mol m⁻³.

**IDP field.** Lateral bending yielded on average higher intradiscal pressure than axial rotation with no notable change between the two considered geometries as depicted by the empirical cumulative distributive functions (Fig. 4B,D). However the kinematics of the ICR drastically impacted the IDP field: the predicted values at the
<table><caption>**Table 4.** Peclet number associated with $O_2$ and lactate at $t = 1$ s for each configuration.</caption><thead><tr><th rowspan="3">Geometry</th><th colspan="6">Lateral bending</th><th colspan="6">Axial rotation</th></tr><tr><th colspan="3">Oxygen</th><th colspan="3">Lactate</th><th colspan="3">Oxygen</th><th colspan="3">Lactate</th></tr><tr><th>LB-1</th><th>LB-2</th><th>LB-3</th><th>LB-1</th><th>LB-2</th><th>LB-3</th><th>AR-1</th><th>AR-2</th><th>AR-3</th><th>AR-1</th><th>AR-2</th><th>AR-3</th></tr></thead><tbody><tr><td>Initial</td><td>743</td><td>903</td><td>817</td><td>1604</td><td>1950</td><td>1628</td><td>552</td><td>720</td><td>886</td><td>1193</td><td>1554</td><td>1913</td></tr><tr><td>Modified</td><td>694</td><td>860</td><td>787</td><td>1499</td><td>1858</td><td>1699</td><td>531</td><td>705</td><td>884</td><td>1147</td><td>1523</td><td>1909</td></tr></tbody></table>

centre of the disc range between -0.064 MPa (LB-3) and 0.251 MPa (LB-1) for the initial geometry and between -0.058 MPa (LB-3) and 0.249 MPa (LB-1) for the modified geometry during lateral bending. On a more global scale, the median IDP value was the lowest with a mobile ICR (LB3: 0.021 MPa, AR3: 0.02 MPa), followed by an ICR at the centre of the disc (LB1: 0.079 MPa, AR1: 0.078) and finally by an off-centre static ICR (LB2: 0.180 MPa, AR2: 0.142 MPa).

**Displacement field.** For a given ICR kinematics, the geometry did not significantly affect the displacement fields as depicted by their empirical cumulative distribution functions (Fig. 4A,C); the relative difference on the maximum of the displacement field between the two geometries was up to 2.8%. More specifically, During LB, the maximum of the displacement field considering a mobile centre of rotation was respectively 0.996 mm and 1.024 mm for respectively the initial and the modified geometry and their location was on the extreme left side of the disc. On the other hand, using a static centre of rotation (LB1, LB2) lead to a maximum of the displacement field in the right side with different magnitudes (LB1: 0.774 mm, LB2: 1.054 mm). Regarding AR, the maxima ranged from 0.600 mm (AR1) to 0.966 mm (AR3), passing through 0.833 mm (AR2). Its location was on the left side of the disc for the assumptions AR1 and AR3 while it was reported on the right side under hypothesis AR2.

**Nutrients concentration fields.** In line with the literature related to the IVD, we define the critical nutrients concentrations as the minimal oxygen ($[O_2]^*$) and maximal lactate concentrations ($[Lact]^*$). For all kinematics scenarios, the modified geometry with the wedge angle had a substantial influence on the critical concentrations when compared to the initial geometry. On average, the minimal oxygen concentration dropped by 32.3% while the maximal lactate concentration increased by 13.4% (Fig. 4E). The ICR kinematics only affected marginally the nutrients fields: considering a static centre of rotation at the disc barycentre (LB1, AR1) or a moving ICR (LB3, AR3) gave the same predictions ($[O_2]^* = 13.51$ mmol m$^{-3}$, $[Lact]^* = 4155$ mmol m$^{-3}$ during lateral bending and $[O_2]^* = 13.56$ mmol m$^{-3}$, $[Lact]^* = 4154$ mmol m$^{-3}$ during axial rotation). However, when considering the ICR in an off-centre location (LB2, AR2), $[O_2]^*$ increased on average by 5.2% and $[Lact]^*$ increased on average by 1%.

**Anisotropic strain energy.** For all the kinematics of interest, the proportion of anisotropic strain energy with respect to the total energy tended to increase with the range of motion (Fig. 4F,G). On average, axial rotation yielded higher proportion (up to 15.3%) than lateral bending. The geometry of the wedge angle had a marginal effect on the maximum of the anisotropic strain energy proportion (an average difference of 3.4% during LB and 5.1% during axial rotation). Spatial analysis of the ratio of the first eigenvalue of $PK_{aniso}$ over the first eigenvalue of $PK_{iso}$ during axial rotation revealed that the fibres in the postero-lateral compartments are the most solicited. However, the right side fibres were more solicited assuming hypothesis AR1 but the left side fibres were the most solicited assuming AR3. In both cases, the maximum of the ratio are comparable (AR1: 2.75, AR3: 2.53).

**Peclet number.** The maximum values (i.e. at $t = 1$ s) for each definition of the ICR and for the two geometries are gathered in Table 4. The reported values are of the same order of magnitude between the geometries and kinematics. A static centre of rotation at the disc barycentre systematically lead to lower maximal values.

## Discussion

The location of the intervertebral instantaneous centre of rotation is a controversial topic. Hence, the purpose of this study was to assess the sensitivity of the main IVD variables to the assumptions made on its centre of rotation kinematics and geometry during lateral bending and axial rotation. The major outcomes are that the kinematics of the centre of rotation considerably affects the pressure and displacement fields as well as the strain energy, while the wedge angle markedly impacts the nutrients' concentrations.

Porcine models are an opportune alternative to human models and are often used to investigate injury mechanisms, surgical approaches or implant tests; more details and references on this topic can be found elsewhere$^{34,35}$. Lumbar porcine discs are indeed similar to human discs in many ways. First, they have similar dimensions, albeit slightly smaller in the anteroposterior axis (average depth of 24.9 mm for the swines versus 34.6 mm for humans)$^{36,37}$. The anatomy of the lumbar unit, especially regarding the posterior elements, are also akin since no statistically significant differences on the spinous and transverse processes lengths and between-facets dimensions were reported$^{38}$. Then, both have analogous concentrations and gradients in terms of water content (67–83% for humans$^{39}$ versus 71–85% for pigs$^{40}$), proteoglycans content (humans: 5–11 µg of uronic acid by mg of wet tissue$^{39}$ versus 1 to 10 in a pig disc$^{41}$) and collagen content (humans: 30–48% of dry weight$^{42}$ versus 20.8–48.6% in pigs discs$^{40}$). Finally, from a chemical aspect, porcine and human discs are closely related as well: the ratio of
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (A) Comparison of the predicted nucleus pressure with previous studies. (1) Schmidt et al., 2006<sup>48</sup>, (2) Heuer et al., 2006<sup>49</sup>, (3) Park et al., 2013<sup>50</sup>, (4) Wilke et al., 2001<sup>51</sup>, (5): Wilke et al., 1996<sup>52</sup>. (B) Total reaction moment over time during LB considering a static (LB-1) and a mobile (LB-3) ICR and using the initial geometry. (C) Total reaction moment over time during AR considering a static (AR-1) and a mobile (AR-3) ICR and using the initial geometry.</figcaption></figure>

chondroitin sulfate to keratan sulfate was shown to be similar (1.75–4 for pigs<sup>41</sup> versus 1.79–2.5 for humans<sup>39</sup>). The main difference that has to be highlighted is the presence of notochordal cells in swine nuclei, even at the adult stage when such cells disappear in human discs during adolescence<sup>43</sup>. If this specificity has to be kept in mind while investigating disc degeneration or therapeutic treatments, it has, to the best of our knowledge, only little to no direct influence on the macroscopic, mechanical response of the disc.

However for the kinematic analysis we aimed to perform in this study, some questions may have arisen as pigs are quadrupeds when humans are bipeds. Previous studies compared human and porcine spines and reported that both spines had similar range of motion during lateral bending and axial rotation<sup>34</sup>: range of motion up to 6° were reported during lateral bending and up to 2° during axial rotation of a swine spine. These values match the upper values of the study of Pearce and Tibrewal<sup>16</sup> we used as reference for our simulations (AR: 3°, LB: 6°). Moreover, in standing situations, the human spine bears the weight of the upper body which is a notable difference with quadrupeds. Despite this distinction, previous work argues that the load exerted on quadrupeds spines are comparable in terms of magnitude and direction to the load put on humans spines, as a result of the muscles forces required for stabilisation<sup>44</sup>. Finally, the parameters (Supplementary Table 1) describing the mechanical behaviour of our porcine model fall in the range of values reported in previous models of human discs<sup>45-47</sup>. We therefore consider our model to be of valuable insights for the study of the human centre of rotation kinematics with the main advantage of having a realistic description, not only of the geometry like many other models, but more importantly of the internal biomechanical variables.

Our results can be validated through comparisons with previous works. Range of motion happens to be a preferential validation metric in the literature but as we used it as input through a prescribed displacement, we will instead use the intradiscal pressure. The pressure computed in the initial unloaded state is in good agreement with the in vivo data previously reported by Wilke and colleagues<sup>53</sup>, whose values in the NP range from 0.10 to 0.25 MPa for different kinds of lying positions. As for the specific validation during lateral bending and axial rotation, Fig. 5A displays a comparison of the mean IDP in the nucleus predicted by our simulations with measures or predictions found in the literature. The reported values show large variations that can be explained by the choice of the methodology. First, the age of the specimens used in the reference studies ranges from 26 to
52 which is most likely to span several disc degeneration grades which is known to alter the IDP. Additionally, the experimental conditions were diverse: the magnitude and type of loading differ between papers. Nonetheless, our mean values (0.22 MPa during AR and 0.31 MPa during LB) are close of the values predicted by the model of Schmidt and coworkers⁴⁸ for their combined moment without follower load conditions which can be compared to our simulations. The in vitro tests of Wilke and colleagues⁵² under pure moment solicitation also tend to validate our results. At last, the notable contrast with the in vivo measurements of Wilke⁵¹ might be explained by the muscles acting on the spine. In his *in vitro* experiment, he indeed showed the significant impact of the muscles on IDP measurements for a L4/L5 disc. Because of inter-subject variability, an accurate prediction of the IDP was proven to be a tough issue: a comparative study showed up to a five-times difference on the prediction between already validated models of the lumbar spines³⁴. As our estimations are in the range of the results presented in this above-mentioned comparative study, we consider our model valid.

The predicted maximum lactate concentration of 4.12 mol m⁻³ is below the value (5.45 mol m⁻³) computed by earlier models⁵⁵. However, both the average and maximal lactate concentration compares with the values predicted by a recent study using patient-specific models for a cohort of ten patients without disc-related disorder³³ and the computed values fall in the range 2–6 mol m⁻³ measured in human discs with different levels of degeneration⁵⁶. As for the oxygen levels, the critical oxygen concentration of 1.36 kPa is slightly above the measured range 0.3–1.1 kPa reported in canine nucleus⁵⁷ but considerably below the average value of 4 kPa predicted by Shalash et al.³³. As the values measured on human discs during spinal surgery show a significant variability (0.67–20 kPa⁵⁶), our predictions still seem acceptable. Finally, the critical oxygen and lactate concentration are found in the antero-lateral part of the disc, at the boundary between the inner annulus and the nucleus, which is consistent with previous findings³³,⁵⁵.

As modeled by the empirical work of Bibby and colleagues²⁷, the IVD metabolism is closely linked to the local concentration of oxygen and lactate through its connection to the medium pH. The concentration fields and the critical concentrations are therefore key features whose sensitivity needs to be investigated. Although the critical concentrations stay in the range of values reported in literature, we do report an average drop of 32.3% of the minimal oxygen concentration in the modified geometry with respect to the initial. As our metabolic model is only diffusive and does not account for any advection, this drop is most likely due to the higher diffusion distance from the cartilaginous endplates which are known to be the major nutrients supply pathways⁵⁷,⁵⁸. In parallel, using the modified geometry, the maximal lactate concentration predicted increases only by 13.4% on average; the relative difference with the drop in oxygen is explained by the different diffusibility of the nutrients. These results highlight the importance of the geometry to compute the nutrients concentration fields. A previous set of studies indeed showed a 3% increase (respectively 22% decrease) in the maximal lactate (respectively minimal oxygen) concentration between a 3D axisymmetric and a 3D realistic model of a lumbar disc⁵⁵,⁵⁹. The significant impact of the disc size was also corroborated and quantified by a recent study³³ which performed a sensitivity analysis to a change on the diffusivity, reaction rates or disc size.

Lateral bending and axial rotation are, by definition, non-uniaxial motions. As a consequence, large pressure gradients arose in our simulations which can impact the fluid flows. Even though we chose a motion duration twice as long as in a previous study¹⁵, we have the same order of magnitude for our Peclet numbers (10⁰ to 10³) which confirms that the advection should be taken into account when modelling physiological rotations of the spine. However, a spatial analysis using local values of fluid velocity and diffusibility showed significant variations up to 60 times lower than the values reported in Table 4. It appeared that the high global value of the Peclet number was the consequence of the high local fluid velocity at the lateral boundaries of the disc. We therefore acknowledge that our numerical values should be interpreted with care, but as a complete assessment of the impacts of the advective term showed that for LB and AR there are no significant differences on the critical concentrations and only moderate differences on the concentration fields near the borders⁶⁰, we consider our predictions reliable. We also investigated the potential sensitivity of the porosity and osmotic pressure fields as they are the variables driving the water flows. Our dynamic simulations did not show differences between the geometries or the assumptions on the ICR. Although not representative of activities of daily living, we also ran stationary simulations for LB-1, LB-2, AR-1 and AR-2 hypotheses, which showed respectively a 12.5% and 25% increase of the maximal osmotic pressure with respect to the initial state during respectively axial rotation and lateral bending. These differences between stationary and dynamic simulations are likely explained by the diffusion being too slow with respect to the motion duration.

While the predicted IDP and displacement fields from our simulations appear to be fairly insensitive to the wedge angle (Fig. 4A–D), the ICR kinematics strongly altered these fields. For both LB and AR, a mobile ICR yields substantially lower predictions of the IDP. A possible explanation may lie in an underlying pressure optimisation mechanism. This idea was first raised using musculoskeletal models: in a sensitivity study of the lumbar spine kinetics to the position of the ICR during flexion, it was suggested that the ICR is moving to minimize joint reaction forces¹², which would subsequently minimize the IDP. However the authors only investigated this idea for motion restricted to the sagittal plane and they acknowledged that further work would be needed to conclude on this interpretation. Our results tend to be in agreement with such a hypothesis: with a mobile ICR, the total moment undergone by the disc over time is up to 19% smaller than with a static ICR set at the centre of the disc (Fig. 5B,C). The evolution of the IDP at the centre of the disc also hints in this direction as at the full range of motion the total pressure is on average 0.326 MPa lower with a mobile ICR than with a static ICR at the disc barycentre. On the other hand, a mobile ICR also increases the maximal values of the displacement field up to 28.7% during LB and 61% during AR.

Furthermore, in our simulations a mobile ICR increased the proportion of anisotropic strain energy, especially during axial rotation where it represents up to 15.3% of the total strain energy, with the proportion increasing with the range of motion. In addition, the analysis of the Piola-Kirchhoff stress tensors associated with the isotropic and anisotropic part of the extracellular matrix (ECM) revealed spatial differences in terms of fibres
solicitation. First, it showed that in the areas of higher deformation, most of the stress is put on the fibres (Fig. 4H,I). The collagen fibres are therefore considerably solicited in activities of daily life requiring end-range rotations of the spine. Finally, the fibres in the postero-lateral compartments of the IVD sustain the highest stress which is consistent with clinical observation reporting these areas as common places of injuries such as annular tears. On the other hand, the results displayed on Fig. 4F,G point out that for a smaller range of motion, equivalent to the first 0.5 s of our simulations, the anisotropic strain energy proportion stays nearly constant for every conditions. It implies that the isotropic part of the ECM is the first to stretch and once a certain threshold is reached, the anisotropic component of the strain energy quickly rises.

In short, we tested how the different definitions of the ICR kinematics of a L4/L5 lumbar segment reported in the literature impact the main biomechanical variables of the IVD during lateral bending and axial rotation. Previous studies proposed sensitivity analyses but only for sagittal motions and/or using arbitrary, non-physiological shifts of the ICR<sup>11,12</sup>. Additionally, previous works predominantly focused on how the rotations alter the mechanical variables, such as intradiscal pressure, facet forces, fibre strains or joint reaction forces<sup>8,12,19,48,61</sup> without looking at the nutrients. On the other hand, papers related to nutrients distribution mostly modeled uni-axial compression<sup>32,59,62</sup>. To our knowledge, we are the first to quantify the differences consequent to the assumption made on the ICR kinematics using finite-element simulations accounting for both the nutrients distributions and the mechanical variables.

Our results find clinical applications, notably for the design and implantation of lumbar discs prostheses. Such prostheses are of particular interest as they are a solution to avoid spinal fusion which is known to significantly decrease the range of motion and have noxious effects on adjacent levels. However, a physiological, moving ICR is required to avoid an overload of the other elements of the spine, such as the facets<sup>63</sup>. Hence our simulations can provide some insights to predict the locations with the highest displacement which can be of help for the conception of new devices.

Even though we wanted to be as close as possible to the physiological reality, some limitations have to be acknowledged and will need to be addressed in future work. Since the L4/L5 lumbar unit is the most prone to degeneration, we chose to focus our study on this level. It can be hypothesized that the behaviour of the other lumbar units would qualitatively be close but our results can not be directly transposed without further validation. Additionally, we used a non-degenerate disc which restricts our findings to healthy discs only. Even though porcine model are appropriate for the study of the spine, future work would greatly benefit from an experimental validation on L4/L5 healthy human discs. For the hypotheses of a moving ICR (LB-3 and AR-3) we used the data reported by Schmidt et al.<sup>9</sup>. However, they computed those trajectories using pure applied moments in the anatomical planes while we made the choice to model the coupled rotations described in the in vivo studies. Nonetheless, the ligaments and posterior elements of their model generate reaction moments and consequently create coupled motion, qualitatively similar to our simulations. The trajectories of the ICR are also constrained in an anatomical plane (frontal for LB and transverse for AR) which is undoubtedly an approximation that does not hold in reality. Further work could focus on more accurate, three-dimensional descriptions of the centre of rotation or alternative ways to model the motion, such as finite helical axis<sup>64</sup>. Finally, it could also be argued that our model lacks of anatomical completeness but as we modeled the same motion for all conditions, it should not bias results. Musculoskeletal simulations entailing more kinematically relevant anatomy (e.g. with muscles and ligaments) could also be of interest and bring more insights on the ICR location.

## Conclusion

We showed manifest evidence that the instantaneous centre of rotation (ICR) kinematics considerably tailor the IVD mechanical response during lateral bending and axial rotation, which emphasizes the importance of patient-specific modeling. Most notably, our work showed that a mobile ICR decreases the intradiscal pressure while increasing the maximal displacement and the stress put on the collagen fibres. On the contrary, it has marginal effects on the nutrients critical concentration; hence an average centre of rotation is accurate enough for nutrients-related simulations but not for the prediction of the displacement and pressure fields. Nonetheless, our results only relate to the L4/L5 level and could benefit from further experimental validation on human healthy discs. Consequently, this work brings new insights on the IVD biomechanics which can notably be of help for the creation of new models or the modeling and design of surgical materials.

## Data availability

Upon reasonable request to roman.allais@centrale-marseille.fr.

Received: 9 September 2022; Accepted: 6 February 2023

Published online: 23 February 2023

## References

1. Seligman, J., Gertzbein, S., Tile, M. & Kapasouri, A. Computer analysis of spinal segment motion in degenerative disc disease with and without axial loading. *Spine* **9**, 566–573. https://doi.org/10.1097/00007632-198409000-00006 (1984).

2. Gertzbein, S. et al. Centrode patterns and segmental instability in degenerative disc disease. *Spine* **10**(3), 257–61. https://doi.org/10.1097/00007632-198504000-00014 (1985).

3. Ogston, N. et al. Centrode patterns in the lumbar spine, baseline studies in normal subjects. *Spine* **11**, 591–595. https://doi.org/10.1097/00007632-198607000-00010 (1986).

4. Pearcy, M. & Bogduk, N. Instantaneous axes of rotation of the lumbar intervertebral joints. *Spine* **13**, 1033–1041. https://doi.org/10.1097/00007632-198809000-00011 (1988).

5. Dimnet, J., Fischer, L., Gonon, G. & Carret, J. Radiographic studies of lateral flexion in the lumbar spine. *J. Biomech.* **11**, 143–150. https://doi.org/10.1016/0021-9290(78)90006-4 (1978).
6. Rolander, S. D. Motion of the lumbar spine with special reference to the stabilizing effect of posterior fusion: an experimental study on autopsy specimens. *Acta Orthopaed. Scand.* **37**, 1-144. https://doi.org/10.3109/ort.1966.37.suppl.90.1 (1966).

7. White, A. & Panjabi, M. Kinematics of the spine. In *Clinical biomechanics of the spine. 2nd ed.*, 109-112 (Philadelphia: JB Lippincott Company, 1990).

8. Rousseau, M. A., Bradford, D. S., Hadi, T. M., Pedersen, K. L. & Lotz, J. C. The instant axis of rotation influences facet forces at l5/ s1 during flexion/extension and lateral bending. *Eur. Spine J.* **15**, 299-307. https://doi.org/10.1007/s00586-005-0935-1 (2006).

9. Schmidt, H., Heuer, F., Claes, L. & Wilke, H.-J. The relation between the instantaneous center of rotation and facet joint forces—A finite element analysis. *Clin. Biomech.* **23**, 270-278. https://doi.org/10.1016/j.clinbiomech.2007.10.001 (2008).

10. Cossette, J., Farfan, H., Robertson, G. & Wells, R. The instantaneous center of rotation of the third lumbar intervertebral joint. *J. Biomech.* **4**, 149-153. https://doi.org/10.1016/0021-9290(71)90025-X (1971).

11. Han, K.-S., Kim, K., Park, W. M., Lim, D. S. & Kim, Y. H. Effect of centers of rotation on spinal loads and muscle forces in total disk replacement of lumbar spine. *Proc. Inst. Mech. Eng. Part H J. Eng. Med.* **227**, 543-550. https://doi.org/10.1177/0954411912 474742 (2013) (PMID: 23637264).

12. Senteler, M., Aiyangar, A., Weisse, B., Farshad, M. & Snedeker, J. G. Sensitivity of intervertebral joint forces to center of rotation location and trends along its migration path. *J. Biomech.* **70**, 140-148. https://doi.org/10.1016/j.jbiomech.2017.10.027 (2018).

13. Heuer, F., Schmidt, H., Klezl, Z., Claes, L. & Wilke, H.-J. Stepwise reduction of functional spinal structures increase range of motion and change lordosis angle. *J. Biomech.* **40**, 271-280. https://doi.org/10.1016/j.jbiomech.2006.01.007 (2007).

14. Rohlmann, A., Zander, T., Rao, M. & Bergmann, G. Realistic loading conditions for upper body bending. *J. Biomech.* **42**, 884-890. https://doi.org/10.1016/j.jbiomech.2009.01.017 (2009).

15. Etienne, M., Deplano, V., Boiron, O. & Tropiano, P. 3D dynamic numerical simulations of intervertebral disc: Bending and twisting. *Comput. Methods Biomech. Biomed. Eng.* **16**, 236-238. https://doi.org/10.1080/10255842.2013.815948 (2013).

16. Pearcy, M. & Tibrewal, S. Axial rotation and lateral bending in the normal lumbar spine measured by three-dimensional radiography. *Spine* **9**, 582-587. https://doi.org/10.1097/00007632-198409000-00008 (1984).

17. Li, G. et al. Segmental in vivo vertebral motion during functional human lumbar spine activities. *Eur. Spine J.* **18**, 1013-1021. https://doi.org/10.1007/s00586-009-0936-6 (2009).

18. Berry, D. B. et al. Lumbar spine angles and intervertebral disc characteristics with end-range positions in three planes of motion in healthy people using upright mri. *J. Biomech.* **89**, 95-104. https://doi.org/10.1016/j.jbiomech.2019.04.020 (2019).

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

20. Xu, M., Yang, J., Lieberman, I. H. & Haddas, R. Lumbar spine finite element model for healthy subjects: Development and validation. *Comput. Methods Biomech. Biomed. Eng.* **20**, 1-15. https://doi.org/10.1080/10255842.2016.1193596 (2017).

21. Chetoui, M., Boiron, O., Dogui, A. & Deplano, V. Prediction of intervertebral disc mechanical response to axial load using isotropic and fiber reinforced Fe models. *Comput. Methods Biomech. Biomed. Eng.* **20**, S39-S40. https://doi.org/10.1080/10255842.2017. 1382850 (2017).

22. Chetoui, M. A., Boiron, O., Ghiss, M., Dogui, A. & Deplano, V. Assessment of intervertebral disc degeneration-related properties using finite element models based on rho-h -weighted MRI data. *Biomech. Model. Mechanobiol.* **18**, 17-28. https://doi.org/10.1007/ s10237-018-1064-1 (2019).

23. Ghiss, M., Giannesini, B., Tropiano, P., Tourki, Z. & Boiron, O. Quantitative MRI water content mapping of porcine intervertebral disc during uniaxial compression. *Comput. Methods Biomech. Biomed. Eng.* **19**, 1079-1088. https://doi.org/10.1080/10255842. 2015.1101072 (2016).

24. Mow, V. C., Kuei, S., Lai, W. M. & Armstrong, C. G. Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiments. *J. Biomech. Eng.* **102**, 73-84. https://doi.org/10.1115/1.3138202 (1980).

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

26. Holzapfel, G. A., Sommer, G., Gasser, C. T. & Regnitig, P. Determination of layer-specific mechanical properties of human coronary arteries with nonatherosclerotic intimal thickening and related constitutive modeling. *Am. J. Physiol.-Heart Circ. Physiol.* **289**, H2048-H2058. https://doi.org/10.1152/ajpheart.00934.2004 (2005).

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

28. Chetoui, M.-A. *Modélisation du Comportement Biomécanique du Disque Intervertebral*. Theses, Ecole Centrale Marseille, École Nationale d'Ingénieurs de Monastir (Tunisie) (2017).

29. Kimura, S., Steinbach, G. C., Watenpaugh, D. E. & Hargens, A. R. Lumbar spine disc height and curvature responses to an axial load generated by a compression device compatible with magnetic resonance imaging. *Spine* **26**, 2596-2600. https://doi.org/10. 1097/00007632-200112010-00014 (2001).

30. De Carvalho, D., Grondin, D. & Callaghan, J. The impact of office chair features on lumbar lordosis, intervertebral joint and sacral tilt angles: A radiographic assessment. *Ergonomics* **60**, 1393-1404. https://doi.org/10.1080/00140139.2016.1265670 (2017).

31. Soukane, D. M., Shirazi-Adl, A. & Urban, J. Analysis of nonlinear coupled diffusion of oxygen and lactic acid in intervertebral discs. *J. Biomech. Eng.* https://doi.org/10.1115/1.2073674 (2005).

32. Jackson, A. R., Huang, C.-Y. & Gu, W. Y. Effect of endplate calcification and mechanical deformation on the distribution of glucose in intervertebral disc: A 3D finite element study. *Comput. Methods Biomech. Biomed. Eng.* **14**, 195-204. https://doi.org/10.1080/ 10255842.2010.535815 (2011).

33. Shalash, W., Ahrens, S. R., Bardonova, L. A., Byvaltsev, V. A. & Giers, M. B. Patient-specific apparent diffusion maps used to model nutrient availability in degenerated intervertebral discs. *JOR Spine* **4**, e1179. https://doi.org/10.1002/jsp2.1179 (2021).

34. Wilke, H.-J., Geppert, J. & Kienle, A. Biomechanical in vitro evaluation of the complete porcine spine in comparison with data of the human spine. *Eur. Spine J.* **20**, 1859-1868. https://doi.org/10.1007/s00586-011-1822-6 (2011).

35. Patel, S. A., Kepler, C. K., Schaer, T. P. & Anderson, D. G. Large animal models of disc degeneration. *Intervertebral Disc*.https:// doi.org/10.1007/978-3-7091-1535-0_18 (2014).

36. Chanchairujira, K. et al. Intervertebral disk calcification of the spine in an elderly population: Radiographic prevalence, location, and distribution and correlation with spinal degeneration. *Radiology* **230**, 499-503. https://doi.org/10.1148/radiol.2302011842 (2004).

37. Wang, J.-L., Tsai, Y.-C. & Wang, Y.-H. The leakage pathway and effect of needle gauge on degree of disc injury post anular puncture: A comparative study using aged human and adolescent porcine discs. *Spine* **32**, 1809-1815. https://doi.org/10.1097/brs.0b013e3181 1ec282 (2007).

38. McLain, R. F., Yerby, S. A. & Moseley, T. A. Comparative morphometry of l4 vertebrae: Comparison of large animal models for the human lumbar spine. *Spine* **27**, E200-E206. https://doi.org/10.1097/00007632-200204150-00005 (2002).

39. Lyons, G., Eisenstein, S. M. & Sweet, M. Biochemical changes in intervertebral disc degeneration. *Biochim. Biophys. Acta (BBA) Gen. Subj.* **673**, 443-453. https://doi.org/10.1016/0304-4165(81)90476-1 (1981).

40. Holm, S., Holm, A. K., Ekström, L., Karladani, A. & Hansson, T. Experimental disc degeneration due to endplate injury. *Clin. Spine Surg.* **17**, 64-71. https://doi.org/10.1097/00024720-200402000-00012 (2004).

41. Kääpä, E. et al. Proteoglycan chemistry in experimentally injured porcine intervertebral disk. *J. Spinal Disord.* **7**, 296-306 (1994).
42. Adams, P., Eyre, D. & Muir, H. Biochemical aspects of development and ageing of human lumbar intervertebral discs. *Rheumatology* **16**, 22–29. https://doi.org/10.1093/rheumatology/16.1.22 (1977).

43. Omlor, G. W. et al. A new porcine in vivo animal model of disc degeneration: Response of anulus fibrosus cells, chondrocyte-like nucleus pulposus cells, and notochordal nucleus pulposus cells to partial nucleotomy. *Spine* **34**, 2730–2739. https://doi.org/10.1097/BRS.0b013e3181b723c9 (2009).

44. Smit, T. H. The use of a quadruped as an in vivo model for the study of the spine-biomechanical considerations. *Eur. Spine J.* **11**, 137–144. https://doi.org/10.1007/s005860100346 (2002).

45. Galbusera, F. et al. Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs. *J. Mech. Behav. Biomed. Mater.* **4**, 1234–1241. https://doi.org/10.1016/j.jmbbm.2011.04.008 (2011).

46. Castro, A., Wilson, W., Huyghe, J., Ito, K. & Alves, J. Intervertebral disc creep behavior assessment through an open source finite element solver. *J. Biomech.* **47**, 297–301. https://doi.org/10.1016/j.jbiomech.2013.10.014 (2014).

47. Cheung, J.T.-M., Zhang, M. & Chow, D.H.-K. Biomechanical responses of the intervertebral joints to static and vibrational loading: A finite element study. *Clin. Biomech.* **18**, 790–799. https://doi.org/10.1016/s0268-0033(03)00142-6 (2003).

48. Schmidt, H. et al. Intradiscal pressure, shear strain, and fiber strain in the intervertebral disc under combined loading. *Spine* **32**, 748–755. https://doi.org/10.1097/01.brs.0000259059.90430.c2 (2007).

49. Heuer, F., Schmidt, H., Claes, L. & Wilke, H.-J. Stepwise reduction of functional spinal structures increase vertebral translation and intradiscal pressure. *J. Biomech.* **40**, 795–803. https://doi.org/10.1016/j.jbiomech.2006.03.016 (2007).

50. Park, W. M., Kim, K. & Kim, Y. H. Effects of degenerated intervertebral discs on intersegmental rotations, intradiscal pressures, and facet joint forces of the whole lumbar spine. *Comput. Biol. Med.* **43**, 1234–1240. https://doi.org/10.1016/j.compbiomed.2013.06.011 (2013).

51. Wilke, H.-J., Neef, P., Hinz, B., Seidel, H. & Claes, L. Intradiscal pressure together with anthropometric data—A data set for the validation of models. *Clin. Biomech.* **16**, S111–S126. https://doi.org/10.1016/S0268-0033(00)00103-0 (2001).

52. Wilke, H.-J., Wolf, S., Claes, L. E., Arand, M. & Wiesend, A. Influence of varying muscle forces on lumbar intradiscal pressure: An in vitro study. *J. Biomech.* **29**, 549–555. https://doi.org/10.1016/0021-9290(95)00037-2 (1996).

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

54. Dreischarf, M. et al. Comparison of eight published static finite element models of the intact lumbar spine: Predictive power of models improves when combined together. *J. Biomech.* **47**, 1757–1766. https://doi.org/10.1016/j.jbiomech.2014.04.002 (2014).

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

56. Bartels, E. M., Fairbank, J. C., Winlove, C. P. & Urban, J. P. Oxygen and lactate concentrations measured in vivo in the intervertebral discs of patients with scoliosis and back pain. *Spine* **23**, 1–7. https://doi.org/10.1097/00007632-199801010-00001 (1998).

57. Holm, S., Maroudas, A., Urban, J., Selsam, G. & Nachemson, A. Nutrition of the intervertebral disc: Solute transport and metabo- lism. *Connect. Tissue Res.* **8**, 101–119. https://doi.org/10.3109/03008208109152130 (1981).

58. Sélard, É., Shirazi-Adl, A. & Urban, J. P. Finite element study of nutrient diffusion in the human intervertebral disc. *Spine* **28**, 1945–1953. https://doi.org/10.1097/01.BRS.0000087210.93541.23 (2003).

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

60. Etienne, M. *Modelisation du Disque Intervertebral. Theses* (Aix-Marseille Université, 2013).

61. Kuo, C.-S. et al. Biomechanical analysis of the lumbar spine on facet joint force and intradiscal pressure—A finite element study. *BMC Musculoskelet. Disord.* **11**, 1–13. https://doi.org/10.1186/1471-2474-11-151 (2010).

62. Huang, C.-Y. & Gu, W. Y. Effects of mechanical compression on metabolism and distribution of oxygen and lactate in intervertebral disc. *J. Biomech.* **41**, 1184–1196. https://doi.org/10.1016/j.jbiomech.2008.02.002 (2008).

63. Choi, J., Shin, D.-A. & Kim, S. Biomechanical effects of the geometry of ball-and-socket artificial disc on lumbar spine: A finite element study. *Spine* **42**, E332–E339. https://doi.org/10.1097/BRS.0000000000001789 (2017).

64. Kettler, A. et al. Finite helical axes of motion are a useful tool to describe the three-dimensional in vitro kinematics of the intact, injured and stabilised spine. *Eur. Spine J.* **13**, 553–559. https://doi.org/10.1007/s00586-004-0710-8 (2004).

## Acknowledgements

This work is funded by the Emplois Jeunes Doctorants 2021 grant 2021_04412 of the Région Provence Alpes Côte d'Azur.

## Author contributions

Conceptualization: R.A., O.B.; methodology: R.A.; numerical modeling: R.A.; analysis and validation: R.A., A.C., A.D.S. and O.B.; writing—draft preparation: R.A.; writing—reviews and editing: R.A., A.C., A.D.S. and O.B.; supervision: O.B. and A.D.S. All authors have reviewed and agreed to the published version of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-29551-7.

**Correspondence** and requests for materials should be addressed to R.A.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023