Effect of tumor shape, size, and tissue transport properties on drug delivery to solid tumors

Abstract

Background

The computational methods provide condition for investigation related to the process of drug delivery, such as convection and diffusion of drug in extracellular matrices, drug extravasation from microvessels or to lymphatic vessels. The information of this process clarifies the mechanisms of drug delivery from the injection site to absorption by a solid tumor. In this study, an advanced numerical method is used to solve fluid flow and solute transport equations simultaneously to investigate the effect of tumor shape and size on drug delivery to solid tumor.

Methods

The advanced mathematical model used in our previous work is further developed by adding solute transport equation to the governing equations. After applying appropriate boundary and initial conditions on tumor and surrounding tissue geometry, the element-based finite volume method is used for solving governing equations of drug delivery in solid tumor. Also, the effects of size and shape of tumor and some of tissue transport parameters such as effective pressure and hydraulic conductivity on interstitial fluid flow and drug delivery are investigated.

Results

Sensitivity analysis shows that drug delivery in prolate shape is significantly better than other tumor shapes. Considering size effect, increasing tumor size decreases drug concentration in interstitial fluid. This study shows that dependency of drug concentration in interstitial fluid to osmotic and intravascular pressure is negligible.

Conclusions

This study shows that among diffusion and convection mechanisms of drug transport, diffusion is dominant in most different tumor shapes and sizes. In tumors in which the convection has considerable effect, the drug concentration is larger than that of other tumors at the same time post injection.

Background

Cancer, the main cause of death in developed countries, is the second leading cause of death in developing countries [1]. Solid tumors account for 85% of human cancers [2]. Chemotherapy is one of the ways widely used for cancer treatment. Based on the findings from clinical applications, most cancer treatments with drugs fail to eliminate solid tumors completely [3]. The computational method can investigate why systemic administration cannot distribute drug uniformly in tumors. The drug exchange between microvessels and extracellular matrices, drug removal by lymphatic system, drug diffusion and convective transport in extracellular matrices should be included by mathematical simulation. Computational fluid dynamics (CFD) can model the whole drug delivery process and clarify the mechanisms of drug delivery from the injection site to absorption by a solid tumor.

Baxter and Jain, based on the theoretical framework in their 1D mathematical method, found the effective factors on drug delivery such as microvessel permeability, interstitial fluid pressure (IFP), and interstitial fluid velocity (IFV) [4–7].

Improving the 1D model of Baxter et al. [5, 6, 8] and Saltzman et al. [9], Wang et al. [10–12] developed a simulation framework of drug delivery to tumors by considering the complex 3D geometry. Wang and Li [10] used modified MRI images for tumor geometry. They considered interstitial fluid flow with blood and lymphatic drainage in their model. Wang et al. [11] studied the effect of elevated interstitial pressure, convective flux, and blood drainage on the delivery of specified solute to brain tumors.

The study of tissue transport property effect on drug delivery is considered in recent studies. Zhao et al. [13] used a 3D computational model to predict the distribution of IFV, IFP, and solute transport through a tumor. Arifin et al. [14] studied the sensitivity of drug distribution to physiochemical properties in realistic models of brain tumors. A specific tumor captured by MRI is used by Pishko et al. [15] for modeling drug distribution in tissue with spatially-varying porosity and vascular permeability. The sensitivity of solute distribution to tumor shape and size is not considered in above mentioned works.

In our previous work [2], tumor shape and size effect on drug delivery is investigated by modeling interstitial fluid flow and assuming that drug particles flow with the interstitial fluid. In the present work, by adding the solute transport equation to the previous developed model in our group [16–20], new governing equations are investigated to find drug concentration in interstitial fluid (DCIF). Solving fluid flow and solute transport equations simultaneously, the effects of tumor shape, size, and tissue transport properties on drug delivery to solid tumor are also investigated.

Spherical and non-spherical tumors and their surrounding normal tissue are modeled with assumption of rigid porous media. The vasculature as a source term and lymphatic vessel as a sink term vary spatially. In the following parts of this paper, the sensitivity analysis provides a better understanding of the effects of tissue transport parameters on drug delivery.

Results

Simulation of interstitial fluid flow for baseline values (Table 1) predicts that IFP has its greatest value in the tumor center. IFP is non-dimensionalized by effective pressure. The effective pressure, P _eff , is a parameter defined by intravascular pressure, plasma osmotic pressure, and interstitial osmotic pressure by Equation (1). The non-dimensionalized pressure is defined by Equation (2)

Table 1 Interstitial transport properties used in numerical simulations

