Diagnosing spinal disorders and their underlying causes has long posed a substantial challenge, as recently evidenced again by a study that revealed significant discrepancies in diagnoses from different MRI centers based on one single MRI (Herzog et al., 2017). To offer individual prevention strategies, estimate personal risks or plan patient-specific therapies, more information needs to be derived from medical images. One important step in this process is the analysis of spinal biomechanics, which can, for instance, be done by experimental measurements or numerical simulations such as multi-body simulation (MBS) and the finite element method (FEM). While experimental setups are often limited to few subjects, numerical simulations can more easily be realized patient-specifically, as well as in large datasets. The latter aspect bears the potential to deliver new pathophysiologic insights (Lerchl et al., 2024), for instance, relating certain spine characteristics to IVD degeneration. This could eventually result in better prediction of disease progression. However, analyzing large datasets requires two key components: patient-specificity and automation.

During the last decade, valuable achievements have been made towards this goal. Automated medical image segmentations have profited from the rise of machine learning (Sekuboyina et al., 2021; Moeller et al., 2024) and started to serve as a basic prerequisite for numerical models. Considering MBS, Lerchl et al. published an automated pipeline for the calculation of compression and shear forces between vertebrae in musculoskeletal models, including individualized body weight, spinal alignment, and the attachment points of ligaments and muscles (Lerchl et al., 2022). Geometries were derived from computed tomography (CT) images using machine learning (Lerchl et al., 2022; Sekuboyina et al., 2021). In contrast to the surface meshes defining rigid bodies in the MBS, FEM meshes require solid 3D meshes to calculate the inner stress and strain in vertebrae or intervertebral discs (IVD). The methods to derive patient-specific FEM meshes from medical imaging data can broadly be divided into voxel-based mesh generation, mesh morphing or surface mesh filling. Multiple approaches have been presented in the past decade, which apply one of these methods in conjunction with different grades of automation (Castro-Mateos et al., 2015; Castro-Mateos et al., 2016; Malandrino et al., 2015; Clouthier et al., 2023; Campbell and Petrella, 2015; Campbell et al., 2016; Caprara et al., 2021; Lavecchia et al., 2018).

In voxel-based models, segmented voxels directly serve as hexahedral elements of an FEM mesh (Castro-Mateos et al., 2016).

Mesh morphing is based on the deformation of an FE template mesh towards landmarks (Campbell and Petrella, 2015; Lavecchia et al., 2018; Caprara et al., 2021; Fasser et al., 2022) or complete surface meshes (Castro-Mateos et al., 2015; Caprara et al., 2021; Fasser et al., 2022) that were identified in individual medical images. The morphing process inherently involves striking a balance between accuracy and mesh quality. Non-conforming anatomies of IVDs and vertebrae can only be adapted by template meshes through large deformations. However, large deformations result in poor mesh quality, likely leading to convergence issues during the final FEM simulation (Malandrino et al., 2015). To minimize the displacements between a template and target mesh, recent approaches used statistical shape modeling (SSM) for their generation (Clouthier et al., 2023; Campbell and Petrella, 2015; Caprara et al., 2021; Castro-Mateos et al., 2015). However, ensuring sufficient mesh quality remains one of the main challenges in morphed patient-specific spine models.

These limitations do not apply to the filling approach, in which surface meshes are derived directly from segmentations in an initial step, which is often done by applying the marching cubes algorithm (Lorensen and Cline, 1987). Acquiring MRI scans involves capturing multiple slices of one to 3 mm thickness from various perspectives, partly leading to limited spatial resolution and voxels with heterogeneous properties. This inherent characteristic is translated directly to the generated surface meshes, limiting their accuracy and demanding subsequent smoothing (Molinari and Falcinelli, 2021). In the literature, limited information is provided on applied smoothing protocols. However, low resolution or artifacts like stair steps have been compensated using overall vertebra smoothing (Mönch et al., 2010; Mönch et al., 2013). This results in a loss of detail, whereby features like vertebra edges or osteophytes are smoothed to the point of disappearance, which we here refer to as oversmoothing. Consequently, calculated FEM stresses might underestimate actual stresses, potentially leading to the overlooking of wedge fracture predictions (Nispel et al., 2024). To avoid over-smoothing and enable accurate biomechanical results, selective smoothing algorithms have been introduced in recent studies (Nispel et al., 2024). Once smoothed, surface meshes are filled with volume elements, which is usually done manually and involves computer-aided design (CAD) software (El Bojairami et al., 2020; Lavecchia et al., 2018; Zadpoor and Weinans, 2015). Although for the morphing approach, efforts have been made towards automated, patient-specific model generation, no algorithm has yet been published to automate the filling approach.

