In silico modelling of the corrosion of biodegradable magnesium-based biomaterials: modelling approaches, validation and future perspectives.

Metallic biomedical implants based on magnesium, zinc and iron alloys have emerged as bioresorbable alternatives to permanent orthopaedic implants over the last two decades. The corrosion rate of biodegradable metals plays a critical role in controlling the compatibility and functionality of the device in vivo. The broader adoption of biodegradable metals in orthopaedic applications depends on developing in vitro methods that accurately predict the biodegradation behaviour in vivo. However, the physiological environment is a highly complex corrosion environment to replicate in the laboratory, making the in vitro-to-in vivo translation of results very challenging. Accordingly, the results from in vitro corrosion tests fail to provide a complete schema of the biodegradation behaviour of the metal in vivo. In silico approach based on computer simulations aim to bridge the observed differences between experiments performed in vitro and vivo. A critical review of the state-of-the-art of computational modelling techniques for predicting the corrosion behaviour of magnesium alloy as a biodegradable metal is presented.

Introduction

Internal fixation in bone fracture management is achieved by using a combination of plates and screws that act together to stabilize the fractured bone and/or loose fragments. Titanium is currently the gold standard biomaterial for internal fixation due to its high mechanical properties and relative inertness in the physiological environment.[1] However, there are a number of limitations associated with titanium implants including stress-shielding,[2] thermal sensitivity,[3-6] and release of toxic metallic ions,[7-10] that in some cases may necessitate revision surgery.[11, 12] Stress-shielding that can result in osteoporosis in the case of permanent metallic implants.[13, 14] Magnesium (Mg) alloys have emerged as a promising alternative metallic biomaterial due to their resorbability and mechanochemical compatibility in vivo,[15-21] with recent clinical trials demonstrating their efficacy in orthopaedic applications.[22-24] However, excessively high rates of Mg corrosion are detrimental to the mechanical performance of the implant due to premature loss of fixation. Moreover, an excessive corrosion rate results in the evolution of hydrogen gas that may accumulate at the implant site if the rate of corrosion is higher than the rate of hydrogen gas diffusion into the surrounding tissues.[25] Fortunately, the corrosion rate and mechanical properties of Mg alloys can be controlled via chemical composition,[26] thermomechanical processing,[27-29] and surface coatings.[30, 31]

Generally, the corrosion rate of Mg alloys is higher in vitro than that observed in vivo due to differences in the type and amounts of inorganic species (Na+, K+, Ca2+, HCO3-, Cl-, HPO42-, SO42-) and/or organic compounds (e.g., proteins) included in the in vitro test design.[32] The above variables may influence the corrosion rate through surface reactions and coatings (e.g., hydroxyapatite formation). Further, the precise dynamics of the local physiological conditions (including pH, ion concentration, surface coatings, etc.) are not well defined for different implant sites. Consequently, it is not known how to precisely mimic the in vivo conditions in an in vitro test.[33-35] Additionally, the stresses and strains introduced during surgical implantation of devices or bone remodelling process are typically not considered during in vitro testing.[17, 29, 36-41] The measurement of stresses and strains within the implanted material is not performed routinely, although these factors may increase the rate and type of corrosion in Mg alloys.[42-45] Finally, the forms of corrosion (e.g., general, pitting, fretting, filiform, crevice) that manifest during in vitro and vivo testing may differ, bringing into question the validity of in vitro testing for experimentally modelling the in vivo conditions.

More recently, in silico approaches to modelling the biodegradation of Mg alloys have been reported as a possible means of bridging the gap between in vitro and vivo experimental findings.[46] The present work is a critical review of the state-of-the-art in terms of corrosion modelling techniques used to predict the biodegradation behaviour of Mg-based biodegradable implants. An emphasis is placed on the theoretical aspects of computational simulation and how these techniques are calibrated and validated via experimental results. The search for peer-reviewed articles from 1991 to 2021 was conducted via the University of Canterbury online library. These databases included Scopus, SciFinder, Materials science & engineering database, IEEE Xplore and Compendex. Google Scholar was also utilized to locate open access articles. Different variations of following terms were used to ensure exhaustively to locate articles specific to this study: corrosion modelling of bioresorbable metals, corrosion modelling of magnesium, finite element modelling of corrosion and so on.

Corrosion of Mg-Based Alloys

Song et al.[36] have reviewed the corrosion mechanisms of Mg alloys in detail. A thin (20-50 nm) surface film of MgO is observed to form on the Mg substrate with exposure to air.[47] The corrosion reaction proceeds with the evolution of H2 gas from cathodic regions and release of Mg2+ ions from anodic regions in an aqueous environment. Mg corrosion may be facilitated by the presence of impurities or secondary phases with differing electrochemical potentials to the matrix phase, rupture in the MgO film, and/or crystallographic misorientation between adjacent grains. The overall corrosion reaction for Mg results in the formation of Mg hydroxide (Mg(OH)2) that forms a partially protective surface layer equations (1-3).

Mg and its alloys are susceptible to micro-galvanic corrosion, pitting corrosion, stress-corrosion cracking and corrosion fatigue. In general, micro-galvanic corrosion results from the Mg matrix and secondary phases acting as anodic or cathodic sites. Impurities such as ferrum, nickel or cupper with a low hydrogen overvoltage can also cause micro-galvanic corrosion.[48] Various studies have also demonstrated that Mg alloys are susceptible to pitting corrosion and transgranular stress corrosion in a physiological environment that may lead to sudden failure of the implant.[36, 39-41] Mg substrates that are coated with a protective layer may also be susceptible to filiform corrosion.[29, 37] Fretting corrosion may occur on mating surfaces such as in a screw and plate fixation system.[38]

The corrosion reactions associated with Mg-based implant devices in the presence of physiological fluids are complex, involving several electrochemical processes resulting from interactions between Mg2+ ions and other ionic species present in the physiological environment.[49] Based on several in vitro corrosion studies, the corrosion of Mg-based bioresorbable materials in a physiological environment would likely be dependent on several factors including the chemical composition,[50] microstructure,[51] surface roughness,[52, 53] mechanical stress,[44] and plastic strain introduced during implantation.[42, 43, 45] The chemomechanical effect describes the influence of chemical reaction (i.e., corrosion) on the mechanical properties of the Mg (e.g., yield strength),[54] while the effect of mechanical loads on corrosion behaviour can be referred to as the mechanochemical effect (e.g., stress corrosion).[55]

