Drug encapsulation into nanocarriers makes it possible to control the pharmacokinetic properties of therapeutics agents, including their release and circulation half-life, while limiting their interactions with healthy tissues [51]. The PBPK modeling approach has existed for many years and has been used to describe time-dependent concentration profiles of substances in various organs of a living system and interspecies scale-up [52]. This approach divides the body into anatomical compartments, interconnected by body fluid systems. The number and nature of the compartments, as well as the complexity of PBPK models, are determined depending on the scientific task and physiological characteristics of the modeled organism. For example, Gilkey et al. considered a five-compartment PBPK model including liver, spleen, kidneys, plasma, and “other” in accordance with the biodistribution of fluorescently labeled NPs in mice used for the controlled delivery of dexamethasone in the therapy of acute lymphoblastic leukemia [53]. While to investigate the in silico effects of NP properties, tumor-related variables, and individual physiological differences on systemic bioavailability, mononuclear phagocytic system sequestration, tumor delivery, and excretion of NPs in rats, Dogra et al. developed a tumor-compartment-bearing PBPK model consisting of 12 compartments of interest: brain, heart, lungs, plasma, liver, spleen, gastrointestinal tract, kidneys, muscle, ‘others’, lymph nodes, and a facultative tumor [54]. Models designed for primates (including humans) tend to have even more compartments [55,56,57,58]. Thus, Perazzolo et al. investigated a whole-body PBPK model for the anti-HIV drug-loaded NPs in nonhuman primates [56]. The model describes the uptake of the injected dose from the subcutaneous site to the adjacent lymphoid depot, passage through the lymph nodes and throughout the lymphatic network, and its subsequent passage into the blood circulation. For this, the model includes 23 compartments: subcutaneous injection site, two adjacent-to-injection lymphoid tissue compartments, thoracic duct, compartments of regional lymph node (cervical node, hilar node, axillary node, mesenteric node, inguinal node), vein, artery, head and neck, lungs, upper body, kidneys, small and large intestines consisting of tissue and mucosa, spleen, liver, lower body and tail, as well as the rest of the body.

Typically, each compartment in PBPK models can be described in two ways. The first is referred to as permeability-limited model (also known as the diffusion-limited or membrane-limited model), in which tissue cell membranes are considered as the main barriers to drug (nanotherapeutic) penetration. The other one is referred to as perfusion-limited model (also known as the flow-limited model) and considers blood perfusion as the only limiting step for drug (nanotherapeutic) penetration through the tissue cell membrane [24,25,59,60]. It is generally accepted that the permeability-limited model is more appropriate for modeling NPs [24,60,61].

It should be noted that the PBPK modeling techniques for NPs and small molecules could be quite different. Traditional route-to-route extrapolation for small molecules is typically performed by using administration of route-specific parameters and keeping other chemical-specific parameters the same [62]. However, unlike small molecules, upon contact with different body fluids following different routes of administration, NPs will be covered by different proteins and other biomolecules, producing different biomolecular coronas and resulting in different patterns of biodistribution. Thus, a recent study by Chou et al. [62] clearly demonstrated that the traditional approach for small molecules is not applicable to NPs, and multiroute PBPK models for NPs should be developed using route-specific data.

