Physiologically Based Multiphysics Pharmacokinetic Model for Determining the Temporal Biodistribution of 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). In the future, this model will incorporate 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. Author summary Nanoparticles (NP) have been used in various drug delivery contexts because they can target specific locations in the body. However, there is a translational gap between animals and humans, so researchers have begun toward computational models to guide in vivo NP experimentation. Here, we present several versions of physics-based multiscale physiologically based pharmacokinetic models (PBPK) for determining NP biodistribution. We successfully developed two versions of ODE-based compartmental models (models A and B) and an ODE-based branched vascular model implemented in MATLAB and Julia and validated models with experimental data. Additionally, we demonstrated using a neural network to solve our ODE system. In the future, this model can integrate different NP surface chemistries, immune system effects, multiscale vascular hydrodynamic effects, which will enhance the ability of this model to guide a variety of in vivo experiments.

The multiscale model must incorporate NP hydrodynamic properties in the vasculature, 48 March 5, 2022 2/34 NP-Endothelial Cell (EC) adhesion properties, and NP subcellular interactions that 49 govern targeted uptake to fully describe the movement and accumulation of NP within 50 the body [6]. In order to create a comprehensive multiscale model, NP behavior must 51 be understood at the system, hydrodynamic, and cell adhesion scales.

52
A previously published multiscale PBPK model has determined binding constants of 53 intracellular adhesion molecule 1 (ICAM1) coated NPs to endothelial cell (EC) surface 54 receptors in mice and humans by utilizing the biophysical properties of the antibody to 55 receptor interactions, and the cell surface [16]. Additionally, this model determines the 56 percent injected dose (%IDG) that distributes to the tissue in given organs at 57 steady-state; non-specific uptake is not accounted for in this model (NP uptake via 58 passive diffusion in the intercellular cleft). The binding constants determined through 59 the previous model [16] model will help to drive the binding characteristics of NP in the 60 multiscale model described in this study. Additionally, the rate at which NP bound to 61 the EC layer are endocytosed into a specific organ tissue must be considered. The 62 concentration of NP retained within the tissue or biodistribution is ultimately most 63 important since this allows for understanding how effectively NP can target tissue. 64 The purpose of this study is to 1) advance existing steady-state multiscale PBPK 65 models [16] to incorporate NP uptake via nonspecific transport, 2) develop novel 66 multiscale PBPK compartmental models to predict temporal effects, and 3) introduce a 67 compartmental branched vascular model that can predict the effect of hydrodynamic 68 interactions that depend on NP size, flow, and vasculature network properties, 4) 69 perform validation with experimental murine biodistribution data, 5) propose efficient 70 solvers for coupled stiff systems that embody the above properties, as well as make the 71 solvers compatible with contemporary machine-learning-based modules (such as neural 72 networks) which can capture and incorporate multiphysics models in the PBPK 73 workflow. 74 Materials and methods 75 Overview 76 Many Monte Carlo-based models have been proposed in previous works and have been 77 integral in understanding multivalent receptor-NP interactions at a molecular level. 78 Agrawal and Radhakrishnan [19] quantitatively characterized NP-endothelial cell 79 interactions by determining the multivalence of NP binding as well as antigen clustering, 80 ultimately providing future models with details about the energetics of the NP binding 81 process. Ramakrishnan [16] used these binding properties to understand the role that 82 expression levels of NP receptors play when targeting live cells, validating with 83 experimental data. While these adhesion-centric models are necessary for understanding 84 NP characteristics at the molecular level, they do not always include the interaction of 85 NP with the vascular network and translate into the pharmacodynamic scales. By 86 considering hydrodynamic parameters of the vascular system, and cell-binding/uptake 87 parameters into a variety of organ tissues (as determined in previous studies such 88 as [16]), a model that is governed by multiple time scales can be created, providing a 89 more physiologically relevant determination of NP biodistribution with temporal 90 resolution. This paper describes three variations of a multiscale pharmacokinetic model: 91 1) a steady-state model that facilitates the translation of the multivalent free energy of 92 adhesion into a biodistribution; 2) two different compartmental models (models A and continuous temporal biodistribution of ICAM-1 targeted NP for a murine model in each 98 case. The multiscale nature of the model described here can allow for customization 99 with various NP sizes, shapes, and uptake parameters, resulting in a customizable 100 predictive platform for different NP chemistries, organisms, and pathophysiologies. 101

102
A steady-state, multiscale, PBPK model was developed by Ramakrishnan [16] to 103 determine the percent injected dose of NP per gram of tissue (%idg) in five organ 104 compartments: lung, heart, kidney, liver, and spleen. However, this model only 105 considers NP uptake via antibody-receptor mediated multivalent interactions, failing to 106 account for receptor-mediated internalization and non-specific transport through pores 107 between the intercellular cleft of adjacent endothelial cells, especially in clearance 108 organs. Here, the existing steady-state PBPK model was modified to incorporate this  The original model [16] is described using Eq (1), where κ p is the non-specific binding of NP, K EC is the association constant for binding 118 of NP to endothelial cell surface receptors, K M is the association constant for NP 119 binding to a macrophage cell, D EC and D M represent the diameter of endothelial cells 120 and macrophage cells respectively, φ M and φ EC represent the concentration of 121 endothelial cells and macrophage cells in the target tissue respectively, L EC,b and L M,b 122 represents the distance from an EC or macrophage surface receptor that the NP can 123 successfully bind respectively, C out represents the injected concentration of NP, and 124 L cap represents the size of the cell-free layer in the capillary in which the NP is perfused. 125 Incorporating the non-specific transport of NP into the tissue, the modified model 126 equation can be described using Eq (2): consists of all organs that are not explicitly included in the model. We have developed 142 two ways in which the organ compartments can be connected by the arteries and veins, 143 shown in Figure 1. Figure 1a, is the model A configuration (an oversimplified version of 144 the circulatory system with 5 organ compartments), while Figure 1b shows the model B 145 configuration model, which is more physiologically relevant (lung circulation has been 146 separated out to keep track of oxygenation and oxygen distribution, and gut and spleen 147 compartments are coupled to the liver compartment, and gut/other compartments were 148 incorporated as well).

149
Each organ compartment in either model configuration consists of three additional 150 compartments: vascular, endothelial cell, and tissue compartments shown in Figure 2.

151
As NPs enter an organ compartment through the vascular compartment, they can either 152 bind to the ICAM-1 endothelial cell surface receptors (K on ), enter the organ tissue 153 compartment via non-specific uptake (K N S ), or be degraded (K deg  The model A configuration can be described using a system of 17 linear ordinary 163 differential equations (ODEs). Model A is described using equations 3-7, where i = 164 lung, heart, kidneys, liver, spleen. The venous compartment is described by: 166 Eq (3) describes the change in concentration of NP in the venous compartment over 167 time, where Q i and Q vein represent the flow of blood through organ compartments and 168 the veins, respectively, V vein denotes the volume of blood in the vein, C vein is the 169 concentration of NP in the veins, and K vein deg is the degradation rate of NP in the veins. 170 The arterial compartment is described by: 172 Eq (4) describes the change in concentration of NP in the arterial compartment over 173 time, where Q art represents the flow of blood through the arteries, V art denotes the 174 volume of blood in the arteries, C art is the concentration of NP in the arteries, and 175 K art deg is the degradation rate of NP in the arteries. For each tissue type i, the vascular 176 compartment is described by: Eq (5)   receptors, and K i N S denotes the rate of nonspecific NP uptake into the organ tissue.

187
The endothelial compartment in each tissue is described by: 189 Eq (6) describes the change in concentration of NP bound to the endothelial cell 190 surface receptors over time, where K i up denotes the rate of uptake of NP into the organ 191 tissue via transcytosis. Each tissue compartment is described by: 193 Eq (7) describes the change in concentration of NP in the organ tissue over time ODEs, however, due to the addition of two organs and the alternative configuration of 199 the model, several equations differ, while the overall structure is the same. Equations liver, spleen, and gut compartments (Q hep = Q liver + Q spleen + Q gut ). C bl denotes the 206 concentration of NP in the vascular compartment of each respective organ (lung, heart, 207 kidneys, liver, spleen, gut, other) March 5, 2022 7/34 Eq (10) describes the change in concentration of NP in the vascular compartment of 214 the lung over time.
Eq (11) describes the change in concentration of NP in the vascular compartment of 216 the heart over time.
Eq (12) Eq (13) describes the change in concentration of the NP in the vascular 221 compartment of the liver over time. Eq (6) and Eq (7) were used in the original model 222 to describe the change in concentration of NP in the endothelial cell compartments and 223 tissue compartments respectively, can also be used in the modified model configuration. 224 Additionally, Eq (14) was used to express the total NP degraded in the system at a 225 given time t, which will be useful to determine mass conservation of the system.

M ol
where t is the total time the model is run, g is the total number of organs in the 227 model, and ∆t is the t step of the model.

228
Eq (15) was used to determine the mass conservation of the system. If the system is 229 closed, the line formed by M ol t total over every time point t of the simulation should have 230 a slope of 0.
Eq (16) describes the mass conservation in molar basis, this equation was used to 232 determine the mass conservation when the system of ODE was set up in terms of molar 233 profiles.

234
Total i ∈ {V asculature, Endothelial, T issue, V ein, Artery, Degraded}. The development of this branched network was modeled after the branching network 252 described in [20], but is a simplified version of their three-dimensional vascular 253 branching network model. The diameter of the daughter vessels is governed by the 254 power law relationship, where the parent vessel is d 0 and daughter vessels are d 1 and d 2 , 255 and k=3, which is typically assumed based on a minimum dissipation principle 256 representing a stable flow condition: Due to the branching nature of the model, we can determine the number of segments 259 (N ) at any generation of branching (n) using the relationship N = 2 n . While the symmetric bifurcation (d i,1 = d i,2 ), the following relationship can be used to relate the 264 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 266 ratio β = 3, and the diameter of the vessel, d.
Equation (21) can be used to determine the total vascular volume of a branching 279 element,

281
It is important to note that the number of generations, volume, and surface area of 282 one branching element is fixed. I.e., the volume of a single branching element does not 283 differ across organs. Each organ will have a different number of branching elements 284 dependent on ϕ i , which is the ratio of total organ vascular volume to that of one 285 branching element.
Eq (22) Eq (23) 319 Equation (24) describes the NP concentration in the lung vasculature in generation 320 n = 17 to n = 31 of the branching network.

321
V bl,32 lung dC bl,32 Equation (25) describes the NP concentration in the lung vasculature in generation 322 n = 32 of the branching network. The heart compartment vascular branching network 323 can be described similarly to the lung compartment, with a series of four equations.
March 5, 2022 11/34 Equation (26) describes the NP concentration in the heart vasculature in generation 325 n = 1 of the branching network. The equations describing the NP concentration in the 326 heart vasculature in generations n = 2 to n = 16 and n = 17 to n = 31 are the same as 327 equations (23) and (24) respectively, but i = lung should be replaced with i = heart. 328 V bl,32 heart dC bl, 32 Equation (27)  = Q spleen C bl,32 spleen + Q gut C bl,32 Equation (28) describes the NP concentration in the liver vasculature in generation 333 n = 1 of the branching network. The equations describing the NP concentration in the 334 liver vasculature in generations n = 2 to n = 16 and n = 17 to n = 31 are the same as 335 equations (23) and (24) Equation (29) describes the NP concentration in the liver vasculature in generation 337 n = 32 of the branching network. The equation describing the NP concentration in the 338 kidneys, spleen, gut, and other compartments are again constructed similarly to the 339 branching equations in the lung and heart compartments, and again consist of a series 340 of four equations.
March 5, 2022 12/34 Equation (31) describes the NP concentration in the liver vasculature in generation 347 n = 32 of the branching network. Next, we can describe the NP concentration that is 348 bound to the endothelial cell layer.
Equation (32) describes the NP concentration bound to the endothelial cells of i = 350 lung, heart, kidney, liver, spleen, gut and other in generation n =1 to n = 32. Finally 351 we can describe the concentration of NP that is distributed to the organ tissue.
Equation (33)  The compartmental model is parameterized with a variety of physiological inputs such 359 as blood, tissue, endothelial cell volumes, and blood flow rates collected from previously 360 published sources ( [21], [22]). Other parameters such as K on and K of f for ICAM 361 antibody-coated NP were computationally determined in a previous study [16], while 362 yet other rate constants were determined via local sensitivity analysis through 363 comparison to existing translational studies [23].

364
The flow rates through each individual organ, Q i , were determined using Table 3 (in 365 supporting information) from [21], which gave blood flow rates through mouse organ 366 vasculature in units of L/g/min, which were ultimately converted into units of L/min 367 for use in the compartmental models by using mouse organ weight data from [22] and 368 the female and male masses were averaged. Additionally, to calculate the total flow, Q, 369 the sum of flow rates across all organs was taken.

370
Table 1 (in supporting information) in [21] listed values for vascular volume in 371 various mouse tissues. The volumes were listed in units of L/g. The mouse organ weight 372 data from [22] was used to determine the tissue volume for the entire organ, V i bl .

373
The volume of blood in the mouse vein was computed using the following equation, 374 V total = V art + V vein + V i bl ,where V total = 2146 µL. The volume of blood in the veins 375 (V vein ) and the arteries (V art ) is considered to be the same. Then the volume of blood 376 in the mouse venous system was computed.

377
Values for interstitial tissue volume and extracellular tissue volume were obtained 378 using Table 2 (in supporting information) from [21]. Values for interstitial tissue volume 379 (given in L/g) were used to determine the volume of each organ in its entirety, V i T . By 380 using the organ weights from [21], as was done to determine flow rates (Q i ), the tissue 381 volume for each organ in the mouse was computed with units of L.

382
The parameter l N C,b denotes the height of the endothelial cell layer, is constant 383 between organs and species and has already been defined by a previous model [16]. A i 384 can be calculated by using known relationships from [16], eg, It is known that . Given ϕ= 0.3, D EC = 5 × 10 −6 , and the V T for each 386 mouse organ previously described, l 2 EC can be calculated. When A i and l N C,b are 387 multiplied together, the resulting value will be representative of the volume of the 388 endothelial cell layer; V i EC .

389
The binding rate of NP to the endothelial cell surface (K i on ) was determined using 390 the relationship D l 2 given in [16]. The unbinding rate of NP from the endothelial cell 391 surface (K i of f ) can be determined by using K i on and K EC from [16] for each organ).

392
However, the log of the K EC had to be taken to accommodate for crowding effects [24]. 393 So, K i of f is calculated as a ratio of K i on to K i EC ; i.e., we assume a diffusion limited The relation between K EC and particle diameter was obtained from Figure 9 of [24]. A 398 linear relationship between log(K eq )) and particle diameter (a) was established using 399 the slope-point form for every organ, the point on the line being the (a, log(K EC )), for 400 a = 800nm, which was obtained using Monte Carlo simulation in [16]. The binding rate 401 is computed as; K on = D l 2 where l = thickness of the glycocalyx layer and D is given by 402 the Stokes-Einstein equation; D = K B T 6πηr , then K of f = Kon log(K EC ) . The values of K on ,