Modelling the Corrosion Rate of Metals

Computational modelling of metal alloys in different corrosive environments has been of research interest over the past couple of decades. The approaches typically used to model corrosion of metals can be broadly classified into three categories: probabilistic modelling, physics-based modelling and finite-element based modelling. A probabilistic modelling approach is adopted to quantify the effect of corrosion on structural integrity and associated uncertainty of the predicted effects due to factors not considered in the model. Probabilistic models are focused on estimating the time to failure of a full-scale industrial component for a specific corrosive environment.[56, 57] Physics-based models on the other hand, focus on a particular corrosion phenomenon such as galvanic corrosion,[58, 59] fatigue cracking,[60] or stress corrosion cracking.[61] These models are ideal to understand the corrosion mechanisms involved during material degradation. However, these models do not predict the overall corrosion rate of a metal.[61] The continuum damage mechanics (CDM)-based finite element method (FEM) modelling approach is a promising approach to model corrosion of metals subjected to complex corrosive environments.[62-65] The CDM-FEM approach has been used to model mass loss and loss of mechanical integrity of a component as a function of time.

Finite Element Modelling of Magnesium Corrosion

It is major challenge to develop computational models that embody all of the chemomechanical (and mechanochemical) effects for predicting corrosion of Mg-based implants in vivo.[66, 67] Probabilistic and physics-based modelling approaches attempt to consider the multiphysics of the corrosion process. However, it is not practical to implement such models to routinely describe the biodegradation of orthopaedic implants in vivo due to the (i) difficulty of carrying out experimental-based calibration procedures for all of the model parameters, and (ii) associated high computational cost of three-dimensional geometries. The FEM is a promising tool since it makes possible the simulation of different corrosion mechanisms simultaneously through the implementation of multiple models. FEM is a general numerical method for solving partial differential equations in the two-dimensional (2D) and three-dimensional (3D) domain. The domain is subdivided into a finite number of smaller parts, referred to as finite elements. The overall corrosion over the entire domain can be modelled as a set of models that predict corrosion of each element as a function of time. It is important to note that FEM is discrete in both position and time and hence it is important to ensure that timestep and length scale associated with each element (characteristic length) is small enough to be able to provide sufficient resolution to capture the underlying physical phenomenon. There is no universal definition of characteristic length. However, characteristic length can be considered to be length of a side cube whose volume is same as the element in 3D elements or length of a side square whose area is same as the element in 2D elements. FEM has been used to model the time-dependent chemomechanical behaviour of Mg-based biodegradable materials.[68] FEM models are relatively straightforward to calibrate with experimental data and may be implemented using commonly-available proprietary software. Additionally, a CDM-based FEM modelling approach has been successfully implemented to model different types of corrosion phenomenon over past couple of decades. FEM-based models designed to simulate the corrosion of Mg alloys mainly fall into one of three categories: phenomenological, physical or hybrid models.

Phenomenological modelling

Phenomenological approaches are based on CDM theory that relates the presence of geometrical discontinuities (e.g., dimensional changes due to loss of material via corrosion) to changes in the macroscopic mechanical properties of the material (e.g., elastic modulus, yield strength, etc.).[69] A damage field (D) is used to describe changes in mechanical properties as a function of position and time. D varies from 0 to 1 where 0 indicates an uncorroded element, while 1 indicates the complete corrosion of the element. The evolution of damage due to uniform corrosion with time is given by equation (4).[68]

where, k (time-1) is a parameter related to the kinetics of uniform corrosion, L is the characteristic length of the particular element, α is a dimensionless corrosion susceptibility parameter (α ≥ 0), and δc is a characteristic dimension of the uniform corrosion process (e.g., thickness of the corrosion film). The value of α is 0 for an element not exposed to the corrosion medium (and greater than 0 if exposed). A constant value of α (typically set to 1) would correspond to a uniform corrosion rate. Technically, the chosen value of α has no impact on uniform corrosion model as it scales (up or down) the calibrated value of k. The role of α in the CDM-based phenomenological corrosion models is to introduce variability to the corrosion rates of different elements, enabling a specific element in the model to undergo uniform or localised corrosion. For example, the value of α is randomly chosen from a Weibull distribution centred around constant value (usually 1) in pitting corrosion models. However, α can be any constant value in the case of uniform corrosion. Various phenomenological corrosion models presented in literature may be further classified as one of four broad categories: uniform, pitting, stress- and strain-corrosion.

Uniform corrosion

The modelling of uniform corrosion assumes that the mass loss occurs homogeneously over the material surface, providing an approximate representation of the micro-galvanic corrosion that occurs in Mg.[68] For example, Gastaldi et al.[68] present a phenomenological model of uniform corrosion for application to coronary stents. The model is based on CDM in which α of all surface elements is given an initial value of 1, while α = 0 for all internal elements. The value of δc can then be selected so that the model closely aligns with the in vitro test data. Gastaldi et al. assigned a value of 100 µm to δc, based on the average thickness of the corrosion layer observed for a Mg alloy (AM60B-F) under dynamic flow conditions in vitro.[70] However, the translation of such in vitro-calibrated models to predictions of corrosion rates in vivo is problematic since the corrosion layer in vivo is like to be different in terms of growth rate due to drastically different corrosion rates reported in vivo compared to in vitro corrosion rates and chemical composition when compared with that observed in vitro.[71, 72] Thus, an in vitro-based approach for the selection of δc would likely result in phenomenological models that are not valid for prediction of corrosion rates in vivo.

Gorgan et al.[73] present a similar phenomenological model to that of Gastaldi et al.[68] However, the approach of Gorgan et al.[73] was to define δc based on the average grain size (17 μm) of the alloy being studied (AZ31). The determination of the average grain size of an equiaxed microstructure is relatively straightforward. However, a microstructure generated by thermomechanical processing (e.g., warm extrusion or rolling) in which the grains become highly elongated may lead to uncertainty in the optimal selection of δc.