$${\mathit{P}}_{\mathit{eff}}={\mathit{P}}_{\mathit{B}}^{\mathit{t}}-{\mathit{\sigma }}_{\mathit{s}}^{\mathit{t}}\left({\mathit{\pi }}_{\mathit{B}}^{\mathit{t}}-{\mathit{\pi }}_{\mathit{i}}^{\mathit{t}}\right)$$

(1)

$${\mathit{P}}_{\mathit{n}}=\frac{{\mathit{P}}_{\mathit{i}}}{{\mathit{P}}_{\mathit{eff}}}=\frac{{\mathit{P}}_{\mathit{i}}}{{\mathit{P}}_{\mathit{b}}^{\mathit{t}}-{\mathit{\sigma }}_{\mathit{s}}^{\mathit{t}}\left({\mathit{\pi }}_{\mathit{b}}^{\mathit{t}}-{\mathit{\pi }}_{\mathit{i}}^{\mathit{t}}\right)}$$

(2)

where t stands for tumor tissue. Parameters used in Equation (1) and (2) are introduced in the method section.

Non-dimensionalized IFP along transverse and vertical lines (shown in Figure 1), are illustrated in Figures 2 and 3. The maximum value of IFP occurs in the tumor center (Figure 2). This value is equal to P _eff , 1.53 kPa, except for the case of λ = 0.1 (Prolate shape). IFP has uniform distribution at tumor region and in the inner boundary (for more detail see Figure 1 and boundary condition section) IFP falls down sharply to around 150 Pa as shown in Figure 2 for all shapes except for prolate one (λ = 0.1). In normal tissue, pressure has uniform distribution close to outer boundary (for more detail, see Figure 1 and boundary condition section) then it decreases smoothly to peripheral pressure. For prolate shape as shown in Figure 2, IFP reduces smoothly from tumor center to the outer boundary, and IFP has a non-zero gradient in whole domain. Figure 2 shows the IFP along vertical line. Results are very similar to IFP results in Figure 2. Only for a prolate shape IFP has different pattern from what is observed in Figure 2. Maximum value of IFP in the prolate shape does not occur in the tumor center and takes place somewhere between the tumor center and tumor periphery.

Figure 1

Schematic of considered geometry and boundary conditions. The transverse and vertical lines are used to show results along them.

Figure 2

Interstitial fluid pressure profile along transverse line. The 1 in “x” axis shows the boundary of tumor (inner boundary).

Figure 3

Interstitial fluid pressure profile along vertical line. The 1 in “x” axis shows the boundary of tumor (inner boundary).

IFV distribution along the vertical and transverse lines is presented in Figures 4 and 5. Maximum value of IFV occurs close to the inner boundary. Also, in normal tissue, IFV reaches zero far from the inner boundary. However, for the prolate shape, velocity is not zero, especially along transvers line (Figure 4).

Figure 4

Interstitial fluid velocity profile along transverse line. The 1 in “x” axis shows the boundary of tumor (inner boundary).

Figure 5

Interstitial fluid velocity profile along vertical line. The 1 in “x” axis shows the boundary of tumor (inner boundary).

DCIF is simulated in two cases of injection. In the first case, the continuous injection which leads to constant plasma concentration (C _p  = constant) is considered and in the second case, the bolus injection in which the plasma concentration decreases with time exponentially (${\mathit{C}}_{\mathit{p}}={\mathit{C}}_{\mathit{p}}^0{\mathit{e}}^{-ln\left(0.5\right)\mathit{t}/\mathit{\tau }}$) is studied, in which τ is the drug half-life in plasma (Table 2). DCIF are non-dimensionalized by C _P for continuous injection and ${\mathit{C}}_{\mathit{P}}^0$ for bolus injection, respectively. Average of non-dimensionalized DCIF for two injection cases in different times is shown in Figures 6 and 7. DCIF of prolate shape (λ = 0.1) has the maximum value and DCIF of oblate shape (λ = 10) has the minimum value. Other considered shapes show the similar transient behavior.

Non-dimensionalized DCIF along two lines (transverse and vertical) is illustrated in Figures 8 and 9. The bolus injection results are presented in 8 hr post injection in which the concentration is maximum based on Figure 7 and for continuous injection is presented at final time of simulation, in 72 hr post injection. The profiles of two types of injections are similar in spite of different values of DCIF. As observed in IFV and IFP profiles, the prolate shape is different from the other shapes in DCIF distribution. A bump is observed in DCIF curves at the inner boundary of all tumor shapes. The normal tissue has uniform distribution of DCIF except near the boundaries (inner and outer). DCIF distribution in normal tissue is the same for all tumor shapes.

Table 2 Parameters of solute transport model used in numerical simulations

Figure 6

Average non-dimensionalized DCIF in tumor region in during time for continuous injection.

Figure 7

Average non-dimensionalized DCIF in tumor region in during time for bolus injection.

Figure 8

