A Predictive Toxicokinetic Model for Nickel Leaching from Vascular Stents

In vitro testing methods offer valuable insights into the corrosion vulnerability of metal implants and enable prompt comparison between devices. However, they fall short in predicting the extent of leaching and the biodistribution of implant byproducts under in vivo conditions. Physiologically based toxicokinetic (PBTK) models are capable of quantitatively establishing such correlations and therefore provide a powerful tool in advancing nonclinical methods to test medical implants and assess patient exposure to implant debris. In this study, we present a multicompartment PBTK model and a simulation engine for toxicological risk assessment of vascular stents. The mathematical model consists of a detailed set of constitutive equations that describe the transfer of nickel ions from the device to peri-implant tissue and circulation and the nickel mass exchange between blood and the various tissues/organs and excreta. Model parameterization was performed using (1) in-house-produced data from immersion testing to compute the device-specific diffusion parameters and (2) full-scale animal in situ implantation studies to extract the mammalian-specific biokinetic functions that characterize the time-dependent biodistribution of the released ions. The PBTK model was put to the test using a simulation engine to estimate the concentration–time profiles, along with confidence intervals through probabilistic Monte Carlo, of nickel ions leaching from the implanted devices and determine if permissible exposure limits are exceeded. The model-derived output demonstrated prognostic conformity with reported experimental data, indicating that it may provide the basis for the broader use of modeling and simulation tools to guide the optimal design of implantable devices in compliance with exposure limits and other regulatory requirements.


INTRODUCTION
In recent decades, there has been a significant surge in the use of implantable devices driven by not only a growing clinical demand but also rapidly evolving technologies in composite and synthetic materials science. 1,2The increasing incidence of chronic illnesses and the expanding elderly demographic are the primary drivers of the global medical implant market, 2 and thus, continual research and development are required to improve their performance and extend their useful life span. 3tents are cylindrically shaped, hollow mesh metal structures used to restore the flow in blood vessels and ducts. 4−7 Although this passivating film insulates the bulk material from the corrosive physiologic fluids, any mechanically induced break or defect in it increases the risk of corrosion. 5While metal alloys might endure significant forward and reverse mechanical strain excursions, the presence of a nonsuperelastic oxide layer could lead to cracking.This creates a pathway for the exposure of metal ion-rich phases to the in vivo environment, ultimately causing elevated ion release and reduced resistance to pitting. 5,8Experimental and clinical data, collectively reported by the U.S. Food and Drug Administration (FDA), suggest that metal implants experience different degrees of wear and corrosion due to the mechanical and biochemical environment at the specific site of implantation. 9−20 The majority of metal alloys used for stent manufacturing contain high levels of nickel, which is metallurgically necessary to impart enhanced mechanical properties but may lead to several adverse health effects if leached in high concentrations. 9o assess general corrosion susceptibility and metal ion release, the FDA recommendations include testing per the American Standard of Testing and Materials (ASTM) F2129-08 (Standard Test Method for Conducting Cyclic Potentiodynamic Polarization Measurements to Determine the Corrosion Susceptibility of Small Implant Devices), 21 the ASTM F3306-19 (Standard Test Method for Ion Release Evaluation of Medical Implants), 22 and the ASTM G31-72(2004) (Standard Practice for Laboratory Immersion Corrosion Testing of Metals). 23−26 Specifically, Sullivan et al. 24 studied the effect of oxide layer composition in nitinol stents manufactured by various processing methods and suggested strong correlations between the layer thickness and cumulative Ni release.Nagaraja et al. 26 tested generic heart valve frames manufactured with different vacuum arc remelting methods and surface finishes and showed that the localized and uniform corrosion performance is maintained in higher microstructural nitinol purity.
In addition to providing a framework for comparing different alloys, designs, or manufacturing processes, this type of testing and data are used to estimate exposure as part of toxicological risk assessment per ISO 10993 (Biological Evaluation of Medical Devices). 27Despite the valuable insights that in vitro test methods can offer regarding the corrosion susceptibility of a particular device, these tests are usually carried out under idealized or hyperphysiological conditions.Consequently, even though these tests facilitate easy comparisons between devices, the degree to which in vitro performance aligns with in vivo corrosion behavior remains uncertain.Interestingly, Nagaraja et al. 28 highlighted the impact of fatigue loading on uniform corrosion for different nitinol stent surfaces, while Sussman et al. 29 demonstrated increased Ni release from nitinol devices exposed to different physiological environments characterized by the level of pH and reactive oxygen species (ROS).These findings emphasize the significance of taking the implantation site into account when designing studies to predict nickel release from medical implants and underscore the importance of the biomechanochemical environment at the device−tissue interface.
The FDA guidance document (Technical Considerations for Non-Clinical Assessment of Medical Devices Containing Nitinol Guidance for Industry and Food and Drug Administration Staff, issued on July 9, 2021) 30 recommends that testing should be performed on as-manufactured devices and in simulated-use conditions.This type of testing should closely mimic physiologic conditions and capture both the initial bolus release of the substance and the longer-term release profile in vitro.Animal studies are often used to evaluate biocorrosion in a relevant anatomical site under simulated clinical conditions. 31However, biocompatibility assessments are not typically designed to evaluate the biological response to mechanical failure, and in some cases, anatomical differences between humans and animal models limit quantitative correlation.Additional studies are required to evaluate the in vivo response to failure modes such as coating defects and the release of wear particles and ions.Although linking the results of in vitro testing to in vivo outcomes is a formidable challenge, available data do suggest that there are at least qualitative consistencies.Data following in situ implantation of stents are scarce in the literature, yet Nagaraja et al. 32 presented one of the very few studies that investigated the relationship between nitinol surface processing, in vivo nickel release, and biocompatibility.Nitinol stents manufactured using different surface treatments, and previously characterized for their in vitro release profile, 24 were evaluated in a porcine stent implantation model.The investigation focused on the concentration of ions in local vascular tissue and in serum and urine samples and highlighted that those stents with nonoptimized surface finishing had significantly greater nickel leaching and gave rise to adverse inflammatory reactions and restenosis compared to polished stents.
While there is qualitative consistency between engineering testing and the behavior inside the body among devices, the quantification of these relationships remains a challenge for corrosion prediction within individual patients.The FDA suggests that a relatively easy and promising approach to obtain such quantitative relationships is the use of modeling and simulation tools. 9,30,33Physiologically based toxicokinetic (PBTK) models bear the capability to connect in vivo ion release with clinically measurable parameters such as ion concentrations in blood or urine.This allows for a direct comparison between the inferred in vivo exposure from these measurements and the outcomes of in vitro testing. 34stablishing such correlations necessitates a combination of computational modeling and in vitro, ex vivo, and in vivo testing and eventually clinical studies to parameterize and validate the model predictions.Although a number of studies have reported data regarding the concentration of different metal ions in serum and urine for healthy adults, 35−38 there is still limited information about the levels of released ions, especially their accumulation in tissues/organs following the in situ implantation of medical devices. 39,40Biodistribution studies in animals are mostly after the injection 41,42 or oral administration 43,44 of metal ions, while the collection of toxicological data from implanted patients has many challenges associated with data heterogenicity and nonroutine follow-up examinations.
Saylor et al. 45 proposed a biokinetic model that estimates nickel release from an implanted device.However, due to the lack of coherent data, the model was parameterized using immersion corrosion testing of nitinol stents 24 and values following implantation of atrium occluders in humans. 46A recent work by some of us 47 presented an expanded and enriched version of that model by adding a separate organ (kidney) compartment to better resemble normal physiology and by introducing time-dependent functions to describe some of the biokinetic parameters.The upgraded model was exercised in conjunction with probabilistic Monte Carlo simulations and exhibited quantitative consistency with nickel levels in serum and urine following the implantation of atrium occluders in humans 46 and data reported by Nagaraja et al. 32 on stent implantation in minipigs.
The present work is, to the best of our knowledge, the first fully comprehensive study, which comprises device in vitro and in vivo testing, along with the formulation, parameterization, and testing of a multicompartment PBTK model and a simulation engine for toxicological risk assessment of vascular stents.The mathematical model describes the transfer of nickel ions from the device to adjacent tissue and circulation and the exchange between blood and the various tissues/organs (see Figure 1 and refer to Section 2.5.2).It is based on a detailed set of transport equations with the diffusion, absorption, distribution, and excretion variables as well as the appropriate initial and boundary conditions.Model parameterization was performed using (1) in-house-produced data from immersion testing to compute the device-specific diffusion parameters and (2) full-scale animal in situ implantation studies to extract the mammalian-specific biokinetic functions that characterize the time-dependent biodistribution of the released ions in tissues/ organs, body fluids, and excreta.The developed in silico tool was then put to the test using a probabilistic Monte Carlo methodology to estimate the concentration−time profiles, along with confidence intervals, of nickel ions leaching from the implanted devices and determine if permissible exposure limits are exceeded.

