Utilizing Computational Modelling to Bridge the Gap between In Vivo and In Vitro Degradation Rates for Mg-xGd Implants

: Magnesium (Mg) and its alloys are promising materials for temporary bone implants due to their mechanical properties and biocompatibility. The most challenging aspect of Mg-based implants involves adapting the degradation rate to the human body, which requires extensive in vitro and in vivo testing. Given that in vivo tests are signiﬁcantly more labour-intensive than in vitro and ethics prohibit direct experiments on animals or humans, attempts are commonly undertaken to infer conclusions on in vivo degradation behavior from in vitro experiments. However, there is a wide gap between these tests, and in vitro testing is often a poor predictor of in vivo outcomes. In the development of biodegradable Mg-based implants, considerable efforts are being made to reduce the overall time and cost of in vitro and in vivo testing. Finding a suitable alternative to predict the degradation of Mg alloys, however, remains challenging. We present computational modelling as a possible alternative to bridge the gap between in vitro and in vivo testing, thus reducing overall cost, duration and number of experiments. However, traditional modelling approaches for complex biodegradable systems are still rather time-consuming and require a clear deﬁnition of the relations between input parameters and the model result. In this study, Kriging surrogate models based on the peridynamic in vitro degradation model were developed to simulate the degradation behavior for two main alloys, Mg-5Gd and Mg-10Gd, for both in vitro and in vivo cases. Using Kriging surrogate models, the simulation parameters were calibrated to the volume loss data from in vitro and in vivo experiments. In vivo degradation of magnesium has one order of magnitude higher apparent diffusion coefﬁcients than in vitro degradation, thus yielding the higher volume loss observed in vivo than in vitro. On the basis of the diffusivity of the Mg 2 + ions modeled under in vitro degradation, Kriging surrogate models were able to simulate the in vivo degradation behavior of Mg-xGd with a ratio between 0.46 and 0.5, indicating that the surrogate-modelling approach is able to bridge the gap between in vitro and in vivo degradation rates for Mg-xGd implants.


