Mathematical Analysis for the Effects of Medicine Supplies to a Solid Tumor

Abstract

Objective: 1. Interpretation of the variations of solute medicine amount in blood vessels and TAF concentration with respect to the flow rates of injected drugs into liver and heart. 2. Description of the alteration of tumor cell density versus the time and radius variations. Methodology: Step 1. Compartmental analysis is adopted for the concentration of chemotaxis caused by injected substances $L$ and $H$ based on the assumption: two different medicines $I_1$ and $I_2$ are injected into heart and liver to recover the functions of each organ, respectively, without any side effects. Step 2. A partial differential equation is derived for the growth of TAF considering the diffusion of TAF and the rate of decay of TAF according to the disturbance of medicine $M$ in blood vessels. Step 3. A partial differential equation is derived for the motion of tumor cells in the lights of random motility and chemotaxis in response to TAF gradients. Step 4. Exact solutions are obtained for the concentration of chemotaxis caused by injected substances $L$ and $H$ under the assumption that the loss of mass is proportional to mass itself. Step 5. Exact solution is obtained for the partial differential equation describing the growth of TAF using the separation of variables. Step 6. A finite volume approach is executed to search approximated solutions due to the complexity of the partial differential equation describing the motion of tumor cells. Results: 1. The concentration of medicine ($M$) decreases as the ratio of flow rate from heart into vessel to flow rate from liver into heart ($\frac{k_1}{k_2})$ increases. 2. TAF concentration increases with the growth of the value of ratio $\frac{k_1}{k_2}$ and TAF shows the smallest concentration when the flow rate of each injected medicine is similar. 3. Tumor cells react highly sensitive as soon as medicine supplies and tumor cell’s density is decreased drastically at the moment of medicine injection. 4. Tumor cell density decreases exponentially at an early stage and the density decrease is developed in a fluctuating manner along the radius. Conclusions: 1. The presented mathematical approach has the potential for the profound analysis of the variations of solute medicine amount in blood vessels, TAF concentration, and the alteration of tumor cell density according to the functional recoveries of liver and heart. 2. The mathematical approach may be applicable in the investigation of tumor cell’s behavior on the basis of complex interaction among five represented organs: kidney, liver, heart, spleen, and lung. A mathematical approach is developed to describe the variation of a solid tumor cell density in response to drug supply. The investigation is progressed based on the assumption that two different medicines, $I_1$ and $I_2,$ are injected into heart and liver with flow rates $k_1$ and $k_2$ to recover the functions of each organ, respectively. A medicine function system for the reactions of tumor angiogenic factors (TAF) to medicine injection is obtained using a compartmental analysis. The mathematical governing equations for tumor cells motion are derived taking into account random motility and chemotaxis in response to TAF gradients and a finite volume method with time-changing is adopted to obtain numerical solutions due to the complexity of the governing equations. The variation of the flow rates $k_1$ and $k_2$ exerts profound influences on the concentration of medicine, and similar flow rate of $k_1$ and $k_2$ produces the greatest amount of medicine in blood vessels and suppresses strong inhibition in TAF movement. Tumor cells react very sensitively to drug injection and the tumor cell density decreases to less than 20% at an early stage of administration. However, the density of tumor cell diminishes slowly after the early stage of sudden change and the duration for complete therapy of tumor cells requires a long time.

1. Introduction

Cancer is the situation of abnormal cell growth in a specific part of the body and spread to other parts beyond boundaries and is becoming the major death cause of the world population. Most tumors undergo insufficient nutrition at an early avascular stage and have difficulty in wastes disposal, and thus the transfer of nutrients and wastes has been progressed through simple diffusion. Tumor cells stimulate endothelial cells located in the interior surface of blood vessels to obtain solutes, such as oxygen, for growth. If the tumor reaches dormant state, the cancer cells release a diffusible chemical substance called tumor angiogenesis factor (TAF). Angiogenesis induced by tumor cells contains cancerous cells that secrete a number of chemicals into surrounding tissues [1,2], and capillary sprouts are formed in response to chemical stimuli. Anderson and Chaplain [3] analyzed the formation of the capillary sprout network in response to chemical stimuli supplied by a solid tumor using both continuous and discrete models. Macklin and Lowengrub [4] investigate the effects of the interaction between the genetic characteristics of the tumor and the tumor micro-environment, and mathematical analysis for the growths of non-necrotic and necrotic tumors with time delays achieved by Cui and Xu [5].

