Our newly developed TO strategy (Smit et al., 2024) (Smit, 2023) to optimise spinal fusion cages was, in this work, slightly modified to improve its effectiveness and to allow the optimisation of titanium and PEEK implant materials. The computational domain Ω comprised an implant domain (Ω_implant), a bone domain (Ω_bone), and a rigid domain (Ω_rigid), with Ω = Ω_implant ∪ Ω_bone ∪ Ω_rigid (Figure 1A). Ω_bone contains the bone structures and patient-specific material properties. Ω_implant was the design domain that was subjected to optimisation in the TO process.

Domain Ω, with a size of 54 mm × 54 mm × 54 mm, was discretised using a structured grid of 1.728 million hexahedral elements with an isotropic element edge length of 0.45 mm. Resampling was performed using SimpleITK version 2.1.1.2 (Lowekamp et al., 2013) to adjust the CT scan resolution to match the resolution of Ω and to position the Functional Spinal Unit (FSU) in the approximate centre of Ω. The vertebrae were partially embedded in Ω_rigid.

Six load cases (axial compression, lateral shear, posterior–anterior shear, flexion, lateral bending, and axial rotation) were used in the TO process. The applied loads were based on the study by Rohlmann et al. (2008) and the OrthoLoad dataset (Bergmann and Damm, 2008) (Figure 1C; Table 1). Loads were applied on nodes at the top of Ω. Nodes on the bottom of Ω were fixed in all directions. The load magnitude for all patients was scaled with respect to the patient’s body weight and the body weight of patient 1 (see Section 2.2), as was used in the subsidence risk assessment (Han et al., 2013) (Zhang et al., 2017).

Similar to the TO formulation that was previously used (Smit et al., 2024) (Smit, 2023), the compliance of the bone–implant system was minimised for all load cases (optimisation objective), but to avoid the overloading of the adjacent vertebrae, the maximum principal strains in Ω_bone were constrained (optimisation constraints). Constraining the minimum principal strains was omitted because our previous work indicated that this constraint was inactive (Smit et al., 2024) (Smit, 2023) because the compression load from the implant is transferred to the shear, orthogonal to the endplate, along the periphery of the implant. In pure shear, the maximum and minimum principal strains have equal magnitudes but opposite signs, indicating that the tensile limit is reached first due to a lower threshold. Shearing is the dominant loading mode in the endplates, which was confirmed by several authors who identified shear failure of the vertebra endplates as the dominant failure mode in subsidence (Cadman et al., 2016); (Au et al., 2011). Furthermore, we use a local volume constraint on Ω_implant to create a porous implant structure. To promote manufacturability, the minimum strut diameter was controlled using robust TO formulation. No other manufacturing-specific constraints were added to make a comparison between titanium and PEEK fair because these materials have different manufacturability requirements.

The compliance minimisation problem with the maximum principal strain constraint and local volume constraint in discrete form is written as min x ∑ k = 1 l c f x e k   i n   Ω subject to (1) g x d 1 ≤ 0.0   i n   Ω b o n e , (2) v x d ≤ 0.0   i n   Ω i m p l a n t , (3)

where l c represents the six load cases with k = 1 … 6 , f is the objective function (Equation 1), g 1 is the constraint on the maximum principal strain response (Equation 2) for the compression load case in the Ω_bone, and v is the local volume constraint (Equation 3) on Ω_implant. The maximum principal strain constraint on the compression load case was used as this was the active constraint (Smit et al., 2024) (Smit, 2023) with a maximum principal strain limited to 0.73% based on the work by Bayraktar et al. (2004). The local volume constraint v x d followed the implementation published by Wu et al. (2018). The local volume constraint was applied on the dilated design x d , with a local volume fraction of 0.35, a filter radius of 3.5, and p-mean penalty of 16. Finally, Poisson´s ratio for Ω_implant was set to 0.3 for titanium and 0.43 for PEEK. The E-modulus for the titanium implant was set to 110 GPa (Niinomi and Nakai, 2011), and for the PEEK implant, the E-modulus was set to 3 GPa (Zhao et al., 2018). The optimisation was run for 24 h on a computational cluster using 1,000 cores in parallel. The final design should at least exhibit a discreteness measure that is lower than 3%. This value is considered to be an acceptable convergence to a solid/void design (Wang et al., 2011).

