Multiphysics pharmacokinetic model for targeted nanoparticles

Nanoparticles (NP) are being increasingly explored as vehicles for targeted drug delivery because they can overcome free therapeutic limitations by drug encapsulation, thereby increasing solubility and transport across cell membranes. However, a translational gap exists from animal to human studies resulting in only several NP having FDA approval. Because of this, researchers have begun to turn toward physiologically based pharmacokinetic (PBPK) models to guide in vivo NP experimentation. However, typical PBPK models use an empirically derived framework that cannot be universally applied to varying NP constructs and experimental settings. The purpose of this study was to develop a physics-based multiscale PBPK compartmental model for determining continuous NP biodistribution. We successfully developed two versions of a physics-based compartmental model, models A and B, and validated the models with experimental data. The more physiologically relevant model (model B) had an output that more closely resembled experimental data as determined by normalized root mean squared deviation (NRMSD) analysis. A branched model was developed to enable the model to account for varying NP sizes. With the help of the branched model, we were able to show that branching in vasculature causes enhanced uptake of NP in the organ tissue. The models were solved using two of the most popular computational platforms, MATLAB and Julia. Our experimentation with the two suggests the highly optimized ODE solver package DifferentialEquations.jl in Julia outperforms MATLAB when solving a stiff system of ordinary differential equations (ODEs). We experimented with solving our PBPK model with a neural network using Julia's Flux.jl package. We were able to demonstrate that a neural network can learn to solve a system of ODEs when the system can be made non-stiff via quasi-steady-state approximation (QSSA). Our model incorporates modules that account for varying NP surface chemistries, multiscale vascular hydrodynamic effects, and effects of the immune system to create a more comprehensive and modular model for predicting NP biodistribution in a variety of NP constructs.

1. Introduction

Clinical medicine has entered an era where nanotechnology is becoming more and more prevalent, with the use of drug-carrying nanoparticles (NP) increasing in recent years. NP can overcome free therapeutic limitations by drug encapsulation, thereby enhancing solubility by promoting transport across cellular membranes (1). Because of this, NPs are increasingly being explored as vehicles for targeted drug delivery to healthy and cancerous tissues, for diagnostic imaging purposes, and to enhance T-cell-based immunotherapies (2, 3). While extensive studies are being performed to determine the efficacy of certain NP based therapeutics in in vitro and in vivo animal models, there exists a translational gap from animal to human studies resulting in only a select few NP being approved for FDA use (4).

One major underlying cause of this translational gap is the difference in the physiology of animal models compared to humans, which have the potential to affect NP behavior and functionality (5). However, the cause of the translational gap that is the motivation for this study is the heterogeneity of NP constructs and experimental models. There are nearly endless NP constructs (e.g., rigid, flexible, spherical, non-spherical, polymeric, DNA-based, etc.), sizes (a few nm to a few microns), and experimental models for translational studies, making NP a complex mix of biology and engineering. The wide array of applications, targets, and physical characteristics of NP significantly impedes the ability of NP to be researched effectively as possible bench-to-bedside therapeutics.

To this end, researchers have begun to turn toward models for adhesion and transport to guide in vivo experimentation and better understand NP targeting behavior and performance in the human body, resulting in more effective and efficient usage for a variety of previously described applications (6–9). However, incorporation of these models in a pharmacokinetic framework has been elusive. Traditional pharmacokinetic (PK) models describe the concentration of drugs in the blood plasma over time using emperical functions (10, 11), while physiology based pharmacokinetics (PBPK) models consist of compartments that represent abstraction of kinetics and transport in and across actual tissues and organ spaces. Existing PBPK models for small molecules or biologics (12–14) use an empirically derived framework for parameterization, resulting in a model that cannot be universally applied with varying NP constructs and experimental settings (15, 16). Additionally, multiphysics aspects, including physiological and hydrodynamic factors governing NP biodistribution and tissue targeting, involve mechanisms operating at multiple lengths and timescales (17). Therefore, a multiscale computationally driven model with physiologically relevant inputs can be utilized to understand the organ-specific biodistribution of NP. The multiscale model must incorporate NP hydrodynamic properties in the vasculature, NP-Endothelial Cell (EC) adhesion properties, and NP subcellular interactions that govern targeted uptake to fully describe the movement and accumulation of NP within the body (6). In order to create a comprehensive multiscale model, NP behavior must be understood at the system, hydrodynamic, and cell adhesion scales.

A previously published multiscale PBPK model has determined binding constants of intracellular adhesion molecule 1 (ICAM1) coated NPs to endothelial cell (EC) surface receptors in mice and humans by utilizing the biophysical properties of the antibody to receptor interactions, and the cell surface (15). Additionally, this model determines the percent injected dose per gram of tissue (%IDG) that distributes to the tissue in given organs at steady-state; non-specific uptake is not accounted for in this model (NP uptake via passive diffusion in the intercellular cleft). The binding constants determined through the previous model (15) model will help to drive the binding characteristics of NP in the multiscale model described in this study. Additionally, the rate at which NP bound to the EC layer are endocytosed into a specific organ tissue must be considered. The concentration of NP retained within the tissue or biodistribution is ultimately most important since this allows for understanding how effectively NP can target tissue.

The purpose of this study is to 1) advance existing steady-state multiscale PBPK models (15) to incorporate NP uptake via nonspecific transport, 2) develop novel multiscale PBPK compartmental models to predict temporal effects, and 3) introduce a compartmental branched vascular model that can predict the effect of hydrodynamic interactions that depend on NP size, flow, and vasculature network properties, 4) perform validation with experimental murine biodistribution data, 5) propose efficient solvers for coupled stiff systems that embody the above properties, as well as make the solvers compatible with contemporary machine-learning-based modules (such as neural networks) which can capture and incorporate multiphysics models in the PBPK workflow.

2. Methods

2.1. Overview

Many Monte Carlo-based models have been proposed in previous works and have been integral in understanding multivalent receptor-NP interactions at a molecular level. Agrawal and Radhakrishnan (18) quantitatively characterized NP-endothelial cell interactions by determining the multivalence of NP binding as well as antigen clustering, ultimately providing future models with details about the energetics of the NP binding process. Ramakrishnan et al. (15) used these binding properties to understand the role that expression levels of NP receptors play when targeting live cells, validating with experimental data. While these adhesion-centric models are necessary for understanding NP characteristics at the molecular level, they do not always include the interaction of NP with the vascular network and translate into the pharmacodynamic scales. By considering hydrodynamic parameters of the vascular system, and cell-binding/uptake parameters into a variety of organ tissues [as determined in previous studies such as Ramakrishnan et al. (15)], a model that is governed by multiple time scales can be created, providing a more physiologically relevant determination of NP biodistribution with temporal resolution. This paper describes three variations of a multiscale pharmacokinetic model: 1) a steady-state model that facilitates the translation of the multivalent free energy of adhesion into a biodistribution; 2) two different compartmental models (models A and B) that are physiologically based and can predict temporal biodistribution; and 3) a compartmental model that incorporates the vascular branching network and can include multiphysics effects such as hydrodynamic interactions, NP size-dependent effects of flow and adhesion, with temporal resolution (6, 19). We describe the steady-state or continuous temporal biodistribution of ICAM-1 targeted NP for a murine model in each case. The multiscale nature of the model described here can allow for customization with various NP sizes, shapes, and uptake parameters, resulting in a customizable predictive platform for different NP chemistries, organisms, and pathophysiologies.

2.2. Steady-state model

A steady-state, multiscale, PBPK model was developed by Ramakrishnan (15) to determine the percent injected dose of NP per gram of tissue (%IDG) in five organ compartments: lung, heart, kidney, liver, and spleen. However, this model only considers NP uptake via antibody-receptor mediated multivalent interactions, failing to account for receptor-mediated internalization and non-specific transport through pores between the intercellular cleft of adjacent endothelial cells, especially in clearance organs. Here, the existing steady-state PBPK model was modified to incorporate this non-specific uptake, to have the modified model better predict in vivo biodistribution data at short times (<30 min, when NP does not internalize substantially) than the original model as determined by the R (square root of the coefficient of determination R²) values. Several versions of the model (incorporating different multiphysics) are considered based on specific adhesion of NP to epithelial cell membrane surfaces, including flat membrane and membrane mimicking live cells. In addition, the presence of resident macrophages and activated macrophages are considered.

The original model (15) is described using Equation (1),

$\begin{array}{l}\%idg=\left\{{\kappa }_p{K}_{EC}{C}_{out}+\frac{{\phi }_{EC}{K}_{EC}{L}_{EC,b}}{D_{EC}}{C}_{out}\right\}\times \frac{L_{cap}}{L_{EC,b}}\\ +\frac{{\phi }_{EC}{K}_M{L}_{EC,b}}{D_M}{C}_{out}\times \frac{L_{cap}}{L_{M,b}},& (1)\end{array}$

where κ_p is the non-specific binding of NP, K_EC is the association constant for binding of NP to endothelial cell surface receptors, K_M is the association constant for NP binding to a macrophage cell, D_EC and D_M represent the diameter of endothelial cells and macrophage cells, respectively, φ_M and φ_EC represent the concentration of endothelial cells and macrophage cells in the target tissue, respectively, L_EC,b and L_M,b represents the distance from an EC or macrophage surface receptor that the NP can successfully bind, respectively, C_out represents the injected concentration of NP, and L_cap represents the size of the cell-free layer in the capillary in which the NP is perfused. Incorporating the non-specific transport of NP into the tissue, the modified model equation can be described using Equation (2):

$\begin{array}{l}\%idg=\left\{{\kappa }_p+\frac{{\phi }_{EC}{K}_{EC}{L}_{EC,b}}{D_{EC}}{C}_{out}\right\}\times \frac{L_{cap}}{L_{EC,b}}\\ +\frac{{\phi }_{EC}{K}_M{L}_{EC,b}}{D_M}{C}_{out}\times \frac{L_{cap}}{L_{M,b}}.& (2)\end{array}$

The model above is expected to approximate the biodistribution at short timescales, defined as the regime in which the observation time is less than the timescale for internalization. The model is easy to compute as it does not involve solving dynamics equations owing to its steady-state nature. Later, we show that the temporal model predicts the same behavior as the steady-state model at these short timescales, thereby validating the steady-state approximation utilized here for short time scales.

2.3. Temporal model

2.3.1. Compartmental model development