The role of the parameter δc is to make the models mesh independent. Mesh independent models predict the same corrosion rate for a particular alloy regardless of the element size in the FEM model which is a user decision based on factors such as geometry of corrosion specimen, computational cost, etc. The approach used to select δc becomes unimportant if the ratio of δc to Le is approximately constant for different specimen geometries. However, this is not always possible due to complex geometry and thus it is critically important to identify which of the two approaches is physically relevant based on in vitro corrosion studies. Furthermore, it is important to identify whether adopted values of δc are physically relevant in vivo, to avoid unintended errors in predicted corrosion rates. At this stage, neither approach for the selection of δc has been validated by experiments carried out in vivo.

Gorgan et al.[73] and Gastaldi et al.[68] also used different methods to compute Le via equations (5) and (6), respectively. Interestingly, Gastaldi et al.[68, 74] evaluated Le for a 3D model with hexahedral elements using an equation that is more appropriate for 2D axisymmetric simulations.

where, Ve is the volume of a hexahedral element.

where, Ae is the cross-sectional area of the element and r is the circumradius.

Both of the models above were implemented in ABAQUS (Dassault Systàmes, Vélizy-Villacoublay, France) in which elements were deleted for D ≥ 1. The above approach based on hexahedral elements could also be used with different types of 2D and 3D elements due to the simple element deletion approach used. However, a major difference between the two above models is the manner in which the corrosion proceeds following deletion of an element. Firstly, the neighbouring elements of the deleted element are identified as the new surface elements provided they lie within a predefined radius equation (7).

where, ρ is the radius of influence from the centroid of the deleted element, ρo is a constant, Le is the characteristic length of deleted element, and Le (max) is the largest characteristic length in the FEM model.

The value of α for all new surface elements was then set to 1 to indicate their exposure to the corrosion medium.[68] Two major drawbacks of using this approach is the (i) dependence of the corrosion rate on the largest characteristic length, and (ii) arbitrary selection of ρo. Firstly, there is no direct correlation between the variance in the grain size (which is a function of thermomechanical processing, alloy composition, etc.) and that of the characteristic length (which is a function of geometry). For example, consider two meshed geometries with the same average element size, where one of the meshes contains a relatively large element. The presence of the anomalously large element would result in a significantly higher predicted corrosion rate compared to a meshed geometry having a more homogenous element size. Secondly, ρo is a constant that represents the corrosion length scale. Consequently, the model would predict that two different elements with same Le would corrode at same rate regardless of their depth from the surface. However, the progressive formation of a passivation layer as corrosion proceeds would be expected to decrease the corrosion rate of elements that lie increasingly further below the original surface elements.

In contrast, Gorgan et al.[73] consider the neighbouring elements to be the common nodes between elements. Additionally, α of each neighbouring element, n, was calculated with equation (8).

where, αe is α of the deleted element ‘e’, αn is α of the neighbouring element ‘n’, and β is a positive dimensionless parameter that captures changes in the corrosion rate as a function of time due to the formation of surface passivation layers. β for each element (βe) is scaled according to equation (9) for each element to eliminate effect of characteristic length of the element on corrosion rate

where, Le is the characteristic length of the element. L0 and β0 are the characteristic length and corresponding value of parameter β respectively. The simulations carried out during calibration deploy a simple geometry (e.g., cubic) where all elements are cubic with L0 as length of each side. The calibration experiments are then performed using specimens with same geometry to identify all model parameters including parameter β0.

Gorgan et al.[73] determined kc and β0 from hydrogen gas evolution experiments. However, the corrosion rate estimated by hydrogen evolution can be a poor estimate of the long-term corrosion rate due to increases in the pH over time if the corrosion medium is not buffered.[74] Gastaldi et al.[68] instead considered kc to be a constant that was determined from an immersion mass loss experiment. Mass loss testing is more reliable method for estimating the long term corrosion rate in a physiological environment as pH control is more practical.[34]

Pitting corrosion

Mg-based biometals are susceptible to pitting corrosion that may lead to premature failure of the implant.[36] Gorgan et al.[73] modified their uniform corrosion model to account for pitting by stochastically assigning the value of α using a Weibull distribution rather than setting α to 1 for every surface element. The probability (P) of α for an element ‘e’ (αe) lying in the range [a, b] can be given by equation (10).

where, f(x) is the probability distribution function given by equation (11).

where, γ and ψ are shape and scaling parameters, respectively. Gorgan et al.[73] implemented their FEM model using ABAQUS and the VUMAT subroutine. Others have also adopted this FEM model, utilising the UMAT subroutine to reduce computational cost.[75, 76] Gorgan et al.[73] proposed that kc, γ, ψ and β0 can be determined via hydrogen gas evolution measurements, while the model was validated by comparing with in vitro test data, namely the ultimate tensile strength (MPa) of corroded specimens and time-to-failure (hour) under constant load corrosion. The specimens were subjected to a constant tensile load while being immersed in the corrosion medium during the constant load corrosion tests.[73] The model gives accurate predictions of the ultimate tensile strength and time-to-failure, although the validation was limited to relatively short term corrosion times (≤ 90 hours). However, the model could potentially overpredict long term mass loss as hydrogen gas evolution measurements were used to calibrate the model parameters. Amerinatanzi et al.[77] also proposed a response surface methodology-based approach to determine the above model parameters that relies on only using immersion mass loss measurements.

Stress corrosion

Gastaldi et al.[68] proposed the first computational model to simulate stress corrosion of Mg. The stress corrosion model also proposes the use of a scalar damage parameter (Dsc), as was utilised in the uniform and pitting corrosion models. Stress corrosion is only considered to occur when the equivalent stress (σeq*) exceeds the threshold stress (σth), where σth is 50% of the value of the yield strength. σeq* was arbitrarily assumed to be equal to the maximum principal stress. Changes in the damage field with time are described by equation (12) when σeq* > σth.[63]