The resulting cage designs from the optimisation process have a voxelised surface because of voxel-based discretisation. The final designs were achieved by applying the “Iso volume” filter to the optimised cage using a contour of 0.5 as the threshold and exporting to an .stl file using vtk 9.1.0 (Schroeder et al., 1996). The resulting file was used in the subsidence risk assessment Section 2.3 and Section 2.4.

Anonymised pre-operative computed tomography (CT) scans with clinical resolution (voxel sizes of approximately 0.30 mm × 0.30 mm × 0.50 mm), from seven female patients, were provided by the Schulthess Klinik, Zürich (Table 2). The patients, who were diagnosed with degenerative spondylolisthesis or degenerative disk disease, were treated with a TLIF procedure. Patients with other spinal diseases such as trauma, infection, or tumour were excluded. The CT images from patients 1 and 2 were reused from a previous study (Smit et al., 2024) (Smit, 2023), and all the patients provided informed consent for the use of their data for research purposes. The CT scans were calibrated using a phantom-less calibration procedure (Lee et al., 2017). The vertebrae were segmented using open-source software MITK-GEM (Pauchard et al., 2016). T-scores, based on DXA scans, were only available for the L4 level for three patients. Bone quantity was thus assessed by calculating the integral volumetric bone mineral density (vBMD) for each vertebra that was modelled, including the endplates and cortical bone in accordance with the work by Kaiser et al. (2020) (Table 2; Figure 2).

The model building process follows the workflow presented by Smit et al. (2024) (Smit, 2023) but is briefly described here for clarity and context. A commonly used, commercially available OTS cage, for a TLIF procedure with a height of 9 mm, was reengineered into a CAD model. The cages were placed anteriorly and across the mid-line in an optimal position, restoring spinal alignment as much as possible. The implant placement was checked by a spine surgeon according to standard surgical fusion techniques for the TLIF approach. The ACD and AMCD design domains were created by matching the top and bottom surfaces of the CAD model to the patient’s endplates. It was assumed that the intervertebral disc and cartilaginous endplates were removed during the surgery. The cages were placed according to standard surgical fusion procedures, and the implant position was checked by a spinal surgeon.

In the optimisation process (Section 2.1), daily living load cases were applied, which cause strains in the bone structures within the linear regime. The optimisation process was followed by the subsidence risk assessment (Section 2.4) that was performed in commercial FE analysis software where the cages were loaded with hyper-physiological loading that can potentially lead to non-linear material response, both in the implant and the bone structures.

In the subsidence risk assessment, 10-node tetrahedral elements (C3D10M) were used to mesh the bone and implant structures with an average edge length of 1.5 mm for the vertebrae and 0.5 mm for the cages. The apparent density-to-modulus relationship from the study by Ouyang et al. (1997) (Table 5) was used to assign material properties to the bone elements in the FE mesh based on an internal calibration of the CT grey levels to the apparent density. Poisson’s ratio was set to 0.3, and the minimum value for Young’s modulus was set to 25 MPa for the bone structures. The subsidence risk assessment was performed in an explicit non-linear FE solver (Abaqus 2021, Dassault Systèmes, Vélizy-Villacoublay, France), assuming non-linear geometry and a general contact formulation. The bone structures were modelled using a rate-independent elasto–plastic material model, which captures the asymmetric tension–compression post-yield behaviour of bone (Figure 3B) (Bayraktar et al., 2004). A compressive follower load (FL) of 310 N for patient 1 was assumed based on the study by Rohlmann et al. (2014) and the OrthoLoad dataset (Bergmann and Damm, 2008). The magnitude of FL for all patients was scaled with the patient’s body weight to the body weight of patient 1 and subsequently distributed on the superior endplate of the superior vertebra such that the direction of the load was as defined by Rohlmann et al. (2009) (Figure 1B). The inferior endplate of the inferior vertebra was constrained in all translations. The models were solved on a computational cluster using 64 cores in parallel. Poisson’s ratio was set to 0.3 and 0.43 for titanium and PEEK, respectively. The material models for titanium and PEEK are shown in Figure 3A.

