Computational modelling of magnesium degradation in simulated body flui under physiological conditions

Berit Zeller-Plumhoff, Tmdur AlBrghtheh, Dniel H?che, Regine Willumeit-R?mer

aHelmholtz-Zentrum hereon GmbH, Institute of Metallic Biomaterials, Max-Planck-Stra?e 1, Geesthacht, 21502, Germany

b Helmholtz-Zentrum hereon GmbH, Institute of Surface Science, Max-Planck-Stra?e 1, Geesthacht, 21502, Germany

Abstract Magnesium alloys are highly attractive for the use as temporary implant materials, due to their high biocompatibility and biodegradability.However,the prediction of the degradation rate of the implants is difficult therefore,a large number of experiments are required.Computational modelling can aid in enabling the predictability, if sufficientl accurate models can be established.This work presents a generalized model of the degradation of pure magnesium in simulated body flui over the course of 28 days considering uncertainty aspects.The model includes the computation of the metallic material thinning and is calibrated using the mean degradation depth of several experimental datasets simultaneously.Additionally, the formation and precipitation of relevant degradation products on the sample surface is modelled, based on the ionic composition of simulated body fluid The computed mean degradation depth is in good agreement with the experimental data(NRMSE=0.07).However, the quality of the depth profil curves of the determined elemental weight percentage of the degradation products differs between elements (such as NRMSE=0.40 for phosphorus vs. NRMSE=1.03 for magnesium).This indicates that the implementation of precipitate formation may need further developments.The sensitivity analysis showed that the model parameters are correlated and which is related to the complexity and the high computational costs of the model.Overall, the model provides a correlating fi to the experimental data of pure Mg samples of different geometries degrading in simulated body flui with reliable error estimation.

? 2021 Chongqing University.Publishing services provided by Elsevier B.V.on behalf of KeAi Communications Co.Ltd.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/)

Peer review under responsibility of Chongqing University

Keywords: Biodegradation; Magnesium; Computational modelling; Corrosion; Uncertainty quantification Kriging.

1.Introduction

Magnesium (Mg) and its alloys are emerging as promising alternatives to permanent implant materials, due to their biodegradability, low toxicity and ability to facilitate bone growth.There are numerous studies investigating the degradation of Mg-alloysin vitroin different cell culture media andin vivoin animal models.Due to the high complexity and sensitivity of the system a gap between the results ofin vitroandin vivostudies remains [1].Computational modelling can be used to facilitate closing the gap and enable predictions based onin vitroexperiments, and thus speed up the development process of new implants if reliable models are developed.

Different kinds of computational models of Mg degradation have been published until now; including physical, phenomenological, analytical and machine learning-based models.A simplifie analytical model based on the mass conservation law was presented by Dahms et al.[2].Chemical and electrochemical processes were not included in this model due to the analytical nature of the approach.Nevertheless, by capturing both long-term(degradation propagation)and shortterm(initiation)characteristics,the degradation rate calculated by the analytical model was found to be in agreement with experimental data for up to 18 days.In other studies, continuum damage (CD) models are used to describe the degradation of Mg phenomenologically.Gastaldi et al.[3]and Grogan et al.[4]used CD models to simulate the degradation of Mg-based stents.The CD approach enabled studying the micro-scale geometric discontinuities induced by the degradation on the mechanical integrity of the overall specimen without the need for the progression of the degradation to be modelled explicitly.Therein lies also the disadvantage of the CD approach,as it fails to incorporate the complex electrochemical processes during degradation, which change (locally) e.g.due to the composition of the degradation medium.Thus, a CD model will require calibration for each change in the corrosive environment or material composition.