where, S and R are experimentally-determined parameters related to the kinetics of the stress corrosion process.

The total damage is considered to be a linear superposition of the damage from both uniform and stress-corrosion. Wu et al.[78] applied the model to estimate the percentage reduction in the outer diameter of an expanded bioresorbable Mg alloy stent; decreases in the expanded diameter of the stent under the pressure of the surrounding artery is known as the recoil. Stress-induced corrosion in addition to uniform corrosion results in faster recoil of the stent according to the model, although the results were not validated experimentally. Oppeel[75] and Debusschere et al.[76] combined the model described by equation (12) with the pitting corrosion model proposed by Gorgan et al.[73] However, the validity of all of these models for the prediction of in vivo biodegradation of Mg alloys is yet to be determined.

Strain corrosion

Galvin et al.[42] proposed a strain-mediated uniform corrosion model that captures the corrosion behaviour of Mg-based coronary stents that have undergone severe plastic deformation. One unique feature of this model is that the corrosion rate parameter kc is a function of both time and strain which eliminates the need for an additional parameter to capture the acceleration or deceleration of corrosion with time for the uniform and pitting corrosion models described earlier. The magnitude of kc of an element at time t was calculated by multiplying the surface exposure parameter (ke) by the strain corrosion parameter (фe) calculated by interpolating between the data in (equation (13)). was obtained experimentally by conducting immersion corrosion tests on pre-strained specimens. The surface exposure parameter ke is another unique characteristic of the model. ke is a scalar derived from a linear approximation of the slope of the curve that describes the rate of diffusion of Mg2+ ions through a brick (or trilinear hexahedron) element where the number of exposed surfaces was varied from 0 to 6 and the mesh is based on hexahedral elements. фe was calculated from a cubic function dependent on the pre-strain (ε) if the strain exceeded 8.4% (equation (14)).

Table 1

фe as a function of pre-strain and immersion corrosion time.

Time (hour)	ε (%)
	0	2.7	8.4
0-24	0.368	2.190	5.830
24-48	1.287	1.265	1.264
48-96	1.000	0.632	0.360
96-144	1.035	1.092	1.149
144-192	0.862	1.149	1.093
192-240	0.875	1.095	1.148

Note: Data were from Galvin et al.[42]

Physical modelling

The chemical interactions that occur between Mg2+ and other ionic species result in the formation (and dissolution) of corrosion products of varying complexity depending on the composition of the corrosion medium.[50] The formation of surface films and their stability plays a crucial role in controlling the kinetics of corrosion of the underlying substrate.[79] The concentration of dissolved O2, inorganic ions and proteins will also influence the corrosion kinetics.[80] However, phenomenological modelling does not aim to capture these intimate chemical interactions, or the effects of alloy composition and microstructure on the corrosion kinetics. In contrast, physically modelling approaches such as activation-controlled corrosion (ACC) and transport-controlled corrosion (TCC) are able to address some of the limitations associated with phenomenological modelling.

Activation-controlled corrosion

ACC modelling assumes that the corrosion process is driven solely by the difference in electrochemical potential between the metal and surrounding physiological medium.[81] ACC models can capture the preferential dissolution of the α-Mg matrix phase based on chemical activity, allowing for the effect of microstructure on the rate of corrosion.[74] FEM is typically used to solve the Laplace equation (equation (15)) with suitable boundary conditions for a scalar field of potential difference (E).

Deshpande et al.[59, 81] adopted an ACC approach to model micro-galvanic corrosion of Mg alloys based on the proportion and distribution of secondary phases. The ACC approached used by Deshpande et al.[59, 81] is based on the unique non-linear polarization behaviour of each individual phase. Different polarisation behaviour of the individual phases results in preferential corrosion of one phase with respect to another. However, the model has not been validated against in vivo or in vitro corrosion studies. Montoya et al.[82] used the ACC approach to model the effect of electrolyte volume on the corrosion rate of an AZ31 implant in vivo. The model predicted that a higher volume of electrolyte would result in higher corrosion rate. The model was partially validated with an in vivo study. The in vivo studies involved intramedullar insertion of AZ31 alloy cylindrical specimens in the femurs of Wistar rats. The effect of presence and absence of fracture in a region surrounding the specimen on the rate of corrosion of implant was evaluated. The corrosion rate was observed to be higher in the specimen attached to the fractured femur, presumably due to increased body fluids surrounding the implant due to the fracture. However, the validation experiments were only qualitative as the corrosion rate measured in vivo was different compared to corrosion rate estimated from the model.

Transport-controlled corrosion

Scheiner et al.[83] demonstrated that the corrosion rate of Mg is governed by activation phenomenon only in the early stages of corrosion due to the formation of a passive or semi-passive surface film. In contrast, TCC models treat the corrosion of Mg as a diffusion problem. Thus, the corrosion kinetics are driven by the difference in Mg2+ ion concentration across the passive surface layer in TCC models.[83, 84] FEM has been used to solve the simplified diffusion equation in TCC models equation (16).

where, C is the concentration of Mg2+ ions, and D is the diffusivity of Mg2+ ions in the physiological medium.[83, 84]

Gorgan et al.[84] and Scheiner et al.[83] used this approach to model the corrosion rate of Mg stents. The velocity (v) of the moving boundary was calculated using equation (17).

where, Csol is the Mg2+ ion concentration in the solid Mg alloy being modelled (Csol is 1735 kg/m3 for pure Mg which is the density of pure Mg), and Csat is the saturation concentration of Mg2+ ions in the physiological medium.[84]

Dahms et al.[85] proposed an analytical model that assumes the long term corrosion rate in vitro to be constant, following the initial rapid corrosion rate equation (18). However, the model does not explicitly account for stress corrosion as the corrosion is assumed to be solely controlled by either an activation or transport phenomenon, as in other ACC- and TCC-based physical models.

where, h is the long-term corrosion rate, h0 is the initial corrosion depth, and h is the corrosion depth after time, t equation (19).

Hybrid modelling

