Enhanced Drug Delivery to Solid Tumors via Drug-Loaded Nanocarriers: An Image-Based Computational Framework

Objective Nano-sized drug delivery systems (NSDDSs) offer a promising therapeutic technology with sufficient biocompatibility, stability, and drug-loading rates towards efficient drug delivery to solid tumors. We aim to apply a multi-scale computational model for evaluating drug delivery to predict treatment efficacy. Methodology Three strategies for drug delivery, namely conventional chemotherapy (one-stage), as well as chemotherapy through two- and three-stage NSDDSs, were simulated and compared. A geometric model of the tumor and the capillary network was obtained by processing a real image. Subsequently, equations related to intravascular and interstitial flows as well as drug transport in tissue were solved by considering real conditions as well as details such as drug binding to cells and cellular uptake. Finally, the role of periodic treatments was investigated considering tumor recurrence between treatments. The impact of different parameters, nanoparticle (NP) size, binding affinity of drug, and the kinetics of release rate, were additionally investigated to determine their therapeutic efficacy. Results Using NPs considerably increases the fraction of killed cells (FKCs) inside the tumor compared to conventional chemotherapy. Tumoral FKCs for two-stage DDS with smaller NP size (20nm) is higher than that of larger NPs (100nm), in all investigate release rates. Slower and continuous release of the chemotherapeutic agents from NPs have better treatment outcomes in comparison with faster release rate. In three-stage DDS, for intermediate and higher binding affinities, it is desirable for the secondary particle to be released at a faster rate, and the drug with slower rate. In lower binding affinities, high release rates have better performance. Results also demonstrate that after 5 treatments with three-stage DDS, 99.6% of tumor cells (TCs) are killed, while two-stage DDS and conventional chemotherapy kill 95.6% and 88.5% of tumor cells in the same period, respectively. Conclusion The presented framework has the potential to enable decision making for new drugs via computational modeling of treatment responses and has the potential to aid oncologists with personalized treatment plans towards more optimal treatment outcomes.

Methodology: Three strategies for drug delivery, namely conventional chemotherapy (one-stage), as well as chemotherapy through two-and three-stage NSDDSs, were simulated and compared. A geometric model of the tumor and the capillary network was obtained by processing a real image. Subsequently, equations related to intravascular and interstitial flows as well as drug transport in tissue were solved by considering real conditions as well as details such as drug binding to cells and cellular uptake. Finally, the role of periodic treatments was investigated considering tumor recurrence between treatments. The impact of different parameters, nanoparticle (NP) size, binding affinity of drug, and the kinetics of release rate, were additionally investigated to determine their therapeutic efficacy.
Results: Using NPs considerably increases the fraction of killed cells (FKCs) inside the tumor compared to conventional chemotherapy. Tumoral FKCs for two-stage DDS with smaller NP size (20nm) is higher than that of larger NPs (100nm), in all investigate release rates. Slower and continuous release of the chemotherapeutic agents from NPs have better treatment outcomes in comparison with faster release rate. In three-stage DDS, for intermediate and higher binding affinities, it is desirable for the secondary particle to be released at a faster rate, and the drug with slower rate. In lower binding affinities, high release rates have better performance. Results also demonstrate that after 5 treatments with three-stage DDS, 99.6% of tumor cells (TCs) are killed, while two-stage DDS and conventional chemotherapy kill 95.6% and 88.5% of tumor cells in the same period, respectively.

