Turing type instability in a diffusion model with mass transport on the boundary

Some reaction-diffusion models describing the cell polarity are proposed, where the system has two independent variables standing for the concentration of proteins in the membrane and the cytosol respectively. In this article we deal with such a polarity model consisting of one equation on a unit sphere and the other one in the ball inside the sphere. The two equations are coupled through a nonlinear boundary condition and the total mass is conserved. We investigate the linearized stability of a constant steady state and provide conditions under which a Turing type instability takes place, namely, the constant state is stable against spatially uniform perturbations on the sphere for all choices of diffusion rates, while unstable against nonuniform perturbations on the sphere as the diffusion coefficient of the equation on the sphere becomes small relative to the one in the ball.


1.
Introduction. There are a huge number of mathematical models in the fields of cellular and molecular biologies to understand complex phenomena exhibited by organisms in vivo. One of the interesting phenomena is the polarity of cells (see [8,16,23,24], and references therein). To understand mechanisms of the polarity, reaction-diffusion models have been proposed by a number of authors in the works [11,15,17,29,30,31,32]. Most of these studies are able to account for the emergence of polarity, in which the reactions of active proteins on the membrane and inactive ones in the cytosol are taken into account and these two types are transformed from one to the other by some action of enzymes. In addition, the active proteins and the inactive ones diffuse in their own domains, slowly in the membrane while faster in the cytosol.
As for rigorous mathematical studies in the cell polarity models we refer to [12,13,15,14,21,22,23,24,25,26]. These papers are concerned with a class of reaction-diffusion systems in the bulk under periodic boundary conditions or the Neumann boundary conditions on the boundary. One of the interesting aspects of the model for the cell polarity is the emergence of wave patterns by a Turing type instability (see [11] and [30]). In fact, it is easily verified that the constant solution, say (u, v) = (u, v), of the model equations under periodic boundary conditions is stable against spatially uniform perturbations, but unstable for spatially nonuniform perturbations provided that the conditions are met for a wave number k (see also [26]). We note that u and v are supposed to stand for the active and inactive proteins on the membrane and in the cytosol, respectively. However, in (1) as it stands, u and v are both defined in the entire bulk interval (0, 2π). The only distinction between the active one (u) residing on the membrane and the inactive one (v) residing in the entire bulk interval is modeled as the difference between the respective diffusion rates. Hence, u moves slower than v, so that d 1 < d 2 is assumed. If one distinguishes the domains of membrane and cytosol, it is more realistic to consider model equations in each domain, the one corresponding to the cytosol, and the other on the boundary supposed to be the inner surface of the membrane. In this setting the linearized stability analysis around the constant steady state solution, together with some numerics displaying spatial patterns, is reported in [31] and [32] (see also [17]).
Taking the model equations in [30] and those of [31] and [32] into account, instead of (1), we consider the following system and develop an instability analysis around a constant steady state. In the system (3)-(4)- (5), Ω is a bounded domain in R n (n = 2, 3) with smooth boundary ∂Ω =: Γ and ∆ stands for the Laplacian in Ω while ∆ Γ stands for the Laplace-Beltrami operator on the boundary. We assume Γ is a closed manifold at first, and later we specify it as the unit ball in R n . As for the relation between our model (3)-(4)- (5) and the models in [31] and [32], refer to §4 (Appendix I) of the present article. On the well-posedness of the initial value problem in various function spaces, consult [1]. In addition one can refer to [33] on the global existence of the solution ( [18], [19] and [20] are helpful, as well). We are interested in the linearized stability of a constant equilibrium of (3)-(4)- (5). We first observe that this system has a conservation property by TURING TYPE INSTABILITY WITH MASS TRANSPORT ON THE BOUNDARY   3815 The system has a positive constant steady state (u, v) = (u c , v c ), defined as a solution of h(u, v) = 0, |Γ|u + |Ω|v = m.
We always assume which biologically means that the protein represented by v in the cytosol is transformed on the boundary (cf. (5)) by the action of some enzymes and contribute to increase the growth rate of the protein u residing on the membrane (cf. (3)). Moreover, in order to avoid the critical situation in which the two curves of (7) are tangent at (u, v) = (u c , v c ), it is natural to set Here we assume that the sign of the left hand side is negative: Geometrically, (9) means that the slope −α/β of the tangent line to the curve h(u, v) = 0 at (u c , v c ), is larger than the slope −|Γ|/|Ω| of the straight line |Ω|v + |Γ|u = m. Note that the condition (9) is met for any α ≤ 0. Although we have no explanation of a biological significance on this assumption, as it will turn out in §2, it is the condition to ensure the stability against spatially uniform perturbations on the boundary Γ of the constant steady state (u c , v c ), so it corresponds to the former condition of (2) for the equations (1) (see Lemma 2.5 and its proof). Under the conditions (8) and (9), we deal with the linearized eigenvalue problem: in the space For any non-zero eigenvalue, the corresponding eigenfunctions (φ, ψ) satisfy the constraint. This is easily proven by integrating the equations (10) and (11) on Γ and on Ω, respectively, and using the boundary condition (12). In the absence of the constraint Γ φ dx Γ + Ω ψ dx = 0, the eigenvalue problem (10)-(11)-(12) with λ = 0 admits the nontrivial solution being the constant vector (φ, ψ) = c(β, −α) (c = 0), which is tangent to the curve defined by h(u, v) = 0 at (u, v) = (u c , v c ). However, by the assumption (9), this vector never satisfies the constraint, hence, it is not an eigenvector of the problem. When α ≤ 0, the constant steady state (u c , v c ) is stable for general domains Ω, which is stated as follows: (8) and (9). If α ≤ 0, then all the eigenvalues λ of (10)-(11)- (12) are negative for any bounded domain Ω with smooth boundary.
Therefore, Turing instability never occurs when α ≤ 0, hence, our main concern in what follows is the case α > 0. We prove Proposition 1 in the next section.
In order to go further in our (in)stability analysis, we specify the domain to be Ω = B n := {x ∈ R n : |x| < 1} and Γ = S n−1 := ∂B n (n = 2, 3). Denote by ∆ S the Laplace-Beltrami operator on S n−1 , that is, ∆ S in the polar coordinate is given by k=0 denote the complete orthogonal system of eigenspaces of −∆ S corresponding to the eigenvalues µ k = k 2 + (n − 2)k, namely, for k ≥ 1, dim V k = 2 for n = 2 (dim V k = 2k + 1 for n = 3), and We expand an eigenfunction (φ, ψ) for (10)-(11)- (12) in terms of eigenspaces {V k (ω)} ∞ k=0 and unknown radial functions W k (r) and unknown constants ξ k as follows.
is met, then any ( , Γ)-mode complex (non-real) eigenvalue has negative real part.
We may regard Theorem 1.1 as a Turing instability theorem for our system (3)-(4)-(5) at the constant equilibrium (u c , v c ). This summarize as a corollary. Corollary 1. When the diffusion rate d 1 is not small (i.e., d 1 > α/µ 1 > 0), the constant equilibrium is stable for any perturbation, in particular, for spatially uniform perturbation on Γ. When the diffusion rate d 1 > 0 of u-component gets smaller and crosses one of the values 0 < α/µ < α/µ −1 < . . . , α/µ 2 < α/µ 1 successively, and moreover the pair (d 1 , d 2 ) crosses the hyperbola from right to left, then the constant equilibrium loses stability in the direction of ( , Γ)-mode perturbations.
Proof. Thisl is proven by the fact that the largest ( , Γ)-mode real eigenvalue λ crosses λ = 0 from left to right.
Remark 2. Now, let us compare our conditions in Theorem 1.1 and the conditions for the Turing type instability for (1). The conditions in (2) read (for α > d 1 µ ) α < β and for some wave number (we replaced k by ). Therefore, in our system (3)-(4)-(5) the condition (16) allowing the Turing type instability is less restrictive on the choice of d 1 and d 2 than the second condition in (19) except for = 1.
Remark 3. We see that β > α in (19) implies that the second condition of (19) is not met for the equal diffusion case: d 1 = d 2 . However, our condition (9) in Theorem 1.1 accommodates the situation α > β. In fact, (9) reads which allows the situation α > β if β is taken, for example, as follows: For the equal diffusion case d : If we set = 1, then holds. In view of µ 1 = 1 (n = 2) and µ 1 = 2 (n = 3), we easily find that the condition (16) is fulfilled for d 1. Theorem 1.1 (2) applies and the constant steady state is unstable in the direction of (1, Γ)-mode perturbations for sufficiently small d > 0 by (16) even if the diffusion coefficients are the same. Moreover, it is possible to take d 1 slightly larger than d 2 so that (16) is still valid.
More generally, if β < µ α, we could obtain the Turing instability in the direction of ( , Γ)-mode perturbations for small equal diffusion However, because of the standing hypothesis (9), the condition α < nβ must be satisfied. Then, only the case = 1 is realized.
There are two possible ways to further reduce the eigenvalue problem (10)-(11)-(12). One way of reduction is this: Consider the following Dirichlet problem for the standing for the eigenvalues of the Laplacian under the zero Dirichlet boundary condition on Γ. Then there is a unique solution f (·; λ, g) ∈ H 2 (Ω). Define the Dirichlet-Neumann map for (∆, Ω) by We can verify that T (λ)g is analytic in λ ∈ C\σ Dir ( ) and for each λ ∈ C\σ Dir ( ), the map T (λ) : H 3/2 (Γ) → H 1/2 (Γ) is a bounded linear operator (see [4], and [3] is also helpful). Utilizing T , we can reduce the eigenvalue problem (10)- (12) to In this direction the linear stability analysis for a related model in the specific domain Ω = B 3 is done in [31] and [32] though there is no explicit formulation by the Dirichlet-Neumann map (see also the work in [27] for a related problem).
The second way to further reduce the problem (10)- (12) is to recast it as follows: In our present study, rather than identifying the parameter values at which the stability properties switch, we reveal how the eigenvalues qualitatively change as the parameters varies. To this end we use a variational approach and some properties of Dirichlet-Neumann map for the analysis. In what follows, we display a clear mathematical argument for the proof of the instability. We finally remark that the possibility of complex eigenvalues with positive real part is discussed neither in [31] nor in [32]. In our study as seen in Theorem 1.1 we provide conditions on nonexistence of complex eigenvalues with positive real part.
When we had almost completed the manuscript, Professor Muratov kindly sent us a recent paper [7] which deals with a system almost similar to our system (3)-(4)-(5). Given a specific kinetic term and reducing the system to a single equation, he and his collaborators constructed spherical cap-solutions and analytically determined the linear stability of these non-uniform steady states. In their model, it is easy to verify that constant steady states are always linearly stable, while in our general case, constant steady states may change stability properties according to the diffusion rates (d 1 , d 2 ) through Turing mechanism.
In §2, we will provide proofs of our results. In §2.1 we prove Proposition 1 and in §2.2 we study stability against (0, Γ)-perturbations. §2.3 is devoted to the investigation of the eigenvalues corresponding to (k, Γ)-perturbations (k ≥ 1). By using the variational characterization of (k, Γ)-eigenvalues, we discuss complex eigenvalues and complete the proof of Theorem 1.1. In §3, to illustrate, we apply the main results to specific examples. In order to relate our study to the problem of pattern formation, we derive a shadow system by reducing original system of two equations into a system of a single equation on Γ and state a spatial pattern for a specific case in §4. We conclude the paper by § §5, 6 and §7 of three Appendices in which technical results are proven for the sake of completeness.
In this paper, we have proven that at the constant steady state (u, v) = (u c , v c ) the system (3)-(4)-(5) exhibits Turing type instability when the conditions (8)- (9) are satisfied. This was established in Theorem 1.1 and Corollay 1 by the careful analysis of eigenvalues of the linearized eigenvalue problem (10)-(11)- (12). The next step may be to develop a weakly nonlinear analysis around the destabilization threshold in order to determine the branch of bifurcating solutions. We are aiming to perform computations to establish such weakly nonlinear analysis. However, it turned out that even formal computations are quite lengthy and making those formal computations rigorous require substantial amount of work. We therefore leave such investigations as future projects and do not pursue in this paper. Proof. Multiply (10) by φ (the complex conjugate of φ) and integrate over S, and then multiply (11) by ψ and integrate over Ω. Then, we obtain the following two identities: When α = 0, it is clear that λ is real and negative. Let α < 0. Multiply the first identity by α and multiply the second identity by β and subtract the latter from the former. We obtain the identity This implies λ is real and negative.
From this point on, therefore, we are concerned with α > 0 throughout.
Proof. Suppose that W and λ are a complex-valued function and a complex number respectively. We assume that the imaginary part of λ is nonzero and then derive a contradiction. Since Dividing, then, the second equation above by λ 2 = 0 and plugging the resulting identity into the first equation, we obtain the identity which implies that (23) cannot hold for nontrivial W (r). Consequently, for the pair (W, λ) to be (0, Γ)-eigenpair, the imaginary part of λ has to vanish, λ 2 = 0, and hence, W must be real valued and the eigenvalues λ must be real.
We have the next lemma.
Proof. The continuity of a simple eigenvalue in the parameter is known (for instance, see [6]). The recent work of [13] could also be of help, where the simplicity of the eigenvalue is not required for the continuity. In order to prove the monotonicity of σ (0) (s)/s, we differentiate the numerator of R[Φ, s] divided by s with respect to s; Indeed, this leads us to the monotonically decreasing ofσ (0) (s)/s (s = 0) by utilizing the Rayleigh quotient (for the details of the argument, consult [13]). In view of σ (0) (0) = 0 and the next lemma,σ (0) (s)/s is defined at s = 0 so that it is continuous at s = 0. This concludes the proof.
(ii) If α > d 1 µ k and if then there is no positive eigenvalue.
then there exists a positive eigenvalue.
To this end we need some preparation. Given η 1 > 0, we letw(r;s, η 1 ) be a solution to In terms of the Dirichlet-Neumann map we have that T k (s)[η 1 ] :=w r (1;s, η 1 ), and then the problem (25)- (27) is reduced to For λ = α − d 1 µ k , the first equation reads Inserting this into the second equation yields Since the operator is linear in η 1 , we consequently obtain the equation in which we simply write T k (λ/d 2 )[1] as T k (λ/d 2 ) = T k (λ/d 2 )[1] being nothing but the eigenvalue of the Dirichlet-Neumann map T (λ/d 2 ) corresponding to the eigenfunction Ψ k ∈ H 3/2 (Γ). We restrict our attention to real λ, and rewrite the equation (32) by replacing λ by s; A real eigenvalue λ of (25)- (27) is characterized as the solution s = λ of (33). For brevity of notation, we define the differential operator A k by and we simply denote byw(r;s) a unique solution of This means that T k (s) =w r (1;s). The next lemma is useful.
The proof of this Lemma will be given §7.1 (Appendix III).
Remark 4. We note that if α − d 1 µ k > 0, then the condition (31) is equivalent to In order to handle the case of complex eigenvalues, we use the following technical inequality.

