Numerical modeling of nanodrug distribution in tumors with heterogeneous vasculature

The distribution and accumulation of nanoparticle dosage in a tumor are important in evaluating the effectiveness of cancer treatment. The cell survival rate can quantify the therapeutic effect, and the survival rates after multiple treatments are helpful to evaluate the efficacy of a chemotherapy plan. We developed a mathematical tumor model based on the governing equations describing the fluid flow and particle transport to investigate the drug transportation in a tumor and computed the resulting cumulative concentrations. The cell survival rate was calculated based on the cumulative concentration. The model was applied to a subcutaneous tumor with heterogeneous vascular distributions. Various sized dextrans and doxorubicin were respectively chosen as the nanodrug carrier and the traditional chemotherapeutic agent for comparison. The results showed that: 1) the largest nanoparticle drug in the current simulations yielded the highest cumulative concentration in the well vascular region, but second lowest in the surrounding normal tissues, which implies it has the best therapeutic effect to tumor and at the same time little harmful to normal tissue; 2) on the contrary, molecular chemotherapeutic agent produced the second lowest cumulative concentration in the well vascular tumor region, but highest in the surrounding normal tissue; 3) all drugs have very small cumulative concentrations in the tumor necrotic region, where drug transport is solely through diffusion. This might mean that it is hard to kill tumor stem cells hiding in it. The current model indicated that the effectiveness of the anti-tumor drug delivery was determined by the interplay of the vascular density and nanoparticle size, which governs the drug transport properties. The use of nanoparticles as anti-tumor drug carriers is generally a better choice than molecular chemotherapeutic agent because of its high treatment efficiency on tumor cells and less damage to normal tissues.


Introduction
Nanodrug carriers are advantageous over conventional molecular medicine in cancer therapy due to their higher tumor selectivity [1]. The therapeutic efficiency of anti-cancer drugs is highly correlated with their spatial and temporal concentration distributions in the tumor, PLOS  which are governed by the tumor environment [2] and the physicochemical properties of a drug. The uniformity of the drug concentration distribution affects the therapeutic effect on the entire tumor, and the cumulative concentration dominates the survival rate of cells. Therefore, the aim of drug delivery is to achieve a high and uniform distribution of the cumulative drug concentration in a tumor. Since drug delivery relies on the vascular system, an abnormal vasculature affects the deposition of drug molecules in a tumor through blood vessels. The presence of the high interstitial pressure in the tumor also hinders the drug delivery [3,4]. The drug molecules are extravasated from blood vessels, and their transport in the interstitium is driven by diffusion and convection effects. Diffusion effect is caused by the concentration difference in the interstitium, while the convection effect is driven by the interstitial pressure gradient.
The concentration difference in the interstitium is mainly the result of the heterogeneous vascular distribution in the tumor [5]. Tumor blood vessels are highly irregular in their structure compared with those in normal tissues. Unlike normal vessels, tumor vessels are dilated and tortuous, and their vascular walls are leaky and more permeable than normal vessels [6][7][8]. Moreover, the vascular distribution of tumor is highly heterogeneous. Tumor angiogenesis starts from the outer region and then spreads into the inner region. The proliferation of tumor cells results in a well-vascularized region in the periphery and a less vascularized region near the tumor center, in which a necrotic core may form, as illustrated by Fig 1. The heterogeneity of the blood vessel network leads to a non-uniformly cumulative concentration distribution of the drug within the tumor. In the tumor, the interstitial pressure is high and the interstitial pressure gradient is near zero due to a less functional lymphatic network. The function of a lymphatic network is to drain excess fluid from tissues to maintain the interstitial fluid balance and to prevent the occurrence of high pressure. However, functional lymphatic vessels can only be found in the tumor periphery, and the lymphatic vessels together with blood vessels at the center of a tumor are compressed by cancer cells and therefore often collapsed [9,10]. As mentioned in the previous paragraph, the tumor vessel walls are leaky and thus fluid can easily leak from blood vessels to tumor tissues. The less functional lymphatic system in a tumor gives rise to the insufficient drainage of fluid, thereby leading to the fluid accumulation in the interstitium and a high interstitial pressure around the center of tumor tissues. On the other hand, the vasculature at the outer region of a tumor can drain the excessive fluid; therefore, the interstitial pressure drops quickly [9,10]. The pressure gradient at the periphery region induces an outward convection, which pushes drug particles away from the tumor.
The tumor cell survival rate can serve as an indicator to evaluate the therapeutic effect and to estimate the probability of tumor recurrence. Putten and Lelieveld reported that there existed a log-linear relationship between the tumor survival rate and the extracellular drug concentration [11]. However, El-Kareh and Secomb also argued that the cell survival rate was more closely related to the intracellular drug concentration than the extracellular concentration [12]. Nevertheless, both works suggested that the tumor survival rate was proportional to the drug concentration whether it is intracellular or extracellular. High drug concentrations generally lead to a lower survival rate as more cells are killed by the drug. In this work, we developed a mathematical tumor model starting from the governing equation describing interstitial pressure distribution introduced by Soltani and Chen [13], and in addition we added the drug transport model to depict the temporal evolution and spatial distribution of drug concentration. The model was applied here to study a subcutaneous tumor with heterogeneous vascular networks considering the physiological characteristics of a tumor. Both nanoparticle carriers and a molecular therapeutic agent were investigated. Their difference in terms of the therapeutic effect and damage to the surrounding normal tissues were compared. The objective is to investigate the influence of tumor characteristics and drug properties on the drug AUC (area under the curve) distributions in the entire tumor by numerical simulation of the proposed mathematical model. In addition, the cell survival rate was estimated to quantify the therapeutic effect.
Since this work is primarily focused on drug delivery, we did not study the detailed tumor growth model, but described it by a simplified equation (Eq (23)). Because little is known about the regrowth of human tumors after chemotherapy [14], several methods have been proposed to model the kinetics of tumor repopulation [15][16][17][18][19]. One of the most popular models is to assume an exponential growth of tumor cells while accounting for a decreased growth rate as the volume of tumor increases [17]. The tumor growth is considered to follow the Gompertz model, which has become a widely accepted growth process for tumor growth in particular [18,20,21]. Consequently, we employed the Gompertz model to simulate the tumor repopulation after chemotherapy.