INTRODUCTION
Encapsulation of drugs in nanoparticles (NPs) can enhance the delivery of therapeutic and diagnostic agents to tumors, and at the same time, reduce their accumulation in healthy tissue (1). In treatments, this issue results in improved local drug concentrations at disease sites while decreasing systemic toxicity. On the other hand, conventional small molecule compounds usually distribute randomly among all tissues (2). However, NPs using advantages such as enhanced permeability and retention effect (EPR) have provided moderate enhancements over small molecule therapeutics regarding patient survival and care (3). A qualitative analysis of the pharmacokinetics of NP-mediated doxorubicin (DOX) in patients with cancer was carried out by Gabizon et al. (4) comparing with an equal DOX dosage delivered in free form. The findings distinctly indicated that NP-mediated delivery is capable of decreasing plasma clearance and enhancing DOX concentration inside tumors (4). In vivo tests on human prostate carcinoma showed that NP-mediated DOX may improve treatment effectiveness because of decreased systemic elimination, enhanced penetration in tumor, and prolonged NP existence with a slow release rate of drug (5). The clinical advantages of NP-encapsulated DOX is also evaluated in patients with metastatic breast cancer (6). The results suggested that NPmediated DOX was capable of significantly reducing cardiotoxicity and improving therapeutic efficacy. Nevertheless, despite multiple promising ideas and experimental outcomes, NPs do not commonly reach tumors efficiently in clinical trials (7). Successful drug delivery depends on the ability of NPs to deep penetration into tissue, achieve an efficient spatial distribution, ensure its proper binding affinity, and adequately release the drugs (8). However, progress in NPs will need a number of technical and physiological barriers to be realized and overcome, such as opsonization and nonspecific protein adsorption, non-specific uptake by cells and organs encompassing the immune system, targeting and penetrating the tumor microenvironment (TME), and gaining access to cancer cells for intracellular drug delivery (9). Physical circumstances influencing the function of NPs facing these obstacles are still understood imperfectly, and the nanomedicine lacks a detailed explanation of physical principles to guide logical NP designs that can overcome biological barriers in TME (1).
For most desirable efficiency, sufficient amounts of the therapeutic agents must reach tumors in order to eliminate cancer cells, and at the same time, they must not induce considerable side effects on normal tissues. Generally, relatively small NPs experience higher transvascular and interstitial transports (10). Investigation of the various functions of NPs with different physiological features in combating biological barriers has demonstrated that NPs require to be dynamically adapted to these obstacles, resulting in the emergence and development of multi-stage DDSs (11). In such a system, the NPs can change their size in response to stimuli at different stages (e.g., initial size of 100 nm and secondary size of 10 nm) (12). Wong et al. (13) suggested a multi-stage system comprising the primary particle containing smaller secondary particles, which subsequently carry the therapeutic agent. Eventually, therapeutic agents contained in the secondary particles must be released to reach cancer cells. If an additional stage was included in the conventional NPs, it could improve drug distribution into the tumor and also tumor penetration, as well as increase treatment efficacy. The secondary particles release their cargo within the tumor, triggered by exposure to external (magnetic and electric fields, acoustic, etc.) or internal stimuli (TME properties such as enzymes, pH, etc.).
Since the introduction of the first model for chemotherapy, the use of mathematical and computational models has become widespread in examining different DDSs (14)(15)(16)(17)(18)(19). These models can provide guidance on required composition and preparation methods for administration of DDSs. Various computational models have been employed for simulation of NSDDSs to examine efficacy, understand biological phenomena, and select optimal treatment plans. Based on the investigated spatial and temporal scales, these models are classified to discrete, continuous, and hybrid models (8). A summary is presented in Table 1 of important studies conducted on employing mathematical modeling of drug-containing NPs for drug delivery to solid tumors .
As seen in Table 1, there exist a number of gaps in the literature regarding usage of drug-containing NPs for delivery of drug to solid tumors. With this motivation, in the present study, the main contribution is to apply a computational model of three DDSs, namely (i) a conventional chemotherapy system (onestage DDS) and (ii) a conventional two-stage DDS containing NP and chemotherapy (two-stage DDS), and (iii) a three-stage DDS containing a primary NP, a secondary NP, and the chemotherapy; on a real image of vascularized tumor. To this end, intravascular and interstitial fluid flows as well as drug transport equations are solved by considering actual conditions in tumor. Then, various parameters, NP size, binding affinity of drug ligands to the receptors of cells, and the kinetics of release rate, are studied to determine their effects on treatment efficacy. NPs employed in this study could for instance be drug-loaded nanocarriers, liposomes, or magnetic NPs. Additionally, the A PK/PD model Authors recommended that shorter injection duration might enhance treatment performance, and explored it by computational studies. The treatment outcome was predicated according to peak intracellular concentration over the whole treatment period. A comparison between bolus injection and continuous infusion demonstrated that duration of infusion had a great effect on the treatment outcome. Optimal duration is dependent on cellular pharmacokinetics. Drug release rate from non-TSLs is an effective parameter so that if this rate is optimized, the efficacy of non-TSLs is slightly fewer than continuous infusion.
-Spatial distribution is not considered.
Zhang et al.  Proposing a computational model for magnetically-assisted drug delivery approach to assess the penetration of drug into peritoneal tumors nodules and improve intraperitoneal (IP) chemotherapy.

A 2D-0D tumor model & FEM
A great enhancement in the intratumoral concentration of magnetic NPs compared to free drugs. The success of magnetic drug targeting in larger tumors (10-20 mm in size) is found to be significantly due to the strength of magnetic field and tumor-magnet distance while these two parameters are less important in small tumors.
-Avascular model; -Lack of capillary network and assumption of uniform vascular density in the tissue.
-Concentration equations are solved ignoring binding affinity and cellular internalization. Stylianopoulos et al (33)./2018 Reengineering the TME to enhance the efficacy of drug delivery from computational modeling to bench to bedside.

