Computational modeling of degradation process of biodegradable magnesium biomaterials

Despite the advantages of using biodegradable metals in implant design, their uncontrolled degradation and release remain a challenge in practical applications. A validated computational model of the degradation process can facilitate the tuning of implant biodegradation by changing design properties. In this study, a physicochemical model was developed by deriving a mathematical description of the chemistry of magnesium biodegradation and implementing it in a 3D computational model. The model parameters were calibrated using the experimental data of hydrogen evolution by performing a Bayesian optimization routine. The model was validated by comparing the predicted change of pH in saline and buffered solutions with the experimentally obtained values from corrosion tests, showing maximum 5% of difference, demonstrating the model's validity to be used for practical cases.


Magnesium biodegradation
Due to their bio-friendly properties, biodegradable metallic biomaterials, including magnesium (Mg), iron (Fe), and zinc (Zn), are regaining attention in recent years [1]. These biomaterials find important applications in the design and manufacturing of supportive implants such as temporary devices in orthopedics and the cardiovascular field [2,3]. In orthopedics, the biodegradable metallic biomaterials are used as fixation devices, providing adequate support in the early stages while being absorbed gradually during the bone healing process [4]. Implants fabricated using Mg and its alloys are being used for such a purpose [5] due to the similarity of the stiffness between natural bone and Mg, which helps to reduce the stress shielding induced by the implanted device. Additionally, Mg is reported to have a non-toxic contribution to the human body's metabolism and the bone healing process, which makes the release and absorption of metallic ions safe and biocompatible [6].
Accumulation of mechanistic understanding of Mg degradation achieved by experimental approaches over the years gradually provided a mechanistic understanding of the biodegradation process. Combining these insights with in silico modeling approaches enables researchers to study the biodegradation properties and behavior of the implant in a virtual environment prior to conducting any in vitro or in vivo tests. When fully validated, computational modeling can (in part) replace certain stages of costly and time-consuming experiments verifying the expected degradation behavior of the designed implants. Additionally, the developed models can be efficiently combined with existing computational models to examine other related phenomena such as tissue growth or mechanical integrity.

Computational modeling of Mg degradation
Previous contributions to the computational modeling of the degradation process include a wide range of different approaches, from the basic phenomenological implementations to comprehensive mechanistic models that take into account various aspects of the degradation and resorption process.
Continuous damage (CD) modeling has always been a common approach for corrosion simulation, but from a physicochemical point of view, it focuses on the mechanical integrity of the degradation and neglects the diffusion process. As a result, its application in the degradation modeling of biomaterials, which includes various fundamental phenomena such as mass transfer through diffusion and reaction, is relatively limited. Despite this issue, a CD model proposed by Gastaldi et al. showed a good performance for simulation of bioresorbable Mg-based medical devices [7], in which geometrical discontinuities were interpreted as the reduction of material.
Alternatively, mathematical modeling using transport phenomena equations has shown great flexibility in capturing different mechanisms involved in the biodegradation process. As an example, in Ahmed et al., a set of mathematical equations in cylindrical and spherical coordinates was derived to model the chemical reactions of Mg degradation [8]. Despite the simplicity of their approach from the computational perspective, their model was able to demonstrate the contribution of various chemical components to the in vitro degradation of Mg. Similarly, Grogan et al. developed a mathematical model based on the Stefan problem formulation in 1D space to correlate the mass flux of metallic ions into the solution to the velocity of shrinkage of the material during degradation [9]. This was done by considering the mass diffusion and change of the concentration of Mg 2+ ions, and then, employing an arbitrary Eulerian-Lagrangian (ALE) approach to extend the model to 3D on an adaptive mesh. A similar approach was taken by Shen et al. to develop a theoretical model of the degradation behavior of Mg-based orthopedic implants showing great consistency with in vitro test results [10].
An ultimate application of the computational modeling of the biodegradation process of biomaterials can be the prediction of how biodegradation affects the shape of the bulk material, medical device, or implant over time. One of the ways to achieve such a prediction is to capture the movement of the corrosion front mathematically using an appropriate method. The level set method (LSM) is a widely used example in this regard, which is an implicit mathematical way of representing the moving interfaces. This approach was used in Wilder et al. to study galvanic corrosion of metals [11]. They employed LSM on an adaptive mesh to track the moving corrosion interface, but their model lacked a thorough validation using experimental data. Gartzke et al. also worked on a simplified representation of the interface movement by developing a mechanochemical model of the biodegradation process, which helped them to study the effect of degradation on the mechanical properties [12]. They performed a basic qualitative validation on the predictions made by the model. Another similar study in this regard is the Sun et al. work [13], in which a detailed mathematical model was derived and validated to study the deposition of corrosion products on the surface of materials. This mathematical approach was also employed in the biomedical field by Bajger et al. to study the mass loss of Mg biomaterials during biodegradation [14]. They used LSM as well as a set of reaction-diffusion equations to track the change of geometry, which can be directly correlated to the loss of material. The derived equations were also able to capture the formation of the corrosion film that decreases the rate of degradation. Another comprehensive mathematical model was developed by Sanz-Herreraa et al. to investigate the role of multiple chemical components involved in the in vitro degradation of Mg implants [15]. One important drawback of this study was its 2D nature. Although the computational model was capable of studying the effect of multiple components, due to the high number of derived equations, it would be difficult to extend and use the same model for real 3D implants. Additionally, a 2D model cannot capture the full phenomenon of corrosion, and as a result, the validation of the model will be more qualitative. It was shown in the study conducted by Gao et al. [16], where they compared the results of a multi-dimension model of the degradation of cardiovascular stents with those of a singledimension model, that the number of considered dimensions had an important effect on the model predictions. In the end, it is worth mentioning that no dedicated experiments were performed in the aforementioned studies to validate the constructed mathematical and computational models.