Introduction
Magnesium (Mg) and its alloys are under increasing investigation as temporary implant materials due to their non-toxicity, biocompatibility, and biodegradability properties [1,2]. In vitro and in vivo experiments are designed and conducted to study the degradation of Mg-based alloys in both cell culture media and animal models [3][4][5][6]. In the course of those studies, a substantial amount of experimentation is required to investigate the degradation rate of such implants in order to subsequently adapt them to be suitable for the human body. The analysis of several studies published in the field reflect the presence of a large gap between the predictions of in vitro and in vivo studies even for identical degradation systems. The sensitivity of biodegradable systems to environmental conditions and their complexity have also contributed to this gap [5][6][7][8][9][10][11]. A detailed comparison of several published studies was conducted by Sanchez et al. in 2015 [7] to examine whether it is possible to correlate the in vitro and in vivo degradation behavior of Mg and Mg-based alloys. As a result, the study concluded that a clear correlation could not be established. This poor correlation between in vitro and in vivo degradation behavior remains a challenge until today, hence extensive studies are still required for the development of biodegradable Mg-based implants. Unfortunately, due to the expensive and time-consuming nature of such investigations, as well as ethical implications when considering animal experiments, an alternate strategy to comprehensive experimental trials is required. Computational modelling evolves as a powerful instrument for resolving this problem.
Mathematical modelling of biodegradable implant materials is a comparatively novel tool aimed at improving the understanding and prediction of degradation processes. In the respective literature, various mathematical models have been proposed to model certain aspects of the biodegradation process of Mg-based implants [12][13][14][15][16][17]. A more in-depth examination of the various modeling approaches and models of Mg-based implants in in vivo and in vitro environments may be found in recent reviews [18][19][20][21]. Due to the complexity of the system and its multiphysics nature, the available models are often limited and calibrated to specific in vitro or in vivo study systems only. Few comparisons between both environments have been made. For example, the model presented by Sanz-Herrera [22] provided qualitative measurements of the degradation of the Mg screw alloy (Mg-Zn-Zr) in vitro and implanted in a Japanese rabbit. The authors found a four-fold higher corrosion kinetics in vivo than in vitro. However, the modelling of in vivo degradation was limited to a 2D axissymetric geometry.
A peridynamic model was recently developed to simulate the in vitro volume loss of Mg-5Gd and Mg-10Gd bone implant screws (in wt.%) in complex cell culture media [14]. Whereas in classical degradation modelling, e.g., by means of the Finite-Element-Method (FEM), additional effort is required to capture the evolution of the corrosion surface over the duration of simulated immersion [23], the corrosion model obtained from the nonlocal theory of peridynamics proposes a formulation using integro-differential equations, thus resulting in the evolution of the corrosion interface as a natural component of the solution of the governing equations [16,24]. In recent years, these models were extended to address several corrosion phenomena, e.g., passivation, salt layer formation, intergranular corrosion, stress-assisted corrosion, stress corrosion cracking and galvanic corrosion [25][26][27][28][29]. Using a peridynamic formulation, in [14] the degradation was modelled as a diffusive process depending on the diffusion of metal ions from the metal phase across a diffusive degradation layer into the surrounding medium. The model was demonstrated to accurately simulate the volume loss of the Mg-based implants in 3D and in turn served as input for the subsequent FEM model to predict the residual strength of the degraded implant. Considering that the peridynamic degradation model is diffusion-driven, by calibrating the model diffusivity parameters to both in vitro and in vivo data, it may be possible to elucidate correlations between the observed degradation behaviour. However, peridynamic models are relatively complex and computationally expensive. For such computation-intensive models, surrogate modelling has become an increasingly popular approach, introducing mathematical models that are more compact and faster to evaluate without compromising the accuracy of the original models.
Several techniques have been proposed to construct surrogate models using a minimum number of original model simulations [30][31][32][33]. Kriging, also known as Gaussian process modelling, has become a popular surrogate modeling technique in recent years, due to its accuracy and simplified mathematical construction. The quality of Kriging surrogate model predictions, however, is heavily influenced by the size and distribution of the input parameters as well as the original model response used to construct the surrogate model [34][35][36][37]. Recently, Zeller-Plumhoff and Albaraghtheh et al. [12] implemented Kriging to calibrate and optimize the key parameters of a complex physical degradation model of pure Mg under physiological conditions. The calibrated model was able to simulate the degradation with an accuracy of 7%, which is equal to the experimental error. Furthermore, the measured normalized root-mean-square deviation (NRMSE) of the simulated mean elemental weight percentage of precipitated salts reflected an acceptable agreement with the experimental data used in the calibration.
In this study, Kriging surrogate degradation models were developed based on in vitro peridynamic models that simulate the degradation of Mg-5Gd and Mg-10Gd bone implant screws. The optimized surrogate models were extended further to simulate the diffusion behavior of the aforementioned implants in vivo. The key model parameters describing the diffusion of Mg 2+ ions within the electrolyte and the porous degradation layer forming on the metal surface were calibrated to in vitro and in vivo experiments thus providing the possibility to compare the diffusive behaviour with respect to the observed volume loss and degradation rates for both cases and alloys.