The application of mathematical models to nanomedicine is based on breaking down their transport in different discrete and simpler phases which are eventually modeled separately. The sum of these contributions will eventually provide an overall picture of the phenomenon with fundamental hints of prediction that will allow better optimization of cancer treatment and, in this case, NP synthesis [63]. The model developer literally dissects the voyage of the particles in different phases depending on the function of the organ in which they circulate, the barrier function of various elements (i.e., tumor vasculature vs. healthy vasculature), and the characteristics of the carrier (shape, size, surface charge, targeting properties, release rate, etc.). All these phenomena and characteristics are described by mathematical functions that summarize the physical and chemical variabilities that characterize the biological barriers and represent the various parameters under consideration and the core of the modeling. The model is then run with appropriate software, which, by combining different parameters and solving respective equations, provides an estimation of particle location in various organs, including the tumor. The simulated dynamics can eventually be compared with experimental data to optimize the reliability of the calibrated model. When close collaboration between wet-lab biologist and computational scientist is not possible, the latter can extract experimental data from the published literature using “auxiliary” software that allows for the digitization (image-to-number conversion) of images and graphs from the scientific literature and then publish the model under copyright permission. All these transitions are eventually followed up with statistical analysis to evaluate the significance of the data and validate the model. After validation, the model can be used to predict other situations, nanocarriers, and being applied to other tumor diseases. However, an extrapolation of the modeling results from one species to another should be conducted very carefully due to parameter uncertainty. For instance, a collective fitting process of parameter values to data can provide multiple parameter sets that give similar model dynamics reproducing experimental observations. In this case, identifiability of fit parameters is the crucial step of the model analysis for reliable application of PBPK models and a more confident translation of their parameters into new experimental settings [64].

It is worth noting that reproducibility and replicability crisis in research areas has also affected systems biology research such as PBPK modeling in cancer nanomedicine [65]. One attempt to overcome the issue was to develop some of the necessary standards and approaches for building models and represent them by expert researchers in the relevant community. In this context, SBML (Systems Biology Markup Language, an XML-based format) [32] as the most widespread language for defining computational biochemical models and SBGN (Systems Biology Graphic Notation) [33] as a standard for graphical representation of molecular networks have been proposed as gold standards for the representation of biological networks and related models [66]. In addition, these efforts focused on the development of tool-independent ways of presenting models to help avoid potential human errors in translation. Note that these initiatives are not widely ingrained in tools for PBPK modeling. However, some of them described below support these critical standards.

MATLAB [67] is regularly used in many scientific fields, including systems biology [68]. In the case of PBPK models and related PBPK analyses, intermediate to advanced programming skills are required, which present a considerable disadvantage for its widespread use in the field [69]. MATLAB provides a range of mathematical and numerical methods for solving PBPK model equations, parameter estimation, and sensitivity analysis [53,54,55,56,61,70,71,72,73,74,75,76]. In particular, several MATLAB toolkits can be used for modeling, simulating, and analyzing PBPK systems, namely SimBiology [57,77,78,79,80], PottersWheel [81,82], and IntiQuan IQM Tools [70]. MATLAB was used to simulate the biodistribution of NPs of different sizes in the range from 46 to 162 nm after intravenous injection into the plasma compartment of rats [54]. The authors showed that the model reproduces the in vivo data from a study of the pharmacokinetics of mesoporous silica NPs [83] and then performed local and global sensitivity analyses to rank the importance of various parameters related to the problem of NP delivery to the tumor. Tumor vascularization (fraction and porosity), tumor blood viscosity, NP size and degradation rate have been shown as the main parameters to consider when calculating NP extravasation in the tumor volume, considering passive targeting as the major targeting mechanism. A similar approach was used to evaluate the biodistribution of dexamethasone in a model of acute lymphoblastic leukemia [53]. The model was based on experimental data obtained after intravenous injection of polymeric NPs. A fluorescent dye was used instead of the drug. The authors noted the difficulties in creating a predictive model, in particular, in the first hours after injections, the model overestimated the blood concentration of the dye. This phenomenon could be explained by particle aggregation and margination [84] in the endothelial wall, and to correct the model, they introduced an additional “other” compartment (i.e., lymphatic system). For longer time points, the model was consistent with the experimental data, demonstrating a rapid accumulation of the particles in the liver and spleen 6 h after their injections. This phenomenon was especially appreciated in consideration that liver and spleen are sites for leukemia blast accumulation and proliferation. MATLAB has also been used to develop a predictive model for tumor targeting of dendrimers functionalized with an insulin-like growth factor 1 peptide analog (NPs for molecular imaging that can bind to mRNA inside cancerous cells) [76]. When fitting the model to the experimental data, the authors found that only 10–20% of the administered NPs were available for transport from the blood to the interstitial tissues and suggested that previous mouse imaging trials used more NPs than necessary. Klapproth et al. [85] used this software to create a model comparing particle biodistribution after intravenous and intratumoral injection. The work focused on superparamagnetic iron oxide NPs for thermal ablation against a potential treatment of brain cancer [86], which is currently being treated with both intravenous and local administration of pharmaceuticals. The model confirmed the experimental data on the dynamic NP redistribution in the organism within 72 h, reaching equilibrium 100 h after local injection. Intravenous injection has demonstrated rapid particle retention in all organs (particularly in the liver and spleen) and subsequent slow release.