Objective of current contribution
The current study focuses on developing a physicochemical model of the biodegradation process of commercially-pure (CP) Mg biomaterials by continuing the work of Bajger et al. [14]. In this model, a set of partial differential equations (PDE) was derived according to the underlying chemistry of biodegradation, described as reaction-diffusion processes taking place at the interface of the biomaterial and its surrounding environment. The formation of a protective layer, effects of the ions in the solution, and the change in the pH due to the corrosion phenomenon were taken into account in the mathematical model. The corresponding computational model was implemented in a fully parallelized manner. Model calibration and validation were executed using data obtained from the immersion tests performed in saline (NaCl) and simulated body fluid (SBF) solutions.

Biodegradation as a reaction-diffusion system
The biodegradation process can be considered as a reaction-diffusion system [17], in which the ions are released due to the chemical reactions on the surface, and the released ions diffuse through the surrounding solution and materials. These ions can interact with other ions and form new compounds [18]. As the reaction-diffusion systems have been studied in science and engineering for a couple of decades, the analogy with a reaction-diffusion system makes it convenient to construct a mathematical model of the biodegradation process based on the well-established transport phenomena equations [19]. From the mathematical perspective, a reaction-diffusion system is expressed by a set of parabolic PDEs that describe the conservation of contributing chemical species in the studied system.

Moving boundary -Stefan problems
Moving boundary problems, also called Stefan problems, are the general class of mathematical problems in which the boundary of the domain should also be calculated in addition to the solution of the other equations [20]. Coupling the reaction-diffusion system of biodegradation with a moving boundary problem constructs a mathematical model in which the change of the domain geometry due to the material loss can be correlated to the underlying reaction and diffusion processes of corrosion. As the geometry can be determined accurately, this approach provides a way to measure the mass loss directly by computing the change in the volume of the material. In such a system, the moving boundary is the material-solution interface (corrosion front).
For a 1D corrosion diffusion system, the position of the diffusion interface can be determined by [20]: in which the s(t) represents the position at any given time, and s 0 is the initial interface position. α coefficient can be calculated using: where [Mg] sol is the concentration in the solid bulk (i.e. materials density) and [Mg] sat is the concentration at which the material is released to the medium.
[Mg] 0 represents the initial concentration of the metallic ions in the medium, which is usually zero for most corrosion cases.
Eqs. 1 and 2 can be used to simply track the movement of the corrosion front, which is the employed method in studies like the Gorgan et al. work [9], but apparently, the real-world corrosion problems are 3D and much more complex than the described system.
As will be described later, Eq. 1 is used strictly for the first time step of the simulations in low diffusion regimes for calculating the initial velocity of the interface. Generally speaking, a more sophisticated approach, such as the level set method, is required for tracking the interface of complex 3D geometries.