Materials and methods
The anti-cancer drug concentration in a tumor is affected by the vascular density, the size of a tumor, and characteristics of a drug such as molecular weight, diffusivity, vascular permeability, uptake rate by tumor cells, and the plasma half-life. Doxorubicin and various sized dextrans were considered and they were assumed to be well circulated in the body. The drug concentration in the plasma will gradually diminish due to the drainage by the lymphatic system and tissue absorption, as well as the clearance by the body.

Geometric configuration
In our model, tumors were assumed to be spherical and symmetrical. Unlike previous studies in Baxter and Jain [22], the vascular distribution of the tumor was modeled to be spatially heterogeneous and the vascular system comprised arterial capillaries and lymphatic vessels.
The value of the current one-dimensional (1D) model relies on its consistent and detailed governing equations for drug transport. In spite of simplification in geometry compared with two, three-dimensional (2D/3D) models [23][24][25], a well-constructed 1D model still can explore important physical mechanisms of drug transportation in tumors, and offer significant tumor treatment advices. A good 3D model sure can do better in prediction, but it takes long computing time and huge amount of resources to run a case, as compared to 1D models. Especially, numerous numerical simulations are required for testing and fitting the model parameters and analyzing their sensitivity afterwards, which will be an exhaustive task for 3D models, but they are much easier to be implemented in 1D models. 1D models do have their limitations compared with 2D/3D models. Usually a good tumor mathematical model is hierarchical with dimensions, and 1D model can provide a good test of mathematical model and identify important parameters efficiently before extending to 2D/3D.
Both an isolated and a subcutaneous tumors can be investigated by this model. The subcutaneous tumor is surrounded by normal tissues while an isolated tumor is not. A tumor region can be divided into a vascularized and necrotic regions. The vascularized region is mostly in the peripheral region of a tumor while the necrotic region comprises the tumor core. Though the proposed model is versatile, only the case of a subcutaneous tumor is presented in the current work.
One of the key factors needs to be considered in tumor modeling is the vasculature in tumor. The vasculature of tumors, featuring heterogeneous distribution of vessel sizes and shapes in 3D computations, is simplified here by the vascular density distribution S/V, i.e., the surface area of blood vessel (S) per unit volume of tumor (V). The vascular density is assumed to be zero in the necrotic core and then increase radially from the necrotic core boundary, reaching highest at the tumor boundary since the vascularized region is mostly in the peripheral region. A typical vascular density distribution in a tumor with a necrotic core is illustrated in Fig 2 in which the radius of the necrotic core is set to R n = 0.4R [8], where R n and R denote the radius of the necrotic core and the tumor, respectively. For normal tissues, vascular density distribution was assumed to be homogeneous and therefore a constant value for S/V, since the vasculature in normal tissues is generally well structured. As vasculature consists of arterial and lymphatic(venous) capillaries, both the arterial and lymphatic vascular density distributions of the tumor, S/V and S L /V, are assumed to be zero in necrotic core and otherwise a sine function here, as shown in Fig 2. Venous vessels have similar functions as lymphatic vessels; therefore, the presence of venous vessels will be omitted in the paragraphs that follow. A major difference between the current work and Baxter and Jain [22] is that they considered S/V distribution to be uniform in both tumor and normal tissues under different constant values, and their vasculature only consisted of arterial capillaries without considering the drainage by venous/lymphatic capillaries.