Physical models based on the Nerst-Planck equation for Mg exposed to aqueous electrolytes were developed for example by Jia et al.[5], Deshpande [6]and Grogan et al.[7].In these models, the geometric change of the computational domain due to the degradation was incorporated by applying an adaptive mesh.Grogan et al.[7]et al.applied this model to a stent geometry and additionally took the formation of a stable degradation layer into account.By doing so,they could sustain the assumption of a transport-controlled degradation across this layer.Sun et al.[8]et al.additionally introduced the time-dependent formation and deposition of corrosion products on the metal surface by including the changes in both electrical potential and mass transfer.Bajger et al.[9]et al.developed the model presented in [7]further by assuming the diffusion as the rate-limiting process and the absence of a protective layer.This model assumed a zero concentration of Mg2+ions and only the presence of Clions in the electrolyte at the beginning of the reaction.The moving interface of the degrading metal Mg specimen was modelled implicitly using the level-set method.The model showed good agreement in terms of hydrogen evolution with published datasets.The physical model presented by Sanz-Herrera et al.[10]et al.incorporated the highest number of chemical species involved in the degradation process so far.This model included a set of seven species in the model;(metallic)Mg,H2,H2O,Cl-,Mg(OH)2and MgCl2,where the concentration gradient of each of these species was modelled using reaction-diffusion equations.This model was reported to be valid for slow (0.05 mm/year) and fast (0.3 mm/year)degradation rates, but mostly qualitative comparisons were provided that were within the order of magnitude of experimental values in the literature.Ahmed et al.[11]et al.presented a simplifie model, which was able to demonstrate the effect of different chemical species on thein vitrodegradation of Mg.This model was simplifie by using the dimensionless presentation of different model parameters,which reduced the mathematical complexity of the model.

A main drawback of most physical models is that they consider only the formation of magnesium hydroxide(Mg(OH)2),which is however not the main degradation product forming during the degradation of Mg in simulated body fluid(SBFs) under physiological conditions [12,13].Instead, the formation of a range of other precipitates must be taken into account.This is pivotal, as the forming precipitates and the evolving pH and H2gas influenc the environment around the implant.Thus, if more complex models should be developed that include the cell-material interaction or any mechanical behaviour, it is important that both the degradation rate of Mg implants is modelled correctly, as well as the formation of the degradation layer that forms in SBF.Additionally, the published computational models are lacking three-dimensional (3D) input data on the degradation layer porosity as a key transport parameter, which limits their validity, as in many cases the degradation of Mg is modelled as a diffusion-limited problem.

A recent model by H?che [14]investigated the microgalvanic degradation in a system of Mg coupled to more noble materials.In addition to the formation of magnesiumhydroxide, different chemical reactions occuring within the immersion medium were included in this model,such as water dissociation.This model considers deposit layer porosity and tortousity and was able to simulate and quantitatively predict the effect of degradation products in a time-dependent manner.In this work, we are presenting a physical model that is an extension of the model presented in [14].Specificall , we are modelling the degradation of pure Mg in SBF over four weeks in a semi-static immersion test.SBF was selected because it mimics the ionic composition of blood plasma more closely than NaCl.At the same time SBF does not contain organic components, which would increase the complexity of the ongoing interactions significantl.We assume a homogeneous degradation during which precipitates form at the Mg interface that do not re-dissolve.These precipitates form a porous degradation layer through which the chemical species diffuse.The degradation layer is assumed non-reactive and forming H2gas bubbles don't interact with any model components.

The model is calibrated using data from weight loss measurements and volume loss measurements based on 3D micro computed tomography (μCT) data, as well as data from energy-dispersive X-ray spectroscopy (EDX) on the formation of certain precipitates.Moreover, experimental data from transmission X-ray microscopy was used to specify the porosity of the degradation layer.The estimation and calibration of different model parameters are one of the main challenges in the model development process due to the complex multiphysics nature of the current model and the high computational costs associated with it.Furthermore, most key parameters have no reported values in literature.For example,there are no reported values for any of the reaction rates of the system [9].Uncertainty quantificatio (UQ) provides the methodologies to characterize the output uncertainties of a mathematical model based on variable inputs and to identify the different sources of the uncertainties within the model itself.For the parameter estimation we are employing the UQ approach by using a Kriging-based metamodeling (or surrogate modelling) methodology.Metamodels are constructed by generating a response of the model using a relatively small set of parameters and simulations.This reduces the computational costs and accelerates the calculation processes, but also allows for further advanced analyses, such as sensitivity analysis, reliability analysis and Bayesian model calibration,therefore enabling a better understanding of the system [15].To quantify the influenc of model parameters over the model outputs, a global sensitivity analysis was performed.

2.Methodology

2.1.Modelling of electrochemistry and chemistry

We are modelling the degradation of pure Mg in SBF in a semi-static immersion test.The Mg specimen is assumed to be of a purity between 99.92% and 99.95%, but the influenc of impurities is not directly modelled.The formulation of SBF that is modelled is taken from Bohner and Lemaitre [16]and contains the following chemical species: Na+, Ca2+, HCO-3,HPO24-and Cl-, in addition to H+and OH-.Table 1 shows the different ionic concentrations in the medium.

