Computational modeling of degradation process of biodegradable magnesium biomaterials
Mojtaba Barzegari, Di Mei, Sviatlana V. Lamaka, Liesbet Geris
CComputational modeling of degradation process of biodegradablemagnesium biomaterials
Mojtaba Barzegari a, ∗ , Di Mei b,c , Sviatlana V. Lamaka b , Liesbet Geris a,d a Biomechanics Section, Department of Mechanical Engineering, KU Leuven, Leuven, Belgium b Institute of Surface Science, Helmholtz-Zentrum Geesthacht, Geesthacht, Germany c School of Materials Science and Engineering & Henan Key Laboratory of Advanced Magnesium Alloy,Zhengzhou University, Zhengzhou, PR China d Biomechanics Research Unit, GIGA in Silico Medicine, University of Liege, Liege, Belgium
Abstract
Despite the advantages of using biodegradable metals in implant design, their uncontrolleddegradation and release remain a challenge in practical applications. A validated computa-tional model of the degradation process can facilitate the tuning of implant biodegradationby changing design properties. In this study, a physicochemical model was developed byderiving a mathematical description of the chemistry of magnesium biodegradation and im-plementing it in a 3D computational model. The model parameters were calibrated usingthe 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 bufferedsolutions with the experimentally obtained values from corrosion tests, showing maximum5% of difference, demonstrating the model’s validity to be used for practical cases.
Keywords:
Degradable metals, Magnesium corrosion, Reaction-diffusion models, Finiteelement method, In silico medicine
1. Introduction
Due to their bio-friendly properties, biodegradable metallic biomaterials, including mag-nesium (Mg), iron (Fe), and zinc (Zn), are regaining attention in recent years [1]. Thesebiomaterials find important applications in the design and manufacturing of supportiveimplants such as temporary devices in orthopedics and the cardiovascular field [2, 3]. Inorthopedics, the biodegradable metallic biomaterials are used as fixation devices, providing ∗ Corresponding author, Tel: (+32) 16 193831
Email addresses: [email protected] (Mojtaba Barzegari), [email protected] (Di Mei), [email protected] (Sviatlana V. Lamaka), [email protected] (Liesbet Geris), [email protected] (Liesbet Geris)
URL: (Liesbet Geris)
Preprint submitted to Corrosion Science February 22, 2021 a r X i v : . [ c s . C E ] F e b dequate support in the early stages while being absorbed gradually during the bone healingprocess [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 reducethe stress shielding induced by the implanted device. Additionally, Mg is reported to havea 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 experimentalapproaches over the years gradually provided a mechanistic understanding of the biodegra-dation process. Combining these insights with in silico modeling approaches enables re-searchers to study the biodegradation properties and behavior of the implant in a virtualenvironment prior to conducting any in vitro or in vivo tests. When fully validated, compu-tational modeling can (in part) replace certain stages of costly and time-consuming experi-ments verifying the expected degradation behavior of the designed implants. Additionally,the developed models can be efficiently combined with existing computational models toexamine other related phenomena such as tissue growth or mechanical integrity. Previous contributions to the computational modeling of the degradation process includea wide range of different approaches, from the basic phenomenological implementations tocomprehensive mechanistic models that take into account various aspects of the degradationand resorption process.Continuous damage (CD) modeling has always been a common approach for corrosionsimulation, but from a physicochemical point of view, it focuses on the mechanical integrityof the degradation and neglects the diffusion process. As a result, its application in thedegradation modeling of biomaterials, which includes various fundamental phenomena suchas mass transfer through diffusion and reaction, is relatively limited. Despite this issue, a CDmodel proposed by Gastaldi et al. showed a good performance for simulation of bioresorbableMg-based medical devices [7], in which geometrical discontinuities were interpreted as thereduction of material.Alternatively, mathematical modeling using transport phenomena equations has showngreat 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 sphericalcoordinates was derived to model the chemical reactions of Mg degradation [8]. Despite thesimplicity of their approach from the computational perspective, their model was able todemonstrate the contribution of various chemical components to the in vitro degradation ofMg. Similarly, Grogan et al. developed a mathematical model based on the Stefan problemformulation in 1D space to correlate the mass flux of metallic ions into the solution to thevelocity of shrinkage of the material during degradation [9]. This was done by consideringthe mass diffusion and change of the concentration of Mg ions, and then, employing anarbitrary Eulerian-Lagrangian (ALE) approach to extend the model to 3D on an adaptivemesh. A similar approach was taken by Shen et al. to develop a theoretical model of thedegradation behavior of Mg-based orthopedic implants showing great consistency with invitro test results [10]. 2n ultimate application of the computational modeling of the biodegradation processof biomaterials can be the prediction of how biodegradation affects the shape of the bulkmaterial, medical device, or implant over time. One of the ways to achieve such a predictionis to capture the movement of the corrosion front mathematically using an appropriatemethod. The level set method (LSM) is a widely used example in this regard, which is animplicit mathematical way of representing the moving interfaces. This approach was used inWilder et al. to study galvanic corrosion of metals [11]. They employed LSM on an adaptivemesh to track the moving corrosion interface, but their model lacked a thorough validationusing experimental data. Gartzke et al. also worked on a simplified representation of theinterface movement by developing a mechanochemical model of the biodegradation process,which helped them to study the effect of degradation on the mechanical properties [12]. Theyperformed a basic qualitative validation on the predictions made by the model. Anothersimilar study in this regard is the Sun et al. work [13], in which a detailed mathematicalmodel was derived and validated to study the deposition of corrosion products on the surfaceof materials. This mathematical approach was also employed in the biomedical field byBajger et al. to study the mass loss of Mg biomaterials during biodegradation [14]. Theyused 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 alsoable 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. toinvestigate the role of multiple chemical components involved in the in vitro degradationof Mg implants [15]. One important drawback of this study was its 2D nature. Althoughthe computational model was capable of studying the effect of multiple components, dueto the high number of derived equations, it would be difficult to extend and use the samemodel for real 3D implants. Additionally, a 2D model cannot capture the full phenomenonof corrosion, and as a result, the validation of the model will be more qualitative. It wasshown in the study conducted by Gao et al. [16], where they compared the results of amulti-dimension model of the degradation of cardiovascular stents with those of a single-dimension model, that the number of considered dimensions had an important effect onthe model predictions. In the end, it is worth mentioning that no dedicated experimentswere performed in the aforementioned studies to validate the constructed mathematical andcomputational models. The current study focuses on developing a physicochemical model of the biodegradationprocess 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 tothe underlying chemistry of biodegradation, described as reaction-diffusion processes takingplace at the interface of the biomaterial and its surrounding environment. The formation ofa protective layer, effects of the ions in the solution, and the change in the pH due to thecorrosion phenomenon were taken into account in the mathematical model. The correspond-ing computational model was implemented in a fully parallelized manner. Model calibration3nd validation were executed using data obtained from the immersion tests performed insaline (NaCl) and simulated body fluid (SBF) solutions.
2. Background theory
The biodegradation process can be considered as a reaction-diffusion system [17], inwhich the ions are released due to the chemical reactions on the surface, and the released ionsdiffuse through the surrounding solution and materials. These ions can interact with otherions and form new compounds [18]. As the reaction-diffusion systems have been studiedin science and engineering for a couple of decades, the analogy with a reaction-diffusionsystem makes it convenient to construct a mathematical model of the biodegradation processbased on the well-established transport phenomena equations [19]. From the mathematicalperspective, a reaction-diffusion system is expressed by a set of parabolic PDEs that describethe conservation of contributing chemical species in the studied system.
Moving boundary problems, also called Stefan problems, are the general class of math-ematical problems in which the boundary of the domain should also be calculated in ad-dition to the solution of the other equations [20]. Coupling the reaction-diffusion systemof biodegradation with a moving boundary problem constructs a mathematical model inwhich the change of the domain geometry due to the material loss can be correlated to theunderlying reaction and diffusion processes of corrosion. As the geometry can be determinedaccurately, this approach provides a way to measure the mass loss directly by computingthe change in the volume of the material. In such a system, the moving boundary is thematerial-solution interface (corrosion front).For a 1D corrosion diffusion system, the position of the diffusion interface can be deter-mined by [20]: s ( t ) = s + 2 α √ t, (1)in which the s ( t ) represents the position at any given time, and s is the initial interfaceposition. α coefficient can be calculated using: α = [Mg] − [Mg] sat [Mg] sol − [Mg] sat (cid:114) Dπ exp (cid:16) − α D (cid:17) erfc (cid:16) − α √ D (cid:17) (2)where [Mg] sol is the concentration in the solid bulk (i.e. materials density) and [Mg] sat is theconcentration at which the material is released to the medium. [Mg] represents the initialconcentration of the metallic ions in the medium, which is usually zero for most corrosioncases.Eqs. 1 and 2 can be used to simply track the movement of the corrosion front, whichis the employed method in studies like the Gorgan et al. work [9], but apparently, thereal-world corrosion problems are 3D and much more complex than the described system.4s will be described later, Eq. 1 is used strictly for the first time step of the simulations inlow 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 theinterface of complex 3D geometries. In the current study, the corrosion front is tracked using an implicit function such thatthe zero iso-contour of the function represents the metal-solution interface. As a commonpractice, this implicit function is expressed as a signed distance function that defines thedistance of each point of space (the domain of interest) to the interface. Such a definitionimplies that the zero iso-contour of the function belongs to the interface. The level setmethod provides an equation to declare such an implicit function, φ = φ ( x , t ) , x ∈ Ω ⊂ R ,which can be obtained by solving [21]: ∂φ∂t + −→ V E · ∇ φ + V N |∇ φ | = bκ |∇ φ | (3)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. Asthe effect of perfusion is neglected in the current study, the term containing the externalvelocity is also eliminated, resulting in the following simplified form of the level set equation: ∂φ∂t + V N |∇ φ | = 0 . (4)By having the normal velocity of the interface (V N ) at each point and solving Eq. 4, theinterface can be captured at the zero iso-contour of the φ function.
3. Materials and methods
The chemistry of biodegradation of Mg depends considerably on the surrounding solutionand the presence of certain ions [18]. In NaCl solutions, the anodic and cathodic reactionsas well as the formation and elimination of side corrosion products can be considered asfollows [1]: Mg + 2H O k → Mg + H + 2OH − k → Mg(OH) + H (5)Mg(OH) + 2Cl − k → Mg + 2Cl − + 2OH − . (6)Reaction 6 is not fully correct from the chemical point of view. In fact, Mg surface isalways covered by MgO layer, and Mg(OH) forms on top of that either at atmosphericconditions or during the immersion. The integrity of this MgO layer is undermined by Cl − ions, leading to an increase in degradation rate:MgO + Cl − + H O k → Mg + Cl − + 2OH − . (7)5lthough Cl − formally does not participate in reaction 7, it reflects the dependenceof Mg corrosion rate on Cl − concentration. This effect on the rate of degradation hasbeen widely expressed as the effect of Cl − on the Mg(OH) in the literature [1, 3]. In thedeveloped model, this effect is used interchangeably by omitting the MgO component, sothe protective film formed on the corrosion interface is assumed to contain Mg(OH) only.Moreover, it has been shown recently that oxygen reduction reaction also takes place duringcorrosion of Mg [22, 23, 24]. However, this is a secondary reaction (complementing waterreduction) 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 chemicalreactions are more complicated in SBF solutions due to the presence of further inorganicions and the formation of a layered precipitate structure [18], but the effect of these ions iscurrently encapsulated in the reaction rates and the diffusion coefficients of the developedmathematical model. The summary of the considered chemistry to develop the mathematicalmodel is depicted in Fig. 1. Figure 1: The chemistry of biodegradation of Mg considered in the current study: 1) Mg oxidation andwater reduction processes accompanied by releasing Mg and OH − ions as well as H gas, 2) formation ofa partially protective precipitation layer, 3) dynamic solubility equilibrium and contribution of Cl − . To keep track of the concentration changes of various contributing chemical components,we define four state variables for the concentration of Mg ions, protective film (Mg(OH) ),chloride (Cl − ) ions, and the hydroxide (OH − ) ions: C Mg = C Mg ( x , t ) , C Film = C Film ( x , t ) C Cl = C Cl ( x , t ) , C OH = C OH ( x , t ) x ∈ Ω ⊂ R , (8)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 ateach point of Ω can be calculated as:pH = 14 + log C OH , (9)6here C OH implies the activity of OH − . By having the definition of the state variables inEq. 8, the biodegradation of Mg described by Eqs. 5 and 6 can be represented as a set ofreaction-diffusion PDEs: ∂C Mg ∂t = ∇ · (cid:0) D e Mg ∇ C Mg (cid:1) − k C Mg (cid:18) − β C Film [Film] max (cid:19) + k C Film C Cl2 (10) ∂C Film ∂t = k C Mg (cid:18) − β C Film [Film] max (cid:19) − k C Film C Cl2 (11) ∂C Cl ∂t = ∇ · ( D e Cl ∇ C Cl ) (12) ∂C OH ∂t = ∇ · ( D e OH ∇ C OH ) + k C Film C Cl2 (13)in which the maximum concentration of the protective film can be calculated according toits porosity ( (cid:15) ) [14]: [Film] max = ρ Mg(OH) × (1 − (cid:15) ) . (14) D e is the effective diffusion coefficient for each component. Due to the formation of theprotective film, the diffusion coefficient is not constant and varies from the actual diffusioncoefficient of the ions to a certain fraction of it. This fraction can be defined as (cid:15)/τ [25, 26],in which (cid:15) and τ are the porosity and tortuosity of the protective film, respectively. Theeffective diffusion coefficient can be then calculated by interpolating the two aforementionedvalues: D ei = D i (cid:18)(cid:18) − β C Film [Film] max (cid:19) + β C Film [Film] max (cid:15)τ (cid:19) . (15)The β coefficient is called momentum here and controls the effect of the saturation term(1 − C Film [Film] max ). The derivation of these equations is discussed in detail in our previous work[27].
In order to take advantage of the level set method for tracking the corrosion front, thevelocity of the interface at each point should be determined. Then, by solving Eq. 4, theinterface is obtained at the points with a zero value of the φ function. The interface velocityin mass transfer problems can be calculated using the Rankine–Hugoniot equation [28], andby considering the transportation of Mg ions, it can be written as: { J ( x, t ) − ([Mg] sol − [Mg] sat ) V( x, t ) } · n = 0 (16)where J is the mass flux at the interface. Rearranging Eq. 16 and inserting the value of thenormal interface velocity into Eq. 4 yields: ∂φ∂t − D e Mg ∇ n C Mg [Mg] sol − [Mg] sat |∇ φ | = 0 , (17)7hich is the final form of the level set equation to be solved. In the case of simulationswith a low diffusion rate, the interface moves slowly in the beginning, which results in alinear degradation, whereas based on the experimental results, the degradation rate is fastat the beginning and slows down eventually [29]. So, to mimic the same behavior in the lowdiffusion 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 interfacecan 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 theinterface at the beginning of the simulation. So, the normal velocity of the interface can bewritten in the general form as:V N ( x, t ) = (cid:40) γ α √ t t = 0 D e Mg ∇ n C Mg [Mg] sol − [Mg] sat t > α value should be calculated from Eq. 2. By selecting γ equal to zero, theStefan formulation can be eliminated, and a value of 1 for γ restores the ideal 1D velocitydefinition. The implementation of boundary conditions is relatively challenging and complex forthe developed model as they should be imposed inside the domain of interest on virtualinterfaces defined by mathematical expressions (i.e. on the moving interface defined by thezero iso-contour of the level set equation). The penalty method was used to overcome thisissue and define the desired boundary conditions on the moving corrosion front.Fig. 2 demonstrates a schematic presentation of the boundary conditions and generalconsiderations of each PDE of the biodegradation mathematical model. This figure is dividedinto 5 different parts, presenting the 5 PDEs of the model. The Mg block is depicted in thecenter, and the interface separates it from the surrounding medium. There is no specificboundary condition for the level set and film formation equations, but in comparison tothe other 3 transport equations, it should be noted that diffusivity is not considered forMg(OH) , which is also reflected in Eq. 11. The level set function φ is defined in a waythat is positive inside and negative outside the solid region. For the Mg ions transportequation, a Dirichlet boundary condition is applied on the mathematical interface to makethe concentration equal to the saturation concentration of Mg ions, a value that wasalready used in Eq. 17. For the Cl − and OH − ions transport equations, a no-flux boundarycondition is applied to the interface by making the diffusion coefficient equal to zero insidethe Mg block, preventing ions to diffuse inside the solid material. 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, leadingto discrete forms of these equations, which were subsequently solved using appropriate linearsolvers. 8 igure 2: A schematic overview of the exposed boundary conditions and constraints required for the simu-lation of each equation of the developed mathematical model for Mg biodegradation.
To discretize the temporal terms of the aforementioned parabolic PDEs, a first-orderbackward Euler finite difference scheme was used, whereas the spatial terms were convertedto a weak form using a standard first-order finite element scheme. Then, the open-sourcePDE solver FreeFEM [30] was used to implement the weak form and obtain a linear sys-tem of equations for each PDE. The obtained linear systems were solved in parallel us-ing the HYPRE preconditioner [31] and the GMRES solver [32] via the open-source high-performance computing (HPC) toolkit PETSc [33]. Additionally, to increase the efficiencyof the computation and decrease the simulation execution time, the computational meshwas decomposed and distributed among available computing resources using the interfaceof HPDDM package in FreeFEM [34]. The details of this implementation are presented inour previous work [27] as well as in the supplementary materials of this paper. A simpleiterative solver based on the Newton method was also developed to solve Eq. 2 to obtainthe value of α parameter if it was required in the simulations.The computational mesh was generated using a set of first-order tetrahedral elementsand was adaptively refined on the metal-solution interface to increase the numerical accuracyof the simulation of the level set equation (Eq. 17). The Netgen mesh engine [35] in theSALOME platform [36] was used to generate the mesh.Similar to the technique employed by Bajger et al. [14], the gradient of concentration ofMg 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: ∇ n C = C ( x + h.n ) − C ( x + 2 h.n ) h x ∈ Ω ⊂ R . (19)9onsidering the adaptively refined mesh, the h value is very small, so the gradient is com-puted at the regions close enough to the interface. In addition to this technique, the masslumping feature of FreeFEM was used to prevent the oscillation of concentration values onthe diffusive metal-medium interface. The degradation rate of CP Mg was evaluated based on the hydrogen evolution testsperformed either in NaCl or SBF solutions with eudiometers. The composition of the elec-trolytes is shown in the following table (Table 1). 0.5 g metallic chips (with a surface area of47 . ± . / g and chip thickness ca. 200 microns) of CP Mg were put in 500 ml electrolytefor 22-24 hours for monitoring the amount of evolved hydrogen. The bulk pH of electrolytesbefore and after corrosion was measured by a pH meter (Metrohm-691, Switzerland). LocalpH was measured by positioning pH microprobes (Unisense, Denmark, pH-sensitive tip size10x50 micron) 50 micron above the surface of Mg and monitoring the pH values either inone spot or by horizontal or vertical line-scans or mapping by following a horizontal grid.The measurements were performed at room temperature of 22 ± ◦ C maintained by thelaboratory climate control system. More detailed information about experimental set upand procedures can be found elsewhere [29, 37].
Table 1: Chemical composition of NaCl and SBF electrolytes used to perform hydrogen evolution tests,weight loss, local and bulk pH measurements.
Concentration/ mM0.85 wt. % NaCl SBFNa + + - 5.0Mg - 1.5Ca - 2.5Cl − − - 4.2HPO − / H PO − - 1.0SO − - 4.2Synthetic pH buffer(i.e. Tris/HCl, HEPES) No NoInitial pH value 5.6-5.9 7.35-7.45 The constructed mathematical model contains some parameters that need to be cali-brated prior to final validation of the model: diffusion coefficient of Mg and Cl − ions( D Mg and D Cl to be inserted into Eq. 15 to get effective diffusion coefficients), the reactionrates of Eqs. 5 and 6 ( k and k ), the momentum parameter, β , for controlling the satura-tion term behavior (in Eqs. 10, 11, and 15), and the γ parameter for the initial interface10elocity (Eq. 18). An inverse problem setup was required to estimate the proper value ofthese parameters.Performing a parameter estimation requires running the computational models severaltimes. Considering the computationally-intensive model of the current study, a sensitivityanalysis was performed prior to the parameter estimation to exclude non-essential param-eters and reduce the time required to complete the inverse problem run. This sensitivityanalysis was accomplished separately in the low diffusion (similar to the SBF solution) andhigh 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 reasonfor 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 moreefficient option for such computational intensive cases in comparison to gradient-based orfully-stochastic methods as it takes into account all the preceding iterations in a probabilitytree [39].The objective function of the optimization problem was the difference between the pre-dicted 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 eachmole of corroded Mg is correlated to one mole of released hydrogen (Eq. 5). The mass losscan be obtained using the following volume integral:Mg lost = (cid:90) Ω + ( t ) [Mg] sol d V − (cid:90) Ω + (0) [Mg] sol d V, (20)where Ω + ( t ) = { x : φ ( x , t ) ≥ } , and then, the amount of produced hydrogen is calculatedusing the ideal gas law: H f = Mg lost Mg mol RTP (21)in which R , P , T , Mg mol are the universal gas constant, the pressure, the medium temper-ature, and the molar mass of Mg, respectively. In order to simulate the developed mathematical model, the experimental setup was re-constructed in silico with some minor differences. As there is no perfusion in the solutionchamber, 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 us-ing small metallic chips, yet, as the biodegradation behavior heavily depends on the exposedsurfaces, we represented these chips by a cuboid with the same surface-to-mass ratio. Byconsidering the approximate surface-to-mass of 50cm / g and the total mass of 0 . × × . , ,
471 elements, resulting in 3 , , Figure 3: Computational representation of the experimental set-up, used to perform parameter estimationand numerical validation of the developed model. a) A cuboid of Mg (60mm × × . Simulations were carried out on the VSC (Flemish Supercomputer Center) supercom-puter. Taking advantage of HPC techniques to parallelize the simulation is an inevitableaspect of such a computational-intensive model, so based on what described in the imple-mentation section, the mesh was decomposed among 170 computing cores, i.e. 24 ,
137 DOFper core (which includes the ghost nodes to satisfy the boundary condition in each sub-mesh). On the VSC supercomputer, we made use of 5 nodes, 36 cores each, each nodeholding 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 calibrationprocess. Afterwards, two full simulations (for the NaCl and SBF solutions) were conductedto calculate the pH changes based on the change of the concentration of OH − ions in themedium. This acted as the validation of the numerical model because no calibration wasperformed on the output of this equation. The pH was calculated using Eq. 9, based onthe solution of Eq. 13 and a reported value of 7 . e × − cm / s (25 . / hour) for thediffusion 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 , Cl − , andOH − ions were set to 0 (no Mg ions at the beginning), 146mM (5 . × − g / mm ), and1 × − g / mm , respectively. The porosity ( (cid:15) ) and tortuosity ( τ ) of the protective film wereconsidered to be 0.55 and 1, respectively [13]. The saturation concentration [Mg] sat was setto the solubility of magnesium chloride in water, which is 134 × − g / mm at 25 ◦ C [41]. Thedensity of Mg ([Mg] sol ) and Mg(OH) were set to 1735 × − g / mm and 2344 × − g / mm ,respectively [14]. A time step convergence study was performed to determine the implicittime step size. Based on the results, a time step with a size of 0 .
025 hours was chosen. Theoverall simulated time is 22 hours in accordance with the experimental design of performedimmersion tests. 12 able 2: The effective parameters as the result of the sensitivity analysis and their corresponding range tobe considered in the Bayesian optimization for parameter calibration
Parameter Optimization rangeLow diffusion(SBF solution) D Mg [0 . , . k , , β [0 . , γ [0 , D Mg [0 . , . k , , β [0 . , 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 parametersobtained for the low diffusion regimes. Similar to the simulation of Mg cuboid, the mesh wasrefined on the metal-medium interface, and it consisted of 1 , ,
439 elements with 246 ,
4. Results
Based on the performed sensitivity analysis, two parameter sets were obtained for thehigh diffusion (in NaCl solution) and low diffusion (in SBF solution) simulations, respec-tively. These parameters are listed in Table 2. According to the results, the reaction rateof Eq. 5 ( k ), which demonstrates the rate of oxidation-reduction, has less contribution tothe process in comparison to the rate of the weakening of the protective film ( k ). Becauseof this, the parameter k was not selected for the parameter estimation. Also, the modelwas sensitive to the effect of parameter k in different ranges of values and not on a specificpoint, and as a result, three constant values were chosen as the delegates of these rangesin 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 ofthe partially protective MgO film, its transport equation (Eq. 12) is purely diffusive anddoes not include any reaction term.The parameter optimization process was performed on the specified range of selectedparameters, 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. Fortwo estimation processes, 120 optimization iterations were taken cumulatively, which took276 hours of simulation execution time using 170 computing nodes for each simulation.13 able 3: Values used to evaluate the model performance, obtained from the output of the optimizationprocess and the literature. Parameter D Mg D Cl D OH k k β γ Unit mm hour mm hour mm hour 1hour mm g hour - -SBF solution 0 . .
05 25 . .
125 0 . . .
05 25 . . Figure 4: Comparing the quantitative output of the model for the rate of degradation and the pH changes inNaCl and SBF solutions with experimentally measured values as well as the simulation results for ion release,mass loss, protective film formation, and pH changes after 22 hours of simulated time: a) calibrated outputof the formed hydrogen gas during the degradation (the SBF curves are overlapped), b) the simulationresults of Mg ions release, c) the simulation results of protective film concentration at the end of thesimulation (the color contour shows the concentration of species, and the gray surface is the zero iso-contourof the level set function, which indicates the surface of the Mg block), d) de novo prediction of the global pHchanges in the medium, showing a good agreement between the model output and the experimental results,e) pH changes in different regions of the medium in NaCl solution, f) pH changes in SBF solution. Fig. 4 shows the model output for the predicted produced hydrogen, protective filmformation, and the pH changes. The graph of the evolved hydrogen is used as input duringthe parameter optimization process, but the pH results are produced by the simulations usingthe optimized parameters to demonstrate the validation of the developed mathematical andcomputational models. The predicted pH result (Fig. 4-d) shows a difference of 5 .
35% forthe simulation in NaCl and 1 .
03% for SBF simulation. Each simulation took about 3 hoursto complete.In Fig. 4, a post-processed view of the final shape of the Mg cuboid in the NaCl solutionis presented, in which the degraded geometry is plotted on the Mg ions (Fig. 4-b) andprotective film concentration (Fig. 4-c) contours. A transparent contour of the pH valuesin the solution is depicted for both the NaCl (Fig. 4-e) and SBF (Fig. 4-f) solutions. Therange of colors is kept equal for both contours to make it easy to compare the change of pHin both solutions. 14 igure 5: The change of concentration for the involved chemical components, Mg , Cl − , OH − ions, andthe protective precipitation structure (which can be correlated to the thickness of the layer) plotted over adiagonal line as shown in the right.Figure 6: A cross-section of the computational mesh and simulation results of the degradation process ofthe use-case screw in SBF solution as well as the mass loss graph over time. The contours display theconcentration of Mg ions on a cross-section view of the medium beside the moving surface of the screwat 1) 1st day (initial state), 2) 6th day, 3) 12th day, and 4) 18th day. The concentration values of the state variables of the derived transport PDEs (Mg ,Cl − , OH − , and Mg(OH) ) are plotted along a diagonal line in the solution container in Fig.5, showing how they change in the zones close to the corrosion surface and far from it. Theanimated output of the degradation of the Mg cuboid is presented as a set of movie filesin the supplementary materials of this paper (post-processed using NVIDIA IndeX softwareon a GPU). The simulation of 42 days (19 ,
200 time steps) of the degradation of the simple screwtook 9 hours to run using 170 computing cores. Fig. 6 depicts the post-processed interfaceand Mg ions release (similar to Fig. 4-b) as well as the mass loss during the degradationof the screw in the SBF solution. 15 . Discussion In this study, a physicochemical model of the biodegradation process of commercially-pure Mg was developed by constructing a mathematical model formulating the mass transferphenomena 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 biodegradationof the Mg in saline (NaCl) and buffered (SBF) solutions, which includes the oxidation ofthe metallic part, reduction of water, changes in pH, and formation of a protective film onthe surface of the scaffold which contributes to a slower rate of degradation. Beside theseaspects, it was also crucial to consider the effect of different ions in the medium on the rate ofdegradation. Additionally, investigating the structural changes of the scaffolds and implantsin practical applications, like resorption of temporary fixation devices, requires tracking themovement of the corrosion surface. This was done by constructing an equation based on thelevel set principle, which captured the movement of the medium-metal interface by definingan implicit surface. The derived equations were coupled and solved using a combination offinite difference and finite element methods. The degradation data to validate the modelswas collected from immersion tests of small Mg chips, reconstructed as a single cuboid in thecomputational study with a similar surface over volume properties. The model parameterswere calibrated using a Bayesian optimization algorithm, and the obtained parameters wereused 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 cap-ture the loss of materials, which makes it impossible to model the formation of new chemicalcompounds or interaction of species [42]. The physical models, like the one developed inthis 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 canbe modeled. Adding an appropriate interface tracking method enables the physical modelsto act like the phenomenological models in capturing the corrosion interface movement. Inthe current study, this has been accomplished using a level set function. Technically speak-ing, this approach has certain benefits over the ALE method, which is the method used byseveral similar studies, including Grogan et al. [9]. In comparison to the ALE method, thelevel set function tracks the interface instead of a Lagrangian mesh, and elements can freelybe marked as solid or liquid. Additionally, employing the ALE method for degradationsimulation requires remeshing the geometry as the interface moves, which is not efficient for3D 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 valuefor the parameters of said models, requiring dedicated experimental input. To overcome thischallenge, an efficient inverse problem was constructed based on the Bayesian optimizationapproach to estimate the unknown parameters. To save time and resources, the parameterestimation process was performed on the most effective parameters, which were selectedbased on a sensitivity analysis. This selection process implied the importance of parametersin high and low diffusion rates. 16he degradation rate is fast at the beginning, but then it slows down due to the for-mation of a partially protective film and also because of the saturation concentration. Thisphenomenon 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 accordingto the Stefan formulation of the moving interface problems. This was controlled by theparameter γ (Eq. 18). It should be noted that the inclusion of the γ parameter is crucialfor 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 requiresa 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 waspossible to get the diffusion coefficient of Mg from the previously conducted experimentsin the literature (similar to what was done for D OH ), we decided to not do so because oftwo reasons: 1) the values reported in the literature are mostly valid for saline solutionsonly, 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 boththe NaCl and SBF solutions. The obtained value of D Mg (0 . / hour) was in linewith the values that Grogan et al. have already suggested (0 . − . / hour)[9], showing that the constructed inverse problem was successful in reproducing previousresults of similar studies. The obtained value of D Mg in the Bajger et al. work [14] is0 . / hour, which is mostly related to the simplicity of the employed parameterestimation 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 surroundingsolution increases less than that in NaCl solution. This is because the Mg(OH) formed inNaCl stabilizes pH at 10.4 [43], while Mg-Ca-P-C containing products stabilize the pH at7.8-8.5 since OH- is consumed for the formation of this product [44, 29]. This phenomenonwas 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 globalpH is considered as the validation criterion, which means that the average value of thesolution pH is calculated using a volume integral and is compared with the ones obtainedfrom the experiments. Fig. 4-d shows that such a prediction has a good agreement withthe experimental data. It is worth noting that the model can be extended to non-pureMg by considering the effect of alloying elements on the reaction rates as well as addingmore terms to the transport equations to capture the electrochemical potential changes,converting the PDEs into the Nernst–Planck equation [45]. By doing so, more complexforms 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 con-tribution of pH changes to the biodegradation mechanism of Mg. Although doing that isrelatively simple and straightforward in the approach taken by this study, it results in non-linear terms in the derived PDEs. This non-linearity inserts another level of complexityto the computational model as the order of the state variables are in the range of 10 − to10 − , which makes it difficult to yield convergence in the iterative non-linear solvers. By17eveloping a robust non-linear solver, this effect can be added simply by including morerelevant terms as the effect of Eq. 13 into Eq. 10.Additionally, buffered solutions and the physiological fluids inside the human and animalbodies contain more ions interacting with more complex chemistry [18]. In this study, thiseffect was encapsulated in a limited number of parameters (such as k and k in Eqs. 10, 11,and 13), but while the results show its success to reproduce experimental observations, itstill needs additional elaboration to be able to capture more chemical interactions. For SBFsolutions, the effect of presented inorganic ions such as HCO − , HPO − / H PO − , and Ca can be added similar to the way the effect of Cl − was considered. Additionally, formulatingthe effect of HPO − that exists in the physiological environments will make the model capableof making more accurate predictions for in vivo studies.Although the pH simulations are not enough experimental input to call the model fullyvalidated, the obtained validation results show that the derived mathematical model andthe corresponding parallelized computational model give a correct in silico representationof the studied process. The performed predictive simulations, including the case study,demonstrate the potential of the developed computational model and software to studythe biodegradation behavior of implants. This can be further combined with other com-putational models to provide a multidisciplinary environment to investigate the mechanicalintegrity of implants or induced neotissue growth for different applications in orthopedicsand tissue engineering.
6. Conclusions
The use of biodegradable metals for designing medical devices and implants has thechallenge of controlling the release and rate of degradation, which is usually investigated byconducting in vitro and in vivo tests requiring conducting multiple experiments for differentscenarios and situations. In this study, we have developed a mathematical model to predictthe biodegradation behavior of commercially pure Mg-based biomaterials, which makes itpossible to study the corrosion of implants and scaffolds in a simulated environment. Despitethe assumed simplifications, the model can serve as an important tool to find the biodegrad-able metals properties and predict the biodegradation behavior of Mg-based implants thatimproves current design workflows.
Acknowledgment
This research is financially supported by the Prosperos project, funded by the InterregVA Flanders – The Netherlands program, CCI grant no. 2014TC16RFCB046 and by theFund for Scientific Research Flanders (FWO), grant G085018N. LG acknowledges supportfrom the European Research Council under the European Union’s Horizon 2020 researchand innovation programme, ERC CoG 772418. The computational resources and servicesused in this work were provided by the VSC (Flemish Supercomputer Center), funded by theResearch Foundation - Flanders (FWO) and the Flemish Government – department EWI.18 eferences [1] Y. Zheng, X. Gu, F. Witte, Biodegradable metals, Materials Science and Engineering: R: Reports 77(2014) 1–34. doi:10.1016/j.mser.2014.01.001 .[2] Y. Chen, Z. Xu, C. Smith, J. Sankar, Recent advances on the development of magnesium alloys forbiodegradable implants, Acta Biomaterialia 10 (11) (2014) 4561–4573. doi:10.1016/j.actbio.2014.07.005 .[3] D. Zhao, F. Witte, F. Lu, J. Wang, J. Li, L. Qin, Current status on clinical applications of magnesium-based orthopaedic implants: A review from clinical translational perspective, Biomaterials 112 (2017)287–302. doi:10.1016/j.biomaterials.2016.10.017 .[4] Y. Qin, P. Wen, H. Guo, D. Xia, Y. Zheng, L. Jauer, R. Poprawe, M. Voshage, J. H. Schleifenbaum,Additive manufacturing of biodegradable metals: Current research status and future perspectives, ActaBiomaterialia 98 (2019) 3 – 22. doi:10.1016/j.actbio.2019.04.046 .[5] U. Riaz, I. Shabib, W. Haider, The current trends of mg alloys in biomedical applications—a review,Journal of Biomedical Materials Research Part B: Applied Biomaterials (dec 2018). doi:10.1002/jbm.b.34290 .[6] Y. Xin, K. Huo, H. Tao, G. Tang, P. K. Chu, Influence of aggressive ions on the degradation behavior ofbiomedical magnesium alloy in physiological environment, Acta Biomaterialia 4 (6) (2008) 2008–2015. doi:10.1016/j.actbio.2008.05.014 .[7] D. Gastaldi, V. Sassi, L. Petrini, M. Vedani, S. Trasatti, F. Migliavacca, Continuum damage modelfor bioresorbable magnesium alloy devices — application to coronary stents, Journal of the MechanicalBehavior of Biomedical Materials 4 (3) (2011) 352–365. doi:10.1016/j.jmbbm.2010.11.003 .[8] S. Ahmed, J. Ward, Y. Liu, Numerical modelling of effects of biphasic layers of corrosion products tothe degradation of magnesium metal in vitro, Materials 11 (1) (2017) 1. doi:10.3390/ma11010001 .[9] J. Grogan, S. Leen, P. McHugh, A physical corrosion model for bioabsorbable metal stents, ActaBiomaterialia 10 (5) (2014) 2313–2322. doi:10.1016/j.actbio.2013.12.059 .[10] Z. Shen, M. Zhao, D. Bian, D. Shen, X. Zhou, J. Liu, Y. Liu, H. Guo, Y. Zheng, Predicting thedegradation behavior of magnesium alloys with a diffusion-based theoretical model and in vitro corrosiontesting, Journal of Materials Science & Technology 35 (7) (2019) 1393–1402. doi:10.1016/j.jmst.2019.02.004 .[11] J. W. Wilder, C. Clemons, D. Golovaty, K. L. Kreider, G. W. Young, R. S. Lillard, An adaptive level setapproach for modeling damage due to galvanic corrosion, Journal of Engineering Mathematics 91 (1)(2014) 121–142. doi:10.1007/s10665-014-9732-3 .[12] A.-K. Gartzke, S. Julmi, C. Klose, A.-C. Waselau, A. Meyer-Lindenberg, H. J. Maier, S. Besdo, P. Wrig-gers, A simulation model for the degradation of magnesium-based bone implants, Journal of the Me-chanical Behavior of Biomedical Materials 101 (2020) 103411. doi:10.1016/j.jmbbm.2019.103411 .[13] W. Sun, G. Liu, L. Wang, T. Wu, Y. Liu, An arbitrary lagrangian–eulerian model for studying the in-fluences of corrosion product deposition on bimetallic corrosion, Journal of Solid State Electrochemistry17 (3) (2012) 829–840. doi:10.1007/s10008-012-1935-9 .[14] P. Bajger, J. M. A. Ashbourn, V. Manhas, Y. Guyot, K. Lietaert, L. Geris, Mathematical modellingof the degradation behaviour of biodegradable metals, Biomechanics and Modeling in Mechanobiology16 (1) (2016) 227–238. doi:10.1007/s10237-016-0812-3 .[15] J. Sanz-Herrera, E. Reina-Romo, A. Boccaccini, In silico design of magnesium implants: Macroscopicmodeling, Journal of the Mechanical Behavior of Biomedical Materials 79 (2018) 181–188. doi:10.1016/j.jmbbm.2017.12.016 .[16] Y. Gao, L. Wang, X. Gu, Z. Chu, M. Guo, Y. Fan, A quantitative study on magnesium alloy stentbiodegradation, Journal of Biomechanics 74 (2018) 98–105. doi:10.1016/j.jbiomech.2018.04.027 .[17] Y. Wang, J. Pan, X. Han, C. Sinka, L. Ding, A phenomenological model for the degradation ofbiodegradable polymers, Biomaterials 29 (23) (2008) 3393–3401.[18] D. Mei, S. V. Lamaka, X. Lu, M. L. Zheludkevich, Selecting medium for corrosion testing of bioab-sorbable magnesium and other metals – a critical review, Corrosion Science 171 (2020) 108722. doi:10.1016/j.corsci.2020.108722 .
19] P. Grindrod, The theory and applications of reaction-diffusion equations : patterns and waves, OxfordUniversity Press, 1996.[20] J. Crank, Free and Moving Boundary Problems, OUP Oxford, 1987.[21] S. O. Ronald Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer New York, 2002.[22] C. Wang, D. Mei, G. Wiese, L. Wang, M. Deng, S. V. Lamaka, M. L. Zheludkevich, High rate oxygenreduction reaction during corrosion of ultra-high-purity magnesium, npj Materials Degradation 4 (1)(dec 2020). doi:10.1038/s41529-020-00146-1 .[23] M. Strebl, M. Bruns, S. Virtanen, Editors’ choice—respirometric in situ methods for real-time moni-toring of corrosion rates: Part i. atmospheric corrosion, Journal of The Electrochemical Society 167 (2)(2020) 021510. doi:10.1149/1945-7111/ab6c61 .[24] E. L. Silva, S. V. Lamaka, D. Mei, M. L. Zheludkevich, The reduction of dissolved oxygen duringmagnesium corrosion, ChemistryOpen 7 (8) (2018) 664–668. doi:10.1002/open.201800076 .[25] P. Grathwohl, Diffusion in Natural Porous Media: Contaminant Transport, Sorption/Desorption andDissolution Kinetics, Springer US, 1998. doi:10.1007/978-1-4615-5683-1 .[26] D. H¨oche, Simulation of corrosion product deposit layer growth on bare magnesium galvanically cou-pled to aluminum, Journal of The Electrochemical Society 162 (1) (2014) C1–C11. doi:10.1149/2.0071501jes .[27] M. Barzegari, L. Geris, Highly scalable numerical simulation of coupled reaction-diffusion systems withmoving interfaces arXiv:http://arxiv.org/abs/2008.11057v1 .[28] S. Scheiner, C. Hellmich, Stable pitting corrosion of stainless steel as diffusion-controlled dissolutionprocess with a sharp moving electrode boundary, Corrosion Science 49 (2) (2007) 319–346. doi:10.1016/j.corsci.2006.03.019 .[29] D. Mei, S. V. Lamaka, J. Gonzalez, F. Feyerabend, R. Willumeit-R¨omer, M. L. Zheludkevich, The roleof individual components of simulated body fluid on the corrosion behavior of commercially pure mg,Corrosion Science 147 (2019) 81–93. doi:10.1016/j.corsci.2018.11.011 .[30] F. Hecht, New development in freefem++, J. Numer. Math. 20 (3-4) (2012) 251–265. doi:10.1515/jnum-2012-0013 .URL https://freefem.org/ [31] R. D. Falgout, U. M. Yang, hypre: A library of high performance preconditioners, in: Lecture Notes inComputer Science, Springer Berlin Heidelberg, 2002, pp. 632–641. doi:10.1007/3-540-47789-666 .[32] Y. Saad, M. H. Schultz, Gmres: A generalized minimal residual algorithm for solving nonsymmetriclinear systems, SIAM Journal on Scientific and Statistical Computing 7 (3) (1986) 856–869. arXiv:https://doi.org/10.1137/0907058 , doi:10.1137/0907058 .URL https://doi.org/10.1137/0907058 [33] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener,V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T.Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Webpage, (2019).URL [34] P. Jolivet, F. Hecht, F. Nataf, C. Prud’homme, Scalable domain decomposition preconditioners forheterogeneous elliptic problems, in: Proceedings of the International Conference on High PerformanceComputing, Networking, Storage and Analysis, SC ’13, Association for Computing Machinery, NewYork, NY, USA, 2013. doi:10.1145/2503210.2503212 .URL https://doi.org/10.1145/2503210.2503212 [35] J. Sch¨oberl, NETGEN an advancing front 2d/3d-mesh generator based on abstract rules, Computingand Visualization in Science 1 (1) (1997) 41–52. doi:10.1007/s007910050004 .[36] A. Ribes, C. Caremoli, Salome platform component model for numerical simulation, in: 31st AnnualInternational Computer Software and Applications Conference - Vol. 2 - (COMPSAC 2007), IEEE,2007. doi:10.1109/compsac.2007.185 .[37] D. Mei, S. V. Lamaka, C. Feiler, M. L. Zheludkevich, The effect of small-molecule bio-relevant organiccomponents at low concentration on the corrosion of commercially pure mg and mg-0.8ca alloy: An verall perspective, Corrosion Science 153 (2019) 258–271. doi:10.1016/j.corsci.2019.03.039 .[38] J. Mockus, Bayesian Approach to Global Optimization, Springer Netherlands, 1989. doi:10.1007/978-94-009-0909-0 .[39] M. Mehrian, Y. Guyot, I. Papantoniou, S. Olofsson, M. Sonnaert, R. Misener, L. Geris, Maximizingneotissue growth kinetics in a perfusion bioreactor: An in silico strategy using model reduction andbayesian optimization, Biotechnology and Bioengineering 115 (3) (2017) 617–629. doi:10.1002/bit.26500 .[40] S. H. Lee, J. C. Rasaiah, Proton transfer and the mobilities of the H+ and OH- ions from studiesof a dissociating model for water, The Journal of Chemical Physics 135 (12) (2011) 124505. doi:10.1063/1.3632990 .[41] J. S. H. Graham C. Hill, Chemistry in Context - Laboratory Manual, Oxford University Press, 2001.[42] M. Abdalla, A. Joplin, M. Elahinia, H. Ibrahim, Corrosion modeling of magnesium and its alloys forbiomedical applications: Review, Corrosion and Materials Degradation 1 (2) (2020) 219–248. doi:10.3390/cmd1020011 .[43] R. J. Santucci, M. E. McMahon, J. R. Scully, Utilization of chemical stability diagrams for improvedunderstanding of electrochemical systems: evolution of solution chemistry towards equilibrium, npjMaterials Degradation 2 (1) (jan 2018). doi:10.1038/s41529-017-0021-2 .[44] S. V. Lamaka, J. Gonzalez, D. Mei, F. Feyerabend, R. Willumeit-R¨omer, M. L. Zheludkevich, LocalpH and its evolution near mg alloy surfaces exposed to simulated body fluids, Advanced MaterialsInterfaces 5 (18) (2018) 1800169. doi:10.1002/admi.201800169 .[45] K. B. Deshpande, Numerical modeling of micro-galvanic corrosion, Electrochimica Acta 56 (4) (2011)1737–1745. doi:10.1016/j.electacta.2010.09.044 ..