Level set method
In the current study, the corrosion front is tracked using an implicit function such that the zero iso-contour of the function represents the metal-solution interface. As a common practice, this implicit function is expressed as a signed distance function that defines the distance of each point of space (the domain of interest) to the interface. Such a definition implies that the zero iso-contour of the function belongs to the interface. The level set method provides an equation to declare such an implicit function, φ = φ(x, t), x ∈ Ω ⊂ R 3 , which can be obtained by solving [21]: in which − → V E is the external velocity field, and V N is the value of the normal interface velocity. The last term is related to the curvature-dependent interface movement and is omitted. As the effect of perfusion is neglected in the current study, the term containing the external velocity is also eliminated, resulting in the following simplified form of the level set equation: By having the normal velocity of the interface (V N ) at each point and solving Eq. 4, the interface can be captured at the zero iso-contour of the φ function.

Underlying chemistry
The chemistry of biodegradation of Mg depends considerably on the surrounding solution and the presence of certain ions [18]. In NaCl solutions, the anodic and cathodic reactions as well as the formation and elimination of side corrosion products can be considered as follows [1]: Reaction 6 is not fully correct from the chemical point of view. In fact, Mg surface is always covered by MgO layer, and Mg(OH) 2 forms on top of that either at atmospheric conditions or during the immersion. The integrity of this MgO layer is undermined by Cl − ions, leading to an increase in degradation rate: Although Cl − formally does not participate in reaction 7, it reflects the dependence of Mg corrosion rate on Cl − concentration. This effect on the rate of degradation has been widely expressed as the effect of Cl − on the Mg(OH) 2 in the literature [1,3]. In the developed model, this effect is used interchangeably by omitting the MgO component, so the protective film formed on the corrosion interface is assumed to contain Mg(OH) 2 only. Moreover, it has been shown recently that oxygen reduction reaction also takes place during corrosion of Mg [22,23,24]. However, this is a secondary reaction (complementing water reduction) contributing to 1-20% of the total cathodic current depending on the conditions. Hence, it is not taken into consideration in this model. Additionally, the involved chemical reactions are more complicated in SBF solutions due to the presence of further inorganic ions and the formation of a layered precipitate structure [18], but the effect of these ions is currently encapsulated in the reaction rates and the diffusion coefficients of the developed mathematical model. The summary of the considered chemistry to develop the mathematical model is depicted in Fig. 1.

Mathematical modeling
To keep track of the concentration changes of various contributing chemical components, we define four state variables for the concentration of Mg 2+ ions, protective film (Mg(OH) 2 ), chloride (Cl − ) ions, and the hydroxide (OH − ) ions: which are indeed 4 scalar functions of space and time. Ω denotes the whole region of interest, including both the Mg bulk and its surrounding medium. By doing this, the value of pH at each point of Ω can be calculated as: where C OH implies the activity of OH − . By having the definition of the state variables in Eq. 8, the biodegradation of Mg described by Eqs. 5 and 6 can be represented as a set of reaction-diffusion PDEs: in which the maximum concentration of the protective film can be calculated according to its porosity ( ) [14]: D e is the effective diffusion coefficient for each component. Due to the formation of the protective film, the diffusion coefficient is not constant and varies from the actual diffusion coefficient of the ions to a certain fraction of it. This fraction can be defined as /τ [25,26], in which and τ are the porosity and tortuosity of the protective film, respectively. The effective diffusion coefficient can be then calculated by interpolating the two aforementioned values: The β coefficient is called momentum here and controls the effect of the saturation term . The derivation of these equations is discussed in detail in our previous work [27].