Interstitial pressure
To reach targeted tumor cells, systemically delivered drug needs to first extravasate from blood vessels, passes through the extracellular matrix, and finally penetrates into intracellular sites.
The extravascular transport of drug carriers relies on convection and diffusion, which are determined by the velocity of the interstitial fluid and the drug concentration distribution. The fluid velocity is proportional to the interstitial pressure gradient as described by Darcy's law, where k (cm 2 /mmHg-s) is the hydraulic conductivity of the interstitium and P i (mmHg) is the interstitial pressure. In the biological tissues with the fluid source and fluid sink in the medium, the continuity equation for the steady-state incompressible flow should be modified as r Áũ ¼ 0 B þ 0 L , where ϕ B and ϕ L (s -1 ) are the fluid source from blood vessels and the sink by lymphatic drainage, respectively. The fluid source term denotes the flow flux out of the vascular wall per unit volume, and is governed by Starling's law as follows: where ϕ B is the volumetric flow rate out of blood vessels per unit volume of the tumor, L p is the hydraulic conductivity of the vascular wall (cm/mmHg-s); the vascular density S/V is the blood vessel surface area per unit volume of the tumor (cm -1 ); σ s is the average osmotic reflection coefficient; P B and P i are pressure in blood vessels and interstitium (mmHg); π B and π i are the osmotic pressure of the plasma and interstitial fluid (mmHg). The lymphatic drainage term is described similar to blood vessel, but without osmotic pressure difference: where −ϕ L is the volumetric flow rate into lymphatic vessels per unit volume of the tumor, L pL is the hydraulic conductivity of lymphatic wall (cm/mmHg-s), the lymphatic vascular density S L /V is the lymphatic vessel surface area per unit volume of tumor (cm -1 ), and P L is the pressure in lymphatic vessels (mmHg). Combining Darcy's Law and the continuity equation with the substitution of ϕ B and ϕ L can result in Eq (4) can be written as [22] and R is the tumor radius. The dimensionless parameter α measures the ratio of the vascular conductance to the interstitial conductance. P ss is the steady-state pressure, a weighted average of P e and P L at which the flow rate outflux from blood vessels to the interstitium equals the influx from interstitium into lymphatic vessels. The effective pressure is defined as P e = P B − σ s (π B − π i ). Note that it seems P ss in the tumor will be a function of r, since S/V and S L /V are generally functions of r, as depicted in Fig 2. However, with S/V and S L /V both being zero in the necrotic core and having sine distribution in other region of tumor, P ss will be reduced to a piecewise constant distribution in tumor. Together with S/V and S L /V being constants in normal tissues, P ss can be more specifically expressed as, Because of symmetry, a pole condition with the pressure gradient (hence the interstitial fluid flow velocity) being zero is implemented at the center of tumor: For a subcutaneous tumor, the boundary conditions at the tumor-normal tissue interface, based on the continuity of pressure and velocity, are as follows, and where k t and k n are the hydraulic conductivity of the interstitium in tumor and normal tissues. Another requested boundary condition for pressure in normal tissue is that meaning the interstitial pressure would approach its P ss in normal tissue when far from the tumor edge. All the equations above describe the hydrodynamics of interstitial fluid in a subcutaneous tumor. Eq (5) will be solved to obtain the interstitial pressure distribution P i . After that, the interstitial fluid flow velocityũ i , flow source ϕ B , and flow sink ϕ L would then be computed, which will be used in the drug transportation model described below. A typical P i and P ss distributions in a subcutaneous tumor would look like Fig 3A while the accompaniedũ i , and ϕ B , ϕ L are shown in Fig 3B and 3C, respectively. The figure shows that the interstitial pressure P i remains constant with the value being tumor's P ss2 within most of the tumor, and has a steep descent across the interface of tumor and normal tissues, and then decays to another constant value of P ss in normal tissues. This causes large outward interstitial fluid flow around the interface, as shown in Fig 3B, and the accompanied drug convection would tend to push drug away from tumor and towards normal tissues. It might degrade the therapeutic effect and be more harmful to normal tissues at the same time. One might wonder why most part of the tumor, especially in the non-vascular necrotic core, experiences high interstitial pressure P ss2 . This is understandable that the fluid inside the necrotic core is trapped there and is forced to have the same pressure as P ss2 . This is mainly because all the vascularized region is mostly in the periphery of tumor and this pressure source forms a large back pressure at the core of tumor due to lack of lymphatic vessels to release the pressure. Because S/V and S L /V in tumor tissues share the same shape of sine distribution here, P ss yields a piecewise constant distribution, as shown in Eq (8). Under this circumstance, the effects of radial distributions of S/V and S L /V on the interstitial pressure happen to be exactly the same as S/V and S L /V specified to be constant distributions throughout the tumor like in [13]. Although having the same interstitial pressure distribution, hence the same interstitial fluid velocity distribution, these two different distributions of S/V and S L /V in [13] and the current model still will have different influence on drug transport through Eqs (14) and (15).