Model Calibration Data
The experimental data used for the calibration of the simulations were obtained by evaluating the relative volume loss (VL in %) of two magnesium gadolinium alloys, Mg-5 wt.%Gd and Mg-10 wt.%Gd screw implants (4 mm length, 2 mm diameter, M2 thread with 0.5 × 0.5 mm slotted screw head) in in vitro and in vivo tests. Details regarding the implant manufacturing and testing were published in [8,38]. In summary, in all cases, non-degraded screws were imaged using micro computed tomography (µCT) to obtain their initial geometry. For the in vitro testing the screws were immersed in complex cell culture medium (α-Minimum Essential Medium or Dulbecco's Modified Eagle's Medium supplemented with 10% Fetal Bovine Serum and 1% Penicillin/Streptomycin) for 1, 2, 3, 4 and 8 weeks under cell culture conditions (37°C, 5% CO 2 , 20% O 2 , and 95% rel. humidity). The immersion medium was changed every 2-4 days to avoid saturation of the medium with ions. To obtain the geometry after degradation the degradation layer was removed using chromic acid and the implants were again measured using µCT. For details concerning methods, such as imaging parameters, see [38]. In vivo tests were carried out by implanting the screws into rat tibia with healing durations of 4, 8 and 12 weeks. Again, µCT imaging, in this case synchrotron radiation-based, was used to study the morphology of the degraded screws. The experimental details can be found in [8]. Importantly, the data provided by (synchrotron radiation-based) µCT yields highly resolved information on the screw morphology and can be used to determine the volume loss of the implants [39].
The volume loss [%] for Mg-xGd is calculated as where V(0) and V(t) are the volume of residual metallic material before and after the degradation tests, respectively, as determined from µCT images. The degradation rate (DR) [mm/yr] is calculated based on the volume loss as where A(0) is the initial surface area of the sample and t is the degradation time.

Peridynamic Model and Implementation
In the present study, the peridynamic degradation model from [14] was employed to simulate the degradation behavior of Mg-based implant screws. Given that computational modelling of the electrochemical degradation processes associated with Mg-based alloys is generally very challenging, we adopt a bi-material solid-liquid interface description to efficiently capture the electrochemically driven dissolution and surface evolution. Con-sidering a generic domain Ω ∈ R n d , where n d ≡ 3 is the spatial dimension, the governing equation of the bi-material peridynamic degradation model is given bẏ where C(x, t) is the (molar) concentration of the Mg-based alloy material associated with point x at time t, H is the spherical integration region or neighborhood with radius δ called horizon, p ∈ (−1, ∞) is a scalar value assumed as p ≡ 0 and c(t) is the proportionality coefficient called micro-diffusivity. The latter is a model parameter that depends on the regions (solid or liquid) in which the interacting points x and x are located within the neighbourhood: x and x in liquid part, where χ 2 s and χ 2 l are the diffusivities of the classical local Fickian diffusion equation within the solid and liquid domains, respectively, V L (t) is the volume loss at time t of the Mg-based implant screw and l > 0 is an analytical non-negative model parameter, which represents the depth of the growing degradation layer and therefore leads to a dampening of the degradation with increasing values [40]. The values for l were found to be 1/4.25 and 1/15.53 for Mg-5Gd and Mg-10Gd, respectively. For the numerical implementation of (3) we employ the standard one-point Gaussian integration scheme and an adaptive multigrid discretization method in space and an implicit Euler time-stepping algorithm [14,41]. The in-house high performance computing (HPC) cluster was utilized to perform the simulation on a computational node, that consists of two sockets, each with a 24-core 2.1 GHz Intel Xeon Scalable Platinum 8160 processor, thus providing a shared memory multiprocessing parallelism on 48 cores.