─ (Review paper)
Authors discussed the mechanics of both solid and fluid components of tumor, focusing on how they prevent the delivery of drug and create an abnormal TME that promotes tumor growth and resistance to treatment. They also provide strategies to re-engineer the TME by normalizing the vessels of tumor and the ECM to enhance the treatment of cancer. Eventually, they summarized different mathematical approaches that have provided insights into the physical obstacles against efficient cancer treatment and suggested novel methods to overcome these impediments. Presenting a mathematical modeling for spatial-temporal distribution of chemotherapy drug in TSL-mediated DDSs A PK/PD model Authors demonstrated that complicated relationships between the related factors (various chemotherapy drugs, release rate constants, and heating duration) and the predicted treatment result, making it difficult to identify the best parameter set. a model-based optimization approach is presented to overcome this challenge. Optimization showed that the best result would be obtained with a low drug release rate at physiological temperature, combined with a moderate to high release rate at mild hyperthermia and 1 h heating post-injection.
-Spatial distribution is not considered; Rezaeian et al.  Developing a mathematical modeling to analyze nanomedicine distributions in solid tumors A PK model Authors quantified the effect of influencing parameters on the efficacy of tumor delivery, the magnitude of heterogeneous distribution, and the EPR effect. They also compared the spatial distributions of the NPs and the free drugs within tumors. The model predicted high degrees of distributional heterogeneity for both NPs and free drugs. They found that diffusion coefficient of NPs was the most efficient factor in decreasing the NPs distributional heterogeneity but it has moderate impact on the free drugs.
-Real image of tumor is not considered; -Spatial distribution is not considered; -Lack of capillary network and assumption of uniform vascular density in the tissue.

Dogra et al. (38)/2019
Overviewing different mathematical modeling about application of nanomedicine in cancer treatment.
─ (Review paper) Authors provided an overview on mathematical modeling works that have been applied towards a better insights of nano-bio interactions for enhancing the efficacy of drug delivery to tumor. This model allows for drug features (e.g., size and binding affinity) to be explicitly defined, thus facilitating investigation into the interaction between the changing TME and cytotoxic and NP drugs. They predict a heterogeneous distribution of NPs after delivery; that NPs need a ECM with high porosity to cause tumor regression; and that transcapillary fluid velocity is dependent on porosity of ECM, and implicitly on the drug size.
-Real image of tumor is not considered. Presenting in silico modeling of nanomedicine for cancer, across scales and transport obstacles.
Authors investigated latest outcomes in multi-scale modeling of NP transport obstacles, as well as existing software packages, with the goal of focusing the wider research community in building a common computational platform that able to overcome some of the current barriers facing effective design of NPs. ─ (Review paper) Investigation of various issues regarding the use of NPs as vehicles of anticancer drug delivery: specifically, administration into the circulation system, transvascular transport, distribution in the extracellular matrix, cellular internalization, and release of drug from NPs.
impact of successive treatment cycles is numerically examined considering tumor recurrence between two consecutive treatments for three investigated DDSs.

MATERIAL AND METHODS
In this section, a detailed description of mathematical models and their assessment are investigated. First, the mathematical equations are presented, including equations governing the interstitial fluid flow, solute transport and cellular uptake of chemotherapeutic agents, and NP delivery system. Then, parameters of model, relevant physical and transport properties, and NP-related calculations are also provided. Finally, numerical methods to assess the mathematical models, model parameters, boundary conditions (BCs), input images as geometry, computational domain, as well as solution strategy are presented.

Delivery Mechanisms
Two well-known delivery mechanisms, namely chemotherapeutic delivery and delivery through drug loaded NP, are considered in the present study. In the following, their mechanisms are described in detail. After intravenous (IV) injection, the chemotherapeutic drugs (here, DOX) are transported to the tumor site through the blood vessels. Subsequently, they pass via the tumor vessel wall and travel the remaining distance from the vessel wall to the cancer cells. In tumor ECM, free molecules of drug agents can bind to receptors of the cells, unbind or get internalized (28,46).
In general, only a small fraction of an injected drug reaches tumor tissue, the remaining being cleared from the body. One possible way to overcome this problem is to target drug delivery by encapsulating the anticancer drug in a nanocarrier (47). In such a system, drugs which encapsulated inside NPs are injected into the blood, and subsequently, NP characteristics cause them to be released in accordance with a pre-designed controllable procedure. NP formulations have benefits over commonly used chemotherapies since they can combine many diagnostic and therapeutic factors (48). They are also associated with significantly lower side effects, owing to their capacity for optional accumulation in tumorous tissue. To eradicate TCs, NPs have to first reach the tumor via the vascular system, then extravasate from the relatively leaky regions of the microvessels into the TME. Subsequently, NPs release their cargo in a controlled manner, while staying put at the ECM. Drug molecules diffuse into the tumor slightly more than into normal tissue and can bind to cancer cells and/or TME components eventually, becoming internalized by cells (28). A schematic of drug delivery mechanisms considered in the current study is shown in Figure 1.