Physical models assume that either the electrochemical potential of different phases or impurity-Mg coupling (activation-controlled models) or the difference in concentration of Mg2+ ions between alloy and surrounding fluid (transport-controlled models) govern the rate of corrosion. Physical corrosion models may attempt to capture the effect of microstructure[81] and Mg2+ ion concentration on the corrosion rate.[84, 86] However, physical corrosion models do not attempt to mimic the formation (or dissolution) of corrosion products due to the complexity and computational costs associated with modelling the physicochemical processes underlying the formation and dissolution of corrosion products.[86] Furthermore, physical corrosion models cannot directly capture the effect of coatings or surface roughness on the corrosion rate of Mg-alloy based biomaterials in contrast to phenomenological corrosion models. In contrast, phenomenological modelling can capture the effect of microstructure, surface coatings due to corrosion product or specific treatments, residual stress or stain, and surface roughness on the corrosion rate of Mg implant devices, provided the models are calibrated using test specimens manufactured to incorporate the above material characteristics. Therefore, a phenomenological model of the in vivo biodegradation of Mg should be highly accurate for a specific material state if the in vitro calibration experiments mimic the corrosion conditions in vivo. However, in vivo conditions are challenging to replicate in a simple in vitro experiment due to the level of complexity and dynamic control of the physiological environment at the implantation site. Further, the application of phenomenological modelling is time consuming as the experimental calibration needs to be repeated whenever there is a change in material characteristics (e.g., alloy composition, microstructure, surface finish, etc.).

Hybrid modelling aims to overcome the above limitations of each modelling method by combining physical and phenomenological models for the prediction of corrosion rate. For example, Bajger et al.[86] modified a physical corrosion model to incorporate the effect of surface passivation on the corrosion rate of Mg device-scale specimens equation (20). Hybrid corrosion models neglect the formation of hydroxyapatite or protein layers that are deposited on the substrate surface in vitro or in vivo, similar to other physical corrosion models. Thus, two additional parameters, K and K, are included to denote the formation and dissolution rates of the passivation layer, respectively, to indirectly consider these effects. The model proposed by Bajger et al.[86] is based on the fact that the diffusion of Cl- into the corrosion film results in the conversion of Mg(OH)2 to Mg(Cl)2, a more soluble corrosion product.[87] The rate of change of the total ion concentration (i.e., [Mg2+] and [Cl-]) within the corrosion layer is governed by equation (21), while the rate of diffusion of Cl- into the corrosion layer is given by equation (22). Essentially, the rate of formation of the corrosion film in this model is driven by the difference between the saturation concentration of Mg2+ and Cl- and the concentration of Mg2+ and Cl- in the corrosion film at a given time. The rate of dissolution of the corrosion film is driven by the difference between the saturation concentration of Mg2+ in the surrounding fluid and concentration of Mg2+ in the corrosion film.[86]

where, C is the concentration of Cl-, F is the concentration of ions in the corrosion film, F is the saturation ions in the corrosion film, and D is diffusivity of Cl- within the corrosion film.

Similarly, Shen et al.[88] incorporate surface passivation into an existing physical model with the inclusion of a Mg(OH)2 film on a Mg substrate. Consequently, the model modifies the velocity (v) with which the outer surface of the model moves inwards to mimic mass loss equation (23). An additional parameter is used to account for the possible presence of porosity and/or other corrosion products. ρMg(OH) and MMg(OH) are density and molecular weight of Mg hydroxide.

where, Ccor is the concentration of Mg2+ ions in the corrosion layer that can be determined from equation (24).

Other corrosion modelling techniques

Artificial neural networks (ANNs) are an eficient numerical modelling tool, inspired by the working principle of the human brain, that is also based on a neuronal system.[89, 90] ANNs can be used to perform computing simulations and to characterize corrosion processes based on experimental observations.[91] ANNs can be designed and trained to estimate corrosion rates of metallic materials, from a set of environmental and material parameters. Material parameters would mainly include chemical composition and microstructural features, while the environmental parameters could include pH, chemical composition, temperature, etc. The ANN algorithm trains itself using the above information to predict corrosion rate of an alloy based on the experimentally determined corrosion rate. Kirkland et al.[92] and Birbilis et al.[93] used ANNs to predict the effect of chemical composition on the mechanical strength and corrosion rate of Mg alloys. Willumeit et al.[94] used ANN to predict the effect of the composition of the simulated body fluid (SBF) on the in vitro corrosion rate of Mg-based alloys in SBF. Xia et al.[95] proposed a model based on ANNs for predicting the corrosion rate and hardness of Mg alloys based on their chemical composition. The ANN approach is promising due to its ability to estimate the synergistic influence of several factors on the corrosion rate of metals. However, a major limitation of the ANN-based corrosion modelling approach is the large amount of data required to train the model before being able to be reliably deployed.

Finally, new corrosion modelling approaches still under development have been based on fuzzy logic-based methods[96, 97] and microstructural parameter-based modelling.[98] These methods have been proposed for the prediction of the corrosion rate of Mg-based alloys and other biometals as a function of their chemical composition and microstructure. However, these models have not been fully validated due to the large amount of experimental data necessary to train these models.

Experimental Calibration of Corrosion Models

Introduction

The determination of the various parameters that appear in the above phenomenological and hybrid corrosion models usually requires a number of experimental measurements - a process also known as model calibration. The inherent limitation of phenomenological, physical and hybrid models is their reliance on replicating in vivo conditions using in vitro experiments to calibrate these models. Potentially, the models could be directly calibrated using in vivo studies. However, the calibration of a model needs to be performed whenever the material is changed, making in vivo calibration impractical from both ethical and financial standpoints.

Phenomenological models are all calibrated in a fundamentally similar approach. Firstly, the mass loss is experimentally determined at several time points using immersion corrosion testing or hydrogen evolution testing. The model is then calibrated by identifying a set of parameters that minimize the root mean squared difference between estimated mass loss and predicted mass loss by the model for a particular set of parameters. Once calibrated the model should ideally be validated using in vivo techniques by comparing the mass loss in vivo and mass loss predicted by the model.

