Meta-Analysis of Nanoparticle Delivery to Tumors Using a Physiologically Based Pharmacokinetic Modeling and Simulation Approach.

Numerous studies have engineered nanoparticles with different physicochemical properties to enhance the delivery efficiency to solid tumors, yet the mean and median delivery efficiencies are only 1.48% and 0.70% of the injected dose (%ID), respectively, according to a study using a nonphysiologically based modeling approach based on published data from 2005 to 2015. In this study, we used physiologically based pharmacokinetic (PBPK) models to analyze 376 data sets covering a wide range of nanomedicines published from 2005 to 2018 and found mean and median delivery efficiencies at the last sampling time point of 2.23% and 0.76%ID, respectively. Also, the mean and median delivery efficiencies were 2.24% and 0.76%ID at 24 h and were decreased to 1.23% and 0.35%ID at 168 h, respectively, after intravenous administration. While these delivery efficiencies appear to be higher than previous findings, they are still quite low and represent a critical barrier in the clinical translation of nanomedicines. We explored the potential causes of this poor delivery efficiency using the more mechanistic PBPK perspective applied to a subset of gold nanoparticles and found that low delivery efficiency was associated with low distribution and permeability coefficients at the tumor site (P < 0.01). We also demonstrate how PBPK modeling and simulation can be used as an effective tool to investigate tumor delivery efficiency of nanomedicines.

Nanomaterials or nanoparticles (NMs or NPs) can be engineered to have different physicochemical and biological properties, such as different shapes, sizes, charges, and surface coatings, to provide a multifunctional platform for diagnosis and targeting therapy of various diseases, including cancer.[1−4] In particular, the enhanced permeability and retention (EPR) effect is purported to be one of the major mechanisms for passive retention of 10–200 nm NPs, due to the impaired lymphatic drainage and relatively high permeability of vascular endothelial cells in tumors.[4−6] Examples include, but are not limited to, self-assembled polymeric micelles[7] and liposomes[8] that can be encapsulated with anticancer drugs to enhance tumor targeting and on-site drug releasing. In contrast, active targeting is a strategy that counts on surface functionalization with a certain ligand,[9] peptide,[10] or coating[11] to facilitate the recognition and binding to the surface of tumors.

Numerous studies have been devoted to the design of nanomedicines with higher therapeutic indices, i.e., higher tumor targeting ability, longer circulating half-lives in blood, tumor accumulation, and lower systemic toxicity in normal tissues.[12,13] Despite these purported advantages, NP translation into clinical applications is limited partly due to a low tumor delivery efficiency of only 0.70% of the injected dose (%ID) as previously reported.[14] Subsequent quantitative analyses further revealed that <0.0014%ID intravenously (IV) administered NPs with active targeting moiety were delivered to targeted cancer cells in the tumor sites with only 2.0% of the cancer cells interacting with NPs.[15] These surprisingly low tumor delivery and cancer cell targeting efficiencies suggest the importance of examining key physicochemical and pharmacokinetic determinants of NP disposition within the tumor microenvironment. Additionally, these low delivery efficiencies were estimated using empirical noncompartmental and nonphysiologically based modeling approaches[14,15] that are inherently unable to be extrapolated to predict time-dependent kinetics of NP distribution to tumors/tissues across species with different physiologies and study designs nor to provide a mechanistic explanation of estimated NP delivery in tumor-bearing animals. These computational limitations in prior pharmacokinetic modeling approaches can be addressed by using a physiologically based pharmacokinetic (PBPK) modeling approach, which is a mechanism-based computational modeling method that simulates absorption, distribution, metabolism, and excretion (ADME) of substances in an organism.[16−19] A number of PBPK models have been developed for different types of NPs loaded with or without anticancer drugs.[17,20−30] However, none of these studies use a PBPK modeling approach to systematically analyze and compare NP disposition across hundreds of disparate data sets from tumor-bearing animals in order to obtain a more insightful and comprehensive conclusion on the key determinants of NP delivery efficiency to tumors.

To improve our understanding and reveal the critical factors in the systemic delivery of NM to tumors, the present study aimed to develop and apply a generic PBPK model for describing NM disposition in tumor-bearing mice and then employ this model to analyze tissue and tumor distributions of various types of NMs in order to identify key determinants of NM tumor delivery efficiency. In addition, we applied the model to simulate short-term (24 h), long-term (168 h), and time-dependent biodistribution to tumors to examine the time dependence of NM tumor delivery efficiency. Finally, we proposed a long-term strategy for future studies of nanomedicines to enhance the design of preclinical trials and to facilitate clinical translation of NMs with higher tumor delivery efficiency and optimal therapeutic index, with a focus on the role of PBPK modeling and simulation in this process.

Results and Discussion

PBPK Model Calibration and Simulation in Tumor-Bearing Mice

This study used the inclusion/exclusion criteria listed in Figure to screen the literature and identified 200 pharmacokinetic studies (376 data sets in total) for subsequent PBPK modeling and simulations with experimental protocols summarized in Table S1 in Supporting Information. A PBPK model for NMs in healthy mice (Figure A) was calibrated with the biodistribution data for venous plasma, lungs, liver, spleen, and kidneys in healthy mice injected IV with 13 nm gold (Au) NPs (AuNPs).[31] The predicted values were in good agreement with reported kinetic profiles for up to 168 h postdosing as indicated by the estimated coefficient of determination (R2) of 0.95 (Figure S1 in Supporting Information). The tumor-bearing PBPK models (Figure B) extrapolated from the PBPK model in healthy mice (Figure A) adequately simulated most of the tumor delivery kinetics for 313 out of 376 total data sets (83%) with R2 ≥ 0.75 or <10% difference in tumor delivery efficiency. This difference was estimated at the last sampling time point and calculated using the PBPK model (DETlast) (Table S2) versus a noncompartmental linear trapezoidal integration method (DETlast_PK) (Table S1).