Governing Equations
The mathematical models for delivery of drugs to solid tumors consist of equations as: i. Mass and momentum conservation for interstitial fluid flow; ii. Mass transport for the free, bound, and internalized drug; iii. Mass transport for nano-encapsulated drug; and iv. Treatment efficacy for intracellular drug concentration.
Detailed descriptions of each aspect of mathematical modeling are provided next.
The dynamics process in drug delivery includes binding and unbinding of drug agent ligands with receptors of the cells at the rates of K ON and K OFF , respectively; exchange of drug between capillary network and interstitium, and influx/efflux of drugs from interstitium to TCs, internalization of drugs to cellular space, and finally cell-killing by drug agents. The cell-killing rate is governed by a model according to the predicted intracellular concentration of anticancer drugs. In the case of two-stage and

Fluid Flow in Interstitium
First, the interstitial fluid flow equations are solved to provide the basic biomechanical environment for transport of drug. Darcy equation, which demonstrates the relationship between interstitial fluid velocity (IFV) and interstitial fluid pressure (IFP), is employed to describe interstitial fluid flow in a porous environment as follows (16): where v i and k are the IFV and the hydraulic conductivity.
Considering the presence of source/sink terms in biological tissues, the continuity equation is modified as (7): in which f B is the rate of fluid flow from the microvessels to the extracellular matrix (ECM) and and f L is the rate of fluid flow from ECM to lymph system, defined as (7): in which S/V demonstrates the surface area per unit volume of microvessels. L PL ( S V ) L and P L are the lymphatic filtration coefficient and lymphatic pressure, respectively.

Transport of Drug in the Interstitium
The comprehensive model for drug transport in the interstitial space includes: -Transport in ECM by diffusion and convection mechanisms, -Transport across microvessels by diffusion and convection mechanisms, and -Binding to cells and internalization.
The drug transport equation for biological tissue can be written as (18): in which Ф B is the drug transport rate through the microvessels into the interstitial space, and Ф L is the drug transport rate from interstitial space into the lymph system. Ф B is defined according to Patlak's model, as (18,49): Pe demonstrates the Peclet number, s f is the coefficient of filtration reflection, P represents the permeability coefficient of microvessels, and C p is the injected drug concentration.
The drug transport rate through lymphatic microvessels has been considered only in healthy tissue as (18,49): A bolus injection of chemotherapeutic agents, representing the initial vascular concentration of the drug, is modeled as (50): where C 0 and k d are the initial concentration and blood circulation decay, respectively.

Convection-Diffusion-Reaction Modeling of Drug Transport in the Interstitium
The drugs exist in different forms as: NPs in the interstitial space (C N ), free drugs in the interstitial space (C F ), bound drugs (C B ), and intracellular drugs (C INT ) (20,49). In tissues, with regards to the existence of mass flow source/sink, a system of equations, which is called CDR equations, is utilized to represent the process of drug delivery. The general block diagram of the model of current study considering NPs and chemotherapeutic drugs for two-stage DDS is shown in Figure 2.

Conventional Chemotherapy
The CDR equations for conventional chemotherapy are as follows (17,41): Drug internalized into the cell (11) in which C F demonstrates free drug concentration, C B bound drug concentration, C INT internalized drug concentration, and C rec is the concentration of cell-surface receptor. In this equations, v is the IFV, D is the diffusion coefficient, K ON is the association rate, K OFF is disassociation rate, K INT , represents the cellular internalization; and j demonstrates the tumor volume fraction available for the drugs.

Two-Stage Drug Delivery System
Systemically administered NPs, as demonstrated in Figures 1  and 2, are transported to tumor sites through the circulation system, undergo transvascular extravasation followed by distribution in the interstitial space, and are finally delivered to c a n c e r c e l l s . F or t r a n s p o r t of a n N P c o n t a i n in g chemotherapeutic agents (C N ), the system of equations is adjusted thus (51): (12) in which K rel , D N , and a are respectively the drug release rate from the carrier, its diffusion coefficient, and the number of chemotherapy molecules contained in the nanocarrier.

