• An automatic mesh morphing procedure was established to transform a structured IVD FE mesh, previously calibrated (Ruiz et al., 2013), into IVD PP models based on segmented medical images.

A cohort of 169 geometrical models of the lumbar spine IVD (Pfirrman scale range from 1 to 4) was acquired through T2-weighted magnetic resonance imaging (MRI) from the former project MySpine (FP7-269909) to generate the morphed FE meshes. These models were developed following the innovative approach introduced by Castro-Mateos et al. (2014), which enables accurate 3D segmentation of IVD annulus and nucleus, with increased resolution as axial and transversal MRI slices (isotropic voxels of 0.68 × 0.68 × 0.68 mm³) are combined. This method is effective for all degrees of IDD, including discs with protrusions or herniations. The segmentation algorithm, employing a feature selector, iteratively deforms an initial shape, which is projected into a statistical shape model space at first and then into a B-Spline space to improve accuracy. Expert clinicians validated the initial manual segmentation and the generated 3D morphologies.

The point cloud (PC) of the AF and NP external surfaces could be extracted from the MRIs, but the segmented models did not incorporate the CEP (Figure 1B). To overcome this limitation, our morphing process estimates the shape of the CEP (presented in Figure 1C). This algorithm tailors an IVD FE generic mesh (previously calibrated through a comprehensive mesh convergence analysis (Ruiz et al., 2013), and originally developed by Noailly et al. (2010), presented in Figure 1A) to the PC of the PP models while preserving the relative dimensions of the elements and the mesh structure at material discontinuities. The entire mesh of the IVD contains 83,481 nodes, and the disc tissues—namely AF, NP, and CEP—were discretized with 19,392 second-order hexahedral elements (20 nodes).

The collagen fibers of the AF were modeled with 7,680 second-order quadrilateral rebar elements (8 nodes) embedded in the hexahedral elements of the annulus to model the fiber-reinforcement of the tissue. These rebars follow the concentric mesh structure of the annulus hexahedral elements to represent the functional mechanical behavior of the annulus structure within the IVD (Noailly et al., 2010).

The primary purpose of the pipeline presented in Figure 1 is to create a FE PP mesh file (.inp) compatible with the ABAQUS mechanical solver, with a mesh structure able to preserve simulation convergence and minimize poro-mechanical instabilities under physiological loading conditions (Ruiz et al., 2013). Furthermore, our algorithm guarantees the non-convexity of the elements, contributing to the simulations’ accuracy.

The simulations adopt the same boundary conditions, material properties, and physiological loads used by Ruiz et al. (2018) for healthy discs, based on in vivo intradiscal pressure measurements by Wilke et al. (1999). Each PP FE model has the same number of elements and nodes and the same connectivity. The sole difference among the PP models lies in the coordinates of the nodes within the FE mesh.

To validate our model’s capacity to represent the mechanical behavior of a normal healthy disc, we compared our simulations with the in vitro vertical creep displacement tests by Malandrino et al. (2015b). We employ the morphing algorithm to obtain the morphology of the IVD used in these experiments and consider its effects on the validation simulation.

Thanks to the consistent topology across the 169 PP morphed models, extracting mechanical responses from targeted zones, such as the NP’s center, the anterior and posterior TZs (regions prone to early degeneration, observed by (Smith et al., 2011) through clinical images), and the CEP—was straightforward and systematic for all models. Accordingly, local mechanical predictions could consistently correlate with various morphological factors, such as heights, areas, and volumes. The relevance of each factor was then ranked by leveraging SHAP (SHapley Additive exPlanations) values (Ning et al., 2022).

The Bayesian Coherent Point Drift ++ (BCPD++) (Hirose, 2021a; Hirose, 2021b) was employed to adapt the structured generic IVD FE mesh (Ruiz et al., 2013) to the PP point cloud (AF and NP). This Bayesian framework allows the inclusion of prior knowledge about the distribution of potential transformations, such as degrees of translation, rotation, and scaling among the points. Overall, the BCPD++ algorithm was chosen to facilitate the automation of transforming FE meshes.

• Algorithm convergence is guaranteed by variational Bayesian inference while introducing motion coherence using a prior distribution of displacement vectors.