The goal of the compartmental model is to develop a basic framework for determining targeted NP biodistribution in a murine model that can act as a predictive model when provided with experimentally and empirically derived parameters. This model consists of five to seven organ compartments (lung, heart, kidney, liver, spleen and in model A; lung, heart, kidney, liver, spleen, gut, and 'other in model B) each of which is interconnected via the arterial and venous compartments. The ‘other' compartment consists of all organs that are not explicitly included in the model. We have developed two ways in which the organ compartments can be connected by the arteries and veins, shown in Figure 1. Figure 1A, is the model A configuration (an oversimplified version of the circulatory system with 5 organ compartments), while Figure 1B shows the model B configuration model, which is more physiologically relevant (lung circulation has been separated out to keep track of oxygenation and oxygen distribution, and gut and spleen compartments are coupled to the liver compartment, and gut/other compartments were incorporated as well).

FIGURE 1

Figure 1. Compartmental model configurations. (A) model A configuration. (B) model B configuration.

Each organ compartment in either model configuration consists of three additional compartments: vascular, endothelial cell, and tissue compartments shown in Figure 2. As NPs enter an organ compartment through the vascular compartment, they can either bind to the ICAM-1 endothelial cell surface receptors (K_on), enter the organ tissue compartment via non-specific uptake (K_NS), or be degraded (K_deg). Then, once bound to the endothelial cell layer, the NP can either unbind from the EC compartment returning to the vascular compartment (K_off), be taken into the organ tissue via transcytosis (K_up), or degraded (K_deg). Once NPs are in the tissue compartment, they can be degraded (K_deg). It is also important to note that NP can be degraded (K_deg) within the arterial and venous compartments as well. The difference in time scales represented in the arteries/veins and the cellular scale compartments contribute to the temporal multiscale nature of the model.

FIGURE 2

Figure 2. Sub-compartments within each organ compartment. Nanoparticles can be taken directly into the tissue via nonspecific uptake (K_NS), bound to endothelial cell surface receptors (K_on), unbound from endothelial cell surface receptors (K_off), or taken into the tissue via endo or transcytosis (K_up).

2.3.2. Model A equations

The model A configuration can be described using a system of 17 linear ordinary differential equations (ODEs). Model A is described using Equations (3–7), where i = lung, heart, kidneys, liver, spleen. The venous compartment is described by:

$\begin{array}{ll}{V}_{vein}\frac{d{C}_{vein}}{dt}=C_{vein}{\displaystyle \sum }{Q}_i-Q_{vein}{C}_{vein}-K_{deg}^{vein}{C}_{vein}{V}_{vein}.& (3)\end{array}$

Equation (3) describes the change in concentration of NP in the venous compartment over time, where Q_i and Q_vein represent the flow of blood through organ compartments and the veins, respectively, V_vein denotes the volume of blood in the vein, C_vein is the concentration of NP in the veins, and $K_{deg}^{vein}$ is the degradation rate of NP in the veins. The arterial compartment is described by:

$\begin{array}{ll}{V}_{art}\frac{d{C}_{art}}{dt}=Q_{art}{C}_{art}-C_{art}{\displaystyle \sum }{Q}_i-K_{deg}^{art}{C}_{art}{V}_{art}.& (4)\end{array}$

Equation (4) describes the change in concentration of NP in the arterial compartment over time, where Q_art represents the flow of blood through the arteries, V_art denotes the volume of blood in the arteries, C_art is the concentration of NP in the arteries, and $K_{deg}^{art}$ is the degradation rate of NP in the arteries. For each tissue type i, the vascular compartment is described by:

$\begin{array}{l}{V}_{bl}\frac{d{C}_{bl}^i}{dt}=Q_{art}{C}_{art}-Q_{vein}{C}_{vein}-K_{on}^i{V}_{bl}^i{C}_{bl}^i+K_{off}^i{C}_{EC}^i{V}_{EC}^i\\ -K_{NS}^i{C}_{bl}^i{V}_{bl}^i.& (5)\end{array}$

Equation (5) describes the change in concentration of NP in the vascular compartment within each organ compartment of the model over time, where $V_{bl}^i$ and $V_{EC}^i$ represent the volume of blood in the vascular and the endothelial cell compartments (which is defined by the product of the length of the endothelial cell receptors and the surface area of the vascular compartment) of the organ, respectively, $C_{bl}^i$ and $C_{EC}^i$ denotes the concentration of NP in the vascular and endothelial cell compartments of the organ, respectively, $K_{on}^i$ denotes the rate of binding of NP to the endothelial cell surface receptors, $K_{off}^i$ denotes the rate of NP unbinding from the endothelial cell surface receptors, and $K_{NS}^i$ denotes the rate of nonspecific NP uptake into the organ tissue. The endothelial compartment in each tissue is described by:

$\begin{array}{l}{V}_{EC}\frac{d{C}_{EC}^i}{dt}=K_{on}^i{C}_{bl}^i{V}_{bl}^i-K_{up}^i{C}_{EC}^i{V}_{EC}^i\\ -K_{off}^i{C}_{EC}^i{V}_{EC}^i-K_{deg}^i{C}_{EC}^i{V}_{EC}^i.& (6)\end{array}$

Equation (6) describes the change in concentration of NP bound to the endothelial cell surface receptors over time, where $K_{up}^i$ denotes the rate of uptake of NP into the organ tissue via transcytosis. Each tissue compartment is described by:

$\begin{array}{ll}{V}_T\frac{d{C}_T^i}{dt}=K_{NS}^i{C}_{bl}^i{V}_{bl}^i+K_{up}^i{C}_{EC}^i{V}_{EC}^i-K_{deg}^i{C}_T^i{V}_T^i.& (7)\end{array}$

Equation (7) describes the change in concentration of NP in the organ tissue over time where $C_T^i$ and $V_T^i$ denote the concentration of NP in the tissue compartment of the organ and the volume of the tissue compartment, respectively.

2.3.3. Model B equations

The model B configuration (Figure 1B) can be described using a system of 23 linear ODEs, however, due to the addition of two organs and the alternative configuration of the model, several equations differ, while the overall structure is the same. Equations (8–13) Describe the modified model configuration:

$\begin{array}{l}{V}_{vein}\frac{d{C}_{vein}}{dt}=Q_{kidney}{C}_{bl}^{kidney}+Q_{hep}{C}_{bl}^{liver}\\ +Q_{other}{C}_{bl}^{other}-Q_{vein}{C}_{vein}-V_{Vein}{K}_{deg}^{vein}.& (8)\end{array}$

Equation (8) describes the change in concentration of NP in the veins over time, where Q_kidney and Q_other denote the flow rate of blood through the kidneys and “other” compartment, respectively, Q_hep denotes the combined flow rate of blood through the liver, spleen, and gut compartments (Q_hep = Q_liver + Q_spleen + Q_gut). C_bl denotes the concentration of NP in the vascular compartment of each respective organ (lung, heart, kidneys, liver, spleen, gut, other)

$\begin{array}{l}{V}_{art}\frac{d{C}_{art}}{dt}=Q_{art}{C}_{bl}^{heart}-C_{art}(Q_{lung}+Q_{heart}\\ +Q_{gut}+Q_{other})-C_{art}{V}_{art}{K}_{deg}^{art}.& (9)\end{array}$

Equation (9) describes the change in concentration of NP in the arteries over time where Q_lung and Q_heart denote the flow rate of blood through the lung and heart compartments, and $C_{bl}^{heart}$ is the concentration of NP in the vascular compartment of the heart.

$\begin{array}{l}{V}_{bl}^{lung}\frac{d{C}_{bl}^{lung}}{dt}=Q_{heart}{C}_{bl}^{heart}-Q_{lung}{C}_{bl}^{lung}-K_{on}^{lung}{V}_{bl}^{lung}{C}_{bl}^{lung}\\ +K_{off}^{lung}{V}_{EC}^{lung}{C}_{EC}^{lung}\\ -K_{NS}{V}_{bl}^{lung}{C}_{bl}^{lung}-K_{deg}^{lung}{V}_{bl}^{lung}{C}_{bl}^{lung}.& (10)\end{array}$

Equation (10) describes the change in concentration of NP in the vascular compartment of the lung over time.

$\begin{array}{l}{V}_{bl}^{heart}\frac{d{C}_{bl}^{heart}}{dt}=(Q_{vein}{C}_{vein}+Q_{lung}{C}_{bl}^{lung})\\ -(Q_{lung}{C}_{bl}^{heart}+Q_{vein}{C}_{bl}^{heart})\\ -C_{bl}^{heart}{V}_{bl}^{heart}+K_{off}^{heart}{C}_{EC}^{heart}{V}_{EC}^{heart}\\ -K_{NS}^{heart}{C}_{bl}^{heart}{V}_{bl}^{heart}{-K}_{deg}^{heart}{C}_{bl}^{heart}{V}_{bl}^{heart}.& (11)\end{array}$

Equation (11) describes the change in concentration of NP in the vascular compartment of the heart over time.

$\begin{array}{l}{V}_{bl}^i\frac{d{C}_{bl}^i}{dt}=Q_i{C}_{art}-Q_i{C}_{bl}^i-K_{on}^i{C}_{bl}^i{V}_{bl}^i+K_{off}^i{C}_{EC}^i{V}_{EC}^i\\ -K_{NS}^i{C}_{bl}^i{V}_{bl}^i-K_{deg}^i{C}_{bl}^i{V}_{bl}^i.& (12)\end{array}$

Equation (12) describes the change in concentration of the NP in the vascular compartment of the kidneys, spleen, gut, and ‘other' compartment over time (where i = kidney, spleen, gut, other).

$\begin{array}{l}{V}_{bl}^{liver}\frac{d{C}_{bl}^{liver}}{dt}=Q_{spleen}{C}_{bl}^{spleen}+Q_{gut}{C}_{bl}^{gut}+Q_{liver}{C}_{art}\\ -Q_{hep}{C}_{bl,}^{liver}-K_{on}^{liver}{C}_{bl}^{liver}{V}_{bl}^{liver}\\ +K_{off}^{liver}{C}_{EC}^{liver}{V}_{EC}^{liver}\\ -K_{NS}^{liver}{C}_{bl}^{liver}{V}_{bl}^{liver}-K_{deg}^{liver}{C}_{bl}^{liver}{V}_{EC}^{liver}.& (13)\end{array}$

Equation (13) describes the change in concentration of the NP in the vascular compartment of the liver over time. Equations (6) and (7) were used in the original model to describe the change in concentration of NP in the endothelial cell compartments and tissue compartments, respectively, can also be used in the modified model configuration.

Additionally, Equation (14) was used to express the total NP degraded in the system at a given time t, which will be useful to determine mass conservation of the system.

$\begin{array}{l}Mo{l}_{deg}^t={\displaystyle \sum _{g=1}^g(C_{bl}^g{V}_{bl}^g{K}_{bl,deg}^g\Delta t+C_{EC}^g{V}_{EC}^g{K}_{EC,deg}^g\Delta t}\\ +C_T^g{V}_T^g{K}_{T,deg}^g\Delta t)\\ +C_{vein}^t{V}_{vein}^t{K}_{vein,deg}\Delta t+C_{art}^t{V}_{art}^t{K}_{art,deg}\Delta t)& (14)\end{array}$