Non-dimensionalized DCIF profile along transverse line for normal. Left for bolus injection and right for continuous injection. The 1 in “x” axis shows the boundary of tumor (inner boundary).

Figure 9

Non-dimensionalized DCIF profile along vertical line. Left for bolus injection and right for continuous injection. The 1 in “x” axis shows the boundary of tumor (inner boundary).

2D contours of DCIF in tumor region for bolus injection in 8 hr post injection are shown in Figure 10 for all tumor shapes. Results show that DCIF for tumors with λ = 0.1 and λ = 10 has different distributions. Generally, in the inner boundary, DCIF has its maximum value. In Figures 11 and 12, Peclet number distribution along two lines for continuous injection is shown. Peclet number demonstrates the ratio of convection to diffusion across the microvessel wall. Results show that in the tumor region Peclet number is zero except for prolate shape. Peclet number for the prolate shape especially along short radius (transverse line) is greater than zero.

Figure 10

2D contour of DCIF in tumor region for bolus injection in 8 hr post injection for all tumor shapes.

Figure 11

Peclet number (ratio of convection to diffusion across the microvessel wall) along transverse line.

Figure 12

Peclet number (ratio of convection to diffusion across the microvessel wall) along vertical line.

Some of tissue transport parameters mentioned (effective pressure, hydraulic conductivity, and tumor size) in Table 1 are investigated for sensitivity analysis. The values of these parameters are selected near the ranges reported in the literature [2, 15, 21].

Figure 13 shows the influence of changing effective pressure, P _eff , on numerical results of tumor shapes with the same volume. Maximum value of IFP in tumor region (Figure 13a) linearly increases when P _eff increases. Only in the prolate shape the pressure increases more than other tumor shapes. Average of IFV on inner boundary also has the same pattern and linearly changes with P _eff (Figure 13b). The changes in average of DCIF in tumor region for two cases of injection are shown in Figure 13c and d. For these two cases of injection, the prolate shape is more sensitive to P _eff changes than other tumor shapes and DCIF changes around 20% when the P _eff changes from 500 Pa to 2500 Pa. The other shapes do not show significant variations in DCIF for these ranges of P _eff .

Figure 13

The effect of different values of effective pressure on IFP, IFV and DCIF. a) The variation of maximum IFP in tumor region. b) The variation of averaged IFV at tumor boundary. c) Average of DCIF in tumor region for continuous injection. d) Average of DCIF in tumor region for bolus injection.

Different tumor sizes are studied with changes in their volume. One of important metric of disease development and response to tumor therapy with drug is volume of tumors [22–25]. To quantify response to therapeutic regimens and also assess disease progression, tumor volume is used as a metric in many studies, such as Char et al. [26], Jensen et al. [27] and Gass et al. [28].

As shown in Figure 14a, IFP has less value than P _eff when the tumor volume is smaller than 1 cm³. When the tumor volume is in the order of 1 cm³, IFP reaches P _eff and by increasing the tumor radius, IFP remains constant. Average of IFV on the inner boundary generally decreases with increasing tumor size, Figure 14b. Only in prolate shape in small radii, IFV increases with the tumor size. As shown in Figure 14c and d, the mean values of DCIF have the greatest value in the smallest tumor. Also, the prolate shape has the highest value of DCIF in all studied tumors.

Figure 15 shows the sensitivity of IFF parameters and DCIF to hydraulic conductivity changes. Results show that in all tumor shapes, if hydraulic conductivity increases, the maximum value of IFP remains constant and then decreases sharply. Average of IFV increases by increasing hydraulic conductivity. DCIF increases by increasing hydraulic conductivity and then reaches a constant value in spite of increasing hydraulic conductivity.

Figure 14

The effect of different values of tumor size on IFP, IFV and DCIF. a) The variation of maximum IFP in tumor region. b) The variation of averaged IFV at tumor boundary. c) Average of DCIF in tumor region for continuous injection. d) Average of DCIF in tumor region for bolus injection.

Figure 15

The effect of different values of hydraulic conductivity on IFP, IFV and DCIF. a) The variation of maximum IFP in tumor region. b) The variation of averaged IFV at tumor boundary. c) Average of DCIF in tumor region for continuous injection. d) Average of DCIF in tumor region for bolus injection.

Discussion

This study presents DCIF, IFP, and IFV in solid tumors surrounded by normal tissue in two types of injection; bolus and continuous one. The model used in this study investigates the effect of two characteristics of tumors on concentration distribution; tumor shape and size.