Figure 1

Procedure, strategies, and inclusion/exclusion criteria for the literature search. Following the literature search from the databases of Cancer Nanomedicine Repository (CNR) and PubMed and application of listed selection criteria, 200 tumor-bearing mouse studies with a total of 376 data sets published from 2005 to 2018 were identified for subsequent PBPK modeling and simulation analyses.

Figure 2

Schematic diagram of PBPK models in (A) healthy and (B) tumor-bearing mice intravenously administered with AuNPs and various inorganic and organic nanomaterials (INMs and ONMs), respectively. Except plasma and brain, each compartment is divided into three major parts: capillary blood, tissue interstitium, and endocytic/phagocytic cells (PCs) or tumor cells (TCs).

Figure shows comparisons of PBPK model simulations with 29 representative data sets of tumor delivery kinetics in tumor-bearing mice following IV injection with different types of inorganic (I) or organic (O) NMs (INMs or ONMs). Additional information associated with performing the PBPK simulations as well as the estimated nanoparticle-specific parameter values are summarized in Tables S1–S2 with details provided in the Supporting Information Excel files E1–E4. Within 24–48 h after IV administration, most tumor kinetic profiles displayed distinctive uptake during the accumulation phases, indicating that NMs were delivered to solid tumors successfully and resided in the tumors for at least 24–48 h and then were eliminated from the tumors (Figure ).

Figure 3

Representative simulation results from the PBPK model in tumor-bearing mice intravenously administered with various types of INMs, including (A) and (B) gold,[167,179] (C) iron oxide,[176] (D) gadolinium (Gd)-calcium phosphate (CaP),[147] and (E) silica[157] NMs as well as ONMs, including (F) liposome,[121] (G) dendrimer,[35] (H) hydrogel,[36] (I) polymeric,[11] (J) single-wall carbon nanotube (SWCNT),[122] (K) ginseng extract,[195] and (L) anticancer drug 10-hydroxycamptothecin (HCPT)[102] NMs. Tumor tissue concentrations as presented in the y-axis are expressed in the units of percent of the injected dose (%ID),%ID/g, or μg/g according to units used in the original articles. R2 is the coefficient of determination. Uppercase letters P and A followed by each legend represent passive and active targeting strategies, respectively. Abbreviations: A, active targeting; AuNP, gold nanoparticle (NP); BSA, bovine serum albumin; F, folate; FA, folic acid; G4 dendrimer, generation 4 polyamidoamine dendrimer; GNC, gold nanocluster; IONP, iron oxide NP; NC/ND/NR, nanocube/nanodisc/nanorod; P, passive targeting; PEG, polyethylene glycol; PSMA, prostate-specific membrane antigen; RGD, arginine-glycine-aspartic acid peptide; SNP, silica NP; Tat, peptide; Zn, zinc(II).

In the present meta-analysis, we established a general tumor-bearing PBPK modeling framework that successfully simulated most of the available published data sets (i.e., 313 out of 376 data sets) for describing NM biodistribution to tissues/tumors in tumor-bearing mice for up to 168 h after systemic administration. These PBPK models allow one to simulate the maximum, short-term (24 h), long-term (168 h), and NM concentration–time profiles in the tumor that reflects the time dependency of tumor delivery efficiency; a perspective typically not available or difficult to obtain from traditional animal studies.

Effect of Parameter Sensitivity on Tumor Delivery

A local sensitivity analysis was performed based on one representative study[32] with sufficient measured data points conducted using AuNPs in tumor-bearing mice for up to 168 h in order to determine the contribution, as well as the importance of each parameter, to the short-term (24 h) and long-term (168 h) dose metrics of tissues and tumor. Specifically, positive values of the normalized sensitivity coefficient (NSC) estimates indicate that an increase in the parameter value would also increase the dose metrics, and vice versa. Table S3 lists all the calculated NSC values for highly influential parameters on short-term and long-term dose metrics of AuNPs distributed to tumors or organs/tissues, in which parameters with |NSC| ≥ 0.3 were considered sensitive.

For physiological parameters, cardiac output (QCC) had more impact on the 168 h dose metrics of liver, spleen, and kidneys than the 24 h dose metrics; while fractional cardiac output to spleen (QSC) was highly influential to the liver dose metrics at 24 h (Table S3). Other physiological parameters such as fractional volume of the body for plasma, liver, spleen, kidneys, and tumor (VPlasmaC, VLC, VSC, VKC, and VTC, respectively) and blood volume fraction of liver tissue (BVL) had significant contribution to the increase in the dose metrics for individual tissues with most |NSC| estimates of >0.5. NM-specific parameters, distribution coefficients for spleen, kidneys, and tumor (PS, PK, and PT) were identified as highly influential parameters that affected both 24 h and 168 h biodistributions in spleen, kidneys, and tumor with estimated |NSC| values ranging from 0.7 to 1. NM-specific parameters associated with endocytic/phagocytic and tumor cell uptake within particular tissues had a greater impact on the dose metrics for liver, spleen, kidneys, and tumor estimated at 168 h than on the 24 h dose metrics.

Our sensitivity analysis suggests that tumor delivery efficiency of NMs was highly sensitive to the following parameters; PT (i.e., tumor tissue:plasma distribution coefficient), Kmax, (i.e., maximum uptake rate constant by tumor cells in the tumor tissue), and VTC (i.e., volume fraction of tumor tissue in the body) (Table S3). We further explored the relationship between tumor delivery efficiency and nanoparticle-specific parameters using AuNPs as an example (66 data sets). The results of rank-sum test (i.e., the estimated AuNPs-related parameters were separated into two groups with 33 data sets in each group according to their ranked values) suggest that lower distribution and permeability coefficients at the tumor site (PT and PATC) would hinder the targeting delivery of AuNPs and therefore lower the delivery efficiency (Table S4). The multivariable linear regression results indicate that PATC is a critical factor in affecting delivery efficiency of AuNPs to the tumor tissue. However, these results were based on estimated parameter values that had a degree of uncertainty. They remain to be optimized using more advanced statistical methods, such as the Bayesian approach with Markov chain Monte Carlo simulation[33] and verified experimentally with in vitro and/or in vivo studies. The cause and mechanism that lead to low tumor delivery efficiency might be due in part to multiple tumor physiological factors that are worth exploring in future studies.

