Due to their advantageous properties, ceramics are used in TDRs (Kölle et al., 2022). The established bioceramic zirconia-toughened alumina material, BIOLOX®delta (CeramTec GmbH, Germany), was chosen, as it offers low wear (Döring et al., 2023), biocompatibility (Asif, 2018), high strength (Kuntz et al., 2009), MRI (Mödinger et al., 2023) and CT compatibility, and years of successful clinical use in hip arthroplasty (Alshammari et al., 2023). The design concept of the bearing geometry is based on toroidal and cylindrical shapes subtracted from the inferior part, and a spherical shape on the superior part of the TDR (Figure 1).

An implant height of 6 mm was chosen, as the mean disc height of C6/C7 was reported between 5.48 (Choi et al., 2017) and 6.4 mm (Yukawa et al., 2012). Appropriate implant height is necessary to enable a sufficient range of motion and avoid an increase of facet joint pressure (Wang et al., 2021). Dimensions of the vertebrae were taken from (Bazaldúa Cruz et al., 2011). A radial clearance of 0.07 mm was chosen. The baseline design, with which the optimization started, had a ball radius of 7.2 mm and torus radii of 3 mm inspired by the Prestige™ LP (Medtronic, United States) (Jung et al., 2013).

A finite element model (FEM) was built in LS-PrePost V4.6.4 (LSTC, United States) containing only the necessary structures (Figure 1) to minimize computational cost. The effects of other structures are represented by applying the kinematics of the intact spine and by using the compressive forces predicted for the intact spine. The material parameters of BIOLOX®delta were provided by the manufacturer (CeramTec GmbH, Germany) and a linear elastic material model was used for the TDR with a density of 4.37 g/cm³, a Young’s modulus of 360 GPa and a Poisson`s ratio of 0.24. The friction coefficients for BIOLOX®delta on BIOLOX®delta in a lubricated condition were chosen as µ = 0.09 based on (Bishop et al., 2013; Haeussler and Pandorf, 2020).

A mesh sensitivity study showed negligible differences (percentage difference for the average absolute deviation of the moment: 0.75% in flexion and 0.88% in extension) between a mesh size of 0.25 and 0.5 mm (between 0.5 mm and 0.125 mm, it is 0.73% in flexion and 0.87% in extension). To reduce simulation time, the optimization was run with the larger mesh size (0.5 mm). As there is typically a hysteresis in the IVD’s moment-rotation curves, that the ex vivo data used in the objective does not contain (as they report only one way, likely neutral to deflection and not back to neutral), the moment was not recorded at the bottom of the TDR but 0.25 mm cranially, so that the TDRs hysteresis would be around the curve from literature. Following the optimization, a simulation of a full flexion/extension-anterior/posterior translation movement cycle, was run with smaller mesh size (0.25 mm) as it led to a smoother moment-rotation curve, and the moment at the bottom of the TDR was reported.

The moment-rotation curves of intact anterior columns were quantified based on the literature (Hartman et al., 2016), in which the kinematics of natural intact specimens were recorded before dissection and replayed. While the study reported rotation angles in relation to an ICR, it did not report the ICR itself. Therefore, the ICR-trajectory from a separate in vivo study (Anderst et al., 2013) was used in the present work. We support combining these datasets by the observation of Jonas et al. (Jonas et al., 2018), who found that ICRs of an in vivo study were comparable to those of an ex vivo cervical study, in which pure moments were applied to mono-segmental specimens. An in silico study (Bayoglu et al., 2019) reported axial compressive loads during movements based on musculoskeletal simulations as a percentage of bodyweight. In the present work, loads applied were scaled to the average bodyweight of an european adult: 70.8 kG (Walpole et al., 2012).

The design objective was to minimize the mean square error (MSE) between the Y-moment-rotation curve of the TDR (Y-axis pointing laterally) and that of C6/C7 anterior columns replaying the kinematics of the intact functional spinal units (FSUs). The ICR-trajectory was used to calculate the displacements (rotational and translational) that were applied to the vertebra superior of the TDR to which also the scaled load curve was applied.

The corresponding optimization problem can be written as (Equations 1a, 1b, 1c): min p  MSE p (1a) s . t .       p ∈ Ω   design space (1b) σ p ≤ σ max   stress constraint (1c) where p is a vector corresponding to the geometrical optimization variables: SphereRadius, AnteriorCylinderLength, PosteriorCylinderLength, AnteriorTorusRadius, PosteriorTorusRadius, Channeldepth (see Figure 1). The design space Ω imposes direct bounds on the geometric variables transcribing constraints regarding manufacturing and anatomy. Additionally, multiple dependent variables and sampling constraints were defined. Finally, σ(p) ≤ σ_max imposes bounds on the maximal principal stresses during the simulations (σ_max = 0.3 GPa, based on the biaxial flexural strength of “extra-high strength” zirconia-toughened alumina ceramics from ISO 6472 Part 2 (2019) and a safety factor of 2) to prevent material failure. The implicit finite element simulations consisted of transient stress initializations followed by simulations of flexion and extension motions coupled with anterior/posterior translation.

Figure 2 shows the optimization setup in LS-OPT V7.0.0 (LSTC, United States). A baseline design and variations of it were created (sampling- and pre-processing stages) and their performance during the simulations (simulation stage) was used to create metamodels (metamodel stage). Based on this, the parameter ranges were narrowed (domain-reduction stage) so that in the next iteration, sampling was done only in the expedient parts of the initial ranges. This way, the design parameters were narrowed down until one of the termination criteria was reached. This optimization strategy is also known as sequential metamodel based optimization with domain reduction. Polynomial metamodels of linear order and D-optimal point selection (maximization of determinant of moment matrix of least squares formulation) of eleven design versions per iteration were selected, as well as the Leapfrog optimizer for constrained minimization algorithm. Termination criteria were set to a design change tolerance of 1e-3, an objective function tolerance of 1e-4 and a maximum of 25 iterations. Reaching any of these criteria was sufficient for termination. The metamodels are generated based on the MSE computed in LS-OPT.

LS-Dyna version R13.0, LS-OPT V7.0.0 and LS-PrePost V4.6.4 (all software from LSTC, United States) were used for the optimization that was run on a computer cluster with 60 cores and 320 GB RAM.

A sensitivity analysis evaluated robustness of the optimized design to uncertainties in the input data. Therefore, the superior-inferior and anterior-posterior trajectories of the ICR were varied within a 95% confidence interval, the friction coefficient and patient-bodyweight were varied between 50% and 150% of the aforementioned values. The maximal MSEs were calculated.

Furthermore, maximal axial rotation and lateral bending based on Hartman et al. (2016) were combined with compressive forces during these movements reported in Bayoglu et al. (2019) and simulated to evaluate whether the stress limits of the ceramic were violated.

It was not possible to apply the vertical load in the universal joint simulator in the same way it was applied in the optimization (application to the vertebra above the implant (optimization) instead of to the centre of the spherical part (experiment)). This affects the moment, as the load application point is deflected more during coupled rotations, leading to a longer lever arm. The load- and motion profiles and application point from the experiment were applied to a FEM (virtual test bed) including the average geometry of the specimens tested, and read out at the position of the loadcell; to compare the measured and simulated responses (Figure 5). The sample geometries were based on geometry measurements using an optical measuring device GOM (ATOS III, Carl Zeiss GOM Metrology GmbH, Germany) prior to any testing.

Load and displacement outputs from the universal joint simulator (mean of cycles 180–189) were used as simulation inputs. This accounts for imperfections of the simulators controls in meeting the desired profiles and smoothing of the profiles for the simulator input. Mesh size (0.0625 mm), initial and maximal timestep (both 1 ms instead of 10 ms) were reduced compared to previous simulations.

Experimental testing of pure flexion/extension under static load was performed on five samples to measure the dynamic friction coefficient of the samples in the tested condition. This friction coefficient was used in the simulation of the experiment.

The optimization reached a termination criterion after 24 iterations, within 5.5 h. The optimized design parameters are given in Table 1, the design is shown in Figure 6. There are clear differences between the parameters of the baseline and optimized design. When performing a coupled flexion/extension-anterior/posterior translation movement, the moment-rotation curve exhibits nonlinearity and hysteresis (Figure 7) as does the moment-rotation curve of an IVD (Newell et al., 2017). Especially in extension, the behaviour of the optimized design matches that of the design objective considerably better than the baseline design (Figure 7). The design optimization achieved a 99% reduction of the design objective (MSE) compared to the baseline design.

In axial rotation and lateral bending, the maximal principal stresses were 0.05 GPa (min: −0.12) and 0.06 GPa (min: −0.10), which is within the bounds of the stress constraints of −0.3 ≤ σ(p) ≤ 0.3 GPa. The sensitivity analysis showed that bodyweight and superior-inferior ICR-trajectory variation have the strongest effect on moment-rotation curves with maximal MSEs of 0.04 and 0.04 N m 2 . For anterior-posterior ICR-trajectory and friction coefficient variation, the maximal MSEs were 0.002 and 1.55e-04 N m 2 . However, a direct comparison of these values is not possible, as for the ICR-parameters, variation was investigated between ±2 SD, and for the other between 50% and 100%. The study from which the bodyweight was taken did not provide SD.

The simulation of the experiment produced moment-rotation curves that are comparable to the experiment (Figure 8). The moment-rotation curves are sensitive to the initial relative positioning of the two implant parts to each other in anterior-posterior direction (Figure 8 - red). Initial placement of the superior implant part 0.45 mm more posterior leads to a better match with the experimental curve and is in line with video data of an initial dry experiment with a test sample when overlaying the video and an animation video of the simulation. The peak moment of the experiment differed by 0.036 Nm (0.92%) in flexion and 0.17 Nm (2.38%) in extension from the simulation of the experiment when taking the 0.45 mm into account (simulation data filtered using moving average filter of 4). However, after the maximal flexion and extension, the simulation had a strong change in the moment without changing the rotation as much as in the experiment. Furthermore, the central parts of the simulation curve are more linear than in the experiment.

Geometric measurements of the samples prior to any testing showed that the main difference between the optimized design and as-manufactured samples was a larger radial clearance of the samples.