Kriging-Based Surrogate Model and Implementation
The implementation of Kriging techniques to construct a surrogate model for studying the degradation of Mg-based implants is described in detail in an earlier publication of ours [12]. In brief, the peridynamic degradation model is treated as a black-box and the Kriging-based surrogate model is built using a data-driven approach, as illustrated by the schematic diagram in Figure 1. The training data for the Kriging model are sampled from the input distribution in the parametric space of the input parameters. The sampled parameter combinations are evaluated using the expensive computational model (blackbox) and the resulting observations are exploited to build the surrogate model on the whole domain of the input distribution. The Latin hypercube sampling (LHS) method is used to draw thirty different sets of samples based on the assumed parameter distribution of the input parameters (χ 2 s and χ 2 l ) over a user-defined parametric space, similar to the sampling domain shown in Figure 2.
The sampling of the input variables ranges should be representative and capture all possible information about the original peridynamic model [42,43]. Due to the fact that the size and distribution of data points used in the numerical experiments have a major impact on the quality of surrogate model predictions [34,36]. In this case, the parametric space was extended from the range of optimized parameters proposed in [14].
M is the Kriging mapping function, which is the realization of the Gaussian process, that is used to calculate the volume loss Y for any input parameter vector X for a finite period of time. β T f (x) is the trend of the Kriging, which is the mean value of the Gaussian process, σ 2 Z(x, ω) is the realization of the stochastic process that is assumed to have a zero-mean, unit-variance Z(x, ω), where σ is the variance of the process and ω is the underlying probability space, which is defined in terms of the correlation function of the stochastic process R(x, x , θ) [44]. The performance of the Kriging-based surrogate model based on the training simulation data generated from the peridynamic model is estimated by the leave-out-error (ε LOO ), which is defined as where N e is the number of of data points considered during Kriging, M (−i) (x i ) designates that the Kriging surrogate model is obtained using all points of the numerical experimental design except x i and M(X ) is the corresponding surrogate model response of the initial numerical experimental design generated using the original model [44]. The generated Kriging-based surrogate model is calibrated based on the experimental data of both in vitro and in vivo degradation tests of Mg-xGd biodegradable implants. The calibration process was performed by dividing the model domain into equal intervals and then re-sampling each of the intervals using LHS. Intervals with higher uncertainty were further divided and sampled. The goodness of the calibration was measured by the means of the mean absolute error (MAE), which is calculated via with N the number of measurement points, y t the (mean) experimental value at time t for the volume loss at point j,ŷ j the Kriging surrogate model response at point j. The Krigingbased surrogate model is implemented using the Kriging module within the UQLab framework [45], which is installed and integrated under BSD 3-clause license in MATLAB 2021b (The MathWorks, Inc., Natick, MA, USA).