Subgroup Analyses and Future Nanomedicine Design from the PBPK Perspective

Among all the estimates of tumor delivery efficiency, including DETlast, DE24 (tumor delivery efficiency estimated at 24 h post-IV administration), DE168 (tumor delivery efficiency estimated at 168 h post-IV administration), DEmax (maximum tumor delivery efficiency post-IV administration), and DETlast_PK; DE168 had the lowest median delivery efficiency of 0.35%ID, suggesting that by 168 h post-systemic administration, most of the NMs delivered to the tumor site had left the tumor site (Figure and Figure S2). There was a slight improvement in the mean DETlast value of 2.23%ID estimated from this study based on the data sets collected from 2005 to 2018 using a physiologically based approach, compared with the previously reported result of 1.48%ID estimated based on the literature from 2005 to 2015 using a nonphysiologically based approach[14] (Figure A and Table ). The median DETlast estimates from Wilhelm et al.[14] and the present study had no apparent difference (0.70% vs 0.76%ID), implying that the studied NMs can be delivered to solid tumors successfully at a median delivery efficiency of ∼0.7%ID. Surprisingly, using the PBPK approach, there was no apparent improvement in the mean tumor delivery efficiencies estimated from data sets after 2015 (2.33%ID) versus before 2015 (2.13%ID) (P = 0.65). There was no significant improvement of median tumor delivery efficiency from recent published data sets (after 2015) (P > 0.15). No statistically significant differences were observed in further comparison of mean delivery efficiency estimated using nonphysiologically and physiologically based approaches (i.e., DETlast_PKvs DETlast) (P = 1.00), indicating that both approaches generate similar tumor delivery efficiency estimates (Table S5).

Figure 4

Subgroup analyses on tumor delivery efficiencies estimated at the last sampling time point according to the original literature (DETlast) using our tumor-bearing PBPK model. Box-and-whisker plots of tumor delivery efficiency data (%ID) for different subgroups: (A) year, (B) targeting strategy, (C) type of nanomaterials (NMs), (D) inorganic NMs, (E) organic NMs, (F) shape, (G) hydrodynamic diameter, (H) ζ potential, (I) tumor model, and (J) cancer type. The boxes represent the 25th to 75th percentiles, and solid lines in the boxes indicate the median values. The pink dashed and solid lines denote the median and mean values of tumor delivery efficiencies derived from a previous study based on 193 published data sets from 2005 to 2015.[14] The green dashed and solid lines stand for the median and mean values of tumor delivery efficiencies derived from the present study based on 376 published data sets from 2005 to 2018.

Table 1

Summary of Tumor Delivery Efficiency Estimates for Different Types of Nanomaterials*

DETlast, DE24, and DE168 represent tumor delivery efficiency estimated at the last sampling time point according to the original pharmacokinetic study at 24 h and 168 h, respectively. DEmax is the maximum tumor delivery efficiency based on individual PBPK simulation.

DETlast_PK is the tumor delivery efficiency estimated at the last sampling time point according to the original pharmacokinetic study using the noncompartmental linear trapezoidal integration method as used by Wilhelm et al., i.e., area-under-the-tumor-concentration-curve (AUC) method.[14]

Median (mean).

Several trends were observed from subgroup analyses by comparing median estimates of tumor delivery efficiencies (DETlast) (Table ). First, active targeting strategies had slightly higher delivery efficiencies of 0.89%ID compared to 0.70%ID of passive targeting strategies (P < 0.01) (Figure B, Table , and Table S5). Second, INMs had superior tumor delivery efficiencies of 1.12%ID in comparison with 0.62%ID of ONMs (P < 0.05) (Figure C, Table , and Table S5). Furthermore, INMs composed of iron oxide and ONMs such as dendrimers tended to have higher delivery efficiencies of 2.80% and 7.96%ID, respectively (Figure D,E and Table ). Third, rod-shaped NMs had the highest median DETlast estimate of 1.62%ID followed by spherical and other shaped NMs of 0.74% and 0.58%ID, respectively (Figure F and Table ). Fourth, NMs with hydrodynamic diameters of <10 nm displayed the highest median DETlast estimate of 1.41%ID, followed by >200 nm, 10–100 nm, and 100–200 nm NMs of 0.94%, 0.75%, and 0.56%ID, respectively (Figure G and Table ). Fifth, NMs with neutral and positive surface charges possessed higher delivery efficiencies of 0.81% and 0.90%ID, respectively, compared to 0.47%ID of negatively charged NMs (Figure H and Table ). Finally, mice with orthotopically inoculated tumors, regardless of allografted or xenografted ones, had relatively higher delivery efficiencies of >1.00%ID, compared with approximate or lower than 0.70%ID for heterotopically inoculated tumor models (Figure I and Table ).

Tables S6 and S7 summarize the results of a one-way analysis of variance (ANOVA) test and simple linear regression for all 376 data sets and the selected 313 data sets that had confidence in the estimation of tumor delivery efficiency. The results for the curated confidently predicted data sets showed that the variables “tumor model” and “cancer type” had similar and significant contributions to the tumor delivery efficiencies of Au NMs and INMs (P < 0.05) (Table S7). By contrast, in addition to variable “cancer type”, variables related to NM physicochemical properties, including “core materials”, “shape”, “hydrodynamic diameter”, and “ζ potential” had significant influence on the overall tumor delivery efficiencies of ONMs (P < 0.05) (Table S7).