With MRI images being more challenging to segment automatically, inter alia due to their subtle contrast variations, automatically generated spine FE models are mainly based on CT. They thus often include an approximation of the IVD shapes, neglecting anatomical characteristics such as bulge (Campbell and Petrella, 2015; Lavecchia et al., 2018; Caprara et al., 2021). In addition, material parameters for the IVD, which vary significantly depending on the degree of degeneration, cannot be extracted from CT images. Recent achievements in fully automated MRI segmentations opened the door to geometrically accurate patient-specific FEM models for vertebrae and IVDs including varying material definitions for soft tissue (Moeller et al., 2024; Graf et al., 2023).

However, the derivation of both IVDs and vertebrae poses another challenge: Contact modeling becomes more complex if the contact surfaces of both bodies are not equivalent or do not share nodes. This diminishes computational efficiency, potentially resulting in less accurate outcomes (Papadopoulos and Taylor, 1992) and hinders coupled MBS and FEM simulations (Nispel et al., 2023). To better understand the challenges of realizing shared contact surfaces, it is important to note that smoothing typically involves specific node translations for each body and can only be applied to bodies separately. Thus, smoothing steps need to be implemented in a specific order and with respect to the preservation of the adjacent body’s node positions (El Bojairami et al., 2020).

In summary, methods either automate the process using the morphing approach or rely on the filling approach in combination with a manual procedure. To the best of our knowledge, there is currently no existing approach that automates the filling approach.

Addressing the need for automated approaches to realize large dataset investigations, we here present the first pipeline to create and simulate patient-specific FEM models of vertebrae and IVDs with shared-node contact surfaces from MRI segmentations in a completely automated manner. Note that the focus of this work is on methodological development, particularly in automating model creation. As such, complex material models or biologically related parameters such as varying loading conditions or patient-specific weights are deliberately excluded from this study. The FEM simulation presented here serves solely to verify that the pipeline produces models capable of converging under simplified conditions, laying the groundwork for future biological model validation.

We implemented an automated pipeline that takes segmented MRI images of the spine as an input and provides FEM simulation results as an output (Figure 1). All steps were automated using Python as a baseline programming language. We used an in-house dataset that included 30 patients. The data contains MRI images of a 1 × 1 mm resolution in the sagittal plane and a 2.5–3.5 mm slice thickness. The implemented pipeline can broadly be divided into the following substeps: Generation of surface meshes from MRI segmentations, smoothing, mesh filling, volume meshing, FEM modeling and FEM simulation (Figure 2). FE models can be used in a plug-and-play manner to simulate either single bodies, functional spinal units (FSU) or complete spine models. To demonstrate the functionality of the FEM models, we included the automated execution of FEM simulations as the final substep. Therefore, we defined an exemplary load and material model. Note that this work focused on the automation aspect and the resulting FEM stresses and displacements were not interpreted biomechanically. Given that a manual approach would not significantly affect the key pipeline steps—mesh smoothing, volume body creation, or FEM meshing—and that the manual process is highly user-dependent, making precise time comparisons challenging, we chose not to include a traditional manual control group in this study. In what follows, we present a detailed description of the substeps.

We implemented a pipeline that is able to automatically calculate FEM results of vertebrae and IVDs in large cohorts. The pipeline is based on an automated segmentation of vertebrae and IVDs in MRI images (Moeller et al., 2024). From this, surface meshes were derived and selectively smoothed to mimic the biological endplate shapes, compensating for image resolution inaccuracies. Hereby, we take advantage of the MRI segmentations by including both bone and soft tissue. Smoothed surface meshes were then automatically transformed to FEM volume meshes, which contain individual node sets of the superior and inferior contact surface to the adjacent vertebra or IVD, respectively. FEM models were available as a combination of nodes and elements. In our study, these models were supplemented with boundary conditions, loading and material definition to create and simulate FEM models. Using this pipeline, we processed all 30 patient MRI scans, resulting in a total of 637 vertebrae FEM simulations and 600 IVD FEM simulations.