Berkeley Madonna [87] is relatively easy for beginners and is sufficient to develop basic PBPK models for routine PBPK analyses (e.g., parameter estimation, route-to-route extrapolation, and Monte Carlo analysis) [69,88]. However, it is not intended for visualization of biochemical systems, requiring specialists in the field of biomedicine and mathematical modeling [89]. Examples of PBPK models of NPs coded in Berkeley Madonna are proposed in [50,62,82,90,91,92,93,94,95,96,97,98]. The model by Cheng et al. was used to analyze 376 experimental data sets on the kinetics of NPs in mice with tumors [50]. The authors confirmed that nanomaterials with a hydrodynamic diameter of less than 10 nm can be delivered to tumors with greater efficiency compared to larger particles. In addition, nanomaterials with a hydrodynamic diameter of over 200 nm have a relatively higher efficiency of delivery to the tumor than particles with a size of about 10−200 nm. The study also showed that rod-shaped NPs show better tumor accumulation than variants with other geometry, including spherical, plate-like, or flake-like shapes. Based on several studies, the authors noted that elongated nanostructures, compared to nanospheres, exhibit greater tumor accumulation and a longer half-life in blood circulation, perhaps because of adherence to endothelial cells lining the blood vessel walls, thus, enhancing the site-specific delivery. Furthermore, they concluded that the administration of nanomaterials with a positive (>10 mV) and almost-neutral (−10 to 10 mV) surface charge provides a similar tumor accumulation, which is higher than for negatively charged NPs. These results were partially confirmed by Zhang et al. [97], who used the Berkeley Madonna software to characterize the biodistribution of spherical and rod-shaped gold nanoprobes and tumor accumulation in a model of lung cancer. They found that while nontargeted rod-shaped NPs showed higher tumor accumulation during the first hours after intravenous injection, similar results were obtained at longer time points between nontargeted rod-shaped and spherical delivery platforms. A clear advantage in tumor accumulation was provided by surface functionalization of the particles with RGD to target αvβ3 integrin-positive cancer cells and tumor angiogenic vessels. In the case of active targeting, the authors found a higher distribution coefficient and a much higher maximum tumor uptake rate constant for rod-shaped particles (compared to spherical particles), which eventually resulted in a more effective inhibition of tumor growth by NP-mediated chemoradiotherapy.

The R language [99] is a powerful high-level programming platform that is used in various fields of study for statistical computing and graphics [100]. The advantages of the R language are that it is freely distributed and can perform all PBPK analyses. However, it requires medium- to high-level programming skills (same as MATLAB) [69]. The R language may be the optimal choice for projects related to Markov chain Monte Carlo analysis or statistical analysis. For example, Cheng et al. used it for normality testing, one-way ANOVA, and simple and multivariate linear regression to analyze data on the efficiency of NP delivery to tumors [50], whereas Chou et al. applied it to optimize the parameters of a multiroute PBPK model constructed for different sizes (1.4–200 nm) of gold NPs in adult rats at different routes of administration (i.e., intravenous, oral gavage, intratracheal instillation, and endotracheal inhalation) [62]. In addition, the R Shiny package can be used to create an interactive web interface for PBPK models [62].

