Chapter 10 Model applications: patient-specific porous acetabular implants

Two patient-specific finite elements (FE) models of the hip joint, i.e., one with implant and one without, were derived from computed tomography (CT) scans of the pelvis, implemented using MARC/MENTAT software (Perez-Boerema et al., in preparation). The elastic modulus of the implant was parametrized, making use of two paraboloids, one for the upper half of the acetabular implant and one for the lower half. The FE model of the hip joint with implant was then used to generate a surrogate model by first applying a principal component analysis (PCA) in order to reduce the size of the output space and subsequently constructing a radial basis function (RBF) surrogate, implemented in MATLAB. The resulting surrogate model was used in a multi-start local optimization routine that relied on sequential quadratic programming to perform the local optimizations. The optimization had as objective minimizing the difference in stresses between the physiological control (no implant) and the case with implant. The coefficients of the paraboloids used to parametrize the elastic modulus were taken to be the design variables. The elastic modulus was limited to values between 1.2 and 110 GPa, while implant stresses were constrained to at most 0.3 times the yield stress. The optimized stiffness distribution was subsequently used to design an acetabular implant with a local porosity chosen to match the optimized global stiffness distribution. To this end, the open-source tool ASLI [207] was developed in-house, facilitating the generation of a skeletal-gyroid-based infill of varying porosity that can be tuned to match a predefined (in this case, the optimized) stiffness.

In this application chapter, we apply the previously developed quantitative mathematical model capturing the biodegradation of magnesium-based implants (Chapter 3). Magnesium (Mg) was selected to start with due to its acceptable mechanical properties, biocompatibility, and contribution to osteoinductivity [5]. The developed model captures the release of Mg ions, the formation of a protective film, the dissolution of this film due to the presence of various ions in the surrounding environment, and the change in pH. This was accomplished by deriving a system of nonlinear time-dependent reaction-diffusion partial differential equations (PDEs) from the underlying oxidation-reduction reactions occurring during the biodegradation process of metallic materials in simulated body fluid. The level-set formalism was employed to track the movement of the biodegradation front between the implant and its surroundings. The equations were solved implicitly using the finite element method for spatial terms (with a 1st order Lagrange polynomial as the shape function) and backward-Euler finite difference method for temporal terms on an Eulerian mesh, implemented in the open-source PDE solver FreeFEM [112]. The details of this model can be found in Chapter 3 and [31]. Two separate simulations were performed using two different diffusion coefficients, as calibrated in our previous work, to model the degradation behavior in non-buffered and buffered solutions.

In order to build the computational model, the resulting surface mesh of the surrogate-based optimization routine, consisting of 5347924 faces, was converted to a volume mesh and embedded in a cubic container that was to act as the environment during the biodegradation simulations. The conversion of the surface mesh was performed using GMSH [87], and the embedding process was carried out using the internal mesh engine of FreeFEM, called BAMG. The mesh was refined on the implant-environment interface to increase the numerical accuracy of the interface tracking model, leading to a final mesh with 45870053 elements. The mesh refinement process was handled using Mmg [63] (and ParMmg [24]). Fig. 10.3 shows the embedded cup model and a cross-section of the final generated mesh, which was refined on the metal-solution interface.

For dealing with a problem of this size and making the model capable of being scaled in large-scale computing systems, the model implementation made use of the HPC techniques available in the FreeFEM language v4.10 and PETSc toolkit v3.16.1 [25]. In this implementation, METIS and ParMETIS graph partitioners [144] were used to decompose the mesh into various partitions, and then the partitioned mesh was distributed among available computing resources using the HPDDM preconditioner [140]. Fig. 10.4 shows the decomposition of the computational mesh. Moreover, the HYPRE BoomerAMG [78] preconditioner and GMRES iterative solver [223] were used to solve the linear systems of equations obtained from applying the finite element discretization on the model. More details of the implementation can be found in Chapter 7 and [30].

The simulation was carried out using 2,000 CPU cores with 16.5 TB of available memory on the Dutch national supercomputer, Snellius. The simulation results, comprising of 95000 files with a total size of 148 GB, were visualized using a parallel client-server remote rendering approach in ParaView server v5.9.1 running on 128 CPU cores on the ARCHER2 supercomputer.

In order to obtain the scaling behavior of the model in an HPC environment, strong and weak scaling tests were performed on the Snellius supercomputer. For the weak scaling test, in addition to the mesh with 45870053 elements, two more mesh files consisting of fewer elements were generated using the aforementioned procedure for embedding and refining the mesh. These models had 15989521 and 29035491 tetrahedral elements, respectively, to which the number of employed CPU cores was adjusted accordingly. The strong scaling test was carried out for all three model sizes by varying the number of employed CPU cores from 60 to 9000.

The optimized design to minimize stress shielding, depicted in Fig. 10.5, has a skeletal-gyroid infill of varying porosity with apparent elastic moduli ranging from 4.7 to 49 GPa and an estimated acetabular stress shielding reduction of about 56% compared to a solid titanium implant.

The varying porosity was adjusted during a surrogate-based optimization process to match the stiffness distribution required to minimize stress shielding during the bone healing process. Considering the target application is bone tissue engineering, the infill used to generate the lattice structure was selected based on triply periodic minimal surfaces (TPMS). TPMS-based lattices have shown good performance for biological processes of cell attachment, migration, and proliferation [212, 113], making them a suitable choice for tissue engineering applications where their aim is to guide tissue growth. This positive influence is routed in the appropriate balance of various structural and mechanical properties required to improve bone tissue formation, such as yield strength, fatigue strength, adequate mechanical stimulus, and permeability [201, 121].

