Asymptotic analysis of a biphase tumor fluid flow. The weak coupling case
AASYMPTOTIC ANALYSIS OF A BIPHASE TUMOR FLUID FLOW. THE WEAKCOUPLING CASE.
CRISTINA VAGHI, SEBASTIEN BENZEKRY, CLAIR POIGNARDTEAM MONC, INRIA, INSTITUT DE MATH´EMATIQUES DE BORDEAUX, CNRS, BORDEAUX INP, UNIV.BORDEAUX, FRANCE
Abstract.
The aim of this paper is to investigate the asymptotic behavior of a biphase tumor fluid flowderived by 2-scale homogenisation techniques in recent works. This biphase fluid flow model accountsfor the capillary wall permeability, and the interstitial avascular phase, both being mixed in the limithomogenised problem. When the vessel walls become more permeable, we show that the biphase fluid flowexhibits a boundary layer that makes the computation of the full problem costly and unstable. In thelimit, both capillary and interstitial pressures coincide except in the vicinity of the boundary where differentboundary conditions are applied. Thanks to a rigorous asymptotic analysis, we prove that the solution tothe full problem can be approached at any order of approximation by a monophasic model with appropriateboundary conditions on the tumor boundary and appropriate correcting terms near the boundary are given.Numerical simulations in spherical geometry illustrate the theoretical results. Introduction
Motivation.
Drug delivery in tumors is affected by the fluid flow phenomena that occur within themalignant tissues, which include blood flow, interstitial convection and transvascular transport [1]. Anefficient quantification of these processes is of great importance to evaluate the drug penetration within thetumor site.Malignant tissues differ from normal tissues for several aspects [5]. Neoplastic vasculature is unevenlydistributed, leaving avascular spaces [15] and vessel walls are leaky and highly permeable [8]. Furthermore,the tumor interstitial matrix is dense and heterogeneous [12]. These features lead to an elevated interstitialfluid pressure (IFP) at the center of the tumor with a sharp drop at the periphery, impacting the transportof fluids and drug [2].Due to the high complexity of the architecture of tumors, two-scale asymptotic analysis has been employedto derive macroscopic models that take into account the microscopic properties of the malignant tissues[16, 13]. In a recent study, we derived several asymptotic models according to the magnitude of the hydraulicconductivity of the interstitium and of the vessel walls [18]. Under the assumption of a periodic structureof neoplastic tissues, tumors are modeled as double porous media and Darcy’s law describes the interstitialfluid flow and the blood transport. The coupling is driven by the permeability of the vessel walls. Abiphase model was derived under the assumption of low permeability of the vessel walls, while a monophasemodel was obtained under the assumption of leaky vessel walls. The former consists of a coupled systemof Darcy’s equation and it is in agreement with previous findings [16]. However, the numerical simulationsof the biphase model are computationally expensive due to the coupling of two elliptic equations. Indeed,the domain discretization has to be thick enough to catch the sharp pressure gradient that occurs at theperiphery of the tumor.Here, we prove that the biphase model can be approximated by the monophase model far from the bound-ary for large values of the hydraulic conductivity of the vessel walls. Moreover, boundary layer correctorscan describe precisely the behavior of the biphase model solutions in a simple and computationally efficientway.Our methodology can be applied to efficiently simulate the fluid transport in tumors, which might giveinsights in the drug delivery process in malignant tissues and can be applied to optimize treatments.1.2.
Model statement and objectives.
Given two Dirichlet boundary conditions π t and π c belongingto H s +1 / ( ∂ Ω) with s ≥ a r X i v : . [ m a t h . A P ] J a n hapman, Shipley et al [17] also derived recently in [18]. In this paper we focus on the case where theporosity tensors involved in [17, 18] are colinear. After simple change of notation, the problem reads in thesmooth domain Ω as − ∇ · ( E∇ p εt ) + α ε ( p εt − p εc ) = 0 , (1a) − ∇ · ( E∇ p εc ) − β ε ( p εt − p εc ) = 0 , (1b)with Dirichlet boundary conditions: p εt | ∂ Ω = π t , p εc | ∂ Ω = π c , on ∂ Ω.(1c)The tensor E is a positive definite tensor satisfying E min | X | ≤ X T E X ≤ E max | X | , ∀ X ∈ R , (2)and α/ε and β/ε are two positive constants that account for the permeability of the capillaries, the volumefraction of the capillary and the interstitium media, ε being a nondimension small parameter. Remark 1 (The weakly coupled case) . Problem (1) is the case where the two phases are weakly coupled.Indeed, performing the change of unknowns q ε = ( p εt − p εc ) / , p ε = ( p εt + p εc ) / , Problem (1) is equivalent to find ( p ε , q ε ) such that − ∇ · ( E∇ q ε ) + α + β ε q ε = 0 , (3a) − ∇ · ( E∇ p ε ) = − α − β ε q ε , (3b) q ε | ∂ Ω = π t − π c , p ε | ∂ Ω = π t + π c , on ∂ Ω , (3c) hence q ε is entirely determined by its boundary condition, while p ε involves q ε .To our opinion, the weakly couple case contains sufficient technical results, and describes already a lot ofapplications (for instance for scalar tensors [1] ) to justify the present paper. We are interested in the asymptotic regime ε tending to 0, which corresponds to leaky vessel walls. Inparticular, we show that p εt − p εc decays exponentially fast from the domain boundary, making appear atypical skin depth effect on the pressure difference. In addition we show that as ε goes 0, Problem (1) canbe approached by the solution ˜ P to the following monophase Laplace problem with well-designed boundarycondition: − ∇ · ( E∇ ˜ P ) = 0 , (4a) ˜ P | ∂ Ω = π t + π c − α − β α + β π t − π c , on ∂ Ω.(4b)The next section is devoted to the prove the well-posedness of the problem, and the preliminary estimatesof the solution ( q ε , p ε ). In particular, we prove the exponential decay from the boundary of q ε . In section 3,the asymptotic expansion of the solution at any order of ( q ε , p ε ) and optimal error estimates are given,proving the efficacy of the methodology. Numerical simulations illustrate the theoretical results in the lastSection . The highly coupled case, where the porosity tensors are not colinear for q ε qnd p ε will be treatedin a forthcoming work. 2. Well-posedness and a priori estimates
Proposition 1.
Let s ≥ and let π t and π c belongs to H / s ( ∂ Ω) , there exists a unique solution ( p εt , p εc ) to Problem (1) in ( H s (Ω)) . In addition, there exists a constant C independent of ε such that (cid:107) p εt − p εc (cid:107) L (Ω) ≤ C (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) (5a) (cid:107) p εt (cid:107) H (Ω) + (cid:107) p εc (cid:107) H (Ω) ≤ Cε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . (5b) roof. The uniqueness of the solution is obvious and left to the reader. Thanks to standard elliptic regular-ity [9, 10], one just has to prove the well-posedness in ( H (Ω)) . To prove the existence, denote by f t ( resp. f c ) a lift of π t ( π c resp. ) which belongs to H s (Ω) defined by − ∇ · ( E∇ f t ) = 0 , in Ω − ∇ · ( E∇ f c ) = 0 , in Ω f t | ∂ Ω = π t , f c | ∂ Ω = π c . It is well-known that there exists a constant C such that (cid:107) f t (cid:107) H (Ω) ≤ C (cid:107) π t (cid:107) H / ( ∂ Ω) , (cid:107) f c (cid:107) H (Ω) ≤ C (cid:107) π c (cid:107) H / ( ∂ Ω) , (6)Then, ( p εt , p εc ) reads ( p εt , p εc ) = ( f t + φ εt , f c + φ εc ) , where the couple ( φ εt , φ εc ) satisfy − ∇ · ( E∇ φ εt ) + α ε ( φ εt − φ εc ) = − α ε ( f t − f c ) − ∇ · ( E∇ φ εc ) − β ε ( φ εt − φ εc ) = β ε ( f t − f c ) , with homogeneous Dirichlet condition: φ εt | ∂ Ω = 0 = φ εc | ∂ Ω , on ∂ Ω,Consider the following bilinear form A ε defined on ( H (Ω)) by: ∀ ( u , v ) = (( u , u ) , ( v , v )) ∈ ( H (Ω)) × ( H (Ω)) , A ε ( u , v ) = β (cid:90) Ω E∇ u ∇ v dx + α (cid:90) Ω E∇ u ∇ v dx + α β ε (cid:90) Ω ( u − u )( v − v ) dx. Thanks to Poincar´e inequality, since E is coercive by (2), the bilinear form A is continuous and coerciveon ( H (Ω)) , and the coercivity constant does not depend on ε . Therefore there exists a unique solution φ ε satisfying for any v ∈ ( H (Ω)) : A ε ( φ ε , v ) = − α β ε (cid:90) Ω ( f t − f c )( v − v ) dx. (7)To prove the a priori estimates, thanks to Dirichlet and Neumann trace theorems, one just has to showthe inequalities on φ ε . Taking v = φ ε in (7), one infers β E min (cid:107)∇ φ εt (cid:107) L (Ω) + α E min (cid:107)∇ φ εc (cid:107) L (Ω) + α β ε (cid:107) φ εt − φ εc (cid:107) L (Ω) ≤ α β ε (cid:107) f t − f c (cid:107) L (Ω) (cid:107) φ εt − φ εc (cid:107) L (Ω) . Then one infers successively that for a constant C independent of ε (cid:107) φ εt − φ εc (cid:107) L (Ω) ≤ C (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) , and (cid:107)∇ φ εt (cid:107) L (Ω) + (cid:107)∇ φ εc (cid:107) L (Ω) ≤ Cε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . Estimates (5a)–(5b) fall then easily thanks to (6) since ( p εt , p εc ) = ( f t + φ εt , f c + φ εc ). (cid:3) Proposition 2.
For any d > , denote by Ω d the inner domain defined by Ω d = { x ∈ Ω : dist( x, ∂ Ω) > d } .There exists ε > , C d and µ > such that for any ε < ε , (cid:107) p εt − p εc (cid:107) H (Ω d ) ≤ C d e − µ/ε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . roof. The proof of the proposition is very similar to the proof of Haddar, Joly, Nguyen in [11] even thoughthe problem is slightly different. We recall here the main ideas for the self-consistency of the paper. Let φ be a smooth non negative function of Ω such that φ ( x ) = (cid:40) , if x ∈ Ω \ Ω d/ , µ, if x ∈ Ω d , where µ is a constant that is fixed later on. Denote by u ε the function of Ω defined by ( p εt − p εc )( x ) = e − φ ( x ) /ε u ε ( x ). It satisfies − ∇ · ( E∇ u ε ) + 1 ε (cid:0)(cid:0) E + E T (cid:1) ∇ φ · ∇ u ε (cid:1) + α + β − E∇ φ · ∇ φ + ε ∇ · ( E∇ φ ) ε u ε = 0 , in Ω ,u ε | ∂ Ω = π t − π c . Multiplying by u ε and integrating by parts lead to (cid:90) Ω E∇ u ε · ∇ u ε dx + 1 ε (cid:90) Ω (cid:0)(cid:0) E + E T (cid:1) ∇ φ · ∇ u ε (cid:1) u ε dx + 1 ε (cid:90) (cid:0) α + β − E∇ φ · ∇ φ + ε ∇ · ( E∇ φ ) (cid:1) u ε dx = (cid:90) ∂ Ω u ε ∂ n u ε dx = (cid:90) ∂ Ω ( p εt − p εc ) ∂ n ( p εt − p εc ) dx. Then using the fact that φ vanishes identically near ∂ Ω, one infers (cid:90) Ω (cid:0)(cid:0) E + E T (cid:1) ∇ φ · ∇ u ε (cid:1) u ε dx = − (cid:90) Ω ∇ · (cid:0)(cid:0) E + E T (cid:1) ∇ φ (cid:1) u ε dx, and Proposition 1 implies that (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ∂ Ω u ε ∂ n u ε dx (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ∂ Ω ( p εt − p εc ) ∂ n ( p εt − p εc ) dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ C ε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . Thus one infers the following estimate (cid:90) Ω E∇ u ε · ∇ u ε dx + 1 ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) (cid:16) α + β − E∇ φ · ∇ φ + ε ∇ · (( E − E T ) ∇ φ ) (cid:17) u ε dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ C ε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . One then just has to choose µ independent of ε such that (cid:107)E∇ φ · ∇ φ (cid:107) ≤ ( α + β ) / , to infer that for ε small enough (cid:13)(cid:13)(cid:13) ( α + β ) / − E∇ φ. ∇ φ + ε ∇ · (( E − E T ) ∇ φ (cid:13)(cid:13)(cid:13) L ∞ (Ω) ≥ α + β . Then (cid:107) p εt − p εc (cid:107) H (Ω d ) = e − µ/ε (cid:107) u ε (cid:107) H (Ω) ≤ Ce − µ/ε ε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) ≤ Ceµ e − µ/ε (cid:0) (cid:107) π t (cid:107) H / ( ∂ Ω) + (cid:107) π c (cid:107) H / ( ∂ Ω) (cid:1) . (cid:3) Asymptotic analysis
Consider here Problem (3) satisfied by q ε = ( p εt − p εc ) / p ε = ( p εt + p εc ) / − ∇ · ( E∇ q ε ) + α + β ε q ε = 0 , in Ω , − ∇ · ( E∇ p ε ) + α − β ε q ε = 0 , in Ω , ith Dirichlet boundary conditions: p ε | ∂ Ω = ( π t + π c ) / , q ε | ∂ Ω = ( π t − π c ) / , on ∂ Ω,3.1.
Geometry.
Let d > ∂ Ω, such that the tubular neighborhood of ∂ Ω of radiusd is parameterized by local coordinates. More precisely, let x T = ( x , x ) be a system of local coordinateson ∂ Ω = { Ψ( x T ) } . Define the map Φ by ∀ ( x T , x ) ∈ ∂ Ω × R , Φ( x T , x ) = Ψ( x T ) − x n ( x T ) , (8)where n is the normal vector of ∂ Ω outwardly directed. Then we assume that O d the tubular neighborhoodof ∂ Ω is parameterized as O d = { Φ( x T , x ) , ( x T , x ) ∈ ∂ Ω × (0 , d) } . The Euclidean metric in ( x T , x ) is given by the 3 × g ij ) i,j =1 , , where g ij = (cid:104) ∂ i Φ , ∂ j Φ (cid:105) : g = 1 , (9a) ∀ α ∈ { , } , g α = g α = 0 , (9b) ∀ ( α, β ) ∈ { , } , g αβ ( x T , x ) = g αβ ( x T ) − x b αβ ( x T ) + x c αβ ( x T ) , (9c)where g αβ = (cid:104) ∂ α Ψ , ∂ β Ψ (cid:105) , b αβ = (cid:104) ∂ α n , ∂ β Ψ (cid:105) , c αβ = (cid:104) ∂ α n , ∂ β n (cid:105) . (9d)3.2. The operator ∇ · ( E∇· ) in local coordinates. Define (cid:101) E the matrix E written in the new basis( ∂ i Φ) i =1 , , . In other words, (cid:101) E = P − E P, (10a)where P is the transfer matrix from the Euclidean basis to ( ∂ i Φ) i =1 , , : P = (cid:0) ∂ Φ , ∂ Φ , ∂ Φ (cid:1) . (10b)We denote by ( g ij ) the inverse matrix of ( g ij ) defined by (9), and by g the determinant of ( g ij ). For anyinteger l ≥ ε : a lij = ( − l ∂ l (cid:32) ∂ i (cid:0) √ gg ij (cid:1) √ g (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x =0 , for ( i, j ) ∈ { , , } , b lαβ = ( − l ∂ l (cid:0) g αβ (cid:1)(cid:12)(cid:12) x =0 , for ( α, β ) ∈ { , } , (11)and we denote by S l∂ Ω , (cid:101) E the sequence of differential operators on ∂ Ω of order 2 defined by S l∂ Ω , (cid:101) E = (cid:88) α,β =1 , (cid:101) E αβ (cid:0) a lαβ ∂ β + b lαβ ∂ α ∂ β (cid:1) . (12)Considering the rescaled local coordinates ( x T , ρ ) = ( x T , x /ε ), the operator ∇ · ( E∇· ) can be expandedin series of ε as follows: ∇ · ( E∇· ) = 1 √ g (cid:88) i,j =1 , , ∂ i (cid:16) √ g (cid:101) E ij g ij ∂ j (cid:17) , ∀ ( x T , x ) ∈ ∂ Ω × (0 , d) , = (cid:101) E ε ∂ ρ − (cid:101) E ε a ( x T ) ∂ ρ + (cid:88) (cid:96) ≥ ( − (cid:96) ε (cid:96) ρ (cid:96) (cid:96) ! L (cid:101) E ,(cid:96) , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , d /ε ) , (13)where the operators L (cid:101) E ,(cid:96) are of first order in ρ and of second order in x T , and given for any (cid:96) ≥ ∀ l ≥ , L (cid:101) E ,(cid:96) = (cid:18) ρ(cid:96) + 1 (cid:101) E a (cid:96) +133 ( x T ) ∂ ρ + S (cid:96)∂ Ω , (cid:101) E (cid:19) , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , d /ε ) . We refer to [14] for the detailed calculation of these expansions in the case of Laplace operator. .3. Asymptotic expansion of q ε . The problem satisfied by q ε is similar to the Helmholtz problem inhigh conductive materials that has been studied extensively in the last decade [11, 7, 4, 3, 6], in differentcontext. We recall here the main results regarding the expansion of q ε , and the proof is given in the Appendixsection refsec:appendix for self consistency of the paper.It is important to note that thanks to the properties of E given in (2), (cid:101) E is strictly positive.Denote by γ the positive parameter such that γ = α + β (cid:101) E , γ > . The following proposition provides the asymptotic expansion of q ε : Proposition 3.
Let d > , such that the mapping Φ defined by (8) is smooth and one-to-one. Define thesmooth cut-off function χ d equal to 1 in Ω \ Ω d , whose support is compactly embedded in Ω \ Ω d .Assume that Ω is a smooth domain, and that π t and π c are smooth functions of ∂ Ω .Then for any k ≥ , there exists C k depending on π t , π c and their derivatives such that (cid:107) q ε − k (cid:88) (cid:96) =0 χ d ε (cid:96) q (cid:96) ( x T , x /ε ) (cid:107) L (Ω) ≤ C k ε k +1 , (cid:107) q ε − k (cid:88) (cid:96) =0 χ d ε (cid:96) q (cid:96) ( x T , x /ε ) (cid:107) H (Ω) ≤ C k ε k +1 / , where the profiles q (cid:96) satisfy: • For k = 0 : − ∂ ρ q + γ q = 0 , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , + ∞ ) , q | η =0 = 12 ( π t − π c ) , q → ρ → + ∞ = 0 , (14) • For k ≥ : − ∂ ρ q k + γ q k = − a ( x T ) (cid:101) E ∂ ρ q k − + k − (cid:88) (cid:96) =0 ( − (cid:96) (cid:101) E ρ (cid:96) (cid:96) ! L (cid:101) E ,(cid:96) ( q k − − (cid:96) ) , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , + ∞ ) , q k | η =0 = 0 , q k → ρ → + ∞ = 0 , (15) where by convention, q n ≡ for n ≤ .More precisely, q ( x T , ρ ) = 12 ( π t − π c ) e − γρ , (16a) and for any k ≥ , q k ( x T , ρ ) = Q k ( x T , ρ ) e − γρ , (16b) where Q k is polynomial of degree k in the variable ρ , which vanishes in ρ = 0 and whose coefficients aresmooth functions of ∂ Ω .Proof. Problem (3a) reads in local coordinates in the vicinity of ∂ Ω as follows: ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , d /ε ) , − ∂ ρ q ε + γ q ε + ε a ( x T ) (cid:101) E ∂ ρ q ε − (cid:88) (cid:96) ≥ ( − (cid:96) (cid:101) E ε (cid:96) +2 ρ (cid:96) (cid:96) ! L (cid:101) E ,(cid:96) q ε = 0 , (17a) q ε | ρ =0 = 12 ( π t − π c ) . (17b) • Formal expansion.Set the Ansatz q ε ( x ) = χ d ( x ) (cid:88) k ≥ ε k q k ◦ Φ − ( x ) , ∀ x ∈ Ω , here q k are defined in ∂ Ω × (0 , + ∞ ).By Proposition 2, the exponential decay to 0 of q ε implies that the above ansatz is consistent in Ω d .Identifying the terms with the same power in ε k implies that the coefficients of the expansions q k satisfythe following inductive elementary problems. • For k = 0: − ∂ ρ q + γ q = 0 , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , + ∞ ) , q | η =0 = 12 ( π t − π c ) , q → ρ → + ∞ = 0 , • For k ≥ − ∂ ρ q k + γ q k = − a ( x T ) (cid:101) E ∂ ρ q k − + k − (cid:88) (cid:96) =0 ( − (cid:96) (cid:101) E ρ (cid:96) (cid:96) ! L (cid:101) E ,(cid:96) ( q k − − (cid:96) ) , ∀ ( x T , ρ ) ∈ ∂ Ω × (0 , + ∞ ) , q k | η =0 = 0 , q k → ρ → + ∞ = 0 , where by convention, q n ≡ n ≤
0. The expression of q is obvious since it satisfies the equation (14).Assuming that equality (16) is satisfied up to a give ( k − ≥
0. Then q k satisfies − ∂ ρ q k + γ q k = − a ( x T ) (cid:101) E ( ∂ ρ Q k − − γ Q k − ) e − γρ + e − γρ k − (cid:88) (cid:96) =0 ( − (cid:96) (cid:101) E ρ (cid:96) (cid:96) ! (cid:18) ρ(cid:96) + 1 a (cid:96) +133 ( ∂ ρ Q k − − (cid:96) − γ Q k − − (cid:96) ) + S (cid:96)∂ Ω , (cid:101) E ( Q k − − (cid:96) ) (cid:19) . (18)One then just has to observe that the righthand side of the above equality reads as R k − e − γρ , where R k − = (cid:80) k − (cid:96) =0 b (cid:96) ( x T ) ρ (cid:96) is polynomial in ρ and with smooth coefficients in ∂ Ω by induction hypothesis.Denoting by ( c (cid:96) ) k(cid:96) =1 the sequence of smooth functions defined on ∂ Ω by c k = b k − γk , for (cid:96) = k − , · · · , c (cid:96) = 12 γ(cid:96) (( (cid:96) + 1) (cid:96)c (cid:96) +1 + b (cid:96) ) , then q k = (cid:80) k(cid:96) =1 c (cid:96) ( x T ) ρ (cid:96) e − γρ satisfies (18) and vanishes at the boundaries ρ = 0 and ρ → + ∞ . • Proof of the estimates.The proof of the estimates is standard and we just recall here the main ingredients. The reader can referto [11] for more details. Thanks to Proposition 2, since (cid:107) q ε (cid:107) H (Ω d ) ≤ e − µ/ε for a specific µ >
0, and since (cid:80) k(cid:96) =0 χ d ε (cid:96) q (cid:96) ◦ Φ − also decays exponentially fast in Ω d one just has to prove the estimate in the tubularneighborhood O d of ∂ Ω, where local coordinates can be used. Using (13), simple a priori estimates showthat (cid:107) q ε − n (cid:88) (cid:96) =0 χ d ε (cid:96) q (cid:96) ( ., ./ε ) ◦ Φ − (cid:107) H (Ω) ≤ C n ε n − / , for any n . Then the proposition is proven by observing that since q k +1 P k +1 e − γρ , with P k +1 polynomial in ρ which vanishes in 0, one the following estimate (cid:107) q ε − k (cid:88) (cid:96) =0 χ d ε (cid:96) q (cid:96) ( ., ./ε ) ◦ Φ − (cid:107) H (Ω) ≤ (cid:107) q ε − k +2 (cid:88) (cid:96) =0 χ d ε (cid:96) q (cid:96) ( ., ./ε ) ◦ Φ − (cid:107) H (Ω) + ε k +1 (cid:107) χ d ( q k +1 + ε q k +2 ) (cid:107) H (Ω) ≤ C k ( ε k +3 / + ε k +1 / ) . (cid:3) .4. Asymptotic expansion of p ε = ( p εt + p εc ) / . We are now ready to approximate the function p ε =( p εt + p εc ) /
2, which is coupled to q ε by Problem (3b). Using the expansion of q ε given by Proposition (3), p ε satisfies formally − ∇ · ( E∇ p ε ) = − χ d α − β ε (cid:88) k ≥ ε k q k ( x T , x /ε ) , in Ω , (19a) p ε | ∂ Ω = ( π t + π c ) / , (19b)where by abuse of notation the mapping Φ is omitted, and where δ k is the Kronecker symbol equal to 1 is k = 0 and 0 elsewhere.Since the each term q k decays exponentially fast as x /ε , a fine mesh is necessary to solve the aboveequality in order to capture the the source term q ε . To prevent this drawback, we propose to determinean asymptotic expansion of p ε which splits between the fast variable x /ε and the slow variable . Moreprecisely, we look for p ε = (cid:88) k ≥ ε k ( P k + χ d ( x ) p k ( x T , x /ε )) , where each function p k describes the behavior near the boundary while P k is a function depending on theslow-variable far from the boundary. As before, the idea is to use the expansion of the operator ∇ · ( E∇· ) inthe local coordinates3.4.1.
Expansion of p ε . Injecting the Ansatz in Problem (3b) we infer that for any k ≥ p k ) k ≥ are defined on Γ × (0 , + ∞ ) vanish as ρ goes to infinity and satisfy ∂ ρ p = α − β (cid:101) E q , (20a)and for k ≥ ∂ ρ p k = α − β (cid:101) E q k + a ( x T ) (cid:101) E ∂ ρ p k − − k − (cid:88) n =0 ( − n (cid:101) E ρ n n ! L (cid:101) E ,n p k − n − . (20b) Proposition 4.
For any k ≥ , there exists a unique function p k defined on Γ × R + , vanishing as ρ tendsto infinity and satisfying Problem (20) . Moreover p k is given by p k = α − β α + β q k . More precisely, p ( x T , ρ ) = α − β α + β π t − π c e − γρ , and for any k ≥ , p k ( x T , ρ ) = α − β α + β Q k ( x T , ρ ) e − γρ , where Q k is polynomial of degree k in the variable ρ , which vanishes in ρ = 0 and whose coefficients aresmooth functions of ∂ Ω . In particular, p k | ρ =0 = δ k α − β α + β π t − π c . Proof.
Uniqueness easily comes from the fact that 2 functions satisfying (20) are equal modulo an affinefunction, which is necessary zero since it vanishes as ρ goes to infinity.To prove that p k equals ( α − β ) q k /γ , one just has to use induction. Indeed it is obviously true for k = 0, by the first equation of (20). Assume that this is true up to the rank k − ≥
0. One has, by definition f q k : α − β α + β ∂ ρ q k = α − β (cid:101) E q k + α − β α + β a ( x T ) (cid:101) E ∂ ρ q k − − α − β α + β k − (cid:88) n =0 ( − n (cid:101) E ρ n n ! L (cid:101) E ,n q k − n − = α − β (cid:101) E q k + a ( x T ) (cid:101) E ∂ ρ p k − − k − (cid:88) n =0 ( − n (cid:101) E ρ n n ! L (cid:101) E ,n p k − n − , by induction. Thus α − β α + β q k and p satisfy the same problem and both functions vanish at infinity, so theyare equal, which ends the proof. (cid:3) Then P is given by − ∇ · ( E∇ P ) = 0 , in Ω ,P | ∂ Ω = (cid:18) π t + π c (cid:19) − α − β α + β π t − π c , which is nothing but the problem (4) satisfied by ˜ P , and for any k ≥ P k = 0. One has the followingproposition. Proposition 5.
Assume that Ω is a smooth domain, and that π t and π c are smooth functions of ∂ Ω .Then for any k ≥ , there exists C k depending on π t , π c and their derivatives such that (cid:107) p ε − P − k (cid:88) n =0 ε n p n ( x T , x /ε ) (cid:107) L (Ω) ≤ C k ε k +1 , (cid:107) p ε − P − k (cid:88) n =0 ε n p n ( x T , x /ε ) (cid:107) H (Ω) ≤ C k ε k +1 / Proof.
The proof is standard and is performed recursively. Using Proposition 3, the expansion of ∇ · ( E∇· ),and by construction for any k ≥ (cid:107) p ε − P − k (cid:88) n =0 ε n p n ( x T , x /ε ) (cid:107) H (Ω) ≤ C k ε k − , then for a given k , one also has easily p ε − P − k (cid:88) n =0 ε n p n ( x T , x /ε ) = p ε − P − k +3 (cid:88) n =0 ε n p n ( x T , x /ε ) + ε k (cid:88) n =1 ε n p k + n ( x T , x /ε ) . By construction, for any n there exists a constant C n which depends on n , and π t , π c and their derivativessuch that (cid:107) χ d p n ( x T , x /ε ) (cid:107) L (Ω) ≤ C n √ ε, (cid:107) χ d p n ( x T , x /ε ) (cid:107) H (Ω) ≤ C n √ ε , which end the proof. (cid:3) Applications and numerical simulations
Expression of the first order of approximation.
Even though the above Propositions 3–5 makesit possible to derive the whole expansion of q ε and p ε , in the applications only the first order terms are oftensufficient. Using (14)–(15), one infers successively q ( x T , ρ ) = ( π t − π c )( x T )2 e − γρ , (21a) q ( x T , ρ ) = ( π t − π c )( x T )2 a ( x T )2 ρe − γρ , (21b) nd thus Proposition 4, one infers p ( x T , ρ ) = α − β γ ( π t − π c )( x T )2 e − γρ , (21c) p ( x T , ρ ) = α − β γ π t − π c a ( x T )2 ρe − γρ , (21d)while ˜ P satisfies (4).4.2. Numerical simulations.
To illustrate our results, we perfom simulation in the simple spherical case,with constant Dirichlet conditions on the boundary of the unit sphere. By symmetry, the problem to solvereads as follows d dr q ε + 2 r ddr q ε − α + β ε q ε = 0 , for 0 < r < , (22a) d dr p ε + 2 r ddr p ε = α − β ε q ε , for 0 < r < q ε | r =1 = π t − π c , p ε | r =1 = π t + π c , dq ε dr | r =0 = dp ε dr | r =0 = 0 . (22c)Denote by γ = (cid:112) α + β . Note that q ε admits the explicit formula q ε = π t − π c γr/ε ) r sinh( γε ) . On the contrary, the function p ε does not have explicit solution and it requires a numerical method to besolved. We used here the full second order finite difference method, by discretizing the flux at the order2 with a decentered stencil in r = 0, and by approaching at the order 2 the operator d dr + 1 r ddr at anypoint r i of the discretization of the segment (0 , r = 0 we multiplyequation (22b) by r and use the following approximation: r i ( d dr f + 1 r ddr f ) | r i ∼ r i +1 f i +1 − r i f i + r i − f i − δr + O ( δr ) , where δr is the path of the discretization. It is worth noting that due to the exponential decay in γ/ε of q ε , solving numerically the problem satisfied by p ε requires a very fine discretization especially the border r = 1, as shown in Figure 1.On the other hand, the solution to (4) and the first order profiles given by (21) provides the approximationat the order ε of both q ε and p ε by denoting by q approx, and p approx, respectively the following functions: q approx, ( r ) = ( π t − π c )( x T )2 e − γ (1 − r ) /ε + ( π t − π c )( x T )4 (1 − r ) e − γ (1 − r ) /ε , ∀ r ∈ (0 , p approx, ( r ) = π t + π c − α − β α + β (cid:18) π t − π c q approx, (cid:19) , ∀ r ∈ (0 , . (23b)Figure 2 illustrates the order of convergence of Propositions 3–5. References [1] L. T. Baxter and R. K. Jain. Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection.
Microvasc. Res. , 37(1):77–104, Jan. 1989.[2] Y. Boucher and R. K. Jain. Microvascular pressure is the principal driving force for interstitial hypertension in solid tumors:Implications for vascular collapse.
Cancer Res. , 52(18):5110–5114, Sept. 1992.[3] G. Caloz, M. Dauge, E. Faou, and V. P´eron. On the influence of the geometry on skin effect in electromagnetism.
ComputerMethods in Applied Mechanics and Engineering , 200(9-12):1053–1068, 2011.[4] G. Caloz, M. Dauge, and V. P´eron. Uniform estimates for transmission problems with high contrast in heat conductionand electromagnetism.
Journal of Mathematical Analysis and Applications , 370(2):555–572, May 2010.[5] V. P. Chauhan, T. Stylianopoulos, Y. Boucher, and R. K. Jain. Delivery of Molecular and Nanoscale Medicine to Tumors:Transport Barriers and Strategies.
Annu. Rev. Chem. Biomol. Eng. , 2(1):281–298, July 2011. .0 0.2 0.4 0.6 0.8 1.0 r q function q , for = 0.1 function q , for = 0.06 function q , for = 0.04 r p function p , for = 0.1 function p , for = 0.06 function p , for = 0.04 Figure 1.
Plots of the pressures q ε (Left) and p ε (Right) for 3 values of ε : 0.1, 0.07 and0.04. As ε goes to zero, a boundary layer appears near the boundary, making the numericalcomputation of the full problem costly and unstable. The following constants have beenchosen: π t = 0 . π c = 1, α = 1, β = 1 .
5, and the discretization step is δr = 10 − . e rr o r qq a pp r o x , theoretical slope L -norm of the error (slope:2.5)theoretical slope H -norm of the error (slope:1.5) 10 e rr o r pp a pp r o x , theoretical slope L -norm of the error (slope:2.4)theoretical slope H -norm of the error (slope:1.5) Figure 2. (Left): Convergence rate of the L -norm and H -norm of q ε − q approx, . (Right):Convergence rate of the L -norm and H -norm of p ε − p approx, . The following constantshave been chosen: π t = 0 . π c = 1, α = 1, β = 1 .
5. In this simple particular case, thenumerical slopes of the errors are slightly better than the theoretical results for the L norm. [6] M. Dauge, P. Dular, L. Kr¨ ahenb¨uhl, V. P´eron, R. Perrussel, and C. Poignard. Corner asymptotics of the magneticpotential in the eddy-current model. Mathematical Methods in the Applied Sciences , 37(13):1924–1955, 2014.
7] M. Dauge, E. Faou, and V. P´eron. Comportement asymptotique `a haute conductivit´e de l’´epaisseur de peau en´electromagn´etisme.
Comptes Rendus Math´ematique , 348(7-8):385–390, Mar. 2010.[8] H. F. Dvorak, L. F. Brown, M. Detmar, and A. M. Dvorak. Vascular permeability factor/vascular endothelial growthfactor, microvascular hyperpermeability, and angiogenesis.
Am. J. Pathol. , 146(5):1029–1039, May 1995.[9] L. Evans and A. M. Society.
Partial Differential Equations . Graduate studies in mathematics. American MathematicalSociety, 1998.[10] D. Gilbarg and N. Trudinger.
Elliptic Partial Differential Equations of Second Order . Grundlehren der mathematischenWissenschaften. Springer Berlin Heidelberg, 2013.[11] H. Haddar, P. Joly, and H.-M. Nguyen. Generalized impedance boundary conditions for scattering by strongly absorbingobstacles: The scalar case.
Mathematical Models and Methods in Applied Sciences , 15(08):1273–1300, 2005.[12] R. K. Jain. Transport of molecules in the tumor interstitium: A review.
Cancer Res. , 47(12):3039–3051, June 1987.[13] R. Penta, D. Ambrosi, and A. Quarteroni. Multiscale homogenization for fluid and drug transport in vascularized malignanttissues.
Math. Models Methods Appl. Sci. , 25(01):79–108, Jan. 2015.[14] R. Perrussel and C. Poignard. Asymptotic expansion of steady-state potential in a high contrast medium with a thinresistive layer.
Applied Mathematics and Computation , 221(0):48 – 65, 2013.[15] E. M. Sevick and R. K. Jain. Viscous resistance to blood flow in solid tumors: Effect of hematocrit on intratumor bloodviscosity.
Cancer Res. , 49(13):3513–3519, July 1989.[16] R. J. Shipley and S. J. Chapman. Multiscale Modelling of Fluid and Drug Transport in Vascular Tumours.
Bull Math Biol ,72(6):1464–1491, Aug. 2010.[17] R. J. Shipley and S. J. Chapman. Multiscale modelling of fluid and drug transport in vascular tumours.
Bulletin ofMathematical Biology