Three-Stage Drug Delivery System
Three-stage NP drug delivery system is an efficient option to overcome the barriers of two-stage DDS in ECM. The system of equations for three-stage DDS, as demonstrated in Figure 1C, is as following (32): (13) in which C N1 and C N2 are the primary and secondary NP concentrations, respectively. K rel1 is the rate constant for the release of the secondary NP from the primary one and K el2 the rate constant for the drug release from the secondary NP. a and b are the number of secondary NPs released by the primary and drug particles released by the secondary NP, respectively.

Fraction of Killed Cells
Despite its cardiotoxicity, the drug DOX is frequently employed in tumor treatment. DOX is a standard-of-care, DNA-damaging agent employed in the treatment of multiple tumors (e.g., bladder, breast, and lung cancers). Using the internalized drug concentration, the efficacy of drugs was calculated according to the empirical equation obtained for DOX (52). The FKCs parameter is defined, as (53): in which S F is the fraction of cells remaining after treatment and w is a fitting parameter specified for DOX based on the results of experiments (54).

Tumor-Cell Survival
The number of TCs after a period of time (n i ), which is obtained through Gompertz's model, depends on intervals between chemotherapy sessions (t) (55). Gompertz equation is a function of three parameters (as demonstrated in Eq.)15(): the number of TCs surviving after the i th treatment (N i ), the number of saturated cells after a very long period (N ∞ ), and eventually the rate of tumor progression (b).
The number of surviving TCs after each therapeutic phase will be examined as a criterion of treatment efficacy evaluation. This criterion is obtained by using the FKCs according to Eq. (14). The number of remaining TCs as a result of the difference between the number of cells after (N i ) and before (N 0 ) treatment will also be examined as a criterion for evaluating treatment efficacy. In the present model, S F is also defined, as (30): The initial numbers of TCs here, adopted from York et al. (55), are N 0 = 5 × 10 9 , N ∞ = 3.1 × 10 12 , and b = 0.0283 month -1 . The number of cells for healthy tissue was considered to be N 1 = 4.64 × 10 12 (55). The cell number is assumed to depend on tumor volume; i.e., when the cell number decreases, the tumor shrinks (R reduces). The ratio of the density of healthy tissue to the density of TC is considered to be 0.2 (56), and it is also assumed that the microvascular density distribution in the computational field does not change after each treatment. In addition, the regrowth of healthy tissue cells was examined by using Eq. (15) with the assumption that the growth rate (b) of healthy tissue is half that of the tumor (30).

Model Parameters
Interstitial and Drug Transport Parameters Tables 2 and 3 demonstrate both the interstitial and DOX drug transport parameters, respectively, including both tumor and healthy tissues. Table 4 represents the baseline state parameters for NSDDS for 20 nm particle size and 200 nm vessel-wall pore (VWP) size.

NP-Related Calculations
The hydraulic conductivity, vascular permeability, and reflection coefficient are calculated for the NPs by using the theory of particle transport through cylindrical pores (28,59): in which g is the fraction of the surface area of a porous vesselwall, r o is the pore radius, h is the viscosity of water at 310 K, and L is the vessel-wall thickness. D o represents a particle diffusion coefficient in a free solution at 310 K, given by the Stokes-Einstein relationship as in (53,59): in which K b , T, and r s are the Boltzmann constant, the temperature, and the diffusing particle radius, respectively. H and W are diffusive and convective hindrance factors, respectively, and related to hydrodynamic and electrostatic interactions. Neglecting electrostatic interactions (E=0), H and W are reduced to (53,59): where F is the partition coefficient and is defined as (53,60): (23) in which l is the ratio of drug particle size to the VWP size, as follows (53,59): K t and K s coefficients in Eqs. (21) and (22) are given by (53,59): Parameter values used in the Eqs. (17) to (25) are demonstrated in Table 5.

Model Geometry and Boundary Conditions
In this study, a tumor model (as shown in Figure 3A) is employed as an input geometry based on real image of a tumor with a capillary network surrounded by healthy tissue, extracted from Roudnicky et al. (60). About the circumstance of this input image, it should be mentioned that this type of tumor is inoculated in a mice by injecting human A431 squamous cell carcinoma (SCC) cells (60). This image was taken about 12 days after tumor inoculation. In fact, 2 days after inoculation, thrombospondin-2 (TSP2), an anti-angiogenic matricellular protein that inhibits tumor growth and angiogenesis, was injected for 10 days and subsequently this image was taken. After image-processing, a computational field is considered with the existence of a tumor in the middle of the domain as well as the parent vessels ( Figure 3B).
Two boundaries are assumed in the computational field: the inner boundary (between the tumor and the normal tissue), and the outer boundary (at the outer edges of computational field). For the inner boundary, the continuity BC is considered for IFV, IFP, as well as concentration and its flux. For the outer boundary, in which the IFP is constant, the Dirichlet BC is used for fluid flow, and the open boundary is employed for concentration (50). The input and output pressure amounts for the parent vessels are selected as: P Inlet,1 = 25 mmHg, P Inlet,2 = 25 mmHg, and P Outlet =10 mmHg, according to realistic physiological conditions reported in literature (50,61).