Stent Manufacturing and Processing.
Custom-made selfexpanding nitinol stents (Admedes GmbH, Pforzheim, Germany), 0.7 mm × 3.3 mm in dimension, were used throughout the study.Two closed-cell designs with diamond-shaped patterns V1 (version 1� double cell) and V2 (version 2�single cell) (shown in Figure 2a−c) were fabricated to evaluate the effect of stress raisers (connective links/struts) on Ni ion release.After laser cutting, mechanical polishing was applied to eliminate the heat-affected zone from all stent samples.Subsequently, the stents were categorized into two groups based on material surface condition and processing steps�heat treatment (HT) and electropolishing (EP)�with a resulting strut thickness of approximately 40 and 20 μm, respectively.Heat treatment was employed to modify the metal surface chemistry and topography by creating a thicker titanium oxide layer and affecting the nickel-rich zone, thus acting as a source of high nickel release to  mimic active in vivo corrosion.Electropolishing is the standard surface finishing method and was used to simulate a reduced level of ion leaching.These two groups cover a broad spectrum of potential manufacturing processes and help define a range of feasible release model parameters.Overall, we have studied three different devices: (1) the double-cell design, heat-treated stent (V1HT), (2) the double-cell design, electropolished stent (V1EP), and (3) the singlecell design, electropolished stent (V2EP).
2.2.Immersion Testing.All stent samples were weighted (n = 3 per group; 5-fold), and the exact values of their resulting volume and active surface area A surf (per stent design/surface treatment method) were calculated based on the nitinol alloy density value (6.5 g/cm 3 ) and three-dimensional (3D) model drawings (provided by the manufacturer).
Immersion testing was conducted on stent specimens (n = 3 per group; V1HT, V1EP, and V2EP), as per ASTM G31:72(2004) 23 and ASTM F3306-19, 22 to study the device-specific cumulative release profile data.All containers and instruments used for the testing, handling, and storage of specimens were acid-washed using a 10% HNO 3 solution (100441 Supelco; Merck KGaA, Darmstadt, Germany) before use.Due to the small surface area of each device (∼11−12 mm 2 ) and the low nickel concentrations in typical test solutions, sample pooling was performed for more accurate estimates.In each case, all 3 stents per group were placed in 15 mL polypropylene tubes filled with 10 mL of phosphate-buffered saline (PBS) (pH 7.4 ± 0.1), resulting in a PBS-volume-to-stent-surface-area ratio of 0.27−0.35mL/mm 2 , depending on A surf value.Although individual sample variability could not be analyzed with this method, a sufficient metal ion mass release (above the nickel detection limit of the analysis system) was required. 25During a 60-day incubation period at 37 °C, sampling was conducted at 10 time points (days 1, 2, 3, 5, 7, 14, 21, 28, 42, and 60) for nickel-ion measurements.At each time point, the total PBS volume was withdrawn from the test container, which was first flushed and then refilled with fresh solution.
The concentration of nickel ions was measured using highresolution inductively coupled plasma mass spectrometry (ICP-MS) (XSERIES2, Elemental Scientific Inc., Omaha, NE).To monitor metallic contamination throughout all steps, nickel levels were also measured in n = 3 control blanks (PBS in tubes with no stent specimen).Samples were acidified with 2% HNO 3 before analysis to ensure stability and comparability with calibration standards.The ICP-MS was tuned to the manufacturer's recommendations and calibrated using blanks.Known standards were analyzed to create a calibration curve with a linearity greater than 0.999.Raw measurements were obtained in units of ppb (ng/mL).The three stent types (V1HT, V1EP, and V2EP) were tested in vivo to establish the time-dependent concentration of the released Ni ions in tissues, organs, and body fluids.Drug treatment and surgical procedures were carried out as per the reported experimental protocol by Kapnisis et al., 48 which describes in detail a method for stent implantation in mouse carotid artery.In brief, following general anesthesia, a small median incision is performed at the ventral neck area, and the left common carotid artery (LCCA) is exposed.The internal carotid artery (ICA) and LCCA are controlled, distally and proximally, respectively, by slings, and the external carotid artery (ECA) is ligated distally.Following blood flow interruption, a small incision is performed on the ECA and a polymeric catheter containing the stent is guided into the desired position in the LCCA.After stent deployment, with a stent-to-artery-size ratio range of 1.1−1.2, the wound is closed and the animals are allowed to recover.