Both immersion mass loss and hydrogen evolution are frequently utilised for calibration of phenomenological models. However, the use of electrochemical test methods for the calibration of corrosion models is not commonly reported due to the indirect nature of the mass loss measurement.[74] Theoretically, in vitro corrosion tests including potentiodynamic polarisation (PDP), electrochemical impedance spectroscopy, and pH or Mg2+ ion concentration monitoring can potentially be used to calibrate phenomenological models, although there are some complicating factors. For instance, the mass loss predicted from PDP assumes that the rate of corrosion remains constant. However, the corrosion rate estimated from PDP can deviate from that observed with in vitro or vivo mass loss experiments. The reason is that there is normally a decrease in corrosion rate with the greater formation of corrosion product in longer term tests. Thus, the corrosion rates measured by short term PDP or electrochemical impedance spectroscopy testing may overestimate the long-term corrosion rate predictions for Mg alloys. In reality, the mass loss directly measured during in vitro tests tends to decrease over time due various factors including formation and build-up of corrosion products, hydroxyapatite precipitation, protein adsorption, changes in pH, and buffer selection. Additionally, there is bone remodelling, and other physiological processes that occur in vivo that can decrease the corrosion rate over time.[16]

Immersion corrosion testing

Immersion corrosion testing of metals (aka mass loss testing) is a standard method for directly measuring mass loss to determine the general corrosion rates of metals. The experimental conditions that may be controlled include test duration, temperature, pH, oxygen concentration, flow rate, and medium composition. Mass loss testing is frequently used to provide a quantitative measure of the rate of corrosion of biometals when immersed in a physiologically-relevant corrosion medium (e.g., SBF). The ratio of the mass loss to surface area is measured as a function of the immersion time to estimate the corrosion rate. However, the measurement of mass loss alone does not indicate the corrosion mechanism, and normally microscopy is needed to determine the nature of the corrosion (e.g., uniform, pitting, etc.). A variety of corrosion media are used for mass loss testing of Mg alloys including sodium chloride (saline) solution, phosphate-buffered saline, Ringer’s solution, Earle’s balanced salt solution, Hank’s balanced salt solution, conventional SBF, cell culture media (e.g., minimum essential medium, Dulbecco’s modified eagle medium, etc.), and cell culture media containing proteins (e.g., fetal bovine serum, bovine serum albumin, etc.).[99] The pH and temperature of the SBF should be maintained at pH7.2 and 37°C, respectively. Mass loss tests are often carried out in accordance with the standard ASTM G31-21.[100] The weight of the Mg specimen is recorded prior to immersion in the corrosion medium. Post-testing the Mg alloy specimen is acid-cleaned to remove corrosion products prior to reweighing.[101] Any model can be calibrated using immersion testing by identifying a set of parameters that minimise the root mean squared difference between mass loss measured from experiments and that predicted by the model. However, the composition of the corrosion medium influences the in vitro corrosion rate, making it difficult to compare test data obtained via different corrosion media.[99] Thus, it is critically important to choose a corrosion media that replicates the in vivo conditions as closely as possible such that the calibrated model can reliably predict the in vivo corrosion rate.

Stress- and strain-mediated immersion corrosion testing

The introduction of applied stresses or pre-straining of specimens to be subjected to immersion corrosion testing adds yet another approach that has been used for the calibration of phenomenological (or hybrid) stress or strain corrosion models. For example, Grogan et al.[73] utilised a constant load corrosion test to validate a pitting corrosion model. A constant tensile stress (10-180 MPa) was imposed on the corrosion specimens during immersion testing. The time-to-failure and mass loss were both compared to the predicted values from the model. The model was able to accurately predict the non-linear reduction in time-to-failure with increasing load. However, orthopaedic devices are typically subjected to a combination of tensile, compressive and shear loads during both the surgical implantation procedure and in situ as a bone fixation device. Further work is required to assess the calibration of phenomenological (or hybrid) corrosion models under different modes and combinations of loading.

Galvin et al.[42] used tensile pre-straining of immersion specimens to calibrate a strain-mediated corrosion model. Specimens were subjected to tensile plastic pre-strains of up to 8.4%.[42] The calibration procedure was similar to that adopted by Grogan et al.[74] in terms of identifying the model parameters. The calibrated model was able to predict an increase in the corrosion rate with an increase in plastic strain. However, once again only the effect of tensile strain on the corrosion rate was considered, although orthopaedic devices are subject to a combination of tensile, compressive and shear loads. Further work is needed to modify the calibration approach to allow for the adjustment of calibration parameters depending on both the magnitude and combination of plastic strain modes.

Hydrogen gas evolution

The dissolution of Mg during the corrosion reaction in an aqueous medium result in the evolution of a stoichiometric equivalent amount of hydrogen gas (H2(gas)) equation (3).[102, 103] Consequently, the corrosion rate of Mg can be determined indirectly by collecting and measuring the volume of evolved H2(gas) with the use of a funnel and burette equation (25).

where VH is the volume of the H2(g), and Patm is atmospheric pressure.

The main limitations associated with the measurement of hydrogen evolution are difficulties with (i) capturing 100% of the evolved H2(gas), and (ii) maintaining a constant physiological pH with a sodium bicarbonate/CO2 buffer within a closed system.[74] Consequently, hydrogen evolution is not frequently reported as a method for calibration of phenomenological (or hybrid) corrosion models. Gorgan et al.[73] determined the corrosion rate of Mg AZ31 alloy specimens in Hank’s balanced salt solution from the volume of evolved H2(gas) to calibrate the corrosion model equation (25). The values of the parameters for uniform (k, β0) and pitting (k, β0, γ, ψ) corrosion models were then estimated by minimising the root mean squared difference between the mass losses determined from hydrogen evolution and that predicted by the models. If the corrosion is relatively uniform, the values of parameters would be, ψ ≌ 1 and γ << 1 which essentially result in value of parameter αe for each element to be almost 1 as in case of uniform corrosion model. Conversely, if the corrosion is highly localised, uniform corrosion model would fail to account for varying corrosion rates across various elements and thus it is unlikely to fit the data.