Drug transport
The flux of drug in the interstitium can be described asÑ ¼ À DrC i þũC i , where C i is the interstitial drug concentration, D is the diffusion coefficient (cm 2 /sec), and −DrC i andũC i denote the diffusion and convection terms, respectively. Diffusion is driven by the concentration gradient and convection is governed by the fluid velocity. Besides these two processes, the drug source from blood vessels, drug drainage by lymphatic vessels and the decay of drug concentration in the interstitium due to uptake of drug by tumor and normal cells immersed in interstitium are also taken into consideration. As a result, the drug delivery equation is expressed as where φ R = R d C i represents the drug uptake with R d (1/s) denoting the uptake rate or sink coefficient; φ SB (mole/cm 3 -s) and φ SL (mole/cm 3 -s) are respectively the solute effluxes from blood vessels and lymphatic vessels per unit volume, which are described by Kedem-Katchalsky equation as where ω is microvascular permeability (cm/s), ΔC = C h − C l (mole/cm 3 ) and " Þ are the concentration difference and average of solutions placed at both sides of the blood vessel membrane, respectively. Here the high concentration C h will be the drug concentration in blood vessels C B and the low concentration C l will be C i . Likewise, the rate of solute transport across the lymphatic vessel φ SL can also be expressed as where ω L is permeability for lymphatic vessels (cm/s), ΔC L (mole/cm 3 ) and " C L (mole/cm 3 ) are the concentration difference and average of solutions placed at both sides of the lymphatic vessel membrane, respectively. The expression of ΔC L and " C L are same as ΔC and " C s above, but with C l = C i and C h being the drug concentration in lymphatic vessels C L . σ L is the average osmotic reflection coefficient for lymphatic vessels and is set to be the same as σ s here. Similarly, the permeability of lymphatic vessels ω L was set to be the same vascular permeability ω in the tissue. As to the boundary/interface conditions for Eq (13), at the center of the tumor, there is a pole condition @C i @r At the interface of tumor and normal tissue, the concentration and its flux needs to be continuous over there: and À D @C i @r where the latter can be further simplified to by Eqs (10) and (17).
As r ! 1 in normal tissues, the boundary condition over there is a convective boundary condition, which is simply Eq (13) with the diffusion term neglected over there.
To close the equations, we still need governing equations for C B and C L , which are governed by an exponential decay, characterized by the half-life decaying rate λ (1/s) and related source/ sink terms from Eqs (14) and (15) based on mass conservation as and where λ is related to the plasma half-life time τ as λ = ln 2/τ. Several things should be noted for Eqs (20) and (21). First, they are "local" compartment models. C B and C L are temporal-spatial variables like C i . They are consistent with the spatial distributions of ϕ B , ϕ L , S/V and S L /V, as shown in Eqs (14) and (15). However, due to circulation constraint to vessels, unlike Eq (13) there will be no diffusion and convection in them. This reduces Eqs (20) and (21) to ODE models. It means C B and C L can be computed locally without coupling to its neighborhood explicitly. Second, C B and C L have the same plasma half-life time λ here. It is because the drug decays in circulation system (consisted of both blood and lymphatic vessels) through the uptake of liver. The initial conditions for C i , C B and C L are specified as in both tumor and normal tissues and The equality of C B and C L initially is due to drug goes directly from blood vessels to lymphatic vessels at the instant of drug injection (not yet drained to the interstitium). The vanishing of C B and C L in the necrotic region is due to the lack of vasculature. Eqs (13), (20) and (21) will be solved for C i (r, t), C B (t) and C L (t) in both tumor and normal tissues with boundary/interface and initial conditions described above. Note that since C i (r, t), C B (t) and C L (t), are proportional to C max , judging from Eqs (13), (20) and (21). We can scale these concentrations by C max , which is equivalent to set C max = 1 all the time.