Solution Strategy
There exist two distinct phases for solving the present problem: steady-state and transient. Calculations related to blood flow in a model with a discretized capillary network in the computational field   After image-processing, geometries were meshed and analyzed utilizing the COMSOL Multiphysics software-version 5.5a. The coupled nonlinear set of the above-mentioned governing equations and also the BCs were assessed through FEM. A segregated approach is applied to solve the equations with the time-step of 0.1 [s] and relative tolerance of 0.001. A sixfold drop of residuals is chosen as the criterion for convergence. For solving drug delivery equations in vascularized tumors, a Core (TM) i24 CPU @ 3 GHz with 32 GB RAM system is used.

RESULTS AND DISCUSSION
Characteristics of tumors are determinant factors in transport of drug and final efficacy of treatment. Since chemotherapy drugs are carried by the circulatory system, the tumor vasculature features play an important act in delivery of drug. In the present study, a tumor model incorporating details of capillary network distribution is employed to evaluate the impact of heterogeneous distribution of capillary network on delivery of drug. This is followed by studies of transport of drug in tumors with microvascular density (MVD). First, the results of chemotherapeutic agent delivery for a case study, extracted from real image of tumor, are proposed, and different parameters are investigated in detail. Subsequently, results of computational modeling for delivery of NPs are presented.

Validation of the Results
In this study, validation is performed with the same governing equations as the literature (28) for the FKCs over time ( Figure 4). As is clear, there is a correspondence between the results of the present study and those in the literature so that by considering the real geometry and physics, the FKCs has a similar trend. However, considering similar conditions, its value has about 4% difference, due to differences in tumor geometry, computational domain, and structure of capillary network. One should note the 2D geometry of the tumor and capillary network investigated in this study, while Stylianopoulos et al. (28) utilized 1D for the capillary network.

Baseline Model Analysis
The spatial-temporal distribution of the non-dimensional concentrations of drug that is taken up by TCs in the tumor and surrounding normal tissue for three investigated DDSs ─singlestage (conventional chemotherapy), two-stage, and three-stage ─ are shown in Figures 5-7. These non-dimensional concentrations ( e C i ) were obtained at a given time by dividing that concentration at any point by the maximum concentration in the entire field. Total concentration is defined as the sum of various drug concentrations. In single-state DDS, free drug concentration reduces over time and the drug gradually turns into a bound drug. Then, the bound drug slowly enters the intracellular space and is consumed there. It is also obvious that drug concentration in the tumor is greater than that in the normal tissue. The reason for this is the greater extravasation rate from microvasculature in the tumor area and additionally the highly-dense MVD in this region. The highest value for total concentration takes place in the tumor zone and this value is several times higher than the concentration in the normal tissue. In two-stage DDS, the loads of NPs are released in the interstitium in a free drug form and the remaining process is similar to singlestate system. In three-stage DDS, the primary NPs release the second NPs, and the remaining process is similar to previous two-stage DDS. e C INT represents the internalization of drugs to the cellular space, determining the succeed of drug delivery in killing cancer cells (i.e., leading to higher FKCs). Another important factor for efficient drug delivery is uniform distribution of drug in tumor, expressing the penetration depth of drugs extravasated from microvascular network. Comparison of the distribution of different drug concentrations in three investigated systems demonstrates that three-stage DDS provides much more uniform drug distribution than the other two ones. After that, two-stage system has more efficient drug distribution compared to singlestage DDS. From Figures 5-7, it is also clear that in conventional chemotherapy, there exist a small concentration of drug in healthy tissue, leading to side effects; whereas, in both NSDDSs, the side effects are negligible.

Analysis of NP Drug Delivery System
Parameters of the TME that inhibit the delivery of NPs into the tumor include the size-dependency of the transport both across the tumor microvessel wall and then via the interstitial space of tumor. Transport across the tumor vessel wall is determined by the relative size of the particle compared to the VWP size. On the other hand, not only the size but also other parameters of drug, such as the drug release kinetics might play a crucial role in the outcome of the therapy. Therefore, one of the main goals of the present study is to determine under what circumstances twostage and three-stage DDSs can be beneficial, relative to one another or to conventional chemotherapy. In the following, the results are presented for two-and three-stage DDSs. Then, a parameter study is carried out to examine the effect of three important parameters ─size of NPs, binding affinity, release rate of drug─ on the real geometry of tumor.