acslX. Many PBPK models for environmental chemicals, drugs, and nanomaterials have been developed using acslX [101], which has been deprecated since November 2015. Lin et al. [69] provided guidance on converting PBPK model code from acslX to alternative modeling tools (Berkeley Madonna, MATLAB, and R) and discussed the advantages and disadvantages of each software package in implementing PBPK models in toxicology. Application of acslX to PBPK modeling of NP delivery can be found in many articles [60,90,91,92,102,103,104,105].

BioUML [106] is an integrated web-based platform for systems biology and data analysis [107], which has been successfully tested for modeling biological systems [108,109]. It allows to translate PBPK models from text formats (e.g., Berkeley Madonna) into SBML and SBGN standards, making these models accessible to a wide range of scientists and providing a quick-start guide to work with them [89]. To represent PBPK models in BioUML, a modular approach is used, according to which the system under study is considered as a set of interconnected subsystems (see the example in the Section 6 below). It is also worth noting that BioUML allows one to conduct an identifiability analysis of fitted parameter values mentioned above as a critical procedure for all PBPK models.

Simcyp Simulator [110] is a software platform for population PBPK modeling and simulation. It links in vitro data to in vivo absorption, distribution, metabolism, excretion (ADME), and pharmacokinetic/pharmacodynamic outcomes to explore clinical scenarios and support drug development decisions, including regulatory submissions and drug labels [111,112,113,114]. This platform contains a library of predefined models and a database of physiological parameters, making it popular with PBPK users [115,116]. However, due to the complex ADME processes of NPs, it may not provide enough flexibility and capability to support complex NP models [26]. Therefore, the use of Simcyp for the study of nanoformulations is rare [117].

GastroPlus [118] is a mechanistically based simulation software package that simulates intravenous, oral, oral cavity, ocular, inhalation, dermal, subcutaneous, and intramuscular absorption, biopharmaceutics, pharmacokinetics, and pharmacodynamics in humans and animals [114,119]. It is widely used for PBPK modeling [116], but, like Simcyp, it is inferior to general purpose programming platforms (Matlab-Simulink, Berkeley Madonna, etc.) in the development of complex PBPK models of NPs [26]. However, GastroPlus has been applied for simulation and in silico prediction of pharmacokinetics and absorption of NPs [120,121,122,123].

PK-Sim [124] is a comprehensive software tool for whole-body PBPK modeling. It enables access to all relevant anatomical and physiological parameters for humans and the most common preclinical animal models that are contained in the integrated database [114,125,126]. PK-Sim offers different customized model structures but allows only minor modifications to them. However, it is fully compatible with the stand-alone general purpose graphical modeling tool MoBi [124], which makes the development of PBPK models (also for NPs) possible [125]. Thus, Aborig et al. used the software to parameterize a PBPK model of the biodistribution of gold NPs obtained by green synthesis in mice [70], while Kullenberg et al. applied it to develop a PBPK model describing the disposition of pegylated liposomal doxorubicin [127]. In addition, the platform includes interfaces to R and MATLAB that are useful to analyze and interpret PK-Sim and MoBi models. Advanced model analysis may, for example, involve statistical analysis of the results obtained or the calculation of local or global sensitivity measures to assess model robustness and quality [125].

In this section, we give a brief description of software that is not widely used for PBPK modeling, but can provide additional tools for mathematical analysis:

Julia [128] is a flexible dynamic programming language appropriate for scientific and numerical computing. It has been shown that in the context of PBPK models, a DifferentialEquations.jl package outperforms MATLAB in solving a stiff system of ordinary differential equations [72].

NONMEM [129] solves pharmaceutical statistical problems in which within-subject and between-subjects variability is taken into account when fitting a PK/PD model to data. Bauer recently created two tutorials to apply this tool to simple and complex scenarios [130,131]. NONMEM allows to build population physiologically based PK models [74].