Nickel Biodistribution Analysis.
Stented mice (n = 3 per stent type and time point) were euthanized under anesthesia at 4, 8, and 18 weeks (±1 day at each time point) after surgery.Harvested tissue/organs (stented arteries, kidneys, liver, lungs, brain, and small intestine) and whole blood, urine, and fecal samples were collected for nickel-ion quantification through ICP-MS.ICP-MS is the most commonly employed technique for detecting metal ions in biological samples and is the preferred analytical method due to its high sensitivity. 49ascular tissue surrounding explanted stents was removed through a digestion process in a 1 M solution of NaOH, according to the protocol described in Kapnisis et al. 48Whole blood was collected via retro-orbital bleeding, and untainted urine specimens were collected via terminal direct bladder puncture using a glass Pasteur pipette.To avoid issues and variations arising from handling errors, the tissuedigested solution, the harvested organs, and the body fluid samples were collected and pooled (n = 3 per stent type and time point) for more accurate estimates.Tissue/organs and body fluid samples were also collected from healthy animals (n = 5 pooled) to establish baseline concentrations of nickel.All tools and containers utilized in both the tissue removal process and the handling of the explants were nonmetallic and acid-washed using a 10% HNO 3 solution prior to use.Nickel levels were also estimated in mouse food and water to calculate the average daily dietary uptake.Biological specimens were digested using a closed-vessel microwave digestion procedure before analysis via ICP-MS as described in Section 2.2.

Stent Surface Evaluation.
Following immersion testing and peri-implant tissue dissolution (for explanted devices), each stent sample was transferred into 1.5 mL Eppendorf tubes filled with ultrapure deionized water for 20 min of ultrasonication (twice) and allowed to air-dry overnight before surface evaluation.Scanning electron microscopy (SEM) analysis was performed prior to and post testing (in vitro and in vivo) using a Quanta 200 scanning electron microscope (FEI, Hillsboro, OR) to characterize and compare the surface topography and features of immersed and explanted to asreceived nontested (control) stents.The abluminal (outer), laser-cut side wall, and luminal (inner) surfaces along the entire length of each sample were inspected for any signs of wear and/or corrosion.Elemental analysis, with energy-dispersive X-ray spectroscopy (EDS), was also performed (prior to and after testing) on the outer stent surfaces to test for altered elemental composition and quantify nickelion release upon active corrosion.The weight % ratio of nickel to titanium (Ni/Ti) was derived from the EDS spectra from tested and control stents of the same surface treatment group.
2.5.Model Formulation.2.5.1.Device Release.By considering the physics-based diffusion model described by Giakoumi et al., 47 the device's cumulative amount of nickel release (M d ) per surface area (A surf ) can be mathematically described as with a rate of where α represents the quantity of nickel connected to the surface per unit area, and τ denotes the characteristic release time.2.5.2.Biokinetics.Giakoumi et al. 47 proposed a compartmental time-variant PBTK model for Ni leaching from cardiovascular stents by extending one proposed by Saylor et al. 45 along with a method for determining the confidence intervals associated with the model predictions.The model demonstrated prognostic conformity with reported data following the implantation of nickel-containing cardiovascular devices in humans 46 and minipigs 32 in the compartments of serum, urine, and local tissue while providing nonvalidated predictions for the compartments of kidney and other tissues.However, the unavailability of measurements from other tissues and organs did not enable further validation and development of the model by separating additional compartments that may be required for the thorough toxicological evaluation of the device.In this work, we present a generalized version of the model proposed by Giakoumi et al., 47 wherein the other tissue compartment is separated into the liver, brain, lungs, gut, and fecal compartments (Figure 1).The model considers a zero-order rate of absorption of dietary nickel from the diet into the gut (k diet ), an irreversible first-order elimination of Ni from the gut through feces (k f ), an irreversible first-order elimination of Ni from the kidney to urine (k u ), a reversible first-order exchange of Ni between blood and local tissues (k ltbl and k bllt ), blood and liver (k livbl and k blliv ), blood and brain (k brbl and k blbr ), blood and lungs (k lubl and k bllu ), blood and kidney (k kbl and k blk ), blood and gut (k gbl and k blg ), and the rate at which nickel is released from the medical device, (M ̇d) (eq 2).A fraction of the nickel (0 ≤ F (t) ≤ 1) is directly released into the blood, whereas the residual nickel is released into the local tissue encircling the device.The evolution equations of the PBTK model are as follows where M lt , M bl , M liv , M br , M lu , M k , M g , M u , and M f are the nickel mass in the local tissue, whole blood, liver, brain, lungs, kidney, gut, urine, and fecal compartments, respectively.C u and C f correspond to the concentration of nickel in urine and feces, respectively, whereas Q u and Q f correspond to the volumetric urine and fecal output rates, respectively.Note that most of the kinetic rates are considered as time-dependent.The initial conditions for all of these quantities at t = 0 are set equal to their control/baseline values This enables the calculation of the respective kinetic rates at t = 0 By implementing non-negative, time-dependent kinetic rates, k i (t), along with a fraction F(t), the model can be considered as a timevariant system of ordinary differential equations (ODEs), for which in general no analytical solution exists, thus providing further flexibility for simulating the dynamics of Ni biodistribution in all of the compartments.Giakoumi et al. 47 introduced a symmetric timedependent function characterized by a Gaussian "pulse" behavior, expressed in the form of where A (nickel mass/day) determines the pulse's maximum height, B (1/day) is the pulse's width, C is the abscissa of the saturation point, and D (nickel mass/day) is the kinetic rate's steady-state value.Finally, fraction F(t) takes the form of The kinetic rate parameters A−D are treated as fitting parameters that can be optimized using the available in vivo experimental data, as discussed in Section 2.6.However, it should be noted that certain parameters can be directly calculated if the control values are experimentally available by utilizing eqs 5a−5h.It is important to highlight that in this particular case, when t = 0, only parameters A, C, and D can be calculated directly instead of being obtained through the fitting process.In this study, it was decided that it was more convenient to calculate the steady-state parameter (D) rather than the transient parameters A and C during the optimization procedure.