Interface movement formulation
In order to take advantage of the level set method for tracking the corrosion front, the velocity of the interface at each point should be determined. Then, by solving Eq. 4, the interface is obtained at the points with a zero value of the φ function. The interface velocity in mass transfer problems can be calculated using the Rankine-Hugoniot equation [28], and by considering the transportation of Mg 2+ ions, it can be written as: where J is the mass flux at the interface. Rearranging Eq. 16 and inserting the value of the normal interface velocity into Eq. 4 yields: which is the final form of the level set equation to be solved. In the case of simulations with a low diffusion rate, the interface moves slowly in the beginning, which results in a linear degradation, whereas based on the experimental results, the degradation rate is fast at the beginning and slows down eventually [29]. So, to mimic the same behavior in the low diffusion regimes, we took advantage of the theoretical Stefan formulation (Eqs. 1 and 2) to push the interface in the first time step. According to Eq. 1, the velocity of the interface can be calculated as (2α/ √ t), but as we are dealing with a 3D model and not a 1D one, we pick a fraction (denoted by γ) of this ideal value to be used as the driving force of the interface at the beginning of the simulation. So, the normal velocity of the interface can be written in the general form as: in which the α value should be calculated from Eq. 2. By selecting γ equal to zero, the Stefan formulation can be eliminated, and a value of 1 for γ restores the ideal 1D velocity definition.

Boundary conditions
The implementation of boundary conditions is relatively challenging and complex for the developed model as they should be imposed inside the domain of interest on virtual interfaces defined by mathematical expressions (i.e. on the moving interface defined by the zero iso-contour of the level set equation). The penalty method was used to overcome this issue and define the desired boundary conditions on the moving corrosion front. Fig. 2 demonstrates a schematic presentation of the boundary conditions and general considerations of each PDE of the biodegradation mathematical model. This figure is divided into 5 different parts, presenting the 5 PDEs of the model. The Mg block is depicted in the center, and the interface separates it from the surrounding medium. There is no specific boundary condition for the level set and film formation equations, but in comparison to the other 3 transport equations, it should be noted that diffusivity is not considered for Mg(OH) 2 , which is also reflected in Eq. 11. The level set function φ is defined in a way that is positive inside and negative outside the solid region. For the Mg 2+ ions transport equation, a Dirichlet boundary condition is applied on the mathematical interface to make the concentration equal to the saturation concentration of Mg 2+ ions, a value that was already used in Eq. 17. For the Cl − and OH − ions transport equations, a no-flux boundary condition is applied to the interface by making the diffusion coefficient equal to zero inside the Mg block, preventing ions to diffuse inside the solid material.

Implementation
To simulate the developed mathematical model, which is comprised of Eqs. 10, 11, 12, 13, and 17, a combination of finite difference and finite element methods was used, leading to discrete forms of these equations, which were subsequently solved using appropriate linear solvers. To discretize the temporal terms of the aforementioned parabolic PDEs, a first-order backward Euler finite difference scheme was used, whereas the spatial terms were converted to a weak form using a standard first-order finite element scheme. Then, the open-source PDE solver FreeFEM [30] was used to implement the weak form and obtain a linear system of equations for each PDE. The obtained linear systems were solved in parallel using the HYPRE preconditioner [31] and the GMRES solver [32] via the open-source highperformance computing (HPC) toolkit PETSc [33]. Additionally, to increase the efficiency of the computation and decrease the simulation execution time, the computational mesh was decomposed and distributed among available computing resources using the interface of HPDDM package in FreeFEM [34]. The details of this implementation are presented in our previous work [27] as well as in the supplementary materials of this paper. A simple iterative solver based on the Newton method was also developed to solve Eq. 2 to obtain the value of α parameter if it was required in the simulations.
The computational mesh was generated using a set of first-order tetrahedral elements and was adaptively refined on the metal-solution interface to increase the numerical accuracy of the simulation of the level set equation (Eq. 17). The Netgen mesh engine [35] in the SALOME platform [36] was used to generate the mesh.
Similar to the technique employed by Bajger et al. [14], the gradient of concentration of Mg 2+ in Eq. 17 was calculated at a distance h in the normal direction from the interface, with h being the smallest element size of the mesh: Considering the adaptively refined mesh, the h value is very small, so the gradient is computed at the regions close enough to the interface. In addition to this technique, the mass lumping feature of FreeFEM was used to prevent the oscillation of concentration values on the diffusive metal-medium interface.