Table 1Concentrations of ionic components in simulated body flui [16].

Table 2Dissolution and dissociation constants as presented in [13].

We are considering the following (electro-)chemical reactions taking place at the implant surface, as suggested in [13].Firstly, the magnesium is dissolving leading to the formation of hydrogen.

In the meantime, a number of precipitates form based on the ions contained in the degradation medium.Additionally, we assume that water, hydrogen carbonate and hydrogen phosphate are dissociating in the media according to [13]:

While the oxygen reduction reaction was shown to be of importance during the degradation of pure Mg in NaCl, it is not incorporated into the model as its influenc diminished after degradation over less than an hour for commercial purity Mg [17].We therefore assume that it has little influenc on the degradation process at the time scales considered in the present model.

The precipitation and dissolution constants (pK=-log10K)at 37°C as presented in[13]and adjusted for the above chemical reactions are shown in Table 2.

2.2.Numerical model

A rectangular shape is used approximate a pure Mg sample that is covered with a degradation layer and degradation medium, see Fig.1.An initial porous degradation layer of 10 nm is assumed to be on the sample surface.This equates to the thin MgO layer that forms on Mg samples immediately after exposure to air following the manufacturing of the sample [18].In the following,ΩMgdenotes the metal domain,ΩDLthe porous degradation layer of porosity?andΩLthe aqueous immersion medium.Fig 1 shows a schematic of the geometry, the included domains and modelled (electro-)chemical processes.

2.2.1.Transport of chemical species

The transport of the chemical speciesiof concentrationciis described as

The firs part of the species flu describes the diffusion of the species within the domain, whilst the second part describesthe migration within the electric fieldzidenotes the species charge (e.g., for Mg2+zMg=2),Vis the electric potential,R≈8.314[J / (mol K)]the gas constant andTthe temperature.F≈96,485[C/mol]is the Faraday constant.The diffusion coefficien is given as

Fig.1.(a)Schematic of the degradation process of pure Mg in simulated body fluid The metallic Mg sample is assumed to be covered by a porous degradation layer from the start, with an initial thickness of 10 nm.The ions diffuse through the liquid medium and the liquid phase of the degradation layer towards the metal interface, where the Mg dissolution is taking place.Precipitates form within the degradation layer and the dissociation of different ionic species occurs in the liquid phase.Not to scale.(b) Slice through μCT data of pure Mg disc degraded in SBF after 3 weeks for comparison.

whereτis the tortuosity of the porous layer,?the porosity andDl,ithe diffusion coefficien of the chemical species in the liquid phase of the porous medium.For simplicity, it is assumed that the porosity can be described by?only, i.e.τ=1.?=1 inΩLtherefore,

Furthermore,

The reaction rateRiof chemical speciesiis described by

withjdenoting the chemical species with whichireacts under equilibrium constantKeq,ljand backward reaction constantklj, i.e.in reactionljof Eqs (3)-(9).However, the formation of precipitates (Eqs.(3) to (9)) is assumed to take place only in the porous degradation layer.

Initial valuesc0,iof ionic speciesiin the liquid phase are determined from the medium composition and given in Table 1.c0,i=0 in the metal domainΩMg.

A constant infl w is define at boundary?ΩL,in, i.e.ci=c0,ito model and approximate the semi-static immersion test.Symmetry is assumed at?Ωsym.For the inner boundary betweenΩMgandΩDL, denoted as?ΩMg, flu discontinuity is assumed, as well as an electrode surface coupling, i.e.,

2.2.2.Formation of precipitates on mg surface

Precipitates are forming on the metal surface depending on the pH as highlighted in Fig.2, according to Eqs.(3)-(9).The formation of precipitateplis given as

with the equilibrium constantKeq,land backward reaction constantklfor reactionlof the Eqs.(3)-(9) andiandjdenoting the two species reacting in the respective reaction.To determine if all chemical reactions presented in Section 2.1 have to be taken into account in the computational model, the freeware Hydra-Medusa [19]is used to determine the fractions of the above precipitates forming in SBF at 37°C depending on the pH.To this end, the solubility constants used in the software were adapted using the constants in Table 2.Fig.2 shows the respective calculated fractions.Within the experimentally observed pH range, marked with vertical lines,neither nesquehonite(Mg(HCO3)OH·2H2O(s))nor portlandite (Ca(OH)2(s)) are predicted form.However,nesquehonite has been observed within the degradation layer experimentally using X-ray diffraction in a previous study[12].By contrast, the formation of portlandite hasn't been reported.Based upon this information, only Eq.(7) is removed from the model.