Parameter values
The material parameters of the tumor and normal tissues used in this work adopted the values published by Baxter and Jain [22], and these values are listed in Table 1. Note that S/V in tumor and normal tissues in [22] were in uniform distribution with constant value of 200 cm -1 and 70 cm -1 respectively. Since in our current model the S/V and S L /V are heterogeneously distributed in tumor, we set the maximum values of S/V to be 200 cm -1 [2,5] and maximum values of S L /V to be 20 cm -1 in tumor tissues, and both to be uniformly distributed in normal tissues with a constant value of 70 cm -1 . As no value of hydraulic conductivity for tumors has been reported in the literature [22], L pL value can only be deduced from the reported lymphatic filtration coefficient L pL S L /V and the estimated S L /V value [26]. S/V and S L /V have the same constant value in normal tissue, but different maximum values in tumor because normal tissue, compared with tumor, has a much more extensively functional lymphatic network, which removes the net fluid filtered from the blood vasculature. Thus the extravasated materials are more quickly re-absorbed by lymphatic vessels in normal tissues [22,27]. Also, unlisted in Table 1, the pressure in lymphatic vessels P L is 5 mmHg smaller than the effective pressure P e according to [22]. Dextrans with the molecular weight of 10 kDa, 70 kDa, and 2 MDa were chosen as the nanoparticle carriers and the characteristics of dextrans in a tumor are shown in Table 2 [28]. Nugent and Jain (1984)  The vascular permeability of nanoparticles in the tumor was estimated to be about 7.8 times that of normal tissues, as measured by Gerlowski and Jain [31]. The counterpart for normal tissues are shown in Table 3. The parameter sink in Tables 2 and 3 is R d in Eq (13), corresponding to the drug uptake in tissues. For comparison with the nanoparticle dextran, doxorubicin was chosen as the conventional molecular anti-cancer drug. The molecular weight of doxorubicin is 544 Da and Table 4 lists its transport parameters utilized in this work. Note that the vascular permeability ω generally has larger value in tumor than in normal tissue. For the reason that drug leaks to both tumor and normal tissues through blood vessels, but blood vessels in tumor are usually much more leaky than those in normal tissues due to the irregularly developed vasculature during angiogenesis.

Tumor cell survival
The cumulative drug concentration can be computed by use of the equations above, which was subsequently used to estimate the cell survival rate of tumor and normal tissues after treatments. The therapeutic effect and the probability of tumor recurrence can be measured based on the tumor cell survival rate. In 1983, Greene et al. reported that the survival rate in a tumor was an exponential function of the extracellular concentration of anti-cancer drug [36]. In 2000, El-Kareh and Secomb argued that it would be more reliable to estimate the cell survival rate with the intracellular concentration in cells [12]. In our model, the cell survival rate (S F ) was defined as S F = N/N 0 , where N 0 and N denote the numbers of cells before and after a treatment. The relationship between the cell survival fraction and the extracellular concentration was defined as: where k s is the dose-survival constant and AUC, abbreviation for area under the curve, is a frequently used pharmacokinetic term denoting the cumulative concentration in the interstitium over time; in other words, the area under the concentration-time curve, which can be described as The dose-survival constant (k s ) employed the value given in Jusko's work [37], which is 4.329 × 10 −3 (1/nM-hr). The cell survival rate after each single treatment can be estimated by Numerical modeling of nanodrug distribution in tumors with heterogeneous vasculature use of Eq (24). In the field of pharmacokinetics, the area under the curve (AUC) is the definite integral (from zero to infinity) in a plot of drug concentration in tissue vs. time. The AUC represents the total drug exposure over time.
When estimating the cell survival rate after multiple treatments, the regrowth of tumor cells after each treatment needs to be considered. The tumor cell growth can be estimated by the three-parameter Gompertzian function [18,19,38]: where N i is the number of viable tumor cells after i-th treatment; n i (t) is the number of tumor cells at time t after i-th treatment; N 1 is the saturated cell number that can be reached after a long period of time and b is related to the initial tumor growth rate [38]. Here we adopted the initial number of tumor cells as N 0 % 5 × 10 9 , N 1 = 3.1 × 10 12 and b = 0.0283 month -1 from the work of Yorke et al. [38]. The cell number for the 1-cm thick normal tissue was assumed N 1 = 4.64 × 10 12 here. We assumed that the number of cells is proportional to the tumor volume, which means that the tumor size shrinks (R decreases) as the number of cells decreases, and assumed that the ratio of normal tissue density to tumor cell density to be 0.2 [39].
Another assumption is that the shape of vascular density distribution in the tumor remains unchanged after each treatment. Moreover, the regrowth of normal tissue cells was also evaluated by use of Eq (26) under the assumptions that the mean growth rate b of normal tissues is half of that in a tumor. The whole simulation procedures are summarized and depicted by the flowchart contained in Fig 4. The flow chart describes the repeated procedures in multiple treatments consisted of the following stages: (1) drug delivery phase: computing C(r, t) and AUC(r), (2) tumor cell killing phase: computing cell survival rate S F and tumor radius after treatment, (3) tumor regrowth phase: computing tumor radius after regrowth.