To adapt the mesh to the PP models and estimate the CEP shape, the BCPD++ algorithm was regulated through our in-house Python script named morph.py (Figure 1C). Additionally, the script substitutes the original hexahedral elements of the BEP (Noailly et al., 2007; Malandrino et al., 2011; Ruiz et al., 2013) from the generic mesh with 3,552 eight-node shell elements in the top and bottom regions of the disc, consequently reducing the computation time of both the morphing and the FE simulations.

In brief, the algorithm has two main parameters that were controlled hereby: λ, which controls the expected length of the displacement vectors, and β, which controls the directional correlation among the displacement vectors. λ and β facilitate the balance of rigid and non-rigid registrations by using large or small values of each parameter, which enables the control of the distortion of the deformation field.

To accurately estimate the CEP, it is crucial to set a specific target thickness range to ensure that modeled healthy discs do not fall below the expected thickness values. Moon et al. (2013) found, via MRI observations, that the central thickness of a healthy disc—considered the minimum thickness across the entire CEP—measured in the sagittal plane averages 0.54 ± 0.12 mm across all lumbar levels. This finding is consistent with our generic model, which specifies a CEP thickness at the center of the disc of 0.545 mm, initially based on histological measurements (Noailly et al., 2007). Consequently, ensuring that the average CEP thickness of the discs generated by the morphing process falls within this range becomes a key objective, particularly for discs not in advanced degeneration stages. Additionally, focusing on the central thickness allows for the creation of CEPs with variable thickness across the tissue, adapting its shape to match the morphology of the AF and NP, thereby introducing a more realistic and nuanced approach to modeling healthy and degenerated discs.

The process of creating the PP models is divided into seven key steps, presented in Figure 1C.

• Step 1: Initial Alignment—The IVD generic FE mesh and the PP point cloud are aligned using a rigid registration process. This step employs a theoretically infinitely large value for both λ and β to achieve a perfect initial alignment without any deformation of the models. Here, λ = 1E9 and β = 1E9 are suggested by Hirose (2021b) to represent the infinity.

Once all the IVD morphed models were prepared, the mesh quality was contrasted with the generic FE mesh of Ruiz et al. (2013). This allowed us to assess whether the BCPD++ algorithm could maintain the proportionality of the nodes’ relative distances without inducing excessive deformation.

Each BCPD++ process was accelerated inside and outside the variational Bayes inference using Nystrom’s method alongside KD-tree search. Furthermore, a downsampling strategy was implemented to manage the number of model points, standardizing on a voxel size of 0.1 across all instances. This approach to acceleration and downsampling adheres to the default parameters as suggested by Hirose (2021b).

The generic mesh of the IVD, including the morphed models (healthy and degenerated disc geometries), uses the same constitutive model, adapted by Ruiz et al. (2018) and Noailly et al. (2010). This model was implemented within an ABAQUS UMAT subroutine, allowing the biomechanical parameters to accurately reflect the properties of healthy disc tissues (Ruiz et al., 2016; Ruiz et al., 2018). This model primarily aims to isolate and examine the effects of disc morphology on biomechanical behavior, ensuring that the constitutive properties remain consistent across healthy and degenerated disc geometries.

The IVD material model considers 1) a solid phase comprising structural macromolecules such as collagen, elastin, and proteoglycans, alongside cells and 2) a fluid phase consisting of water and solutes (Malandrino et al., 2015a). The biphasic-swelling (BS) theory, as detailed by Mow et al. (1980), Mow et al. (1989), delineates both the equilibrium and transient mechanics of charged soft tissues in IVDs. This theory presents each tissue as a composite material featuring a charged solid porous phase saturated by interstitial fluid, thereby enabling the simulation of fluid pressurization and movement within the disc.

This study characterizes the behavior of the entire disc through an osmo-poro-hyper-viscoelastic model. This comprehensive model integrates the constitutive tissue of the bony endplate, treated as a linear poroelastic material (Malandrino et al., 2011) through shell elements. Furthermore, the annulus, nucleus, cartilage, and transition zone, represented with second-order hexahedral elements, employ the BS model used to simulate poromechanical interactions within a poro-hyperelastic matrix saturated with intra- and extra-fibrillar fluid (Wilson et al., 2005b), including the Donnan osmotic pressure gradient effects (Urban and Maroudas, 1981). In addition, the model considers viscoelastic collagen fibers present only in the AF Wilson et al. (2006a) as rebar elements.