Fig.2.Fractions of precipitates forming in SBF at 37°C depending on pH as computed from [19].The figur is taken and adapted from [12]under a CC BY 4.0 license.Vertical lines indicate experimentally measured minimum and maximum pH values.The fractions are with respect to the pecipitates based on Mg and Ca, respectively.

2.2.3.Magnesium degradation

As previous experiments have shown there is no easily quantifiabl correlation between the electrochemical driven anodic dissolution related degradation of Mg and the ionic composition of the degradation medium [12].Thus, we are mathematically uncoupling the formation of precipitates and the degradation of Mg.As the degradation behaviour shown in [20,21]followed a linear trend following an initial period of low degradation the degradation of metallic Mg (Eq.(1))is modelled as:

which equates to the effective rate of Mg2+ion release at?ΩMg,tinitdenotes the starting time of degradation.

2.3.Experimental data for model input and calibration

Many model parameters are taken from the literature, such as the precipitate dissolution constants in Table 2.The initial concentration of the ionic species and their concentration at the boundaries is given by the SBF solution itself as stated in Table 1.Table 3 contains all other parameters involving the ionic species and forming precipitates.

Table 3Model parameters as taken from the literature.

For simplicity we will assume in the beginning that the current density in the electrolyteilis constant as measured experimentally for the corrosion potential.Ecorrof Mg in SBF wasdetermined with respect to SHE as -1.66 V/cm2with a correspondingilof 31.6·10-6A/cm2[22].The backward reaction constants of Eqs.(10)-(12) are all set to 1.4·104m3mol-1s-1for simplicity.

2.3.1.Degradation layer porosity

The porosity of the degradation layer evolving in SBF is determined using transmission X-ray microscopy, as published in [20].In brief, pure magnesium (99.92%) discs were degraded in SBF in a semi-static immersion test for 1-4 weeks.Following degradation, the samples were dried and imaged using μCT to determine the overall degradation morphology.Small samples were then extracted from the degradation layer using focussed-ion-beam (FIB) milling.These were imaged tomographically using TXM at a resolution down to 37 nm.The porosity of the degradation layer was determined based on the 3D image data.As the samples taken from 1 and 2 week specimen were imaged at higher resolution, we use the mean calculated porosity of 18.67% for these time points as input into the model.Fig.3 gives an overview of the experimental data used in the current model for input and data calibration.

2.3.2.Model calibration data

The calibration of the model is performed by evaluating the mean degradation depth of the Mg specimen using both data from μCT imaging and weight loss data.Data for 1-4 weeks is taken from [20],while additional μCT data for 6-10 days is taken from [12]and from newly measured data.For a description of the data acquisition for the previously unpublished data please see the appendix.Moreover, the data from[20]and [12]differ, as discs of 99.92% purity were used and wires of 99.95% purity in the experiments, respectively.By calibrating the model to the data based on both geometries and materials simultaneously, it is possible to assess the generality of the model with respect to metal impurity within a certain range.

Fig.3.Overview of the datafl w in the current model.The experimentally determiend porosity of the degradation layer is given as input to the computational degradation model.The model simulates the Mg metal thinning and the formation of precipitates on the metal surface in an uncoupled manner.Based on the determined concentration changes, the Kriging meta-model determined the mean degradation depth and elemental wt% of the elements in the degradation layer, based on the different parameter settings.The results are compared with the experimental values on the basis of the NRMSE,which is minimized.

To calibrate the model for the formation of precipitates in the degradation layer, calibration data from EDX measurements is taken from [20]and [12].In both cases, the mean elemental weight percentage (wt%) of each element in the degradation layer is calculated and used for the calibration.Again, by calibrating to experimental data from two experiments the model generality can be assessed with respect to the deviation of experimental results.

2.4.Implementation