403
K of f , and log(K EC ), are reported in Table 1, 2, and 3, in supporting information.  Larger K i deg resulted in a smaller maximum biodistribution and a quicker decay in the 413 biodistribution curve. The results of an example local sensitivity analysis for the spleen 414 is shown in Figure 4. The parameters determined via the local sensitivity analysis that 415 were used in the compartmental models can be found at the SI GitHub link. While a 416 local sensitivity analysis was performed on certain parameters, it is important to note 417 that no extensive linear optimization or global sensitivity analysis was done on these 418 parameters; this will be done in future iterations of this work. link. It is logical to assume that non-specific uptake should be greatest in the liver and 422 gut because they are considered clearance organs. However, the local sensitivity analysis 423 revealed that the highest non-specific uptake rate was in the lungs. This discrepancy is 424 due to the presence of potentially invalid experimental data. It was apparent that no 425 matter how high the liver K N S value was set during the sensitivity analysis, the 426 biodistribution produced by the model was always significantly smaller than the 427 experimental data set. In fact, forcing a match to the liver data can only be realized at 428 the expense of violating the conservation of mass, which indicates an experimental error 429 in the reporting of the liver data. It is important to note that earlier iterations of the 430 model not described in the paper used other experimental data sets for validation, and 431 the liver biodistribution output matched the experimental data well. So, we conclude 432 that the experimental data reported by [23] for the liver is invalid. As for the gut and 433 other K N S values, they were entirely arbitrary. Since the experimental data set used 434 for validation did not report NP biodistribution data for gut or other, we could not 435 perform a valid sensitivity analysis on this parameter for these compartments. To ensure that the model output can be considered accurate, the model was validated 452 using data collected from an experimental study [23]. This study used a variety of sizes 453 of PEGylated gold nanoparticles and sampled the biodistribution at various time points 454 (5, 30, 60, 120 minutes) in five locations in a mouse. Using WebPlotDigitizer  [23]) and converted to %ID/g values using the relationship 457 presented in Figure 6 (b) of [23]. Then, using given organ weights from [21], %ID/t 458 (percent of injected dose in whole tissue) was calculated. These %ID/t values 459 determined experimentally in [23] were then compared to the model output to ensure 460 the compartmental model is indeed predictive. Single time points for the liver and 461 spleen were reported reliably, so these data sets were used for model validation. The stiffness ratio characterises the stiffness of the system. For a system of linear 477 ordinary differential equations d y dt = Ay, The stiffness ratio is defined as the ratio of the 478 largest and smallest eigenvalue of the matrix A. Quasi Steady State Approximation (QSSA) 485 The system of ODEs described by our PBPK model is a stiff system. The stiffness in organs. This reduced system consists of ODEs for concentration profiles in vein, artery 500 and organ tissue (9 in total). The system is much less stiff than the complete systems 501 because K on and K of f not do occur in the reduced system of ODEs. Functions f and g 502 consists of linear combinations of terms N ss and N nss for different organs, because the 503 original system of ODEs is a linear system. Then the QSSA involves setting: It has been shown in [25] that Neural Networks can be used to solve partial differential 510 equations; we use the same protocol to solve our system of ODEs using a neural 511 network. Our system being a stiff one makes it ideal for testing the performance of a 512 new method and comparing it against highly optimized ODE solvers. It has been 513 shown [26] that neural networks can be used to solve ODEs when the system of ODEs is 514 nonstiff. We used QSSA to make our PBPK model nonstiff (stiffness ratio ∼ 10 4 ) and 515 trained a neural network to learn the solution to the ODEs.
The derivative of the output from the neural network can be approximated by a finite 520 difference scheme or can be obtained using autograd functionality of a neural network. 521 Here we use a first-order finite difference scheme. 522 The loss function for this problem is: The first term in the loss function is the squared error between LHS and RHS of the 524 ODEs, computed using forward pass through the Neural Networks and the second term 525 is the squared error between true initial condition and predicted initial condition. The 526 solution NN(t, θ * ) is a unique solution to the system of ODEs, where, Eq (41) is the optimization problem aiming to minimize the above defined loss in terms 528 of neural network parameters θ.