Two-Stage Drug Delivery System
The effect of release rate of drug for two sizes of NPs are demonstrated in Figure 8. Overall, FKCs for NPs with the size of 20nm is higher than that of 100nm, in all investigate release rates of drug. The main reason is the greater blood half-life of smaller NPs compared to larger ones. Therefore, a 20nm NP is expected to circulate in the blood for a longer time in comparison with a 100nm particle, further improving efficacy of delivery systems. Moreover, it is demonstrated that slower and continuous release of the chemotherapeutic agents from NPs have better treatment results compared to faster release rate.
NP delivery systems may have other functions in addition to acting as carriers of drug. One of these functions is to control the release rate of therapeutic agents from the NPs (2). Release rate from drug-loaded nanocarriers directly determines the toxicity and anticancer activity of a DDS. Depending on various factors such as the NP formulation, fabrication method, surrounding environment, etc., release rate can vary across a wide range. For instance, stealth liposomes can provide sustainable release over

Three-Stage NP Drug Delivery System
In addition to enhancing the tumoral penetration depth, multistage NP delivery systems can provide further tunability in the spatial delivery control to solid tumors (8,28). Due to the physiological obstacles that a chemotherapeutic agent must encounter, a multi-stage DDS can improve treatment efficacy by altering its physical properties (including shape, size, flexibility, charge, and/or surface coating) to suit the transport across each obstacle (8). As interest in increasingly complicated DDSs grows, we face a corresponding challenge to set the model parameters. In this study, we had presented two common scenarios for multi-stage DDS: (i) a 100-nm particle, which released secondary 10-nm particles; and (ii) a 20-nm particle, which released secondary 5-nm particles. Based on Figure 9 for the three-stage DDSs, we assumed multiple different possible states for release rate constants (i.e., K el1 and K el2 ) of the secondary particle and the drug. The results enable a set of observations as to how release rate constants affect the efficacy of drugs. The four most predominant are: (i) overall, the second scenario (NP1 = 20nm and NP2 = 5nm) has higher treatment efficacy, almost in all investigated states, compared to the first scenario (NP1 = 100nm and NP2 = 10nm); (ii) in low binding rates, the high release rates have better performance; (iii) in moderate and high binding rates, the NP release must have high release rates and the drug release must have the lower release rates; (iv) the least treatment efficacy occurs when both release rates are slow.
A high binding affinity of drug agents to receptors of cancercell leads to a higher drug concentration in the intracellular space of tumor, and it simultaneously decreases the drug concentration in tumor tissue. In other words, with a higher binding affinity, more drugs are taken up by cancer cells, reducing the drug level in the interstitium. Binding affinity plays different roles in the delivery of smaller or larger particles because in smaller particles the diffusion mechanism is more dominant than convection, while in larger particles the convection via vessel-walls is more dominant in the transport across blood vessels. It should be mentioned that a very high binding affinity of the NPs leads to aggregation nearby the vessels that NPs are extravasated from. On the other hand, for transport from microvessels to tissue, NP efficacy depends on the kinetics of drug release from particles (28). Eventually, the released drug may also rapidly bind to the cells, causing heterogeneous and incomplete distribution within tumor (62). After simultaneous investigation of the impacts of NP size, binding affinity, and the release rate of drugs, it is predicted that NP penetration from tumor microvessels would be affected (29). Indeed, there is a competition between diffusion in the ECM (which depends on NP size) and binding affinity and/or the release rate of drugs. However, the one of the main conclusions is that all these three parameters should be considered together and studying their impact alone does not provide the right policy on optimal design of NSDDSs.