The computational model is implemented in Comsol Multiphysics 5.6 (COMSOL AB, Stockholm, Sweden).The model geometry is a quasi-1D geometry.A schematic of the model geometry is given in Fig.1, however, the Mg domain measures 0.2 mm in length, the initial degradation layer is 10 nm thick and the electrolyte domain is 0.2 mm long.The geometry is 1 μm high.The mesh is based on rectangular finit elements with a fi ed height of 1 μm and a variable thickness, to refin the mesh near the degradation layer.Fig.4 shows the extents of the initial mesh.

The arbitry Lagrangian-Euler (ALE) method is used to capture the dynamic deformation of the model geometry and its moving boundaries.The ALE includes two frames:the reference frame, which has a fi ed coordinates frame, and the spatial frame, which has coordinates moving with time, subject to boundary conditions[6,7].Thus,utilizing ALE enables explicitly tracking the moving interface of the degrading metal and it produces a smooth deformation of the mesh subjected to the constraints placed on the boundaries.

The normal velocity of the metal interface is calculated as

from which the mean degradation depth(MDD)is determined as a function of the degradation rate (DR) [29]:

Note that no displacement is considered for all other boundaries.

The mean elemental weight percentage of each element is calculated based on the concentration of each precipitated speciescprecjand the molar mass.For example, for elementi, we initially calculate the overall mass concentration in the degradation layer as:

withMelemithe molar mass of the element, andNprecj,elemithe stochiometric factor for elementiin precipitatej.The elemental weight percentage of the same elementiis then

for elements C, Ca, P, H, O and Mg.H cannot be measured experimentally, therefore no calibration data is available for its wt%.The presence of H needs to be taken into account to calculate the mean element wt% of all other elements though.

2.4.1.Parameter estimation under uncertainty

Traditionally,the optimal model parameters are determined by running the computationally expensive simulation for a large number of possible parameters, followed by an error estimation.By employing metamodelling, the computational cost can be reduced while still covering the same range of possible parameter values.In this case, Kriging is used to quantify and optimize the reaction ratesklforl=3,...,9,kdegandtinit, by obtaining the optimal reaction rate constants to calibrate the model output with respect to the experimental data.In Kriging, the response of the model is constructed as the summation of the mean structure and a zero-mean stationary Gaussian stochastic process [15,30].

for anyx=(kdeg,tinit,k3,k4,k6,k7,k8,k9).The firs term of the right-hand side of the Eq.(24) is the mean value of the Gaussian process, which is also known as the trend,wherePis the number of the random arbitrary functionsf fj;j=1,...,Pandβis the corresponding coefficientβj;j=1,...,P, which is obtained from the generalized least-squares method.The second term of the right-hand side of Eq.(24) is a realization of the stochastic process that is assumed to have a zero-mean,unit-varianceZ(x,ω),whereσisthe variance of the process andωis the underlying probability space, which is define in terms of the correlation function of the stochastic processR(x,x′,θ)[15,31].

Fig.4.Initial mesh used to mdel the magnesium corrosion.The overall extents of the mesh are denoted in the graph and the different regions are denoted in reference to Fig.1 with several stages of zoom to enable the visualisaton of the mesh within the initial 10 nm thick degradation layer..

The mathematical presentation of the Kriging model to calculate the degradation rate constant would be

and for the mean elemental weight percentage:

for elementsibeing C, Ca, P, O and Mg andjinkjreferring to the number of the precipitation reaction in Eqs.(3)-(9).

The Kriging algorithm was implemented in MATLAB(R2020b,The MathWorks, Inc., Massachusetts, United States)using the Kriging module of UQLab [32], a MATLAB-based Uncertainty Quantificatio framework.A wrapper script is used to convert data-format into a compatible format between COMSOL (viaLiveLink) (true model output) and UQLab(data needed to create the input module).Later the output data of the true model was used to create the experimental design, which represent the initial known function values for the objective function in the Kriging metamodel.The observations of the true model are obtained by running the initial input samples in COMSOL through LiveLink and report the model response.The MDD and the elemental wt% of each chemical element in the true model are calculated then by the wrapper script.To describe the similarity between observationsxand the new pointsx′generated by the Kriging metamodel, the Mate'r 5/2 correlation function was used, which is given by [33]

where the hyper-parameterθwas determined using the Maximum-likelihood (ML) estimation method.

