Bifurcation analysis of a tumor-model free boundary problem with a nonlinear boundary condition

In this paper we study existence of nonradial stationary solutions of a free boundary problem modeling the growth of nonnecrotic tumors. Unlike the models studied in existing literatures on this topic where boundary value condition for the nutrient concentration \begin{document}$ \sigma $\end{document} is linear, in this model this is a nonlinear boundary condition. By using the bifurcation method, we prove that nonradial stationary solutions do exist when the surface tension coefficient \begin{document}$ \gamma $\end{document} takes values in small neighborhoods of certain eigenvalues of the linearized problem at the radial stationary solution.

1. Introduction. In this paper we study the following multi-dimensional free boundary problem of a 2-system of elliptic equations modeling the stationary state of a living tumor: x ∈ Ω, −∆p = g(σ), x ∈ Ω, −∂ n σ = h(σ), x ∈ ∂Ω, p = γκ, x ∈ ∂Ω, ∂ n p = 0, x ∈ ∂Ω. (1) Here Ω is the domain in R n occupied by the tumor whose boundary ∂Ω is a priori unknown and has to be determined together with the unknown functions σ = σ(x) and p = p(x), which represent the nutrient concentration in the tumor region and the pressure between tumor cells, respectively. Moreover, f , g and h are given functions, with f (σ) representing the consumption rate of nutrient by tumor cells, g(σ) representing the proliferation rate of tumor cells when the nutrient concentration is at level σ, and h being a function measuring the combined effect of the strength of nutrient supply in the host tissue of the tumor and the ability that the tumor receives nutrient from its surface. Besides, κ is the mean curvature of the tumor surface ∂Ω whose sign is designated by the convention that for the sphere it is positive, and γ is a positive constant called surface tension coefficient, which measures the strength of the surface tension of the tumor surface. Finally, ∂ n is the derivative in the direction of the outward normal n of ∂Ω. In practice n = 3, but for mathematical interest in this paper we consider the general case n 2. The above model is developed from recent work of Friedman and Lam [13], where the time-dependent version of (1.1) with f , g and h being linear functions was studied. In our recent work [20], we studied the time-dependent version of (1.1) with f , g and h being general nonlinear functions satisfying some biologically meaningful conditions (see the assumptions (A1)-(A5) below) and proved that under those conditions, (1.1) has a unique radial solution (σ s (r), p s (r), Ω s ) (with Ω s = B(0, R s )) which corresponds to the radial stationary solution of the time-dependent version of the above model. Moreover, we also proved that there exists a threshold value γ * > 0 for the surface tension coefficient γ, such that for γ > γ * this radial stationary solution is asymptotically stable under small non-radial perturbations, whereas for γ < γ * it is unstable with respect to non-radial perturbations. The aim of this paper is to study existence of nonradial solutions of the problem (1.1) by using the method of bifurcation analysis. We note that in the case that the nonlinear boundary condition −∂ n σ| ∂Ω = h(σ) is replaced by the Dirichlet boundary condition σ| ∂Ω =σ, whereσ is a positive constant reflecting the concentration of nutrient in the host tissue of the tumor, the above bifurcation problem has been intensively studied during the past twenty years, cf. [1,5,9,11,14,17,19] and references therein. Recently, Huang, Zhang and Hu [15] studied the special case of the problem (1.1) that f, g and h are linear functions, and Zhuang and Cui [6] extended their results to the case that f, g are nonlinear functions and h(σ) is still a linear function: h(σ) = β(σ −σ), whereσ is as above and β is another positive constant measuring the ability that the tumor attracts blood vessels from its host tissue. In this paper we aim at further extending the results of [6] and [15] to the general case that f, g and h are all nonlinear functions. We note that nonlinear functions are more realistic in practice due to complexity of real world environment, and linear functions are just approximations of nonlinear functions for the purpose to simplify mathematical treatment. For bifurcation analysis of other related tumor models, we refer the reader to see [7,10,12,18,21,22] and the references cited therein.
In the following, we always denote r = |x| and ω = x/r (for x = 0). As in [20], we assume that the functions f, g, h satisfy the following conditions: , and there existsσ > 0 such that h(σ) = 0; (A5)σ > σ. Note that the assumption (A5) ensures that the tumor can receive sufficient nutrient from its host tissue so that it can be alive; without this assumption the tumor will die of hunger and the problem (1.1) does not have a radial solution, cf. [23], for instance. The assumption (A1) is for mathematical simplicity. It is possible to replace this assumption with certain weaker assumptions. The assumptions (A2)-(A4) are from biological considerations, cf. [3,4,20,23] for explanation.
Here we particularly point that later on we always put Y k1 (ω) = Z n k (ω), where Z n k (ω) denotes the zonal harmonic of degree k (see [16] and Sections 4, 5 below for details). In particularly, for n = 2, by identifying the point ω = (cos θ, The main result of this paper is as follows: Theorem 1.1. Assume that conditions (A1)-(A5) are satisfied and let {γ k } k≥2 be the sequence of eigenvalues of the problem (3) mentioned above. There exists an integer k * ≥ 2 such that the following assertions hold: (1) In the case n = 2, for any k ≥ k * the problem (1) has a branch of nontrivial solutions bifurcating from the radial solution (σ s (r), p s (r), Ω s ) at γ = γ k , with the bifurcation solution having the following form: where a k is a nonzero real constant, u k (r), v k (r) are some smooth functions (see (31) for explicit expressions), and ε represents a small parameter. The second equation is the polar coordinate equation of the free surface ∂Ω.
(2) In the general case n ≥ 3, at least for every even k ≥ k * the problem (1) has a branch of nontrivial solutions bifurcating from the radial solution (σ s (r), p s (r), Ω s ) at γ = γ k , with the bifurcation solution having the following form: , v k (r) and ε are similar as above. Again, the second equation is the polar coordinate equation of the free surface ∂Ω.
Organization of the rest part is as follows. In section 2 we compute the linearization of the problem (1) and study the eigenvalue problem (3) by regarding γ as a real parameter. We shall prove existence of a sequence of eigenvalues {γ k } ∞ k=2 . In section 3 we transform the problem (1) into a functional equation in certain Banach space by using the technique of Hanzawa transformation. After making these preparations, in the last two sections we give the proof of Theorem 1.1.
2. Linearization and eigenvalues. By Theorem 1.1 of [20], we know that the problem (2) has a unique solution (σ s (r), p s (r), R s ). In this section we compute the linearization of the problem (1) at this spherically symmetric solution (σs(r), ps(r), Rs) and determine eigenvalues of the linearized problem (3).
First of all, as in [5,6,11,14,15,19], by considering a solution (σ, p, Ω) of (1) of the form , with a small real parameter ε and using a similar argument as in the proofs of Lemma 3.1 of [5], Lemma 2.1 of [6] and Lemma 5.2 of [20], we easily obtained the following preliminary result: Lemma 2.1. The linearization of the problem (1) at the radial solution (σ s (r), p s (r), Ω s ) is the system (3).
Comparing the system (3) with the system (1.4) of [6], we see that by letting β = h σ s (R s ) , (3) coincides with (1.4) of [6]. Thus, all results obtained in [6] concerning the system (1.4) there apply to the present system (3) by replacing β there with h σ s (R s ) . For completeness, in what follows we present a skeleton of the process of solving (3).
Note that since all known functions in (3) are smooth, by some standard results in the regularity theory of elliptic partial differential equations we see that the solution (ϕ, ψ, ρ) of (3) belongs to C ∞ (Ω s ) × C ∞ (Ω s ) × C ∞ (S n−1 ). Hence we can expand ϕ, ψ, ρ via spherical harmonics {Y kl (ω)} as follows: Substituting these expressions into (3) and comparing coefficients of every Y kl (ω), we get Note that we naturally have the following boundary conditions: u kl (0) = v kl (0) = 0 for all k ∈ Z + and l = 1, 2, · · · , d k , which arise from smoothness of ϕ(x) and ψ(x) at x = 0.
As in [5,6], for each k ∈ Z + = {k ∈ Z : k 0} we use the notationū k (r) to denote the unique solution of the following initial value problem: It is not difficult to see that the unique solution of (6) 1 subject to the boundary condition (6) 3 is given by u kl (r) = µ k c kl r k u k (r), where .
Again as in [5,6], for each k ∈ Z + we use the notationv k (r) to denote the unique solution of the following initial value problem: In the same way, we can easily obtain the unique solution of (6) 2 subject to the boundary condition (6) where Substituting the above expressions (7), (8) of v kl into (6) 5 , we see that the system (6) reduces into the following equation for c kl : where (k = 0, 1, 2, · · · ). Hence, (ϕ, ψ, ρ) given by (5) is a nontrivial solution of (3) if and only if there exist k, l such that (9) has a nontrivial solution c kl . Clearly, this is possible if and only if there exists a nonnegative integer k such that γ is a solution of the equation a k (γ) = 0. By (10), for every k = 2, 3, · · · , the above equation has a unique solution which we denote as γ k , having the following expression: where In fact, it is clear that the following relations hold: Similarly as in Lemma 4.1 of [23], for k = 0, 1 we have, for any γ > 0, In summary, we have proved the following lemma: Since by letting β = h σ s (R s ) , (3) coincides with (1.4) of [6], by Lemma 2.3 of [6] we have: Lemma 2.3. There exists an integer k * 2 such that γ k is strictly monotone decreasing for k k * .

Reduction into a functional equation.
In this section we reduce the problem (1) into a functional equation of the form F(ρ, γ) = 0, where ρ represents a function defined in S n−1 describing the free boundary ∂Ω, and [ρ → F(ρ, γ)] is an operator mapping some neighborhood of the origin of C m+µ (S n−1 ) into C m−3+µ (S n−1 ) (for each fixed γ > 0), where m is a positive integer not less than 2 and 0 < µ < 1. We note that unlike the problem treated in [5,6] where both boundary conditions for σ and p are linear, in the present situation we meet with a nonlinear boundary condition for σ. Fortunately, the elliptic boundary value problem with a nonlinear boundary condition has been successfully treated in [20]. We cite the result of [20] that we need as follows: Lemma 3.1. Assume that Conditions (A1)-(A4) are satisfied and let γ > 0 be fixed. Let m, µ be as above. For any bounded domain Ω in R n which is C m+µdiffeomorphic to the unit ball in R n , the problem has a unique solution σ ∈ C m+µ (Ω). Moreover, the map Ω → σ is smooth.
To reduce the problem (1) into an equation of the form F(ρ, γ) = 0 as mentioned above, we need the so-called Hanzawa transformation which we introduce as follows. First we choose an even function χ ∈ C ∞ (R) such that it satisfies the following conditions: Here δ represents a small number such that 0 < δ < 1 In what follows, as before we use the notation (r, ω) to denote the polar coordinate of a point x in R n \{0}, so that r = |x| and ω = x/r. For ρ ∈ O δ , we denote Using the above notations, we see that Ω s = Ω ρ=0 and ∂Ω s = ∂Ω ρ=0 . Given ρ ∈ O δ , the Hanzawa transformation Ψ ρ : R n → R n is defined as followed: x, otherwise.
From (15), it is not difficult to verify that Ψ ρ is a C m+µ -diffeomorphism of R n onto itself. Besides, it is also easy to see that and Ψ ρ (Ω s ) = Ω ρ and Ψ ρ (∂Ω s ) = ∂Ω ρ .
Now we substitute the above expression into the equation (20) 5 . It follows that, by letting where π * denotes the pull-back of the natural projection π : S n−1 → ∂Ω s , i.e., π(ω) = R s ω for ω ∈ S n−1 , the problem (20) reduces into the following equation containing the unknown ρ only: From (19), (21) and (22) we see that the following relation holds: In summary, we have proved the following preliminary result: We denote by D ρ F(0, γ) the Frëchet derivative of F(ρ, γ) at ρ = 0. To compute D ρ F(0, γ), we need to linearize the system (20) at ρ = 0 and next reduce it into a functional equation as we have just done for the system (20). Since (20) comes from (1) by making the Hanzawa transformation, a simple analysis shows that the linearization of the system (20) is not other but just the system (3). Hence, from the analysis in Section 2 we obtain the following preliminary result: Lemma 3.4. For any γ ∈ R, the operator D ρ F(0, γ) is a Fourier multiplier in the sense that it has the following expression: where c 0 is as in (13) and a k (γ) (k = 0, 1, 2 · · · ) are as in (10).
From the expressions (12) and (25) we immediately obtain the following result: Corollary 3.5 The following relation holds: if γ = γ k for some k 2, 4. Bifurcation in case n = 2. In this section we present the proof of the assertion (1) of Theorem 1.1. We shall use the Crandall-Rabinowitz bifurcation theorem to prove the desired result. For reader's convenience, let us first recall this theorem as follows: Theorem 4.1 (Crandall-Rabinowitz [2]). Let X, Y be real Banach spaces and let F (x, λ) be a C p (p 3) map of a neighborhood (0, λ 0 ) in X × R into Y . Suppose that the following conditions are satisfied: . 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 -curves Γ 1 and Γ 2 which intersect only at the point (0, λ 0 ); Γ 1 is the curve (0, λ) and Γ 2 can be parameterized as follows: Throughout this section, we shall use (r, θ) to denote the polar coordinate system of R 2 . Hence we always identity the point ω = (cos θ, sin θ) with θ ∈ [0, 2π). In this way, instead of Y kl (ω) used before we use a different notation Y kl (θ) for spherical harmonics on the circle S 1 , and correspondingly we use the notation ∆ θ to replace ∆ ω . Recall that for n = 2 we have d 0 = 1 and d k = 2, k = 1, 2, 3, · · · , and the standard basis of spherical harmonics is as follows: Recall also that in the polar coordinate system the Laplacian operator in R 2 has the following expression: Proof. Theorem 1.1 (1): Let k * be as in Lemma 2.3, and let k be an arbitrary integer such that k k * . For an integer m 2 and a number µ ∈ (0, 1) we define X m+µ k = closure of span{cos(jkθ) : j = 0, 1, 2, · · · } in C m+µ (S 1 ); and we simply write X and Y for the spaces X m+µ k and X m−3+µ k , respectively. For any ρ ∈ O δ , we know that the outward normal derivative of the surface r = R s +ρ(θ) is given by and the mean curvature of this surface has the expression By using Lemma 5.1 of [5] By using these facts and some similar argument as in the proof of the assertion A ∈ C ∞ (C m+µ 2π,l,0 , C m+µ 2π,l,0 ) in Lemma 5.2 of [5] we see that the following relations hold: Using these facts and once again Lemma 5.1 of [5] we finally obtain the following relation: In what follows we use the same notation F to denote its restriction to (O δ ∩X)×R. Since ρ = 0 corresponds to the unique radial solution (σ s , p s , Ω s ) of the problem (1), from the reduction of the equation (23) we see that F(0, γ) = 0 for all γ > 0.

5.
Bifurcation in case n ≥ 3. In this section we give the proof of the assertion (2) of Theorem 1.1.
In what follows we redenote θ 1 as θ. We use the notation Z n k (ω) to denote the zonal spherical harmonics of degree k on S n−1 with pole e = (1, 0, · · · , 0), i.e., where C λ k represents the ultraspherical polynomial of degree k and index λ, i.e., and c nk is the normalization factor. From (32) and (33) we see that the zonal spherical harmonics Z n k (ω) are independent of the variables θ 2 , · · · , θ n−1 , and they are polynomials of cos θ with all exponents being even numbers when k is an even integer. Therefore, for any ρ ∈ X, ρ is a function of the variable θ only and satisfies the conditions ρ(θ) = ρ(π − θ) and ρ (0) = ρ (π) = 0.
A significant difference between the cases n ≥ 3 and n = 2 is that in the latter case the closure of the linear space spanned by all Z 2 jk (θ) = cos(jkθ) (j = 0, 1, 2, · · · ) in C m+µ (S 1 ) is an algebra for any k ≥ 2. In the case n ≥ 3, however, this is no longer true. Hence there are some changes in the following deduction.