Validation of Corrosion Models

In vivo testing (i.e., animal models and clinical testing) is critical to understanding the longer term biological response to Mg-based implants, including bone healing/remodelling, bone-implant compatibility, bone-implant interfacial properties, inflammatory responses, etc.[104] However, there are ethical concerns and tremendous costs associated with carrying out in vivo testing. Thus, it is perhaps not surprising that the majority of mathematical models described in the literature for the prediction of Mg corrosion are not yet validated by in vivo test results.[82] In vitro corrosion testing provides a range of systematic test methods that permit a deeper understanding of how corrosion rates (and associated changes in mechanical properties) of Mg-based biomaterials are influenced by the properties of the material and corrosion medium. The majority of computer simulation-based corrosion modelling methods has the objective of correlating findings from in vitro testing with that observed from in vivo testing. However, it has been widely demonstrated that the in vitro and vivo corrosion rates of Mg differs significantly due to the complexity of mimicking the physiological environment.[17] Thus, a robust validation of any corrosion model still mandates the need for in vivo test results. In the following sections, a range of measurement techniques that have the capacity to evaluate biodegradation rates in vivo is presented, providing a possible framework for future experimental validation of mathematical models that aim to predict the biodegradation rates of Mg and its alloys in vivo. Micro-computed tomography (micro-CT or µCT) scanning is a chief method for the measurement of mass loss in Mg-based orthopaedic implants in vivo.[50] However, a range of in vivo techniques are being developed that could be prove to be valuable for validation of corrosion models, including µCT, electrochemical monitoring, and microsensor measurements.[50]

Micro-computed tomography

Micro-CT uses X-rays to illuminate a rotating sample, capturing several 2D images from different angles. The 2D images are reconstructed to create a high-resolution 3D image of the structure, providing volumetric information. Changes in the volume of the object of interest over time can then be related to the in vivo corrosion rate. Several researchers have used µCT to quantify the rate of corrosion in vivo.[105-108] µCT is a non-destructive technique with high temporal and spatial resolution that can be used to measure the biodegradation of the implant device, while also monitoring the bone remodelling in situ.[109-117] The measurement of the volumetric (or mass) loss with µCT at various timepoints could be used to validate the predicted mass loss from the corrosion model. The resolution of µCT (typically ∼25 µm) is sufficient to accurately measure mass loss of the implant. However, there are certain limitations associated with µCT including the (i) use of ionising radiation, which restricts the number of allowable timepoints for in situ monitoring of the implant in vivo, and (ii) directly determining changes in the mechanical integrity of the implant.

Electrochemical-based monitoring of biodegradation

Doepke et al.[118] developed a system that permits real-time monitoring of OH-, Mg2+ and H2(gas) concentrations in the corrosion medium during an in vitro immersion test. The rate of mass loss (or corrosion rate) is obtained by measuring the concentration of H2(gas) or Mg2+ in the corrosion medium. Although the system was originally developed for in vitro testing, the system could be potentially adapted to monitor in vivo corrosion using a vascular bioreactor.[119] Unfortunately, the method is not well suited to continuous, long-term monitoring due to limitations on the electrode life (∼minutes to days) and the invasive nature of the technique, limiting the total number of measurement timepoints.[50] Electrochemical-based monitoring techniques provide insight into the rate of corrosion and gas formation, although the measurements are not fully quantitative as unlike similar in vitro techniques it is impossible to definitively quantify the H2(gas) or Mg2+ ions released during the corrosion process.

Microsensor-based monitoring of biodegradation

The H2 concentration in vivo can be measured in the immediate vicinity of the implant at different time points using various microsensor-based techniques that are described below. The corrosion rate at different timepoints is estimated from the measured H2 gas concentration, which will be invaluable for the validation of the various proposed corrosion models.

Zhao et al.[120, 121] developed highly sensitive electrochemical microsensor mounted inside a 16-gauge needle to measure H2(gas) evolution in vivo. A strong correlation was reported between the measured mass loss of the ex-implant after 4 weeks and concentration of H2(gas) measured after 1 week of implantation. However, this method does not allow for continuous monitoring of H2(gas) the needle cannot be left permanently in position in an animal model. Thus, several measurements at different timepoints are required to build up an understanding of the long-term corrosion rate in vivo. One major advantage of needle-mounted sensors is the high-resolution measurement of H2 at multiple locations surrounding the implant providing more detailed spatially-sensitive information on the corrosion rate. Visual biosensors have also been recently developed for the transdermal monitoring of the H2(gas) concentration in vivo.[122] These visual sensor consisted of a thin film of H2 sensitive material coated on a flexible plastic sheet that was pressed against the mouse skin directly above the implant. A major advantage of this method is the non-invasive nature of the measurement. However, the diffusivity of H2(gas) via the skin is relatively low, rendering this method only effective for corrosion monitoring of subcutaneous implants. The result obtained from visual biosensors used with subcutaneous implants are comparable to that from needle-mounted microsensors.[122] Boutry et al.[123] propose the development of biodegradable wireless implantable sensors to monitor the corrosion of bioresorbable implants. Wireless sensors will allow non-invasive, continuous monitoring of the corrosion rate, although these devices are still very much in the development phase.

Microdialysis-based monitoring of biodegradation

Microdialysis is a minimally-invasive sampling technique used for continuous measurement of analyte concentrations (e.g., Mg2+) in any tissue. Microdialysis-based monitoring techniques have been coupled with plasma mass spectrometry and a micro-flow capillary injection system to continuously monitor concentration of Mg2+ and H2(gas) for estimating the corrosion rate.[124, 125] However, microdialysis systems are still in the development phase due to a limited probe life, potential risk of biofouling, and qualitative nature of the data. Thus, further development is required before this technique can validate the results obtained from the computational modelling techniques.

Future Outlook for Mathematical Modelling, Calibration and Validation

A summary of the various modelling techniques, along with calibration and validation methods, is presented in . Corrosion modelling along with its calibration and validation is a rather complex process. This section highlights the future outlook for developing future corrosion models, calibration methods and validation techniques.