Table and Tables S8 and S9 list the multivariable linear regression models constructed for Au NMs, INMs, ONMs, and all NMs as well as the related statistical criteria for determining the significance and goodness-of-fit of each model. By excluding 63 data sets with weakly estimated (low confidence) tumor delivery efficiencies (R2 < 0.75 or >10% differences by comparing DETlast with DETlast_PK), the full regression models based on the selected confidently predicted 313 data sets with the same variables being included can better describe the response, e.g., log(DETlast), for gold NM (R2 = 0.62; P < 0.01), INMs (R2 = 0.51; P < 0.001), ONMs (R2 = 0.43; P < 0.001), and all NMs (R2 = 0.36; P < 0.001) (Table ), compared to the full regression models based on all the 376 data sets (Table S8). Overall, the estimated best regression models for the curated 313 confident data sets revealed that “cancer type”, “tumor model”, and “ζ potential” were critical factors in determining tumor delivery efficiency, regardless of the type of NMs studied (Table ).

Table 2

Multiple Linear Regression Results of Selected Models for the Log-Transformed Tumor Delivery Efficiency Estimated at the Last Sampling Time Point (log(DETlast)) Based on the 313 Confidently Predicted Datasetsa

**P < 0.01, and ***P < 0.001. R2, coefficient of determination; Adj-R2, adjusted R2; AIC, Akaike information criterion; BIC, Bayesian or Schwarz Bayesian information criterion; Type, inorganic or organic nanomaterials (INMs or ONMs); MAT, core material of INMs or ONMs; TS, targeting strategy; CT, cancer type; TM, tumor model; log(HD), log-transformed hydrodynamic diameter; ZP, ζ potential.

Based upon the Welch’s t-test, active targeting had statistically significant higher mean tumor delivery efficiency than passive targeting (Table S5). In contrast, the rank-sum test results suggest the median tumor delivery efficiency of active targeting was not statistically superior compared to passive targeting (Table S5). While most of the calibrated pharmacokinetic studies report that active targeted NMs have statistically significant higher tumor delivery efficiency, anticancer drug efficacy,[34−37] tumor cell toxicity,[8,10,35−42] and tumor growth inhibition[8,10,35−43] compared to passive targeting, some studies indicate either the opposite or no significant differences depending on the targeting strategy.[44−48] Comparing INMs with ONMs, both mean and median tumor delivery efficiencies of INMs were statistically significant higher than ONMs (Table S5). However, it must be noted that statistically significant differences do not always translate to biologically significant events or improved clinical outcome, and vice versa.[49−51] Biologically significant events or clinical improvements depend on multiple factors. Besides the delivery efficiency of the NP to the tumor site, other critical factors include the therapeutic potency and efficacy of the drug, the duration of exposure of the NP-carried drug at a concentration above the therapeutic threshold, as well as the toxicological profile of the drug.

This study confirms that smaller NMs with <10 nm hydrodynamic diameter can be delivered to tumors at a higher delivery efficiency compared to larger NMs (>10 nm). Significantly, there were Au[52,53] and dendrimer[54−56] NMs with measured hydrodynamic size <10 nm and estimated tumor delivery efficiency >5.00%ID (Table S1). Another noteworthy finding is that NMs with hydrodynamic diameter >200 nm have a relatively higher tumor delivery efficiency compared to NMs with ∼10–200 nm in hydrodynamic diameter. Several studies using >200 nm NMs with higher tumor delivery efficiency (>5.00%ID) support this observation (Table S1).[37,43,57−59]

Our results suggest that in addition to size, the geometry of NMs can modulate tumor uptake and in vivo disposition of NMs. In particular, a tumor delivery efficiency of 1.62%ID was observed for rod-shaped NMs which was greater than spherical (0.74%ID), plate-, or flake-shaped NMs (0.46%ID) as well as other geometries (0.58%ID). Several studies show that elongated nanostructures, compared to nanospheres, exhibit greater tumor accumulation and longer half-lives in blood circulation perhaps due to adherence to the endothelial cells lining the blood vessel walls that enhance the site-specific delivery.[43,60,61]

Our analyses show that administered NMs with positive (>10 mV) and near-neutral (−10 to 10 mV) surface charges have similar tumor deliver efficiencies (0.81–0.90%ID), but are higher than the negatively charged NMs (0.47%ID). Surface properties, such as charge (ζ potential), play a critical role in the type and magnitude of biomolecule (e.g., proteins, lipids, carbohydrates) adsorption that results in the formation of a biocorona, which in turn influences the pharmacokinetics, biodistribution, and cellular uptake of systemically administered NMs.[62−69] Surface charge of administered NMs can affect not only the above-mentioned site-specific extravasation[10,70−72] but also subsequent cancer cell type-specific internalization, e.g., higher internalizations of positively charged NMs have been reported in several cancer cell types, compared with their negatively charged counterparts.[73−76] Not surprisingly, our regression analyses reveal that in addition to NM ζ potential, heterogeneity in tumor physiology (i.e., cancer type) is an important factor in determining efficiency of tumor cell targeting, which is consistent and supported by previous studies.[7,73,77]

Similar to the nonphysiologically based study by Wilhelm et al.[14] that covered the time period from 2005 to 2015, the current PBPK model-based study suggests that there was no statistically significant improvement in NM tumor delivery efficiency before and after 2015 (up to year 2018). Compared to the results reported by Wilhelm et al.,[14] several key findings and major differences from the present study have been summarized in Table . Briefly, the current study employs a physiologically based approach (i.e., PBPK modeling) that is capable of estimating time-dependent tumor delivery efficiency and generates models that are extrapolatable across administration routes and doses as well as species; a computational approach that is not possible using traditional noncompartmental pharmacokinetics as used by Wilhelm et al.[14] Besides “cancer type” and “tumor model” as critical factors affecting tumor delivery efficiency reported by the previous study,[14] we have identified another critical factor “ζ potential”. These differences may be due in part to the incorporation of recent published data sets after 2015, different categorization methods used in subgroup and statistical analyses, and different tumor delivery efficiency estimates between the previous study[14] and the present study. Other statistical methods, including principal component analysis-based regression[78] or random forest regression[79] may be worth exploring in the future to determine the potential key determinants of tumor delivery efficiency.