where t is the total time the model is run, g is the total number of organs in the model, and Δt is the t step of the model.

Equation (15) was used to determine the mass conservation of the system. If the system is closed, the line formed by ${Mol}_{total}^t$ over every time point t of the simulation should have a slope of 0.

$\begin{array}{l}Mo{l}_{total}^t={\displaystyle \sum _{g=1}^g}(C_{bl}^g{V}_{bl}^g+C_{EC}^g{V}_{EC}^g+C_T^g{V}_T^g)+C_{vein}^t{V}_{vein}\\ +C_{ar{t}_{art}^{tV}}+{Mol}_{deg}^t.& (15)\end{array}$

Equation (34) describes the mass conservation in molar basis, this equation was used to determine the mass conservation when the system of ODE was set up in terms of molar profiles.

$\begin{array}{ll}\begin{array}{c}\text{Total Mass (t)}={\displaystyle \sum _i}{N}_i(t),\hfill \\ i\in \left\{Vasculature,Endothelial,Tissue,Vein,Artery,Degraded\right\}.\end{array}& (16)\end{array}$

2.4. Vasculature branching

2.4.1. Branched model development

The purpose of developing a branched model was to create a more detailed and physiologically relevant version of the basic compartmental model. Additionally, utilizing this branched model will ultimately increase the specificity of uptake rate constants for NP of various sizes. Furthermore, the branched model will enable the inclusion of key margination and hydrodynamic interactions whose effects are determined by the flow rate and blood hematocrit concentration. The branched network consists of a branched vascular tree that begins at the main arteries and veins and bifurcates into the capillary beds, connecting the arterial and venous branching networks (Figure 3). While asymmetric branching patterns characterize typical vasculature networks, for simplicity, the branching model described here will consist of identical daughter vessel segments at each generation of bifurcation. Daughter vessels will continue to bifurcate, their diameters following the power-law relationship until the diameter of the vessel approaches the size of a red blood cell (the smallest vessel in our network). Below the development of the branched network is discussed in greater detail.

FIGURE 3

Figure 3. Branched Network. An example branching network with i = 3 bifurcations and n = 6 generations.

The development of this branched network was modeled after the branching network described in Yang and Wang (20), but is a simplified version of their three-dimensional vascular branching network model. The diameter of the daughter vessels is governed by the power law relationship, where the parent vessel is d₀ and daughter vessels are d₁ and d₂, and k=3, which is typically assumed based on a minimum dissipation principle representing a stable flow condition:

$\begin{array}{ll}{d}_0=d_1^k+d_2^k,& (17)\end{array}$

where d₀ > d₁ = d₂

Due to the branching nature of the model, we can determine the number of segments (N) at any generation of branching (n) using the relationship N = 2ⁿ. While the number of segments increases as the number of generations increase, the diameter of the branches decreases. Suppose that the diameter of the parent vessel at a given bifurcation (i) is d_{i, 0}, and the diameter of the daughter vessels are d_{i, 1} and d_{i, 2}. Given a symmetric bifurcation (d_{i, 1} = d_{i, 2}), the following relationship can be used to relate the diameter of the parent vessels to that of the daughter vessels.

$\begin{array}{ll}{d}_{i,2}=d_{i,1}=\sqrt[k]{0.5•d_{i,0}^k}.& (18)\end{array}$

Then, the length of each vessel can be determined using the know length to diameter ratio β = 3, and the diameter of the vessel, d.

$\begin{array}{ll}L(d)=\beta •d.& (19)\end{array}$

This network will begin at the diameter of the main artery d₀ = 600 μ m (main artery diameter in mouse), continue to bifurcate until the diameter of the daughter vessels reaches the diameter of a red blood cell (d = 15 μ m). This lower-bound in vessel diameter results in a network of 16 bifurcations and 32 generations of branching for each branching element. A schematic of this branching network can be seen in Figure 3.

After constructing this branched network, it is important to determine the surface area and volume of the vascular network element as a whole. Supposing d_i, l_i, and N_i are the diameter, length, and the number of daughter vessels in a given generation i, respectively, Equation (19) can be used to describe the total vessel surface area in a branching element.

$\begin{array}{ll}SA=2{\displaystyle \sum _{i=0}^i}\pi d_i{l}_i{N}_i.& (20)\end{array}$

Equation (20) can be used to determine the total vascular volume of a branching element,

$\begin{array}{ll}V=2{\displaystyle \sum _{i=0}^i}\pi {\left(\frac{d_i}{2}\right)}^2{l}_i{N}_i.& (21)\end{array}$

It is important to note that the number of generations, volume, and surface area of one branching element is fixed. I.e., the volume of a single branching element does not differ across organs. Each organ will have a different number of branching elements dependent on ϕ_i, which is the ratio of total organ vascular volume to that of one branching element.

2.4.2. Branched model equations

The system of ODEs used to describe the branched compartmental model consists of 457 equations. The organs are organized in the model B configuration, and therefore the equations in the branched model take on a similar form when describing the transport between organs. Therefore, the equations describing the concentration of NP in the arteries and veins over time are the same as in the modified compartmental model configuration, Equations (8) and (9), for veins and arteries, respectively.

Equations (21–32) describe the rest of the branched vascular compartmental model that bifurcates from the diameter of the main vein to the size of a red blood cell and branches back out to the diameter of the main artery. Each organ compartment contains four types of equations to describe the concentration of NP in the vasculature: 1) one equation describing NP concentration in the first generation (n = 1), 2) series of equations describing NP concentration in the branching network until the diameter of the branch is that of a red blood cell (n = 2 to n = 16), 3) series of equations describing NP concentration in branching network until diameter of branch branches back out to the generation before the diameter of the vasculature is that of the main artery (n = 17 to n =31), 4) one equation describing the NP concentration in the last generation (n = 32).

$\begin{array}{l}{V}_{lung}^{bl}\frac{d{C}_{lung}^{bl,1}}{dt}=Q_{heart}{C}_{heart}^{bl,32}-Q_{lung}^1{C}_{lung}^{bl,1}{N}^1\\ -{\phi }_{lung}{C}_{lung}^{bl}{V}_{lung}^{bl,1}(K_{lung}^{on}+K_{lung}^{NS})\\ -{\phi }_{lung}{C}_{lung}^{bl}{V}_{lung}^{bl,1}{K}_{deg}^{lung}\\ +K_{lung}^{off}{\phi }_{lung}{C}_{lung}^{EC}{V}_{lung}^{EC,1}.& (22)\end{array}$

Equation (22) describes the NP concentration in generation n = 1 in the vascular branching network of the lung. $C_{lung}^{bl,1}$ denotes the NP concentration of the lung vasculature in generation n = 1 of the branching network,$Q_{lung}^1$ represents the flow rate through generation n = 1 of the lung vascular network, $C_{heart}^{bl,32}$ denotes the NP concentration of the heart vasculature in generation n = 32 (the concentration of NP exiting the heart compartment), N¹ represents the number of branches in generation n = 1, φ_lung represents the total number of branching elements in the given organ compartment, $V_{lung}^{bl,1}$ denotes the vascular volume of a single branching element in generation n = 1, and $V_{lung}^{EC,1}$ denotes the volume of the lung endothelial cells of a single branching element in generation n = 1.

$\begin{array}{l}{V}_{lung}^{bl,n}\frac{d{C}_{lung}^{bl,n}}{dt}=Q_{lung}^{n-1}{C}_{lung}^{bl,n-1}{N}^{n-1}-Q_{lung}^n{C}_{lung}^{bl,n}{N}^n\\ -{\phi }_{lung}{C}_{lung}^{bl,n}{V}_{lung}^{bl,n}(K_{lung}^{on}+K_{lung}^{NS}+K_{lung}^{deg})\\ +K_{lung}^{off}{\phi }_{lung}{C}_{lung}^{EC,n}{V}_{lung}^{EC,n}.& (23)\end{array}$

Equation (23) describes the NP concentration in the lung vasculature in generation n = 2 to n = 16 of the branching network. $Q_{lung}^{n-1}$, $C_{lung}^{bl,n-1}$, and Nⁿ⁻¹ denote the flow rate, NP concentration, and number of branches in the lung vasculature of the previous generation of the branching network.

$\begin{array}{l}{V}_{lung}^{bl,n}\frac{d{C}_{lung}^{bl,n}}{dt}=Q_{lung}^{n-1}{C}_{lung}^{bl,n-1}{N}^n-Q_{lung}^n{C}_{lung}^{bl,n-1}{N}^{n+1}\\ -{\phi }_{lung}{C}_{lung}^{bl,n}{V}_{lung}^{bl,n}(K_{lung}^{on}\\ +K_{lung}^{NS}+K_{lung}^{deg})\\ +K_{lung}^{off}{\phi }_{lung}{C}_{lung}^{EC,n}{V}_{lung}^{EC,n}.& (24)\end{array}$

Equation (24) describes the NP concentration in the lung vasculature in generation n = 17 to n = 31 of the branching network.

$\begin{array}{l}{V}_{lung}^{bl,32}\frac{d{C}_{lung}^{bl,32}}{dt}=Q_{lung}^{31}{C}_{lung}^{bl,31}{N}^{32}-Q_{lung}^{32}{C}_{lung}^{bl,32}\\ -{\phi }_{lung}{C}_{lung}^{bl,32}{V}_{lung}^{bl,32}(K_{lung}^{on}+K_{lung}^{NS}+K_{lung}^{deg})\\ +K_{lung}^{off}{\phi }_{lung}{C}_{lung}^{EC,32}{V}_{lung}^{EC,32}& (25)\end{array}$

Equation (25) describes the NP concentration in the lung vasculature in generation n = 32 of the branching network. The heart compartment vascular branching network can be described similarly to the lung compartment, with a series of four equations.