GNU MCSim [132] is specifically designed to conduct Monte Carlo stochastic simulations that can be used for population variability analysis when the distribution of each parameter of the PBPK model is known [69,133]. In particular, the software was applied to estimate parameters in several PBPK models of NPs [58,105,134].

Phoenix WinNonlin [135] is one of the most common applications used in the industry for the analysis of PK/PD data. The platform is suitable for modeling, simulations, and analysis of PBPK models of NPs as well [71,136].

GraphPad Prism [137] is a commercial scientific software for 2D plotting and statistical analysis [138] which can be applied to evaluate the overall quality of fit between simulated and experimental data [102,103] or to compare delivery efficiency of different types of NPs [50].

Crystal Ball [139] is the spreadsheet-based application for predictive modeling, forecasting, simulation, and optimization which is suitable for statistical analysis (e.g., Monte Carlo simulations [74,98] or regression analysis [95]) to parametrize and validate PBPK models.

TableCurve 2D [140] is an automated curve fitting software for engineers and researchers that includes a wide range of linear and nonlinear models [141]. It can be used to fit the experimental data on the biodistribution of NP in PBPK simulations [98].

ADAPT [142] is a computational modeling platform developed for PK/PD applications. Mager et al. used it to estimate the unknown parameters of the PBPK model for composite nanodevices [82].

NanoSolveIT [143] is a new informatics system for in silico nanosafety assessment [144] that, in particular, can be applied for PBPK modeling of nanomaterial biodistribution [145].

SPSS [146] is a platform for advanced statistical analysis. Zazo et al. used it in the development of a PBPK model of a drug delivery system based on gold NPs for stavudine biodistribution [136].

Microsoft Excel [147] is not a specialized tool for creating PBPK models; however, it can successfully code and analyze them [148,149]. Statistical analysis and calculations with this software can also be used in modeling NP kinetics [105,127].

Statistics Calculator [150] was developed by StatPac, Inc. Despite the modest design, the program works correctly and can be used by researchers for statistical analysis of experimental data when creating PBPK models [151].

Minitab [152] is a statistics package for identifying trends in data, finding and predicting patterns, uncovering hidden relationships between variables, and visualizing data. An example of a nanoparticle study in which the software was used for two- and three-way ANOVA of the antitumor effect and cellular uptake ratio of free and liposomal doxorubicin in human cancer cell lines (HepG2, Huh7, SNU449, and MCF7) is given in [127].

COMSOL Multiphysics [153] is a software package that provides fully coupled multiphysics and single-physics modeling capabilities. This tool is very useful for numerically solving complex physical equations. Thus, Chen et al. applied it to solve the advection-diffusion equations simulating the concentration of superparamagnetic iron oxide NPs coated by gold and conjugated with polyethylene glycol in the cerebral blood and brain tissue in mice [79].

OriginPro [154] is a data analysis and graphing software that includes tools for peak fitting, surface fitting, statistics, and signal processing [155]. It can be utilized for parameter estimation of PBPK models [80].

WebPlotDigitizer [156] is an open source and cross-platform (web and desktop) semiautomated tool for extraction of the numerical data from engineering images of data visualizations. It works with a wide variety of charts (x-y, histograms, polar, ternary, maps, etc.) and can be used to digitize experimental data on the biodistribution of NPs published as graphs [50,60,72,102].

PlotDigitizer [157] is an open source Java program that converts information from 2D plots or graphs to standard x-y values (tabular format). For instance, the program can be used to digitize experimental data from the literature, after which PBPK models can be calibrated using the obtained quantitative values [70].

UN-SCAN-IT [158] is a software that automatically converts graph images to their underlying (x, y) data [159]. It works with most image formats (JPG, TIFF, GIF, BMP, PNG, etc.) and can integrate peak areas, smooth data, take derivatives, rescale graphs, and export (x, y) data for use in other programs. Lee et al. used the software to obtain the PK data for quantum dots from the published figures [160].