The maximum factor of fracture risk (FFR) was used as a surrogate for subsidence risk. The maximum FFR is quantified as the degree of overloading of the vertebrae using a maximum principal strain limit of 1.5% and a minimum principal strain limit of −2.0% (Soyka et al., 2016). For each element “ i ” in the bone domain, the FFR ( F F R i , Equation 4) and maximum FFR (Equation 5) (Soyka et al., 2016) were calculated using F F R i = ⁡ max ε max ⁡ , i 1.50 % , ε min ⁡ , i − 2.00 % , (4) M a x i m u m   F F R = ⁡ max   F F R i , (5) where ε max ⁡ , i is the maximum principal strain of element i , in the centroid of the element, and ε min ⁡ , i is the minimum principal strain of element i , in the centroid of the element.

Credibility assessments were performed following the ASME V&V40–2018 standard. Note that the following applies to the subsidence risk assessment of Section 2.3, 2.4. The process starts with the definition of the question of interest (QoI). The next step was to formulate the context of use (CoU) including the function and scope of the computer model that was used in answering the QoI (Table 3).

The model risk describes the possibility that the use of the computer model leads to a decision causing harm to the patient. The model risk is composed of the “model influence” and the “decision consequence.” The model influence describes how much influence the results of the computer model have on the decision maker.

The model influence is defined to be MEDIUM because the decision of the AMCD has lower subsidence risk than the ACD or OTS, for a specific patient, is informed by the results from the computer model, in combination with the observations and judgment of the spinal surgeon.

The decision consequence is the significance of harm to the patient in the case of an incorrect decision.

The decision consequence is defined to be LOW because a worst-case incorrect decision would lead to a patient being treated by an OTS cage, which is the current standard practise.

The selection of computer model credibility activities is based on the credibility requirements for the computational model that should correspond with the model risk that was defined (Table 4). Model credibility is defined as the trust in the predictive capability of the computer model for the CoU, according to the ASME V&V40–2018.

The ASME V&V40 standard suggests several credibility factors that should be evaluated to demonstrate that the overall credibility of the computational model agrees with that of the model risk (Table 8).

A mesh convergence study of the bone structures was performed by calculating changes in the maximum FFR for gradually decreasing element edge length. Mesh convergence of the implant structures was assessed by decreasing the element edge length of the implant mesh while maintaining an element edge length of 1.5 mm for the bone structures. The maximum FFRs were calculated, as well as the percentage difference in the maximum FFR between different meshes ( c u r r e n t − p r e v i o u s p r e v i o u s * 100 ). The simulations were run using patient 2 and the ACD in titanium. Furthermore, for the mesh convergence study, a homogeneous Young’s modulus of 350 MPa for the bone structures was used.

The sensitivity analysis of the maximum FFR with respect to the load magnitude was performed by calculating the variation in the maximum FFR while applying variations (2.5% and 25% change), positive or negative, to the applied load magnitude. The sensitivity analysis of the maximum FFR to different apparent density-to-modulus relationships was performed by calculating the maximum FFR using two alternative relationships, from the study by Kopperdahl and Keaveny (1998) and by Morgan et al. (2003), with all other model parameters kept the same (Table 5). The sensitivity of the maximum FFR with respect to the constitutive model for the subsidence risk assessment was investigated by additionally calculating the maximum FFR with a linear constitutive model and a symmetric post-yield behaviour for the bone structures with all other model inputs equal. The sensitivity studies were conducted for patient 2 with the ACD in titanium and the AMCD of patient 6 in PEEK. Percentage differences were calculated compared to the base model, and the maximum percentage differences were reported. In our previous study, mechanical tests were conducted, and the computer model was subsequently validated against them (Smit et al., 2024) (Smit, 2023). A more in-depth explanation of the credibility factors is provided in the work by Viceconti et al. (2016) and Aldieri et al. (2023).

For statistical analysis, Spearman’s rank correlation coefficient analysis was used to calculate the correlation coefficient r with (−1 < r < 1) between two variables. A Mann–Whitney U test was used to compare two samples with statistical significance assumed when p < 0.05. The small sample size in this study does not allow for reliable normality testing (Taylor, 2017); therefore, statistical methods were selected that allow for non-normally distributed data and are robust on small datasets, e.g., the median was reported. SciPy version 1.11.1 (Virtanen et al., 2020) and Seaborn version 0.11.2 (Waskom, 2021) were used to implement the statistical methods and data visualisations.