529
The procedure described in Algorithm 1 aims to solve Eq (35) and (36) iteratively. 530 Notations used in the algorithm mean the same as described in these equations.  Table 5 in supporting information. 546 This suggests that the incorporation of non-specific uptake into the model results in a 547 much more physiologically relevant model. We also compare the results obtained from a 548 30-minute simulation of an unbranched model with experimental data (reported at a 549 30-minute time-stamp). We observe that by performing local sensitivity analysis for 550 area and volume of the endothelial layer and organ tissue respectively, we obtain a high 551 R for lungs compared to all other models. The excellent agreement of the steady-state 552 data with the 30-minute time data from the transient model suggests that the 553 steady-state model is a good approximation for time scales.   [23]. 568 The ODE solvers in Julia's DifferentialEquations.jl package does not require 569 user-defined timestep and exact mass conservation was obtained by using any of the 570 available stiff solvers. MATLAB's ode15s solver requires a user-defined timestep (∆t). 571 Timestep was chosen so the system remains stable, and so the system exhibits mass March 5, 2022 19/34 seconds, the system became unstable and mass was not conserved. On the other hand, 575 if ∆t was decreased, the model was unable to complete running due to the computing 576 constraints of MATLAB. It is important to note that mass conservation is only accurate 577 to the order of ∆t but integration is valid to a higher order. So, to run the model for a 578 time scale similar to that of experimental studies while maintaining stability and 579 conservation in the system, a ∆t of 0.001 for models A and B is necessary. The development of the branched model is incredibly important because it provides a 583 more accurate representation of the circulatory system; specifically, the surface area to 584 volume ratio of the blood vessels in the branched model is more representative of in vivo 585 mouse models. The framework of the branched model allows for the incorporation of 586 more specific hydrodynamic interactions and margination allowing for a physics-based 587 prediction of NP biodistribution, enabling us to account for characteristics such as NP 588 size, shape, and surface chemistry in the model, following the theory in [27].