$\begin{array}{l}{V}_{heart}^{bl,1}\frac{d{C}_{heart}^{bl}}{dt}=Q_{vein}{C}_{vein}+Q_{lung}{C}_{lung}^{bl,32}-Q_{heart}^1{C}_{heart}^{bl,1}{N}^1\\ -{{\phi }_{heart}C}_{heart}^{bl,1}{V}_{heart}^{bl,1}(K_{heart}^{on}+K_{heart}^{NS}+K_{heart}^{deg})\\ +K_{heart}^{off}{\varphi }_{heart}{C}_{heart}^{EC,1}{V}_{heart}^{EC,1}.& (26)\end{array}$

Equation (26) describes the NP concentration in the heart vasculature in generation n = 1 of the branching network. The equations describing the NP concentration in the heart vasculature in generations n = 2 to n = 16 and n = 17 to n = 31 are the same as equations (22) and (23), respectively, but i = lung should be replaced with i = heart.

$\begin{array}{l}{V}_{heart}^{bl,32}\frac{d{C}_{heart}^{bl,32}}{dt}=Q_{heart}^{31}{C}_{heart}^{bl,31}{N}^{32}-Q_{heart}{C}_{heart}^{bl,32}+Q_{vein}{C}_{vein}\\ {-\phi }_{heart}{C}_{heart}^{bl,32}{V}_{heart}^{bl,32}(K_{heart}^{on}+K_{heart}^{NS}+K_{heart}^{deg})\\ K_{heart}^{off}{\varphi }_{heart}{C}_{heart}^{EC,32}{V}_{heart}^{EC,32}.& (27)\end{array}$

Equation (27) describes the NP concentration in the heart vasculature in generation n = 32 of the branching network. The equations describing the concentration of NP in the liver compartment are again constructed similarly to the branching equations in the lung and heart compartments and again consist of a series of four equations.

$\begin{array}{l}{V}_{liver}^{bl,1}\frac{d{C}_{liver}^{bl,1}}{dt}=Q_{spleen}{C}_{spleen}^{bl,32}+Q_{gut}{C}_{spleen}^{bl,32}+Q_{other}{C}_{other}^{bl,32}\\ -Q_{liver}{C}_{liver}^{bl,32}{N}^1-{\phi }_{liver}{C}_{liver}^{bl,1}{V}_{liver}^{bl,1}\\ (K_{liver}^{on}+K_{liver}^{NS}+K_{liver}^{deg})\\ +K_{liver}^{off}{\varphi }_{liver}{C}_{liver}^{EC,1}{V}_{liver}^{EC,1}.& (28)\end{array}$

Equation (28) describes the NP concentration in the liver vasculature in generation n = 1 of the branching network. The equations describing the NP concentration in the liver vasculature in generations n = 2 to n = 16 and n = 17 to n = 31 are the same as Equations (22) and (23) respectively, but i = lung should be replaced with i = liver.

$\begin{array}{l}{V}_{liver}^{bl,32}\frac{d{C}_{liver}^{bl,32}}{dt}=Q_{liver}^{31}{C}_{liver}^{bl,31}{N}^{32}-Q_{liver}^{32}{C}_{liver}^{bl,32}\\ -{\varphi }_{liver}{C}_{liver}^{bl,32}{V}_{liver}^{bl,32}(K_{liver}^{on}+K_{liver}^{NS}+K_{liver}^{deg})\\ +K_{liver}^{off}{\phi }_{liver}{C}_{liver}^{EC,32}{V}_{liver}^{bl,32}.& (29)\end{array}$

Equation (29) describes the NP concentration in the liver vasculature in generation n = 32 of the branching network. The equation describing the NP concentration in the kidneys, spleen, gut, and other compartments are again constructed similarly to the branching equations in the lung and heart compartments, and again consist of a series of four equations.

$\begin{array}{l}{V}_i^{bl,1}\frac{d{C}_i^{bl,1}}{dt}=Q_{vein}{C}_{vein}-Q_i^{bl,1}{C}_i^{bl,1}{N}^1\\ -{{\phi }_{i}C}_i^{bl,1}{V}_i^{bl,1}\left(K_i^{on}+K_i^{NS}+K_i^{deg}\right)\\ +K_i^{off}{\phi }_i{C}_i^{EC,1}{V}_i^{EC,1}.& (30)\end{array}$

Equation (30) describes the NP concentration in the i = kidneys, spleen, gut, and other compartment's vasculature in generation n = 1 of the branching network. The equations describing the NP concentration in the liver vasculature in generations n = 2 to n = 16 and n = 17 to n = 31 are the same as equations (22) and (23), respectively, but i = lung should be replaced with i = kidneys, spleen, gut, other.

$\begin{array}{l}{V}_i^{bl,32}\frac{d{C}_i^{bl,32}}{dt}=Q_i^{bl,31}{C}_i^{bl,31}{N}^{32}-Q_i{C}_i^{bl,32}\\ -{{\phi }_{i}C}_i^{bl,32}{V}_i^{bl,32}(K_i^{on}+K_i^{NS}+K_i^{deg})\\ +K_i^{off}{\phi }_i{C}_i^{EC,32}{V}_i^{EC,32}.& (31)\end{array}$

Equation (31) describes the NP concentration in the liver vasculature in generation n = 32 of the branching network. Next, we can describe the NP concentration that is bound to the endothelial cell layer.

$\begin{array}{l}{V}_i^{EC,n}\frac{d{C}_i^{EC,n}}{dt}=K_i^{on}{{\phi }_{i}C}_i^{bl,n}{C}_i^{bl,n}\\ -{{\phi }_{i}C}_i^{EC,n}{V}_i^{EC,n}(K_i^{up}+K_i^{off}+K_i^{deg}).& (32)\end{array}$

Equation (32) describes the NP concentration bound to the endothelial cells of i = lung, heart, kidney, liver, spleen, gut and other in generation n =1 to n = 32. Finally we can describe the concentration of NP that is distributed to the organ tissue.

$\begin{array}{l}{V}_i^{T,n}\frac{d{C}_i^{T,n}}{dt}={\displaystyle \sum _{n=1}^{32}}\left(K_i^{up}{\phi }_i{C}_i^{EC,n}{V}_i^{EC,n}+K_i^{NS}{{\phi }_{i}C}_i^{bl,n}{V}_i^{bl,n}\right)\\ -K_i^{deg}{\phi }_i{C}_i^{T,n}{V}_i^{T,n}.& (33)\end{array}$

Equation (33) describes the NP concentration in the tissue of i = lung, heart, kidneys, liver, spleen, gut, and other in generation n = 1 to n = 32, where the sum of the antibody-receptor mediated endocytosis ($K_i^{up}$) and non-specific uptake ($K_i^{NS}$) is summed over all generations of branching.

2.5. Parameterization

2.5.1. Parameters for compartmental model

The compartmental model is parameterized with a variety of physiological inputs such as blood, tissue, endothelial cell volumes, and blood flow rates collected from previously published sources (21, 22). Other parameters such as K_on and K_off for ICAM antibody-coated NP were computationally determined in a previous study (15), while yet other rate constants were determined via local sensitivity analysis through comparison to existing translational studies (23).

The flow rates through each individual organ, Q_i, were determined using Supplementary Table 3 from Diehl and Morse (21), which gave blood flow rates through mouse organ vasculature in units of L/g/min, which were ultimately converted into units of L/min for use in the compartmental models by using mouse organ weight data from Boswell et al. (22) and the female and male masses were averaged. Additionally, to calculate the total flow, Q, the sum of flow rates across all organs was taken.

Supplementary Table 1 in Diehl and Morse et al. (21) listed values for vascular volume in various mouse tissues. The volumes were listed in units of L/g. The mouse organ weight data from Boswell et al. (22) was used to determine the tissue volume for the entire organ, $V_{bl}^i$.

The volume of blood in the mouse vein was computed using the following equation, $V_{total}=V_{art}+V_{vein}+\sum V_{bl}^i$,where V_total= 2146 μL. The volume of blood in the veins (V_vein) and the arteries (V_art) is considered to be the same. Then the volume of blood in the mouse venous system was computed.

Values for interstitial tissue volume and extracellular tissue volume were obtained using Supplementary Table 2 from Diehl and Morse (21). Values for interstitial tissue volume (given in L/g) were used to determine the volume of each organ in its entirety, $V_T^i$. By using the organ weights from Diehl and Morse (21), as was done to determine flow rates (Q_i), the tissue volume for each organ in the mouse was computed with units of L.

The parameter l_NC,b denotes the height of the endothelial cell layer, is constant between organs and species and has already been defined by a previous model (15). A_i can be calculated by using known relationships from Ramakrishnan et al. (15), e.g., It is known that ${\varphi }_{EC}=\frac{V_{EC}}{V_T}$ and ${\varphi }_{EC}=\frac{l_{EC}^2{D}_{EC}}{V_T}$. Given ϕ= 0.3, $D_{EC}=5\times {10}^{-6}$, and the V_T for each mouse organ previously described, $l_{EC}^2$ can be calculated. When A_i and l_NC,b are multiplied together, the resulting value will be representative of the volume of the endothelial cell layer; $V_{EC}^i$.

The binding rate of NP to the endothelial cell surface ($K_{on}^i$) was determined using the relationship $\frac{D}{l^2}$ given in (15). The unbinding rate of NP from the endothelial cell surface ($K_{off}^i$) can be determined by using $K_{on}^i$ and K_EC from Ramakrishnan et al. (15) for each organ). However, the log of the K_EC had to be taken to accommodate for crowding effects (24). So, $K_{off}^i$ is calculated as a ratio of $K_{on}^i$ to $K_{EC}^i$; i.e., we assume a diffusion limited on-rate, where the diffusion occurs through the glycocalyx. The off-rate is computed based on the on-rate and the equilibrium constant.

2.5.2. Parameterization of NP-size dependent branched model

The relation between K_EC and particle diameter was obtained from Figure 9 of McKenzie et al. (24). A linear relationship between log(K_eq)) and particle diameter (a) was established using the slope-point form for every organ, the point on the line being the (a, log(K_EC)), for a = 800 nm, which was obtained using Monte Carlo simulation in Ramakrishnan et al. (15). The binding rate is computed as; $K_{on}=\frac{D}{l^2}$ where l = thickness of the glycocalyx layer and D is given by the Stokes-Einstein equation; $D=\frac{K_{B}T}{6\pi \eta r}$, then $K_{off}=\frac{K_{on}}{log(K_{EC})}$. The values of K_ON, K_OFF, and log(K_EC), are reported in Supplementary Tables S1–S3, in supporting information. K_UP K_NS