Initial tests are performed to specify the testing ranges for each parameter.These tests are performed using the rate equation and guided by the experimental data (see Section 2.3.2).The maximum and minimum limits of the test intervals of each parameter were provided to the input module.Within the input module of the Kriging model, the Latin hypercube sampling (LHS) strategy is used to draw samples within the user-define range.Random test sets within the input intervals were created using the LHS sampling strategy and tested based on the convergence of the model, the fina input ranges are given in Table 4.Initially, 80 different sets of each tested parameters were generated using the random sampling strategy of LHS.These sets were randomly created under the assumption of a uniform distribution.After generating the test samples, a random selection of 25 sets out of the 80 sets was selected and runviaLiveLink COMSOL to generate the true model response, which is needed to create the Kriging metamodel.In the current case, the true model object was directly define as the Comsol model (viaLiveLink), hence the experimental design test-matrices will be automatically de-fine within the UQLab through the wrapper function based on both model object and the input module.

Table 4Tested ranges for the model parameters and optimal parameters as estimated by Kriging.LOO-error for each parameter is given.As multi-output Kriging was performed for k3-k9, one LOO-error is given for these.

The performance of the Kriging metamodel was estimated by the leave-one-out error (LOO-error) given by

withNKthe number of data points considered during Kriging and whereM(-i)(xi)designates that the Kriging metamodel is obtained using all points of the experimental design exceptxiandM(X)is the corresponding metamodel response of the initial experimental design [15,31].

The optimized value of each parameter based on the Kriging metamodel was found by the means of the normalized root mean squared error (NRMSE):

withNthe number of measurement points,ytthe (mean)experimental value at timetfor outputiat pointj, i.e.either MDD or wt%, ?yjthe Kriging metamodel response at pointj,andymaxandyminthe maximum and minimum experimental values.

2.5.Sensitivity analysis of model parameters

In order to quantify the relative importance of each input parameter on the variability of the model output a global sensitivity analysis was performed.Global sensitivity analysis considers the inputs of the mathematical model as random variables and computes the relative importance of each parameter over the model output.In the metamodel we assumed that the model parameters were uniformly distributed.If all the input parameters for the sensitivity analysis are assumed independent and the parameter vector is given byXi=(kdeg,tinit,k3,k4,k6,k7,k8,k9)ifori∈{1...,P}, wherePis the number of metamodel evaluations.Xuis a subset of metamodel inputs, i.e.Xu=(Xi1,Xi2,...,Xik)foru={i1,i2,...,ik}?{1,2,...,P} with corresponding outputsM(Xu).The metamodel resultsM(X)can decomposed into contributions from each subset of input parameters according to the variance decomposition analysis(ANOVA), with recursively define contributions [34,35]

The Sobol' index for the corresponding subsetuof metamodel runs was calculated as the variance of the subset over the total variance:

Consequently, the total Sobol' index for the entire set was calculated as:

In order to reduce the computational cost of the system, the sensitivity analysis was performed using the Kriging metamodels.

3.Results and discussion

3.1.Metamodel response

The visualization of the Kriging metamodel response at 28 days for different degradation rate constants(kdeg)is shown in Fig.5.The y-axis is the MDD estimated by Kriging and the x-axis represents the different values of kdeg.The observation points are the true model response, i.e.the Comsol simulation output, where no Kriging estimation is done.The distribution over functions enables the Kriging metamodel to predict the mean degradation depth values at different time intervals and generate the confidenc interval, which is important in systems with uncertain input parameters and parametric variability.The uncertainty in the model predictions is lower within the 95% confidenc intervals, where the probability of findin the optimized value of the desired parameters is higher.Furthermore, large confidenc intervals are related to systematic uncertainty, which indicates that assumptions underlying the physical model may need to be modifie lateron [36].

Fig.5.The mean degradation depth estimated by the Kriging metamodel at different values of the degradation rate constant kdeg [mol/m2s]at t = 28 days.The blue line presents the Kriging interpolation of the response at different values of kdeg, with 95% confidenc intervals and no calculations at the observed points.The observations are the COMSOL model output.(For interpretation of the references to colour in this figur legend, the reader is referred to the web version of this article.)

The computational time of the Kriging metamodel was within the range of 6 - 12 second for each training matrix(25 test sets) for both MDD and wt% calculations.Within COMSOL, the average computational time for each training set of the same matrices is approximately 5-23 minutes.Thereduction in the computational time didn‘t affect the accuracy of the system as can be seen in Table 4, which represents the optimized values for the model‘s parameters, with the input range and LOO eror.