The total stress tensor σ _tot is expressed as the superimposition of the effective stress σ _eff (defined in Section 2.3.1) of the solid skeleton within the pores, a fluid pore pressure component p, and Darcy’s law: σ tot = σ eff − p I (2) q = κ ∇ p (3) Where I is the identity tensor, q is the fluid mass flow to the spatial gradient of pore pressure ∇p, and κ is the hydraulic permeability tensor of the tissue. Also, the fluid flow can be expressed by: q = u f n f (4) Where u _f is the pore fluid velocity, and n _f represents the total water fraction, i.e., the porosity of the medium.

Due to the fixed charges, the cation concentration inside the tissue is higher than in the surrounding body fluid (Wilson et al., 2005b). This excess of ion particles within the matrix creates the Donnan osmotic pressure, Δπ, which drives the fluid flow, causing the swelling of the tissue (Urban et al., 1979). Incorporating the osmotic pressure into Eq. 2, where Schroeder et al. (2007) adapted this equation to the IVD, the hydrostatic fluid pressure p is defined as: p = u w + Δ π (5) u _w is the water chemical potential, linked with the pore pressure degree-of-freedom generated by the interstitial fluid permeation effects through the permeability (introduced in Section 2.3.5) by applying Darcy’s law to describe a relationship between fluid flow and the swelling pressure. Therefore, fluid flow between the different tissues of the model depends on the tissue-specific mappings of permeability. Δπ represents the osmotic pressure gradient generated by the difference between the internal and external salt concentrations (more details in Section 2.3.3).

This simulation was conducted following a free swelling step of 17 h (as per the BS theory outlined in Section 2.3.3). No external loads were applied on the morphed IVD to simulate free swelling, which stood for an initialization of the expected equilibrium osmotic pressurization of the disc, according to the Eq. 8. The pre-swelling initial values for fixed charge density and water content, denoted by c _f* and n _f* in Table 1, were set to gradually reach the c _f,0 and n _f,0 values proposed in the literature by the end of the swelling process (Ruiz et al., 2018).

Three daily load cycles were then simulated to identify which morphological factors significantly affected the mechanical response during average human activities. Accordingly, each cycle started with an 8-h resting period (creep step) under 0.11 MPa compression (load step of 10 s), which simulated overnight best rest. It was followed by 16 h of average day activity (creep step) under a load of 0.54 MPa compression (load step of 10 s), as illustrated in Figure 2A. These load values were selected based on average activity and resting intradiscal pressures as measured in vivo by Wilke et al. (1999), as proposed by Ruiz et al. (2016). The compressive loads were applied to the BEP shell elements (as mentioned in Section 2.2) of the top, while the nodes of the caudal BEP remained fully constrained (U _x = U _y = U _z = 0). The simulation also accounted for atmospheric external pressure (Figure 2B).

Given the incorporation of swelling into our simulations, it is necessary to detect potential deviations between the real geometry from the MRI models and the geometry following the morphing process, which includes swelling. To address this, we selected three IVD models with varying mid-heights (MH) for detailed analysis (See Table 2). Post-swelling, these were used to re-evaluate the Similarity Score (refer to Section 2.2 for details).

To discern possible mechanical changes because of the morphological variability in regions likely relevant to IDD pathophysiology, we focused on particular regions susceptible to being altered in early IDD. On one hand, the center of the nucleus pulposus is a natural candidate region, according to Pfirrmann’s grading. On the other hand, our mesh structure contains a region, the transition zone, that is worth exploring. This zone, as defined in our FE meshes, emerged out of a need for computational stability to ensure the FE mesh convergence and cope with the negative effects of weak discontinuities between the nucleus pulposus and the annulus fibrosus elements (Ruiz et al., 2013). Interestingly, though, a transition between the nucleus and the annulus exists, as revealed by several quantitative MRI, synchrotron imaging, cell phenotypes, and structural and composition measurements through the IVD (Marchand and Ahmed, 1990; Bruehlmann et al., 2002; Disney et al., 2023; Kapoor et al., 2023). Hence, the definition of this region in the FEM results in a more realistic description of the IVD, in contrast to an abrupt change of material properties from the NP to the AF.