2.5.3. Local sensitivity analysis

Parameters that do not have an explicit value stated in previously published literature were subject to local sensitivity analysis. These parameters included $K_{deg}^i,K_{NS}^i,K_{up}^i$; where i = lung, heart, kidneys, liver, spleen, gut, and other. To perform the local sensitivity analysis the values of $K_{deg}^i,K_{NS}^i,K_{up}^i$ were changed incrementally so the output curve closely resembled that of an experimental data set (discussed in Section 2.5.1). Larger $K_{NS}^i$and $K_{up}^i$ values resulted in a steeper initial biodistribution curve. Generally, $K_{NS}^i$and $K_{up}^i$ had a similar magnitude result in the biodistribution curve. Larger $K_{deg}^i$ resulted in a smaller maximum biodistribution and a quicker decay in the biodistribution curve. The results of an example local sensitivity analysis for the spleen is shown in Figure 4. The parameters determined via the local sensitivity analysis that were used in the compartmental models can be found at the SI GitHub link.

FIGURE 4

Figure 4. Local sensitivity analysis: Plot for the K_deg, K_up, and K_NS parameters in the spleen compartment of the model. The blue line is the model output using optimized parameters obtained via the sensitivity analysis. Other lines show the model output with unoptimized K_deg, K_up, and K_NS parameters.

Following the sensitivity analysis in model A, model B, and the branching model, values for $K_i^{deg},K_i^{NS},$ and $K_i^{up}$ were determined and are available at the SI GitHub link. It is logical to assume that non-specific uptake should be greatest in the liver and gut because they are considered clearance organs. However, the local sensitivity analysis revealed that the highest non-specific uptake rate was in the lungs. This discrepancy is due to the presence of potentially invalid experimental data. It was apparent that no matter how high the liver K^NS value was set during the sensitivity analysis, the biodistribution produced by the model was always significantly smaller than the experimental data set. In fact, forcing a match to the liver data can only be realized at the expense of violating the conservation of mass, which indicates an experimental error in the reporting of the liver data. It is important to note that earlier iterations of the model not described in the paper used other experimental data sets for validation, and the liver biodistribution output matched the experimental data well. So, we conclude that the experimental data reported by Dong et al. (23) for the liver is invalid. As for the gut and other K^NS values, they were entirely arbitrary. Since the experimental data set used for validation did not report NP biodistribution data for gut or other, we could not perform a valid sensitivity analysis on this parameter for these compartments.

2.5.4. Global sensitivity analysis

We performed global sensitivity analysis to determine the combined effect of model parameters on specific quantities of interest. In particular, we determine the most significant parameters for maximum uptake of NPs inside organ tissue and the mean value of NPs over time in the endothelial compartment. Maximum uptake of NPs inside organ tissue is an essential target for better design of NPs for specific targeting, and the mean of NPs bound to the endothelial layer can help us understand the role of non-specific targeting.

We utilized Julia's GlobalSensitivity.jl package to perform Sobol sensitivity analysis (25). Sobol sensitivity captures the effect of variance in model inputs on model outputs in terms of Sobol indices. Namely, we look at first-order and total-order Sobol indices. First-order indices tell us about a single input parameter's contribution to the model output, and the total-order indices include the higher-order contributions of an input parameter.

Figure 5 contains the top four first-order and total-order Sobol indices for maximum NP uptake by organ tissue. It can be seen that NP uptake in organs other than the kidney and liver is affected mainly by parameters corresponding to the same organ, whereas in clearance organs such as the liver and kidney, the dominant parameters do not correspond to the same organ.

FIGURE 5

Figure 5. Global Sensitivity Analysis: Total order and First-order Sobol indices of K_deg, K_up, and K_NS parameters, for global sensitivity analysis performed on maximum NP uptake in organ tissue.

Figure 6 contains the top four first-order and total-order Sobol indices for the mean value of NPs in the endothelial compartment over time. The specific uptake depends on the binding of NPs to the endothelial layer, and the dominant parameters for this process provide insight into if the organ of interest is suitable for targeted drug delivery. If organ-specific parameters do not dominate the binding of NPs to the endothelial layer, then the optimization of targeted agents for that organ will rely on systemic factors and parameters of other organs. By way of this reasoning, the lungs can be an ideal targets, whereas it can be hard to design NPs for specific targeting of the heart, liver and kidney.

FIGURE 6

Figure 6. Global Sensitivity Analysis: Total order and First-order Sobol indices of K_deg, K_up, and K_NS parameters, for global sensitivity analysis performed on mean value of NPs bound to endohtelial layer.

2.6. Model validation

2.6.1. Computational performance metrics

The unbranched model was solved using MATLAB's ode15s solver for stiff systems and also using various stiff solvers from Julia's DifferentialEquations.jl Library. Both A and B models was run from 0 to 10,000 s (2.78 h) to reflect the time scale of a given experimental data set. The branched model was run from 0 to 10,000 s on Julia using the same solvers from DifferentialEquations.jl but only from 1 ms on MATLAB due to computational constraints. The simulations were performed on a 2020 MacBook Pro with 1.4 GHz Quad-Core Intel core i5 processor and 8 GB 2133 MHz LPDDR3 memory. The timestep of integration was optimized considering stability and the conservation of mass as two metrics. While the stiff solvers in Julia's DifferentialEquations.jl package do not require a user-defined timestep, for models A and B the timestep of integration, Δt, was 1 ms, and was 1 ns for the branched model, in MATLAB. All plots were generated through Julia.

2.6.2. Comparison to experimental data

To ensure that the model output can be considered accurate, the model was validated using data collected from an experimental study (23). This study used a variety of sizes of PEGylated gold nanoparticles and sampled the biodistribution at various time points (5, 30, 60, and 120 min) in five locations in a mouse. Using WebPlotDigitizer (https://apps.automeris.io/wpd/), attenuation values for 100 np NP were extracted from Figure 5 in Dong et al. (23) and converted to %ID/g values using the relationship presented in Figure 6b of Dong et al. (23). Then, using given organ weights from Diehl and Morse (21), %ID/t (percent of injected dose in whole tissue) was calculated. These %ID/t values determined experimentally in Dong et al. (23) were then compared to the model output to ensure the compartmental model is indeed predictive. Single time points for the liver and spleen were reported reliably, so these data sets were used for model validation. While this paper only describes validation with one experimental study, it is important to note that previous iterations of the model used a variety of other experimental studies for validation.

2.6.3. Mass conservation analysis

A conservation analysis was employed in every iteration of the model to ensure that each model described can be considered a closed system. The conservation curve of the model is the sum of NP that have been degraded (Equation 14) as well as those previous time points in every compartment of the model (arterial, venous, organs, vasculature, endothelial cell, and tissue) Equation (15). Mass is conserved in the model when the line formed by all data points from Equation (15) over time t has a slope of 0. Equation (16) is used for mass conservation analysis when we setup the system on a molar basis and the mass is conserved where the moles of NPs in all compartments add up to give initial value of NPs in venous compartment.

2.6.4. Stiffness of the unbranched model

The stiffness ratio characterizes the stiffness of the system. For a system of linear ordinary differential equations $\frac{d\overline{y}}{dt}=A\overline{y}$, The stiffness ratio is defined as the ratio of the largest and smallest eigenvalue of the matrix A.

$\begin{array}{ll}\text{Stiffness Ratio }=\frac{\left|{\lambda }_{max}\right|}{\left|{\lambda }_{min}\right|}.& (34)\end{array}$

Supplementary Table S6 shows the dependence of stiffness ratio on model parameters, based on the stiffness ratio of our system, we chose ode15s in MATLAB and QNDF, Rodas4, KenCarp4, TRBDF2 and RadauIIA5 from Julia's DifferentialEquations.jl package to solve our system. All the mentioned solvers are stiff ODE solvers.

2.7. Quasi steady state approximation

The system of ODEs described by our PBPK model is a stiff system. The stiffness in the system arises because of the large differences in order of magnitude of model parameters; in our system we have very fast binding and unbinding of NPs to endothelial layer and very slow uptake of NPs into tissue. The model can be made non-stiff by omiting the ODEs containing parameters at both time scales i.e. very fast (K_on, K_off) and very slow (K_UP, K_NS). A more qualitative explanation can be given by looking at the results from branched model simulation; it can be seen in the Supplementary Figure S1 that concentration profiles in the vasculature and endothelial layer quickly drop to zero after the initial spike. Therefore, the model can be approximated by omitting the differential equations related to the vasculature and endothelial layer. Consider the two subsystems

$\begin{array}{ll}\frac{d{N}_{ss}}{dt}=f(N_{nss},N_{ss})& (35)\end{array}$

and,

$\begin{array}{ll}\frac{d{N}_{nss}}{dt}=g(N_{nss},N_{ss}).& (36)\end{array}$

where N_ss denotes the moles of NP particles in vasculature and endothelial in all organs and N_nss denotes the moles of NPs in vein, artery and tissue compartment in all the organs. This reduced system consists of ODEs for concentration profiles in vein, artery and organ tissue (9 in total). The system is much less stiff than the complete systems because K_on and K_off not do occur in the reduced system of ODEs. Functions f and g consists of linear combinations of terms N_ss and N_nss for different organs, because the original system of ODEs is a linear system. Then the QSSA involves setting:

$\begin{array}{ll}\frac{d{N}_{ss}}{dt}=0.& (37)\end{array}$

We then use the values of N_ss obtained by solving equation (37) as constants in equation (36), this is no different than solving a system of ODEs in Julia with an extra step of solving a system of linear (Equation 37). Since the reduced system is nonstiff equation (36) can be solved with any nonstiff/explicit ODE solver.

2.8. Neural networks to solve ODEs

It has been shown in Raissi et al. (26) that Neural Networks can be used to solve partial differential equations; we use the same protocol to solve our system of ODEs using a neural network. Our system being a stiff one makes it ideal for testing the performance of a new method and comparing it against highly optimized ODE solvers. It has been shown (27) that neural networks can be used to solve ODEs when the system of ODEs is nonstiff. We used QSSA to make our PBPK model nonstiff (stiffness ratio ~10⁴) and trained a neural network to learn the solution to the ODEs.

Consider a system of first-order ODEs; $\frac{dy}{dt}=f(y)$, where both f and y are vectors of same size. Hypothesize that the solution can be approximated using a neural network, with trainable parameters θ, i.e.

$\begin{array}{ll}\text{y}(t)=\text{N}\text{N}(t,\theta ).& (38)\end{array}$

The derivative of the output from the neural network can be approximated by a finite difference scheme or can be obtained using autograd functionality of a neural network. Here we use a first-order finite difference scheme.

$\begin{array}{ll}\frac{d\text{y}(t)}{dt}=\frac{\text{d NN}(t,\theta )}{dt}=\frac{\text{N}\text{N}(t+ϵ,\theta )-\text{N}\text{N}(t,\theta )}{ϵ}.& (39)\end{array}$

The loss function for this problem is:

$\begin{array}{ll}L(\theta )={\displaystyle \sum _{t_i\in t}{\left[\frac{\text{d}NN({\text{t}}_{\text{i}},\theta )}{dt}-f(NN(t_i,\theta ))\right]}^2}+[NN(0,\theta ))-y_{ic}{]}^2.& (40)\end{array}$