Table 3

Summary* of the Differences and Major Findings between Current Study and the Previous Study by Wilhelm et al.(14)

DETlast_PK is the tumor delivery efficiency estimated at the last sampling time point according to the original pharmacokinetic study using the noncompartmental AUC approach as used by Wilhelm et al.[14] DE24, DE168, and DETlast represent tumor delivery efficiency estimated at 24 h, 168 h, and the last sampling time point, respectively, according to the original pharmacokinetic study. DEmax is the maximum tumor delivery efficiency based on the individual PBPK simulation.

The mean value of 1.48%ID was calculated from 193 included data sets based on the reported tumor delivery efficiencies by Wilhelm et al.[14] and the Cancer Nanomedicine Repository (CNR) database. There were 238 data sets in the CNR database at the time of the present study, 232 of which were analyzed and reported in the Wilhelm et al.[14] paper. One additional study containing 6 data sets was uploaded into the CNR database after the publication of the Wilhelm et al.[14] paper. After excluding studies due to lack of information and/or did not meet the criteria for subsequent PBPK modeling and simulation as described in Figure , there were 193 data sets from the CNR database being included and analyzed.

In the past several decades, many small molecule cancer drugs have been approved and used to treat different cancers successfully.[80] In general, the potential weaknesses of cancer chemotherapy include its low bioavailability, high-dose requirements, adverse side effects, low therapeutic indices, development of multiple drug resistances, and nonspecific targeting.[81] In the field of small molecule cancer chemotherapy, the pharmacokinetic focus is on improving overall drug bioavailability. For most of the approved cancer chemotherapeutic drugs, bioavailability data via oral route are available. However, the specific information on the “delivery efficiency” of a small molecule cancer drug to the tumor site is lacking, which is quite different from the delivery efficiency to the tumor site of NPs that is commonly reported in NM literature. In theory, we can calculate the “delivery efficiency” of a small molecule drug to the tumor site provided that we have the complete pharmacokinetic data in the plasma and tumor. However, this would be an entirely new investigation whose goals would be to confirm whether NP-based drug formulations outperform conventional small molecular chemotherapy.

Our PBPK model successfully simulates the majority of the tumor delivery kinetics for various types of NMs in tumor-bearing mice. Yet, there are several challenges and limitations to implement this PBPK computational framework in order to more appropriately describe NM disposition to the tumor and tissues in tumor-bearing animals. First, the model fails to capture the kinetics in the tumor for several NMs (63 data sets). This may be a result of many factors, including an insufficient number of and improperly timed experimental data points, unobvious uptake or release phases of the tumor kinetic profiles, a sudden increase in the tumor uptake phase and/or abrupt decrease following IV injection, or extended retention of delivered NMs until later time points. The exact mechanisms of these phenomena remain to be investigated, prohibiting our inclusion of these mechanisms in our general model. Importantly, the present study has employed a strict standard to evaluate model simulation results via both qualitative (i.e., visual inspection of the simulated and measured kinetic profiles) and quantitative evaluation (i.e., linear regression analysis and quantitative comparison between DETlast and DETlast_PK). In addition, a review of the curated confidently estimated data sets indicates that such a mechanism-based PBPK modeling approach could be used as a screening tool to select well-behaved data sets for further analyses. Finally, these analyses were conducted in mice, and integrating data sets from other laboratory animals would ultimately help improve NM designed for humans.[23]

Some limitations and challenges to the current approach remain to be addressed. Our PBPK model does not simulate the process of biocorona formation that dynamically alters the NMs’ properties and disposition[64,68] within an animal as we have previously reported.[66,82,83] PBPK models incorporating biocorona formation kinetics of NMs would better describe NM-specific tumor dosimetry and remain to be implemented. Moreover, the tumor vasculature, extracellular matrix, and microenvironment are highly variable in animals and humans bearing solid tumors. It is necessary to develop a nanomedicine-specific and tumor-bearing PBPK model in humans incorporating drug/nanomedicine pharmacodynamics[84,85] to describe individual tumor progression,[86] anticancer drug releasing pharmacokinetics,[19,22] and therapeutic response after receiving chemotherapy, i.e., PBPK/PD, to facilitate designs of the next-generation nanomedicines with optimal therapeutic index.[87−91]

To address these aspects, we propose a long-term integrative computational strategy from the PBPK perspective that could aid in the design of nanomedicine studies (Figure ). Specifically, by implementing our tumor-bearing PBPK models for a specific type of NMs coupled to constructed regression models, researchers may be able to estimate the tumor delivery efficiency of administered NMs and/or propensity of loaded anticancer drugs. By comparing to the specific type of NMs in our database, nanoscientists performing fundamental bench work would be able to modify their synthesized NMs with desired physicochemical properties and optimize the therapeutic index of NMs with higher drug loading capacity and tumor delivery efficiency with minimum systemic toxicity (left panel in Figure ). The present study provides a ready-to-use PBPK modeling framework calibrated by hundreds of pharmacokinetic data sets in mice bearing various tumor types. These informative and organized data sets for NMs that have been tested in the literature are provided in the Supporting Information Excel files E1–E4. The representative model code and all model parameters are provided in Supporting Information. This PBPK modeling framework may be extrapolated to rats, dogs, monkeys, and humans to gain more insight into determining the optimal dose with minimum side effects and systemic toxicity for designing personalized and optimized cancer therapy for preclinical trials as well as for individual cancer patients undergoing clinical trials (right panel in Figure ).