Numerical method
Here we employ method of lines (MOL) with the semi-discretization in space at first for Eq (13) by multi-block Chebyshev pseudospectral method. Together with interface/boundary conditions, it then would result in a set of ordinary differential algebraic equations (ODAEs) in time. We can then solve this ODAE set by a highly efficient variable-step-variable-order (VSVO) ODAE solver like ode15s in Matlab. Due to the high order of accuracy intrinsic to pseudospectral method, far less grid points are needed to achieve quite accurate numerical results compared with traditional finite difference/element methods. Also, since ode15s has an error estimate at each step, it can automatically adjust its time step to avoid instability and gain the optimal efficiency in time integration.

Results and discussions
Our model can compute the spatial distribution of the interstitial pressure and the cumulative drug concentration featured by AUC. Furthermore, the spatially-averaged AUCs resulted from different anti-cancer agents can be calculated and hence the therapeutic effect of nanoparticle carriers and molecular chemotherapeutic agents can be evaluated via the cell survival rate. A 2-cm (diameter) subcutaneous tumor with a 0.8-cm necrotic region, i.e., R = 1 cm and R n = 0.4R, embedded in a 1 cm thick normal tissues was used in our current study.  blood nor lymphatic vessels are present and the interstitial pressure is flat there as well, the convection and source/sink terms vanish in Eq (13). The only way for drug to penetrate the necrotic region is through diffusion, which is generally weak there. Though most tumor cells are already deceased in the necrotic region, it does not mean the amount of the drug penetration is insignificant in treating a tumor since there are still some surviving tumor stem cells hiding inside. In the viable tumor region, the drug level increases fast in the beginning due to double sources from the permeation terms of φ SB and φ SL in Eq (13). The first term in Eq (14) acts as the source while the first term in Eq (15) acts as sink due to ϕ B and ϕ L , as shown in Fig  3C. These two terms roughly cancel each other except around the interface of tumor and normal tissues. The second term in Eqs (14) and (15) is the major source term for drug transported from vessels (blood and lymphatic vessels) into the interstitium, and the driving force is the concentration difference between vessels and interstitium. This driving force is modulated by the vessel permeability (ω and ω L ) and vascular density (S/V and S L /V). Also the increase of C i is biased towards the interface due to convection, as shown in Fig 3B. Accompanying the increase of C i naturally comes the decay of both C B and C L . Here C B decays much faster than C L in the viable tumor region chiefly due to S/V >> S L /V.
In normal tissues, unlike the viable tumor region, the interstitial concentration C i increases very slowly in the beginning due to small values of φ SB and φ SL . This is a result of the much smaller vascular permeability (ω and ω L ) in normal tissues compared with its counterpart in tumors, as shown in Table 4. In addition, smaller vascular densities (S/V and S L /V) in normal tissue than those in tumor contribute to this slow increase as well. Drug penetrating into normal tissue through convection (driven by a higher interstitial pressure in a tumor) and diffusion (due to faster growth of drug concentration in the viable tumor region in the beginning) also causes the increase of drug concentration in normal tissue, but only limited to normal tissues near the interface. Generally, we wish AUC for normal tissues to be as small as possible to avoid harming normal tissues. As time goes, C i changes from increasing to decreasing and dies out eventually due to continuous weakening of driving force and the uptake of drug by cells in interstitium.
The concentration contours of 10 kDa, 70 kDa, and 2 MDa dextrans as well as that of doxorubicin are displayed in subfigures A-D of Fig 6. It can be seen from Fig 6A and 6D that concentrations of small particles like 10 kDa dextran and molecular drug doxorubicin have larger contour values, but decay quickly over time. On the contrary, larger nanoparticles like 70 kDa and 2 MDa dextrans have small contour values, but decay much more slowly. This is mainly controlled by vascular permeability. As shown in Tables 2 and 4, vascular permeability monotonically decreases with molecular weight. Generally speaking, a larger vascular permeability causes drug to leak more quickly from vessels into interstitium, result in a higher level of interstitial concentration, but at the same time decay more quickly. On the contrary, smaller vascular permeability causes drug to leak slowly from vessels into interstitium, result in smaller levels of concentration, but sustain longer before vanishing. Both vascular permeability and density are smaller in normal tissue than in tumor; consequently, less drug is able to leak from vessels into the interstitium. As a result, drug concentration levels are generally lower in normal tissue than in tumor. Note that the concentration contours of 2 MDa dextran exist almost strictly in the viable tumor region, which indicates almost no drug transport from tumor into normal tissue. This can be attributed to its small concentration difference between tumor and normal tissues owing to the smallest vascular permeability among all drugs. Besides, its smallest diffusion coefficients among all drugs also hinders the diffusion of drug from tumor to normal tissues. It naturally implies that 2 MDa dextran kills tumor cells most and at the same time harms normal tissue least. Actually, the therapeutic effect of drug should be evaluated by its AUC distribution, i.e., the time integral of the local concentration. Fig 6 shows the contours of drug concentration C i (r, t) computed from Eq (13), and its time integral from zero to infinity will generate the spatial distribution of AUC like Fig 7A. Therefore, both the drug concentration level and its duration in the interstitium are important contributing factors to AUC. Judging from 2 MDa and 70 kDa dextrans having larger spatial distributions of AUC in tumor tissues than 10 kDa dextran and doxorubicin, as shown in Fig 7A, it appears that the duration time affects the AUC more significantly than the momentary concentration levels here.
The spatial distributions of AUC are shown in Fig 7A and the average values in different regions are listed in Table 5. This gives us a more precise therapeutic effect prediction than Fig  6. All anti-cancer drugs shared a similar distribution pattern in AUC. Compared with the normal tissue region, their AUCs in the tumor vascular region were usually high, but extremely low in the necrotic core. This low AUC in the necrotic region is simply because there is no vascular vessel there to transport drugs. Therefore, very limited amount of drug is transported into a necrotic region solely via diffusion from a vascular region. Also, this low AUC in the necrotic region might imply insufficient dosage to kill the tumor stem cells hiding in it. Note Numerical modeling of nanodrug distribution in tumors with heterogeneous vasculature that smaller drugs penetrate deeper into normal tissues and affect more normal tissues while larger nanoparticles display a shallower penetration. This has been analyzed via Fig 6. The survival rate of normal tissue will be particularly affected by the degree of penetration.
As shown in Table 5, the average AUCs of larger-sized dextrans were much higher in the vascular region and lower in normal tissue than those of doxorubicin, which suggests that larger-sized dextrans generally have a better therapeutic effect and less side effect than doxorubicin. This shows that using larger-sized nanoparticles as drug carriers is more effective at treating a tumor and less harmful to the surrounding normal tissues than molecular chemotherapeutic agents. Based on average AUCs, 2 MDa dextran has the best therapeutic effect to the vascular tumor region and second least side effect to normal tissues, and this generally indicates that larger nanoparticles are the better choice for treating a well-vascularized tumor.
Besides calculating the average AUC of drugs, we estimated the cell survival rate based on Eq (24) and used it to quantify the therapeutic effect after treatments, as shown in Fig 7B. Consistently, 2 MDa dextran has the best performance with lowest survival rate distribution in tumor and second highest survival rate distribution in normal tissue.
The tumor cell responses to all four kinds of drug after multiple treatments were also investigated and the spatially-averaged survival rates of viable tumor region and normal tissue are depicted in Fig 8. The survival rate drops after each treatment followed by regrowth based on Eq (26) Fig 8A, which was about 73.79% after one doxorubicin treatment. The overall tumor cell survival rate was about 36.15% after four doxorubicin treatments. In addition to the effectiveness of killing cancer cells, the side effects to normal tissues during chemotherapy also needs to be addressed. Therefore, the cell survival rates of normal tissue were also investigated. From Fig 8B, the Numerical modeling of nanodrug distribution in tumors with heterogeneous vasculature normal tissue survival rate after the first treatment of 10 kDa dextran delineated by the red curve was about 91.56%, and about 70.49% of normal tissue cells were left after the whole therapy program. At the same time a large fraction of normal tissue was killed by doxorubicin treatment (the survival rate after the first treatment was about 77.15%) and only 35.69% of normal tissue cells were left after the whole therapy program. It is worth noting that the normal tissue survival rate was almost the same as the tumor survival rate in doxorubicin treatment. In other words, doxorubicin kills cancer and normal cells without distinguishment. These results show that 10 kDa dextran, though less efficient at killing tumor cells, was less harmful to normal tissues compared with doxorubicin. Fig 8C and 8D show the tumor and normal tissue survival rates after multiple treatments of 70 kDa and 2 MDa dextrans. For multiple treatments, larger-sized drugs like 70 kDa and 2 MDa dextrans are generally more effective in killing tumor cells compared with smaller-sized drugs like doxorubicin and 10 kDa dextran, and have a normal tissue survival rate between doxorubicin and 10 kDa dextran, which can be comprehended by AUC and S F in the normal tissue region shown in Fig 7. Effects of different treatment plans can be compared by the current model as well, as depicted by Fig 9. In the first therapy program, the maximum concentration (C max ) in the blood vessel was assumed to be 1000 nM and the total number of treatments was four times with 10 days interval between two treatments. In the second therapy program, the maximum concentration (C max ) in the blood vessel was cut into half to be 500 nM and 8 treatments were administered with 5 days interval between two treatments. The tumor and normal tissue cell survival rates resulted from these two programs were delineated in Fig 9. The result shows that the performance differences between 4 longer interval treatments and 8 shorter interval treatments are negligible, judging from both the tumor and normal tissue survival rates. Though the dose per treatment in the second plan is half of the first one and seems to be less toxic to the patient, it does not end up with a higher survival rate of normal tissue than the first treatment. The figure of the tumor and normal tissue cell survival rates after multiple treatments provides important prognostic information regarding the survival rate and percentage of the The overall survival rate after 10 kDa dextran treatments. The survival rates of tumor cells (blue curve) and normal cells (red dashed curve). A: The maximum concentration (C max ) in the blood vessel was assumed to be 1000 nM and four treatments were administered with a time interval between two treatments of 10 days. B: The maximum concentration (C max ) in the blood vessel was assumed to be 500 nM and the total number of treatments was eight times with a time interval between two treatments of 5 days. https://doi.org/10.1371/journal.pone.0189802.g009 Numerical modeling of nanodrug distribution in tumors with heterogeneous vasculature cell regrowth. It is of great help to evaluate the efficiency of a treatment plan. Note that in these multiple-treatment simulations, S/V and S L /V have been fixed values as tabulated in Table 1. However, angiogenesis after each treatment may actually be different and heterogeneity in the spatial distribution of blood and lymphatic vessels may occur, which means S/V and S L /V may vary after each treatment. Not just angiogenesis, even geometry of tumor would be different after each treatment. All these effects are not considered here, but will be studied in the future work by extending the current model to higher dimensions.