Probabilistic Modeling.
According to the methodology described by Giakoumi et al., 47 the implementation of probabilistic modeling techniques is considered more appropriate for evaluating the PBTK model prediction uncertainties associated with such a complicated system as the human body.This can be achieved by considering the fundamental physiological hypothesis of any PBTK model, which dictates that each kinetic rate must comprise the same values regardless of the device used but depending only on the mammalian-specific (device host) absorption, distribution, and excretion properties.Consequently, with the parameterization of the kinetic rates' parameters for each device (see Section 2.6), it is possible to construct the probability density function (PDF) for each parameter.
The construction of the PDFs of the kinetic rates' parameters, e.g., f Akd i for parameter A of the kinetic rate k i , can be achieved with the implementation of point estimation techniques, such as the maximum likelihood estimator (MLE), 50 by considering each kinetic rate's parameters as a random variable, e.g., A, with observed values, e.g., A 1 , A 2 , ... , A n , in a random sample of size n (the number of tested devices).With the assumption that the kinetic rates' parameters follow a log−normal or normal distribution based on the domain of the particular kinetic rates' parameter, two unknown PDF's parameters remain to be determined: the mean μ and the variance σ 2 for a normal distribution (or the mean θ or variance ω 2 for a log− normal distribution, the two types of distributions used in this paper), which are considered as optimization parameters that should maximize a function.This function, called MLE, takes the form of the likelihood function of the sample.For example, assuming that the parameter A of the kinetic rate k i follows a normal distribution, the likelihood function is defined as The likelihood functions are maximized at the maximum likelihood estimators μ̂A kd i and σ̂A kd i 2 .
With the determination of the estimators for all of the kinetic rates' parameters and the construction of their PDF, the PBTK model can be solved using probabilistic methods, such as the Monte Carlo (MC) probabilistic method, 51 where each kinetic rate parameter's value is randomly drawn from its PDF through a series of random virtual experiments (virtual experiments that result in a negative concentration in any compartment at any time are interrupted and discarded).This will produce the biodistribution profile in various tissues/organs and body fluid compartments, which can be depicted in a Whisker plot featuring the median value, the interquartile range (IQR), and the upper−lower Whisker range (maximum = Q3 + 1.5 × IQR and minimum = Q1−1.5 × IQR) of the data. 50.6.Model Parameterization: Optimization.For the determination of the kinetic rates' parameters for each tested device, an optimization problem must be constructed.In this work, we choose to use as the objective function the root-mean-square error (RMSE) between the prediction of the PBTK model (M i (t)) and the experimental data j (j ∈ [1, N i ], where N i is the total number of experimental data in compartment i).Since the solution of the PBTK model provides the prediction for each compartment i, with a total of P compartments (P ≥ 2), and assuming that there are experimental data available for each compartment, then the objective function can be considered as a multiobjective function, defined as min(RMSE , RMSE , , RMSE ) , where k is the vector that includes the kinetic rates' parameters (A, B, C, D) and the F(t) parameters (A, B) to be optimized (as defined by Giakoumi et al. 47 ), and K is the feasible set of k.Note that the optimal value of k is the one that aligns with the solution's Pareto front, of the multiobjective optimization problem, as determined by a preference function.
One way is to define the preference function as implemented by Giakoumi et al. 47 Briefly, by encoding the preferences as a vector of weights 52 and by using the weighted sum method, 53 the multiobjective function is converted to a single objective function.Since currently there is no indication in the literature regarding the importance of any specific compartment relative to the rest of the compartments, the same trade-off weight factors 52,53 can be used, redefining the optimization problem as As such, we may now use the single objective optimization algorithm patternsearch from the Global Optimization Toolbox available in MATLAB. 54Note that eq 8 can be interpreted as the average error of the model output and the experimental data in units of mass and can thus be used to assess the "goodness" of the fit or to compare the various predictions along different devices.
The feasible set K 47 is determined by the domain and linear inequalities related to the kinetic rates' parameters as A (9b) whereas for the fraction F(t), the domain and linear inequalities are given as For the establishment of a well-defined objective function, it is important to acknowledge the possibility of having experimental data in the acquired sample that, due to the complexity of their acquisition, are sources of experimental errors.To account for such errors, we assign a weight factor to every experimental data point j in compartment i (w ij ), indicating the level of confidence in the data.A numerical value for these weight factors can be assigned by considering limiting conditions, such as (1) The Ni mass in any compartment cannot be less than its control value (M i1 ), described as M ij ≥ M i1 ,j > 1.
(2) The Ni mass in any compartment cannot exceed the total cumulated release of the device M d (eq 1) plus the control value, described as If these conditions are met, then a 100% confidence is assigned (w ij = 1); otherwise, a confidence less than 100% must be assigned, which in this study is quantified by using a Gaussian radial basis function (RBF). 55The RBF is a kernel-type function that quantifies the similarities between data.In this work, we choose where 2 ln( ) , with ( ) The parameter σ i c refers to the violation of the control value (first limiting condition), and the parameter σ i r refers to the violation of the total cumulated device release (second limiting condition).If these parameters are not preset, they can be estimated by considering the equalities in eqs 11a and 11b after determining w ij (e.g., from eq 11a w ij = M ij /M i1 , which defines the maximum percentage reduction on the experimental data j in the compartment i to stop violating the limiting conditions).However, note that there is only one parameter σ i for each compartment i, and therefore if more than one data violates the conditions, then the determined w ij will not have the same σ i parameter value despite referring to the same compartment.In this situation, the optimized σ i is derived through an optimization procedure similar to the one applied for the in vitro parameters (with the starting point of the optimization procedure to be the mean of the different resulting σ i parameters in compartment i) and then applied to recalculate the w ij of these data via eq 10.Therefore, the level of confidence for each data defined by eq 10 can be included in the objective function defined in eq 8 through the redefinition of RMSE i where the term ∑ j N j w ij normalizes the data among compartments.
Finally, according to Giakoumi et al., 47 the in vitro parameters of eq 2 can be obtained by a standard nonlinear optimization algorithm.
2.7.Tolerable Intake Analysis.We have recently 47 proposed a method for performing a tolerable intake (TI) analysis based on FDA recommendations 27 of a TI value for parenteral (nonoral) exposure to nickel of 0.5 μg/kg/day. 27,30As per the method, the device's nickel release (eq 2) and release rate profile predicted in each compartment (dM i (t)/dt) are evaluated against the TI value for nickel to assess if the design and physicochemical characteristics of a device are toxicologically safe or unsafe.We hypothesized that the TI value per compartment can be evaluated via TI i = ⟨y i ⟩ × 0.50 μg/(kg.day),where ⟨y i ⟩ = ⟨M i ⟩/⟨M total ⟩ represents the time-averaged mass fraction of nickel in the i th compartment, which can be estimated stochastically using available experimental data and the mass balance equation.For the newly proposed model (Figure 1), the following mass balance equation can be defined by considering the control/baseline mass values, the Ni diet intake, and the total devise release