Therefore, our analysis focused on the finite element simulations of the local tissue mechanics in these corresponding volumes that show signs of alterations in early-degeneration stages positioned along the sagittal plane of the disc: the transition zones (Posterior Transition Zone, PTZ; Anterior Transition Zone, ATZ), and the Center of the Nucleus Pulposus (CNP). Furthermore, we analyzed the mechanical features calculated over the central regions of the top and bottom CEP surfaces (TC and BC, respectively) to underscore the potential significance of the CEP in different disc morphologies. For each volume (PTZ, CNP, ATZ), average mechanical responses were calculated over 27 nodes contained in the volume (Figure 3E). For each surface (TC and BC), average mechanical values were computed over 529 nodes belonging to the surface (Figure 3D).

The selection of 27 nodes within each volume is grounded on the definition of a minimal-size transition zone to avoid numerical instabilities in our generic model, as defined by Ruiz et al. (2013). The TZ, characterized by a composition of 5 second-order hexahedral elements in thickness, opts for an internal volume excluding the nodes at the nucleus or annulus interface. Given the second-order nature of the hexahedral elements, the thickness includes 11 nodes, with those at positions 3, 6, and 9 being specifically chosen. Adopting a 27-node volume for the CNP facilitates a more coherent comparison by leveraging the vertex nodes of the 8 central hexahedral elements within the nucleus. Meanwhile, the CEP surface selection aims to minimize the influence of the AF, focusing solely on its relation to the NP.

Despite numerous studies using Pfirrmann’s classification to assess varying degrees of degeneration and morphological measurements, the specific morphological factors associated with IDD remain elusive. Previous studies, including observational and ex vitro research, have focused on the significance of disc height reduction at the posterior, middle, and anterior sections of IDD. These studies, however, have often presented inconclusive outcomes regarding the morphological attributes of IDD (Urquhart et al., 2014; Teichtahl et al., 2015; Bach et al., 2019; Kizilgöz and Ulusoy, 2019; Fleps et al., 2024; Tavana et al., 2024). Hence, it becomes necessary to explore the entire morphology of the discs, combining all the possible factors to understand the interplay between them.

As the morphological variations are reflected by different node coordinates in each model and the direct mechanical and tissue property environment of each region is kept consistent with the entire set of PP models, 20 morphological factors were extracted from finite element (FE) model meshes, with a consistent topology across all models. This enhances the precision of measuring heights and distances, both observable in the sagittal and coronal planes and measurable in clinical images, particularly in sagittal lumbar MRI. Further morphological features, such as the volumes, can only be calculated after 3D image segmentations but are expected to be relevant to the functional mechanics of the IVD, according to the known importance of volumetric mechanical deformations in highly hydrated materials. Therefore, we propose these features as possible novel biomarkers in the study of IDD.

Figure 4 provides an in-depth visual guide to these factors. More specifically, the explored morphological factors of the PP IVD models were: Posterior Height (PH); Middle Height (MH); Anterior Height (AH); Cartilage Endplate Top Height (CTH); Cartilage Endplate Bottom Height (CBH); Sagittal Distance (SD); Coronal Distance (CD); Nucleus Pulposus Perimeter (NPP); Sagittal Area of both the Anterior Zone and Posterior Zone of Annulus Fibrosus (AFAZSA and AFPZSA); Nucleus Pulposus Sagittal Area (NPSA); Cartilage Endplate Top Sagittal Area (CTSA); Cartilage Endplate Bottom Sagittal Area (CBSA); Volumes of Annulus Fibrosus Nucleus Pulposus, and both Cartilage Endplates (AFV, NPV, CTV, and CBV); the Wedge Angle (α), defined as the angle formed by the line connecting the upper heights of PH and AH with the line connecting the base of PH and AH. Two different ratios were evaluated as well, signifying the IVD’s intrinsic asymmetry when the sagittal and coronal plane are observed: the Posterior and Anterior Height Ratio (PAHR)—representing the ratio of PH to AH, Sagittal and Coronal Distance Ratio (SCDR)—indicating the disc’s ellipticity: PAHR = PH AH and SCDR = SD CD (21)

A brief description of the morphological features and their location in Figure 3 is detailed in Table 3.

To refine the selection of mechanical variables for analysis, our approach extends beyond principal stresses and hydrostatic pressure, elements previously examined Tavana et al. (2024); Fleps et al. (2024), to explore potential non-linear interplays with IVD morphology. Moreover, we have incorporated mechanical variables known to affect indirect mechanotransduction phenomena, such as the different strain-dependent or stress parameters related to fluid and proteoglycan contents possibly related to IDD through IVD nutritional aspects (Roberts et al., 1996; Magnier et al., 2009; Travascio et al., 2009; Shirazi-Adl et al., 2010; Malandrino et al., 2015a; Ruiz et al., 2016; De Geer, 2018).