Experimental setup
The degradation rate of CP Mg was evaluated based on the hydrogen evolution tests performed either in NaCl or SBF solutions with eudiometers. The composition of the electrolytes is shown in the following table (Table 1). 0.5 g metallic chips (with a surface area of 47.7±5.0cm 2 /g and chip thickness ca. 200 microns) of CP Mg were put in 500 ml electrolyte for 22-24 hours for monitoring the amount of evolved hydrogen. The bulk pH of electrolytes before and after corrosion was measured by a pH meter (Metrohm-691, Switzerland). Local pH was measured by positioning pH microprobes (Unisense, Denmark, pH-sensitive tip size 10x50 micron) 50 micron above the surface of Mg and monitoring the pH values either in one spot or by horizontal or vertical line-scans or mapping by following a horizontal grid. The measurements were performed at room temperature of 22 ± 2 • C maintained by the laboratory climate control system. More detailed information about experimental set up and procedures can be found elsewhere [29,37].

Parameter estimation
The constructed mathematical model contains some parameters that need to be calibrated prior to final validation of the model: diffusion coefficient of Mg 2+ and Cl − ions (D Mg and D Cl to be inserted into Eq. 15 to get effective diffusion coefficients), the reaction rates of Eqs. 5 and 6 (k 1 and k 2 ), the momentum parameter, β, for controlling the saturation term behavior (in Eqs. 10, 11, and 15), and the γ parameter for the initial interface velocity (Eq. 18). An inverse problem setup was required to estimate the proper value of these parameters.
Performing a parameter estimation requires running the computational models several times. Considering the computationally-intensive model of the current study, a sensitivity analysis was performed prior to the parameter estimation to exclude non-essential parameters and reduce the time required to complete the inverse problem run. This sensitivity analysis was accomplished separately in the low diffusion (similar to the SBF solution) and high diffusion (similar to NaCl solution) regimes.
After determining the essential parameters to include, a Bayesian optimization approach [38] was used to construct the inverse problem and calibrate the parameters. The reason for choosing a Bayesian approach was to minimize the number of optimization iterations, in each of which the simulation should run once. The Bayesian optimization is a more efficient option for such computational intensive cases in comparison to gradient-based or fully-stochastic methods as it takes into account all the preceding iterations in a probability tree [39].
The objective function of the optimization problem was the difference between the predicted and experimentally obtained values of evolved hydrogen. In the computational model, the evolved hydrogen can be computed directly at any time through the mass loss as each mole of corroded Mg is correlated to one mole of released hydrogen (Eq. 5). The mass loss can be obtained using the following volume integral: where Ω + (t) = {x : φ(x, t) ≥ 0}, and then, the amount of produced hydrogen is calculated using the ideal gas law: in which R, P , T , Mg mol are the universal gas constant, the pressure, the medium temperature, and the molar mass of Mg, respectively.