Stent Surface Characterization and Ni Ion
Quantification.Stent surface processing resulted in visually different colors for the two types (HT vs EP), indicating unique oxide layers.HT stents were dark in appearance, whereas EP stents had a bright, shiny metallic appearance.According to the manufacturer, the resulting oxide layer thickness could range between 5 and 100 nm for EP-and HT- treated surfaces, respectively.Yet, the samples' surface topography did not allow us for any precise oxide layer thickness measurements, neither via X-ray reflectometry (XRR) nor ultraviolet−visible (UV−vis) spectroscopy (data not shown).
Nontested (control) HT stents had inherent manufacturing surface features and imperfections that pre-existed the in vitro and/or the in vivo testing (see Figure 2d), while EP stents had smooth consistent surfaces throughout.The evaluation also confirmed that the surface roughness decreases after the various processing steps, from HT to EP (Figure 2d−f).Next, the 60 day immersion testing revealed no evidence of significant uniform corrosion for any group when compared to preimmersion samples (see Figure 2g−i).On the contrary, explanted stents from mice, HT, in particular, displayed rougher surfaces predominately on the abluminal (outer) part compared to their nonimplanted controls (see Figure 2j−l).Some microcracks were observed with "bubbling" of the oxide indicating potential subsurface corrosion (see Figure 2j).In any case, the EDS analysis did not expose any differences in the Ni/Ti ratio between tested (immersed or explanted) and nontested (control) stents since the levels of Ni ion leach were typically below the 1−2% by weight EDS detection limit (data not shown).
The concentration of nickel ions in all test solutions was quantified through ICP-MS.The resulting detection limit for Ni was below 1 ppb, and the reproducibility observed after triplicate measurement per sample resulted in low standard deviation values.Nickel levels in control blanks (PBS in tubes with no stent specimen) averaged 5.7 ppb and were deducted from all measurements from immersion testing.The cumulative Ni release profile for the V1HT, V1EP, and V2EP devices, following the 60-day immersion testing, is shown in Figure S1.The detected concentration in each step was multiplied by the solution volume per specimen (i.e., 3.33 mL) to calculate the total mass of nickel (in ng) leached from each device.Regardless of the geometric design (V1 or V2), HT stents exhibited the greatest amount of Ni release throughout 60 days of immersion (Figure S1), confirming that electropolishing is a much more effective surface processing method for controlling ion leaching from metallic biomaterials.Interestingly, the V1 design demonstrated higher cumulative release than the V2 design when comparing the EPtreated stents.This is most probably because the double-cell diamond-shaped pattern (V1) results in a larger number of stress raisers (connective links and struts) that promote stressinduced corrosion.Nickel levels were measured in harvested tissue/organs (stented arteries, kidneys, liver, lungs, brain, and small intestine), blood, urine, and fecal samples and were processed to determine the time-dependent concentration profile of implant-leached ions in the mammalian living system (see Figures 3−5).In general, HT stents demonstrated a significantly higher Ni ion release, most notably in the periimplant tissue, compared to EP stents with varying biodistribution patterns over the 18-week implantation period.These data were subsequently used to extract the devicespecific release parameters and the mammalian-specific biokinetic parameters during the model training phase and also to validate the model-derived predictions during the prognostic phase.

