A physiologically based pharmacokinetic model for polyethylene glycol-coated gold nanoparticles of different sizes in adult mice.

Nanoparticles (NPs) are widely used in various fields of nanomedicine. A systematic understanding of NP pharmacokinetics is crucial in their design, applications, and risk assessment. In order to integrate available experimental information and to gain insights into NP pharmacokinetics, a membrane-limited physiologically based pharmacokinetic (PBPK) model for polyethylene glycol-coated gold (Au) NPs (PEG-coated AuNPs) was developed in mice. The model described endocytosis of the NPs in the liver, spleen, kidneys, and lungs and was calibrated using data from mice that were intravenously injected with 0.85 mg/kg 13 nm and 100 nm PEG-coated AuNPs. The model adequately predicted multiple external datasets for PEG-coated AuNPs of similar sizes (13-20 nm; 80-100 nm), indicating reliable predictive capability in suitable size ranges. Simulation results suggest that endocytosis of NPs is time and size dependent, i.e. endocytosis of larger NPs occurs immediately and predominately from the blood, whereas smaller NPs can diffuse through the capillary wall and their endocytosis appears mainly from the tissue with a 10-h delay, which may be the primary mechanism responsible for the reported size-dependent pharmacokinetics of NPs. Several physiological parameters (e.g. liver weight fraction of body weight) were identified to have a high influence on selected key dose metrics, indicating the need for additional interspecies comparison and scaling studies and to conduct pharmacokinetic studies of NPs in species that are more closely related to humans in these parameters. This PBPK model provides useful insights into the size, time, and species dependence of NP pharmacokinetics.