The first term in the loss function is the squared error between LHS and RHS of the ODEs, computed using forward pass through the Neural Networks and the second term is the squared error between true initial condition and predicted initial condition. The solution NN(t, θ^*) is a unique solution to the system of ODEs, where,

$\begin{array}{ll}{\theta }^{*}={\displaystyle \underset{\theta }{\text{argmin}}}L(\theta ).& (41)\end{array}$

Equation (41) is the optimization problem aiming to minimize the above defined loss in terms of neural network parameters θ. The procedure described in Algorithm 1 aims to solve Equations (35) and (36) iteratively. Notations used in the algorithm mean the same as described in these equations. We solved our PBPK model with a neural network using Julia's Flux.jl package (35).

ALGORITHM 1

Algorithm 1. Physics Informed Neural Network with QSSA

3. Results

3.1. Biodistribution from steady state model agrees with murine in vivo data

Using the steady-state model equation (Equation (2) in Methods), we plot the %idg values in each of the five organs over three different antibody concentrations on the surface of NP (41, 100, and 162 antibodies coating the surface of the NP), four different theoretical membrane properties (flat, membrane, membrane & present macrophages, membrane & present and active macrophages), and five different organs (lung, heart, kidney, liver, spleen). Additionally, the in-vivo experimental %idg data is plotted, in Figure 7 and R values representing the correlation between the lung %idg values of the given model and experimental lung values, as well as R values that show the correlation between the %idg values of all organs of the model output and the experimental data are given. It is important to note that the total R values across all organs are much higher in the modified model discussed in this paper compared to the original model, comparisons of R values in both models is shown in Supplementary Table S5. This suggests that the incorporation of non-specific uptake into the model results in a much more physiologically relevant model. We also compare the results obtained from a 30-min simulation of an unbranched model with experimental data (reported at a 30-min time-stamp). We observe that by performing local sensitivity analysis for area and volume of the endothelial layer and organ tissue, respectively, we obtain a high R for lungs compared to all other models. The excellent agreement of the steady-state data with the 30-min time data from the transient model suggests that the steady-state model is a good approximation for time scales.

FIGURE 7

Figure 7. Modified model from Ramakrishnan et al. (15) that includes non-specific uptake. % idg reported for several membrane conditions, antibody concentrations, and organ combinations.

3.2. Compartmental model B provides a better representation of the temporal biodistribution observed in experiments

Both compartmental models (A and B) were solved using a variety of stiff solvers available in Julia's DifferentialEquations.jl package and the MATLAB ode15s solver (for stiff systems) for a time period of 10,000 s (2.78 h). The output and corresponding conservation analysis graphs of models A and B are shown in Figures 8, 9, respectively.

FIGURE 8

Figure 8. Output of model A: (A) the normalized concentration of NP in each of the five organ compartments and the arteries and veins of the original model, and (B) the conservation figure.

FIGURE 9

Figure 9. Output of model B: (A) the normalized concentration of NP in each of the seven organ compartments and the arteries and veins of the model, and (B) the conservation figure.

The normalized root mean squared deviation values (NRMSD) were calculated comparing the model output to the experimental data set where the experimental data set provided time-series data (kidney, liver, spleen). The NRMSD values in model A were 1.0503, 37.4273, and 0.1025 for the kidney, liver, and spleen, respectively. The NRMSD values in model B were 0.6735, 20.1909, and 0.0614 for the kidney, liver, and spleen, respectively. The NRMSD values were lower in model B than model A suggesting the modified model (model B) more accurately represents the experimental data set (23).

The ODE solvers in Julia's DifferentialEquations.jl package does not require user-defined timestep and exact mass conservation was obtained by using any of the available stiff solvers. MATLAB's ode15s solver requires a user-defined timestep (Δt). Timestep was chosen so the system remains stable, and so the system exhibits mass conservation. The Δt of 0.001 s was chosen for both compartmental models A and B, for the simulation time of 10,000 s. If the Δt was increased beyond 0.001 s, the system became unstable and mass was not conserved. On the other hand, if Δt was decreased, the model was unable to complete running due to the computing constraints of MATLAB. It is important to note that mass conservation is only accurate to the order of Δt but integration is valid to a higher order. So, to run the model for a time scale similar to that of experimental studies while maintaining stability and conservation in the system, a Δt of 0.001 for models A and B is necessary.

3.3. Branched model predicts a delayed temporal response compared to the lumped compartmental model

The development of the branched model is incredibly important because it provides a more accurate representation of the circulatory system; specifically, the surface area to volume ratio of the blood vessels in the branched model is more representative of in vivo mouse models. The framework of the branched model allows for the incorporation of more specific hydrodynamic interactions and margination allowing for a physics-based prediction of NP biodistribution, enabling us to account for characteristics such as NP size, shape, and surface chemistry in the model, following the theory in Jabeen et al. (28). Additionally, the construction of the branched model will allow for K_on to change depending on vessel diameter and blood flow rate. This cannot be done in models A and B without empirical measurements.

The branched model was solved using Julia's DifferentialEquation.jl (29) package. Namely, the stiff solvers, QNDF, Rodas4, KenCarp4, TRBDF2 and RadauIIA5 were able to solve for t = [0, 10⁴]s. Whereas the most efficient stiff solver in MATLAB (ode15s) was unable to solve the system of ODEs. Figure 10 shows the comparison between the molar profiles for the branched and unbranched models.

FIGURE 10

Figure 10. NP concentration vs. Time for the branched and unbranched models: For the organs which have φ>1 the onset is quicker. For lungs φ~1 and for spleen φ < 1.

3.4. Effect of nanoparticle size on biodistribution

The branching model construction allows for exploration of the effect of differing NP sizes on biodistribution. The effect of nanoparticle size on biodistribution was also explored using the branched model and nanoparticle diameters of 4, 15, 50, 79, and 100 nm. In Figure 11, it can be seen that the 5-min simulation result from our branched model is in good agreement for 5-min experimental data from Dong et al. (23) for the kidney. Even though we were not able to get an exact match for Liver and Spleen, we observed similar trends between experimental data and simulation data.

FIGURE 11

Figure 11. Effect of NP size: Y axis represents relative concentration of nanoparticles at 5 min w.r.t initial concentration.

3.5. Model reduction and performance of newer solvers

3.5.1. Quasi steady state approximation for the unbranched model

The unbranched model consists of 23 ordinary differential equations, i.e., two equations for vein/artery, seven equations for branched vasculature, seven equations for branched endothelial layer, and seven equations for organ tissues. The reduced system after using QSSA includes only nine equations. The system of linear equations for vasculature and endothelial layer (14 equations) is solved using the similar approach.

Figure 12 depicts the comparison between the unbranched model's QSSA solution and the unbranched model's complete solution. The difference in the solutions can be explained based on the vasculature and endothelial profiles for the unbranched model. The difference in profiles for some organs (Kidney, Liver, Gut, Others) is due to the non-zero gradient in profiles of NP bound to endothelial and vasculature for large times, as shown in Supplementary Figure S2.

FIGURE 12

Figure 12. Comparison between complete model and QSSA for unbranched model.

3.5.2. QSSA for the branched model

The full branched model consists of 457 ordinary differential equations, i.e., two equations for vein and artery, 32 x 7 equations for branched vasculature, 32 x 7 equations for branched endothelial layer, and seven equations for organ tissues. The reduced system after using QSSA has nine equations. The system of linear equations for vasculature and endothelial layer (64*7 equations) is solved using the backslash (A\b) operator in Julia, and the values of y are updated in this way at every iteration. However, this makes the task computationally more intensive and takes longer than solving for the complete model, but the objective is to reduce the system's stiffness. After using QSSA the system of ODE was solved using Tsit5 solver in Julia, which is a nonstiff solver. Figure 13 shows the comparison between QSSA and the complete solution of the branched model. The difference in initial onset can be explained from Supplementary Figure S1, where we can see the gradient is zero for most of the time but there is a spike initially. QSSA for the branched model is a better approximation compared to QSSA in the unbranched model because of the more rapidly vanishing gradients in the former's case.

FIGURE 13

Figure 13. Comparison between complete model and QSSA for branched model.

3.5.3. Coupling QSSA and neural networks

Using the methods described in the above sections, we solved the concentration profile in the tissue, vein and artery compartments. The neural network predicts these nine outputs, which are used to update the steady-state solution of the remaining 14 equations iteratively.

The neural network was trained for 30,000 epochs on a CPU. The neural network architecture consists of an input layer to which discretized time is given as input, two hidden layers each followed by a hyperbolic tangent activation and an output layer consisting of 9 outputs followed by Sigmoid activation. The input to the neural network is a batch of equally spaced numbers between [0, 1], and the first-order derivative of the neural network is scaled with maximum time ($t_{max}=10^3$) before computing the loss function.

Figure 14 depicts the comparison between the Neural Network solution and the solution obtained using a nonstiff ODE solver (Tsit5 in Julia). This method demonstrates the ability of Neural Networks to solve ODEs. However, we tested this implementation of a simple feedforward neural network to solve the system of ODEs using QSSA. For the full model (without QSSA), we found that the ODE solvers from Julia's DifferentialEquations.jl library are more adept and a neural network fails to solve the full system.

FIGURE 14

Figure 14. Comparison between Neural Networks and Tsit5 with QSSA for unbranched model.

4. Discussion