Simulation setup
In order to simulate the developed mathematical model, the experimental setup was reconstructed in silico with some minor differences. As there is no perfusion in the solution chamber, the mixing effect was neglected, so, as can also be seen in the mathematical model, the advection terms were not considered. Furthermore, the experiments were conducted using small metallic chips, yet, as the biodegradation behavior heavily depends on the exposed surfaces, we represented these chips by a cuboid with the same surface-to-mass ratio. By considering the approximate surface-to-mass of 50cm 2 /g and the total mass of 0.5g, the chips were replaced by a cuboid with the size of 60mm × 21mm × 0.2mm, which approximately has the same ratio, surface area, volume, and mass. Also, the solution chamber with a capacity of 500ml was represented by a cubic container with an edge size of 80mm. Fig. 3 depicts the constructed geometry as well as the computational mesh generated to represent it. The mesh is refined on the interface and contains 18,049,471 elements, resulting in 3,077,227 degrees of freedom (DOF) for each PDE. Simulations were carried out on the VSC (Flemish Supercomputer Center) supercomputer. Taking advantage of HPC techniques to parallelize the simulation is an inevitable aspect of such a computational-intensive model, so based on what described in the implementation section, the mesh was decomposed among 170 computing cores, i.e. 24,137 DOF per core (which includes the ghost nodes to satisfy the boundary condition in each submesh). On the VSC supercomputer, we made use of 5 nodes, 36 cores each, each node holding CPUs with a clock speed of 2.6 GHz, with 960 GB of the total available memory.
The OH − transport equation (Eq. 13) was not solved during the parameter calibration process. Afterwards, two full simulations (for the NaCl and SBF solutions) were conducted to calculate the pH changes based on the change of the concentration of OH − ions in the medium. This acted as the validation of the numerical model because no calibration was performed on the output of this equation. The pH was calculated using Eq. 9, based on the solution of Eq. 13 and a reported value of 7.00e × 10 −5 cm 2 /s (25.2mm 2 /hour) for the diffusion coefficient of OH − ions (D OH to be used in Eq. 15) in aqueous solutions [40].
According to the experimental setup, the initial concentration of the Mg 2+ , Cl − , and OH − ions were set to 0 (no Mg 2+ ions at the beginning), 146mM (5.175 × 10 −6 g/mm 3 ), and 1 × 10 −7 g/mm 3 , respectively. The porosity ( ) and tortuosity (τ ) of the protective film were considered to be 0.55 and 1, respectively [13]. The saturation concentration [Mg] sat was set to the solubility of magnesium chloride in water, which is 134×10 −6 g/mm 3 at 25 • C [41]. The density of Mg ([Mg] sol ) and Mg(OH) 2 were set to 1735 × 10 −6 g/mm 3 and 2344 × 10 −6 g/mm 3 , respectively [14]. A time step convergence study was performed to determine the implicit time step size. Based on the results, a time step with a size of 0.025 hours was chosen. The overall simulated time is 22 hours in accordance with the experimental design of performed immersion tests.

Case study
To further investigate the predictions of the current model on more complex shapes, the biodegradation of a simple screw was studied in the SBF solution using the parameters obtained for the low diffusion regimes. Similar to the simulation of Mg cuboid, the mesh was refined on the metal-medium interface, and it consisted of 1,440,439 elements with 246,580 DOFs for each PDE. All the simulation parameters and materials properties were identical to the simulation of biodegradation in the SBF solution, and the target was to simulate 42 days (1008 hours) of the process. This was selected as a sufficiently long simulated time to observe the effects of biodegradation on larger time scales.

Optimization results
Based on the performed sensitivity analysis, two parameter sets were obtained for the high diffusion (in NaCl solution) and low diffusion (in SBF solution) simulations, respectively. These parameters are listed in Table 2. According to the results, the reaction rate of Eq. 5 (k 1 ), which demonstrates the rate of oxidation-reduction, has less contribution to the process in comparison to the rate of the weakening of the protective film (k 2 ). Because of this, the parameter k 1 was not selected for the parameter estimation. Also, the model was sensitive to the effect of parameter k 2 in different ranges of values and not on a specific point, and as a result, three constant values were chosen as the delegates of these ranges in the optimization process. The model was not sensitive to the diffusion rate of Cl − ions, which was also expected because although Cl − has an important role in the weakening of the partially protective MgO film, its transport equation (Eq. 12) is purely diffusive and does not include any reaction term.
The parameter optimization process was performed on the specified range of selected parameters, while the rest of the parameter values were obtained from the literature [14,40]. Table 3 shows the output of this process, which was used to simulate the full model. For two estimation processes, 120 optimization iterations were taken cumulatively, which took 276 hours of simulation execution time using 170 computing nodes for each simulation.   Fig. 4 shows the model output for the predicted produced hydrogen, protective film formation, and the pH changes. The graph of the evolved hydrogen is used as input during the parameter optimization process, but the pH results are produced by the simulations using the optimized parameters to demonstrate the validation of the developed mathematical and computational models. The predicted pH result (Fig. 4-d) shows a difference of 5.35% for the simulation in NaCl and 1.03% for SBF simulation. Each simulation took about 3 hours to complete.