589
Additionally, the construction of the branched model will allow for K on to change 590 depending on vessel diameter and blood flow rate. This cannot be done in models A 591 and B without empirical measurements.

593
Namely, the stiff solvers, QNDF, Rodas4, KenCarp4, TRBDF2 and RadauIIA5 were 594 able to solve for t = [0, 10 4 ]s. Whereas the most efficient stiff solver in MATLAB 595 (ode15s) was unable to solve the system of ODEs. Figure 8 shows the comparison 596 between the molar profiles for the branched and unbranched models.

598
The branching model construction allows for exploration of the effect of differing NP 599 sizes on biodistribution. The effect of nanoparticle size on biodistribution was also 600 explored using the branched model and nanoparticle diameters of 4nm, 15nm, 50nm, 601 79nm, and 100nm. In Fig 9, it can be seen that the 5-minute simulation result from our 602 branched model is in good agreement for 5-minute experimental data from [23] for the 603 kidney. Even though we were not able to get an exact match for Liver and Spleen, we 604 observed similar trends between experimental data and simulation data.  solving for the complete model, but the objective is to reduce the system's stiffness.

627
After using QSSA the system of ODE was solved using Tsit5 solver in Julia, which is a 628 nonstiff solver. Figure 11 shows the comparison between QSSA and the complete 629 solution of the branched model. The difference in initial onset can be explained from 630 Figure S1 in supporting information, where we can see the gradient is zero for most of 631 the time but there is a spike initially. QSSA for the branched model is a better 632 approximation compared to QSSA in the unbranched model because of the more rapidly 633 vanishing gradients in the former's case.  The original steady-state model from [16] was modified by adding a nonspecific 677 uptake term that represents NP uptake via passive diffusion through the intracellular 678 cleft. The addition of nonspecific uptake increased the predictive ability of the model, 679 as evidenced by higher R values when compared to an experimental data set than the 680 original model. This suggests that incorporating a nonspecific uptake term into future 681 models is necessary to increase physiological relevance and predictive ability. Next shortening the length of the plateau, leading to a quicker decrease in NP concentration 713 in a given compartment.

