Original article Biomath 3 (2014), 1411111, 1–12 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Regular and Discontinuous Solutions in a Reaction-Diffusion Model for Hair Follicle Spacing Peter Rashkov FB Mathematik und Informatik Philipps-Universität Marburg 35032 Marburg, Germany Email: rashkov@mathematik.uni-marburg.de Received: 21 July 2014, accepted: 11 November 2014, published: 24 November 2014 Abstract—Solutions of a model reaction-diffusion system inspired by a model for hair follicle initiation in mice are constructed and analysed for the case of a one-dimensional domain. It is shown that all regular spatially heterogeneous solutions of the problem are unstable. Numerical tests show that the only asymptotically stable weak solutions are those with large jump discontinuities. Keywords-dynamical systems; reaction-diffusion equation; stationary solutions; weak solutions I. INTRODUCTION A parabolic reaction-diffusion system is pro- posed in [14] to model the WNT signaling path- way in primary hair follicle initiation in mice. The authors in [14] use a modified version of the well-known activator-inhibitor (Gierer-Meinardt) model [3], [4] with saturation and without source terms. An important characteristic of the model is that both species share the same (up to scaling) non-linear production term for both activator and inhibitor. A modified version of this model was stud- ied in [12] as a proxy to reduce the parameter complexity and to capture the dynamics of the original model. Global existence of solutions of both the original and the modified systems was demonstrated by estimating time-independent up- per bounds for the solutions. A parameter space analysis indicated the range of the existence of Turing patterns. It is demonstrated that het- erogeneous solutions arise not only because of diffusion-driven instability, but also due to conver- gence to far-from-equilibrium solution branches. This short note compares stationary solutions in the singularly perturbed problem (letting the inhibitor’s diffusion rate tend to 0) and the reduced problem (setting the inhibitor’s diffusion rate equal to 0) based on the modified equations from [12] for the case of a one-dimensional domain. Stability of the stationary solutions is analysed. We show that all strictly positive, spatially hetero- geneous, regular solutions of the reduced problem are unstable. Furthermore, for some parameter Citation: Peter Rashkov, Regular and Discontinuous Solutions in a Reaction-Diffusion Model for Hair Follicle Spacing, Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 1 of 12 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... values, the spatially homogeneous solution is also unstable. The only asymptotically stable solutions are weak solutions where the activator exhibits jump discontinuities. This work is organised as follows: first we define the model system and its reduced variant as a coupled ODE-reaction-diffusion system. Then we restrict our attention to a one dimensional domain and convert the problem to an auxiliary two-point boundary value problem. Energy meth- ods are employed to construct the regular and weak stationary solutions. Finally we establish the stability properties of the different solutions. II. THE MODEL PROBLEM Let Ω ∈ Rn be a bounded domain with suf- ficiently regular boundary ∂Ω. Consider the fol- lowing problem that describes the spatio-temporal dynamics of two interacting species ut = � 2∆u + f(u,v), vt = d∆v + g(u,v), ∂xu(·,x) = ∂xv(·,x) = 0, x ∈ ∂Ω, (1) with nonlinearities f,g given by f(u,v) = ρu u2 v(1 + κu2) −µuu, g(u,v) = ρv u2 v(1 + κu2) −µvv. (2) The functions u = u(t,x),v = v(t,x) describe the concentrations of the species at x ∈ Ω for time t > 0. The initial conditions u(0, ·),v(0, ·) are sufficiently smooth so that the second derivatives in space are well-defined. The model parameters ρu,ρv,µu,µv,κ have the following physical interpretation. κ is saturation parameter for the production law for u and v, which is scaled respectively by ρu,ρv. µu,µv denote the decay rates of u and v. The diffusion constants �,d describe the diffusion speeds in the domain Ω. The model equations are based on the equations proposed in [14] to model hair follicle spacing in mice. Of particular interest are the properties of the non-negative stationary solutions of (1), i.e. those pairs (u,v) such that ut = vt = 0. These are those pairs (u,v) solving the problem of two coupled elliptic PDEs 0 = �2∆u + f(u,v), 0 = d∆v + g(u,v), ∂xu(·,x) = ∂xv(·,x) = 0, x ∈ ∂Ω. (3) A. Diffusion-driven instability The mechanism of diffusion-driven (or Turing) instability has been used used in mathematical and biological models to motivate the emergence of patterns and forms (spatial heterogeneities) in development processes. The classical form of the mechanism is described by a reaction-diffusion model system with two morphogens that react and diffuse in the domain producing heterogeneous spatial patterns [10]. Let us recall the conditions for diffusion-driven instability of a steady state (û, v̂) of (3). The Jacobian of the reaction-kinetic system ut = f(u,v),vt = g(u,v) evaluated at this steady state is J = ( fu fv gu gv ) . (4) From the definition of diffusion-driven instability, in the absence of diffusion d = 0, the steady state (û, v̂) must be locally unstable to spatially inhomogeneous perturbations ũ(t,x) = u(t,x) − û, ṽ(t,x) = v(t,x) − v̂, but locally stable to spa- tially homogeneous perturbations ũ(t) = u(t) − û, ṽ(t) = v(t) − v̂. Under the Ansatz ũ(t,x) = ueikxeλt, ṽ(t,x) = veikxeλt, where u,v are scalars, ũ(t,x), ṽ(t,x) will be growing in time t if the eigenvalue λ associated to the wave number k > 0 satisfies Reλ > 0. For the spatially homogeneous perturbation (k = 0) the eigenvalue λ associated to the wave number k = 0 satisfies Reλ < 0. Both conditions can be written in terms of the dispersion relation M(λ,k) relating the eigen- value λ to the wavenumber k. The conditions for Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 2 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... diffusion-driven instability can be formulated as follows. First, we require the equation M(λ, 0) = 0 to have solutions with Reλ < 0. Second, there must exist at least one k0 > 0 such that the equa- tion M(λ,k0) = 0 has a solution with Reλ > 0. In particular, this implies that the diffusion constants must be different d 6= �2. In particular biological applications the diffu- sion constant �2 may be so small to be negligible or u may not diffuse at all. The reduced problem with � = 0 may also exhibit diffusion-driven insta- bility, but the above conditions may have different significance. In particular, the properties of the stationary solutions cannot be derived from a linear stability analysis of the spatially homogeneous steady state because these solutions are far-from- equilibrium solutions. We remark that spatial heterogeneities arise also in reaction-diffusion systems where the nonlinear- ities do not even allow the existence of a spatially homogeneous steady state (û, v̂) [12]. In this note, our attention is restricted to the one- dimensional case Ω = [0, 1]. For simplicity we set ρu = ρv = 1 in (2). We begin by providing some a priori estimates for the solutions u,v > 0 of the stationary prob- lem (3). B. A priori estimates Throughout the rest of the discussion we let Ω = [0, 1]. Lemma 1. Let f,g be given by (2). Assume that the pair (u,v) with positive u,v ∈ C2(Ω) solves (3). Then max Ω u ≤ e √ µu/� min Ω u, max Ω v ≤ e √ µv/d min Ω v. Proof: For shortness we shall prove this for v. The computations for u > 0 are analogous. Since we are looking for a solution (u,v) > 0, we can divide both sides of the equation for v in (3) by v and integrate over Ω, µv = d ∫ 1 0 ∆v v dx + ∫ 1 0 u2 v2(1 + κu2) dx. Integration by parts using the boundary condition ∂xv(0) = ∂xv(1) = 0 gives µv = d ∫ 1 0 ( ∂xv v )2 dx + ∫ 1 0 u2 v2(1 + κu2) dx, whence ∫ 1 0 ( ∂xv v )2 dx ≤ µv d . Choose xmin,xmax ∈ Ω such that v(xmin) = min Ω v, v(xmax) = max Ω v. Then log max v − log min v = ∫ xmax xmin ∂xv v dx ≤ |xmin −xmax|1/2 · ∣∣∣∣∣ ∫ xmax xmin ( ∂xv v )2 dx ∣∣∣∣∣ 1/2 ≤ √ µv d , implying max v ≤ e √ µv/d min v. Theorem 1. Let (u,v) ∈ C2(Ω) solve (3). Then there exist monotone functions ψ↑,ψ↓,α such that α,ψ↑ are monotone increasing, ψ↓ is monotone decreasing, and ψ↑(�) < u(x) < ψ↓(�), (5) α(min Ω u) ≤ v(x) ≤ α(max Ω u). (6) The functions ψ↑,ψ↓,α are independent of d. Proof: Solving the equation g(u,v) = 0 (2), v can be expressed in terms of u as v = α(u) :≡ ( u2 µv(1 + κu2) )1/2 . α(u) is a monotone increasing function of u. Observe that for a fixed u g(u,v) < 0 iff v > α(u) and g(u,v) > 0 iff v < α(u). Choose xmin,xmax ∈ Ω such that v(xmin) = min Ω v, v(xmax) = max Ω v. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 3 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... Because ∆v(xmin) ≥ 0, ∆v(xmax) ≤ 0 the equal- ity d∆v + g(u,v) = 0 implies g(u(xmin),v(xmin)) ≤ 0, g(u(xmax),v(xmax)) ≥ 0. Therefore, min v = v(xmin) ≥ α(u(xmin)), min v = v(xmax) ≤ α(u(xmax)), and due to monotonicity of α we obtain min Ω v ≥ α(min Ω u), max Ω v ≤ α(max Ω u), whence α(min Ω u) ≤ v(x) ≤ α(max Ω u), x ∈ Ω. (7) From the equation f(u,v) = 0 we express v = β(u) :≡ u (1 + κu2)µu , for u > 0. Observe that for a fixed u > 0 f(u,v) < 0 iff v > β(u) and f(u,v) > 0 iff v < β(u). Choose xmin,xmax ∈ Ω such that u(xmin) = min Ω u, u(xmax) = max Ω u. Because ∆u(xmin) ≥ 0, ∆u(xmax) ≤ 0 the equal- ity �2∆u + f(u,v) = 0 implies f(u(xmin),v(xmin)) ≤ 0, f(u(xmax),v(xmax)) ≥ 0. Therefore, v(xmin) > β(u(xmin)), v(xmax) < β(u(xmax)). We obtain the following estimates max Ω v(x) > β(min Ω u), min Ω v(x) < β(max Ω u). (8) Combining estimates (7) and (8) we have β(min Ω u) < α(max Ω u), β(max Ω u) > α(min Ω u). Note that α is monotone increasing in u, so the estimates in Lemma 1 transforms these inequalities to β(min Ω u) < α(e √ µu/� min Ω u), β(max Ω u) > α(e− √ µu/� max Ω u). Therefore, min u ≥ ζ1, where ζ1 is the solution of β(z) = α(e √ µu/�z), (9) and max u ≤ ζ2 , where ζ2 is the solution of β(z) = α(e− √ µu/�z). (10) Hence we obtain the functions ψ↑,ψ↓ by setting ψ↑(�) = ζ1(�), ψ↓(�) = ζ2(�), and applying Lemma 2 completes the proof. Remark: The functions ψ↑,ψ↓ may have jump discontinuities. To estimate the behaviour of ζi, i = 1, 2, we need Lemma 2. The equation (9) has as unique solution ζ1 = 0 if µv < µ2u. Otherwise, the equation (9) has a unique non-negative solution ζ1(�) which is an increasing function of �. Furthermore, as � → 0, ζ1(�) → 0, and as � →∞, ζ1(�) → 1κ( µv µ2u − 1). The equation (10) has a unique non-negative solution ζ2(�) which is a decreasing function of �. Furthermore, as � → 0, ζ2(�) → ∞, and as � → ∞, ζ2(�) → 0 if µv < µ2u and ζ2(�) → 1κ( µv µ2u − 1) else. Proof: Note that α(0) = β(0) = 0, and α′ > 0, while β has a maximum at z = 1√ k . Let s = e± √ µu/�,C = µv µ2u . The problem reduces to solving the equation β(z) = α(sz) or C (1 + κz2)2 = s2 1 + κs2z2 . This is a quadratic in z2, s2 −C + (2 −C)κs2z2 + κ2s2z4 = 0. (11) This quadratic has real solutions if its discriminant D = κ2s2(C2s2 + 4C − 4Cs2) ≥ 0. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 4 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... First, note that s = e √ µu/� ≥ 1, so by Descartes’ rule of signs (11) has no positive so- lutions for C < 1. For 1 ≤ C < 2 (11) has a positive solution iff C > s2. The positive solution of (9) is the positive square root ζ1 = √ Z where Z(s) = C − 2 2κ + √ C2s2 + 4C − 4Cs2 2κs . For C ≥ 2 (11) has at least one positive solution iff D > 0. Observe that as � → ∞,s → 1 so the positive solution ζ1 tends to the square root of lim s→1 Z(s) = C − 1 κ . Second, note that s = e− √ µu/� ≤ 1, so D > 0 (solutions are always real). For C < 1, Descartes’ rule of signs shows (11) has a positive solution iff s2 < C. For C ≥ 1 Descartes’ rule of signs shows that (11) has always one non-negative solution. Hence we conclude that for µv µ2u < 1, ζ1 = 0, and for µv µ2u ≥ 1, the solution ζ1 of (9) is an increasing function of �, and lim �→0 ζ1(�) = 0, lim �→∞ ζ1(�) = √ 1 κ ( µv µ2u − 1). Furthermore, ζ2(�) →∞ as � → 0, while lim �→∞ ζ2(�) = 0 for µv µ2u < 1, but lim �→∞ ζ2(�) = √ 1 κ ( µv µ2u − 1) for µv µ2u ≥ 1. In the singular-perturbation limit � → 0, the problem (3) will exhibit spike solutions with the spikes in u having small support in Ω. We refer to [2], [15] for construction and analytic properties of such solutions. Fig. 1 shows a typical spike solution. Fig. 1. A pattern with spikes. Parameter values are � = 0.01,d = 0.1,µu = 1,µv = 1.2. C. Reduced problem We show that the reduced problem (� = 0) ad- mits another class of solutions. By setting � = 0 in (3) the resulting reduced problem is an algebraic- PDE system 0 = f(U,V ), 0 = d∆V + g(U,V ), 0 = ∂xV, x ∈ ∂Ω. (12) In the following, we shall characterise solutions of (12) and their stability properties. Let us recall some basic definitions. A solution (U,V ) to the problem (12) is called regular if there exists a function h ∈ C1(R) such that the a solution U(x) of f(U,V ) = 0 given by U(x) = h(V (x)) for all x ∈ Ω. If the function h is not unique, there is more than one way to choose the solution for U in the first equation in (12), so the problem may have only piecewise continuous solutions U on Ω. Hence it is convenient to study such solutions in a weak sense. The weak solution (U,V ) of (12) belongs to the class L∞(Ω)×H1(Ω) and satisfies 0 = f(U,V ), a.e. x ∈ Ω, d〈∇V,∇ψ〉 = 〈g(U,V ),ψ〉, ψ ∈ H1(Ω), (13) where 〈·, ·〉 denotes the H1-scalar product. Proposition 1. Suppose the problem (3) has a spa- tially homogeneous steady state (û, v̂). Diffusion- Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 5 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... driven (Turing) instability at (û, v̂) occurs if u is self-activating there, in other words, fu(û, v̂) > 0. Proof: The dispersion relation is a quadratic in λ. Solving the dispersion relation M(λ, 0) = 0 for λ we obtain the conditions trJ = fu + gv < 0, (14) det J = fugv −fvgu > 0, (15) so that Reλ1,2 < 0. The derivatives are evaluated at (û, v̂). Next we solve the dispersion relation M(λ,k) = 0 for λ. By Vieta’s formulae λ1 + λ2 = −dk2 + trJ, λ1λ2 = −fudk2 + det J, whence λ1 + λ2 < 0, for all k > 0. If fu ≤ 0, λ1λ2 > 0, so Reλ1,2 < 0 for all k > 0, and no diffusion-driven instability would be possible. This proves the claim. III. AUXILIARY PROBLEM In this section, we consider an auxiliary elliptic problem when at least one of U,V > 0 on some subinterval of Ω. Suppose that the equation f(U,V ) = 0 can be solved (not necessarily) uniquely on a subset I ⊂ Ω. Let U(x) = h(V (x)),x ∈ I, with h ∈ C1(R). Then every regular solution of (12) on I satisfies the elliptic problem 0 = d∆V + φ(V ), x ∈ I, (16) with φ(V ) = g(h(V ),V ). The solutions of (16) can be constructed using an energy method for two-point boundary value problems. There are two cases for the function φ, depending on whether U = 0 on I or U > 0 on I. Let us consider each case separately. If U(x) = 0,x ∈ I, we have h ≡ 0, so φ = −µvV almost everywhere on I. The problem (16) is reduced to the elliptic problem (17), 0 = d∆V −µvV, x ∈ I, (17) which has a non-trivial solution only under Dirich- let or Robin boundary conditions. Next we classify the solutions when h 6≡ 0 on I. We solve formally for U in (12), U = hi(V ) = 1 ± √ 1 − 4κµ2uV 2 2κµuV , i = 1, 2. (18) and use U2 V (1 + κU2) = µuU, x ∈ I. When U > 0, the equation f(U,V ) = 0 may have locally at most two solutions. Then (16) becomes 0 = d∆V + µuhi(V ) −µvV, (19) so we must solve a two-point boundary-value problem in two cases i = 1, 2 (for each solution branch for U), 0 = ∆V + 1 d (µuhi(V ) −µvV ), x ∈ I. (20) We set φi(y) = 1 d (µuhi(y) −µvy) = 1 d ( 1 ± √ 1 − 4κµ2uy2 2κy −µvy ) (21) with φ1 denoting the choice of positive square root and φ2 denoting the choice of negative square root in (21). The auxiliary elliptic problem is thus formu- lated: Solve for V = y(x) such that 0 = y′′ + φi(y), x ∈ I, i = 1, 2. (22) Recall that only solutions y ∈ (0, 1 2µu √ k ) are considered in order for the square root in (21) to be real-valued. Problem (22) can be rewritten as the equivalent system (23) of first-order equations y′ = z, z′ = −φi(y). (23) Note that y < 1 2µu √ κ in order for the square root in (21) to be well-defined in R. Hence, without loss of generality we may assume that V < 1 2µu √ κ on I. Else, we restrict the domain to a subset {x : V (x) < 1 2µu √ κ }⊂ I. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 6 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... We employ an energy formulation to describe the solutions y,z of (23). We set E as the total energy, U as the potential energy. The first integral of (23) is z2 2 + U(y) = E, U′ = φi. (24) To see this, differentiate the left-hand side of (24) and apply (22) and (23),( z2 2 + U(y) )′ = zz′ + U′y′ = zz′ + φiz = 0. Let U have a local minimum at y0. Choose total energy E such that U(0) > E > U(y0). According to [1, Satz, p.92], for E > U(y0), the equation (24) defines a closed smooth curve in the (y,z)- plane, which is symmetric with respect to the y- axis. Then there exists some t > 0 such that z(0) = z(t), corresponding to a solution of prob- lem (22) under homogeneous Neumann boundary conditions on (0, t). As in [1] we express the solution y of the ODE using (24) as y′ = ± √ 2(E −U(y)). Choosing the positive value of y′ (corresponding to a monotone increasing y), rearranging the above as 1 = y′√ 2(E −U(y)) (25) and integrating both sides of (25) over x, we obtain for every L > 0, L = ∫ L 0 y′(x)√ 2(E −U(y(x)) dx = ∫ y(L) y(0) dy√ 2(E −U(y)) . (26) Let 0 < y1 < y0 < y2 < 12µu √ κ be such that U(y1) = U(y2) = E, but U′(yi) 6= 0, i = 1, 2. Then the integral I(E) = ∫ y2 y1 dy√ 2(E −U(y)) := L 2 . (27) is convergent [1, p.93]. Then (26) defines im- plicitly a continuous solution y of (22) such that y(0) = y1,y( L 2 ) = y2. This solution can be con- tinued periodically on R and the periodic function has a period L 2 . Therefore, every such closed curve for suitable L will correspond to a solution of the system (23) under homogeneous Neumann boundary condi- tions on (0, L 2 ). The properties of the solutions of (23) will depend, therefore, on the properties of the integral (27), I(E). The following lemma is a modification of a well-known result. For the idea of proof we refer to [11, Lemma 3.1] or [5, Lemma 5.3-5.5]. Lemma 3. Let U have a local minimum at x0 and a local maximum at 0. Suppose 0 < x1 < x0 < x2 are such that U(x1) = U(x2) = E,U′(x1),U′(x2) 6= 0. Then I(E) is a continuous function in E, and lim E→U(x0) I(E) = πU′′(x0) , lim E→U(0) I(E) = ∞. The extrema of the potential energy Ui depend on the zeros of the functions φi. Let us examine φi’s zeros for i = 1, 2. Note that in order for the square root in (21) to be real-valued, y > 0 is such that 1 − 4κµ2uy2 ≥ 0, so we search for zeros in the interval (0, 1 2µu √ k ). Lemma 4. The equation φ1(y) = 0 has no solutions in (0, 1 2µu √ k ). For µv < µ2u or µv > 2µ 2 u, φ2(y) = 0 has no solution in (0, 1 2µu √ k ). For µ2u < µv < 2µ 2 u, the solution of φ2(y) = 0 is y0 = √ µv−µ2u κµ2v . Proof: After rearrangement of the terms and squaring both sides, we obtain 1 − 4κµ2uy2 = (1 − 2κµvy2)2, and after cancellation of y2 from both sides, κµ2vy 2 + µ2u − µv = 0. For µv < µ2u, the left- hand side is strictly positive, hence neither φi has a positive root. If µv > µ2u, a direct computation shows that the solution of the above quadratic is y0 = √ µv−µ2u κµ2v . Yet, y0 is a root of φ2 only. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 7 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... When µv < 2µ2u, 1 − 2κµvy2 > 1 − 4κµ2uy2 ≥ 0, so φ1(y) = 1 d ( 1 − 2κµvy2 + √ 1 − 4κµ2uy2 2κy ) > 0. Therefore, for µ2u < µv < 2µ 2 u, φ1(y) = 0 has no solution y ∈ R+. For µv > 2µ2u, the number y0 lies outside the domain of definition of the square root, (0, 1 2µu √ k ). Hence none of φi, i = 1, 2 has a zero in (0, 1 2µu √ k ). Now we are able to characterise the extrema of the potential energies U1,U2 on (0, 1 2µu √ k ). Lemma 5. Let µu,µv > 0. The potential energy U1 has a maximum at 1 2µu √ k and no local minima. • for µv < µ2u, U2 has a minimum at 0; • for µv ∈ (µ2u, 2µ2u), U2 has a local maximum at 0 and a local minimum at y0 = √ µv−µ2u κµ2v ; • for µv > 2µ2u, U2 has a maximum at 0. Proof: Lemma 4 implies that φ1 never changes sign on the interval (0, 1 2µu √ k ). Thus, it is clear that U1 has no local extrema in (0, 1 2µu √ k ). In fact, lim y→0 U1(y) = −∞, and U1 has a maximum at 1 2µu √ k . Furthermore, U2 has an extremum at y0 =√ µv−µ2u κµ2v , see Lemma 4. Note that for this y0,√ 1 − 4κµ2uy20 = ∣∣∣∣1 − 2µ2uµv ∣∣∣∣ = 2µ2uµv − 1. Next, we compute U′′2 (y0) = φ′2(y0). Note that φ′2(y) = 1 d ( 1 2κy2 √ 1 − 4κµ2uy2 − 1 2κy2 −µv ) . If µv ∈ (µ2u, 2µ2u), U′′2 (y0) = φ′2(y0) = 2µv(µv −µ2u) d(2µ2u −µv) > 0. Hence, U2 has a local minimum at y0. If µv = 2µ2u, the point y0 = 1 2µu √ κ coincides with the endpoint of the interval, so it is of no interest. We compute by L’Hôpital’s rule lim y→0 φ2(y) = 0, showing 0 is an extremum for U2. Next we ex- amine the extremum properties of 0 by applying again L’Hôpital’s rule lim y→0 φ′2(y) = 1 d (µ2u −µv). Therefore, 0 is a maximum for U2 if µ2u < µv, and a minimum if µ2u > µv. This completes the proof. Lemma 5 and [1, Satz, p.92] allow us to relate the existence of regular solutions of problem (12) on Ω to the extrema of the potential energies asso- ciated to the auxiliary problem (22). We conclude that no regular solutions (U,V ) can be constructed using the potential energy U1 because it does not have local minima. The only possibility to construct regular solutions (U,V ) is by using the potential energy U2 when µ2u < µv < 2µ2u. The following Lemma provides an important property of the integral I(E) which will be em- ployed in the construction of regular solutions of (12). Lemma 6. I(E) is monotone in E on (U2(y0),U2(0)). Proof: Using reasoning as in [5, Lemma 5.5] it is enough to show that U′′′2 ≤ 0 on (0, 12µu√κ). Then we estimate U′′′2 (y) = φ′′2 (y) = 2 y3 (1 − (1 − 4µ2uκy2)− 1 2 ) − 4κµ 2 u y (1 − 4κµ2uy2)− 3 2 , but (1 − 4µ2uκy2)− 1 2 ≥ 1 so U′′′2 (y) < 0. IV. STATIONARY SOLUTIONS Combining the results on the auxiliary problem, we can classify the regular and the weak solutions of the problem (22). Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 8 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... A. Regular stationary solutions We first remark on the ‘trivial’, constant solution to (12). If µ2u < µv, there is a constant solution of (12), û = √ µv −µ2u κµ2u , v̂ = √ µv −µ2u κµ2v . (28) For µ2u < µv < 2µ 2 u we can construct a regular solution using the auxiliary problem and the potential energy U2 as follows. Proposition 2. Let Lmin := min E I(E) = π√ U′′2 (√ µv−µ2u κµ2v ) and set N = max{n ∈ N : NLmin ≤ 1}. The two-point boundary value problem (20) with homogeneous Neumann boundary conditions has the following solutions • the spatially homogeneous solution V = v̂ given in (28). • if Lmin ≤ 1, there exists a unique mono- tone increasing solution V↑(x), and a unique monotone decreasing solution V↓(x) = V↑(1 −x) • for all 2 ≤ n ≤ N there exists a unique n-periodic solution Vn,↑ which is monotone increasing on (0, 1 n ), as well as a unique n-periodic solution Vn,↓ which is monotone decreasing on (0, 1 n ). Proof: The validity of the first claim is ob- vious. The remaining claims use the properties of the integral I in Lemma 3 and 6. These imply that for all n ≤ N, there exists a unique energy level En : I(En) = 1n, corresponding to a unique mono- tone increasing solution Vn,↑(x) on (0, 1 n ), and a monotone decreasing solution Vn,↓ = Vn,↑( 1 n −x) on (0, 1 n ). If n ≥ 2, either solution can be extended to the entire domain Ω by the folding principle [1], [11]. B. Weak stationary solutions The previous section showed that for values of µu,µv such that µv < µ2u (12) has only weak solutions. These solutions are constructed piecewise using the auxiliary problem for each branch of U = hi(V ). Then U ∈ L∞(Ω) and V is continuous on Ω. Start on the y-axis at x = 0 and begin tracing along any admissible trajectories defined by • (y,z) : z = ± √ 2E −U0(y), • (y,z) : z = ± √ 2E −U1(y), or • (y,z) : z = ± √ 2E −U2(y). Here U0(y) = −µvd y 2 is the potential energy associated with the problem (17). The potential energies U1,2 are represented in closed form (up to a constant) by U1(y) = 1 dκ log y + 1 2dκ √ 1 − 4κµ2uy2 − µv 2d y2 − 1 2κ log ( 1 2µu √ k + √ 1 4µ2uκ −y2 ) , U2(y) = − 1 2dκ √ 1 − 4κµ2uy2 − µv 2d y2 + 1 2dκ log ( 1 2µu √ k + √ 1 4µ2uκ −y2 ) . Continue tracing until returning to the y-axis at x = 1 (Fig. 4). In this way we obtain a partitioning of the domain Ω into subintervals Ii. On each Ii the solution V is given by the y-coordinate of the admissible solution trajectories in the (y,z)-space. Of course, under this construction U may be discontinuous in Ω as on each subinterval Ii U is given by a different branch hi. However, note that V := y belongs to C1(Ω) because by construction z = y′ is continuous at the intersection of such trajectories. On Fig. 2 and Fig. 3 are plotted several solution curves corresponding to an energy level E for the different cases of potential energies Ui. Note, for example, that a weak solution can be traced by following a trajectory given by U0,U2 (in that order) when µv < µ2u or U2,U1 (in that order) when µv > 2µ2u. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 9 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... U1 E z y y z y y U0 E Fig. 2. Solution curves in (y,z)-space associated to potential energies U1 and U0. Fig. 3. Solution curves in (y,z)-space associated to potential energies U2 for the different cases (from left to right) µv > 2µ2u,µ 2 u < µv < 2µ 2 u,µv < µ 2 u. V. STABILITY OF STATIONARY SOLUTIONS For a stationary solution (U,V ) we can establish linear stability in the classical sense: for small perturbations ũ(t,x) = u(t,x) − U(x), ṽ(t,x) = v(t,x)−V (x), the behaviour of the solutions u,v of (3) is governed locally by the linear approxi- z yV (0) V (1) Fig. 4. A weak solution on Ω following different energy trajectories. mation. By setting an Ansatz ũ(t,x) = u(x)eλt, ṽ(t,x) = v(x)eλt, for local stability we have to consider the sign of Reλ. λ is an eigenvalue of the linearised differen- tial operator L = diag(0,d∆) + J, where J = ( 2U V (1+κU2)2 −µu − U 2 V 2(1+κU2) 2U V (1+κU2)2 − U2 V 2(1+κU2) −µv ) . If every eigenvalue λ of L has negative real part, the stationary solution (U,V ) is locally stable. Note that the spectrum of L need not be discrete. In the one-dimensional case we formulate re- sults on the stability of spatially nonuniform sta- tionary solutions of the system. Easy linear stabil- ity analysis leads to Proposition 3. Let µ2u < µv < 2µ 2 u, the constant solution (û, v̂) is locally asymptotically unstable. When µv > 2µ2u, the constant solution (û, v̂) is locally asymptotically stable. Proof: A linearisation of the right-hand side of (12) at (û, v̂) gives the following: L = diag(0,d∆) + ( 2µ3u µv −µu −µv 2µ3u µv −2µv ) . Note that the saturation parameter κ does not influence the stability of the constant solution. L has constant coefficients and its spectrum can be computed by matrix eigenvalue analysis. For values µv > 2µ2u, the spectrum of L lies entirely in the left half-plane. Hence, the constant solution (û, v̂) is locally asymptotically stable. For values µv < 2µ2u, L has positive eigenval- ues. This proves the claim. The following result establishes the local insta- bility of spatially heterogeneous regular solutions. Theorem 2. Let µ2u < µv < 2µ 2 u. Any spatially heterogeneous regular solution (U,V ) of (12) on the interval Ω is unstable. Proof: For the proof we use a result on the spectrum of L established in [6, Corollary 2.7], Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 10 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... namely the effect of autocatalysis in the nonlinear- ity f which implies that the spectrum of L com- prises eigenvalues with positive real part. This fact implies linear instability of the stationary solution (U,V ). Note that f(0,V ) = 0 for all V ∈ R. Let (U,V ) be a regular solution of (12). Then due to the energy method construction, U = h2(V ) > 0 is continuous on I, and so is fu(U,V ) = 2U V (1 + κU2)2 −µu. Furthermore, using the formula (18) we see that U = h2(V ) < κ −1/2. We estimate, first, using the relation f(U,V ) = 0 on I, fu(U,V ) = 2µu 1 + κU2 −µu < µu Second, fu(U,V ) = 2µu 1 + κU2 −µu > 2µu 1 + 1 −µu = 0. The result of [6, Corollary 2.7] implies that the linearised differential operator L has eigenvalues with positive real part. Hence, (U,V ) is an unsta- ble solution. The result of Theorem 2 is also a consequence of the results in [7] that establish instability of heterogeneous solutions for semilinear diffusion equations. However, in this particular case insta- bility of the regular solution can be established by direct computation. Theorem 3. Let µv < 2µ2u. Any weak solution (U,V ) of (12) such that U < κ−1/2 on Ω is unstable. The proof is identical to that of Theorem 2. Theorem 4. Let n ∈ N, 0 ≤ x1 < x2 < ... < xn ≤ 1. The pair (U,V ) defined by U(x) > 0,x = xj,U(x) = 0,x 6= xj, and V a solution to (17) with Robin boundary conditions, is a solution to (13). Moreover, any such (U,V ) is locally asymptotically stable. Proof: For any (U,V ) of the given type, the linearised operator L looks almost everywhere like L = diag(0,d∆) + ( −µu 0 0 −µv ) . Therefore, the spectrum of L is bounded away from 0 in the left halfplane almost everywhere. Hence, any stationary solution (U,V ) is locally asymptotically stable. VI. DISCUSSION Theorem 3 and Theorem 4 imply that the locally stable weak solutions must necessarily exhibit large jump discontinuities of amplitude at least κ−1/2. In the literature such solutions are said to exhibit striking patchiness [8], [9] or transition layers [13]. The results show that all spatially heteroge- neous, strictly positive, regular solutions of (12) over a one-dimensional domain Ω are unstable. This is in contrast to the singularly perturbed sys- tem (3) where stable spike solutions exist [15]. The asymptotically stable solutions, which are different from the homogeneous steady state (û, v̂) (28) are discontinuous solutions that exhibit large transition layers. Fig. 5 shows such a pattern. In contrast to the results in [5], [6], the prob- lem (12) does not fulfil the autocatalysis condition in [6, Corollary 2.7], whose estimates do not hold for all bounded weak solutions of (12). That explains the existence of stable weak solutions with large transition layers. Fig. 5. A discontinuous pattern with large amplitude. Parameter values are µu = 1,µv = 1.2,κ = 0.1,d = 0.1. Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 11 of 12 http://dx.doi.org/10.11145/j.biomath.2014.11.111 P. Rashkov, Regular and Discontinuous Solutions in a Reaction ... ACKNOWLEDGMENT The author acknowledges the support of the Centre for Synthetic Microbiology in Marburg, promoted by the LOEWE Excellence Program of the state of Hesse, Germany. Current address: College of Life and Environmental Sciences, University of Exeter, Stocker Road, Exeter EX4 4QD, United Kingdom. REFERENCES [1] V. I. Arnold, Gewöhnliche Differentialgleichungen. VEB Deutscher Verlag der Wissenschaften, Berlin, 1979. [2] A. Doelman, R. A. Gardner, and T. J. Kaper. Large stable pulse solutions in reaction-diffusion equations. Indiana Univ. Math. J. 50 (1) 443–507, 2001. [3] A. Gierer and H. Meinhardt. A theory of biological pattern formation, Kybernetik 12, 30–39, 1972. [4] A. J. Koch and H. Meinhardt. Biological pattern forma- tion: from basic mechanisms to complex structures. Rev. Mod Phys. 66 (4), 1481–1507, 1994. [5] A. Marciniak-Czochra, G. Karch, and K. Suzuki. Un- stable patterns in reaction-diffusion model of early car- cinogenesis. J. Math. Pures Appl. 99, 509–543, 2013. DOI:10.1016/j.matpur.2012.09.011 [6] A. Marciniak-Czochra, G. Karch, and K. Suzuki. Unsta- ble patterns in autocatalyctic reaction-diffusion systems. Preprint. [7] H. Matano. Asymptotic behavior and stability of solutions of semilinear diffusion equations. Publ. RIMS, Kyoto Univ. 15, 401-454, 1979. [8] M. Mimura and J.D. Murray. On a diffusive prey-predator model which exhibits patchiness. J. Theor. Biol. 75, 249– 262, 1979. [9] M. Mimura, M. Tabata, and Y. Hosono. Multiple solu- tions of two-point boundary value problems of Neumann type with a small parameter. SIAM J. Math. Anal. 11 (4), 613–631, 1980. [10] J. D. Murray, Mathematical Biology. Springer, New York, 1993. [11] Y. Nishiura. Global structure of bifurcating solutions of some reaction-diffusion systems. SIAM J. Math. Anal. 13 (4), 555–593, 1982. [12] P. Rashkov. Remarks on pattern formation in a model for hair follicle spacing. preprint. [13] K. Sakamoto. Construction and stability analysis of transition layer solutions in reaction-diffusion systems. Tôhoku Math. J. 42, 17–44, 1990. [14] S. Sick, S. Reinker, J. Timmer, and T. Schlake. WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism. Science 314 (5804), 1447– 1450, 2006. DOI:10.1126/science.1130088 [15] F. Veerman and A. Doelman. Pulses in a Gierer- Meinhardt equation with a slow nonlinearity. SIAM J. Dyn. Sys. 12(1), 28–60, 2013. DOI:10.1137/120878574 Biomath 3 (2014), 1411111, http://dx.doi.org/10.11145/j.biomath.2014.11.111 Page 12 of 12 http://dx.doi.org/10.1016/j.matpur.2012.09.011 http://arxiv.org/abs/1301.2002 http://dx.doi.org/10.1126/science.1130088 http://dx.doi.org/10.1137/120878574 http://dx.doi.org/10.11145/j.biomath.2014.11.111 Introduction The model problem Diffusion-driven instability A priori estimates Reduced problem Auxiliary problem Stationary solutions Regular stationary solutions Weak stationary solutions Stability of stationary solutions Discussion References