637 vertebrae and 600 IVDs were modeled and simulated using our automated pipeline, with an average duration of 4.4 min per vertebra and 31 s per IVD. We evaluated the quality and performance of the pipeline by investigating the quality of generated meshes, failed attempts and processing times.

With the predefined mesh quality criteria, more than 98% of generated models achieved good quality. This was also reflected in the low amount of failed attempts - only ten vertebrae and seven IVDs raised an error within the modeling and simulation process. A detailed analysis of these failed attempts revealed insights into the failure mechanisms in our pipeline. For all cervical IVDs, errors were due to self-intersecting faces or non-manifold edges in the surface mesh, which occurred in regions of small segmentation volumes. Failure mechanisms in thoracic IVDs were connected either to a hole in the raw mesh or small segmentation volumes. In one case, a hole led to a T-vertex in the presmoothed mesh, which caused the self-intersection of faces in the adaptively smoothed IVD. In the only failed lumbar IVD, T-vertices were created for an unidentified reason regarding the good quality of the raw mesh. We additionally monitored all failed vertebrae for apparent quality issues, and found similar failure mechanisms as in the IVDs. In six of the ten vertebrae, self-intersecting faces were found, which led to failing FEM mesh generation. In the remaining cases, mostly cervical vertebrae, the segmentation mask contained two disconnected volumes. Our pipeline was able to generate two disconnected surface mesh parts from these segmentation masks. However, only one volume was converted to a solid mesh. We consider this an issue of the segmentation mask and not primarily a limitation of the here presented pipeline.

The processing time for the FEM simulation of vertebrae was approximately tenfold higher compared to that for IVDs. It is worth noting that meshing and smoothing times are directly proportional to the number of elements in the mesh, which is roughly 4.5 times greater in vertebrae models than in IVD models.

Assessing the efficiency of our approach versus traditional manual methods is challenging due to the lack of qualitative literature data on the time required for manual implementation. Manual processes typically involve segmentation software for surface mesh derivation, followed by volumetric model generation and FEM simulations using separate tools, a process that can take up to several days (Caprara et al., 2021). In contrast, our automated pipeline drastically reduces this timeframe. With no user interaction and reliance solely on an MRI image, our method achieves biomechanical analysis of a vertebra or IVD within just about 4.3 min or 27 s, respectively. In addition, it is insusceptible to variability due to manual steps.

One limitation of our pipeline is the homogeneous, isotropic modeling approach, that is in conflict with recent IVD and vertebra models, which contain anisotropic material definitions in the case of vertebrae (Molinari and Falcinelli, 2021; Fleps and Morgan, 2022), and biphasic, hyperelastic or fiber-reinforced components in the case of the IVD (Gruber et al., 2024; Dreischarf et al., 2014).

For the IVD, these more complex modeling approaches go hand in hand with anisotropic meshing that comprises, in the best case, hexahedral elements (Dreischarf et al., 2014). The presented approach does not enable a structured mesh including multiple phases, such as e.g. a nucleus pulposus cavity or annulus fibrosus fibers. Besides the lack of detail this brings to the resulting mechanics, calculation times would likely increase when using a more realistic material definition. However, experiments indicate that already in mildly degenerated discs, fluid content decreases significantly (Galbusera et al., 2014; Iatridis et al., 1998; Iatridis et al., 1997; Johannessen and Elliott, 2005). As a result, tissue behaved almost linearly and a change from anisotropic to isotropic nature was reported (Iatridis et al., 1998). For lumbar IVDs, a slowly increasing loss of flexibility correlates with an increasing grade of degeneration (Kettler et al., 2011). Interestingly, studies have found an abrupt decrease in flexibility for mildly degenerated thoracic discs with Pfirrmann grade (Pfirrmann et al., 2001) of 1 (Liebsch and Wilke, 2022) as opposed to healthy discs. Precisely, the tetrahedral representation of IVDs, as it results from our pipeline, might be sufficient for analyzing further degenerated IVDs with the advantage of low complexity and comparably short computation time. In addition, biomechanical FEM studies of multiple FSUs that investigate the holistic behavior of the spine, like the one conducted by El Bojairami et al. (2020), may profit from the reduced complexity of the IVD mesh. To particularly investigate the advantages of a more complex IVD model in their use case, El Bojairami et al. implemented an additional two-phase fluid-structure model besides their isotropic one. And extracted the hydrostatic pressure from an enclosed, hydrostatic pressure element within the nucleus pulposus cavity (El Bojairami et al., 2020). The comparison between the two modeling approaches (isotropic vs. complex) showed that both methods produced results that overlapped to a great extent, with a maximum discrepancy of approximately 4% at certain flexion angles. Moreover, tetrahedral elements allow for a more straightforward implementation of shared-node contact surfaces, which is beneficial for computation times and numerical stability (El Bojairami et al., 2020). A recent observation in FEM simulations of degenerated discs even found that tetrahedral elements were more stable and accurate in the results (Fasser et al., 2022). Altogether, these findings may support employing tetrahedral elements and isotropic material properties for the IVD in limited circumstances, particularly when considering factors such as computational efficiency and simplicity of implementation. However, multiple studies have analyzed the benefits of biphasic or fiber-reinforced material models for the IVD (Gruber et al., 2024; Stadelmann et al., 2018; Knapik et al., 2022; Sciortino et al., 2023). A differentiated consideration of modeling approaches is crucial, and their application should be carefully selected depending on the specific research question.