In clinical settings, the use of nanotechnology, including drug-carrying nanoparticles (NP), has increased in recent years. However, the range of NP applications, target, and physical characteristics significantly impede the ability of NP to be researched effectively as bench-to-bedside therapeutics. To this end, researchers have begun to turn toward physiologically based pharmacokinetic (PBPK) models to guide experimentation and better understand the targeting behavior of various nanoparticle compositions in the human body. A multiscale computationally driven model with physiologically relevant inputs can be utilized to determine organ-specific biodistribution since the physiological and hydrodynamic factors governing NP biodistribution and tissue targeting involve mechanisms that operate at different timescales. NP behavior must be understood at every level to create a comprehensive multiscale model. This includes the binding landscape of a NP in the presence of an endothelial cell layer. A previous multiscale PBPK model has determined binding constants of intracellular adhesion molecule 1 (ICAM1) coated NPs to endothelial cell surface receptors in mice and humans by utilizing the biophysical properties of the antibody to receptor interactions, and the cell surface (15). However, the nonspecific uptake, or uptake via passive diffusion in the intercellular cleft, is not accounted for in that model. The purpose of this current study was to 1) modify an existing steady-state PBPK model (15) to incorporate NP uptake via nonspecific transport, 2) develop a novel multiscale PBPK compartmental model to predict temporal effects, and 3) introduce a compartmental branched vascular model that can predict the effect of NP size, 4) perform validation with experimental murine biodistribution data.

The original steady-state model from Ramakrishnan et al. (15) was modified by adding a nonspecific uptake term that represents NP uptake via passive diffusion through the intracellular cleft. The addition of nonspecific uptake increased the predictive ability of the model, as evidenced by higher R values when compared to an experimental data set than the original model. This suggests that incorporating a nonspecific uptake term into future models is necessary to increase physiological relevance and predictive ability. Next, a novel multiscale PBPK compartmental model was created to predict the continuous temporal biodistribution of NP in five to seven organs. Two versions of this model were created, model A and model B, which differ based on the number of organ compartments represented (model A: 5, model B: 7), and the way flow is routed through the model (Model B is more physiologically relevant). Models A and B were represented with a system of ODEs (Model A: 17 ODEs, Model B: 23 ODEs) solved using stiff solvers from Julia's DifferentialEquations.jl package and MATLAB's ode15s solver, then validated with experimental data. The predictive ability of Model B was greater than Model A as evidenced by the normalized root mean squared deviation (NRMSD) analysis. Finally, a branched model was developed to create a more detailed and physiologically relevant version of the basic compartmental model while still maintaining the simplistic compartmental foundation. The branched network consists of a branched vascular tree that begins at the main arteries and veins and bifurcates into the capillary beds, connecting the arterial and venous branching networks. This branching model was represented with 457 ODEs which were solved using Julia. The branching framework allowed for customized output based on NP size. K_on and K_off values can be calculated and are dependent on NP size. Biodistribution for NP size of 4, 15, 50, 79, and 100 nm was plotted.

The compartmental models' A and B differed slightly in their composition, with model B more physiologically relevant. Model A consists of five organ compartments (lungs, heart, kidneys, liver, spleen), while model B includes the addition of two additional compartments (gut, other). Additionally, the lung circulation has been separated out to keep track of oxygenation and oxygen distribution, and gut and spleen compartments are coupled to the liver compartment in model B. A local sensitivity (Figure 4) analysis was performed on several parameters: K_deg, K_up, and K_NS to determine the values of these parameters that would best fit the experimental data set available. K_up and K_NS had a similar effect on the model output across organ compartments. When increased, the slope of the biodistribution curve over time would increase, with a higher plateau concentration. When K_deg is increased, the slope of the biodistribution curve over time would decrease, typically decreasing the plateau concentration and also shortening the length of the plateau, leading to a quicker decrease in NP concentration in a given compartment.

The multiscale PBPK model framework presented in this paper presents a significant advance because the predictive ability is purely mechanism-based. Multiscale physics-based modeling allows for the system's behaviors and interactions to be completely described mathematically, rather than relying on empirical observations and data to make predictions. Typical PBPK models (7, 9) are generally empirically based and do not describe the entire behavior of the system. Creating a purely physics-based PBPK model allows for more customization. For example, in this case, it allows NP composition to be varied by changing certain model parameters to reflect differing NP surface chemistry or size. This is advantageous since NP exist in many forms with various surface chemistries, compositions, and sizes and allows for further model customization in the future.

To continue to evolve this model to allow for customization beyond NP size, additional modules can be added in the future. Incorporating a module that defines internalization rates of varying NP surface chemistries is vital to extending this model to other NP besides ICAM coated. Additionally, hydrodynamic interactions in the bloodstream, immune system effects, and separating various types of NP degradation can be added to the model to increase the physiological relevance and translational potential.

The branching model results in an incredibly large and stiff system of ODEs. To attempt to combat these issues, the stiff MATLAB ode solver, ode15s, was used. MATLAB produced biodistribution graphs, but only from 0 to 1 ms. This is because a small time-step needed to be used to ensure stability within the model. If the time-step was increased beyond 1 ns, the system became unstable, and if the run time was increased beyond 1 ms with the 1 ns time-step, MATLAB became unresponsive. So, it is clear that there are computing power limitations in MATLAB solving large stiff ODE systems. We used stiff equation solvers from DifferentialEquations.jl package in Julia to address this issue. All the stiff solvers from the package successfully solved the system large times.

The physiologically relevant PBPK model can produce a correct output describing the biodistribution of NP. While the neural net ODE solver was successfully demonstrated here, the full power of neural networks can be realized by embedding the multiphysics in the neural networks. Eventually, this computationally driven PBPK model could be used for developing a Neural Network to create a reduced but accurate model for determining NP biodistribution, with more efficiency than the PBPK compartment model. An input to the neural network could be characteristics of the NP such as NP size, vesicle cargo, and concentration of surface proteins. The discrete NP concentrations in each organ as determined by the neural network could be trained against the continuous PBPK compartmental model output (described in this study) at a given time point. Utilizing ML techniques in this model will allow for a much more efficient and automated predictive model for determining NP biodistribution. However, it is necessary to construct an informative PBPK compartmental model to ultimately use in the training process of this Neural Network (26, 27).

4.1. Limitations and future work

• While our model does not include charge effects, parameter sensitivity is a simple way to address this before an in-depth study. Indeed, we have carried out local and global sensitivity analyses on many of the parameters to identify the sensitive parameters. The model is context-specific, so if implemented for a different type of nanoparticle, all that would be required are changes to model parameters.

• We chose to study nanoparticles in a size range of between 10 and 100 nm because experiments reporting temporal distribution (and not just single time points) were available in this size range. The studies that reported larger particles tended to focus on effects at single time points. In principle, our studies can describe larger nanoparticles, but we need to know how the parameters such as K_on, K_off, and K_up change for these larger particles. We note that the dependence of these parameters on size for the 10–100 nm range was available through our earlier studies (30); these studies also reported results for larger particles such as 500 nm. Therefore, in principle, our studies can describe larger nanoparticles so long as we know how the parameters such as K_on, K_off, and K_up change for these larger particles. The smaller size range we considered is more appropriate for translational applications as described in Anselmo and Mitragotri (31), which reports that the nanoparticles used are primarily in the range we have considered.

• We used the Matlab ODE 15s solver in this study due to its ease of implementation in Matlab. We also tried other solvers, such as ODE 45s and ODE 23tb but each posed problems (inefficiency). We tried using the Matlab ODE 45s solver in earlier iterations of this project. However, the ODE 45s solver was significantly slower to solve the system than the ODE 15s solver. Solving was never actually completed with the ODE 45s solver due to its inefficiency. We want to emphasize that it was not our intent to compare Matlab and Julia systematically. We undertook a limited comparison presented in this paper out of necessity and concluded that in the situations where Matlab failed, Julia was able to provide a viable and computationally tractable solution. However, a more systematic comparison and depth analysis of performance have been reported in the literature, and our observations are consistent with these reports (32).

• We believe our model presents a minimal framework to interpret salient features of systemic NP transport. This model can be applied to any other species by changing the model parameters and validating against experimental data for the particular species. However, the model is (woefully) inadequate to capture the full complexity of the physiological system. Hence it is our philosophy that we will explore the dynamic range of the model predictions through sensitivity (local and global), evolvability, and robustness analysis (35), to the extent that the model results align with physiological measurements. Our model offers one plausible interpretation; we can confirm this by additional constraints such as exploring the effects of critical parameters in the model and experiments to ensure the model predicts the correct responses for the right reasons. Additionally, validation can be performed at multiple scales by carrying out independent experiments (33). These multiple steps of statistical analysis, sensitivity analysis, and multiscale experimental validation must be done before determining if the minimal framework of the model is entirely adequate for physiological comparison. At this point, clinical adaptation or at least adaptation in translational settings can be attempted. If the model fails these tests in a given scenario, it is most often because one or more crucial effects are missing from the minimal framework. At this point, the model needs to be expanded, and further validation of the expanded model is warranted before further use.

• We attempted the simpler model A to determine if the simple wiring of the organs performs as well as the more physiologically correct wiring. While model B outperforms model A, model A gets most of the constraints correct, implying that most effects governing tissue accumulation depend on the flow rates and the compartmental volumes. The flow rates, for the most part, are only partially sensitive to changes in wiring. The purpose of the stripped-down versions of the model is to assess the accuracy of prediction vs. model complexity to determine what level of complexity of the model provides an acceptable description of the results. Moreover, this acceptable limit is set by a user threshold of how much error tolerance we can accommodate in our predictions.

• Several papers in the journals described by the reviewer either do not report temporal data, or they are reported for nanoparticles which are not targeted; also, very few studies do all this and study the effect of nanoparticle size. A physics-based (or data-driven) understanding of the dependence of different nanoparticle architecture on key parameters is required before they can be factored into our model. Hence we chose to focus on (limit ourselves to) solid spherical nanoparticles in the current study. We hope to explore more detailed comparisons of other classes of nanoparticles in the future. We present this perspective under the limitations in the discussion section. One drawback of physics-based models is that while the models are generalizable, they are predicated on physics being known for different classes of particles. So far, we have explored the effect of the particle architecture for three classes—rigid, flexible, and semiflexible (34), and hope to cover other particles in the future. At this time, we will systematically explore the literature, uncover the temporal data reported across nanoparticle classes, and subject them to comparison to our model. However, we have not done this step for the current version of our study.

