Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 1, pp. 271-279, Warsaw 2014 MATHEMATICAL MODELING OF FLUID FLOW IN BRAIN TUMOR Daniel N. Riahi, Ranadhir Roy University of Texas-Pan American, Department of Mathematics, Edinburg, USA e-mail: driahi@utpa.edu We consider the problem of fluid flow in a brain tumor.We develop a mathematical model for the one-dimensional fluid flow in a spherical tumor where the spatial variations of the interstitial velocity, interstitial pressure and the drug concentration within the tumor are only with respect to the radial distance from the center of the tumor. The interstitial ve- locity in the radial direction and the interstitial pressure are determined analytically, while the radial variations of two investigated drug concentrations were determined numerically. We calculated these quantities in the tumor, in a corresponding normal tissue and for the concentrations also in the cavity that can exist after the tumor is removed.We determine, in particular, theway the interstitial pressure and velocity vary, which agreeswith the expe- riments, as well as the way one drug concentration changes in the presence or absence of a second drug concentration within the tumor. We find that the amount of drug delivery in the tumor can be enhanced in the presence of another drug in the tumor, while the ratio of the amount of one drug in the tumor to its amount in the normal tissue can be reduced in the presence of the second drug in the tumor and the tissue. Key words: tumor, brain tumor, spherical tumor, drug concentration, fluid flow 1. Introduction It is well documented that cancer is the second leading cause of death in the North American Continent (Jain, 2005). There have been a number of studies of tumors in general (Baxter and Jain, 1989, 1990, 1991;Wang andLi, 1998; Zhao et al., 2007; Soltani andChen, 2011) and brain tumors in particular (Saltzman andRadomsky, 1991; Tan et al., 2003; Teo et al., 2005). Despite many experimental or computational investigations on the tumors that have been carried out in the past including those referred to here, there is still little information available about the mechanism that operates for drug interactions when more than one drug concentrations are transported in the patient’s body to effectively reduce the negative effects of malignant tumors. It is known that the most effective way to cure cancer due to a brain tumor is a successful removal of the tumor by surgery and then delivering in an efficientway certain anticancer drugs. However, the presence of drugs in the patient’s brain can lead to negative side effects and the interactions of drugsmay negatively affect the efficiency of the drug delivery and its usefulness. In addition, anticancer drugs also can affect the normal tissues in the brain. So an important studycanbe tounderstandthe effects of drug concentration in thebrain tumor in thepresenceof another anticancer drug, how the interaction of, say, two drugs can affect the overall usefulness of each of such drugs and the negative effects that anticancer drugs can have on the normal tissues as well as on the cavity after the tumor removal by surgery. In the present study, we develop a mathematical model of a brain tumor by approximately simplifying its shape into a sphere, and we first consider steady state features of the fluid flow within the tumor and in the corresponding spherical normal tissue and for the drug concentra- tions also in the cavity that can exist after the tumor is removed by surgery. Such steady state features of the flowquantities, which canbe considered as time average of the flowquantities, are 272 D.N. Riahi, R. Roy reasonable as the leading first approximation on the assumption that the rate of the growth of the tumor is sufficiently small. Next, we consider an approximate analysis that assume small and non-zero growth rate of the tumor.We determine important quantities such as interstitial velo- city and pressure as functions of distance from the center of tumor or as function of time, and, in addition, we determine the radial variations of drug concentration in the presence or absence of another drug within such systems. This is a first type of theoretical investigation to obtain information about the mechanism of the drug interactions that have not yet been investigated sufficiently by other investigators so far aswe referred to earlier in this section. Furthermore, our present investigation can help to improve understanding of the transportmechanisms in tumors which can consequently help to improve drug delivery schemes to tumors. We found some interesting results. In particular, we detected inter relations between these variations and the physical systems aswell as somemechanical cause for the decrease or increase of the drug effectiveness within the patient’s brain which can be useful as further guiding tools for the specialists and doctors to improve chances for the patient’s health recovery. We also found that increasing the value of the interstitial pressure can be due to the increased tumor size, which is in agreement with some experimental results (Raju et al., 2008). 2. Mathematical formulation We first provide a motivation for our present investigation in this paper. Let us first consider a spherical tumor with radius R∗(t∗) and a spherical coordinate system whose origin lies at the center of the tumor andwith radial r∗-axis directed positively outward from the origin. Here t∗ is time variable. We assume that growth rate of the tumor is sufficiently small (∂R∗/∂t∗ ≪ 1), so that leading order of the radius is constant. We consider the governing equations for the fluid flow in the tumor or the corresponding normal tissue or the cavity after tumor removal (Tan et al., 2003; Soltani and Chen, 2011), which are basically the Darcy-momentum equation, the mass continuity equation in the presence of sinks and sources, which is appropriate for the interstitial fluid (Baxter and Jain, 1989). The sources and sinks are based on Starling’s law and the fluid productions by cells due to metabolism (Tan et al., 2003) and the equations for the concentrations of two drugs. These equations and the corresponding boundary conditions are given below µ k u ∗ =∇P∗ ∇·u∗ = a1−a2P∗ ( ∂ ∂t∗ +u∗ ·∇ ) C1 =D1∇2C1−a3C1−a4C2 ( ∂ ∂t∗ +u∗ ·∇ ) C2 =D2∇2C2−a5C1−a6C2 ∂P∗ ∂r∗ =0 at r∗ =0 P∗ =PB∗ ∧ Ci =C∗Bi ∧ ∂Ci ∂r∗ =T∗i at r ∗ =R∗(t∗) i=1,2 (2.1) where u∗ is the interstitial fluid velocity vector, P∗ is the interstitial pressure, P∗B is a non-zero constant, µ is dynamic viscosity, k is permeability, a1 and a2 are constants whose expressions are given in Tan et al. (2003) with a1 due to hydraulic conductivity Kv of the macro-vascular wall, vascular pressure, osmotic reflection coefficient for the plasma protein, osmotic pressures of theplasmaand interstitial pressures and exchange area S/V of thebloodvessels perunit volume of tissues and a2 =KvS/V . In addition, Ci (i=1,2) are the concentrations of two drugs with diffusion coefficients Di,C ∗ Bi and T ∗ i are constants, and ai (i=3, . . . ,6) are constant coefficients of the source or sink terms in the concentration equations, which can bedue to presence of either Mathematical modeling of fluid flow in brain tumor 273 one or two drugs. As explained in Tan et al. (2003), such source and sink terms can be due to chemical elimination by drug degradation in the cavity ormetabolic reactions in the tumor and the normal tissue as well as by the drug gain from the blood capillaries in the tumor and tissue. The expression for the time-dependent radius of the tumor is considered to be in form R∗(t∗)=R0exp{ε[1− exp(−α∗t∗)]} (2.2) where ε is a non-dimensional small parameter (ε≪ 1) representing the order of magnitude of the relative growth rate of the tumor, R0 is the constant radius of tumor in the absence of its growth and α∗ is a constant of order one quantity per unit time whose value can affect the relative growth rate of the tumor. So, α∗ε represent the inverse of a time scale and provide the relative rate of growth of the tumor initially (t′ = 0). We have assumed that the tumor grows very slowly in time, so that (1/R∗)∂R∗/∂t∗ = o(ε)≪ 1. We consider the one-dimensional case of the flow system where the spatial variation is only along the radial direction and the velocity vector has only anon-zero component along the radial direction. We make governing system (2.1) and (2.2) dimensionless by using length scale R0, time scale R20/ν, where ν ≡ µ/ρ is kinematic viscosity and ρ is fluid density, velocity scale ν/R0 and pressure scale ν 2ρ/k, so that the dimensional forms of the variables andR are related to their non-dimensional counterpart by (r∗, t∗,u∗,P∗,R∗)= ( R0r, R20 ν t, ν R0 u, ν2ρ k P,R0R ) (2.3) It should be noted that our non-dimensional procedure makes the quantities of interest just dimensionless, so that the results can be applicable to a wide range of cases and biomedical problems in the subject area of this paper. We now apply a perturbation expansion in powers of small ε by assuming that the growth rate of increase of the tumor with respect to time is sufficiently small and write in the form (u,P,C1,C2,R)= [u0(r),P0(r),C10(r),C20(r),1] +ε[u1(r),P1(r),C11(r),C21(r),1][1−exp(−αt)]+o(ε2) (2.4) where α=α∗R20/ν is a non-dimensional constant. In the first part of the present study, we assume that the growth rate of the tumor is sufficiently small (ε≪ 1), so that to the leading order approximation, we determine the steady state features of the flow quantities. Thus, we consider (2.4) to the lowest order in ε and use it together with (2.3) in (2.1) and (2.2) to find the following simplified form of the system that will be investigated in this paper dP0 dr =u0 du0 dr = b1− b2P0 u0 dC10 dr =L1 [ 1 r2 d dr ( r2 dC10 dr )] − b3C10− b4C20 u0 dC20 dr =L2 [ 1 r2 d dr ( r2 dC20 dr )] − b5C10− b6C20 dP0 dr =0 at r=0 P0 =PB ∧ Ci0 =CBi ∧ dCi0 dr =Ti at r=1 i=1,2 (2.5) where bi = aiR 2 0/ν (i=1,3,4,5,6), b2 = a2νR 2 0/k, PB =P ′ Bk/(ρν 2) and Ti =T ′ iR0. In the second part of our study, we take into account some approximate analysis of (2.4) by keeping terms in the order ε (ε ≪ 1) for the radial velocity and the interstitial pressure, and 274 D.N. Riahi, R. Roy thus determine the solutions for these quantities, which can provide us the effect of the tumor growth on the flow. Thus, using (2.4) in the order ε in (2.1)-(2.3), we find the following system to be solved. The results for the approximate unsteady state will be studied in this paper dP1 dr =u1 du1 dr =−b2P1 dP1 dr =0 at r=0 P1 =− dP0 dr at r=1 (2.6) 3. Results and discussion Using (2.5)1,2 and (2.5)5,6, we find the following analytical solutions for the velocity andpressure P0(r)= b1 b2 + ( PB − b1 b2 )cos( √ b2r) cos √ b2 u0(r)=− √ b2 ( PB − b1 b2 )sin( √ b2r) cos √ b2 (3.1) Using (3.1)2 in (2.5)3,4, we solve numerically these second order differential equations subjec- ted to the boundary conditions given in (2.5)6 by an efficientRunge-Kutta fourth ordermethod. We calculated the main quantities, which are the velocity, pressure and the concentrations of two drugs, by using the numerical values for the non-dimensional constant coefficients bi (i=1, . . . ,6), and for the non-dimensional parameters Li and the boundaryvalues PB, CBi and Ti (i = 1,2). These numerical values were chosen based on the dimensional values of the cor- responding quantities, which were collected mainly from the relevant literature on biomedical applications (Jain and Baxter, 1988; Tan et al., 2003; Soltani and Chen, 2011). In particular, the values of the coefficients b1 − b3 and the diffusivity parameter L1 used in this paper are directly based on dimensional values given in Tan et al. (2003). The non-dimensional values of the constant coefficients bi (i = 1, . . . ,6), the boundary values PB, CBi and Ti (i = 1,2) and the diffusion parameters Li (i=1,2) are given in Table 1. Table 1.Non-dimensional values of the quantities that are used in the calculation Quantities Tumor values Normal tissue Cavity values values b1 −4.436064 ·10−11 −5.84496 ·10−12 b2 8.44 ·10−11 1.08 ·10−11 b3 3.32 ·10−4 0.424 ·10−4 0.424 ·10−4 b4 0.03324 ·10−4 0.00424 ·10−4 0.00424 ·10−4 b5 3.32 ·10−4 0.424 ·10−4 0.424 ·10−4 b6 0.03324 ·10−4 0.00424 ·10−4 0.00424 ·10−4 CB1 0.1 ·10−3 0.1 ·10−3 0.1 ·10−3 CB2 0.1 ·10−3 0.1 ·10−3 0.1 ·10−3 L1 1.4175 ·10−5 0.525 ·10−6 0.3003 ·10−5 L2 1.4175 ·10−8 0.525 ·10−8 0.3003 ·10−8 PB 1.84 ·10−3 1.84 ·10−3 1.84 ·10−3 T1 0.1 ·10−2 0.1 ·10−2 0.1 ·10−2 T2 0.1 ·10−2 0.1 ·10−2 0.1 ·10−2 For the two drug concentrations we chose etanidazole, which is used for cancer patients to regulate the level of oxygen concentration in the tissue in order tomake radiotherapy applicable Mathematical modeling of fluid flow in brain tumor 275 to destroymalignant cells, andCisplatinwhich is a kindof chemotherapy drug.Figure 1 presents interstitial radial velocity of the flowversus the radial variable for both tumor andnormal tissue. As can been seen from this figure, u0 < 0 which is due to the fact that the flow enters into the tumor or tissue from the boundary surface. The velocity profile appears to be linear in the domain for r (0 < r < 1), and its magnitude for the tumor is higher than that for the tissue. It is also seen that |u0|→ 0 as r→ 0 which is reasonable since the high interstitial pressure at the center put the motion into rest there. Fig. 1. Radial velocity ×10−10 versus radial variable for tumor and normal tissue Figure 2 presents interstitial pressure versus the radial variable for both tumor and normal tissue. It can be seen from this figure that the pressure increases as the radius decreases leaving maximum pressure at the center or tumor or tissue. The value of the pressure in the tumor is higher than that in the tissue, and, in addition, the radial rate of change of the pressure in the tumor is higher than the corresponding one in the tissue. Our results for both interstitial pressure and velocity agree qualitatively with the available experimental results (Jain, 1987; Baxter and Jain, 1989; Boucher et al., 1990; Jain et al., 2007). Fig. 2. Instestitial pressure versus radial variable for tumor and normal tissue It should be noted that even though the velocity and pressure profiles are given in terms of trigonometry functions as can be seen from (2.6)1,2, the shape of the velocity profile in Fig. 1 is linear in contrast to one for the pressure shown in Fig. 2. But the linearity of the velocity profile is due to relatively very small values of the velocity that appears linear over the short radial distance from 1 to 0. In the absence of any drug interaction and in the presence of only one drug concentration, we did some calculations, some of which are presented in Figs. 3 and 4. These figures present respectively the concentrations of etanidazole and cisplatin versus the radial variable and for the cases of tumor, normal tissue and cavity.We can see from these figures that both concentrations of these drugs increase with decreasing radius, and the value of concentration in the tumor is less that in the tissue, and the value of the concentration in the tissue is less than that in the cavity. Similarly, the radial rate of change of the concentration is highest in the cavity and is the lowest in the tumor. For a given value of the radius, the concentration for Cisplatin is found to be higher than that for the etanidazole, which is due to smaller diffusion parameter L2 as 276 D.N. Riahi, R. Roy compared to the value of the diffusion parameter L1 (Landau and Lifshitz, 1987). The results for the radial variations of the concentrations indicate that the concentrations close to r = 1 have a smaller radial rate of increase as compared to their rate of increase near the center. We also calculated the drug concentration for different initial loadings and found that the values of the concentration increase with an increase in such loading. Fig. 3. Etanidazole concentration versus radial variable in the absence of the other drug Fig. 4. The same as in Fig. 3 but for concentration of cisplatin We also calculated the concentrations of drugs when both drugs are present, and some of our results for such a case are shown in Figs. 5 and 6. These figures present respectively the concentrations of etanidazole and cisplatin versus the radial variable for the cases of tumor, normal tissueandcavity. It canbe seen fromthesefigures thatboth concentrations again increase with the decreasing radial variable, and the radial rate of change of each concentration is highest in the cavity and lowest in the tumor. Similar to the case in the absence of drug interaction presented inFigs. 3 and 4, we can see fromFigs. 4 and 5 that the radial rate of change of each of these concentrations is higher close to the center, and the amount of concentration of cisplatin is again higher than that for etanidazole. Comparing the results in Figs. 5 and 6 to those in Figs. 3 and 4, we find that as a result of presence of both concentrations, the values of the concentrations increase either in the tumor system, in the tissue or in the cavity. This result indicates that in the actual drug delivery mechanism to the patients, these can be a case where providing a secondary drug to the patient helps improving the amount of the primary drug to the patient’s tumor, thereby improving the patient’s health conditions. It should be noted from Figs. 4 and 6 the inflexion type points in the concentration profiles for the cisplatin in the tissue. These points are due to the fact that the diffusivity parameter as given in Table 1 for the cisplatin has an intermediate value as compared to the corresponding values for the tumor, where the diffusivity parameter is larger and the profile appears nonlinear, and for the cavity, where the diffusivity parameter is smaller and the profile appears to be close to a linear profile. This result appears to be due to the behaviour of the cisplatin concentration and the value of its diffusivity parameter where its profile tends to move away from the linear form and approaches a nonlinear form, and hence, as a result, some inflexion points appear in its profile. Mathematical modeling of fluid flow in brain tumor 277 Fig. 5. The same as in Fig. 3 but in the presence of cisplatin Fig. 6. The same as in Fig. 4 but in the presence of etanidazole We should also refer to the efficacy of the drug delivery system (Tan et al., 2003) that for a given radius it can depend, in particular, on the ratio of the values of the concentration for the tumor to that for the normal tissue. This ratio referred to as Therapeutic Index (TI) in relevant biomedical literature (Tan et al., 2003). So, TI is ameasure of the efficiency of the drug delivery, so that a higher value of TI indicates that a higher amount of the drug delivers to the tumor than to the normal tissue. In our results, for the concentrations of etanidazole and cisplatin and the results presented in Figs. 3-6, we conclude that for either etanidazole or cisplatin, TI ismore in the presence of only one drug in the system. In order to determine some effect of the tumor growth rate on the flow, we first determine the solutions to (2.6), which are given below P1 = ( PB − b1 b2 ) √ b2 tan √ b2 cos( √ b2r) cos √ b2 u1 =−b2 ( PB − b1 b2 ) tan √ b2 sin( √ b2r) cos √ b2 (3.2) Using (3.1) and (3.2) in (2.4) and taking into account the results in Table 1, we generated data for interstitial radial velocity and pressure for both tumor and tissue versus time t for fixed values of r, α (α > 0) and ε (ε ≪ 1) as well as versus r for fixed values of t, α and ε. We found that the radial velocity for both tumor and tissue is negative and itsmagnitude decreases with the decreasing radial variable, and the magnitude of the radial velocity for the tumor is higher than that for the tissue. The interstitial pressure is found to be positive for both tumor and tissue and increases with the decreasing radial variable. At time increases, where the radius of the tumor increases with t as seen in (2.4), the interstitial radial velocity for both tumor and tissue increases inmagnitude. Similarly, the interstitial pressure increaseswith time as the tumor grows in time, which appears to agree, in particular, with some experimental results (Raju et al., 2008). 4. Conclusions We investigated the problem of flow in a brain tumor by developing a theoretical model for a spherical tumor or tissue under the assumption that the growth rate of the tumor is sufficiently 278 D.N. Riahi, R. Roy small. We, thus, first calculated the steady results for the main quantities such as interstitial pressure, unidirectional velocity and concentrations of two different drugs. We found that the interstitial pressure decreases with the increasing radius, while the magnitude of the velocity increases with the radial distance from the center of the tumor or tissue. These results appear to agree qualitatively with the available experimental results. We also calculated the values of concentrations of two drugs (etanidazole and cisplatin) in the presence or absence of the second drug and found, in particular, that the presence of both drugs can improve the amount of the drug in the flow system. This result also indicates that non-trivial cases can be investigated theoretically to discover ways that improve drug delivery to the patient and thereby improve health condition of such a patient. Our results about the influence of the tumor growth on the flow indicate that as time increases the magnitude of the interstitial radial velocity increases with the increasing radius of the tumor. The value of the interstitial fluid pressure is also found to increase with the increasing size of the tumor, which agrees with the experimental results. The present paper can be considered as the first part of a two-part investigation of our theoretical and modeling development of the flow in a spherical brain tumor or normal tissue. In the present part, we carried out investigation and calculations for the steady case for the flow and drug delivery as well as for the effect of a non-zero growth rate of the tumor or tissue on the flow, under the assumption that the growth rate of the tumor is sufficiently small. In the second part of our investigation that we plan to complete in the near future, we will present results for unsteady cases, where the growth rate of the tumor is not assumed to be sufficiently small, and, in particular, we will determine results for drug concentration in the tumor or tissue affected by the growth rate of the tumor and how such a growth rate can be controlled by special drug delivery mechanisms. In addition to the time-dependent extensions that we referred to in the previous Section, the theoreticalmodel developed and investigated here can be further extended and involved to apply to non-unidirectional tumors and tissues in terms of properties and geometrical configurations. Another extension can be prediction of the effect that the rate of drug delivery to the tumor site can have on the tumor growth for the actual operating and drug conditions in the medical cases of brain tumor patients. Acknowledgements This researchwas supported by FRC-UTPA grant during the academic year 2011-2012. References 1. Baxter L.T., Jain R.K., 1989, Transport of fluid and macromolecules in tumors. (I) Role of interstitial pressure and convection,Microvascular Research, 37, 77-104 2. Baxter L.T., Jain R.K., 1990, Transport of fluid and macromolecules in tumors. (II) Role of heterogeneous perfusion and lymphatics,Microvascular Research, 40, 246-263 3. Baxter L.T., Jain R.K., 1991, Transport of fluid and macromolecules in tumors. (III) Role of binding andmetabolism,Microvascular Research, 41, 5-23 4. Boucher Y., Baxter L.T., Jain R. K., 1990, Interstitial pressure gradients in tissue-isolated and subcutaneous tumors: Implication for therapy,Cancer Research, 50, 4478-4484 5. Jain R.K., 1987, Transport of molecules in the tumor interstitium: A review, Cancer Research, 47, 3039-3051 6. Jain R.K., 2005, Normalization of tumor vasculature: An emerging concept in antiangiogenic therapy, Science, 307, 58-62 7. Jain R.K., Baxter L.T., 1988, Mechanisms of heterogeneous distribution of monoclonal anti- bodies and other macromolecules in tumors: significance of elevated interstitial pressure, Cancer Research, 48, 7022-7032 Mathematical modeling of fluid flow in brain tumor 279 8. Jain R.K., Tong R.T., Munn L.L., 2007, Effect of vascular normalization by antiangiogenic therapy on interstitial hypertension, peritumor edema, and lymphaticmetastasis,Cancer Research, 67, 2729-2735 9. Landau L.D., Lifshitz E.M., 1987,Fluid Mechanics, 2nd edition, PergamonPress, Oxford 10. Raju B., Haug S.R., Ibrahim S.O., Hereraas K.J., 2008, High interstitial fluid pressure in rat tongue cancer is related to increased lymph vessel area, tumor size, invasiveness and decreased body weight, Journal of Oral Pathology and Medicine, 37, 3, 137-144 11. Saltzman W.M., Radomsky M.L., 1991, Drugs released from polymers: diffusion and elimina- tion in brain tissue,Chemical Engineering Science, 46, 2429-2444 12. Soltani M., Chen P., 2011, Numerical modeling of fluid flow in solid tumors, Plos One, 6, 6, e20344 13. Tan W.H.K., Wang F.J., Lee T., Wang C.H., 2003, Computer simulation of the delivery of etanidazole to brain tumor fromPLGAwafers:Comparisonbetween linear anddouble burst release systems,Biotechnology and Bioengineering, 82, 278-288 14. Teo C.S., Tan W.H.K., Lee T., Wang C.H., 2005, Transient interstitial fluid flow in brain tumors: Effects on drug delivery,Chemical Engineering Science, 60, 4803-4821 15. WangC.H., Li J., 1998,Three dimensional igg delivery to tumors,Chemical Engineering Science, 53, 3579-3600 16. Zhao J., SalmonH., Sarntinoranont M., 2007, Effect of heterogeneous vasculature on inter- stitial transport within a solid tumor,Microvascular Research, 73, 224-236 Manuscript received November 15, 2012; accepted for print October 28, 2013