Apart from the IVD model complexity, vertebrae models in our pipeline are also limited to isotropic, tetrahedral meshes, which lack the distinction between cortical and cancellous bone. However, the inclusion of shell elements could easily be realized similarly to the approach of Imai et al., who attached a 0.4 mm thick layer of shell elements to the outer surface of the tetrahedra (Imai, 2015). Nevertheless, the application of modeling approaches should also be carefully selected based on the specific research question. By providing the pipeline to the public and ensuring the output of a standardized mesh representation as an .inp file, research groups are enabled to further process and adapt meshes of their own datasets towards their requirement.

In comparison to morphing approaches for FEM model generation (Campbell and Petrella, 2015; Lavecchia et al., 2018; Caprara et al., 2021; Fasser et al., 2022; Castro-Mateos et al., 2015; Caprara et al., 2021; Fasser et al., 2022), the here presented filling approach represents more details of the segmentation, which unfortunately includes low-resolution features like stair steps as well. To counter this, we implemented a tailored smoothing algorithm Nispel et al. (2024). However, since the focus was placed on smoothing the interface surfaces, the circumferential surfaces of vertebrae and IVDs were only considered during the pre- and postprocessing of smoothing, resulting in stair-step artifacts in the final meshes. This may be viewed as a limiting factor of our approach. With advancements in faster MRI scanning methods, higher-resolution scans will become more readily available, potentially reducing these stair-step artifacts in the circumferential surfaces. Considering the final FEM mesh quality, we do not expect any improvement in the FEM mesh with the enhancement of MRI resolution, as the size of the FEM elements is much smaller than the stair-step artifacts.

Finally, given the MRI segmentation in our pipeline, hard tissue prediction is not as stable as in CT image segmentation. Potential errors in correctly predicting detailed vertebrae characteristics like osteophytes could be propagated to the FEM results and may affect the result quality negatively.

The goal of this work was to advance the automated generation of patient-specific FEM models derived directly from MRI segmentations. In comparison to recent approaches, our approach is based on the MRI image itself and does not rely on templates like SSMs. The process of morphing a mesh is thereby obsolete, which constitutes the advantage of having fewer processing steps and higher flexibility in representing geometries that do not fit into statistical norms, like osteophytes, strongly deformed vertebrae or extreme bulges in IVDs.

By realizing surfaces with shared nodes and elements for adjacent vertebrae and IVDs, we spare the time-consuming and unstable process of implementing penalties during contact modeling. In addition, shared interface nodes are advantageous in the framework of exchanging load and displacement data in coupled MBS and FEM simulations.

Manual implementations of full spine models may be subject to automation efforts in the future. Constitutive models of the IVD and vertebrae can be advanced by including locally varying, individualized MRI- or CT-derived material parameters. To include both detailed bone geometries with, e.g., osteophytes and IVD deformities such as bulges, a co-registration of MRI and available CT data could be beneficial.

For the IVD, a differentiation between nucleus pulposus and annulus fibrosus, together with a fiber-reinforced or biphasic implementation could be beneficial in answering specific research questions. Vertebrae might benefit from a differentiation between cortical shell and trabecular bone. Prospectively, combining our pipeline with validated models, as demonstrated in our group’s recent work Gruber et al. (2024), will enable large cohort studies to gain insight into the causes of spinal disorders such as degeneration or back pain.