• Our model is minimal, and within this minimal framework, we have backed up the model parameters and their dependence on essential effects such as size and antibody density on previous physics-based simulations. In this realm, we have considered cell-membrane-NP interactions in K_on, K_off, and K_up. These bundle the different uptake mechanisms and rates. The other mechanisms of non-specific uptake are bundled into K_ns. In principle, systematic studies of each of these uptake routes while blocking others can be done and would be insightful to pursue in future studies.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material. The MATLAB and Julia code (which includes data used to run the models) for this study can be found in the Multiscale PBPK Nanoparticle Biodistrubtion GitHub repository (https://github.com/emmamglass/Multiscale-PBPK-Nanoparticle-Biodistribution-Model).

Author contributions

EG, SK, and RR coordinated and led the research, analyzed the data, and wrote the manuscript. CE, SF, and AM contributed to crucial sections and specific aspects of the codes and analysis. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the NSF-funded Science and Technology Center, Center for Engineering Mechanobiology (CEMB), award number CMMI-1548571, by the Penn Engineering Blair Scholar Fellowship Fund, the NIH Physical Sciences in Oncology/ Cancer Systems Biology Consortium Summer Internship, and the William and Mary Charles Center Honors Research Fellow Program. This study has received funding from the National Institutes of Health under U01CA250044.

Acknowledgments

We thank David Cormode for providing us with experimental validation data, and Greg Conradi Smith for aiding in the continuation of this work at The College of William and Mary.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmedt.2022.934015/full#supplementary-material

References

1. Blanco E, Shen H, Ferrari M. Principles of nanoparticle design for overcoming biological barriers to drug delivery. Nat Biotechnol. (2015) 33:941–51. doi: 10.1038/nbt.3330

PubMed Abstract | CrossRef Full Text | Google Scholar

2. He H, David A, Chertok B, Cole A, Lee K, Zhang J, et al. Magnetic nanoparticles for tumor imaging and therapy: a so-called theranostic system. Pharm Res. (2013) 30:2445–58. doi: 10.1007/s11095-013-0982-y

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Anselmo A, Mitragotri S. Nanoparticles in the clinic. Bioeng Transl Med. (2016) 1:10–29. doi: 10.1002/btm2.10003

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Mitragotri S, Lammers T, Bae YH, Schwendeman S, De Smedt S, Leroux JC, et al. Drug delivery research for the future: expanding the nano horizons and beyond. J Control Release. (2017) 246:183–4. doi: 10.1016/j.jconrel.2017.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hua S, de Matos MBC, Metselaar JM, Storm G. Current trends and challenges in the clinical translation of nanoparticulate nanomedicines: pathways for translational development and commercialization. Front Pharmacol. (2018) 9:790. doi: 10.3389/fphar.2018.00790

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Eckmann DM, Bradley RP, Kandy SK, Patil K, Janmey PA, Radhakrishnan R. Multiscale modeling of protein membrane interactions for nanoparticle targeting in drug delivery. Curr Opin Struct Biol. (2020) 64:104–10. doi: 10.1016/j.sbi.2020.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Li M, Al-Jamal KT, Kostarelos K, Reineke J. Physiologically based pharmacokinetic modeling of nanoparticles. ACS Nano. (2010) 4:6303–17. doi: 10.1021/nn1018818

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lin Z, Monteiro-Riviere NA, Riviere JE. A physiologically based pharmacokinetic model for polyethylene glycol-coated gold nanoparticles of different sizes in adult mice. Nanotoxicology. (2016) 10:162–72. doi: 10.3109/17435390.2015.1027314

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yuan D, He H, Wu Y, Fan J, Cao Y. Physiologically based pharmacokinetic modeling of nanoparticles. J Pharm Sci. (2019) 108:58–72. doi: 10.1016/j.xphs.2018.10.037

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Rao GG, Landersdorfer CB. Antibiotic pharmacokinetic/pharmacodynamic modelling: MIC, pharmacodynamic indices and beyond. Int J Antimicrob Agents. (2021) 58:106368. doi: 10.1016/j.ijantimicag.2021.106368

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Snoek E, Jacqmin P, Sargentini-Maier ML, Stockis A. Modeling and simulation of intravenous levetiracetam pharmacokinetic profiles in children to evaluate dose adaptation rules. Epilepsy Res. (2007) 76:140–7. doi: 10.1016/j.eplepsyres.2007.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Sharma RP, Schuhmacher M, Kumar V. The development of a pregnancy PBPK model for bisphenol a and its evaluation with the available biomonitoring data. Sci Total Environ. (2018) 624:55–68. doi: 10.1016/j.scitotenv.2017.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Husoy T, Martinez MA, Sharma RP, Kumar V, Andreassen M, Sakhi AK, et al. Comparison of aggregated exposure to di(2-ethylhexyl) phthalate from diet and personal care products with urinary concentrations of metabolites using a PBPK model - Results from the Norwegian biomonitoring study in EuroMix. Food Chem Toxicolo. (2020) 143:111510. doi: 10.1016/j.fct.2020.111510

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Li M, Nguyen L, Subramaniyan B, Bio M, Peer CJ, Kindrick J, et al. PBPK modeling-based optimization of site-specific chemo-photodynamic therapy with far-red light-activatable paclitaxel prodrug. J Control Release. (2019) 308:86–97. doi: 10.1016/j.jconrel.2019.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ramakrishnan N, Tourdot R, Eckmann D, Ayyaswamy P, Muzykantov V, Radhakrishnan R. Biophysically inspired model for functionalized nanocarrier adhesion to cell surface: roles of protein expression and mechanical factors. R Soc Open Sci. (2016) 3:160260. doi: 10.1098/rsos.160260

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ayyaswamy PS, Muzykantov V, Eckmann DM, Radhakrishnan R. Nanocarrier hydrodynamics and binding in targeted drug delivery: challenges in numerical modeling and experimental validation. J Nanotechnol Eng Med. (2013) 4:011001. doi: 10.1115/1.4024004

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Radhakrishnan R. A survey of multiscale modeling: foundations, historical milestones, current status, and future prospects. AIChE J. (2020) 67:17026. doi: 10.1002/aic.17026

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Agrawal N, Radhakrishnan R. The role of glycocalyx in nanocarrier-cell adhesion investigated using a thermodynamic model and monte carlo simulations. J Phys Chem C Nanomater Interfaces. (2007) 111:15848–56. doi: 10.1021/jp074514x

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Radhakrishnan R, Yu HY, Eckmann DM, Ayyaswamy PS. Computational models for nanoscale fluid dynamics and transport inspired by nonequilibrium thermodynamics1. J Heat Transfer. (2016) 139:033001. doi: 10.1115/1.4035006

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yang J, Wang Y. Design of vascular networks: a mathematical model approach. Int J Numer Method Biomed Eng. (2013) 29:2534. doi: 10.1002/cnm.2534

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Diehl L, Morse M. A Comparison of Selected Organ Weights and Clinical Pathology Parameters in Male and Female CD-1 and CByB6F1 Hybrid Mice 12-14 Weeks in Age. Charles River, MA: Spencerville (2015).

22. Boswell CA, Mundo EE, Ulufatu S, Bumbaca D, Cahaya HS, Majidy N, et al. Comparative physiology of mice and rats: radiometric measurement of vascular parameters in rodent tissues. Mol Pharm. (2014) 11:1591–8. doi: 10.1021/mp400748t

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Dong Y, Hajfathalian M, Maidment P, Hsu J, Naha P, Si-Mohamed S, et al. Effect of gold nanoparticle size on their properties as contrast agents for computed tomography. Sci Rep. (2019) 9:14912. doi: 10.1038/s41598-019-50332-8

PubMed Abstract | CrossRef Full Text | Google Scholar

24. McKenzie M, Ha S, Rammohan A, Radhakrishnan R, Natesan R. Multivalent binding of a ligand-coated particle: role of shape, size, and ligand heterogeneity. Biophys J. (2018) 114:1830–46. doi: 10.1016/j.bpj.2018.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang XY, Trame MN, Lesko LJ, Schmidt S. Sobol sensitivity analysis: a tool to guide the development and evaluation of systems pharmacology models. CPT. (2015) 4:69–79. doi: 10.1002/psp4.6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys. (2019) 378:686–707. doi: 10.1016/j.jcp.2018.10.045

CrossRef Full Text | Google Scholar

27. Ji W, Qiu W, Shi Z, Pan S, Deng S. Stiff-PINN: physics-informed neural network for stiff chemical kinetics. J Phys Chem A. (2021) 125:8098–106. doi: 10.1021/acs.jpca.1c05102

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Jabeen Z, Yu HY, Eckmann DM, Ayyaswamy PS, Radhakrishnan R. Rheology of colloidal suspensions in confined flow: treatment of hydrodynamic interactions in particle-based simulations inspired by dynamical density functional theory. Phys Rev E. (2018) 98:042602. doi: 10.1103/PhysRevE.98.042602

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Rackauckas C, Nie Q. Differentialequations.jl-a performant and feature-rich ecosystem for solving differential equations in julia. J Open Res Softw. (2017) 5:15. doi: 10.5334/jors.151

CrossRef Full Text | Google Scholar

30. McKenzie M, Ha SM, Rammohan A, Radhakrishnan R, Ramakrishnan N. Multivalent binding of a ligand-coated particle: role of shape, size, and ligand heterogeneity. Biophys J. (2018) 114:1830–46.

PubMed Abstract | Google Scholar

31. Anselmo AC, Mitragotri S. Nanoparticles in the clinic: an update. Bioeng Transl Med. (2019) 4:e10143. doi: 10.1002/btm2.10143

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Rackauckas C. A Comparison Between Differential Equation Solver Suites. In Matlab, R, Julia, Python, C, Mathematica, Maple, and Fortran. The Winnower (2018).

33. Liu J, Weller GER, Zern B, Ayyaswamy PS, Eckmann DM, Muzykantov VR, et al. Computational model for nanocarrier binding to endothelium validated using in vivo, in vitro, and atomic force microscopy experiments. Proc Natl Acad Sci USA. (2010) 107:16530–5. doi: 10.1073/pnas.1006611107

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Farokhirad S, Kandy SK, Tsourkas A, Ayyaswamy PS, Eckmann DM, Radhakrishnan R. Biophysical considerations in the rational design and cellular targeting of flexible polymeric nanoparticles. Adv. Mater. Interfaces (2021) 8:2101290. doi: 10.1002/admi.202101290

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Innes M, Saba E, Fischer K, Gandhi D, Rudilosso MC, Joy NM, et al. Fashionable modelling with flux. CoRR. (2018) abs/1811.01457. doi: 10.48550/arXiv.1811.01457

CrossRef Full Text | Google Scholar