Conclusions
The proposed model was employed to compute the interstitial pressure, dynamics of drug transport, AUC distributions of drugs of various sizes, and tumor/normal tissue cell survival rates in a subcutaneous tumor. The model is useful to quantify the therapeutic effect. It has been found that AUCs of all drugs are generally high in the tumor vascular region, low in normal tissues, and extremely low in the necrotic region. Larger nanoparticles like 2 MDa and 70 kDa dextrans display high AUCs in the tumor vascular region and, on the contrary, smaller nanoparticle drugs like 10 kDa dextran and molecular anti-cancer agents like doxorubicin show low AUC there. AUC shows that the delivery of all kinds of anti-cancer agents in the necrotic region was insufficient to kill the tumor stem cells hiding in it. Therefore, other treatment methods such as surgical removal and thermal therapy may be used to enhance the effectiveness of cancer treatment. The distribution of cell survival rates demonstrated that the side effect to normal tissues by use of dextrans was limited to a small range while doxorubicin caused damage extensively. Overall, larger-sized nanoparticles were found to deliver better therapeutic effects to the tumor region with limited toxicity to the surrounding tissues, as compared with the molecular anti-cancer agents such as doxorubicin.
In addition to estimating the delivery of drugs in tumors, the current model can be of help in the treatment planning. By estimating the cell survival rate after each treatment and regrowth between treatments, our method can be used to compare treatment plans with different parameters like dosage and frequency of treatment cycle.

Future works
In the future, we plan to extend the present model to a more practical and complete one by the modification of the following three modules. (1) Geometry and hydrodynamics module: considering tumor with a geometry of higher dimensions; outlining explicit distributions of leaky blood and lymphatic vessels through angiogenesis into the model while still treating tumor interstitium as a porous media and modeling it by Darcy's law. (2) Drug transport module: modifying the current drug transport model for the tumor environment mentioned in (1).
(3) Time evolution of the spatial distribution of tumor cell density module: deriving a set of differential equations describing the rate of tumor cell density changes at least by its growth with nutrition supplied, death caused by the drug uptake in tumor, and natural death. By considering the above modified modules, simulations of this new model will shed more light on significant mechanisms of drug delivery inside tumor and, at the same time, offer a better and more practical evaluation of tumor treatment.