Results of high IFP in tumors are discussed in our previous studies [2, 16, 17] and in the experimental results of Arifin et al. [29] and Huber et al. [30]. Maximum value of IFP in spherical tumors (1 cm radius) for baseline values in Table 1 is 1529.5 Pa which is close to the studies of Jain et al. [31], Chauhan et al. [32], and Arfin et al. [33]. The current results are verified by experimental data of IFP measured by Nielsen et al. [34]. In their work, the wick-in-needle technique is used to measure IFP in two types of tumors with the same volume as tumors in the current study. They found IFP in the range of 1400 Pa to 1600 Pa.

IFV on tumor boundary in spherical tumors with baseline values is around 0.05 μm/s which is at the same order of the prediction of Jain et al. [31] and experimental observation of Hompland et al. [35]. Also, the profile of drug concentration for simulation with baseline values for spherical tumor in different times (Figures 6 and 7) agrees well with Baxter and Jain’s predictions [8]. Results show that DCIF for prolate shape (λ = 0.1) always has the greatest value. Results of transient DCIF (Figures 6 and 7) show that drug delivery is much easier in the prolate tumors. Also, the oblate shape (λ = 10) has the most resistance to drug delivery. Following of this section, the reason of these phenomena is investigated.

The uniform values of IFP for all tumor shapes except prolate one (Figures 2 and 3) is equal to P _eff . Non-uniform IFP in prolate shape results in the maximum value of DCIF among other tumor shapes with the same volume. Equation (17) is able to legitimize this behavior of tumor shapes. The source term (the last term in the right hand side of Equation (17)) includes diffusion and convection terms. The convection term depends on the differences between IFP and P _eff based on starling’s law. This pressure difference is close to zero (Figures 2 and 3) for all tumor shapes except for the prolate one and therefore the convection term only for prolate shape has non-zero value. Because of the non-zero values of Peclet number for prolate shape in tumor region (Figures 11 and 12), the concentration for this shape is affected not only by diffusion rate from vessels but also by convection rate from vessels. The convection rate leads to higher level of DCIF in prolate shape than that of other tumor shapes (Figures 6, 7, 8, 9 and 10). The non-zero values of Peclet number for prolate shape are seen in an image based work of Zhao et al. as well [13]. In normal tissue, Peclet number shows that drug delivery from microvessel to tissue is done by both mechanisms of transfer, convection and diffusion.

The other effect of uniform pressure is zero pressure gradients in all tumor shapes except for prolate one. Because of zero IFP gradient and based on Darcy’s equation, IFV is close to zero in tumor tissue except for prolate shape. The close to zero value of IFV is predicted in a few numerical studies such as Welter and Rieger [36]; Roy and Riahi [37]; and experimental results of Hompland et al. [35]. Zero IFV results in a negligible convection transport (the second term of the right hand side of Equation (17)) and consequently the convection transport does not affect drug distribution; and the diffusion transport (the first term of the right hand side of Equation (17)) is the only reason of drug transport in all tumor shapes except for prolate one. Therefore, non-zero IFV in prolate shape, also seen in Zhao et al. [13], is another reason of higher DCIF values in this tumor shape with respect to other shapes.

The sharp pressure gradient (Figures 2 and 3) and highest value of IFV (Figures 4 and 5) in the inner boundary for all tumor shapes increases drug transport and make a bump at this boundary for DCIF.

The sensitivity analysis of effective pressure shows that the effective pressure does not have too much effect on DCIF. In all tumor shapes, DCIF for a wide range of effective pressure changes smoothly; however, in prolate shape, this change is sharper and increases by effective pressure (Figure 13c and d). As mentioned, in the tumor region for all tumor shapes except for prolate one the convection rate, which depends on P _eff , has a negligible role on drug distribution; therefore, P _eff cannot have significant effect on drug concentration. However, in prolate shape, increasing P _eff increases pressure difference between IFP and P _eff and consequently the convection rate from vessels in the tumor region; therefore, DCIF in prolate shape is sensitive to P _eff changes.

The tumor volume shows more effects on IFP, IFV, and DCIF than other investigated parameters such as effective pressure and hydraulic conductivity. The increasing tumor volume increases significantly IFP (Figure 14a). The dependency of IFP to tumor volume is observed in experimental study of Gutmann et al. [38], Hompland et al. [39], and Leguerney et al. [40], as well. When the tumor volume is in the order of 1 cm³, the sensitivity of IFP to tumor size decreases. The independency of IFP to tumor volume in in large tumors is observed in the study of Leguerney et al. [40], as well. In their work, IFP changes very slowly with tumor volume. Since the high IFP is introduced as the main barrier of drug delivery to tumors, IFP increasing with tumor volume leads to DCIF decrease in these tumors. This reduction of DCIF with tumor size is observed in Au et al. [41], as well. Since IFP reaches around the effective pressure with increasing tumor volume, the convection rate is vanished and the diffusion rate reaches a constant value, and consequently the sensitivity of DCIF to tumor size reduces. Lower IFP in the small tumor sizes leads to increase the convection rate of source term in Equation (17). Therefore, it is expected to have a better drug distribution in small tumors.