Figure 5

Proposed long-term strategy in facilitating the design of future nanomedicines and translation from preclinical studies to clinical applications, and the role of PBPK modeling and simulation approach in this field.

In view of the generally slow progress and limited success in translating nanomedicine into clinical applications, integrative PBPK modeling and systems biology analyses are needed and would be facilitated if the authors always reported physicochemical (e.g., size, shape, ζ potential, etc.) and biological (e.g., pharmacokinetic, biodistribution, and (cyto)toxicity data) properties of their synthesized NPs and details in experimental protocols (e.g., dose in the units of mg or mg/kg and tumor concentration in the units of μg/g tumor,%ID, or %ID/g tumor) as well as raw concentration–time data.[92] Only extensive interdisciplinary communication of reproducible nanomedicine design aiming at solving a clinical problem, in combination with open source repository documenting transparent experimental details provided by the authors that meet the “minimum information reporting” standard,[92] would possibly realize successful clinical translation.[89]

Conclusions

While one can be impressed with the dramatic advancement of nanotechnology and nanomedicine over the past decades, we must acknowledge that our current progress in translating cancer nanomedicine research into clinical application was slow before 2015, and in the ensuing years after 2015 (up to 2018). This continues to bring into question the validity of the EPR effect which is based on preferential tumor retention. A thorough understanding of NM-tumor interaction is necessary for facilitating clinical development of cancer nanomedicines with NM toxicity screening and safety evaluation, followed by manufacture of favorable, contaminant-free, reproducible, and scalable NMs. To accomplish this goal, a modeling framework that factors in species-specific physiological disposition is essential to extrapolate both from in vitro to in vivo as well as from animals to humans. This study explores tumor delivery efficiencies using data from hundreds of nanoplatforms in such a PBPK modeling and simulation framework to identify influential factors influencing tumor delivery kinetics. We hope this study and the computational approach facilitate the future design of cancer nanomedicines and improve clinical translation from bench work to bedside.

Materials and Methods

Data Source and Model Structure