In this context, we selected 11 key mechanical variables, modeled as the response of the 169 PP IVD models to the imposed mechanical loads as target variables. These variables were derived from averages across the three 27-node volumes (PTZ, CNP, and ATZ) and two 529-node surfaces (TC and BC), as introduced previously. They include the indirect mechanotransduction-related variables: Fluid Volume Ratio (n _w ); Pore Fluid Effective Velocity (‖ v _f ‖); Hydrostatic Pressure (p); Principal Stress I (σ _I); Principal Stress III (σ _III); Water Content (n _f ); Fixed Charge Density (c _f ); Extrafibrillar Fixed Charge Density (c _f,exf); Extrafibrillar Water Content (n _f,exf); Swelling/Osmotic Pressure (Δπ); Void Ratio (e). A description of the variables related to the constitutive model through the equations is detailed in Table 4.

We adopted three machine-learning regression models, i.e., LR, SVR, and XGBoostR, each trained using the disc morphological features to predict the target simulated variables. We partitioned our IVD dataset, consisting of features (X) and targets (y), 60% for training (X _train, y _train), 20% for validation (X _val, y _val), and the remaining 20% for testing (X _test, y _test). We chose the model and the optimal hyperparameters by maximizing the highest R-squared (r ²) value while ensuring a low mean squared error (MSE).

The implementation relies on the sklearn library in Python. For each of the five regions of interest (PTZ, CNP, ATZ, TC, and BC), we employed the following steps.

A total of 169 PP IVD FE models were successfully created out of the 169 segmented volumes. Each FE model was devoid of zero-volume elements (ZVe = 0), and all models shared the same components and connectivity as the generic FE mesh model, described in Figure 1A. The mesh quality of the generic model and the average mesh quality of all morphed discs, in terms of element aspect ratio and angular distortions, were compared by employing the ABAQUS 2020 mesh quality check functions (Figure 4A). The average aspect ratio increased by 9.4% in the morphed models, with 580 out-of-range elements and a standard deviation (std) of 274.7, compared to 530 in the generic mesh, out of a total of 19,392 hexahedral elements in each model. Small and large angles escalated by 3 (std of 14.57) and 58 (std of 118.42), respectively, in the morphed models.

None of these increases led to error elements, convergence issues, or negative Jacobians during the simulations. Considering all the out-of-range elements of the three mesh quality categories within the created cohort, it represents only 3.3% of the 19,392 elements.

The morphological measurements of the resulting morphed models are listed in Table 6 and detailed in Section 4.2. Here, we can observe that the average CEP height at the disc’s center for the top (CTH) and bottom (CBH) zones were 0.534 ± 0.148 and 0.493 ± 0.151, respectively. This agrees with that observed by Moon et al. (2013), although the out-of-range minimum and maximum values are also shown among the generated models corresponding to discs in advanced stages of degeneration.

The similarity score (Figure 4B) between each morphed model and its AF and NP counterpart surfaces in the segmented geometrical models was computed using the Hausdorff distance (as detailed in Section 2.2), as seen in step 7 of Figure 1C). On average, the AF and NP yielded a similarity score of 92.06% and 92.10%, respectively. Moreover, the morphing algorithm consistently achieved a similarity score of at least 85%.

No model was discarded for further morphing or during the simulations since all of them met the quality criteria: ZVe = 0 & S% ≥ 85%.

Figure 5 presents the time history of the vertical displacements of the cranial BEP during the 3-h creep experiment reported in Malandrino et al. (2015b) alongside the simulations of the PP model performed in this work. The simulated creep response was similar to the experimentally measured one (Malandrino et al., 2015b), with a relative error at the end of the simulation of 5.20% between our simulation and the in vitro measurements.

Our morphing process successfully produced 169 PP IVD models, with the mesh quality remaining consistent with the generic FE mesh model, even if minor increases in out-of-range elements occurred. Notably, the increase in aspect ratio criterion and small and large angles did not affect the convergence or the simulations, demonstrating the overall quality of the generated mesh and the ability of the algorithm to restrict the excessive growth of hexahedral edges. Precise deformation control was reflected by the similarity scores for AF and NP, which were 92.06% and 92.10%, respectively (Figure 4B). The BCPD++ algorithm maintained the relative distances and positions between FE mesh nodes, which prevented any single edge of a hexahedral element from growing disproportionately, achieving a patient-personalized FE model of the IVD efficiently, based on segmented medical images, with full respect to important tissue-specific mesh structures (Ruiz et al., 2013). Arguably, this approach does not allow the modeling of structural tissue damage, such as fissures or annular tears. Yet, the effect of these damages can be captured with proper homogenized continuum models and parameters, as proposed by Malandrino et al. (2015b).