714
The multiscale PBPK model framework presented in this paper presents a significant 715 advance because the predictive ability is purely mechanism-based. Multiscale 716 physics-based modeling allows for the system's behaviours and interactions to be 717 completely described mathematically, rather than relying on empirical observations and 718 data to make predictions. Typical PBPK models ( [8], [10]) are generally empirically 719 based and do not describe the entire behaviour of the system. Creating a purely 720 physics-based PBPK model allows for more customization. For example, in this case, it 721 allows NP composition to be varied by changing certain model parameters to reflect 722 differing NP surface chemistry or size. This is advantageous since NP exist in many  The branching model results in an incredibly large and stiff system of ODEs. To 733 attempt to combat these issues, the stiff MATLAB ode solver, ode15s, was used. because a small time-step needed to be used to ensure stability within the model. If the 736 time-step was increased beyond one nanosecond, the system became unstable, and if the 737 run time was increased beyond one millisecond with the one nanosecond time-step, 738 MATLAB became unresponsive. So, it is clear that there are computing power 739 limitations in MATLAB solving large stiff ODE systems. We used stiff equation solvers 740 from DifferentialEquations.jl package in Julia to address this issue. All the stiff solvers 741 from the package successfully solved the system large times.

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

757
In this study, we successfully advanced existing steady-state multiscale PBPK 758 models [16] to incorporate NP uptake via nonspecific transport, as evidenced by an 759 ability to better predict in vivo biodistribution data (determined by R). Specifically, the 760 total R values across all organs are much higher in the modified model presented in this 761 paper than the original steady-state model.