Degradation prediction
In Fig. 4, a post-processed view of the final shape of the Mg cuboid in the NaCl solution is presented, in which the degraded geometry is plotted on the Mg 2+ ions (Fig. 4-b) and protective film concentration (Fig. 4-c) contours. A transparent contour of the pH values in the solution is depicted for both the NaCl (Fig. 4-e) and SBF (Fig. 4-f) solutions. The range of colors is kept equal for both contours to make it easy to compare the change of pH in both solutions.

Example application
The simulation of 42 days (19,200 time steps) of the degradation of the simple screw took 9 hours to run using 170 computing cores. Fig. 6 depicts the post-processed interface and Mg 2+ ions release (similar to Fig. 4-b) as well as the mass loss during the degradation of the screw in the SBF solution.

Discussion
In this study, a physicochemical model of the biodegradation process of commerciallypure Mg was developed by constructing a mathematical model formulating the mass transfer phenomena as well as tracking the location of the surface of the implant during degradation. For the mass transfer model, the equations were derived from the chemistry of biodegradation of the Mg in saline (NaCl) and buffered (SBF) solutions, which includes the oxidation of the metallic part, reduction of water, changes in pH, and formation of a protective film on the surface of the scaffold which contributes to a slower rate of degradation. Beside these aspects, it was also crucial to consider the effect of different ions in the medium on the rate of degradation. Additionally, investigating the structural changes of the scaffolds and implants in practical applications, like resorption of temporary fixation devices, requires tracking the movement of the corrosion surface. This was done by constructing an equation based on the level set principle, which captured the movement of the medium-metal interface by defining an implicit surface. The derived equations were coupled and solved using a combination of finite difference and finite element methods. The degradation data to validate the models was collected from immersion tests of small Mg chips, reconstructed as a single cuboid in the computational study with a similar surface over volume properties. The model parameters were calibrated using a Bayesian optimization algorithm, and the obtained parameters were used to simulate the pH changes in NaCl and SBF solutions.
The developed model falls in the categories of physical models of the corrosion process, which provide more insights of the process in comparison to the phenomenological models. The reason is that the phenomenological models focus on the elimination of elements to capture the loss of materials, which makes it impossible to model the formation of new chemical compounds or interaction of species [42]. The physical models, like the one developed in this study, are capable of capturing the underlying chemical interactions. By doing this, processes like the effect of coating, the formation of a protective layer, and pH changes can be modeled. Adding an appropriate interface tracking method enables the physical models to act like the phenomenological models in capturing the corrosion interface movement. In the current study, this has been accomplished using a level set function. Technically speaking, this approach has certain benefits over the ALE method, which is the method used by several similar studies, including Grogan et al. [9]. In comparison to the ALE method, the level set function tracks the interface instead of a Lagrangian mesh, and elements can freely be marked as solid or liquid. Additionally, employing the ALE method for degradation simulation requires remeshing the geometry as the interface moves, which is not efficient for 3D models and is limited to the available features of the employed numerical solver.
One of the challenging aspects of validating physical models is getting the correct value for the parameters of said models, requiring dedicated experimental input. To overcome this challenge, an efficient inverse problem was constructed based on the Bayesian optimization approach to estimate the unknown parameters. To save time and resources, the parameter estimation process was performed on the most effective parameters, which were selected based on a sensitivity analysis. This selection process implied the importance of parameters in high and low diffusion rates.
The degradation rate is fast at the beginning, but then it slows down due to the formation of a partially protective film and also because of the saturation concentration. This phenomenon is well captured by the model at high diffusion rates, but in low diffusion rates (in SBF solution), this effect can be reproduced by pushing the corrosion front according to the Stefan formulation of the moving interface problems. This was controlled by the parameter γ (Eq. 18). It should be noted that the inclusion of the γ parameter is crucial for short-term simulations only, helping the model mimic the chemical behavior correctly. Defining and considering γ is necessary because from the computational costs perspective, performing the parameter calibration on simulations with thousands of time steps requires a lot more resources and time.
The degradation of the CP Mg was assumed to be mostly diffusion-based. As a result, the value of D Mg plays an important role in the behavior of the model. Although it was possible to get the diffusion coefficient of Mg 2+ from the previously conducted experiments in the literature (similar to what was done for D OH ), we decided to not do so because of two reasons: 1) the values reported in the literature are mostly valid for saline solutions only, and 2) the reported values were not in a good agreement with one another [9,13]. Thus, the diffusion coefficient was obtained using the parameter estimation process for both the NaCl and SBF solutions. The obtained value of D Mg (0.06273mm 2 /hour) was in line with the values that Grogan et al. have already suggested (0.010575 − 0.50575mm 2 /hour) [9], showing that the constructed inverse problem was successful in reproducing previous results of similar studies. The obtained value of D Mg in the Bajger et al. work [14] is 0.00066mm 2 /hour, which is mostly related to the simplicity of the employed parameter estimation method as well as having a 2D model instead of a 3D one.
In the in vitro biodegradation of Mg-based biomaterials, the local pH of the surrounding solution increases less than that in NaCl solution. This is because the Mg(OH) 2 formed in NaCl stabilizes pH at 10.4 [43], while Mg-Ca-P-C containing products stabilize the pH at 7.8-8.5 since OH-is consumed for the formation of this product [44,29]. This phenomenon was captured in Eqs. 13 and 9 to calculate pH based on the concentration of OH − ions, showing the local pH changes at any location ( Fig. 4-e,f). In the current study, the global pH is considered as the validation criterion, which means that the average value of the solution pH is calculated using a volume integral and is compared with the ones obtained from the experiments. Fig. 4-d shows that such a prediction has a good agreement with the experimental data. It is worth noting that the model can be extended to non-pure Mg by considering the effect of alloying elements on the reaction rates as well as adding more terms to the transport equations to capture the electrochemical potential changes, converting the PDEs into the Nernst-Planck equation [45]. By doing so, more complex forms of the corrosion process, such as galvanic corrosion, can be predicted by the model.
One of the biggest simplifications of the current study was made by ignoring the contribution of pH changes to the biodegradation mechanism of Mg. Although doing that is relatively simple and straightforward in the approach taken by this study, it results in nonlinear terms in the derived PDEs. This non-linearity inserts another level of complexity to the computational model as the order of the state variables are in the range of 10 −6 to 10 −10 , which makes it difficult to yield convergence in the iterative non-linear solvers. By developing a robust non-linear solver, this effect can be added simply by including more relevant terms as the effect of Eq. 13 into Eq. 10.
Additionally, buffered solutions and the physiological fluids inside the human and animal bodies contain more ions interacting with more complex chemistry [18]. In this study, this effect was encapsulated in a limited number of parameters (such as k 1 and k 2 in Eqs. 10, 11, and 13), but while the results show its success to reproduce experimental observations, it still needs additional elaboration to be able to capture more chemical interactions. For SBF solutions, the effect of presented inorganic ions such as HCO − 3 , HPO 2− 4 /H 2 PO − 4 , and Ca 2+ can be added similar to the way the effect of Cl − was considered. Additionally, formulating the effect of HPO 2− 4 that exists in the physiological environments will make the model capable of making more accurate predictions for in vivo studies.
Although the pH simulations are not enough experimental input to call the model fully validated, the obtained validation results show that the derived mathematical model and the corresponding parallelized computational model give a correct in silico representation of the studied process. The performed predictive simulations, including the case study, demonstrate the potential of the developed computational model and software to study the biodegradation behavior of implants. This can be further combined with other computational models to provide a multidisciplinary environment to investigate the mechanical integrity of implants or induced neotissue growth for different applications in orthopedics and tissue engineering.

Conclusions
The use of biodegradable metals for designing medical devices and implants has the challenge of controlling the release and rate of degradation, which is usually investigated by conducting in vitro and in vivo tests requiring conducting multiple experiments for different scenarios and situations. In this study, we have developed a mathematical model to predict the biodegradation behavior of commercially pure Mg-based biomaterials, which makes it possible to study the corrosion of implants and scaffolds in a simulated environment. Despite the assumed simplifications, the model can serve as an important tool to find the biodegradable metals properties and predict the biodegradation behavior of Mg-based implants that improves current design workflows.