As mentioned in Section 4.2, most of the research related to morphology and degeneration focuses on the diversity of disc heights. For instance, mid-height (MH) discs ranging from 5 to 16 mm have been used to determine relationships between IDD and MH (Teichtahl et al., 2015). In the present study, we examined the MH, which ranged from 5 to 18 mm. Other studies have correlated the three heights of the IVD with lumbar disc herniation (Kizilgöz and Ulusoy, 2019), with corresponding values of posterior, middle, and anterior as PH: 2–13 mm, MH: 3–15 mm, and AH: 3–18 mm, respectively. Our cohort includes discs with values of PH: 4–14 mm, MH: 5–18 mm, and AH: 6–24 mm, thus representing a diverse array of shapes that cover an extensive range of morphologies in healthy and degenerated discs. This further validates the method’s capability to preserve vital geometric characteristics.

The average CEP height at the center of the disc for the top (CTH) and bottom (CBH) zones measured 0.534 ± 0.148 mm and 0.493 ± 0.151 mm, respectively. These averages align with the findings by Moon et al. (2013) of 0.54 ± 0.12 mm, showcasing the robustness of our modeling approach. Additionally, the range of values captured includes minimum and maximum measurements corresponding to the variation expected in discs exhibiting advanced degeneration, reflecting the diversity of IVD conditions in the population.

The importance of our pipeline and methodology becomes particularly evident when compared to other techniques, such as the one proposed by Fleps et al. (2024), which also used segmented medical images to achieve personalized IVD model geometries. However, these only captured the outer surface of the AF, neglecting to personalize the Nucleus Pulposus individually. Our findings highlight the critical role of the NP in influencing indirect mechanotransduction responses. Specifically, parameters crucial for nutrient transport, including charge densities, water content, and osmotic pressure (represented by c _f , n _f , c _f,exf, n _f,exf, and Δπ), as mentioned in Section 4, especially in the posterior and anterior zone (PH and AH) and both cartilaginous regions, top and bottom (TC and BC), cannot be accurately replicated with the same level of personalized precision using such external-surface-focused techniques. Furthermore, our morphing technique has been uniquely tested on a significantly larger sample of IVD morphologies [n = 169 vs. n = 22 in Fleps et al. (2024)].

One of the primary challenges we faced in this study was the inability to replicate the morphology of the CEP due to the inherent limitation of MRI scans in delineating the external surface of this thin tissue. To overcome this obstacle, we implemented a model customization strategy. Specifically, we adjusted the dimensions of the CEP to fill the volume between the NP’s cranial and caudal surfaces and the external AF. This adaptation involves manipulating the λ variable in the third step of our morphing process, as explained in Section 2.2. Since the AF/IVD and the NP external surface are directly provided by image segmentations, a certain level of personalization of the CEP is uniquely achieved for each set of segmented disc tissue surfaces, although direct verification or validation of the FE CEP morphology was not possible. Remarkably, the algorithm allows the parameterization of the CEP thickness, in that case, adjusting the cranial and caudal external surfaces of the NP coherently. Though this feature was not used in the current study, it paves the way to systematically explore the relative importance of the CEP through models and simulations. Such explorations might support simulation-based identifications of advanced early image biomarkers associated with the risk of IDD or targeting of CEP engineering-based new regenerative strategies.