Nowadays, an active analysis of cancer cells is comprehensively performed to understand the development processes of tumor cells and develop anticancer drugs. Successful treatment of tumor depends on the efficient deliveries of anticancer drugs. Tumor cell propagation creates an abnormal vascular structure causing an elevated interstitial fluid velocity, and chaotic vascular structure appeared commonly in advanced tumor disturbs the efficient drug delivery to solid tumors [6]. Vascular normalization using anti-angiogenic factors, thus, reduces the interstitial fluid velocity and strengthens the vessel perfusion and improves both drug and nutrients delivery to tumor cells [7,8]. By dint of complex and unpredictable actions among various factors in tumor cell growth mathematical approach has become an essential method in the study of drug delivery. E.M. Kashkooli et al. [9] investigates the effect of vascular normalization in executing efficient drug delivery through computational and mathematical modeling approaches.

Moreover, the extensive dosage of chemical medicines to control cancer growth is a discussing issue in diversified perspectives. Magni et al. [10] explored the cancer growth dynamics in response to anticancer agent administration in xenograft models, using an ordinary differential equation. Jackson and Byrne [11] used partial differential equations to investigate vascular tumor behaviors against chemotherapeutic treatment. Maeda et al. [12] suggested a tumor-targeted delivery strategy responding to anticancer medicines for patients who are in the stage of aggravated cancer growth. They studied the pathophysiological mechanisms for various tumor blood vessels containing properties such as enhanced permeability, and different endurance and structure. The impact of apheresis platelet supernatants of different storage periods on tumor cell growth was explored by Fei et al. [13]. They optimized transfusion timing and provided a reference for reducing the risks of platelet transfusion in cancer patients. Hang et al. [14] found that surrounding neighbor cells in tumor induce non-cell autonomous autophagy, and tumor growth can be suppressed significantly by genetic ablation of autophagy induction. A mathematical approach was developed to investigate the growth of both tumor cells and microorganisms expressing the influences of migrating cells, proliferating cells, and nutrient by Alpna [15].

Meanwhile, metabolism of human body is in close contact with the balanced functions of five representative organs: kidney, liver, heart, spleen, and lung. The kidney is the first generated organ in an unborn child and takes charge of water and salinity in human body. The kidney is deeply related to joints in our body and the function of joints in the knee and waist reacts sensitively to kidney conditions [16,17,18]. The liver is generated with the assistance of the matured kidney. The liver purifies corrupted blood and supports the health of tendons in our body [19,20]. The heart organ is born with the help of the liver and plays important roles in maintaining healthy blood and blood vessels [21]. The spleen is generated by the virtue of the heart’s assistance and has jurisdiction for the preservation of resilient skin in the human body [22]. The spleen assists in the creation of the lung, which is the last generated among the five representative organs, and the lung takes charge of the skin pores [10,23]. The interactions of five organs are critically key to preserving good physical shape, and metabolism of the human body is controlled by the co-existence, conflict, and contradiction of five organs. Harmonious interactions of the five organs promote normal metastasis and cell divisions, and active cell divisions assist in maintaining good metabolism, while conflict and contradiction in the functions of the five organs interrupt solute delivery from blood vessels to cells [24]. Insufficient delivery disturbs normal cell divisions, which is a cause of tumor cell growth.

Overdose of chemical medicine interrupts the balanced functions of the five organs and physical condition gets worse because of abnormal metabolism. Thus, it is required to control of the amount of medicines injected maintaining the harmonious balance of five organs during treatment. Due to complex interactions among the five organs, the present investigation is restricted to the variations of tumor cell density when medicines are injected into two organs, the heart and liver. Since the functions of liver are to purify corrupted blood and assist heart functions, only co-existence interactions between heart and liver organs are taken into account for the simple analysis. The complex behaviors of tumor based on the co-existence interactions among five organs will be followed. Mathematical formulations are deployed based on the following assumptions: (i) Two different medicines $I_1$ and $I_2$ are injected into heart and liver to recover the functions of each organ, respectively. (ii) A compartmental analysis is developed for solute transport from liver to blood vessels and no side effect appears in any compartment. (iii) Two factors, random motility and chemotaxis in response to TAF gradients, are considered in tumor cell motion.

2. Mathematical Modelling

2.1. Tumor Cell Mobility

Owing to the facts that the liver purifies corrupted blood and blood flows into heart, liver state deeply influences the performance of heart functions. The injection of substances such as medicine into the liver recovers liver function, and recovered liver activates heart functions. The following assumptions: (i) substances $L$ and $H$ injected into liver and heart, respectively, (ii) heart generates a solute $M$ to interfere with the development of TAF in the blood vessels is applied for mathematical modeling (see Figure 1). Compartmental analysis adopts for the concentration of chemotaxis caused by injected substances $L$ and $H$ under the assumption that the loss of mass is proportional to mass itself. The mathematical equations for the compartmental analysis are

$$\frac{d\left[L\right]}{d\tau }=I_2-k_2\left[L\right]$$

(1)

$$\frac{d\left[H\right]}{d\tau }=I_1-k_1\left[H\right] + k_2\left[L\right]$$