3.
Applications. As mentioned in Introduction the model equations (1) are proposed in [11] and [30] with the following specific functions h(u, v): where a 1 , a 2 and b are positive constants. For these cases we verify the conditions in Theorem 1.1. To this end we show that there are positive equilibria at which α = h u > 0, β = h v > 0 and (9) hold. Then we obtain the instability stated in the theorem for sufficiently small d 1 .
Next consider the case (46). We solve the equations Put w = u + v. Then (48) turns to be By g 2 (w) = 1 − a 2 w (a 2 w + 1) 3 , and κ > 1, there is a unique solution (w c , v c ) of (49) if (1/a 2 ) < m/κ|Ω| = m/|Γ|. We see The last inequality is realized by taking m large, since w c → ∞ as m → ∞, hence, g 2 (w c ) → 0 for large m 1. In this case, our system (3)-(4)-(5) with nonlinearity h in (46) is monostable with respect to (0, Γ)-perturbations. 4. Shadow systems. In this section we discuss the shadow system obtained by taking the limit d 2 → ∞ in (3) with the constraint Plug (51) into (50), then we obtain Thus, the limiting problem turns out to be the scalar equation (52) for u with a nonlocal term. We note that the limiting equation of ξ is equivalent to the equation (51). Indeed, differentiating (51) with respect to t yields Here, we consider the specific nonlinearity h(u, v) = −g(u) + v of (45) in the domain Ω = B 2 . Then, the limiting problem (52) have a more concrete form: and the stationary problem of (52) is the Euler-Lagrange equation of E(u).
As seen in the previous section the constant solution becomes unstable for an appropriate m. Then it is known that there exist a stable unimodal stationary solution to (54) and no stable solutions with other spatial patterns (see [26]).
5. Appendix I. We state how our model is related to those in [32]. We also verify that our conditions are consistent with those in [32] in the stability for the spatially uniform perturbation on Γ. The reduction argument below is mathematically not rigorous, but it supports that our model is worth studying.
The model equations with the same diffusion rate on Γ in [32] reads Putting U := u + v, we obtain the equations for (u, U, V ) as In order to simplify the system we assume that the quasi steady state approximation works for the first equation, and that f (u, U − u) = 0 is solved as u = ϕ(U ). We insert it into q as h(U, V ) := q(ϕ(U ), U − ϕ(U ), V ) and get to This system is nothing but (3)-(4)- (5). We note that if q(u, v, V ) = h(u + v, V ), then (60)-(62) immediately follows from (56)-(57)-(58)-(59) by putting U = u + v, namely, the quasi steady sate assumption is not necessary. We verify that the assumptions (8) and (9) are consistent with those in [32]. First β = h v (u c , v c ) > 0 of (8) coincides with the assumption q V > 0 in (2.5) in [32]. We consider the next condition of Proposition 2.1 in [32] 1 which is a necessary and sufficient condition for the stability of the constant solution under spatially homogeneous perturbations on Γ = S 2 , though in [32] the possibility of complex eigenvalues with positive real part is not discussed. From the equation , at the constant steady state. Then compute α, β as Hence the condition (8) Following the computation of [32], we see Since they assume f v > 0, q v ≤ q u and q v < 0 (see (2.1) and (2.5) of [32]), the condition (63) yields f v − f u > 0. In the sequence, (63) implies (8) This concludes the equivalence of (8) and (63). We remark that the quasi-steady state approximation does not work when the Turing instability is discussed, therefore it is reasonable that the conditions of Theorem 1.1 do not coincide with those in [32]. We emphasize that they don't discuss the case of complex eigenvalues.
6. Appendix II. Although the differentiability for a simple eigenvalue in a parameter is well studied, for completeness of the present article we give a proof of the differentiability of the eigenvalue σ (0) (s) and a corresponding eigenfunction W 0 (r; s) of (24) with respect to s. Moreover, we only consider it in a neighborhood of s = 0 and set d 2 = 1 for the sake of clarity.
We look for a solution in the form where W c is any positive constant. Plug it in the equation (24). Then =σW c (r 2 − 1)/2n + sσ There are ε 1 > 0 and M > 0 such that F maps X into itself provided |s| < ε 1 . Moreover, we can see that there is a constant L > 0 such that holds. This implies that there is a unique fixed pointŴ (r;σ, s) of F for every small |s|. By the continuous differentiability of F(Ŵ 1 ,σ, s) in the arguments and the uniform contraction we see that the fixed pointŴ (r;σ, s) is continuously differentiable in (σ, s) (for instance, see [10]). Inserting this solutionŴ (r;σ, s) into (65) yields (r;σ, s)r n−1 = 0.
7. Appendix III. In this section, we will give the proof of technical Lemma 2.6 used in §2 together with a variational approach to Porposition 2 (iii).
Applying this inequality to the above inequality yields kw(1) 2 ≤ 1 0 w 2 r r 2 dr + k(k + 1) This concludes the proof of Lemma 2.7.

7.3.
A variational treatment of Proposition 2 (iii). We are able to show that a variational approach works for the proof of Proposition 2 (iii). We introduce the following eigenvalue problem with a parameter s = α − d 1 µ k : Let ω (k) (s) be the largest eigenvalue ω (k) (s) of (67). Then the following variational characterization works: We easily see that ω (k) (s) is continuous in s, which is proved by utilizing the variational characterization (this is well known, but the readers may refer to [13]). From the next lemma, Proposition 2 (iii) immediately follows. Indeed, ω (k) (s * k ) in the lemma gives the desired eigenvalue.