762
Secondly, we developed novel multiscale PBPK compartmental models to predict 763 temporal effects of NP distribution. We developed Model A (5 organ compartments, 764 simple configuration) and Model B (7 organ compartments, more physiologically 765 relevant configuration) and solved with both MATLAB ode15s solver and a variety of 766 stiff solvers available in Julia's DifferentialEquations.jl package, for 10,000 seconds (2.78 767 h). to determine NP biodistribution temporally, and model outputs were compared to 768 in vivo experimental data. According to NRMSD, Model B more accurately represents 769 the experimental data set used for comparison [23]. 770 Thirdly, we introduced a compartmental branched vascular model that can predict 771 the effect of hydrodynamic interactions that depend on NP size, flow, and vasculature 772 network properties. This branched vascular model is the most physiologically relevant of 773 the models described, as it provides the most accurate representation of a circulatory 774 system. The branched model framework allows for incorporating specific hydrodynamic 775 interactions and marination, which allows a more physics-based prediction of 776 biodistribution, which allows us to consider more specific NP characteristics, like NP 777 size. We explored the effect of NP size (4nm, 15nm, 50nm, 79nm, 100nm) on 778 biodistribution and determined that a five-minute model simulation agrees well with 779 five-minute experimental data. 780 Finally, we proposed efficient solvers for the stiff systems that make up the models 781 described above and make the solvers compatible with machine learning modules like 782 neural networks that can incorporate multiphysics models into PBPK workflow. The vasculature at large times. The branched model was also reduced from 457 equations to 787 9 equations using QSSA. Again, the outputs of the reduced and full branched models 788 are similar, and differences in initial onset can be explained by the gradient being zero 789 for most of the time except for an initial spike in the vasculature and endothelial. We 790 concluded that QSSA for the branched model is better than QSSA in the unbranched 791 model due to more rapidly vanishing gradients in the former. Finally, QSSA was coupled 792 to a neural network and determined that the output of the neural network was the same 793 as that of the ODE solver, demonstrating the ability of a neural network to solve ODEs. 794 This model can continue to evolve to allow for customization beyond NP size by 795 adding modules in future model iterations. For example, incorporating a module that 796 describes internalization rates of varying NP surface chemistries to describe 797 biodistribution of other NP beyond those that are ICAM coated. Additionally, 798 incorporating immune system effects and more specific hydrodynamic interactions into 799 the model would increase NP's physiological relevance and translational potential to be 800 utilized more frequently in the clinic for targeted therapies.  function of NP diameter. Table 2: K of f values for different organs as a function of NP 806 diameter. Table 3: log(K EC ) values for different organs as a function of NP diameter. 807 Table 4: K deg , K U P , K N S , rates in each organ compartment.  The MATLAB and Julia code (which includes data used to run the models) for this 811 study can be found in the Multiscale PBPK Nanoparticle Biodistrubtion GitHub