(2)

$$\frac{d\left[M\right]}{d\tau }=k_1\left[H\right]$$

(3)

The rectangular bracket in the above equations implies mass or concentration of each substance, $\tau$ time variable, and $k_1$ and $k_2$ reaction rate coefficients.

A tumor is centered on (0,0) with radius $\frac{R}{10}$ and a circle of radius $R$ is chosen for the domain. Let $c$ be the concentration of TAF in the domain, then the growth of TAF can be modeled by [25]

$$\frac{\partial c}{\partial \tau }=\nabla · \left(D_c\nabla c\right) + g\left(M,c\right)$$

(4)

The $D_c$ is the TAF diffusion coefficient, and $g$ is the rate of decay of TAF according to the disturbance of $M$. Due to the assumption that the chemotaxis solute in the blood vessels interferes with the development of TAF into the surrounding tissue, the function $g\left(M, c\right) =-Mc$ and the boundary condition $c\left(\tau ,R\right) =0$ are selected for mathematical approach. In the motion of tumor cells, random motility and chemotaxis in response to TAF gradients are essential factors to be considered. Let $\varnothing$ be the tumor cell’s density per unit area in the domain, then the variation of tumor cells can be expressed with

$$\frac{\partial \varnothing }{\partial \tau }=\nabla · \left(D_{\varnothing }\nabla \varnothing \right) -\nabla ·\left(\mu \varnothing \nabla c\right)$$

(5)

The first term in the right side of Equation (5) represents diffusion motility of cells and the second term is chemotaxis due to TAF gradients. The $D_{\varnothing }$ is the diffusion coefficient of cell and $\mu >0$ is the chemotactic coefficient. No flux condition is assumed on the boundary of domain, which yields the boundaries of the form

$$\overrightarrow{n}· \left(-D_{\varnothing }\nabla \varnothing \right) + \mu \varnothing \nabla c)=0$$

(6)

The $\overrightarrow{n}$ is an appropriate outward unit normal vector on the boundary of domain. With the consideration of the initial concentration of injected medicine

$$I_1\left(0\right) =1-e^{t_1}$$

$$I_2\left(0\right) =1-e^{t_2},$$

(7)

Medicine functions are represented with

$$I_1\left(\mathsf{\tau }\right) =1-e^{\tau -t_1}$$

$$I_2\left(\tau \right) =1-e^{\tau -t_2}$$

(8)

Note that each medicine affects until time level $t_1$ and $t_2$, respectively. The solutions for the differential Equations (1)–(3) then are the followings:

$$L\left(\tau \right) =\frac{1}{k_2}-\frac{e^{\tau -t_2}}{1+k_2}+c_1{e}^{-k_2\tau }$$

(9)

$$H\left(\tau \right) =\frac{2}{k_1}-\frac{e^{\tau -t_1}}{1+k_1}+c_2{e}^{-k_1\tau }+\frac{c_1{k}_2{e}^{-k_2\tau }}{k_1-k_2}-\frac{k_2{e}^{\tau -t_2}}{\left(1+k_1\right)\left(1+k_2\right)}$$

(10)

$$M\left(\tau \right) =\frac{2\tau }{k_1}-\frac{e^{\tau -t_1}}{1+k_1}-\frac{c_2{e}^{-k_1\tau }}{k_1}-\frac{c_1{e}^{-k_2\tau }}{k_1-k_2}-\frac{k_2{e}^{\tau -t_2}}{\left(1+k_1\right)\left(1+k_2\right)}+c_3$$

(11)

The initial conditions $L\left(0\right) =0$, $H\left(0\right) =0$, and $M\left(0\right) =0$ yield the values of integral constants $c_1$, $c_2$, and $c_2$;

$$c_1=-\frac{1}{k_1}+\frac{e^{-t_2}}{1+k_2}$$

$$c_2=-\frac{2}{k_2}+\frac{e^{-t_1}}{1+k_1}+\frac{1}{k_1{k}_2-k_2^2}+e^{-t_2}\left(\frac{1}{\left(1+k_2\right)\left(k_2-k_1\right)}+\frac{k_2}{1+k_1+k_2+k_1{k}_2}\right)$$

$$c_3=\frac{e^{-t_1}}{k_1}+\frac{e^{-t_2}\left(-1+k_1+k_1{k}_2-k_2^2\right)}{k_1\left(k_1-k_2\right)\left(1+k_2\right)}+\frac{k_1-k_1^2-2k_1{k}_2+2k_2^2}{k_1^2\left(k_1-k_2\right)k_2}$$

For the simplicity of TAF concentration, mobility the movement of TAF is assumed to be radially symmetric and the diffusion coefficient $D_c$ is constant. By the setting the $r=\tilde{r}\sqrt{D_c}$ and dropping the tilde for clarity Equation (4) is expressed as