3.2.Calibrated degradation model

The mean degradation depth of the model is calculated using the optimized value of the degradation rate constant(kdeg),see Table 4, obtained based on the Kriging estimations.The simulation output in comparison to the experimental data is shown in Fig.6a.The simulated MDD is in agreement with the experimental μCT imaging and weight loss data for 1-4 week with an overall NRMSE of 0.07.Detailed NRMSE values for each experimental data set can be found in Table B.5.The inset in Fig.6a shows the confidenc interval for MDD estimated by the Kriging model at approximately 14 days of degradation,which reflect the high certainty within the model predictions.The assumption of a linear degradation rate appears to approximate the degradation behaviour well,although it overestimates the degradation depth at earlier time points,both for datasets based on disc and wire degradation.Nevertheless, the variability of the experimental data is high at all time points, suggesting that a linear approximation is suffi cient.Future experiments should aim at collecting additional long-term data to verify or refute this approach.At the current stage the model also does not differentiate between different levels of Mg purity, which influence the MDD.Future work should include the development of microgalvanic cells for lower-purity Mg samples in SBF [20], e.g.by modelling the degradation of Mg using a non-linear law and adding an accelarating term dependent on the impurity level of the Mg in Eq.(19).

For comparison,we have computed the degradation rate for wires and discs respectively, based on the computed MDD,as shown in Fig.6b.The results differ based on the specimen geometry, because the degradation rate is define as

and therefore depends both on the volumeVand the surface areaAof the sample, which are shape-dependent.In addition to the degradation of the pure Mg metal, the model estimated the mean elemental wt% of each element of the precipitated degradation products.The optimized values of the reaction rate constants of each of the precipitation reactions are given in Table 4.

Fig.7 presents the solution of the model based on the optimized parameters for each element.The model response is found to fi the EDX measurements within the standard deviation of the experiments for Ca(NRMSE=0.53),C(NRMSE=0.59) and P (NRMSE =0.40).For O (NRMSE =0.99) and Mg (NRMSE = 1.03) the values differ more strongly.This could be due to the formation of other precipitates, which are not included in the current stage of the model, such as MgO or mixed MgxCay-P precipiates.Moreover, the difference in EDX measurements between different experimental conditions is larger than that for MDD.The 8-day measure-ment based on wire degradation differs strongly from that based on the degraded discs for 1-4 weeks.These differences may be due to different immersion volume to sample surface(VA) ratios in the experiments and other experimental uncertainties, such as the exact medium composition prepared by two different operators.Additionally, the current model has strongly simplifie boundary conditions that do not take into account the change of the immersion medium every 2-4 days,as performed in the experiments.This simplificatio will be valid as long as the depletion of ions due to their precipitation during the degradation process is neglectable, which is the case for high VA ratios and sufficientl frequent medium changes.Therefore, future models should entail a more detailed approach to integrate the experimental boundary conditions.

Fig.6.The mean degradation depth (a) and corresponding degradation rate (b) calculated based on the simulation and experimental data obtained using μCT.μCT measurements were performed at 7, 14, 21 and 28 days [20], 8 days [12]and at 6, 7, 9 and 10 days (Appendix A).The sample geometry differed between wires and discs.

Finally, Fig.8 shows the pH values determined in the liquid phase based on the OH-concentration for different time points in comparison to experimental minimum and maximum values from [12].A direct comparison with respect to time points is not sensible, due to the aforementioned discrepancy in boundary conditions.Nevertheless, the figur shows that the simulated pH is in the range of the experimental data.

Fig.7.The mean elemental weight percentage (wt%).calculated based on the simulation and measurements using energy-dispersive X-ray spectroscopy.EDX measurements were performed at 7,14,21 and 28 days [20]and 8 days [12].

Fig.8.pH value in liquid phase over the arc length of the domains determined by the simulation in comparison to experimental mimimum and maximum data [12].The simulated pH values are in the range of the experimental data.The sudden drop of the pH value is due to the liquid-metal interface, which moves over time.

Fig.9.Total Sobol indices of the calibrated model parameters calculated following the sensitivity analysis.

3.3.Sensitivity analysis of model parameters