Another limitation of our study pertains to the NP-AF boundary transition zone (TZ), which, despite its relevance introduced in the Methods and noted by Ruiz et al. (2013), is difficult to define precisely in terms of material properties and size. Yet, the delineation of this region enhances the realism of the IVD model by providing a gradual material property transition from the NP to the AF, as opposed to an abrupt shift, as suggested by measurements of gradual shifts in cell phenotypes (Bruehlmann et al., 2002), MRI signal (Kapoor et al., 2023), and fiber structure (Disney et al., 2023). Arguably, this region’s local properties and exact size are difficult to assess. Even so, the TZ is physically interesting, as it represents a singular volume of increased radial compression at the periphery of the inner IVD, where the NP material is pressed against the confining annulus material because of the lateral expansion of the nucleus under external mechanical loads. This mechanism potentially influences the transport of nutrients to the cells in the TZ, as suggested by FE models and simulations (Ruiz et al., 2016; Ruiz et al., 2018). Such regional particular tissue physics coincides with the reported local emergence of relatively early signs of tissue disorganization in IVD (Leung et al., 2006; Smith et al., 2011). Moreover, the used segmentations, refined through image fusion (Castro-Mateos et al., 2014), may not have fully captured the shape of the complex NP-AF boundary, even with an improved resolution of 0.68 mm isotropic voxels. Acknowledging this limitation might be crucial to accurately represent this intricate boundary in future models, though some segmentations did represent inward bulging of the AF in degenerated discs ¹ .

Additionally, it is essential to recognize our models’ limitations related to the IVD-vertebrae interface. Arguably, the detailed mechanics of the BEP-CEP system cannot be fully captured with the current meshes, as these do not include full vertebrae.

Finally, the validation of the simulations against the research by Malandrino et al. (2015b), with a relative error of just 5.20%, demonstrated the ability to capture the real morphology with the morphing process, and the osmo-poro-hyper-viscoelastic model (referenced in Section 2.3) and the mechanical and tissue properties (listed in Table 1) in representing the time-dependent behavior of a healthy IVD under compression, as simulated hereby to represent average load cycles of daily life.

The mechanical simulation results stabilized by the second day. Examination of fixed charge densities and water content (c _f,0 and n _f ) across the morphed PP models confirmed that second-day equilibrium was achieved regardless of the particularities of the modeled morphologies. This verification of the consistency of model initialization was crucial to comparatively analyze the mechanical behaviors of the different disc models after the simulated swelling when reaching the c _f,0 and n _f,0 proposed values in Table 1.

The analysis revealed that the post-swelling Similarity Scores of the models with mid heights ranging from 8 to 17 mm remained relatively stable, as detailed in Table 7. This stability suggests that despite its potential to modify the morphed shape, the swelling process does not significantly impact the geometry. This resilience in shape demonstrates the simulations’ effectiveness in capturing the true geometry of the IVD for both the AF and NP throughout the present PP modeling process.

A key finding of this research was the correlation between the mechanical responses within specific local volumes in the IVD (Posterior Transition Zone—PTZ, Center of the Nucleus Pulposus—CNP, and Anterior Transition Zone—ATZ) and the corresponding local morphological characteristics of the disc. Overall, results suggested that local disc shapes preponderantly altered local disc mechanics and shall be considered in personalized modeling. The mechanical variables in the CEP were significantly influenced by the middle height (MH) and the AF and NP volumes (AFV and NPV), except the fluid velocity (‖ v _f ‖) that depended more on the CEP sagittal areas (CTSA and CBSA). This demonstrates the impact of the morphology of specific disc regions on biomechanical variables in other regions, probably caused by particular volume distributions, i.e., the mechanical variables in the reduced volume of the CEP are impacted by the deformation of larger volumes from other regions. Likewise, the ATZ, a region reported to be affected (Section 4.1) already in early stages of IDD (Leung et al., 2006; Smith et al., 2011), was more influenced by adjacent morphological features, such as the NP and CEP sagittal areas (NPSA and CTSA) than by its own height.

The local mechanical variables reported to impact indirect mechanotransduction phenomena, possibly involved in early IDD (see Section 4.3), are highly influenced by the corresponding local morphology. Such outcome suggests that indirect mechanotransduction is affected by disc height and local morphological factors in each evaluated zone, as we explain above. Meanwhile, maximum and minimum principal stresses (σ _I and σ _III) are respectively impacted by different morphological features. On one hand, the MH, PH, and the nucleus sagittal area (NPSA) predominantly affected the maximum principal stress, σ _I. On the other hand, the minimum principal stress, σ _III, was more influenced by ellipticity-related variables (ellipticity, SCDR, and sagittal distance, SD) and the posterior and anterior Annulus Fibrosus sagittal areas (AFPZASA and AFAZSA). These insights suggest that in-depth evaluations of IDD require a nuanced analysis that considers local morphological features. Specifically, the assessment of tensile stresses demands a detailed consideration of local axial morphologies, and the compressive stresses call for an exploration of horizontal dimensions related to the ellipticity of the disc to fully understand the mechanical interplays that might contribute to IDD.