$$\frac{\partial c}{\partial \tau }={\nabla }^{2}c-M\left(\tau \right)c$$

(12)

The separation of variables for Equation (12) by proposing a solution of the form

$$c=T\left(\tau \right)R\left(r\right)$$

Provides

$$\frac{1}{r}\frac{R^{\prime }\left(r\right)}{R\left(r\right)}+\frac{R^{″}\left(r\right)}{R\left(r\right)}=M\left(\tau \right)+\frac{T^{\prime }\left(\tau \right)}{T\left(\tau \right)}$$

(13)

Since the right-hand side of Equation (13) has only one dependent variable $\tau$ and the other terms has no $\tau$ dependence, the right-hand side is assumed to be equal to some unknown constant, that is,

$$M\left(\tau \right)+\frac{T^{\prime }\left(\tau \right)}{T\left(\tau \right)}=-\gamma$$

(14)

Then the following ordinary differential equation appears

$$r\frac{d}{dr}\left(r\frac{d{R}^{\prime }\left(r\right)}{dr}\right)+ r^2\gamma R\left(r\right) =0$$

(15)

The solutions of the ordinary differential Equations (14) and (15) are

$$\begin{gathered} T\left(\tau \right) =c_4{e}^{-\gamma \tau -\int M\left(\tau \right)d\tau }, \\ R\left(r\right) =c_5{J}_0\left(\sqrt{\gamma }r\right)+c_6{Y}_0\left(\sqrt{\gamma }r\right), \end{gathered}$$

(16)

where {$J_{\ast }, Y_{\ast }\}$ are Bessel functions of order $\ast$. Since $Y_{\ast \left(\sqrt{\gamma }r\right)}$ is unbounded at $r=0$, the setting $c_6=0$ removes an infinite value in the solution. Thus, the general solution of Equation (12) is, by superposition

$$c\left(\tau , r\right) =c_7{J}_0\left(\sqrt{\gamma }r\right)e^{-\gamma \tau -\int M\left(\tau \right)d\tau }$$

(17)

To determine the integral constant $c_7$ in the solution (17), an initial TAF concentration field is adopted of the form

$$c\left(0, r\right) = \left\{\begin{array}{c}1, 0\le r\le \frac{R}{10}\\ J_0\left(\sqrt{\gamma }r\right), \frac{R}{10}\le r\le R \end{array}\right.$$

The value $\gamma ={\left(\frac{2.40483}{R}\right)}^2$ is determined to satisfy the boundary condition $c\left(\tau , R\right) =0$ and the integral coefficient $c_7=\frac{1}{J_0\left(\sqrt{\gamma }\frac{R}{10}\right)}$ is obtained under the consideration of the continuity of TAF concentration at $r=\frac{R}{10}$.

2.2. Discretization for Finite Volume Method

The diffusion coefficient $D_{\varnothing }$ is assumed to be constant and let $r=\tilde{r}\sqrt{D_c}$. Equation (5) in polar coordinates, after dropping the tilde, can be expressed as

$$\frac{d\varnothing }{d\tau }=D\left[\frac{{\partial }^2\varnothing }{\partial r^2}+\frac{1}{r}\frac{\partial \varnothing }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2\varnothing }{\partial {\theta }^2}\right]-{\mu }_0\left[\frac{\partial c}{\partial r}\frac{\partial \varnothing }{\partial r}+\frac{1}{r^2}\frac{\partial c}{\partial \theta }\frac{\partial \varnothing }{\partial \theta }\right]-{\mu }_0\varnothing {\nabla }^{2}c,$$

(18)

where $D=\frac{D_{\varnothing }}{D_c}$ and ${\mu }_0=\frac{\chi }{D_c}$. The assumption of radial symmetricity for the movement of

TAF in Equation (18) is reduced to

$$\frac{d\varnothing }{d\tau }=D\left[\frac{{\partial }^2\varnothing }{\partial r^2}+\frac{1}{r}\frac{\partial \varnothing }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2\varnothing }{\partial {\theta }^2}\right]-{\mu }_0\frac{\partial c}{\partial r}\frac{\partial \varnothing }{\partial r}-{\mu }_0\varnothing {\nabla }^{2}c$$

(19)

Since Equation (19) is too involved to obtain the analytical solution, a finite volume approach is executed to search approximated solutions [26]. The typical arrangement of the finite volume cells on the polar grid is presented in Figure 2.

After the multiplication of $\mathrm{r}$ on both sides in Equation (19), the integration over the control volume and over a time interval from $\tau$ to $\tau +∆\tau$ provides