The resulting AMCDs show a “box-like” design and are hollow in the inside with patient-specific internal strut structures, similar as found by Smit et al. (2024) and (Smit, 2023). The outside shape conforms to the endplate shapes of the vertebrae. The designs satisfy all the constraints and reached a discreteness parameter lower than 3%. The morphological details of the AMCDs are listed in Table 6. The porosity and mean pore size diameter are very similar for the titanium and PEEK implants, the reason being that the local volume constraint, which controls the porosity and pore size distribution, is an active constraint. Nevertheless, the internal structures and the exact location of the pores are patient-specific. The fact that the local volume constraint is effective in controlling the implant morphology, independent of other constraints, enables it to tune the porous internal structure of the implant to be favourable for bone in-growth.

Figure 4 highlights elements with FFR larger than 1, which indicates mechanical overloading, for patient 1 for all three implant types. The red elements are associated with the titanium implant and the green elements with the PEEK implant. Figure 4 shows that the titanium implant leads to a higher volume of overloaded elements than the PEEK implant material by 15%, 5%, and 300%, for the OTS cage, ACD, and AMCD, respectively. Furthermore, the two separate concentrations in the superior vertebra compared to the evenly distributed overloaded elements in the lower vertebra are explained by the contact areas between the endplate and the OTS cage (Figure 4A). The ACD and AMCD were endplate conforming and show that the regions of overloaded elements were reduced (Figures 4B, C).

The subsidence risk assessment provided the maximum FFR data for the three implant groups in PEEK and titanium implant materials for all patients (Figure 5A). In general, the subsidence risk for ACDs was reduced, compared to the OTS cages, and the subsidence risk for AMCDs was reduced, compared to the ACD, for all patients.

Comparison of the maximum FFR results between PEEK and titanium implant materials, for the three implant types, showed that there was no statistical difference between the PEEK and titanium groups with p-values of 0.71, 0.95, and 0.53 for the OTS, ACD and AMCD implant types, respectively. The negative correlations between the maximum FFR and vBMD for the different implant types are given in Figure 5A. As an illustration, the linear correlations between the maximum FFR and vBMD indicate that, for a patient with a mean vBMD of 0.25 g/cm³, the maximum FFR values would be 14.8, 3.8, and 2.0 for the OTS cage, ACD, and AMCD, respectively.

A comparison of the implant type and implant materials is given in Figure 5B. This shows that the median maximum FFR reduction between the OTS PEEK to AMCD PEEK, OTS Ti to AMCD Ti, and ACD PEEK to AMCD PEEK is statistically significant with 94%, 89%, and 75%, respectively. The reduction in the maximum FFR for the OTS cage to the ACD, ACD to AMCD, and OTS cage to AMCD is not strongly significantly correlated with vBMD, with p-values of 0.18, 0.07, and 0.09 for the PEEK implant material and 0.12, 0.15, and 0.09 for titanium, respectively.

An element edge length of 1.5 mm was chosen for the bone structures because subsequent decreases in element edge length produced percentage differences in the maximum FFR lower than the 5% threshold (Figure 6A). For the implant structures, an element edge length of 0.5 mm was chosen, and the mesh convergence study indicated that variations had little influence on the maximum FFR (Figure 6B).

Table 7 shows the results of the sensitivity studies. The sensitivity of the maximum FFR with respect to variations in load magnitude, different density-to-Young’s modulus relationships, and different constitutive models is illustrated by reporting the maximum percentage difference derivation from the base case. Patient 2 is loaded in the linear loading regime and shows little variation in the maximum FFR for the applied variations. On the other hand, patient 6 shows a larger variation in the maximum FFR because this patient is loaded in the non-linear loading regime.

The credibility factors were assessed after the outputs of the verification and validation (V&V) activities were obtained. The robustness of the activities and level of rigour were evaluated according to the guidelines and examples from the ASME V&V40–2018 standard. This assessment is summarised in columns “Rigour selected” and “Achieved credibility” in Table 8.

The quantities of interest considered in the verification and validation study were the maximum FFR and the force–displacement curve in the mechanical tests. Therefore, we consider the quantities of interest relevant but not identical since the maximum FFR cannot be calculated directly in mechanical testing of samples. In the validation study, the range of applied forces is physiological and, therefore, deemed relevant for the CoU. The mechanical tests did not include the bone structures; therefore, partial overlap between the validation and CoU is concluded.