Results and Discussion
Surrogate models were trained to predict the volume loss based on the computational experimental design (ED), which was collected from random samples of the target parameters using LHS. The ranges for the input parameters were set similarly to the ranges assumed during the parameter optimization for the PD models. Thus, the ranges were assumed between 1 × 10 −17 m 2 /s to 1 × 10 −12 m 2 /s and 1 × 10 −12 m 2 /s to 1 × 10 −7 m 2 /s for χ 2 s and χ 2 l , respectively. The developed surrogate models can capture the QoI for both implant alloys, Mg-5Gd and Mg-10Gd, for in in vitro degradation sufficiently well, as can be inferred from the performance of the surrogate model with respect to the peridynamic model shown in Figure 3. It should be noted that the response of the PD model for the in vivo case was assessed based on pure numerical findings of the PD model, the ranges of the numerical test were confirmed based on the experimental data from [38]. The error of the surrogate models, represented by ε LOO was between 2.1 × 10 −3 and 5.5 × 10 −3 for Mg-5Gd and 1.7 × 10 −3 and 3.1 × 10 −3 for Mg-10Gd. Furthermore, the Kriging surrogate model reduces the computational time required to run the training matrix, which contains the training sets of parameters, down to 23.6 s in order to simulate the entire domain over the respective simulated immersion time compared to approximately 3 h for the peridynamic model with the same set of parameters. Generally, surrogate models have a wider range of applications due to their ability for mapping the response domain with respect to the full input distribution. The model key parameters χ 2 s [m 2 /s] and χ 2 l [m 2 /s] were calibrated and optimized using the predictions of the surrogate model and later the model response was validated against VL determined from µCT measurements obtained for in vitro degradation [38] and for in vivo degradation [8]. The optimized values of χ 2 s and χ 2 l for Mg-5Gd and Mg-10Gd for both in vitro and in vivo cases are given in Table 1. The PD model responses were collected over 56 days (8 weeks), so the Kriging models of the in vitro case were initially trained and calibrated over 56 days. Surrogate model responses were numerically extrapolated to 100 days in order to capture the time points from in vivo experiments. The diffusivity of Mg 2+ ions was found to be one order of magnitude higher in the in vivo case with respect to in vitro, which is expected based on the reported degradation behaviour of the current implants by [8,38]. The visualization of the calibrated surrogate models' responses in comparison to the experimental data for 100 days of degradation for Mg-5Gd and Mg-10Gd implants for both in vitro and in vivo are shown in Figure 4a. For the in vitro degradation the simulated VL are in agreement with VL determined from µCT measurements (scatter plots in Figure 4a,b [38], with MAE of 0.03 for Mg-5Gd and 0.08 Mg-10Gd. The same observation is also valid for the simulated in vivo VL, where the simulated values agree with VL determined from µCT measurements [8] (scatter plots in in Figure 4a,b), with MAE values of 0.31 and 0.44 for Mg-5Gd and Mg-10Gd, respectively. The accuracy of the calibration process of the key parameters of the surrogate models was slightly compromised by the uncertainties associated with the experimental data. The relatively high standard deviation of the experimental measurements, as can be seen in Figure 4a, was one of the main challenges of the calibration process. There were also limitations in the number of available data, i.e., the in vitro degradation was limited to five time points between one week and eight weeks, while the in vivo degradation was limited to three time points: four, eight, and twelve weeks. The calibration process would be improved by including more degradation time points, which would also reduce numerical overfitting and underfitting. In addition, it is important to increase the number of samples at each point in order to maintain lower standard errors. The degradation rate calculated based on the volume loss. The VL was determined from µCT measurements as published in [8,38].
In general, based on in vitro surrogate models for the same implants, the current approach of surrogate modeling provided an opportunity to model the in vivo degradation behaviour of Mg-5Gd and Mg-10Gd. It has also been shown that the surrogate models are capable of simulating the degradation of Mg-xGd implants for over 100 days, for both in vitro and in vivo degradation; however, further experiments are needed to validate the results. For comparison the degradation rates (DR) of Mg-5Gd and Mg-10Gd were calculated in both in vitro and in vivo systems, Figure 4b. The surrogate model-based DR were found to follow the experimental data for both the cases.
As the simulation results for the DR suggest, calculating a ratio between DR for in vitro and in vivo as performed by [7] will yield different values, depending on the degradation time. By contrast, using computational modelling, it is possible to derive a time-independent ratio between the DR based on the diffusivity of the Mg 2+ ions. For Mg-5Gd, the mean ratio of χ 2 s and χ 2 l between in vitro and in vivo experiments was calculated as 0.5, while for Mg-10Gd the mean ratio was calculated as 0.46. By establishing such ratios for different material systems, it will be possible to predict the degradation behaviour in vivo based on that observed in vitro in the future, thus reducing the effort required by performing a large number of experiments.

Conclusions
In this study, Kriging surrogate models were developed to simulate degradation behavior as a function of Mg 2+ ion diffusivity. Peridynamic in vitro degradation models were used to generate observations in order to construct Kriging surrogate models. These models were shown to simulate the in vitro degradation behaviour of Mg-5Gd and Mg-10Gd implants with MAE-values of 0.03 and 0.08, respectively. The Kriging surrogate models with optimized key parameters, χ 2 s and χ 2 l , were able to simulate in vivo degradation behavior for Mg-5Gd and Mg-10Gd at MAE values of 0.31 and 0.44, respectively. Additionally, the computational approach employed was able to significantly reduce the computational time and cost as compared to the original peridynamic degradation model. More importantly, the proposed approach provided time-independent ratios based on the variation in Mg 2+ ion diffusivities within the degradation system. Using the timeindependent ratios established between in vitro and in vivo experiments by means of Kriging surrogate models, we may in the future be able to infer in vivo degradation rates based on in vitro experiments.

Data Availability Statement:
The experimental data for the validation can be found in [8,38]. The data sets generated and analyzed during the current study are available from the authors on reasonable request.