Figure 1

An overview of the potential of currently-available experimental in vitro and in vivo techniques for the calibration and validation of mathematical corrosion models. µCT: micro-computed tomography; H2(g): hydrogen gas.

Mathematical modelling

There are several limitations that the current modelling approaches need to address. For instance, Witte et al.[126] demonstrate that Mg-based implants can enhance osteoblastic response used µCT scanning. Few groups have attempted to model effect of bone remodelling on mechanical performance of the biodegradable Mg-based implants.[127, 128] However, no such model has been validated. No research group has attempted to model effect of bone remodelling on corrosion rate (or vice versa). Ideally, the corrosion model should predict the effect of bioresorbable implants on osteoblastic response and tissue differentiation and in turn predict effect of the osteoblastic response and tissue differentiation on the mechanochemical performance of the implant. Additionally, it is important to ensure that bone interface has regained sufficient strength before implant loses its mechanical integrity.[14] However none of the current modelling approach is focused on modelling the bone healing process alongside implant corrosion modelling. It is also important to predict potential inflammatory response through modelling in addition to information about mass loss.[129, 130] No current model predicts the potential inflammatory response due to corrosion. Additionally, no current model predicts rate of accumulation or dissipation of metallic ions from surrounding tissue which is important to understand the risk of potential inflammatory reaction.[131]

An ideal corrosion model should combine different features and address the limitations of various modelling approaches discussed above. An ideal corrosion model should be able to predict mass loss and loss in mechanical integrity of the implant alongside predicting the changes in bone interface strength throughout the process of bone healing. The model should also be able to predict gas formation, potential risk of inflammatory reactions and their impact on bone healing. Additionally, the model should be able to predict the potential risk of crevice corrosion and environmentally assisted corrosion. The model should be easy to calibrate, accurate and repeatable. An ideal corrosion model should also be extensively validated using in vivo experiments.

Model calibration

The corrosion models are calibrated using in vitro corrosion monitoring techniques. The corrosion rate observed in vitro strongly depends on the composition of the SBF.[14, 132] Currently, it is impossible to replicate the exact in vivo conditions while performing in vitro corrosion experiments due to several factors such as risk of biofouling and the dynamic nature of the composition of the physiological environment. Thus, further development is necessary to develop in vitro corrosion monitoring techniques that can replicate in vivo conditions as closely as possible. The calibration process in itself should be further developed to indirectly account for various factors that would influence the corrosion rate in vivo but impossible to replicate in vitro such as, dynamic corrosion environment, bone remodelling process, etc.

Model validation

In vitro studies indicate that pitting corrosion can cause a non-linear loss of mechanical integrity.[73] The non-linear loss of the mechanical properties of implants is able to be captured by pitting corrosion models.[73] Loss of mechanical integrity implant can lead to premature implant failure. Thus, it is essential to quantify the loss in the mechanical integrity of the implant in vivo. However, current in vivo corrosion monitoring techniques are limited in their ability to quantify the loss of mechanical strength of implant predicted by pitting corrosion models. Ultrasound velocity and attenuation measurement across the fracture interface is a promising method to measure the bone interface strength.[133] However, as of now this method has only been successfully demonstrated to measure strength across fractures in long bones[133] Additionally, in vitro studies also indicate that, the applied stress and induced strain can enhance the rate of corrosion. However, some recent in vivo studies show that stress do not alter in vivo corrosion behaviour.[134] Nonetheless, further studies are essential to quantify the effect of stress and strain on corrosion rate in vivo. However, no current in vivo corrosion monitoring technique can be used to quantify the applied stress and induced strain in the implant and their impact on corrosion rate and hence further development is essential.

Currently, the state-of-the-art technique for assessing mass loss is µCT by precisely monitoring volumetric changes in situ. However, µCT is expensive and monitoring of Mg biodegradation with µCT is restricted by the number of allowable time points.[50] Several in vivo corrosion monitoring techniques are still in the early stages of development, although in the future it may become possible to utilise a combination of the above discussed techniques to validate mathematical corrosion models. An ideal validation method should be able to accurately measure the mass loss of the implant throughout the process of bone healing and simultaneously quantify the loss of the mechanical integrity of the implant based on a non-invasive approach.

Conclusions

Phenomenological corrosion models are relatively easy to develop and can be deployed to predict the rate of corrosion when an implant is subjected to plastic strain or stress. However, phenomenological corrosion models are specific to the alloy under investigation. Moreover, separate calibrations may be required for the same alloy if there are changes in the corrosion environment (e.g., electrolyte composition, local pH, temperature, etc.).

Physical corrosion models are capable of modelling the effect of local electrolyte composition and microstructure on corrosion performance. The TCC models assume that the corrosion process is driven by the difference in concentration of ions across the partially passive layer of corrosion products. This assumption limits the practical application of these models since the corrosion process is governed by several factors including alloy composition and microstructure, presence of impurities, surface passivation, etc. ACC models can be used to predict the effect of microstructure on the rate of corrosion. However, the corrosion of an implant in vivo is not solely governed by the activation phenomenon. Thus, ACC models do not mimic the long-term changes in the corrosion process. Hybrid corrosion models can provide benefits associated with both phenomenological and physical corrosion models while overcoming several limitations. These models can consider the effects of the passivating surface films that form on the corrosion rate. However, these models are difficult to calibrate due to a large number of parameters. Additionally, hybrid models need further modification to predict stress-mediated (or strain), filiform, or crevice corrosion.

A limitation associated with the majority of models reviewed is that they need to be calibrated for a specific alloy composition, alloy microstructure and surface finish. Furthermore, the models are not able to simulate filiform or crevice corrosion following calibration of a specific alloy. ANN-based corrosion models have the potential to eliminate the need for time-consuming calibration experiments. However, ANN-based models mandate a large number of experimental results to train the models for reliability.

Amongst models proposed thus far, phenomenological models calibrated using in vitro corrosion experiments performed under the appropriate conditions for the particular application would provide the most accurate estimation of in vivo corrosion behaviour. However, there is major scope for improvement in these models by adopting features of other modelling approaches.