Fig.9 shows the global sensitivity analysis of each of the optimized parameters.The mean degradation depth is highly influence by the degradation reaction rate constant (kdeg),while the sigmoid step initialization of the degradation (tinit)has only half the influenc on the mean degradation depth.Due to numerical decoupling of the degradation model from the precipitation, the mean degradation depth only depends on these two parameters.By contrast, the salt precipitation is more complex.For example, the amount of precipitated carbon depends on all parameters that the model is calibrated for.Similarly, the amount of precipitated calcium depends not only on reaction ratesk6andk7, which are directly related to Ca precipitation, but also onk4which models the formation of Mg3(PO4)2·8H2O(s), due to the competing reaction with phosphate ions for the formation of hydroxyapatite.The dependence of the elemental weight percentages on the reaction ratesk3,...,9is highly correlated, therefore, the summation of the Sobol indices is higher than one [35,37].The greater the correlation between parameters, the harder it will be to distinguish the influenc of each individual parameter on the quantity of interest.As a result, a simplificatio of the precipitation reactions will be difficult

Sensitivity analysis in general is an expensive approach that requires a large number of model runs under different conditions [35].To overcome this problem, the total Sobol indices were calculated based on the Kriging metamodel, which reduced the computational cost of the sensitivity analysis to 40 - 60 seconds.

4.Conclusion

Overall, the presented model simulates the anodic dissolution related degradation process of pure Mg in SBF with an error of 7%, which is similar to the experimental uncertainty.The model simulations provide the mean elemental wt% for each of the precipitated elements on the implant surface for the firs time.This information is crucial to enable the prediction of the mechanical behaviour of the degraded specimen in the future.Furthermore, it will be of particular importance when modelling the interaction with cells who are sensitive to the implant surface mechanical properties.The presented model will be valid for Mg degradation in a medium with similar ionic composition in terms of the modelled precipitation and for solutions with a similar degradation profile As the model has been calibrated for two shapes of Mg specimens with a purity between 99.92-99.95% it will be valid for pure Mg within this purity range.Thus,the degradation model presented needs to be enhanced in generality in the future by the formulation of a more general degradation reaction.The implementation of uncertainty quantificatio through the Kriging metamodel significantl reduced the computational load in the parameter estimation phase of the model and enhances the reliability of the model output.The performed sensitivity analysis reflecte the correlation between the different parameters, which explains the complexity of the presented model.This was the case in particular for the salt formation on the sample surface.

Data availability

The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Declaration of Competing Interest

The authors declare that they have no known competing financia interests or personal relationships that could have appeared to influenc the work reported in this paper.

CRediT authorship contribution statement

Berit Zeller-Plumhoff:Conceptualization, Supervision,Writing - original draft.Tamadur AlBaraghtheh:Formal analysis, Writing - original draft.Daniel H?che:Conceptualization, Writing - original draft.Regine Willumeit-R?mer:Supervision, Writing - review & editing.

Acknowledgment

The authors TAB and BZP acknowledge funding from the Helmholtz-Incubator project Uncertainty Quantification This research was supported in part through the Maxwell computational resources operated at Deutsches Elektronen-Synchrotron, Hamburg, Germany.We thank Dr.Julian Moosmann for his support at the P05 beamline and Helmholtz-Zentrum Hereon for providing beamtime.We thank Dr.Damar Wicaksono (HZDR) for the fruitful discussions regarding the UQLab software.

Appendix A.Determination of mean degradation depth for 6-10 days

Pure Mg (99.95%) wires with 1 mm diameter were degraded in SBF for 6-10 days.The degradation experiments were performed as described in [12,20].Briefl , the wires were cut into rods of 2 cm length and cleaned and sterilized in a series of n-hexane, acetone, 100% ethanol and 70%ethanol in an ultrasonic bath.They were then laid into 12-well plates and immersed in 3 ml SBF solution.The well plates were stored in an incubator to ensure physiological conditions(5%CO2,37°C).The medium was changed every 3-4 days to avoid a saturation of ions.Following degradation, the wires were rinsed in double distilled water and ethanol and then dried.They were then imaged using μCT at the P05 beamline at PETRA III, DESY, with a voxel size of 0.64 μm or 1.28 μm.The data was segmented and analysed as described in [12], and the mean degradation depth was then determined.

Appendix B.Model calibration errors per experiment data set

Table B.5 contains the detailed information on the NRMSE divided into the experiments and for the overall data.

Table B.5NRMSE of the model output parameters w.r.t to the different experimental data sets and overall.