Application of the Model.
The methodology described in Sections 2.5.3 and 2.6 was employed to parameterize the developed PBTK model, using the in vitro immersion data (indicated in Figure S1) and the in vivo biodistribution measurements (indicated in Figures S2 and  S3), thus enabling the generation of probabilistic timedependent profiles for nickel levels across the various model compartments.The predictive value of the PBTK model was tested and validated against experimental data (independent training and testing sets) in an attempt to provide a straightforward link between device release characteristics and ion concentration−time profiles in tissues/organs, body fluids, and excreta.
The fitting process and model parameterization using the resulting in vitro and in vivo data can be found in the Supporting Information.Measurements from the V1HT and V1EP devices were used for model training, while data from the evaluation of the V2EP stent type were used as the evaluation (testing) subset.The physiological parameters considered for a 45 g CD1 mouse (Section 2.3.1) were as follows: V b = 3.6 mL, 56 Q u = 1 mL/day, 57 V u = 0.15 mL, 58 Q f = 0.094 g/day (in-house estimation), and k diet = 1250 ng/d (inhouse estimation).The fitted results for devices V1HT and V1EP are depicted in Figures S2 and S3, respectively.We have omitted graphs for the gut and fecal compartments due to high dispersion in the experimental data but mainly because the detected mass did not differ significantly from the control values (M g control = 2180 ng and C f control = 11 593 ng/mL) and was attributed merely to dietary intake and not the device release.
Figures 3 and 4 present the model-derived estimates after conducting 5000 Monte Carlo experiments for devices V1HT and V1EP, respectively.Each parameter was derived from its PDF, as described in Section 2.5.3.It is important to note that these outcomes are not predictive in the dimension of α and τ, as the sampling was performed using parameters trained from the same data set.Nonetheless, the results demonstrate a highly satisfactory description provided by the probabilistic model, as the majority of the data points fall within the uncertainty range, most notably for the V1EP device.This can be observed by the RMSE of each device and for each compartment, as presented in Table 1.The median is represented by red continuous lines, the interquartile range (IQR) is depicted with a dark red shading, and the upper− lower whisker range is represented by a lighter red shading.Furthermore, it is worth mentioning that the data points falling within the numerical uncertainty range are consistent with the resulting weight factors calculated based on Section 2.6, as presented in Table S2 of the Supporting Information.Figure 5 illustrates the predictions obtained by the model, along with a comparison against the in vivo data from V2EP, which were not included in the data set used for training the model.While the prediction does not exhibit the same level of satisfaction as for the previous devices, it is noteworthy that the majority of the data points that fall outside the uncertainty range were detected below the control/baseline levels.Consequently, these data points do not indicate a violation of toxicological safety, suggesting that deviations from the expected range are acceptable.

Toxicological Risk Assessment.
The PBTK modelderived outputs were employed to assess whether nickel exposure levels exceed permissible limits and determine the toxicological safety of implants with different surface finishes.
In Figure S4 (see the Supporting Information), the timedependent device release rate (M ̇d), was compared to the total TI value, revealing that all devices surpass the predefined exposure threshold for approximately the first 20−70 h after implantation.
For evaluating the toxicological exposure in individual tissue/organs and body fluid compartments, the release rate profiles (dM i /dt) for all compartments in the model were compared against the compartmental TI i values, as defined in Section 2.7.Figures 6−8 present the results of this analysis for the V1HT, V1EP, and V2EP devices, respectively, which exhibit similar profiles, with the only discernible difference being an increased response and wider numerical uncertainty before day 1 for the V1HT and V1EP devices.When focusing the evaluation on each separate compartment, local tissue is the only compartment that in all cases does not exceed the permissible limits.On the other hand, the nickel release rate in the blood exceeds the permissible limits during the first few hours post implantation and over again after day 10, with higher values denoted for the V1HT device.The liver, brain, and kidney compartments exhibit similar behavior as in blood but with the early-stage threshold violation expanding until day 1.Conversely, the lung compartment reveals a toxicological risk between days 5 and 30, whereas urine increases a concern only within the first 24 h, with the level and timespan of the predicted violations decreasing between the V1HT, V1EP, and V2EP devices.

