Nonlinear Analysis: Real World Applications 35 (2017) 483–502
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa
Bifurcation for a free-boundary tumor model with angiogenesis Yaodan Huang a , Zhengce Zhang a,∗ , Bei Hu b,c a
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, 710049, China Institute for Mathematical Sciences, Renmin University, Beijing 100872, China c Department of Applied Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA b
article
info
Article history: Received 9 August 2016 Received in revised form 3 December 2016 Accepted 7 December 2016
Keywords: Bifurcation Free boundary problem Stationary solution Symmetry-breaking Angiogenesis
abstract We consider a free boundary tumor model with vasculature which supply nutrients to ∂σ the tumor, so that ∂n +β(σ−σ) = 0 holds on the boundary, where a positive constant β is the rate of nutrient supply to the tumor and σ is the nutrient concentration outside the tumor. The tumor cells proliferate at a rate µ. We show that for each µ2 , µ4 , µ6 , µ8 , · · ·, symmetry-breaking stationary solutions bifurcate from the radially symmetric stationary solutions. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Over the last few decades, an increasing number of mathematical models describing solid tumor growth have been studied and developed; these models are classified into discrete cell-based models and continuum models. At the tissue level, continuum models provide a very good approximation. These models incorporate a system of partial differential equations, where variables such as cell density (or volume fractions), nutrient (i.e., oxygen and glucose), etc. are tracked. Modeling, mathematical analysis, numerical simulations were carried out in numerous papers, such as [1–45] and the references cited there. Lowengrub et al. [46] provided a systematic survey of tumor model studies. Angiogenesis is an important physiological process through which new blood vessels form from pre-existing vessels. It is an essential process in wound healing. However, in a live tissue, angiogenesis is a process that tumor cells secret cytokines that stimulate the vascular system to grow toward the tumor. That is to say, ∗ Corresponding author. E-mail addresses:
[email protected] (Y. Huang),
[email protected] (Z. Zhang),
[email protected] (B. Hu).
http://dx.doi.org/10.1016/j.nonrwa.2016.12.003 1468-1218/© 2016 Elsevier Ltd. All rights reserved.
484
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
the tumor possesses its own vasculature, then nutrient that is essential to the survival of the tumor may be supplied via the capillary network. The tumor core is nonnecrotic and no chemical inhibitor is present. Analysis of this model can lead to a better understanding of the growing mechanism of the tumor and may also assist in the treatment of cancer. We denote the concentration of nutrients by σ, distributed uniformly throughout the tumor, then σ satisfies a reaction–diffusion equation ∂σ = D∆σ + Γ (σB − σ) − λσ ∂t
in Ω (t),
(1.1)
where Ω (t) is the tumor domain at time t with a moving boundary ∂Ω (t), D represents diffusion coefficient, Γ (σB −σ) denotes the nutrient supplied by the vasculature with Γ being the transfer rate of nutrient in blood tissue and σB being the concentration of nutrients in the vasculature, and λσ describes the consumption of the nutrient by cells accompanied with λ being the rate of consumption. This model was proposed in [6] to describe the evolution of tumor growth. By appropriate change of variables and non-dimensional scales (cf. [28], see also [6,10]), (1.1) can be rewritten as c
∂σ = ∆σ − σ ∂t
in Ω (t),
(1.2)
where c is a small parameter. The tumor tissue with angiogenesis receives nutrients through its own vasculature, so it is natural to assume that σ satisfies the boundary condition (see [28]): ∂σ + β(t)(σ − σ) = 0 ∂n
on ∂Ω (t),
(1.3)
where n is the outward normal, β(t) is a positive-valued function which may vary in time and σ is the nutrient concentration outside the tumor. Angiogenesis results in an increase in β(t); on the other hand, if the tumor is treated with anti-angiogenic drugs, β(t) will decrease and the starved tumor will shrink. The pressure p stems from the transport of cells which proliferate or die. As in [10,31], we shall make the assumption that the tumor tissue behaves like a porous medium. Let V be the velocity field of the tumor cell movement, then V = −∇p by Darcy’s law. Combining with the conservation of mass, divV = µ(σ − σ ), we get − ∆p = µ(σ − σ )
in Ω (t),
(1.4)
where µ is a parameter expressing the “intensity” of the expansion by mitosis and σ is a threshold concentration. As in [4,8], the cell-to-cell adhesiveness leads to the boundary condition: p = κ on ∂Ω (t),
(1.5)
where κ is the mean curvature. Furthermore, assuming the continuity of the velocity field up to the boundary of the domain, we derive Vn = V · n = −∇p · n = −
∂p ∂n
on ∂Ω (t),
(1.6)
where Vn is the velocity of the free boundary in the direction n. In this paper, we assume that β(t) ≡ β is a positive constant.
(1.7)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
If the solution is radially symmetric, then Ω (t) = {r < R(t)}, by (1.4) and (1.6), 1 dR(t) ∂p = =− −∆pdV dt ∂r r=R(t) 4πR2 (t) Ω(t) R(t) 1 = µ(σ − σ ) · 4πr2 dr 4πR2 (t) 0 R(t) µ (σ − σ )r2 dr. = 2 R (t) 0
485
(1.8)
In this case, it is proved in [28] that, under the assumptions (1.7) and σ < σ,
(1.9)
there exists a unique radially symmetric stationary solution to (1.2)–(1.5) and (1.8). In fact, if the boundary condition (1.3) is replaced by σ = σ (which is formally the case β(t) = ∞), it was proved in [22] that a branch of symmetry-breaking stationary solutions bifurcates from the radially symmetric stationary solutions for each µn (RS ) (n ≥ 2) with free boundary r = RS + εYn,0 (θ) + O(ε2 )
(1.10)
where RS is the radius of the radially symmetric stationary solution and Yn,0 is the spherical harmonic of order (n, 0). For the two dimensions, this was done in [30]. In the sequel, Friedman and Hu [23] determined the threshold value µ∗ (RS ) for which the radial stationary solution changes from linear stability to linear instability under small radial perturbations, and also proved the spherical stationary solution is nonlinearly stable for µ(RS ) < µ∗ (RS ) in [24]. It was proved in [23] that there is a critical number R (R ≈ 0.62207 if σ = 1) such that µ∗ (RS ) < µ2 (RS ) if RS < R and µ∗ (RS ) = µ2 (RS ) if RS > R. Since µ2 (RS ) is the smallest among all µn (RS ) (n ≥ 3), the bifurcation branches stemming from µn (RS ) (n ≥ 3) are all unstable. Biologically, these unstable bifurcation branches are unlikely be observed in “stationary” mode. Thus the bifurcation branch near µ = µ2 (RS ) is much more significant than any other bifurcation branches. Indeed, it was shown in [27] that it is possible to have both stable and unstable bifurcation branches stemming from µ2 (RS ). For the general case that the nutrient consumption rate and the cell proliferation rate are not linear, Cui and Escher have proved that there exists a branch of non-radial solutions bifurcating from the radially symmetric solution in the neighborhood of each γn (surface tension coefficient) with n > n∗∗ in [16] and found a critical value γ ∗ which separates instability from stability for the radially symmetric equilibrium with respect to non-symmetric perturbations in [15]. In the presence of inhibitor, the asymptotic stability of radially symmetric stationary solution was studied in [11,12,39]. In a more recent article [37], Wang established for each µn (n > n∗∗ ), the existence of a family of non-radially symmetric stationary solutions, bifurcating from a unique radially symmetric stationary solution. Bifurcation from radially symmetric stationary solution into non-radially stationary solution results also extend to the model of tumor growth with a necrotic core [34], and tumor growth in fluid-like tissue which is modeled by Stokes equation [25]. Furthermore, for the tumor growth modeled by Stokes equation, a critical value N∗ (RS , γ) is found in [26] such that the spherical stationary solution is linearly stable if µ/γ < N∗ (RS , γ) and linearly unstable if µ/γ > N∗ (RS , γ) where N∗ (RS , γ) ≤ minn≥2 Mn (RS ). Numerical simulations for bifurcation branches outside a small neighborhood of the bifurcation point were carried out in [33,35,36]. In addition, other models have already been studied and developed in many papers. Escher and Matioc [18–20] studied the growth of nonnecrotic tumors in three different regimes of vascularization and established that for the rate of mitosis Gk (k > k1 ), a series of non-radially symmetric equilibria bifurcate from a unique radially symmetric stationary solution. For a multi-layer tumor model in the absence and presence of inhibitor, the existence of flat solution and non-flat solution and the asymptotic stability of
486
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
flat equilibrium were analyzed in [14,42–44]. If the boundary are nonlinear with Gibbs–Thomson relation, Zhou and Wu [45] showed that γn (n > n∗∗ ) is a bifurcation point of the flat stationary solution. Recently, Wu [38] determined a threshold value θ∗ for which there are two radial stationary solutions for 0 < σ /σ < θ∗ . He proved there exists a family of bifurcation branches of solutions for each radial stationary solution corresponding to γn (n > n∗∗ ) as well. In this paper, we shall work with the boundary condition (1.3) (with the assumption (1.7)) stemming from angiogenesis. We shall establish Theorem 1.1. µn (RS ) given by µn (RS ) =
n 3n(n − 1)(n + 2) RS + β + RS Pn (RS ) · , 2 σ RS4 (n + βRS )P1 (RS ) − (1 + βRS )Pn (RS )
where Pn (RS ) is defined by (2.12), is a monotonically increasing sequence: µ2 (RS ) < µ3 (RS ) < µ4 (RS ) < · · · < µn (RS ) < µn+1 (RS ) < · · · . Furthermore, for every even n ≥ 2, µ = µn (RS ) is a bifurcation point of the symmetry-breaking stationary solution of the system (1.2)–(1.6) with free boundary r = RS + εYn,0 (θ, ϕ) + o(ε).
(1.11)
Our proof of the monotonicity of µn (RS ) is quite technical and somewhat lengthy. However, we expect, just like the case β(t) = ∞, the bifurcation stemming from µ2 (RS ) is the most significant one biologically. Establishment of monotonicity of µn (RS ) enables us to include the bifurcation result for µ = µ2 (RS ). The structure of this paper is as follows. In Section 2, we give some preliminaries which will be needed in the sequel and then we carry out our proof of Theorem 1.1 in Section 3. 2. Preliminaries In this section, for convenience, we collect some properties of the spherical harmonic function Yn,m (θ, ϕ) and the modified Bessel function In (ξ) which will be used in the sequel. We then solve the radially symmetric stationary solution of the problem (1.2)–(1.6) explicitly. Moreover, we state the Crandall–Rabinowitz theorem, which is crucial in our proof. The family of functions {Yn,m } forms a complete orthonormal basis for L2 (Σ ) (cf. [47]), ∆ω Yn,m = −n(n + 1)Yn,m , (2.1) ∂ ∂ ∂2 sin θ ∂θ + sin12 θ ∂ϕ where Σ is the unit sphere and ∆ω = sin1 θ ∂θ 2 is the Laplace operator on Σ . The Laplace operator in R3 can then be rewritten as ∆=
∂2 2 ∂ 1 + + 2 ∆ω . 2 ∂r r ∂r r
(2.2)
The Bessel function In (ξ) given by In (ξ) =
n 2k ∞ ξ 1 ξ 2 k!Γ (n + k + 1) 2
(2.3)
k=0
satisfies (cf. [22,23]) I1/2 (ξ) =
2 sinh ξ, πξ
(2.4)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
I3/2 (ξ) =
2 πξ
−
sinh ξ + cosh ξ , ξ
487
(2.5)
n In (ξ) = In−1 (ξ), n ≥ 1, ξ n In′ (ξ) − In (ξ) = In+1 (ξ), n ≥ 0, ξ 2n In−1 (ξ) − In+1 (ξ) = In (ξ), n ≥ 1, ξ In′ (ξ) +
(2.6) (2.7) (2.8)
In−1 (ξ)In+1 (ξ) < In2 (ξ),
ξ > 0, 2 In−1 (ξ)In+1 (ξ) > In2 (ξ) − In (ξ)In+1 (ξ), ξ > 0, ξ ∞ Γ (m + n + 2k + 1)(ξ/2)m+n+2k . Im (ξ)In (ξ) = k!Γ (m + k + 1)Γ (n + k + 1)Γ (m + n + k + 1)
(2.9) (2.10) (2.11)
k=0
In addition, the function Pn (ξ) defined by Pn (ξ) =
In+3/2 (ξ) , ξIn+1/2 (ξ)
n = 0, 1, 2, 3, . . .
(2.12)
satisfies (by [23, Lemmas 2.1 and 2.3]), for ξ real and positive and for any n ≥ 0, 1 , n+1 (ξ) + (2n + 3) 1 Pn (0) = , 2n + 3 d 1 2n + 3 Pn (ξ) = − Pn (ξ) − ξPn2 (ξ), dξ ξ ξ Pn (ξ) =
ξ2P
(2.13) (2.14) (2.15)
Pn (ξ) > Pn+1 (ξ),
(2.16)
d Pn (ξ) < 0. dξ
(2.17)
A radially symmetric stationary solution of the system (1.2)–(1.6) satisfies (r = |x|) 2 σ ′′ (r) + σ ′ (r) = σ(r), 0 < r < RS , r ∂σ + β(σ − σ) = 0, ∂r r=RS 2 p′′ (r) + p′ (r) = −µ(σ(r) − σ ), 0 < r < RS , r
(2.18) (2.19) (2.20)
and p(RS ) =
1 , RS
p′ (RS ) = 0.
(2.21)
By rescaling if necessary, we assume without loss of generality that σ = 1,
(and therefore σ < σ = 1).
(2.22)
A unique solution of (2.18)–(2.19) is given explicitly in [28] by 1
RS2 I1/2 (r) β σS (r) = , β + g(RS ) r 12 I1/2 (RS )
0 < r < RS ,
(2.23)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
488
where g(RS ) = coth RS −
I3/2 (RS ) 1 = = RS P0 (RS ). RS I1/2 (RS )
(2.24)
Let q = p + µσS , then we can solve q in the form q=
1 µ σ r 2 + C1 , 6
where C1 is a constant. Combining with the boundary condition of p in (2.21), we derive 1 pS (r) = −µσS (r) + µ σ r 2 + C1 , 6 with C1 =
1 RS
+
µβ β+g(RS )
0 < r < RS ,
(2.25)
− 16 µ σ RS2 , and RS satisfies β g(RS ) 1 . = σ β + g(RS ) RS 3
(2.26)
By [28, Theorem 3.1], for 0 < β < ∞ and 0 < σ < 1, there is a unique solution RS for Eq. (2.26), producing a unique solution for (2.18)–(2.21). Next, we state a very useful theorem, the Crandall–Rabinowitz theorem, which is critical in studying bifurcation. Theorem 2.1 (See [48, Theorem 1.7]). Let X, Y be real Banach spaces and F (x, µ) a C p map, p ≥ 3, of a neighborhood (0, µ0 ) in X × R into Y . Suppose (i) (ii) (iii) (iv)
F (0, µ) = 0 for all µ in a neighborhood of µ0 ; Ker[Fx (0, µ0 )] is a one-dimensional space, spanned by x0 ; Im[Fx (0, µ0 )] = Y1 has codimension 1; [Fµx ](0, µ0 )x0 ̸∈ Y1 .
Then (0, µ0 ) is a bifurcation point of the equation F (x, µ) = 0 in the following sense: in a neighborhood of (0, µ0 ) the set of solutions of F (x, µ) = 0 consists of two C p−2 smooth curves Γ1 and Γ2 which intersect only at the point (0, µ0 ); Γ1 is the curve (0, µ) and Γ2 can be parameterized as follows: Γ2 : (x(ε), µ(ε)), |ε| small, (x(0), µ(0)) = (0, µ0 ), x′ (0) = x0 .
3. Proof of Theorem 1.1 In order to tackle the existence of symmetry-breaking stationary solutions of the system (1.2)–(1.6), we would like to apply the Crandall–Rabinowitz theorem. For simplicity, we denote RS by R. Consider a family ϕ) where R(θ, ϕ) = εS(θ, ϕ). Let (σ, p) be the of domains with perturbed boundaries ∂Ωε : r = R + R(θ, solution of ∆σ = σ
in Ωε ,
∂σ + β(σ − 1) = 0 on ∂Ωε , ∂n ∆p = −µ(σ − σ ) in Ωε , p = κ on ∂Ωε .
(3.1) (3.2) (3.3) (3.4)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
489
We define F by µ) = ∂p , F (R, ∂n ∂Ωε
(3.5)
µ) = 0. so that (σ, p) is a symmetry-breaking stationary solution if and only if F (R, We introduce the Banach space: ∈ C l+α (Σ ), R is π-periodic in θ, 2π-periodic in ϕ}, X l+α = {R X2l+α = {closure of the linear space spanned by {Yj,0 , j = 0, 2, 4, . . .} in X l+α }. It will be shown that F (·, µ) maps from a neighborhood of 0 of X l+3+α to X l+α and that the map clearly maps X l+3+α into X l+α (see (3.28) and (3.29)). By [22, p190], Condition A is satisfied for {Yj,0 ; j = 2, 4, 6, . . .}. µ) also maps X l+3+α into X l+α . We choose M = X 3+α (Σ ), N = X α (Σ ). Therefore, F (R, 2 2 2 2 In order to apply the Crandall–Rabinowitz theorem, we need to compute Fr´echet derivatives of F . To do that, we follow the procedure outlined in [27] to compute the expansion of (σ, p) of order ε as follows: σ = σS + εσ1 + O(ε2 ), 2
p = pS + εp1 + O(ε ).
(3.6) (3.7)
Now we proceed to formally compute σ1 and p1 explicitly. We use ⃗er , ⃗eθ , ⃗eϕ to denote the unit normal vectors in the r, θ, ϕ directions, respectively, ⃗er = (sin θ cos ϕ, sin θ sin ϕ, cos θ), ⃗eθ = (cos θ cos ϕ, cos θ sin ϕ, − sin θ), ⃗eϕ = (− sin ϕ, cos ϕ, 0). In addition, the unit outer normal vector n is given in [27]: 1 εSθ εSϕ n= ⃗er − ⃗eθ − ⃗eϕ R + εS (R + εS) sin θ 1 + |ε∇ω S|2 /(R + εS)2 = ⃗er −
εSϕ εSθ ⃗eθ − ⃗eϕ + O(ε2 ), R R sin θ
(3.8)
and 1 1 ∂ϕ . ∇ = ⃗er ∂r + ⃗eθ ∂θ + ⃗eϕ r r sin θ Clearly, ∂σ1 = ∇σ1 · n ∂n r=R+εS r=R+εS εSϕ ∂∇σ1 εSθ ⃗eθ − ⃗eϕ + O(ε2 ) = ∇σ1 + εS · ⃗er − ∂r r=R R R sin θ r=R ∂σ1 = + O(ε). ∂r r=R Furthermore, from (3.2) and (3.6), we have ∂σ1 ∂σ1 + βσ1 =ε + βσ1 + O(ε2 ) ε ∂r ∂n ∂BR ∂Ωε ∂(σ − σS ) = + O(ε2 ) + β(σ − σS ) ∂n ∂Ωε ∂σ ∂σS (R + εS) = + βσ − + βσS (R + εS) + O(ε2 ) ∂n ∂n ∂Ωε
(3.9)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
490
∂σS (R) ∂σS (R) ∂ 2 σS (R) =β− εS + β + βσS (R) + εS ∂r ∂r2 ∂r 2 ∂σS (R) ∂ σS (R) +β = −ε S + O(ε2 ) (by (2.19)) ∂r2 ∂r
+ O(ε2 )
, −ελS + O(ε2 ),
(3.10)
where λ=
∂σS ∂ 2 σS +β ∂r2 ∂r
.
(3.11)
r=R
By (2.6) and (2.7), a direct computation from (2.23) shows that I3/2 (R) β ∂σS = , ∂r r=R β + g(R) I1/2 (R) and β ∂ 2 σS 2 I3/2 (R) = 1− . ∂r2 r=R β + g(R) R I1/2 (R) Therefore, λ=
β β + g(R)
I3/2 (R) β 2 I3/2 (R) +β 1 − 2P0 (R) + βRP0 (R) . 1− = R I1/2 (R) I1/2 (R) β + RP0 (R)
(3.12)
We now expand the boundary condition for p. The mean curvature expansion is computed in [27,49]: ε 1 1 κ= − 2 S(θ, ϕ) + ∆ω S(θ, ϕ) + O(ε2 ). (3.13) R R 2 Combining with (3.7) and (2.21), we then obtain εp1 (R) = εp1 (R + εS) + O(ε2 ) = p(R + εS) − pS (R + εS) + O(ε2 ) = κ − pS (R) − p′S (R)εS + O(ε2 ) ε 1 = − 2 S(θ, ϕ) + ∆ω S(θ, ϕ) + O(ε2 ). R 2
(3.14)
We thus derive the following equations of order ε: ∆σ1 = σ1
in BR ,
∂σ1 + βσ1 = −λS(θ, ϕ) on ∂BR , ∂r ∆p1 = −µσ1 in BR , 1 1 p1 = − 2 S(θ, ϕ) + ∆ω S(θ, ϕ) R 2
(3.15) (3.16) (3.17) on ∂BR .
(3.18)
We set S(θ, ϕ) = Yn,m (θ, ϕ) in (3.15)–(3.18). Using a separation of variables, we seek a solution of the form σ1 (r, θ, ϕ) = σ 1 (r)Yn,m (θ, ϕ),
(3.19)
p1 (r, θ, ϕ) = p1 (r)Yn,m (θ, ϕ).
(3.20)
We first solve σ1 . Substituting (3.19) into (3.15)–(3.16) and combining with (2.1), we derive 2 ′ n(n + 1) σ 1′′ (r) + σ (r) − σ 1 (r) = σ 1 (r), r 1 r2 ∂ σ1 + β σ1 = −λ. ∂r r=R
0 < r < R,
(3.21) (3.22)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
491
Eq. (3.21) admits a solution of the form (r) 1 In+1/2 . σ 1 (r) = C 1 r2 Substituting this into (3.22) and using also (2.7), we obtain 1
1 = C
−λR 2 . nR−1 In+1/2 (R) + In+3/2 (R) + βIn+1/2 (R)
It follows that 1
In+1/2 (r) −λR 2 . σ 1 (r) = 1 nR−1 In+1/2 (R) + In+3/2 (R) + βIn+1/2 (R) r2
(3.23)
Hence, the solution of problem (3.15)–(3.16) is given by 1
σ1 (r, θ, ϕ) =
In+1/2 (r) −λR 2 Yn,m (θ, ϕ). 1 −1 nR In+1/2 (R) + In+3/2 (R) + βIn+1/2 (R) r2
(3.24)
Next, we solve p1 . Let η = p1 + µσ1 , then η satisfies ∆η = 0. We look for a solution of the form η(r, θ, ϕ) = η(r)Yn,m (θ, ϕ). By (2.1), η(r) satisfies n(n + 1) 2 η = 0, η′′ + η′ − r r2 and we solve this problem in the form 2 rn . η(r) = C 2 rn Yn,m (θ, ϕ). Combining this with (2.1), (3.18) and (3.20), we Namely, p1 (r, θ, ϕ) = −µ σ1 (r)Yn,m (θ, ϕ) + C derive 1 n(n + 1) µλ 1 C2 = n − 2 1 − − , R R 2 hn (R) + β where hn (R) =
In+3/2 (R) n n + = + RPn (R). R In+1/2 (R) R
(3.25)
So we can solve p1 in the form 3 p1 (r, θ, ϕ) = −µσ1 (r, θ, ϕ) + C 3 = − 12 1 − with C R
n(n+1) 2
−
r n R
Yn,m (θ, ϕ)
(3.26)
µλ hn (R)+β .
We next want to show that the formal expansions carried out above for (3.6) and (3.7) are actually rigorous. In fact, σ is defined only on Ωε , σS is defined everywhere on R3 and σ1 is defined only on BR , so we shall first transform these three functions into the same domain Ωε by the Hanzawa transformation which is a diffeomorphism defined by r = r′ + χ(R − r′ )εS(θ, ϕ),
θ = θ′ , ϕ = ϕ′ ,
492
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
namely (r, θ, ϕ) = Hε (r′ , θ′ , ϕ′ ), with ∞
χ∈C ,
3 0, if |z| ≥ δ0 4 , χ(z) = 1, if |z| < 1 δ0 4
k d χ C3 dz k ≤ δ k , 0
where δ0 is positive and small. We point out that Hε maps BR onto Ωε while keeping the ball {r < R − 43 δ0 } fixed and the inverse transformation Hε−1 maps Ωε onto BR . Setting σ 1 (r, θ, ϕ) = σ1 (Hε−1 (r, θ, ϕ))
in Ωε ,
(3.27)
then σ, σS and σ 1 are all defined on the same domain Ωε . Now, we rigorously establish the expansions in (3.6) and (3.7) by estimating the O(ε2 ) terms in the C 3+α -norm. Lemma 3.1. ∥σ − σS − ε σ1 ∥C 3+α (Ω ε ) ≤ C|ε|2 ∥S∥C 3+α (Σ ) ,
(3.28)
∥p − pS − ε p1 ∥C 1+α (Ω ε ) ≤ C|ε|2 ∥S∥C 3+α (Σ ) ,
(3.29)
where C is independent of ε and S for 0 < |ε| < 1, ∥S∥C 3+α (Σ ) ≤ 1. Proof. Recall that 1 1 ∂ϕ ∇ = ⃗er ∂r + ⃗eθ ∂θ + ⃗eϕ r r sin θ and n=
1
1 + |ε∇ω S|2 /(R + εS)2 = ⃗er + O(|εS| + |ε∇ω S|).
⃗er −
εSϕ εSθ ⃗eθ − ⃗eϕ R + εS (R + εS) sin θ
For σ 1 , we have −∆ σ +σ 1 = f in Ωε , 1 ∂ σ1 + β σ1 = −λS + k. ∂n ∂Ωε Since f came from various term of the Hanzawa transform after applying the Laplace operator, f involved second order derivatives of εS. Therefore, ∥f∥C 1+α (Ω ε ) ≤ C4 |ε| ∥S∥C 3+α (Σ ) . Similarly, since k involved first order derivatives of εS, ∥ k∥C 2+α (Σ ) ≤ C5 |ε| ∥S∥C 3+α (Σ ) . Let ψ = σ − σS − ε σ1 , then it satisfies −∆ψ + ψ = −εf in Ωε , ∂ψ ∂σS + βψ =β− + βσS + ελS − ε k. ∂n ∂n ∂Ωε ∂Ωε
(3.30)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
493
We now proceed to compute, ∂σS + ελS k1 , β − + βσS ∂n ∂Ωε 1 εSθ εSϕ = β − ∇σS (R + εS) · ⃗er − ⃗eθ − ⃗eϕ R + εS (R + εS) sin θ 1 + |ε∇ω S|2 /(R + εS)2 − βσS (R + εS) + ελS 1 ∂r σS (R + εS) − βσS (R + εS) + ελS =β− 1 + |ε∇ω S|2 /(R + εS)2 = β − ∂r σS (R + εS) − βσS (R + εS) + ελS + O(|ε∇ω S|2 ) = β − ∂r σS (R) − ∂rr σS (R)εS − βσS (R) − β∂r σS (R)εS + ελS + O(|εS|2 + |ε∇ω S|2 ) = O(|εS|2 + |ε∇ω S|2 ), where we made use of (2.19) and (3.11). In a similar manner, we derive ∥ k1 ∥C 2+α (Σ ) ≤ C6 |ε|2 ∥S∥C 3+α (Σ ) . Combining these estimates, we find that ∥ k1 − ε k∥C 2+α (Σ ) ≤ C7 |ε|2 ∥S∥C 3+α (Σ ) . By the Schauder estimates, we derive (3.28). And in the similar manner, we can obtain (3.29). Hence, the lemma is complete. The estimates we obtained are uniformly valid for |ε| small, ∥εS∥C 3+α (Σ ) small, ∥S∥C 3+α (Σ ) bounded. Thus the procedures carried out above give the Fr´echet derivatives of F , which we now proceed to compute explicitly. Since ∂pS (R) = 0, ∂r we have 2 ∂ pS (R) ∂p1 (R, θ, ϕ) ∂(pS + ε p1 ) =ε S+ + O(|ε|2 ∥S∥C 3+α (Σ ) ), (3.31) ∂n ∂r2 ∂r r=R+εS(θ,ϕ) and so, by (3.29), ∂p =ε ∂n r=R+εS(θ,ϕ)
∂ 2 pS (R) ∂p1 (R, θ, ϕ) S+ ∂r2 ∂r
+ O(|ε|2 ∥S∥C 3+α (Σ ) ).
(3.32)
µ) → F (R, µ) from X l+3+α to X l+α is bounded if l = 0, and The relation (3.32) shows that the mapping (R, the same argument shows that the same is true for any l ≥ 0. A similar argument shows that this mapping µ); furthermore ∂F (R, µ)/∂ R (or, ∂F (R, µ)/∂µ) is obtained by solving a is Fr´echet differentiable in (R, linearized problem about (R, µ) with respect to R (or, µ). By using the Schauder estimates, we can then µ) to any order. further obtain differentiability of F (R, µ) in R, at (0, µ), is given by From (3.32) it follows that the Fr´echet derivative of F (R, ∂ 2 pS (R) ∂p1 (R, θ, ϕ) FR (0, µ) S(θ, ϕ) = S(θ, ϕ) + . 2 ∂r ∂r As a matter of fact, from (2.20), (2.21) and (2.23), we obtain ∂ 2 pS (R) β = −µ(σ (R) − σ ) = −µ − σ , S ∂r2 β + g(R)
(3.33)
(3.34)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
494
and from (3.26), for S(θ, ϕ) = Yn,m (θ, ϕ), n(n + 1) ∂p1 (R, θ, ϕ) ∂σ1 (R, θ, ϕ) 1 µλ n = −µ + − 2 1− − Yn,m (θ, ϕ). ∂r ∂r R 2 hn (R) + β R
(3.35)
By (3.16) and (3.24), we have ∂σ1 (R, θ, ϕ) = −βσ1 (R, θ, ϕ) − λYn,m (θ, ϕ) ∂r βλIn+1/2 (R) Yn,m (θ, ϕ) − λYn,m (θ, ϕ) = nR−1 In+1/2 (R) + In+3/2 (R) + βIn+1/2 (R) βλ = Yn,m (θ, ϕ) − λYn,m (θ, ϕ) In+3/2 (R) n + + β R In+1/2 (R) =
−λhn (R) Yn,m (θ, ϕ), hn (R) + β
(3.36)
where hn (R) is from (3.25). Hence, after substituting into (3.35), ∂p1 (R, θ, ϕ) µλhn (R) 1 µλ n n(n + 1) = Yn,m (θ, ϕ) + − 2 1 − − Yn,m (θ, ϕ) ∂r hn (R) + β R 2 hn (R) + β R n n(n + 1) µλ µλ n = −1 + hn (R) − Yn,m (θ, ϕ) R3 2 hn (R) + β hn (R) + β R n n(n + 1) µλ n = − 1 + h (R) − Yn,m (θ, ϕ) n R3 2 hn (R) + β R In+3/2 (R) µλ n n(n + 1) − 1 + Yn,m (θ, ϕ). = R3 2 hn (R) + β In+1/2 (R) Substituting (3.34) and (3.37) into (3.33), we derive FR (0, µ) Yn,m = (An (R) − µBn (R))Yn,m ,
(3.37)
(3.38)
where n An (R) = 3 R Bn (R) =
n(n + 1) −1 2
=
n(n − 1)(n + 2) , 2R3
In+3/2 (R) β λ −σ − . β + g(R) hn (R) + β In+1/2 (R)
(3.39) (3.40)
The equation
FR (0, µ) Yn,m (θ, ϕ) = 0
(3.41)
is then reduced to An (R) − µBn (R) = 0. We now establish: Lemma 3.2. Bn (R) > 0 for any n ≥ 2 and B1 (R) = A1 (R) = 0, A0 (R) = 0, B0 (R) ̸= 0, R ∈ (0, ∞). Proof. By (2.26) and (2.24), we get β σ R σ = = . β + g(R) 3 g(R) 3P0 (R)
(3.42)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
495
From (3.40): In+3/2 (R) β λ −σ − β + g(R) hn (R) + β In+1/2 (R) 1 − 2P0 (R) + βRP0 (R) RPn (R) β β = −σ − (by (3.12), (2.12)) β + g(R) β + g(R) hn (R) + β 1 − 2P0 (R) + βRP0 (R) RPn (R) 1 σ n −3− = (by (3.42), (3.25)) 3 P0 (R) P0 (R) R + RPn (R) + β σ (A + 1 + βR)RPn (R) = A− n , 3 R + β + RPn (R)
Bn =
where A =
1 P0 (R)
(3.43)
− 3. By (2.13), the expression A is simplified to A = R2 P1 (R).
(3.44)
Substituting (3.44) into (3.43), we derive 2 R P1 (R) + 1 + βR Pn (R) σ R Bn = RP1 (R) − n 3 R + β + RPn (R) =
σ R (n + βR)P1 (R) − (1 + βR)Pn (R) · > 0, n 3 R + β + RPn (R)
since, by (2.16), P1 (R) > Pn (R) for n > 1.
(3.45)
We denote the solution of An (R) − µBn (R) = 0 by µn (R). Using (3.39) and (3.45), we obtain µn (R) =
n 3n(n − 1)(n + 2) An (R) R + β + RPn (R) = · . Bn (R) 2 σ R4 (n + βR)P1 (R) − (1 + βR)Pn (R)
(3.46)
We proceed to establish our main lemma. Lemma 3.3. µn (R) is monotonically increasing: µ2 (R) < µ3 (R) < µ4 (R) < µ5 (R) < · · · .
(3.47)
Proof. By (3.46) and Lemma 3.2, we only need to establish the inequality n 3n(n − 1)(n + 2) R + β + RPn (R) · 2 σ R4 (n + βR)P1 (R) − (1 + βR)Pn (R)
<
n+1 3(n + 1)n(n + 3) R + β + RPn+1 (R) · , 4 2 σR (n + 1 + βR)P1 (R) − (1 + βR)Pn+1 (R)
(3.48)
which is equivalent to n
+ β + RPn (R) (n + 1 + βR)P1 (R) − (1 + βR)Pn+1 (R) R n+1 < (n + 1)(n + 3) + β + RPn+1 (R) (n + βR)P1 (R) − (1 + βR)Pn (R) . R
(n − 1)(n + 2)
The inequality (3.49) can be rewritten as (n + 1)(n + 3)[RP1 (R) − RPn (R)] − (n − 1)(n + 2)[RP1 (R) − RPn+1 (R)] β 2 n+1 + (n + 1)(n + 3) + RPn+1 (R) [RP1 (R) − RPn (R)] + (n + 1)(n + 3)[nP1 (R) − Pn (R)] R
(3.49)
496
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
n − (n − 1)(n + 2) + RPn (R) [RP1 (R) − RPn+1 (R)] − (n − 1)(n + 2) (n + 1)P1 (R) − Pn+1 (R) β R n+1 + (n + 1)(n + 3) + RPn+1 (R) [nP1 (R) − Pn (R)] R n − (n − 1)(n + 2) + RPn (R) (n + 1)P1 (R) − Pn+1 (R) > 0, R which simplifies to (3n + 5)RP1 (R) + (n − 1)(n + 2)RPn+1 (R) − (n + 1)(n + 3)RPn (R) β 2 + (3n + 5)(2n + 1)P1 (R) + (n − 1)(n + 1)(n + 2)Pn+1 (R) − (n + 1)(n + 2)(n + 3)Pn (R) + (n + 1)(n + 3)R2 P1 (R)Pn+1 (R) − (n − 1)(n + 2)R2 P1 (R)Pn (R) − (3n + 5)R2 Pn (R)Pn+1 (R) β n(n − 1)(n + 2) (n + 1)2 (n + 3) n(n + 1)(3n + 5) P1 (R) + Pn+1 (R) − Pn (R) + R R R + n(n + 1)(n + 3)RP1 (R)Pn+1 (R) − (n − 1)(n + 1)(n + 2)RP1 (R)Pn (R) − (3n + 5)RPn (R)Pn+1 (R) , Q1 (n, R)β 2 + Q2 (n, R)β + Q3 (n, R) > 0.
(3.50)
We now proceed to show that the quantities within each pair of braces in (3.50), which are denoted by Q1 (n, R), Q2 (n, R), Q3 (n, R) respectively, are positive, and therefore (3.50) must be positive, which will complete our proof. We shall divide our proofs into several lemmas. Lemma 3.4. Let C1 (n) > 0, C2 (n) > 0, and fn (x) = C1 (n)Pn+1 (x) − C2 (n)Pn (x).
(3.51)
If fn (0) > 0, then fn (x) > 0 for x ∈ (0, ∞). Proof. Using (2.15), we derive 1 2n + 5 1 2n + 3 ′ 2 2 fn (x) = C1 (n) − Pn+1 (x) − xPn+1 (x) − C2 (n) − Pn (x) − xPn (x) x x x x 2n + 3 =− fn (x) − x C1 (n)Pn+1 (x) − C2 (n)Pn (x) Pn+1 (x) x C (n) − C (n) 2C (n) 1 2 1 + C2 (n)Pn+1 (x) − C2 (n)Pn (x) Pn (x) + − Pn+1 (x) x x 2n + 3 =− + xPn+1 (x) fn (x) + xC2 (n)Pn (x) Pn (x) − Pn+1 (x) x C1 (n) − C2 (n) − 2C1 (n)Pn+1 (x) + x , E(n, x)fn (x) + J(n, x), (3.52) 2n+3 where E(n, x) = − x + xPn+1 (x) . By (2.16), Pn (x) − Pn+1 (x) > 0. Since we have assumed fn (0) > 0, we deduce, by (2.14) and (2.17), 2C1 (n) C1 (n) − C2 (n) − 2C1 (n)Pn+1 (x) ≥ C1 (n) − C2 (n) − 2n + 5 C1 (n) C2 (n) = (2n + 3) − 2n + 5 2n + 3
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
497
= (2n + 3)fn (0) > 0. It follows that J(n, x) > 0.
(3.53)
Hence fn (x) = e
x 0
E(n,s)ds
fn (0) + e
x 0
E(n,s)ds
x
e
−
r 0
E(n,s)ds
J(n, r)dr > 0,
for x > 0,
(3.54)
0
which completes the proof.
Lemma 3.5. The quantity Q1 (n, R) in (3.50) is positive. Proof. We need to show (3n + 5)RP1 (R) + (n − 1)(n + 2)RPn+1 (R) − (n + 1)(n + 3)RPn (R) > 0
(3.55)
or, by the definition of Pn (R) from (2.12), 2 (3n + 5)I5/2 In+1/2 In+3/2 + (n − 1)(n + 2)I3/2 In+1/2 In+5/2 > (n + 1)(n + 3)I3/2 In+3/2 .
(3.56)
Since, by (2.10), 2 In+3/2 In+5/2 , R the second term in (3.56) will be replaced by this inequality. It suffices to show that 2 In+1/2 In+5/2 > In+3/2 −
(3n + 5)I5/2 In+1/2 > (3n + 5)I3/2 In+3/2 + (n − 1)(n + 2)
2 I3/2 In+5/2 . R
(3.57)
By (2.11), (3.57) is equivalent to ∞
n+3+2k (3n + 5)Γ (n + 2k + 4) R k!Γ (k + 7/2)Γ (n + k + 3/2)Γ (n + k + 4) 2 k=0 n+3+2k ∞ 1 (3n + 5)Γ (n + 2k + 4) (n − 1)(n + 2)Γ (n + 2k + 5) R + . > k!Γ (k + 5/2) Γ (n + k + 5/2)Γ (n + k + 4) Γ (n + k + 7/2)Γ (n + k + 5) 2 k=0
Since all terms from both sides are positive, it suffices to show that the coefficients corresponding to R 2k+n+3 will satisfy the inequality, i.e., 2 (3n + 5)Γ (n + 2k + 4) k!Γ (k + 7/2)Γ (n + k + 3/2)Γ (n + k + 4) 1 (3n + 5)Γ (n + 2k + 4) (n − 1)(n + 2)Γ (n + 2k + 5) > + . k!Γ (k + 5/2) Γ (n + k + 5/2)Γ (n + k + 4) Γ (n + k + 7/2)Γ (n + k + 5) This inequality is equivalent to 5 3 (3n + 5) n + k + n+k+ (n + k + 4) > (n − 1)(n + 2) k + 52 (n + 2k + 4) 2 2 5 5 n+k+ (n + k + 4). + (3n + 5) k + 2 2 After simplification, this inequality is reduced to 1 (6n3 + 44n2 + 95n + 10n2 k + 2nk 2 + 37nk + 2k 2 + 29k + 60) > 0, 2 which is obviously true for any n ≥ 1, k ≥ 0. Therefore, this lemma is completed.
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
498
Lemma 3.6. The quantity Q2 (n, R) in (3.50) is positive. Proof. By (2.16), P1 (R) > Pn (R), we obtain Q2 (n, R) = (3n + 5)(2n + 1)P1 (R) + (n − 1)(n + 1)(n + 2)Pn+1 (R) − (n + 1)(n + 2)(n + 3)Pn (R) + (n + 1)(n + 3)R2 P1 (R)Pn+1 (R) − (n − 1)(n + 2)R2 P1 (R)Pn (R) − (3n + 5)R2 Pn (R)Pn+1 (R) > (n − 1)(n + 1)(n + 2)Pn+1 (R) + (3n + 5)(2n + 1) − (n + 1)(n + 2)(n + 3) Pn (R) + (n + 1)(n + 3)R2 P1 (R)Pn+1 (R) − (n − 1)(n + 2)R2 P1 (R)Pn (R) − (3n + 5)R2 Pn (R)Pn+1 (R) = (n − 1) (n2 + 3n + 2)Pn+1 (R) − (n2 + n − 1)Pn (R) (n + 1)(n + 3) 1 (n − 1)(n + 2) 1 1 2 + (3n + 5)R P1 (R)Pn (R)Pn+1 (R) − − 3n + 5 Pn (R) 3n + 5 Pn+1 (R) P1 (R) , (n − 1)J1,n (R) + (3n + 5)R2 P1 (R)Pn (R)Pn+1 (R)J2,n (R). Since n2 + 3n + 2 n2 + n − 1 − 2n + 5 2n + 3 2 2n + 10n + 11 = > 0, (2n + 5)(2n + 3)
(n2 + 3n + 2)Pn+1 (0) − (n2 + n − 1)Pn (0) =
by Lemma 3.4, J1,n (R) is positive. So we only need to show that (n + 1)(n + 3) 1 (n − 1)(n + 2) 1 1 − − 3n + 5 Pn (R) 3n + 5 Pn+1 (R) P1 (R) 1 1 1 = (bn + 1) − bn − > 0, Pn (R) Pn+1 (R) P1 (R)
J2,n (R) =
for any n > 1 and R ∈ (0, ∞), where bn =
(n−1)(n+2) . 3n+5
(3.58)
Using (2.14), we derive
(n + 1)(n + 3)(2n + 3) (n − 1)(n + 2)(2n + 5) − −5 3n + 5 3n + 5 4n2 + 2n − 6 = >0 3n + 5
J2,n (0) =
for any n > 1. Let ψn (R) =
1 Pn (R) ,
ψn′ (R) = −
then by (2.15), 1 Pn2 (R)
Pn′ (R) = −
2n + 3 1 2 ψn (R) + ψn (R) + R. R R
It follows that ′ J2,n (R)
(3.59)
1 2 2n + 3 = (bn + 1) − ψn (R) + ψn (R) + R R R 2n + 5 1 2 5 1 2 ψn+1 (R) + R − − ψ1 (R) + ψ1 (R) + R − bn − ψn+1 (R) + R R R R 2n + 3 1 2 =− (bn + 1)ψn2 (R) − bn ψn+1 (R) − ψ12 (R) + (bn + 1)ψn (R) R R 2n + 5 5 − bn ψn+1 (R) − ψ1 (R) R R 2 1 2 2 =− (bn + 1)ψn (R) − bn ψn+1 (R) − (bn + 1)ψn (R) − bn ψn+1 (R) − J2,n (R) R 5 2n − 2 2n + J2,n (R) + (bn + 1)ψn (R) − bn ψn+1 (R) R R R
(3.60)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
1 2 5 = J2,n (R) − [(bn + 1)ψn (R) − bn ψn+1 (R)] + J2,n (R) R R R 2 bn (bn + 1) 2 2 + ψn (R) − ψn+1 (R) + (bn + 1)(n − 1)ψn (R) − bn nψn+1 (R) R R R , W1 (n, R)J2,n (R) + W2 (n, R).
499
(3.61)
Let Hn (R) = Then, recalling bn =
(n−1)(n+2) 3n+5
2 2 (bn + 1)(n − 1)ψn (R) − bn nψn+1 (R). R R
and ψn (R) =
1 Pn (R) ,
2 n−1 1 (n + 1)(n + 3)Pn+1 (R) − n(n + 2)Pn (R) R 3n + 5 Pn (R)Pn+1 (R) 1 2 n−1 n (R). , H R 3n + 5 Pn (R)Pn+1 (R)
Hn (R) =
(3.62)
Since 2 n (0) = (n + 1)(n + 3) − n(n + 2) = 2n + 8n + 9 > 0, H 2n + 5 2n + 3 (2n + 5)(2n + 3)
(3.63)
by Lemma 3.4, we have n (R) > 0. H
(3.64)
(3.62) implies that W2 (n, R) > 0. Thus J2,n (R) = e
R 0
W1 (n,s)ds
J2,n (0) + e
R 0
W1 (n,s)ds
R
e
−
r 0
W1 (n,s)ds
W2 (n, r)dr > 0,
(3.65)
0
which completes the proof.
Lemma 3.7. The quantity Q3 (n, R) in (3.50) is positive. Proof. A direct computation shows, Q3 (n, R) =
n(n + 1)(3n + 5) n(n − 1)(n + 2) (n + 1)2 (n + 3) P1 (R) + Pn+1 (R) − Pn (R) R R R + n(n + 1)(n + 3)RP1 (R)Pn+1 (R) − (n − 1)(n + 1)(n + 2)RP1 (R)Pn (R)
− (3n + 5)RPn (R)Pn+1 (R) 1 = n(n + 1)(3n + 5)P1 (R) + n(n − 1)(n + 2)Pn+1 (R) − (n + 1)2 (n + 3)Pn (R) R + R2 n(n + 1)(n + 3)P1 (R)Pn+1 (R) − (n − 1)(n + 1)(n + 2)P1 (R)Pn (R) − (3n + 5)Pn (R)Pn+1 (R) . Since n(n + 1)(3n + 5)P1 (R) > n(n + 1)(3n + 5)Pn (R) > (n + 1)2 (n + 3)Pn (R),
(n ≥ 1),
we can easily deduce that the first bracket is positive. For the second bracket, we utilize the inequality P1 (R) > Pn (R) to replace Pn (R) by P1 (R) for the last term, n(n + 1)(n + 3)P1 (R)Pn+1 (R) − (n − 1)(n + 1)(n + 2)P1 (R)Pn (R) − (3n + 5)Pn (R)Pn+1 (R)
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
500
> P1 (R) n(n + 1)(n + 3)Pn+1 (R) − (n − 1)(n + 1)(n + 2)Pn (R) − (3n + 5)Pn+1 (R) = (n − 1)P1 (R) (n2 + 5n + 5)Pn+1 (R) − (n2 + 3n + 2)Pn (R) . It is not difficult to verify n2 + 5n + 5 n2 + 3n + 2 − 2n + 5 2n + 3 2n2 + 6n + 5 = > 0, (2n + 3)(2n + 5)
(n2 + 5n + 5)Pn+1 (0) − (n2 + 3n + 2)Pn (0) =
so that by Lemma 3.4, the second bracket is also positive. Combining these, we find that Q3 (n, R) > 0.
Combining Lemmas 3.5–3.7, we have completed the proof of Lemma 3.3.
Now, we are ready to prove our main result Theorem 1.1. Proof of Theorem 1.1. We have computed [FR (0, µ)]Yn,0 = (An − µBn )Yn,0 ,
n = 0, 1, 2, 3, . . . ,
therefore the kernel in the space M Ker [FR (0, µ)] = 0,
if An − µBn ̸= 0 for any n = 0, 2, 4, 6, . . . ,
and Ker [FR (0, µ)] =
span{Yj,0 }
if µ = µn .
{j:µj =µn }
Since we have proved Lemma 3.3 that µj (R) are all distinct, we actually have Ker [FR (0, µ)] = span{Yn,0 }
if µ = µn .
Hence, dim(Ker [FR (0, µn )]) = 1. On the other hand, [FR (0, µn )]Yk,0 = (Ak − µn Bk )Yk,0 , where Ak − µn Bk ̸= 0, for k ̸= n. This implies that, for even n ≥ 2, Y1 = Im [FR (0, µn )] = span{Y0,0 , Y2,0 , Y4,0 , . . . , Yn−2,0 , Yn+2,0 , . . .}. Since Y1
span{Yn,0 } = N , we conclude codim Y1 = 1.
Finally, it is clear that [FµR (0, µn )]Yn,0 = −Bn Yn,0 ̸∈ Im [FR (0, µn )]. By the Crandall–Rabinowitz theorem, for even integer n ≥ 2, µn are bifurcation points of the stationary solution of the system (1.2)–(1.6). Since the perturbation is non-radially symmetric, this produces a family of non-radially symmetric stationary solutions.
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
501
4. Conclusion and biological interpretation In this paper, we have studied bifurcation solutions of a free boundary problem modeling the growth of solid tumors with angiogenesis. By taking the cell proliferation rate µ, characterizing the ‘aggressiveness’ of the tumor, as a bifurcation parameter, we show that µ = µ2 (RS ) is the first bifurcation point of the symmetry-breaking stationary solution. In fact, µ = µ2 (RS ) is the most significant one biologically (see [22–24,27]). Note that bifurcation solutions with free boundary r = RS + εYn,0 (θ, ϕ) + o(ε) may be regarded as protrusions, which are associated with the invasion of tumors into their surrounding stroma. It would be interesting to find out the angiogenesis effects on tumor growth. The parameter β represents the rate of angiogenesis. By (1.3), if β increases, the tumor will receive sufficient nutrients from vasculature ( σ < σ), our bifurcation result implies that the tumor may probably loose its spherical shape, develop protrusions and become invasive. In particular, when β = ∞, which indicates that the tumor is completely surrounded by blood vessels and thus regarded as cultured tumor in vitro, is the Dirichlet boundary condition studied intensively in many papers which confirm that this is the case. The implication of our result is that controlling angiogenesis has a positive impact in limiting the ability of tumors invasion, which depends on the growth of protrusions generated by bifurcation solutions. In medical science, it is common to use anti-angiogenic drugs to reduce the invasion of tumors into normal tissue in therapies. Acknowledgments We would like to thank the referees for their valuable recommendations. This work was supported by the National Natural Science Foundation of China (Nos. 11371286, 11401458). Huang is very grateful to Renmin University for providing help during the completion of this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
J.A. Adam, A simplified mathematical model of tumor growth, Math. Biosci. 81 (1986) 229–244. J.A. Adam, A mathematical model of tumor growth: III. Comparison with experiment, Math. Biosci. 86 (1987) 213–227. J.A. Adam, N. Bellomo, A Survey of Models for Tumor-Immune System Dynamics, Birkh¨ auser, Boston, 1997. H.M. Byrne, The importance of intercelluar adhesion in the development of carcinomas, IMA J. Math. Appl. Med. Biol. 14 (1997) 305–323. N.F. Britton, M.A.J. Chaplain, A qualitative analysis of some models of tissue growth, Math. Biosci. 113 (1993) 77–89. H.M. Byrne, M.A.J. Chaplain, Growth of nonnecrotic tumors in the presence and absence of inhibitors, Math. Biosci. 130 (1995) 151–181. H.M. Byrne, M.A.J. Chaplain, Growth of necrotic tumors in the presence and absence of inhibitors, Math. Biosci. 135 (1996) 187–216. H.M. Byrne, M.A.J. Chaplain, Modelling the role of cell–cell adhesion in the growth and development of carcinomas, Math. Comput. Modelling 24 (1996) 1–17. H.M. Byrne, M.A.J. Chaplain, Free boundary value problems associated with the growth and development of multicellular spheroids, European J. Appl. Math. 8 (1997) 639–658. V. Cristini, J. Lowengrub, Q. Nie, Nonlinear simulation of tumor growth, J. Math. Biol. 46 (2003) 191–224. S. Cui, Analysis of a mathematical model for the growth of tumors under the action of external inhibitors, J. Math. Biol. 44 (2002) 395–426. S. Cui, A. Friedman, Analysis of a mathematical model of the effect of inhibitors on the growth of tumors, Math. Biosci. 164 (2000) 103–137. S. Cui, A. Friedman, Analysis of a mathematical model of the growth of necrotic tumors, J. Math. Anal. Appl. 255 (2001) 636–677. S. Cui, J. Escher, Well-posedness and stability of a multi-dimensional tumor growth model, Arch. Ration. Mech. Anal. 191 (2009) 173–193. S. Cui, J. Escher, Asymptotic behaviour of solutions of a multidimensional moving boundary problem modeling tumor growth, Comm. Partial Differential Equations 33 (2008) 636–655. S. Cui, J. Escher, Bifurcation analysis of an elliptic free boundary problem modelling the growth of avascular tumors, SIAM J. Math. Anal. 39 (2007) 210–235.
502
Y. Huang et al. / Nonlinear Analysis: Real World Applications 35 (2017) 483–502
[17] J. Escher, A.-V. Matioc, Analysis of a two-phase model describing the growth of solid tumors, European J. Appl. Math. 24 (2013) 25–48. [18] J. Escher, A.-V. Matioc, Bifurcation analysis for a free boundary problem modeling tumor growth, Arch. Math. 97 (2011) 79–90. [19] J. Escher, A.-V. Matioc, Well-posedness and stability analysis for a moving boundary problem modelling the growth of nonnecrotic tumors, Discrete Contin. Dyn. Syst. Ser. B 15 (2011) 573–596. [20] J. Escher, A.-V. Matioc, Radially symmetric growth of nonnecrotic tumors, Nonlinear Differential Equations Appl. 17 (2010) 1–20. [21] A. Friedman, B. Bazaliy, Global existence and asymptotic stability for an elliptic–parabolic free boundary problem: an application to a model of tumor growth, Indiana Univ. Math. J. 52 (2003) 1265–1304. [22] M.A. Fontelos, A. Friedman, Symmetry-breaking bifurcations of free boundary problems in three dimensions, Asymptot. Anal. 35 (2003) 187–206. [23] A. Friedman, B. Hu, Bifurcation from stability to instability for a free boundary problem arising in a tumor model, Arch. Ration. Mech. Anal. 180 (2006) 293–330. [24] A. Friedman, B. Hu, Asymptotic stability for a free boundary problem arising in a tumor model, J. Differential Equations 227 (2006) 598–639. [25] A. Friedman, B. Hu, Bifurcation for a free boundary problem modeling tumor growth by Stokes equation, SIAM J. Math. Anal. 39 (2007) 174–194. [26] A. Friedman, B. Hu, Bifurcation from stability to instability for a free boundary problem modeling tumor growth by Stokes equation, J. Math. Anal. Appl. 327 (2007) 643–664. [27] A. Friedman, B. Hu, Stability and instability of Liapunov-Schmidt and Holf bifucation for a free boundary problem arising in a tumor model, Trans. Amer. Math. Soc. 360 (2008) 5291–5342. [28] A. Friedman, K.-Y. Lam, Analysis of a free-boundary tumor model with angiogenesis, J. Differential Equations 259 (2015) 7636–7661. [29] A. Friedman, F. Reitich, Analysis of a mathematical model for the growth of tumors, J. Math. Biol. 38 (1999) 262–284. [30] A. Friedman, F. Reitich, Symmetry-breaking bifurcation of analytic solutions to free boundary problems: an application to a model of tumor growth, Trans. Amer. Math. Soc. 353 (2001) 1587–1634. [31] H.P. Greenspan, Models for the growth of a solid tumor by diffusion, Stud. Appl. Math. 51 (1972) 317–340. [32] H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theoret. Biol. 56 (1976) 229–242. [33] W. Hao, J.D. Hauenstein, B. Hu, A.J. Sommese, A three-dimensional steady-state tumor system, Appl. Math. Comput. 218 (2011) 2661–2669. [34] W. Hao, J.D. Hauenstein, B. Hu, Y. Liu, A.J. Sommese, Y. Zhang, Bifurcation for a free boundary problem modeling the growth of a tumor with a necrotic core, Nonlinear Anal. RWA 13 (2012) 694–709. [35] W. Hao, J.D. Hauenstein, B. Hu, Y. Liu, A.J. Sommese, Y. Zhang, Continuation along bifurcation branches for a tumor model with a necrotic core, J. Sci. Comput. 53 (2012) 395–413. [36] W. Hao, J.D. Hauenstein, B. Hu, T. McCoy, A.J. Sommese, Computing steady-state solutions for a free boundary problem modeling tumor growth by Stokes equation, J. Comput. Appl. Math. 237 (2013) 326–334. [37] Z.J. Wang, Bifurcation for a free boundary problem modeling tumor growth with inhibitors, Nonlinear Anal. RWA 19 (2014) 45–53. [38] J. Wu, Stationary solutions of a free boundary problem modeling the growth of tumors with Gibbs–Thomson relation, J. Differential Equations 260 (2016) 5875–5893. [39] J. Wu, S. Cui, Asymptotic behaviour of solutions of a free boundary problem modelling the growth of tumours in the presence of inhibitors, Nonlinearity 20 (2007) 2389–2408. [40] J.P. Ward, J.R. King, Mathematical modelling of avascular-tumour growth, IMA J. Math. Appl. Med. Biol. 14 (1997) 39–69. [41] J.P. Ward, J.R. King, Mathematical modelling of avascular-tumour growth II: Modelling growth saturation, IMA J. Math. Appl. Med. Biol. 16 (1999) 171–211. [42] F. Zhou, S. Cui, Bifurcation for a free boundary problem modeling the growth of multi-layer tumors, Nonlinear Anal. 68 (2008) 2128–2145. [43] F. Zhou, J. Escher, S. Cui, Well-posedness and stability of a free boundary problem modeling the growth of multi-layer tumors, J. Differential Equations 244 (2008) 2909–2933. [44] F. Zhou, J. Escher, S. Cui, Bifurcation for a free boundary problem with surface tension modeling the growth of multi-layer tumors, J. Math. Anal. Appl. 337 (2008) 443–457. [45] F. Zhou, J. Wu, Stability and bifurcation analysis of a free boundary problem modelling multi-layer tumours with Gibbs–Thomson relation, European J. Appl. Math. 26 (2015) 401–425. [46] J.S. Lowengrub, H.B. Frieboes, F. Jin, Y.-L. Chuang, X. Li, P. Macklin, S.M. Wise, V. Cristini, Nonlinear modelling of cancer: bridging the gap between cells and tumours, Nonlinearity 23 (2010) 1–91. [47] G.B. Arfken, H.J. Weber, F.E. Harris, Mathematical Methods for Physicists, seventh ed., Academic Press, 2013. [48] M.G. Crandall, P.H. Rabinowitz, Bifurcation from simple eigenvalues, J. Funct. Anal. 8 (1971) 321–340. [49] A. Friedman, F. Reitich, Nonlinear stability of a quasi-static Stefan problem with surface tension: a continuation approach, Ann. Sc. Norm. Super. Pisa Cl. Sci. 30 (2001) 341–403.