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.

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 Di erentialEquations.jl in Julia outperforms MATLAB when solving a sti system of ordinary di erential 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-sti via quasi-steady-state approximation (QSSA). Our model incorporates modules that account for varying NP surface chemistries, multiscale vascular hydrodynamic e ects, and e ects of the immune system to create a more comprehensive and modular model for predicting NP biodistribution in a variety of NP constructs. KEYWORDS biodistribution, physiologically based pharmacokinetic, quasi steady-state approximation, neural network, vascular branching

. Introduction
Clinical medicine has entered an era where nanotechnology is becoming more and more prevalent, with the use of drugcarrying 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)(7)(8)(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)(13)(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 steadystate 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.

. . 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) 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 receptormediated internalization and non-specific transport through pores between the intercellular cleft of adjacent endothelial cells, especially in clearance organs. Here, the existing steadystate PBPK model was modified to incorporate this nonspecific 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 2 ) 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), 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): 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.
. . Temporal model . . . 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 . /fmedt. . 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). 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.

. . . 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)(4)(5)(6)(7), where i = lung, heart, kidneys, liver, spleen. The venous compartment is described by: (3) 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 vein deg is the degradation rate of NP in the veins. The arterial compartment is described by: 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 art deg is the degradation rate of NP in the arteries. For each tissue type i, the vascular compartment is described by: 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 i bl and C i EC denotes the concentration of NP in the vascular and endothelial cell compartments of the organ, respectively, K i on denotes the rate of binding of NP to the endothelial cell surface receptors, K i off denotes the rate of NP unbinding from the endothelial cell surface receptors, and K i NS denotes the rate of nonspecific NP uptake into the organ tissue. The endothelial compartment in each tissue is described by: Equation (6) describes the change in concentration of NP bound to the endothelial cell surface receptors over time, where K i up denotes the rate of uptake of NP into the organ tissue via transcytosis. Each tissue compartment is described by: Equation (7) describes the change in concentration of NP in the organ tissue over time where C i T and V i T denote the concentration of NP in the tissue compartment of the organ and the volume of the tissue compartment, respectively.

. . . 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)(9)(10)(11)(12)(13) 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) V art dC art dt = Q art C heart bl − C art (Q lung + Q heart 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 heart bl is the concentration of NP in the vascular compartment of the heart.
Equation (11) describes the change in concentration of NP in the vascular compartment of the heart over time.
Equation (12) 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.
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 t total over every time point t of the simulation should have a slope of 0.
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.

. 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 powerlaw 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.
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 0 and daughter vessels are d 1 and d 2 , and k=3, which is typically assumed based on a minimum dissipation principle representing a stable flow condition: where d 0 > d 1 = d 2 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 n . 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.
Then, the length of each vessel can be determined using the know length to diameter ratio β = 3, and the diameter of the vessel, d.
This network will begin at the diameter of the main artery d 0 = 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 lowerbound 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.
Equation (20) can be used to determine the total vascular volume of a branching element, 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.

. . . 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)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(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). = Q heart C bl,32 Equation (22) describes the NP concentration in generation n = 1 in the vascular branching network of the lung. C bl,1 lung denotes the NP concentration of the lung vasculature in generation n = 1 of the branching network, Q 1 lung represents the flow rate through generation n = 1 of the lung vascular network, C bl,32 heart denotes the NP concentration of the heart vasculature in generation n = 32 (the concentration of NP exiting the heart compartment), N 1 represents the number of branches in generation n = 1, ϕ lung represents the total number of branching elements in the given organ compartment, V bl,1 lung denotes the vascular volume of a single branching element in generation n = Q n−1 lung C bl,n−1 lung N n−1 − Q n lung C bl,n lung N n −ϕ lung C bl,n lung V bl,n lung (K on Equation (23) describes the NP concentration in the lung vasculature in generation n = 2 to n = 16 of the branching network. Q n−1 lung , C bl,n−1 lung , and N n−1 denote the flow rate, NP concentration, and number of branches in the lung vasculature of the previous generation of the branching network.
Equation (24) 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.
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. V bl,32 heart dC bl, 32 heart dt = Q 31 heart C bl,31 heart N 32 − Q heart C bl,32 heart + Q vein C vein −ϕ heart C bl,32 heart V bl,32 heart (K on heart + K NS heart + K deg heart ) K off heart φ heart C EC,32 heart V EC,32 heart .
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.
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) 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.
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. 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.

. 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 i bl . The volume of blood in the mouse vein was computed using the following equation, V total = V art + V vein + V i bl ,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 i T . 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), . Given φ= 0.3, D EC = 5 × 10 −6 , and the V T for each mouse organ previously described, l 2 EC 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 i EC . The binding rate of NP to the endothelial cell surface (K i on ) was determined using the relationship D l 2 given in (15). The unbinding rate of NP from the endothelial cell surface (K i off ) can be determined by using K i on 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 i off is calculated as a ratio of K i on to K i EC ; 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.

. . . 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 = D were changed incrementally so the output curve closely resembled that of an experimental data set (discussed in Section 2.5.1). Larger K i NS and K i up values resulted in a steeper initial biodistribution curve. Generally, K i NS and K i up had a similar magnitude result in the biodistribution curve. Larger K i deg 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.
Following the sensitivity analysis in model A, model B, and the branching model, values for K deg i , K NS i , and K up i 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.

. . . 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 totalorder 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 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 organspecific 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.

. . Model validation . . . 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 userdefined 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.

. . . 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

. . . 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.

. . . Sti ness of the unbranched model
The stiffness ratio characterizes the stiffness of the system.

. . 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 .
/fmedt. .   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: 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.

. . 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 4 ) and trained a neural network to learn the solution to the ODEs.

Consider a system of first-order ODEs;
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. NN(t, θ ).

y(t) =
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.
The loss function for this problem is: (40) 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, 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). , 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.

. . 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.
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.

. . 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 4 ]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.

. . E ect 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 Frontiers in Medical Technology frontiersin.org . /fmedt. . Spleen, we observed similar trends between experimental data and simulation data.
. . Model reduction and performance of newer solvers . . . 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.    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.
. . . 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, .
/fmedt. . 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.

. 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  (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 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).

. . Limitations and future work
• While our model does not include charge effects, parameter sensitivity is a simple way to address this before an indepth 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 . /fmedt. . 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.