In General, the PR% for the volumes and surfaces (Figures 8, 9) captured less than 100% of the SR, indicating that regression-based ranges were included within the simulation-based-ranges. However, in the ATZ, the velocity was particularly sensitive to the cartilage sagittal area (CTSA) and to the ellipticity of the disc (SD, CD, SCDR), so much that every top 5 features exceeded 100% of the SR. This emphasizes the importance of the cartilage endplate in critical areas in the early stages of IDD, such as the anterior zone of the IVD.

Additionally, our influential morphology results (Figure 6) align with prior experimental findings that emphasize the central role of the mid-height (MH), influencing the center of the nucleus (CNP) behaviors (Teichtahl et al., 2015; Kizilgöz and Ulusoy, 2019). However, our analysis suggests that the influence is not necessarily reflected in the magnitude of the mechanical response. The MH highly influences the behavior of the CNP, but the magnitude (Figure 8) of the hydrostatic pressure, stresses, water contents, charged densities, and the osmotic pressure do not exceed 50% magnitude variation. Nevertheless, the cumulative impact of all morphological features on both stress variations (Figure 10B) exceeded 100%, implying the necessity for a holistic morphological approach to understanding CNP mechanics.

These calculations reveal the challenge of isolating a single morphological factor as a single leading controller of the mechanical response of the disc. First, we point out possible non-linear relationships between morphology and mechanics. For example, while mid-height (MH) may be a predominant factor, its increase does not necessarily correlate linearly with the disc mechanical behavior (Figure 10A). Second, the influence of morphology is not defined by a single factor (as observed in Figures 6, 7); adjacent morphologies can significantly alter the expected mechanical outcomes by their influence and magnitudes (Figures 8, 9), thereby complicating the ability to attribute mechanical behavior to a singular morphological aspect. This complexity accentuates the need for a comprehensive customization approach of the IVD 3D morphology, including the AF and the NP, in FE simulations for achieving more accurate and realistic simulations.

Some of our morphological features, while typically not assessed in radiographic evaluations, influenced significantly the predicted internal biomechanics of the IVD and could serve, thus, as novel and potentially crucial image-based biomarkers in a clinical context. For instance, the sagittal areas of the annulus and nucleus (AFPZSA, AFAZSA, and NPSA) are readily identifiable on MRIs. Other features, such as the volumes of the annulus, nucleus, and CEPs (AFV, NPV, and CTV and CBV), as well as the nucleus pulposus perimeter (NPP), can be quantified through 3D segmentation (refer to Figure 3 for morphological factors). Current progress in medical image analysis and processing let us envision the forthcoming availability of automated tool capable of generating 3D lumbar spine models from 2D images, e.g., as proposed for osteoporosis assessments (Lopez Picazo et al., 2018).

Although the present study uniquely revealed the intricate interactions among intervertebral disc morphology and internal biomechanics, it has some limitations. Our boundary conditions represent a normal subject with diverse physiological physical activity during the day and lying at night, as reflected by in vivo measurements of intradiscal pressures (Wilke et al., 1999). Although this should be studied following a person’s typical activities, the current approach serves as a standard loading scheme to evaluate a large number of PP models. Since the objective was to focus only on morphology, material properties specific to each degree of degeneration have not been considered. This limitation shall be addressed in future studies, in particular, seizing the opportunity to incorporate Pfirrmann grade-specific material properties (Malandrino et al., 2015b; Barthelemy et al., 2016; Ruiz et al., 2016). Performing morphing for both the AF and NP was challenging due to its difficulty. However, we have achieved it with acceptable levels of similarity to the real model without damaging the mesh, becoming one of the strengths of our work.

The tensile and the compressive stresses are, respectively, the maximum and the minimum principal stresses, which are the local eigenvalues of the stress tensor. Yet, we disregard the analysis of the orientations of these eigenvectors because regardless of the latter, the stress is, indeed, suffered by the materials. Looking at the local orientations of the eigenvectors would be interesting to analyze for future work.

Although the generated IVD files (.inp) are ABAQUS files, these are in ASCII format and can be converted so that models can be reused with alternative open-source software suites, such as FEBio (Maas et al., 2012), MoFEM (Kaczmarczyk et al., 2020), High-Performance Computing (HPC), among others.