Fig. 10.6 shows the mass loss during the biodegradation for simulations performed in buffered (low degradation rate) and non-buffered (high degradation) solutions. Mass loss is one of the commonly-used indicators for the degradation rate, demonstrating that the biodegradation rate is much higher in the non-buffered solution. Although evaluated on different geometries with certain effect on the results, loss of material over time was found to be comparatively in line with the values obtained in Chapter 3 and [31], with the mass loss in the saline solution being 6 to 10 times more than that of buffered solution in evaluated time points before 21 hours. It shows that scaling the model in an HPC environment does not affect the quantitative predictions made by the model. The developed mechanistic model of the biodegradation process includes a level-set equation correlating the loss of material to the velocity at which the implant interface shrinks. The aforementioned agreement of results on a large scale (a model with $\sim$45M elements) shows that the interface tracking formulation behaves efficiently even in problems with a high level of details. This verifies the performance of the developed biodegradation model, which has never been tested before in such high resolutions.

From a qualitative point of view, visualizing the biodegradation results over time, depicted in Fig. 10.7, shows that the acetabular implant degrades faster in the regions with higher porosity, i.e., the regions with more exposed surface area to the environment. These are the regions designed to have lower stiffness, resulting in high porosity after being filled with TPMS unit cells with a lower volume fraction. From the induced bone formation perspective, the implant parts with low stiffness disappear during the bone healing process, a fact that can be further taken into account in the optimization procedure.

Fig. 10.8 demonstrates a similar visualization but with a zoomed view on the surface of the implant being plotted along with a cross-section of the medium showing the concentration of released Mg ions. The movement of the corrosion interface is formulated based on the release of material ions, and as a result, the material loss rate is higher in regions with higher ions concentration. The concentration is directly correlated to the exposed surface area, meaning that a higher surface-to-volume ratio results in a higher material release. This explains why the regions with higher porosity degrade faster in the current biodegradation model.

Regarding the clinical interpretation of the obtained results, it is worth reminding that the simulated cup is a simplified version of the acetabular component of the implant. In real-case scenarios, this component contains a patient-specific part on the back that can be manufactured from biodegradable materials (shown in Fig. 1). In order to simplify the geometry construction and the optimization procedure, the shape of the patient-specific component was neglected and was assumed to be the cup instead (Fig. 10.2), which should be created from non-degradable materials in real applications to provide enough support for fixing the implant behind the polyethylene liner. As a result, the simulation results cannot be used directly for clinical interpretation. However, by only changing the geometry construction from the CT scans performed for the surrogate-based optimization model, the whole workflow can be applied for the patient-specific design of acetabular implants in more clinical-relevant cases.

The results of the strong scaling tests are plotted in Fig. 10.9, which shows the solution time of individual equations of the biodegradation model in a single time step versus the varying number of CPU cores. The comprising equations are the transport equation for Mg ions (Mg equation), hydroxide ions (OH equation), and chloride ions (Cl equation), as well as the film formation (film equation) and the derived level-set equation. These equations are detailed in Chapter 3 and [31]. The run time of a single time step was measured for three model sizes of $\sim$16M (small), $\sim$29M (medium), and $\sim$45M (big) elements.

As can be seen in the strong scaling results, different PDEs show different scaling behavior, the most distinct of which belongs to the OH equation used to calculate the pH change in the surrounding environment. The different scaling behavior is rooted in the boundary conditions used for each equation. They impact the formation of the linear system of equations and appear in the scaling results due to the penalization technique used to implement boundary conditions.

Fig. 10.10 shows the results of the strong and weak scaling tests side by side. As can be seen in the weak scaling result, the model shows an ideal scaling behavior in environments with less than 1000 cores. The system shows non-ideal behavior with more CPU cores than this, as can be observed in strong scaling results, in which adding more CPU cores does not cause a linear drop in the execution time (wall time). The parallel efficiency graphs plotted in Fig. 10.11 show this behavior more clearly, which are calculated from the strong scaling results.

The strong scaling test and parallel efficiency plot show that the maximum size for the distributed computing environment is 2000 CPU cores, above which no significant improvement in the execution time was observed in all three tested problem sizes, and the efficiency drops below 50%. It was expected to observe different scaling behavior for models depending on the size of the problem as in problems with smaller sizes, the costs associated with inter-nodes communication become effective faster in comparison to bigger models. However, no significant difference was observed here, and all three problem sizes show very similar behavior in the parallel efficiency plots (Fig. 10.11), demonstrating the non-ideal behavior of the developed model in terms of efficient usage of available computing resources.

Although the executed performance analysis showed poor parallel efficiency in higher number of CPU cores, it was necessary to reveal the weaknesses of the developed model in this regard. The weak and strong scaling results and the parallel efficiency graphs clearly indicate that the code can be improved from the performance point of view. Given the current state of HPC for biomedical-related simulations and computational modeling works, even taking 10K CPU cores is not unrealistic. The first part of improvement seems to be related to assembling the linear system, in which the schemes used for numerical integration, the order of elements used for various involved physics, and the method for applying boundary conditions play an essential role. Secondly, the configuration of employed preconditioners and iterative solvers can be modified for individual components according to the numerical computing requirements, which can affect their scaling behavior (Fig. 10.9). Thirdly, the load balancing issues caused by the mathematical problem formulation, in which the system has nothing to solve in regions still in the solid block, seem to contribute to the poor parallel performance. This can be improved by revising the mathematical and computational problem definition from the load-balancing perspective. However, making a more tangible conclusion on the performance bottlenecks of the developed model requires further analysis of the results for the parallel speedup, similar to the analysis performed in Chapter 7 and [30]. Furthermore, it should be taken into account that such performance analysis has an inherent simplification of more complex technical aspects affecting the simulation execution time, such as load balancing, network communication, and per-node limitations in memory access, which may differ in each HPC environment due to the configurations made by the system maintainers.