$$\begin{gathered} {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}r\frac{d\varnothing }{d\tau }dVd\tau ={{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}D\left[\frac{\partial }{\partial r}\left(r\frac{\partial \varnothing }{\partial r}\right)+\frac{1}{r}\frac{{\partial }^2\varnothing }{\partial {\theta }^2}\right]dVd\tau \\ -{{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}{\mu }_0\left[r\frac{\partial c}{\partial r}\frac{\partial \varnothing }{\partial r}\right]dVdr-{{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}{\mu }_{0}r\varnothing {\nabla }^{2}c dVd\tau \end{gathered}$$

(20)

Using the first order backwards differencing in time and second order central differencing in space, the following are derived:

$$\begin{gathered} {{\displaystyle \int }}_{CV}^{}\frac{\partial }{\partial r}\left(r\frac{\partial \varnothing }{\partial r}\right)dV=\int \left[{\left(r\frac{\partial \varnothing }{\partial r}\right)}_n-{\left(r\frac{\partial \varnothing }{\partial r}\right)}_s\right]d\theta =r_n∆{\theta }_n{\left(\frac{\partial \varnothing }{\partial r}\right)}_n-r_s∆{\theta }_s{\left(\frac{\partial \varnothing }{\partial r}\right)}_s \\ =r_n∆{\theta }_n\left[\frac{{\varnothing }_N-{\varnothing }_P}{{\delta }_{r_{PN}}}\right]-r_s∆{\theta }_s\left[\frac{{\varnothing }_P-{\varnothing }_S}{{\delta }_{r_{SP}}}\right] \\ =\frac{r_n∆{\theta }_n}{{\delta }_{r_{PN}}}{\varnothing }_N+\frac{r_s∆{\theta }_s}{{\delta }_{r_{SP}}}{\varnothing }_S-(\frac{r_n∆{\theta }_n}{{\delta }_{r_{PN}}}+\frac{r_s∆{\theta }_s}{{\delta }_{r_{SP}}}){\varnothing }_P, \\ {{\displaystyle \int }}_{CV}^{}\frac{1}{r}\frac{\partial }{\partial \theta }\left(\frac{\partial \varnothing }{\partial \theta }\right)dV=\frac{1}{r_w}∆r_w{\left(\frac{\partial \varnothing }{\partial \theta }\right)}_w-\frac{1}{r_e}∆r_e{\left(\frac{\partial \varnothing }{\partial \theta }\right)}_e \\ =\frac{1}{r_w}\frac{∆r_w}{{\delta }_{{\theta }_{WP}}}{\varnothing }_W+\frac{1}{r_e}\frac{∆r_e}{{\delta }_{{\theta }_{PE}}}{\varnothing }_E-(\frac{1}{r_e}\frac{∆r_e}{{\delta }_{{\theta }_{PE}}}+\frac{1}{r_w}\frac{∆r_w}{{\delta }_{{\theta }_{WP}}}){\varnothing }_P, \\ {{\displaystyle \int }}_{CV}^{}{\mu }_0\frac{\partial c}{\partial r}r\frac{\partial \varnothing }{\partial r}dV={\mu }_0\left[∆r∆\theta {\left(\frac{\partial c}{\partial r}r\right)}_P\frac{1}{{\delta }_{r_{NS}}}{\varnothing }_N-∆r∆\theta {\left(\frac{\partial c}{\partial r}r\right)}_P\frac{1}{{\delta }_{r_{NS}}}{\varnothing }_S\right], \\ {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}r\frac{\partial \varnothing }{\partial \tau }dVd\tau =r_P\left({\varnothing }_P-{\varnothing }_P^0\right)∆\theta ∆r, \\ {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}\frac{\partial }{\partial r}\left(r\frac{\partial \varnothing }{\partial r}\right)dVd\tau =\frac{r_n∆{\theta }_n}{{\delta }_{r_{PN}}}∆\tau {\varnothing }_N+\frac{r_s∆{\theta }_s}{{\delta }_{r_{SP}}}∆\tau {\varnothing }_S-\left(\frac{r_n∆{\theta }_n}{{\delta }_{r_{PN}}}+\frac{r_s∆{\theta }_s}{{\delta }_{r_{SP}}}\right)∆\tau {\varnothing }_P, \\ {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}\frac{1}{r}\frac{\partial }{\partial \theta }\left(\frac{\partial \varnothing }{\partial \theta }\right)dVd\tau =\frac{∆r_w}{r_w}\frac{∆\tau }{{\delta }_{{\theta }_{WP}}}{\varnothing }_W+\frac{∆r_e}{r_e}\frac{∆\tau }{{\delta }_{{\theta }_{PE}}}{\varnothing }_E-\left(\frac{∆r_e}{r_e}\frac{∆\tau }{{\delta }_{{\theta }_{PE}}}+\frac{∆r_w}{r_w}\frac{∆\tau }{{\delta }_{{\theta }_{WP}}}\right){\varnothing }_P, \\ {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}{\mu }_0\frac{\partial c}{\partial r}r\frac{\partial \varnothing }{\partial r}dVd\tau ={\mu }_0∆\theta ∆r∆\tau {\left(\frac{\partial c}{\partial r}r\right)}_P{\varnothing }_N-{\mu }_0∆\theta ∆r∆\tau {\left(\frac{\partial c}{\partial r}r\right)}_P{\varnothing }_S, \\ {{\displaystyle \int }}_{\tau }^{\tau +∆\tau }{{\displaystyle \int }}_{CV}^{}{\mu }_{0}r\varnothing {\nabla }^{2}cdVd\tau ={\mu }_0∆\mathsf{\theta }∆\mathrm{r}∆\mathsf{\tau }{\left(r{\nabla }^{2}c\right)}_P{\varnothing }_P. \end{gathered}$$

The discrete equation, thus, can be written as

$$a_p{\varnothing }_p= \left(a_N+a_{\mu N}\right){\varnothing }_N+\left(a_S+a_{\mu s}\right){\varnothing }_s+a_E{\varnothing }_E+a_W{\varnothing }_W+b,$$

where

$$\begin{gathered} a_P^0=\frac{r_P∆r∆\theta }{∆\tau },a_N=\frac{D{r}_n∆{\theta }_n}{{\delta }_{r_n}},a_S=\frac{D{r}_s∆{\theta }_s}{{\delta }_{r_s}},a_E=\frac{D∆r_e}{r_e{\delta }_{{\theta }_e}}, \\ {a}_W=\frac{D∆r_w}{r_w{\delta }_{{\theta }_w}}, a_{\mu N}=-\frac{{\mu }_0∆\theta ∆r}{{\delta }_{r_{NS}}}{\left(\frac{\partial c}{\partial r}r\right)}_{P,}{a}_{\mu S}=\frac{{\mu }_0∆\theta ∆r}{{\delta }_{r_{NS}}}{\left(\frac{\partial c}{\partial r}r\right)}_{P,} \\ {a}_{\mu P}^0={\mu }_0∆r∆\theta {\left(r{\nabla }^{2}c\right)}_P,b=a_P^0{\varnothing }_P^0, a_P=a_P^0+a_{\mu P}^0+a_N+a_S+a_E+a_W. \end{gathered}$$

The ${\varnothing }_P^0$ represents the density of cell at time $\tau$ and the density at time $\tau +∆\tau$ is not superscripted. Note that ${\delta }_{r_{NS}}=2{\delta }_r$, ${\delta }_{r_{PN}}={\delta }_{r_{SP}}={\delta }_r$, ${\delta }_{{\theta }_{EW}}=2{\delta }_{\theta }$, ${\delta }_{{\theta }_{PW}}={\delta }_{{\theta }_{PE}}={\delta }_{\theta }$.

The boundary condition (Equation (6)) is written by

$$\begin{gathered} \overrightarrow{n}∆\left[-\frac{D_{\varnothing }}{\sqrt{D_c}}\nabla \varnothing +\frac{\chi }{\sqrt{D_c}}\varnothing \nabla \mathrm{c}\right]=\overrightarrow{n}∆\left[-\frac{D_{\varnothing }}{\sqrt{D_c}}\left(\frac{\partial \varnothing }{\partial r}{u}_r+\frac{1}{r}\frac{\partial \varnothing }{\partial \theta }{u}_{\theta }\right)+\frac{\chi }{\sqrt{D_c}}\varnothing \left(\frac{\partial c}{\partial r}{u}_r+\frac{1}{r}\frac{\partial c}{\partial \theta }{u}_{\theta }\right)\right] \\ =u_r∆\left[\left(-\frac{D_{\varnothing }}{\sqrt{D_c}}\frac{\partial \varnothing }{\partial r}+\frac{\chi }{\sqrt{D_c}}\varnothing \frac{\partial c}{\partial r}\right)u_r+\left(-\frac{D_{\varnothing }}{\sqrt{D_c}}\frac{1}{r}\frac{\partial \varnothing }{\partial \theta }\right)u_{\theta }\right] \\ =-\frac{D_{\varphi }}{\sqrt{D_c}}\frac{\partial \varnothing }{\partial r}+\frac{\chi }{\sqrt{D_c}}\varnothing \frac{\partial c}{\partial r}=0, \end{gathered}$$

where $u_r$ and $u_{\theta }$ are directional unit vectors. Thus, the finite volume discretized form of the boundary condition is

$$-\frac{D_{\varphi }}{\sqrt{D_c}}\frac{\partial \varnothing }{\partial r}+\frac{\chi }{\sqrt{D_c}}\varnothing \frac{\partial c}{\partial r}=\frac{D_{\varnothing }}{\sqrt{D_c}}\frac{1}{{\delta }_{r_{NS}}}{\varnothing }_N+\frac{D_{\varnothing }}{\sqrt{D_c}}\frac{1}{{\delta }_{r_{NS}}}{\varnothing }_S+\frac{\chi }{\sqrt{D_c}}{\left(\frac{\partial c}{\partial r}\right)}_P{\varnothing }_P=0$$

2.3. Validation of Numerical Results

For validation of numerical results, the results of Kashkooli et al. [27] are compared with the outcomes of the present study (see Figure 3). SR in Figure 3 represents survival rate. Kashkooli et al. considered a three-stage drug delivery system containing a primary nanoparticle, a secondary nanoparticle, and the chemotherapy on a real image of a vascularized tumor. As shown in Figure 3, the survival rate of tumor cells is close to 10% in the present study, while Kashkooli et al. displayed around 30% at an early stage. Its value is 20% different and the difference is due to the different transport method of the drug, boundary conditions, tumor geometry, and capillary domain. The decreasing patterns are similar except for the oscillation of [27], which is due to the growth rate.

3. Results and Discussions

Based on the discretization of the finite volume method in Section 2.2, approximated solutions for the tumor cells density $\varnothing$ are obtained combing the exact solutions of medicine $M$ in blood vessels and the TAF concentration $c.$ An initial value ${\varnothing }_0=1000$ is taken for the tumor cells density, and the values of parameters in Reference [3] displayed in

Table 1. Values of parameters.

$\chi$	${\mathit{D}}_{\mathit{c}}$	${\mathit{D}}_{\varnothing }$
$2.6\times {10}^3 {\mathrm{cm}}^2 {\mathrm{s}}^{-1}{ \mathrm{M}}^{-1}$	$2.9\times {10}^{-7} {\mathrm{cm}}^2{ \mathrm{s}}^{-1}$	$1.0\times {10}^{-10} {\mathrm{cm}}^{2 }{\mathrm{s}}^{-1}$

are applied for numerical solutions. The step sizes of each variable are $∆\mathsf{\theta }=\frac{\pi }{10}$, $∆r=\frac{R}{20}$, $∆T=1$, and the value $R=\frac{4}{\sqrt{D_c}}$ is used.

The values of parameters are adopted because the values are obtained in the investigation of the formation of blood vessels from a pre-existing vasculature. The representative values $\frac{k_1}{k_2}=2, 3, 4, 5$ are chosen to investigate the ratio variation effects of flow rates $k_1$ and $k_2$. The TAF concentrations are described in Figure 4 for the time values $t_1=t_2=10$ and the ratio value $\frac{k_1}{k_2}=2$. Exponential decay develops in TAF concentration with the increase of $r$, which means that tumor angiogenic factors are condensed around the center. By that behavior, TAF reacts sensitively to the injected medicines and TAF concentration decreases exponentially as the time value $\tau$ grows explaining that administered meditation is mostly fatal to tumor angiogenic factors. After the time value $\tau \approx 10$, the TAF concentration is close to zero.

Figure 5 exhibits the volume of medicine $M$ in blood vessels for the time values $t_1=t_2=10$. The medicine solute amount grows linearly until around $\tau =10$ and decreases exponentially after the time value $\tau \approx 10$. The largest medicine concentration appears near the time $\tau =10$ and is close to zero near the time $\tau =13$. The concentration of $M$ decreases as the ratio $\frac{k_1}{k_2}$ increases. The result explains that the medicine solute amount shrinks with the reduction of medicine flow rate into liver when the flow rate into heart is fixed, whereas the medicine solute amount expands as the medicine flow rate into heart decreases for the fixed flow rate into liver. The rate of decline is getting smaller with the increase of the ratio $\frac{k_1}{k_2}$, which implies that similar flow rate of each injected medicine may produce the greatest amount of medicine in blood vessels.

The influences of the variation of ratio $\frac{k_1}{k_2}$ to TAF concentration are presented through Figure 6 and Figure 7, which describe the responses of tumor cells to medicines. The description according to the time value variable $\tau$ is executed at $r=\frac{R}{2}$ when $t_1=t_2=10$ (see Figure 6). TAF concentration increases with the growth of the value of ratio $\frac{k_1}{k_2}$. The result interprets that the increase of the medicine flow rate into liver yields a profound effect to tumor angiogenic factors when the flow rate into heart is fixed, whereas TAF reacts less to the flow rate increment into heart with a fixed flow rate into liver. When the flow rate of each injected medicine is similar, TAF shows the smallest concentration. The influences of the ratio $\frac{k_1}{k_2}$ variation to tumor angiogenic factors along the radius are displayed in Figure 7. The fixed time value $\tau =5$,$t_1=t_2=10$ are chosen. TAF concentration decreases along the radius and a parabolic shape appears in the diminution description. A similar phenomenon occurs in TAF concentration expression, that is, TAF concentration diminishes along the radius with a decay of the ratio $\frac{k_1}{k_2}$ value.

Approximate solutions for the tumor cells density $\varnothing$ are obtained based on Equations (18) and (19), and Figure 8 and Figure 9 depict the tumor cells density $\varnothing$ for the values $t_1=10$, $t_2=10$, and $\frac{k_1}{k_2}=2$. The representative radius values r$=\frac{R}{10}, \frac{R}{5},$ and $\frac{2R}{5}$ are chosen to investigate the density of tumor cells with time variable $\tau$ (see Figure 8) and the difference of initial amount is due to the distance from the center of tumor cell density. Tumor cells react sensitively as soon as medicines are injected and the mobility of tumor cells are disturbed immediately. At the time $\tau =1$, tumor cells density $\varnothing$ decreases to around 20% of the initial amount for the value r$=\frac{R}{10}$, while the density $\varnothing$ is close to zero for the value r$=\frac{2R}{5}$. After the rapid change at the moment of medicine injection, the variation of tumor cells density is trivial until around the time $\tau =20$, and the density $\varnothing$ diminishes with the increment of time value $\tau$. The result evinces that the tumor cells show hypersensitivity reaction at the moment of medicine injection. Figure 9 depicts tumor cell density in respect of radius $r$ for the representative time values $\tau =0, 1, 20,$ and $22$. Initial tumor cell density ${\varnothing }_0=1000$ and angle value $\theta =0$ are taken for the numerical approach. Tumor cell density decreases exponentially at the time value $\tau =0$, which implies that the reaction of tumor cells is highly sensitive to the injected medicines. Owing to the repellence of tumor cells to the injected medicines, the decrease of tumor cell density is developed in a fluctuating manner along the radius. Similar aspects appeared in the tumor cell density variation for time values $\tau =1, 20,$ and $22$, and the tumor cells density approaches zero near the boundary of domain.

The variation of tumor cell density is analyzed through a mathematical approach considering the recovery of liver and heart functions by drug injection. The functional recoveries of liver and heart produce medicine in blood vessel to disable the activity of TAF. The control of TAF activity influences the motility and chemotaxis of tumor cells and the proliferation is restricted. The present work suggests a strong method to reduce the tumor cell density at an early stage of medicine injection, and the survival rate of tumor cells is lower in comparison to previous research [27]. In addition, the flow rate from each organ plays an important role in determining the medicine concentration in blood vessel even though only two organs, liver and heart, are taken into account. The extension to the functional recoveries of five organs: kidney, liver, heart, spleen, and lung, evokes a great curiosity in the study of tumor cell density characteristics. Even so, the consideration of complicated reactions among the five organs in deciding medicine concentration in blood vessel may be an interesting and challenging work.

4. Conclusions

In the present study, a mathematical model of drug supplies is proposed to investigate the characteristics of tumor cell density. A compartmental analysis is developed for solute transport from liver and heart into blood vessel ignoring any side effect in any compartment. Random motility and chemotaxis in response to TAF gradients are considered in tumor cell motion. The mathematical approach has been developed to interpret the variations of solute medicine amount in blood vessel and TAF concentration with respect to the flow rate of injected drugs into liver and heart and describe the alteration of tumor dell density along the time and radius variables. The main results obtained are as follows: (i) the concentration of $M$ decreases with the growth of the ratio $\frac{k_1}{k_2}$, (ii) TAF concentration increases with the growth of the value of ratio $\frac{k_1}{k_2}$ and TAF shows the smallest concentration when the flow rate of each injected medicine is similar, (iii) tumor cells react highly sensitive as soon as medicine supplies and tumor cell density is decreased drastically at the moment of medicine injection, (iv) tumor cell density decreases exponentially at an early stage and the decrease of tumor cells density is developed in a fluctuating manner along the radius.

Throughout the results obtained the flow rate of medicine injected into representative organ is a crucial factor to control the concentration of $M$ and TAF concentration and determine the amount of medicine in blood vessels. Since similar flow rate of each injected medicine produces the greatest amount of medicine in blood vessels, harmonious balanced function between liver and heart may contribute to control the proliferation of tumor cells. Tumor cells react immediately after medication and the density decreases exponentially at the time of drug injection. However, the tumor cell density is diminished in a fluctuating manner along the radius implying that a sufficient period is required for complete treatment. The presented mathematical approach has the potential for the profound analysis of the variations of solute medicine amount in blood vessels, TAF concentration, and the alteration of tumor cell density according to the recovery of liver and heart functions. The mathematical approach may be applicable in the investigation of tumor cells’ behavior on the basis of complex interaction among the five represented organs: kidney, liver, heart, spleen, and lung.