Results show that IFV, IFP, and DCIF are sensitive to tissue hydraulic conductivity changes. The hydraulic conductivity is appeared only in Darcy’s law. This parameter has a direct effect on IFF and indirect effect on DCIF. Theoretical analysis shows that the $\frac{{\mathit{P}}_{\mathit{i}}}{{\mathit{P}}_{\mathit{eff}}}$ in tumor region is proportional to $1-\frac{1}{\sqrt{\mathit{\kappa }}}/sinh\left(\frac{1}{\sqrt{\mathit{\kappa }}}\right)$[8] (κ is hydraulic conductivity, see material section). In low values of κ, the $\frac{1}{\sqrt{\mathit{\kappa }}}/sinh\left(\frac{1}{\sqrt{\mathit{\kappa }}}\right)$ is negligible and P _i is close to P _eff (Figure 15a). Increasing κ, increases the $\frac{1}{\sqrt{\mathit{\kappa }}}/sinh\left(\frac{1}{\sqrt{\mathit{\kappa }}}\right)$ and leads to sharp decrease in IFP. This dependency is also observed by McCarty and Johnson [42]. For high values of κ, IFP is very low and negligible in comparison to P _eff . IFP reduction from effective pressure increases IFV around 5 times (Figure 15b). The hydraulic conductivity affects DCIF by convection rate from vessels (as mentioned this value depends on difference between IFP and P _eff ). In low values of κ, since IFP is equal to P _eff , the effect of convection rate is not significant and DCIF remains constant. Increasing κ increases the convention rate and consequently DCIF.

When hydraulic conductivity increases two to three orders of magnitude, the mean values of DCIF are two times greater than the average of DCIF for baseline values in Table 1. However, after a specific value of hydraulic conductivity, DCIF changes smoothly and reaches a constant value because IFP is very low and convection rate only depends on P _eff .

Conclusions

A numerical approach which couples the mathematical model of the lymphatic system and the interstitial flow with the mathematical model of solute transport demonstrates that DCIF is affected by two transport mechanisms, convection and diffusion.

Drug convection and drug transport from microvessel to tumor are blocked by high interstitial pressure (IFP) which is uniformly distributed in most part of the tumor. The large pressure gradient results in an outward convective flow that washes out the drug extravasated from microvessels at the tumor periphery. This study shows that when there is IFP gradient in the tumor region instead of uniform IFP distribution which occurs in some tumor shapes, DCIF is greater than that of the uniform one.

As the effects of osmotic and intravascular pressure on convection rate are negligible in most of tumor shapes, the dependency of DCIF to these parameters is very low.

The hydraulic conductivity which is another considered parameter in sensitivity analysis has significant effect on drug distribution since it increases the convection rate from vessels.

Method

Mathematical model of interstitial fluid transport

This section introduces the mathematical model of interstitial fluid transport in macroscopic scale [2, 16]. Since normal and tumor tissue have characteristics the same as porous media, fluid flow behavior is defined by coupling the fluid flow governing equations. The mass balance or continuity equation for steady state incompressible fluid in the porous media with source and sink of mass is [16]:

$$\nabla .{\mathit{v}}_{\mathit{i}}={\mathit{\varphi }}_{\mathit{B}}-{\mathit{\varphi }}_{\mathit{L}}$$

(3)

where

v _i : The interstitial fluid velocity,

ϕ _B : The source term, extravasation from microvessels, and

ϕ _L : The drainage term, elimination by lymphatic system.

In biological tissues, the fluid source is evaluated through Starling's law as follows [16]:

$${\mathit{\varphi }}_{\mathit{B}}=\frac{{\mathit{L}}_{\mathit{P}}\mathit{S}}{\mathit{V}}\left({\mathit{P}}_{\mathit{B}}-{\mathit{P}}_{\mathit{i}}-{\mathit{\sigma }}_{\mathit{s}}\left({\mathit{\pi }}_{\mathit{b}}-{\mathit{\pi }}_{\mathit{i}}\right)\right)$$

(4)

where

P _i : Interstitial fluid pressure,

P _B : Blood pressure in microvessel,

$\frac{\mathit{S}}{\mathit{V}}$ : The surface area per unit volume of tissue for transport in the interstitium,

π _B : Microvessel oncotic pressure,

π _i : Interstitial oncotic pressure,

L _p : The hydraulic conductivity of the microvessel wall, and

σ _s : Osmotic reflection coefficient.

and the lymphatic system is related to the pressure difference between the interstitium and the lymphatic vessels and is considered only for normal tissues [43]:

$${\mathit{\varphi }}_{\mathit{L}}\left(\mathit{r}\right)=\left\{\begin{array}{ll}\frac{{\mathit{L}}_{\mathit{PL}}{\mathit{S}}_{\mathit{L}}}{\mathit{V}}\left({\mathit{P}}_{\mathit{i}}-{\mathit{P}}_{\mathit{L}}\right)\hfill & \mathit{Normal}\mathit{tissue}\hfill \\ 0\hfill & \mathit{Tumor}\mathit{tissue}\hfill \end{array}\right.$$

(5)

where

ϕ _L : The volumetric flow rate into the lymphatic,

$\frac{{\mathit{L}}_{\mathit{PL}}{\mathit{S}}_{\mathit{L}}}{\mathit{V}}$: The lymphatic filtration coefficient, and

P _L : The hydrostatic pressure of the lymphatic.

The momentum balance equation in its general form can be written as Equation (6) [44]:

$$\begin{array}{ll}\mathit{\rho }\left(\frac{\partial {\mathrm{v}}_{\mathit{i}}}{\partial \mathit{t}}+\left({\mathrm{v}}_{\mathit{i}}.\nabla \right){\mathrm{v}}_{\mathit{i}}\right)=\nabla ·[-{\mathit{P}}_{\mathit{i}}+\mathit{\mu }\left(\nabla {\mathrm{v}}_{\mathit{i}}+{\left(\nabla {\mathrm{v}}_{\mathit{i}}\right)}^{\mathit{T}}\right)& -\frac{2\mathit{\mu }}{3}\left(\nabla ·{\mathrm{v}}_{\mathit{i}}\right)]\\ -\left(\frac{\mathit{\mu }}{\mathit{K}}\right){\mathrm{v}}_{\mathit{i}}+\mathit{F}\end{array}$$

(6)

where

K: The permeability of the porous medium,

ρ : The density,

μ: The viscosity, and

F: The volume forces.

Since interstitial fluid is a Newtonian fluid and has low velocity through the tissues, Equation (6) in the steady state condition is simplified to Darcy’s law [16]:

$$\nabla {\mathit{P}}_{\mathit{i}}=-\left(\frac{\mathit{\mu }}{\mathit{K}}\right){\mathrm{v}}_{\mathit{i}}$$

(7)

The K/μ is defined as hydraulic conductivity of the interstitium, κ:

$${\mathrm{v}}_{\mathit{i}}=-\mathit{\kappa }\nabla {\mathit{P}}_{\mathit{i}}$$

(8)

Combination of momentum equation (Equation (8)) and the continuity equation (Equation (3)) results in

$$-\nabla ·\left(\mathit{\kappa }\nabla {\mathit{P}}_{\mathit{i}}\right)={\mathit{\varphi }}_{\mathit{B}}-{\mathit{\varphi }}_{\mathit{L}}$$

(9)

When κ is constant, the interstitial pressure can be expressed by

$$-\mathit{\kappa }{\nabla }^2{\mathit{P}}_{\mathit{i}}={\mathit{\varphi }}_{\mathit{B}}-{\mathit{\varphi }}_{\mathit{L}}$$

(10)

By substituting Equations (4) and (5) into Equation (10):

$$-\mathit{\kappa }{\nabla }^2\mathit{p}=\left\{\begin{array}{ll}\frac{{\mathit{L}}_{\mathit{P}}\mathit{S}}{\mathit{V}}\left({\mathit{P}}_{\mathit{B}}-{\mathit{P}}_{\mathit{i}}-{\mathit{\sigma }}_{\mathit{s}}\left({\mathit{\pi }}_{\mathit{b}}-{\mathit{\pi }}_{\mathit{i}}\right)\right)-\frac{{\mathit{L}}_{\mathit{PL}}{\mathit{S}}_{\mathit{L}}}{\mathit{V}}\left({\mathit{P}}_{\mathit{i}}-{\mathit{P}}_{\mathit{L}}\right)\hfill & \mathit{Normal}\mathit{tissue}\hfill \\ \frac{{\mathit{L}}_{\mathit{P}}\mathit{S}}{\mathit{V}}\left({\mathit{P}}_{\mathit{B}}-{\mathit{P}}_{\mathit{i}}-{\mathit{\sigma }}_{\mathit{s}}\left({\mathit{\pi }}_{\mathit{b}}-{\mathit{\pi }}_{\mathit{i}}\right)\right)\hfill & \mathit{Tumor}\mathit{tissue}\hfill \end{array}\right.$$

(11)

Macroscopic solute transport

Molecular transport in tumors is based on the conservation laws for mass and momentum. The interstitial transport of drug is governed by the convection diffusion equation; therefore, the general equation for the mass balance of solutes can be written as:

$$\frac{\partial \mathit{C}}{\partial \mathit{t}}=\nabla ·\left({\mathit{D}}_{\mathit{eff}}·\nabla \mathit{C}\right)-\nabla ·\left({\mathit{v}}_{\mathit{i}}\mathit{C}\right)+\left({\mathrm{\Phi }}_{\mathit{B}}-{\mathrm{\Phi }}_{\mathit{L}}\right)$$

(12)

where

C : The solute concentration based on tissue volume,

Φ _B : The rate of solute transport per unit volume from microvessel into the interstitial space,

Φ _L : The rate of solute transport per unit volume from the interstitial space into lymphatic vessels, and

D _eff : The effective diffusion tensor.

For an isotropic and uniform diffusion in tissues, equation (12) can be written as:

$$\frac{\partial \mathit{C}}{\partial \mathit{t}}={\mathit{D}}_{\mathit{eff}}{\nabla }^2\mathit{C}-\nabla ·\left({\mathit{v}}_{\mathit{f}}\mathit{C}\right)+\left({\mathrm{\Phi }}_{\mathit{B}}-{\mathrm{\Phi }}_{\mathit{L}}\right)$$

(13)

The solute transport rate across the lymphatic vessels can be considered as [15]

$${\mathrm{\Phi }}_{\mathit{L}}=\left\{\begin{array}{cc}\hfill {\mathit{\varphi }}_{\mathit{L}}\mathit{C}\hfill & \hfill \mathit{Normal}\mathit{Tissue}\hfill \\ \hfill 0\hfill & \hfill \mathit{Tumor}\mathit{Tissue}\hfill \end{array}\right.$$

(14)

The solute transport rate across the microvessel can be represented by Patlak equation [45]:

$${\mathrm{\Phi }}_{\mathit{B}}={\mathit{\varphi }}_{\mathit{B}}\left(1-{\mathit{\sigma }}_{\mathit{f}}\right){\mathit{C}}_{\mathit{P}}+\frac{\mathit{PS}}{\mathit{V}}\left({\mathit{C}}_{\mathit{P}}-\mathit{C}\right)\frac{\mathit{Pe}}{{\mathit{e}}^{\mathit{Pe}}-1}$$

(15)

where

$$\mathit{Pe}=\frac{{\mathit{\varphi }}_{\mathit{B}}\left(1-{\mathit{\sigma }}_{\mathit{f}}\right)\mathit{V}}{\mathit{PS}}$$

(16)

ϕ _B : The fluid flow rate per unit volume of tissue across the microvessel wall,

σ _f : The filtration reflection coefficient,

P : The microvessel permeability coefficient,

S/V : The microvessel surface area per unit volume of tissue,

C _p : Solute concentration in the plasma.

By substituting Equations (15) and (14) in to Equation (13):

$$\left\{\begin{array}{l}\begin{array}{l}\frac{\partial \mathit{C}}{\partial \mathit{t}}={\mathit{D}}_{\mathit{eff}}{\nabla }^2\mathit{C}-\nabla ·\left({\mathit{v}}_{\mathit{i}}\mathit{fC}\right)\\ +\left({\mathit{\varphi }}_{\mathit{B}}\left(1-{\mathit{\sigma }}_{\mathit{f}}\right){\mathit{C}}_{\mathit{P}}+\frac{\mathit{PS}}{\mathit{V}}\left({\mathit{C}}_{\mathit{P}}-\frac{\mathit{C}}{{\mathit{K}}_{\mathit{AV}}}\right)\frac{\mathit{Pe}}{{\mathit{e}}^{\mathit{Pe}}-1}-{\mathit{\varphi }}_{\mathit{L}}\mathit{C}\right)\mathit{Normal}\mathit{region}\end{array}\hfill \\ \begin{array}{l}\frac{\partial \mathit{C}}{\partial \mathit{t}}={\mathit{D}}_{\mathit{eff}}{\nabla }^2\mathit{C}-\nabla ·\left({\mathit{v}}_{\mathit{i}}\mathit{fC}\right)\\ +\left({\mathit{\varphi }}_{\mathit{B}}\left(1-{\mathit{\sigma }}_{\mathit{f}}\right){\mathit{C}}_{\mathit{P}}+\frac{\mathit{PS}}{\mathit{V}}\left({\mathit{C}}_{\mathit{P}}-\frac{\mathit{C}}{{\mathit{K}}_{\mathit{AV}}}\right)\frac{\mathit{Pe}}{{\mathit{e}}^{\mathit{Pe}}-1}\right)\mathit{Tumor}\mathit{region}\end{array}\hfill \end{array}\right.$$

(17)

Boundary conditions

A tumor surrounded by normal tissue is considered in this study. The tumor shape is considered ellipsoid. A 2D section of the geometry and boundaries are shown in Figure 1. The three boundaries are indicated for geometry: a) the center of tumor, b) the boundary between tumor and normal tissue, named inner boundary, c) the boundary at the outer edge of geometry, named outer boundary. The appropriate boundary conditions are implemented for Equations (11) and (17). The no flux boundary condition is applied at the center of the tumor; i.e. [16],

$$\nabla {\mathit{P}}_{\mathit{i}}=0\mathit{for}\mathit{r}=0$$

(18)

$${\mathit{D}}_{\mathit{eff}}\nabla \mathit{C}+{\mathit{v}}_{\mathit{i}}\mathit{C}=0\mathit{for}\mathit{r}=0$$

(19)

The continuity of pressure and velocity for Equation (11) and concentration and its flux for Equation (17) are considered as appropriate boundary conditions for inner boundary:

$$\begin{array}{l}{\left.-{\mathit{\kappa }}_{\mathit{t}}\nabla {\mathit{P}}_{\mathit{i}}\right|}_{{\mathit{\Omega }}^{-}}={\left.-{\mathit{\kappa }}_{\mathit{n}}\nabla {\mathit{P}}_{\mathit{i}}\right|}_{{\mathit{\Omega }}^{+}}\\ {\left.{\mathit{P}}_{\mathit{i}}\right|}_{{\mathit{\Omega }}^{-}}={\left.{\mathit{P}}_{\mathit{i}}\right|}_{{\mathit{\Omega }}^{+}}\end{array}$$

(20)

$$\begin{array}{l}{\left.\left({\mathit{D}}_{\mathit{eff}}^{\mathit{t}}\nabla \mathit{C}+{\mathit{v}}_{\mathit{i}}\mathit{C}\right)\right|}_{{\mathit{\Omega }}^{-}}={\left.\left({\mathit{D}}_{\mathit{eff}}^{\mathit{n}}\nabla \mathit{C}+{\mathit{v}}_{\mathit{i}}\mathit{C}\right)\right|}_{{\mathit{\Omega }}^{-}}\\ {\left.\mathit{C}\right|}_{{\mathit{\Omega }}^{-}}={\left.\mathit{C}\right|}_{{\mathit{\Omega }}^{+}}\end{array}$$

(21)

where Ω⁻ and Ω⁺ indicate the tumor and normal tissue at the inner boundary.

For outer boundary, that the interstitial pressure is constant; the Dirichlet type of boundary condition must be applied [16]:

$${\mathit{P}}_{\mathit{i}}=\mathit{\text{for outer region}}$$

(22)

And for concentration, the open boundary condition is used in the outer region [46]. The Open Boundary is used to set up mass transport across boundaries where both convective inflow and outflow can occur and defined by Equation (23):

$$-\mathbf{n}·\nabla \mathrm{C}=0$$

(23)

where n is the normal vector.

Geometry

To investigate the effect of tumor shape on drug delivery, 5 different shapes are considered. The 3D fundamental shape of tumors is assumed to be an ellipsoid in different studies. The assumption of considering ellipsoid tumor shape is investigated in many research such as breast cancer [47], prostate cancer [48, 49] cervical cancer [50], glioma cancer [51], and others [52]. Based on this mentioned reason, ellipsoidal shapes of tumors are considered in this study. Different shapes are produced by changing ratio of two radii of ellipsoid shown in Figure 1. This ratio (λ = b/a) is changed from 0.1 (prolate) to 10 (oblate). In all shapes shown in Figure 16 the volume of the tumor is remained constant and equal to the volume of spherical tumor with radius R. The baseline value of R is 1 cm. The values of R used in sensitivity analysis are changed from 0.1 cm to 2.5 cm. This range is obtained from the literature and is close to the experimental observations [40].

Figure 16

Different shapes of solid tumors considered in this study. The five ellipsoid shapes changed from oblate to prolate are studied.

Model parameterization

The interstitial transport properties for normal and tumor tissue are listed in Table 1. These values are used as baseline and some of them are investigated and changed in specified ranges for sensitivity analysis.

The parameters of solute transport model taken from Baxter and Jain [8] are listed in Table 2. Although, the numerical model is applicable for any type of solute, in present study the properties of Fragment antigen-binding (F(ab’)₂) as a sample is used.

Numerical solution

A systematic flow chart is illustrated in Figure 17 to clarify the computational techniques involved in this paper. The criterion for the convergence of iterative solution based on element–based finite volume (EB-FV) method is to reduce the residual by 6 orders of magnitudes. The details of numerical method is mentioned in our previous works [2, 16, 53].

Figure 17

Algorithm of numerical simulation used for calculating interstitial fluid flow parameter (IFV, IFP) and solute transport parameter (DCIF).