Introduction
Gold (Au) nanoparticles (NPs) are widely used in various fields of nanomedicine, including as medicinal or diagnostic agents for tumors or rheumatoid arthritis, and potentially as carriers for the delivery of drugs, oligonucleotides, and peptides (Arvizo et al., 2012). Their widespread use results in intentional high-dose exposure for diagnostic or therapeutic purposes, and may increase the likelihood of unintentional exposure during the manufacturing process and in the environment (Khlebtsov & Dykman, 2011).
Animal studies show that AuNP exposure results in acute inflammation in the liver (Cho et al., 2009a), changes in the expression of genes involved in regulating the cell cycle, apoptosis, inflammation, oxidative stress, and metabolism in the liver and spleen (Balasubramanian et al., 2010;Cho et al., 2009b;Geffroy et al., 2012). The degree of toxicity of AuNPs is associated with their pharmacokinetics and target organ dosimetry, which depends on surface coating, dose, size, and species (Landsiedel et al., 2012;Riviere, 2009). In order to properly assess the potential toxicity of AuNPs, it is crucial to have a comprehensive quantitative tool that can predict the pharmacokinetics of various types of AuNPs in different species.
One approach that can be used to integrate available pharmacokinetic data of AuNPs and, hence gain deeper insights into the target tissue dosimetry and species differences is physiologically based pharmacokinetic (PBPK) modeling (Li et al., 2010). PBPK models have been commonly developed for drugs to help in the design of therapeutic regimens (Jones et al., 2013;, to aid predictive risk assessment for environmental toxicants (Andersen et al., 1987;Fisher et al., 2011;Lin et al., 2013), as well as to predict drug withdrawal time in food-producing animals (Buur et al., 2008;Leavens et al., 2014). In the field of nanomedicine, researchers have begun to employ PBPK models to simulate the pharmacokinetics of NPs, including quantum dots (Lee et al., 2009;Lin et al., 2008), carbon (Pery et al., 2009), silver (Bachler et al., 2013;Lankveld et al., 2010), titanium dioxide (Bachler et al., 2014), polyacrylamide (Li et al., 2014), poly(lactic-co-glycolic) acid (Li et al., 2012) NPs, and gold/dendrimer composite nanodevices (Mager et al., 2012). These models greatly improve our understanding of the pharmacokinetics of NPs; some of these models have been extrapolated to humans and proved to be helpful in risk assessment (Bachler et al., 2013(Bachler et al., , 2014. PBPK models have not yet been reported for AuNPs. Among the available nine published models, only three of them attempt to simulate the endocytosis of NPs (Bachler et al., 2013(Bachler et al., , 2014Li et al., 2014), a key process affecting their biodistribution. Size dependence of endocytosis has not been assessed. Additionally, although validation with external datasets is an important aspect for developing reliable PBPK models, only three of the available models have been validated thus far (Bachler et al., 2013(Bachler et al., , 2014Li et al., 2012).
In order to improve our understanding of the pharmacokinetics of AuNPs and considering the above-mentioned data gaps, the objective of this study was to develop a PBPK model for PEG-coated AuNPs. PEG-coated AuNPs were selected because PEG coating is one of the most commonly used surface modifications for NPs. Systemic organ toxicity of PEG-coated AuNPs with varied sizes has been reported (Cho et al., 2009a,b). Pharmacokinetic data of PEG-coated AuNPs of different sizes (4-100 nm) in rodents are available (Cho et al., 2009a(Cho et al., , 2010Liu et al., 2013;). These multiple datasets made it possible to assess the role of particle size in cellular uptake and biodistribution of AuNPs, and to evaluate the model with independent data not used in the model calibration.

Data source for model calibration
The PBPK model was calibrated with the experimental data from Cho et al. (2010). Briefly, 6-week old male BALB/c mice were injected in the tail vein with 100 mL (0.85 mg/kg) of PEG-coated AuNPs (4 nm, 13 nm, or 100 nm). At 0.5, 4, 24, 168 h, 1, 3, and 6 months after injection, animals (n ¼ 9 each time point in each group) were sacrificed and their plasma and tissue (e.g. liver, spleen, kidneys, and lungs) samples were collected and analyzed for Au concentrations with inductively coupled plasma-mass spectroscopy (ICP-MS). Data were extracted using the WebPlotDigitizer (version 2.6, Austin, TX, http://arohatgi.info/WebPlotDigitizer/). The 13 nm and 100 nm datasets were used to calibrate the models for smaller and larger AuNPs, respectively. The 4 nm dataset was not modeled because there are no independent data of similar sizes for validation. The present model focused on data up to 7 d after injection because there are sufficient short-term kinetic data to develop a reliable PBPK model and the information about the long-term kinetics of AuNPs is very limited. In order to inform long-term kinetic and PBPK modeling studies, preliminary long-term kinetic simulations up to 6 months post-exposure in Cho et al. (2010) were performed and are discussed in the Supporting material.

Model structure
The model structure was based on published NP PBPK models (Bachler et al., 2014;Li et al., 2014) and on the extractable data from Cho et al. (2010). It consisted of seven compartments: plasma (or blood), liver, spleen, kidneys, lungs, brain, and rest of body ( Figure 1). As shown in other NP PBPK models, the biodistribution of AuNPs was assumed to be governed by two processes: (1) the ability to cross the capillary membrane of the organs and (2) to be taken up via endocytosis (Bachler et al., 2013(Bachler et al., , 2014. The transportation of AuNPs from the blood to the tissue compartments was described with both perfusion-limited (also termed flow-limited) and membrane-limited (also termed diffusion-limited) models. Simulation results using different modeling approaches were compared to decide the final model structure.
The parameter ''distribution coefficient'' was used to account for the uneven distribution of AuNPs between interstitial fluid and plasma in each organ. Inclusion of this parameter was necessary because the composition of AuNP-protein biocorona differs due to different protein compositions between these two locations (Li et al., 2014;Nel et al., 2009). In support of our approach, the majority of available NP PBPK models also included this parameter (Li et al., 2012(Li et al., , 2014Lin et al., 2008). However, some of these models termed this parameter as the ''partition coefficient'' (Li et al., 2012(Li et al., , 2014, which is a parameter commonly used in PBPK models for small molecules. It should be noted that the ''partition coefficient'' represents the ratio of concentrations of a dissolved compound between two phases at thermodynamic equilibrium. Unlike small molecules, NPs do not form solutions, but colloidal dispersions, which are thermodynamically unstable, and thus are not applicable to ''partition coefficient'' (Praetorius et al., 2014). Therefore, the authors used the term ''distribution coefficient'' in the present model. Further discussion on the role of ''distribution coefficient'' in the NP PBPK modeling can be found in the Supporting material.
The authors used the term ''endocytosis'' as an operational term to represent different cellular uptake mechanisms of NPs, including phagocytosis (actin-dependent, ''professional'' endocytosis), macropinocytosis (receptor-independent), and receptormediated (clathrin-, caveolin-, and raft-dependent) endocytosis (Dykman & Khlebtsov, 2014). Another term ''phagocytosis'' was used to present these cellular uptake processes in several recently published NP PBPK models (Bachler et al., 2013(Bachler et al., , 2014Li et al., 2014), but the authors contend that ''endocytosis'' is a more representative term as it also includes professional phagocytosis. In addition, we used ''phagocytic cells (PCs)'' as an operational term for a wide variety of cells that can engulf NPs, including those that other investigators have called reticuloendothelial system cells (Cho et al., 2010), mononuclear phagocyte system cells (Bachler et al., 2014), or organ-specific cells, such as liver Kupffer cells, splenic macrophages, lung macrophages, and kidney mesangial cells (van Furth et al., 1972). Hence, a subcompartment of PCs was included in organs where these cells are primarily located (i.e. liver, spleen, kidneys, and lungs). Available biological evidence suggests that NPs can be taken up via multiple endocytic pathways (Canton & Battaglia, 2012;Dykman & Khlebtsov, 2014). From the computational perspective, the authors classified these pathways into two general pathways based on the cellular location of cellular uptake: (1) directly from the blood to endothelial cells or macrophages located at the tissue capillary blood vessels, and (2) if the size is smaller than the upper limits of pore size of capillary walls, diffused into tissue and further taken up by tissue macrophages (Dykman & Khlebtsov, 2014). However, it will substantially increase the model's complexity and uncertainty if we simulate two uptake pathways because there are no sufficient data granular enough to separate out mechanisms and to validate this approach. Therefore, based on the approaches used in existing NP PBPK models, we tried to describe endocytosis of 13 nm and 100 nm AuNPs either from the blood (Bachler et al., 2014) or from the tissue (Li et al., 2014). We found that for 13 nm NPs, it is better to describe endocytosis from the tissue pathway ( Figure 1A), while it is more appropriate to simulate endocytosis of 100 nm NPs from the blood pathway ( Figure 1B). This model can be easily modified to integrate both uptake pathways once additional data become available.
Metabolism of PEG-coated AuNPs was considered negligible for the time-scale of this study because of their reported high stability both in vivo and in vitro (Manson et al., 2011;). Dissolution was not considered because there is no evidence of AuNP dissolution in the body . Excretion was through both biliary and urinary pathways based on the experimental data (Cho et al., 2010).

Main mathematical description of the model
Endocytosis of NPs was described with both a linear equation and the Hill function. Simulation results using different approaches were compared. More accurate predictions were obtained with the Hill function and it was used in the final model as shown in the following equation: where T (h) is the simulation time, K up t (per h) is the uptake rate parameter of NPs by PCs in the organ t at a particular time T, K max t (per h) is the maximum uptake rate parameter in the organ t, K 50 t (h) is the time reaching half of K max t , and n t (unitless) is the Hill coefficient. The rate of change in the amount of AuNPs in the PCs equals the uptake rate from the tissue for 13 nm NPs (or from the blood for 100 nm NPs) minus the release rate from PCs back to the tissue (or blood). The actual uptake rate (mg/h) depends on the amount that is available for uptake and the rate parameter as a function of exposure time. For example, the equation describing these processes for 13 nm NPs is shown in the following equation: where R up t (mg/h) is the uptake rate of the NPs from the tissue to PCs in the organ t, K up t (per h) is the uptake rate parameter of the NPs, A tissue t (mg) is the amount of the NPs in the tissue sub-compartment in the organ t, R release t (mg/h) is the release rate of the NPs from PCs to the tissue, K release t (per h) is the release rate constant of the NPs, A pc t (mg) is the amount of the NPs in the PCs in the organ t, and R pc t (mg/h) is the rate of change in the mass of the NPs in the PCs. In a membrane-limited model, each compartment contains both a vascular space and a tissue space, each of which should be described separately. As mentioned before, the mass transfer between the two sub-compartments was modeled on the basis of two primary kinetic processes, i.e. the ability of AuNPs to diffuse across the capillary wall of the organs and the extent of endocytosis from the tissue (or blood). For instance, the equations simulating these processes for 13 nm AuNPs are shown in the following equation: where R blood t (mg/h) is the rate of change in the amount of the NPs in the capillary blood sub-compartment of the organ t, Q t (L/h) is the blood flow to the organ t, C a (mg/L) is the concentration of the NPs in the arterial blood (or plasma), CV t (mg/L) is the concentration of the NPs in the venous blood of the organ t, PA t (L/h) is the permeability area cross product between the capillary blood and the tissue of the organ t (PA t is approximated as the product of permeability coefficient between capillary blood and tissue [PAC t : unitless] and regional blood flow [Q t : L/h]), C tissue t (mg/L) is the concentration of the NPs in the tissue sub-compartment, P t (unitless) is the tissue:plasma distribution coefficient for the organ t, and R tissue t (mg/h) is the rate of change in the mass of the NPs in the tissue subcompartment.
Equations describing the kinetics of 100 nm AuNPs in the organ and other key equations are provided and explained in the Supporting material. The PBPK model was coded using the acslX version 3.0.2.1 (Aegis Technologies Group, Inc., Huntsville, AL). The model code in CSL file is provided in the Supporting material.

Model parameterization
All physiological parameter values were from the literature and are given in Table S1 (Supporting material). NP-specific parameters were estimated by using both the Nelder-Mead maximum log likelihood estimation method in acslX (AEgis Technologies Group, Inc, Huntsville, AL) and a manual approach to obtain a visually reasonable fit to the experimental data. The following parameters were optimized: elimination rate constants for the liver and kidney, distribution coefficients and permeability coefficients, and endocytosis-related parameters. Parameter values for different types of NPs are usually different in PBPK models, but in order to create the most parsimonious model it is recommended to use a minimum number of estimated parameters or to apply a generic value for different parameters (Li et al., 2014). Therefore, in the present model, endocytosis-related parameters for each size (13 nm or 100 nm) were optimized separately against respective datasets (Cho et al., 2010) because these parameters were determined to be the most influential in the kinetics of NPs, while other less influential parameters, such as the distribution and permeability coefficients, were kept the same between the 13 nm and 100 nm models, unless it was required to explain the kinetics in a particular compartment. Values of all NP-specific parameters are listed in Table 1. More details regarding the step-by-step parameterization process refer to the Supporting material.

Model evaluation with independent data
The predictive capability of the 13 nm and 100 nm models was evaluated with four and one external datasets, respectively, for PEG-coated AuNPs of similar sizes that were not used in the model parameterization (Cho et al., 2009a;Liu et al., 2013; (Table S2, Supporting Material). The criteria of a validated model were based on World Health Organization (WHO) guidelines (WHO, 2010), namely if the simulations were generally within a factor of two of the experimental values, the model was considered reasonable. The overall goodness-of-fit between predicted and measured values was further analyzed with linear regression.

Sensitivity and uncertainty analyses
The impact of parameter uncertainty on model outputs was evaluated semiquantitatively based on the approach used in Teeguarden et al. (2005). First, we conducted normalized sensitivity analyses to identify highly influential parameters by examining the effects of 1% change of each parameter (p) on key dose metrics, including area-under-the-concentration curve (AUC) in the plasma, liver, and spleen, using the following equation: The normalized sensitivity coefficients (NSCs) were calculated for both 13 nm and 100 nm AuNPs to compare the difference between sizes, and for both 24 h and 168 h to capture the earlier and later kinetic phases, respectively. We also performed local sensitivity analyses by multiplying each highly influential parameter by a factor of 2 or 0.5 and compared the predicted concentrations of AuNPs in a given compartment (e.g. liver) over time (Li et al., 2014). Next, the uncertainty of model parameters with high impacts (i.e. |NSC| ! 0.5) on at least one of the selected dose metrics was evaluated qualitatively based on the following criteria (Teeguarden et al., 2005): Low (L): data were directly available for the parameter in the correct species or verified through successful use in published PBPK models. Medium (M): scaled value from a different species with a high possibility that scaling holds across the species. High (H): data were not directly available from the parameter. Appropriate assumptions had to be made. Values were from Li et al. (2014) except for the brain and the rest of body. The latter was estimated, together with K bile and K urine , by visually fitting to the measured concentration in the plasma from Cho et al. (2010). The value for the brain was set the same as the rest of body. Another exception is the permeability coefficient in the spleen for the 13 nm gold nanoparticles, which had to be increased to fit the measured high concentrations in the spleen from Cho et al. (2010). c These parameters were estimated by visually fitting to measured concentrations in respective organs from Cho et al. (2010). For a detailed description of the step-by-step parameterization, refer to Section 2 in the Supporting material.

Model development and calibration
Both perfusion-limited and membrane-limited models were constructed to determine which was the better model to simulate the pharmacokinetics of different size PEG-coated AuNPs. It was found that both models adequately simulated the 100 nm dataset, but the membrane-limited model performed much better than the perfusion-limited model for the 13 nm AuNPs (Table S3, Supporting Material). To use a consistent model structure for varied sizes of AuNPs, the membrane-limited model was chosen as the final model. Membrane-limited model predictions of concentrations of different sizes of AuNPs in various tissues were compared with the experimental data from Cho et al. (2010) (Figure 2). The model accurately predicted the measured concentrations of 13 nm AuNPs in the plasma, liver, spleen, kidneys, and lungs at all time points. Also, it properly predicted most of the data for 100 nm AuNPs in the plasma, liver, spleen, kidneys, and lungs, except for a slight ($2-fold) overestimation of the concentrations in the plasma, liver, and spleen at 4 h after injection (Figure 2). The overall regression coefficient between measured and predicted data is R 2 ¼ 0.97 ( Figure S1, Supporting Material), suggesting high goodness-of-fit of model calibration results. However, despite the adequate overall predictions, it should be noted that the model predicted a rapid increase followed by a rapid decrease of plasma concentrations of 100 nm AuNPs within 4 h after injection. However, there is a lack of experimental data at this early time point to confirm this, which requires further experiments to either validate this prediction or revise the model accordingly.

Model evaluation with independent data
The validity of the 13 nm PBPK model was first evaluated with data from Cho et al. (2009a) where mice were injected intravenously (iv) with 13 nm PEG-coated AuNPs (0.85 or 4.26 mg/kg). As shown in Figure 3 for the 0.85 mg/kg dataset, the model accurately predicted the concentrations in plasma and the amounts in spleen, and slightly underestimated (52-fold) the amounts in liver. A similar precision was obtained when comparing model simulations with the 4.26 mg/kg dataset ( Figure S2, Supporting material), indicating that the 13 nm model can simulate the pharmacokinetics of PEG-coated AuNPs at an acceptable accuracy across the dose range between 0.85 and 4.26 mg/kg.
To further evaluate the applicable size domain of the 13 nm model, model simulations were compared with other experimental data for similar sizes of PEG-coated AuNPs.  injected 111 In-labeled 20 nm PEG-coated AuNPs iv to healthy (blood pharmacokinetics) or tumor-bearing mice (tissue distribution). Liu et al. (2013) exposed tumor-bearing mice to 16 nm PEG-coated AuNPs via iv injection. Compared with the data from , the model successfully predicted the concentrations in liver and kidneys, but somewhat underestimated ($1.5-fold) the concentration in the blood and overestimated the spleen concentration by $10-fold ( Figure S3, Supporting material). Additionally, the concentrations in the spleen, kidneys, and lungs from Liu et al. (2013) were predicted adequately by the model, but the blood levels were slightly overestimated ($2-fold), and the liver concentrations were underestimated by $3-fold ( Figure S4, Supporting material). It should be noted that this model was calibrated with plasma data from healthy mice, whereas  and Liu et al. (2013) both measured the concentrations in the blood. It was shown that the concentrations of Au in the red blood cells and serum were different in mice after iv injection with 198 Au-radiolabelled mono-sulfonated triphenylphosphine-stabilized AuNPs and this difference was size dependent (Hirn et al., 2011). Therefore, the slight underestimation or overestimation of blood concentrations may be due to differential retention of different sizes of AuNPs between red blood cells and plasma. In addition, the overestimation of spleen concentration and underestimation of liver concentration may be attributable to the differences in the animal (healthy versus tumorbearing mice) and/or the type of AuNPs (13 nm versus 16 nm or 20 nm, non-radiolabelled versus radiolabelled).
The 100 nm model simulations were compared with the experimental data from  who iv injected 111 In-labeled 80 nm PEG-coated AuNPs to healthy or tumorbearing mice (Figure 3). The model precisely predicted the blood kinetics at earlier time points ( 1 h), but underestimated it at the later time points (!2 h) after injection. The liver and spleen levels were underestimated by the model, but the difference was very minor ($2-fold) and is considered acceptable according to the WHO model validation criteria (WHO, 2010). This slight underestimation may also be due to the differences in the animal and/or the size of AuNPs, as explained above.
Overall, there was a good correlation (R 2 ¼ 0.85) between model estimates and selected experimental data for model evaluation ( Figure S5, Supporting material), indicating that the present 13 nm PBPK model is applicable to the size range from 13 nm to 20 nm, and the 100 nm model has a suitable size range of 80-100 nm.

Model application: predictions of endocytic rates of different sizes of nanoparticles
The model-derived estimates of endocytic rates for 13 nm and 100 nm AuNPs in the liver and spleen of mice after iv injection (0.85 mg/kg) are shown in Figure 4. Endocytosis of the 100 nm AuNPs occurred immediately following injection, whereas endocytosis of the 13 nm AuNPs was delayed by $10 h, probably due to delayed activation of PCs by small sizes of AuNPs. The endocytic rates for 13 nm AuNPs were lower ($18-fold in the liver and $3-fold in the spleen) than those for 100 nm AuNPs.

Sensitivity and uncertainty analyses
Physiological and NP-specific parameters that had significant influence on selected key dose metrics are shown in Table 2. Body weight (BW) and plasma volume fraction (VplasmaC) had a high impact on all selected dose metrics regardless of the NP size or time frame of the study. On the other hand, cardiac output (QCC) greatly influenced the dose metrics of 13 nm, but not 100 nm, AuNPs. In addition, vascular space in the liver (BVL) had minimal effects on all dose metrics of 13 nm NPs, but it greatly affected liver dose metrics of 100 nm NPs during both earlier and later phases of the study. Similarly, vascular space in the spleen (BVS) considerably impacted spleen dose metrics of 100 nm NPs throughout the 168 h simulation period, and it also affected the spleen dose metric of 13 nm NPs during the later phase, but not the earlier phase. These data suggest that certain physiological parameters play an important role in the pharmacokinetics of AuNPs.
Similar to physiological parameters, the sensitivity of NPspecific parameters on the pharmacokinetics of AuNPs depends on the NP size and the time (Table 2). Specifically, distribution coefficient in the liver (PL) had a high impact on the 24 h liver AUC for 13 nm AuNPs, but had minimal effects on the 168 h liver AUC. On the other hand, PL had no influence on all selected dose metrics for the 100 nm AuNPs. Similar effects were observed for the distribution coefficient in the spleen (PS). These results suggest that the pharmacokinetics of 13 nm AuNPs is mainly controlled by transcapillary transport during the earlier phase, but is primarily regulated by endocytosis once it is activated during the later phase. On the contrary, 100 nm AuNPs are primarily regulated by endocytosis throughout the entire simulation period.
Among the three endocytosis-related parameters for each organ t, K max t represents the maximum uptake rate, whereas K 50 t and n t determine the time and the extent of activation of endocytosis, respectively. Our data showed that liver 24 h or 168 h AUC for 100 nm AuNPs was highly sensitive to K max liver , but not to the other two parameters, while liver 24 h or 168 h AUC for the 13 nm AuNPs were highly sensitive to K 50 liver and n liver , but not K max liver (Table 2). Similar effects of endocytosis-related parameters on spleen dose metrics were observed.
Local sensitivity analysis demonstrated that, among the highly sensitive NP-specific parameters, K 50 liver and n liver had the greatest influence on the 13 nm model output (i.e. liver concentration) and the effects were time dependent (Figure 5), while all other parameters had time-dependent moderate impacts on the model output. On the other hand, K max liver had the greatest impact on the selected 100 nm model output, followed by K bile , whereas all other parameters had low impacts.
The uncertainty designation for all highly sensitive physiological parameters and for PS was low, for PL it was medium, whereas for endocytosis-related parameters and biliary clearance rate it was high (Table 2).

Discussion
In the field of nanomedicine, there is a need to develop PBPK models that can integrate available data to gain deeper insights and make more definitive conclusions on the pharmacokinetics of  Cho et al. (2009a) and . Comparison of model predictions (solid lines) versus measured concentrations in the plasma (A) or measured amounts in the liver (B) and spleen (C) of healthy mice after iv injection with 0.85 mg/ kg 13 nm PEG-coated gold nanoparticles (AuNPs) (Cho et al., 2009a), or versus measured concentrations (means ± SD) of 80 nm PEG-coated AuNPs in the blood of healthy mice (D), or in the liver (E) and spleen (F) of tumor-bearing mice after iv injection ). The dose used in  was set to be 2 mg/kg based on Khlebtsov & Dykman (2011). NPs since there are hundreds of isolated pharmacokinetic experimental studies, but only a few integrating PBPK models . The present study successfully developed a PBPK model for different sizes of PEG-coated AuNPs in mice after iv injection. This model has been validated with multiple external datasets, indicating robust predictive capability. Besides being the first PBPK model for PEG-coated AuNPs, the mathematical description of endocytosis of NPs, the consideration of size dependence, and the model structure are all novel within the context of NP PBPK modeling.
Existing PBPK models simulate endocytosis of NPs using a linear equation, assuming one uptake pathway, i.e. either from the tissue (Li et al., 2014) or from the blood (Bachler et al., 2013(Bachler et al., , 2014, without considering size-dependent uptake. However, PCs can clear particles either from the blood, lymph, or tissue and endocytosis of NPs is dose and size dependent (Dykman & Khlebtsov, 2014). The present model tested both uptake pathways and assessed both the linear equation and the Hill function. We found that a linear equation could not simulate the rapid increase in the concentration of 13 nm PEG-coated AuNPs in the liver and spleen at 7 d after iv injection. Therefore, we used the Hill function to simulate NP endocytosis, which, we believe, is more physiologically realistic, as multiple in vitro studies have shown that uptake of NPs is saturable, and it is more appropriate to be described using the Hill function (Dykman & Khlebtsov, 2014;Zhang & Monteiro-Riviere, 2009;Zhang et al., 2011). PEG coating is considered as a useful strategy in the design of long circulating NPs because endocytosis of PEG-coated NPs might be avoided (Bazile et al., 1995) or can be substantially decreased (Dykman & Khlebtsov, 2014;Gref et al., 2000). It deserves further investigation to compare the endocytic rates of PEGcoated and non-PEG-coated AuNPs of similar sizes and doses within a PBPK model. However, currently there are no sufficient data that allow is to develop a PBPK model for non-PEG-coated AuNPs . Our modeling strategy for endocytosis of PEG-coated AuNPs can be extrapolated to other types of NPs, including non-PEG-coated AuNPs.
It is important to note that for larger PEG-coated AuNPs (80-100 nm), it is more appropriate to simulate endocytosis from the blood; while for smaller PEG-coated AuNPs (13-20 nm), AUCCA, AUCCL, and AUCCS represent area-under-the-concentration curve of gold nanoparticles (AuNPs) in the plasma, liver, and spleen, respectively. UN, uncertainty designation. L, M, and H stand for low, medium, and high uncertainty, respectively. QCC, cardiac output; BW, body weight; VLC, volume fraction of liver in the body; VplasmaC, volume fraction of plasma in the body; BVL, volume fraction of blood in liver; BVS, volume fraction of blood in spleen; PL, liver:plasma distribution coefficient; PS, spleen:plasma distribution coefficient; K max_liver or K max_spleen , maximum uptake rate constant in the liver or spleen; K 50_liver or K 50_spleen , time reaching half maximum uptake rate in the liver or spleen; n liver or n spleen , Hill coefficient in the liver or spleen; K bile , biliary excretion rate constant. Only parameters with at least one absolute value of NSC greater than 0.5 are presented. NSC values with absolute values higher than, equal to, or around to 0.5 are considered as highly sensitive and are highlighted in bold.  it is better to describe endocytosis from the tissue. Model simulations suggest that endocytosis of 100 nm PEG-coated AuNPs occurs rapidly and at much higher rates than 13 nm PEG-coated AuNPs, which has a $10 h delay, indicating size and time dependence (Figure 4). These results are consistent with previous studies that show plasma circulation time of PEG-coated AuNPs decreases with increasing sizes, e.g. plasma half-lives were 31.9 h and 7.3 h for 15 nm and 100 nm PEG-coated AuNPs, respectively (Chou & Chan, 2012;. Our results suggest that the difference in the pharmacokinetics between varied sizes of PEG-coated AuNPs may be mainly because of the different cellular uptake pathways and rates. By using a more detailed and physiologically realistic mathematical representation of endocytosis, the present study provides valuable insights into the size-and time-dependent pharmacokinetics of PEG-coated NPs. However, it should be noted that the present model was built on the arbitrary classification of two general pathways based on the cellular location. Other mechanisms, such as different endocytic pathways based on different molecular mechanisms (e.g. receptor-mediated versus receptor-independent pathways) with differential rates, may also contribute to what we observed. Therefore, in order to reveal the complete pharmacokinetic profiles of NPs and to further optimize PBPK model structure in vivo, future in vivo pharmacokinetic studies coupled to in vitro cellular uptake studies in major phagocytic cell types are warranted. In addition, it is recognized that many other factors, including the surface coating and biocorona formation, affect the pharmacokinetics of NPs . Therefore, whether this size-and time-dependent cellular uptake mechanism applies to other types of NPs or not remains to be investigated. Sensitivity analysis results show that endocytosis-related parameters, distribution coefficients, and biliary and urinary excretion rates are highly influential parameters on selected critical dose metrics; the extent of sensitivity of each parameter is size and time dependent (Table 2). Among all these highly sensitive parameters, K 50 t and n t show the greatest impact on the dose metrics for 13 nm PEG-coated AuNPs, whereas selected dose metrics for 100 nm PEG-coated AuNPs were sensitive mostly to K max t ( Figure 5). These results suggest that pharmacokinetics of these PEG-coated AuNPs may be regulated primarily via endocytosis. Specifically, endocytosis of 100 nm large PEG-coated AuNPs may be triggered rapidly, so its extent of uptake primarily depends on the maximum uptake rate (K max t ). On the other hand, endocytosis of smaller 13 nm PEG-coated AuNPs may be delayed and mainly dependent on the extent and when endocytosis occurs, which is determined by K 50 t and n t . Therefore, different pharmacokinetic behaviors of varied sizes of PEG-coated AuNPs may be mainly due to different time courses of endocytosis. This model-based conclusion remains to be verified with additional experimental studies that actually measure the cellular uptake rates of AuNPs in various cell types.
The uncertainty analysis data suggest that uncertainty in the endocytosis-related parameters could add high uncertainty in the model predictions. Endocytosis-related parameters had to be estimated from a single study (Cho et al., 2010) due to lack of more time-intensive in vivo data for PEG-coated AuNPs. These parameters may be further refined once more comprehensive pharmacokinetic data for PEG-coated AuNPs become available. Alternative, endocytosis-related parameters derived from in vitro uptake assays using different types of PCs might be able to serve as input parameters to refine the PBPK model. Nevertheless, the successful validation of the present model with multiple independent datasets could increase confidence of these estimated parameters.
The sensitivity analysis reveals several physiological parameters that have high impact on selected critical dose metrics, and this effect is also size dependent (Table 2). These results are important because it suggests that minor differences in these parameters will result in similar extent of changes in the model outputs. Hence, when extrapolating animal study results to humans, it is critical to consider the difference in the values of these parameters between humans and laboratory animals, including mice, rats, dogs, and pigs. As shown in Table S4 (Supporting material), pigs would appear to be closest to humans regarding these highly sensitive physiological parameters, while the mouse values could be up to 3-fold different from humans, such as the blood volume fraction in the liver (BVL). Besides the general physiological parameter differences, there are differences in the anatomical architecture of vascular system and the immune function of macrophages among different species that need to be considered when performing interspecies extrapolation. Specifically, the mouse spleen capillary is non-sinusoidal with predominant blood flow through open-circulation routes where NP filtration is primarily regulated by the function of barrier cells; whereas in rats and humans, the spleen capillary is sinusoidal, where most of the blood flows through the open-circulation route with filtration at interendothelial cell slits (Cesta, 2006). The liver capillaries of mice, rats, and humans are all sinusoidal with open fenestrae, but the number fenestrae/mm 2 and the diameter differ among these species, e.g. the average number of fenestrae/mm 2 is 14, 20, and 20 in mice, rats, and humans, respectively (Braet & Wisse, 2002). There are also species differences in the plasma membrane receptor expression and functionality of macrophages (Gordon & Taylor, 2005), which play a significant role in the NP opsonization and endocytosis. However, despite these mechanistic differences, PBPK models provide the most promising framework to integrate these cellular and systemic physiological factors into whole-animal predictive models.
The model simulation results and the above-mentioned species differences suggest that pharmacokinetics of NPs could be very different between mice and humans, not only because of longer circulation time allowing for further biocorona formation in humans as predicted in our earlier studies (Riviere, 2013;Sahneh et al., 2015) but also due to inherent physiological, anatomical, and immune differences between the two species. Therefore, one should be very cautious in the interpretation of NP pharmacokinetic studies in mice and in the scaling from mice to humans. However, it should be noted that the majority of available pharmacokinetic studies for PEG-coated AuNPs are in mice ; this is the reason why a PBPK model was developed first in mice as a basis for future extrapolation to other species. Notably, this approach was successfully applied to extrapolate NP behavior to humans as demonstrated in a recently published PBPK model for titanium dioxide NPs in mice (Bachler et al., 2014). Since there is no one animal model that is anatomically, physiologically (especially cardiovascular system), and immunologically similar to humans, further pharmacokinetic studies in rats and pigs that are physiologically closer to humans will be needed for interspecies extrapolations using mechanistic PBPK modeling to elucidate the pharmacokinetics of this type of NPs in humans.

Conclusions
A PBPK model for 13 nm and 100 nm PEG-coated AuNPs was successfully developed. The model was validated with multiple independent datasets for the NPs of similar sizes, indicating reliable predicting capability. Model simulations provide insights into the pharmacokinetics of the NPs in vivo, the influence of physiological and physicochemical parameters, and the design of future studies. This model predicted a time-and size-dependent endocytosis profile for AuNPs that needs additional in vivo pharmacokinetic studies coupled to in vitro cellular uptake studies to verify or to inform and optimize NP PBPK model structure in vivo. The several highly sensitive physiological parameters highlight the needs for interspecies comparison and scaling studies and an argument for conducting NP pharmacokinetic studies in species that are physiologically more closely related to humans. This PBPK model describes the model structure selection and parameterization in great detail (Supporting material), and thus provides a general NP PBPK model development strategy and a sound foundation for extrapolating to other types of NPs and to other species.