DISCUSSION
−63 The use of computational modeling and simulation (CM&S) in regulatory submissions is rapidly increasing and together with bench, nonclinical in vivo, and clinical studies can also be used to evaluate the toxicological safety and effectiveness of medical devices. 64,65The FDA 66 and the EU Reference Laboratory for Alternatives to Animal Testing (EURL ECVAM) 67 promote and facilitate the use of nonanimal methods in testing and research that have the potential to provide timelier and more predictive information to assess certain aspects of regulated products.
For many years, CM&S studies have been employed to support engineering design and structural analysis 68−70 as they have the potential to streamline the development process and alleviate challenges related to premarket device evaluation.−73 To ensure the accuracy and precision of results for supporting regulatory submissions, it is essential to develop processes and approaches that foster consistency in the execution and review of CM&S. 33,74While prior art provides the means of modeling the diffusion of substances in mammalian compartments and comprehensively considers the release from implantable devices, 45,75−80 it does not provide a practical and cost-effective way of guiding the optimal design of implantable devices to ensure compliance with a permissible exposure limit of different mammalian compartments to concentrations of released substances including intervals of confidence.There is therefore the need for a robust and validated modeling technology that overcomes the abovementioned limitations so that device manufacturers can optimize device design and ensure that the release of substances will remain below acceptable thresholds.
Herein, a multicompartment PBTK model is presented along with a probabilistic Monte Carlo methodology for toxicological risk assessment of vascular stents.In this study, we have modified the PBTK model proposed by Giakoumi et al., 47 considering the major organs (brain, lungs, and liver) associated with a toxicological concern as well as blood and excreta (urine and feces) as separate compartments.To the best of our knowledge, this is the very first study that reports on the in vivo concentration of nickel ions in all of these compartments.The PBTK model was trained (parameterized) by using only a subset of the available in vivo data (those obtained using the devices V1HT and V1EP; Figures 3 and 4) and inputs using in vitro nickel release measurements; the trained model was then employed to make a prediction for a third device (V2EP, Figure 5) following the leave-one-out cross-validation (LOOCV) approach. 81lthough the latter prediction (V2EP) did not fully meet the specified performance requirements after the training model, it was found to be in good agreement with the majority of the data points that lie beyond the uncertainty range.These are the ones below the control level, which are expected to be erroneous and can also explain why the RMSE of V2EP is very similar to that of V1HT.The physical limitations of the data below the control level play a detrimental role here, resulting in smaller possible disparities between the prediction profile and the outlier data.In contrast, the difference between the prediction profile and the outlier data, which exceeds the total cumulative release, appears to have no limitations and is therefore expected to have significantly higher disparities.Consequently, employing RMSE as a quantitative measure for the model's predictive performance may not be precise in this scenario, underscoring the necessity for future studies to investigate more representative methods to evaluate the prediction performance numerically.
The optimization procedure presented in this work represents a significant improvement over the previous methodology, 47 primarily due to the availability of experimental data in all compartments, allowing for a more comprehensive development of the methodology.However, this advancement has also increased the complexity of the optimization procedure, especially by incorporating a larger sample of experimental data from different compartments that are also subjected to experimental errors and uncertainties.In such cases, having a robust and well-defined objective function becomes crucial for the effective guidance of the implemented algorithm toward the desired basin of attraction points while avoiding stagnation at irrelevant points that do not represent the physiology of the organism under consideration.To simplify the optimization process and address the impact of experimental uncertainties, we take two key steps.First, we transform the model from a multiobjective function to a single objective function 52 by considering an equalization of the importance of each compartment with weight factors (eq 8) 53 and therefore providing a much more simplifying objective function for the algorithm to handle.Second, the methodology accounts for experimental errors in the sample data by integrating them as an intrinsic part of the objective function using kernel-type function 55 (eq 10), which imparts smoothness to the objective function.The fact that the kinetic rates are considered time-dependent also provides enhanced flexibility to handle in vivo data from various mammals even in cases of contradictory behavior between the experimental measurements.
In Section 3.2, Figures 3 and 4 demonstrate that the integration of experiment errors as an intrinsic part of the objective function, based on eq 10 and Table S2 in the Supporting Information, was a successful implementation, with the majority of data points falling within the numerical uncertainty range consistently aligned with the resulting weight factors.However, in certain compartments, like the brain, significant dispersion in data among the devices was observed, not entirely consistent with the high-value weight factors.This discrepancy arises from high dissimilarities between data from devices V1HT and V1EP in the brain compartment, affecting the fitted parameters k brbl and k blbr .Nonetheless, the model's accurate fit to the brain compartment, as shown in Figures S2d  and S3d of the Supporting Information, supports the overall success of the implementation of the weight factors in the optimization procedure.
The methodology presented here can be further applied to increase the level of confidence in experimental data from in vitro testing setups using a probabilistic stochastic (Monte Carlo) methodology.Nagaraja and Pelton 25 studied the nickel elution resistance of electropolished nitinol ocular microstents via immersion testing and showed that the cumulative nickel release profiles of different sets of the same tested samples varied over a 63-day immersion duration.Predictive stochastic simulations could be used as an additional tool in regulatory testing to estimate the median values and the dispersion from multiple experiments and therefore upgrade the credibility of standard guides for the corrosion and elution resistance of medical implants.It is important to emphasize that the proposed methodology is generic and applicable to any simulation engine that requires training with hard-to-obtain data characterized by high uncertainty and experimental errors, providing a valuable approach for modeling complex systems under such conditions.We note that it should be possible to extend the framework presented in this paper to other device applications and metal ions.For instance, cobalt−chromium alloys, widely utilized in both vascular and orthopedic device applications, could be considered, provided there is adequate data for parameterization and, if needed, adjustment of the suggested multiscale model.
As made evident from the literature survey, this is the first fully comprehensive study that reports on a multicompartment PBTK model that includes the most important tissues/organs and the necessary in vivo data to parameterize it.Other studies, for simplicity, have employed PBTK models that include only very few of the organs in separate compartments, and all of the rest are included in an "other tissues" compartment.This then poses the question of whether our multicompartment PBTK model can be "simplified" by lumping together some tissue compartments.It can be easily shown mathematically that this can indeed be done (see the Supporting Information) since the time-averaged mass fractions of the tissue compartments that are lumped are known.However, the reverse (splitting) cannot be done even if the time-averaged mass fractions are assumed.This strongly exemplifies the necessity to acquire measurements from as many organs and tissues as possible.
Herein, the model outputs aligned with existing in vivo data, indicating the potential use of in vitro nickel release testing for estimating both local and systemic exposure.Yet, a constraint of the current model lies in assuming a static local implant environment.The local biomechanical environment plays a decisive role in the degradative performance of an implant. 28he complex and varying in vivo settings comprising dynamic geometries, cyclic and dynamic loading profiles, and harsh chemical conditions (chloride, dissolved oxygen, and pH levels) can affect the corrosion susceptibility that results from surface damage or a change in ion diffusion kinetics.Hence, it is crucial to comprehend and assess the diverse interactions that an implantable device is anticipated to undergo in vivo in order to identify the requisite testing for establishing a reasonable assurance of safety.
In addition, vascular devices undergo dynamic physiological changes, owing to inflammation and gradual tissue coverage.This leads to alterations not only in local biochemistry but also in the evolving geometry at the implant site over time.As the device becomes enveloped by the endothelium, the initially assumed model geometry, which considered nitinol wires to be only partially embedded in peri-implant tissue, is no longer a valid approximation.However, incorporating the rate of coverage, as now captured by F(t), facilitates a more accurate consideration of the impact of this process on local nickel accumulation within the model.Moreover, in certain cases, peri-implant tissue nickel concentrations are more sensitive to release and may surpass local thresholds for toxic effects, even at relatively low release rates compared to the chronic threshold for typical device sizes.In other words, the variability in blood nickel levels within the patient population may obscure the detection of nickel release, which could cause adverse effects locally.This also implies that local effects due to nickel release are more likely than systemic effects in most vascular device applications.Nonetheless, it should be emphasized that the results presented in this study were based on a specific device geometry assumption.
Future work will emphasize advancing the technology using inputs from a larger animal model, like swine, which is a wellestablished model, that will allow us to resolve size-related issues that previously hindered sample collection and detection limits.Data from a wide range of potential manufacturing processes, in terms of bulk material, geometric design, total active surface area, and surface processing methods, are required for extensive model parameterization.A wide data set will also facilitate the implementation of perturbation theory and machine learning (PTML) methods early in the device design cycle to optimize device-specific characteristics to minimize ion leaching.Large animal models offer important translational features and the ability to apply human-like settings, which increases the chances of bench findings being translated into effective clinical tools.Species-specific physiology between the different laboratory animals and further on to human trials is, of course, challenging because of the uncertain comparability of physiological processes.Different knowledgedriven approaches for considering physiological and biochemical differences between them should be evaluated by systematically incorporating specific model parameter domains of a target species into the PBTK model of a reference species.Cross-species extrapolation and ultimately full model validation with coherent clinical data will enable such novel technology to exploit the increased market acceptance and benefit fully from the entry into this rapidly expanding sector.