Treatment Evaluation
In this section, the treatment outcomes of three investigated drug delivery systems are calculated by considering 1-month drug-free breaks between treatments. As demonstrated in Figure 10, threestage system has better treatment outcome than two-stage and one-stage (i.e., conventional chemotherapy) systems, implying superiority based on higher rate of killing TCs and shorter treatment time, simultaneously. For the conventional chemotherapy, the respective efficacies of the first, second, and third stages of treatment in reducing the TCs are about 36.78%, 29.67%, and 24.41%. This rate is 20.24% for the fourth stage, 14.9% for the fifth stage, 5.9% for the sixth stage, 4.65% for the seventh stage, 3.42% for the sixth stage, 2.73% for the ninth stage, and 2.11% for the tenth stage. Tumor regrowth between the first and second stages of treatment, between the second and third stages, and up to the end of the process in the tenth stage are calculated to be 13.28%, 10.57%, 8.06%, 5.59%, 3.4%, 2.85%, 2.4%, 2.17%, and 2.03%, respectively. At the last stage of treatment, it is obvious that the regrowth of TCs is about 2.03%, while the treatment efficacy is about 2.11%, implying that for this size of tumor using conventional chemotherapy alone is not the best choice and need more cycles of treatments, so an adjuvant therapy would be suggested at this size of tumor. In general, after seven treatments, 5.54% of the initial tumor cells remain for conventional chemotherapy.
Results indicate that after 5 treatments with three-stage system, 99.6% of TCs are killed, while two-stage and one-stage system respectively kill 95.6% and 88.5% of TCs in the same period. It should be mentioned that the rate of killing TCs for the two-stage system reaches 98.9% after 7 cycles of treatment. This rate is also 94.46% for one-stage system after 10 treatment cycles. As a consequence, the treatment efficacies for both the NSDDSs have shown significant improvements compared to conventional chemotherapy. Results of the case study show that the survival rates of TCs after several treatment cycles offer significant prognostic insight about the survival rate and the cell regrowth percentage. Being able to evaluate the efficacy of a treatment scenario using multi-scale computational modeling is very helpful in clinical situations.

CONCLUSION
In the present work, a mathematical model of drug transport is proposed to predict delivery of chemotherapeutic drugs, either in their free form or encapsulated in NPs. Parameters of the model describing physiological and biological characteristics of tissue and drug are extracted from experimental data in the literature. The model has been applied to realistic tumor reconstructed from actual image of vascularized tumor to: (i) understand the transport steps of non-encapsulated and encapsulated drugs, and (ii) elucidate the impact of drug properties on drug delivery and treatment. The main findings are as follows: • Spatiotemporal drug distribution illustrates that the complexity of the capillary network is the main factor for non-uniform drug distribution in the tumor. • The FKCs of tumor for two-stage DDS with smaller size of NPs (20nm) is higher than that of larger ones (100nm), in all investigate release rates. Slower and continuous release of the chemotherapeutic agents from NPs have better treatment outcomes in comparison with faster release rate. • For three-stage DDSs, in intermediate and higher binding affinities, it is desirable for the secondary particle to be released with faster rate, and the drug with slower rate. In lower binding affinities, the high release rates have better performance • Three-stage system has better treatment results relative to two-stage and one-stage systems, reaching 99.6% effectiveness in killing TCs after 5 treatments; while two-stage and onestage system respectively kill 95.6% and 88.5% of TCs in the same period.
Overall, a mathematical framework has been developed for drug delivery in both free and encapsulated forms to provide qualitative and mechanistic understanding of transport of drug in solid tumors. The effect of treatment outcomes considering tumor recurrence between the two presented models is highly complementary to PK/PD models as it considers both the spatial and temporal scales at the same time, while PK/PD models merely consider the temporal scales. Moreover, different equation parameters in the present study have physiological and bio-chemical implications (e.g., intravascular pressure, microvascular density, microvascular diameter, permeability, diffusion coefficient to name a few), whereas PK/PD methods might involve a parameter to express several parameters.
It should be mentioned that anti-angiogenic agent (here, TSP2) has effects on real tumor geometry and changes of permeability, but we have no discussion about TSP2 in the present study. On the other hand, if we have three images for three different states (negative control, treatment with TSP2 and without treatment), the mathematical modeling can include the new geometry of microvascular networks. This is an interesting suggestion for future studies. Additionally, this method can be applied on each images of tumor because the geometry of tumor can be obtained after image-processing; however, in this study, we have chosen an actual image from the literature due to the lack of experimental set-up in our group.
Protection systems for NPs such as PEGylation as well as adjustable parameters of NP design, all try to enhance drug halflife in the circulation system, which can increase drug delivery into the tumor and reduce side effects to healthy tissue. Furthermore, NPs can be manipulated to release their loaded drugs if exposed to a specific external (e.g., ultrasound or magnetic field) or internal (e.g., pH or enzyme) stimulus. Multi stimuli-responsive DDSs ─combination of internal and external stimuli─ have not only succeeded in targeted drug delivery but also in the multi-modal cancer diagnosis and treatment. In addition to enhancing efficacy of encapsulation, these systems may enhance the drugs half-life. Our presented approach has the potential of modeling these targeting systems. For example, by simultaneously solving bio-heat and wave equations with the presented equations, the drug delivery through thermo-sensitive NPs can be modeled.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article. Further inquiries can be directed to the corresponding author.