In an elegant analysis, Wilhelm et al.[14] estimated tumor delivery efficiencies (%ID) using a noncompartmental linear trapezoidal integration pharmacokinetic method, which estimates the area of delivery efficiency under the concentration–time curve of the tumor (AUC). This statistical moment approach was applied to over a hundred publications which identified key factors that influenced tumor delivery.[14] The relationships between low tumor delivery efficiencies (<1.0%ID) and multiple factors were reported, including physicochemical properties of NMs, targeting strategies, cancer types, etc. In addition, they constructed the open source database of Cancer Nanomedicine Repository (CNR) that included the tumor delivery efficiency, related biological information, and physicochemical properties of delivered NMs from 118 publications (http://inbs.med.utoronto.ca/cnr/). On the basis of their work, we included published tumor delivery studies since 2015.

The criteria for whether or not to include the latest publications in the current analysis are summarized in Figure . Specifically, relevant tumor delivery studies were selected from the databases of CNR (studies published between 2005 and 2015) and PubMed (from January first, 2015 to September fourth, 2018) for further computational analyses, including PBPK model calibration and simulation, sensitivity analysis, and multivariable linear regression analysis. The literature search was conducted using the following keywords: nanoparticle delivery, nanomaterial delivery, biodistribution, pharmacokinetics, mice, rats, and “tumor or tumour”. In brief, this study only included published articles suitable for PBPK modeling, i.e., studies reported tumor concentration data in more than or equal to 3 sampling time points in units of μg/g tumor, %ID, and %ID/g tumor from tumor-bearing rodents following IV administration in convertible dose units of mg or mg/kg. IV is the major route of administration for nanomedicines and provides a better estimate of the actual dose delivered to the systemic circulation by eliminating potentially confounding absorption kinetics. Tumor delivery studies conducted with tumor-bearing rats were excluded due to insufficient number of studies compared to those in tumor-bearing mice (i.e., 10 out of 393 studies). Our analysis also excluded studies of which the corresponding authors did not respond to our requests for essential data and experimental details that are needed for PBPK analysis. In total, there were 376 data sets from 200 NP tumor delivery studies in tumor-bearing mice that were included for PBPK simulations and analyses (Table S1).[7−11,22,32,34−48,52−60,70,71,73,76,93−257]

There were two phases in PBPK model calibration and simulation in the present analysis: Phase I was to establish a permeability-limited PBPK model to simulate the biodistribution of AuNPs to organs/tissues following IV administration in healthy mice based upon our previously developed PBPK model (Figure A), and phase II was to extrapolate the healthy mouse PBPK model to a tumor-bearing PBPK model. This generic model structure was then employed to simulate the delivery of various INMs or ONMs to solid tumors after IV injection in tumor-bearing mice (Figure B). Specifically, the concentration data in healthy mice after IV injection of AuNPs were obtained from Cho et al. for PBPK model calibration.[31] In brief, 6-week old male BALB/c mice were injected with polyethylene glycol (PEG)-coated AuNPs (4, 13, or 100 nm) via the tail vein at 0.85 mg/kg body weight. Mice (n = 9/group) were euthanized at 0.5, 4, 24, 168 h and 1, 3, and 6 months post-injection at each sampling time point, and their plasma and tissue samples (e.g., lungs, liver, spleen, and kidneys) were collected and analyzed for Au concentrations using inductively coupled plasma-mass spectroscopy (ICP-MS). The 13 nm AuNP mouse pharmacokinetic data were selected to calibrate the healthy mouse PBPK model because of the sufficient sampling time points.

This model structure was based upon our recently developed PBPK model for AuNPs[20] with minor modifications to incorporate muscle and tumor compartments, respectively, for describing the biodistribution of NMs following IV administration in healthy and tumor-bearing mice. Specifically, the permeability-limited PBPK model for healthy mice contained eight compartments, including plasma, lungs, liver, kidneys, spleen, brain, muscle, and remaining tissues (i.e., pooled other tissues) (Figure A). Since tumor cells were inoculated subcutaneously into nude mice according to most tumor delivery studies, this present study included muscle as an additional compartment in the PBPK schematic (Figure B). Moreover, the present model considered uneven distribution between capillary blood and tissue, membrane-limited transcapillary transport, as well as nonlinear endocytic uptake and first-order exocytic release of administered AuNPs in order to describe the permeability-limited pharmacokinetics and tissue distribution of NMs in healthy mice. Besides plasma and brain, using standard PBPK modeling practice, all compartments were divided into three subcompartments: capillary blood, tissue interstitium, and endocytic/phagocytic cells (PCs) (Figure A). Similarly, to describe the NM distribution to tissues as well as to tumor microenvironment in tumor-bearing mice, a tumor compartment was subcompartmentalized as capillary blood, tumor tissue interstitium, and tumor cells (TCs) (Figure B). Pharmacokinetic data from both healthy and tumor-bearing mice were extracted using the WebPlotDigitizer (Version 4.1, Austin, TX, https://apps.automeris.io/wpd/) and provided in Supporting Information Excel files E3 and E4. Other experimental details extracted from original pharmacokinetic studies were documented and/or tabulated in Supporting Information Excel files E1 and E2 as well as Table S1.

Main Mathematical Description of the Model

Based on our earlier study,[20] endocytic uptake of NMs was more accurately described by the Hill function. Therefore, we employed the Hill function to describe the endocytic/phagocytic or tumor uptake of NMs as expressed in the following equation:where t is the simulation time (h), Kup,(t) is the uptake rate constant by endocytic/phagocytic or tumor cells in the tissue interstitium i at a particular time t (per h), Kmax, is the maximum uptake rate constant (per h), t50, represents the time reaching half-maximum uptake rate (h), and n is the Hill coefficient (unitless).

Furthermore, the distribution of NMs among subcompartments of capillary blood, tissue interstitium, and PCs/TCs in lungs, spleen, kidneys, muscle, remaining tissues, and tumor tissues are described using the following equations, respectively:where ABlood, and AT, represent amounts of NMs in the capillary blood and interstitium of tissue i (μg), APC, and ATC are amounts of NMs being taken up by PCs or TCs (μg), Q is the regional blood flow to tissue i (L/h), Ca is the concentration of NMs in the arterial plasma (μg/L), CV is the concentration of NMs in the venous plasma in tissue i (μg/L), PA is approximated to the product of permeability coefficient between capillary blood and tissue membrane (PAC, unitless) and regional blood flow Q (L/h) of tissue i (L/h), CT, is the concentration of NMs in the tissue interstitium i (μg/g), P is the tissue:plasma distribution coefficient for tissue i (unitless), and Kre, is the release rate constant of NMs by PCs or TCs to the tissue interstitium i (per h). Similar equations were used to simulate the distribution of NMs among different subcompartments in the liver. The only difference in the liver compared to other organs is that the liver has Kupffer cells which are liver macrophages that can directly phagocytize NMs from the capillary blood subcompartment since they are located in the liver sinusoids and directly exposed to blood.[258−260]

PBPK Model Calibration and Evaluation

All simulations and model calibration were performed using Berkeley Madonna (Version 8.3.23.0, University of California at Berkeley, CA, USA) to obtain visually reasonable fits to the pharmacokinetic data from healthy and tumor-bearing mice. The model for healthy mice was calibrated with the pharmacokinetic data in mice up to 168 h after IV administration with 13 nm PEG-coated AuNPs.[31] To develop the most parsimonious model following IV administration, most physiological parameters were kept consistent with the literature.[261,262] For physicochemical parameters, distribution, and permeability coefficients, values obtained from our previous study[20] were used as references for further optimization for the present PBPK models in healthy mice using visual fitting and the Curve Fitting Module in Berkeley Madonna. Similarly, other physicochemical parameters were estimated using visual fitting and the Curve Fitting Module, including cellular uptake and release rates by PCs of the lungs, liver, spleen, kidneys, muscle, and remaining tissues as well as excretion rate constants to establish the PBPK model in healthy mice.

The healthy mouse model was expanded to include a tumor compartment to describe NM biodistribution following IV injection in tumor-bearing mice. Except for tumor-related parameters, all physiological as well as NM-specific parameters for other organs/tissues remained the same as those used in the healthy mouse PBPK model. The tumor-related parameters included fractional blood flow to tumor (QTC), fractional tumor weight (VTC), fractional blood volume in the tumor tissue (BVT), distribution coefficient (PT), permeability coefficient to the tumor tissue (PATC), maximum NM uptake rate constant for TCs (Kmax,T), time reaching half-maximum NM uptake rate for TCs (t50,T), Hill coefficient for the uptake by TCs (nT), and NM release rate constant from TCs to the tumor interstitium (Kre,T). Specifically, VTC was estimated based on the original tumor-bearing study. QTC and BVT were assigned with very small numbers initially (∼0.02–0.03) as the blood flow inside and around the tumor only accounts for a small fraction of the cardiac output and optimized afterward using visual fitting and the Curve Fitting Module. Other tumor-related parameters (PT, PATC, Kmax,T, t50,T, nT, and Kre,T) were obtained similarly. The calibrated PBPK models in tumor-bearing mice were then used to predict short-term (24 h) and long-term (168 h) delivery efficiency and kinetics for each NM to the tumor. The PBPK model example code from a representative study is provided in the Supporting Information.

The performance of the PBPK models for both healthy and tumor-bearing mice was evaluated by comparing model simulations with measured pharmacokinetic data based on the criteria described in the World Health Organization (WHO) guideline.[263] Specifically, if the simulated tissue distribution matched the measured kinetic profile visually and the simulated values were within a factor of 2 of the measured values, the model was considered acceptable. In addition, the goodness of fit between log-transformed values of measured and predicted tissue distribution (i.e., μg/g tissue, %ID, or %ID/g tissue) were further evaluated by analyzing coefficient of determination (R2) using simple linear regression with an R2 ≥ 0.75 of predicted over measured values considered adequate.

Sensitivity Analysis

Local sensitivity analyses were conducted to identify highly influential parameters (e.g., physiological and NM-specific parameters) governing the overall tissue distribution and tumor delivery efficiency after IV administration of various NMs. One representative study by Karmani et al.,[32] with an adequate simulation of AuNP distribution to tumor compartment for up to 168 h (R2 > 0.95), was selected to estimate the dose metrics for target organs or tissues. Specifically, each parameter (p) was increased in increments of 1% and the corresponding AUC (%ID × h) of NMs in the venous plasma, liver, spleen, kidneys, and tumor were computed at 24 h and 168 h after IV administration. Highly influential parameters identified by 24 and 168 h AUC estimates were compared to determine the differential effects of parameters on the short-term and long-term internal dose metrics, respectively. NSC was calculated by dividing the relative change in AUC (dAUC/AUC) with the relative change in each parameter (dp/p).[20] Parameters with at least one absolute value of the calculated NSC around or greater than 0.3 (i.e., |NSC| ≥ 0.3) were considered sensitive.

Subgroup Statistical Analyses

To compare with the delivery efficiency results reported by Wilhelm et al. from 2005 to 2015,[14] the present study used different approaches to calculate the delivery efficiency for each data set, including (1) the maximum tumor delivery efficiency of the simulated delivery efficiency kinetics using the calibrated PBPK models (DEmax, %ID); (2) AUC-based tumor delivery efficiency estimated based upon PBPK simulations, i.e., AUCs of tumor delivery efficiency estimated at 24 h, 168 h, and the last sampling time point (Tlast) in accordance with the original study (%ID × h) and then divided by time (h) to generate DE24, DE168, and DETlast (%ID), respectively; and (3) the tumor delivery efficiency estimated at Tlast based upon original tumor pharmacokinetic data (DETlast_PK, %ID) using a noncompartmental linear trapezoidal integration method as used by Wilhelm et al.[14] Subsequently, all tumor delivery efficiency data were categorized into varied subgroups, including the year of publication, physicochemical properties of administered NMs (i.e., shape, size, and ζ potential), targeting strategy, type and core materials of delivered NMs, method of inoculation, as well as inoculated cancer cell types.[14] The ζ potentials measured at pH 7.4 of < –10 mV, −10 to 10 mV, and >10 mV were categorically defined as negative, neutral, and positive, respectively.[14] Following subgrouping, unpaired parametric t-test assuming unequal variances (Welch’s t-test) and nonparametric rank-sum test (Mann–Whitney test) were implemented to examine, respectively, whether a significantly higher mean and median tumor delivery efficiency can be observed by comparing two subgroups. Welch’s t-test and Mann–Whitney test were performed using GraphPad Prism (Version 6.05, GraphPad Software Inc., La Jolla, CA).

Multivariable linear regression was used to determine the potential effects of various physicochemical characteristics, including the hydrodynamic size, ζ potential, type of NMs, targeting strategies, cancer types, and tumor models on tumor delivery efficiency. Tests for normality were performed first for both delivery efficiencies with or without log transformation to identify dependent variables (i.e., delivery efficiencies) with (or at least more like) normal distributions. Based on the results from the normality test, the log-transformed delivery efficiencies were employed for a one-way ANOVA, simple linear regression, and multivariable linear regression analyses. One-way ANOVA and simple linear regression were performed prior to multivariable linear regression to examine the significance (P < 0.05) of categorical and continuous variables, respectively. Statistical analyses, including normality test, one-way ANOVA, simple and multivariable linear regression, were conducted using R language (Version 3.6.0).

Specifically, variables including physicochemical properties of delivered NMs (e.g., log-transformed hydrodynamic diameter (log(HD)), ζ potential (ZP), and shape) as well as type of NMs (ONMs or INMs) (Type), core material of ONMs or INMs (MAT), targeting strategy (TS), cancer type (CT), and tumor model (TM) were considered for one-way ANOVA, simple, and multivariable linear regression analyses. In total, we report 80 regression models from 5 log-transformed tumor delivery efficiencies (DETlast, DE24, DE168, DEmax, and DETlast_PK) for 4 different types of NMs (Au NMs, INMs, ONMs, and all NMs). For a specific response of delivery efficiency from a particular type of NMs, we have included 4 regression models, i.e., full, best, full confident, and best confident models. The full and best regression models were based on all the 376 data sets. The full confident and best confident regression models were computed based on delivery efficiency data that were confidently predicted (only 313 data sets), e.g., the simulations visually matched the measured tumor distribution kinetics reported by original literature and with either R2 ≥ 0.75 or <10% difference by comparing DETlast to DETlast_PK. This difference was estimated at Tlast (up to 168 h) and calculated using the PBPK model versus a nonphysiologically based linear trapezoidal integration method. The best and best confident regression models were determined using a stepwise approach whenever the values of Akaike information criterion (AIC) reached the smallest among all established regression models. Other critical statistical criteria, including R2, adjusted R2 (Adj-R2), P-value, and Bayesian information criterion (BIC) obtained from simulation outputs were used to help determine the appropriateness in selected regression models for predicting tumor delivery efficiency. Moreover, to explore the potential association of NM tumor delivery efficiency with tumor microenvironment and NM properties, a multivariable linear regression analysis was conducted to explore the relationship between log-transformed tumor delivery efficiency (DETlast) and 4 nanoparticle-specific parameters at tumor site (PT, PATC, Kmax,T, and Kre,T) using AuNPs as a case study. This analysis was done only for AuNPs to avoid confounding by different types of NMs.