CONCLUSIONS
A novel multicompartment and time-variant physiologically based toxicokinetic (PBTK) model was developed and parameterized with in vitro immersion tests and in vivo implantation studies.The model predictions were complemented with a probabilistic Monte Carlo approach and a simulation engine to estimate the concentration−time profiles, along with confidence intervals, of implant-leached ions in all of the major organs associated with a toxicological concern as well as in peri-implant tissue, body fluids, and excreta.However, more robust in vivo implantation studies from large animals are required to further improve the feasibility and reliability of the model.Despite its limitations, the proposed model may provide the basis for the broader use of prognostic tools to guide the optimal design of implantable devices in compliance with exposure limits and other regulatory requirements.
■ ASSOCIATED CONTENT financial interests or personal relationships that could have appeared to influence the work reported in this paper.

■ ACKNOWLEDGMENTS
This work was cofunded by the European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation (Project: POST-DOC/0916/ 0237).

Figure 1 .
Figure 1.Diagram depicting the absorption, diffusion, distribution, and excretion principles of the proposed PBTK model.

Figure 2 .
Figure 2. Representative scanning electron microscopy (SEM) images of control (prior testing) vs tested (in vitro and in vivo) stents.(a−c) The two stent designs, V1 and V2, after HT and EP surface processing (at ×60 mag).(d−f) Higher magnification (×2000) images showing inherent manufacturing features and imperfections in (d) HT stents.(g−i) Stent surface appearance and morphology post the 60 day immersion testing (at ×2000 mag).(j−l) Examples of different corrosion morphologies observed in explanted stents (at ×2000 mag).

Figure 3 .
Figure 3. Model-derived output of the nickel concentration−time profiles, released from V1HT stents, in the various tissue/organ and body fluid compartments.

Figure 4 .
Figure 4. Model-derived output of the nickel concentration−time profiles, released from V1EP stents, in the various tissue/organ and body fluid compartments.

Figure 5 .
Figure 5. Model-derived prediction of the nickel concentration−time profiles, released from V2EP stents, in the various tissue/organ and body fluid compartments.

Figure 6 .
Figure 6.Model-derived compartmental release rates (dM i /dt) for V1HT.To enable toxicological risk assessment, the horizontal dashed line corresponds to the compartmental (depicted based only on the median y i value obtained from 5000 Monte Carlo computational experiments) TI i per body mass, as defined in Section 2.7.

Figure 7 .
Figure 7. Model-derived compartmental release rates (dM i /dt) for V1EP.To enable toxicological risk assessment, the horizontal dashed line corresponds to the compartmental (depicted based only on the median y i value obtained from 5000 Monte Carlo computational experiments) TI i per body mass, as defined in Section 2.7.

Figure 8 .
Figure 8. Model-derived compartmental release rates (dM i /dt) for V2EP.To enable toxicological risk assessment, the horizontal dashed line corresponds to the compartmental (depicted based only on the median y i value obtained from 5000 Monte Carlo computational experiments) TI i per body mass, as defined in Section 2.7.
EU of the European Parliament, National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals).Male CD1 mice, weighing 45 ± 5 g (10−12 weeks old), chosen for their relatively large body size, were used throughout the study.

Table 1 .
Root-Mean-Square Error of Each Compartment (RMSE i ) for the Median Output of the Probabilistic PBTK Model for the Devices V1HT, V1EP, and V2EP