Spatiotemporal patterns and multiple bifurcations of a reaction- diffusion model for hair follicle spacing

1.   Introduction

The development of regularly arranged body parts has attracted the attention of a huge number of experimental biologists and theoreticians alike, and one of the crucial issues is getting to know the underlying mechanism in the formation of epidermal appendages such as feathers and hairs. Many theoretical models, including mathematical models of coupled reaction-diffusion equations, have been used to describe the formation of animal pigmentation patterns and distribution.

In ^[1], Sick et al. proposed a system of reaction-diffusion equations to model the influence of the Wnt signaling pathway in primary hair follicle initiation in mice. It is suggested that Wnt and Dkk have a primary influence on the spacing patterns of hair follicles in mice. The interactions between Wnt and Dkk are modeled by using a modified Gierer-Meinhardt reaction-diffusion (activator-inhibitor) model. The Wnt-Dkk interaction is modeled by the following 2-component reaction-diffusion equations:

where $\Omega$ is an open bounded domain in ${\mathbf R}^N$, $N\geq1$, with the smooth boundary $\partial\Omega$; $u = u(x, t)$ and $v = v(x, t)$ stand for the concentrations of the activator Wnt and the inhibitor Dkk at time $t$ and position $x\in\Omega$, respectively; $d_1 > 0$ and $d_2 > 0$ are the diffusion rates of the activator Wnt and the inhibitor Dkk, respectively; $\gamma$ and $k$ are non-negative saturation parameters. The constants $\rho_1$ and $\rho_2$ scale the production rates of the activator Wnt and the inhibitor Dkk, respectively; The negative terms denote that both chemicals decay linearly with constants $\mu$ and $\lambda$; $u_0, v_0\in C^2(\Omega)\cap C^0(\overline{\Omega})$ and the Neumann boundary conditions indicate that there is no flux of the chemical substances of $u$ and $v$ on the boundary.

System (1.1) is a modification of the homogeneous Gierer-Meinhardt model of the following form:

where $u = u(x, t)$ and $v = v(x, t)$ stand for the concentrations of the activator and the inhibitor at time $t$ and position $x\in\Omega$, respectively. The non-negative functions $\sigma_1(x)$ and $\sigma_2(x)$ are background source terms for the activator and the inhibitor, respectively. The exponents $p, q, r$ and $s$ are non-negative.

In ^[2], global existence of the solutions of System (1.1) is considered by calculating the uniform bounds. Analysis of the conditions for the emergence of spatially heterogeneous solutions is performed by using a limiting form of the original reaction-diffusion system. The conditions for pattern formation given in ^[1] are improved by including those subregions in the parameter space where far-from-equilibrium heterogeneous solutions occur.

To simplify the analysis of the Turing mechanism, in ^[1], a modification of System (1.1) is considered. In ^[1], the parameter $\gamma$ is chosen so that $\gamma\approx0$. Then, they obtained the following reaction-diffusion equations:

As argued in ^[2], System (1.3) captures the dynamics of (1.1) quite well, and, due to the smaller number of parameters, the conditions for emergence of Turing patterns are easier to analyze.

In ^[3], Rashkov considered the Turing instability of the positive constant equilibrium solution, as well as the existence and stability of both the regular and the discontinuous non-constant stationary solutions of System (1.3). Veerman and Doelman ^[4] showed that all of the positive non-constant regular equilibrium solutions (in the one-dimensional spatial domain) are unstable. The results are quite different from the corresponding singularly perturbed system, which has stable non-constant spike solutions.

We would also like to mention that a similar but different model used to describe the hair follicle bulb was proposed by Mooney and Nagorcka ^[5,6,7]. In particular, Nagorcka and Mooney ^[7] described a theoretical mechanism for cell differentiation based on the substances $X$ and $Y$, which constitute a reaction-diffusion system, and the substance $Z$, which diffuses radially outward from the dermal papilla through the bulb. The model describes the reaction and diffusion of morphogens $X$ and $Y$, and it is defined by

where $u$ and $v$ can be seen as reduced or dimensionless concentrations of morphogens of $X$ and $Y$, respectively, with positive parameters $a, b, c, e$ and $r$ and diffusion coefficients $d_{1}$ and $d_{2}$. The rates $a, b, c$ and $e$ have units of $time^{-1}$, and $r$ and $p(>1)$ are dimensionless quantities. The terms $-au + bv^p/(1+v^p)$ and $-cu+e(v^p+r)/(1+v^p)$ represent the net rates of the productions $u$ and $v$, respectively, where the reaction-diffusion system constitutes an activator-inhibitor pair with the morphogens X the inhibitor and Y the activator. For System (1.4), Yi et al. ^[8] provided some global analyses of the model that were dependent upon some parametric thresholds/constraints. They found that, when one of the dimensionless parameters is less than one, the unique positive equilibrium is globally asymptotically stable. On the contrary, when this threshold is greater than one, the existence of both steady-state and Hopf bifurcations can be observed under further parametric constraints. In ^[9], Yang and Ju considered the Turing instability of the spatially homogeneous periodic solutions of System (1.4).

In this paper, we are mainly concerned with the dynamics of both the reaction-diffusion equations (1.3) and its corresponding ODE system. The highlights of the paper are as follows:

1. Consider the dynamics of the ODEs. By using the decay rate $\lambda$ of the inhibitor $v$ as the main bifurcation parameter, we are able to show the existence of Hopf bifurcating periodic solutions. That is, when the decay rate $\lambda$ crosses the critical value $\lambda_0$, the densities of the activator $u$ and the inhibitor $v$ will undergo temporal oscillations. Moreover, we are able to prove that, once the Hopf bifurcating periodic solutions occur, they must be unstable. This is one of the novel points of this paper.

2. Consider the dynamics of the PDEs. Three kinds of results are obtained. First, under suitable conditions on the decay rate $\lambda$ of the inhibitor $v$ and the diffusion rates ($d_1$ and $d_2$), $(u_ \lambda, v_ \lambda)$ will always be stable in the reaction-diffusion equations. In this case, the dynamics of the PDEs can be accurately determined by the dynamics of the ODEs. This is the so-called lumped parameter phenomenon ^[10]. Second, we are able to derive precise conditions on $\lambda$, as well as the diffusion rates ($d_1$ and $d_2$) such that Turing instability (see ^[11]) of $(u_ \lambda, v_ \lambda)$ occurs. In this case, Turing patterns (spotted or striped patterns) of the reaction-diffusion system can be observed. Finally, under suitable conditions, non-constant Hopf bifurcating periodic solutions (both spatially homogeneous and spatially non-homogeneous) and non-constant, positive, steady-state bifurcating solutions can be obtained. These results are new and cannot be found in the existing literature. This is the other novel point of this paper.

The structure of this paper is as follows. In Section 2, we mainly study the dynamic behavior of the corresponding ordinary differential equation system of System (1.3), including the existence, stability and instability of the positive steady-state solution and the existence, stability of the Hopf bifurcating periodic solution. In Section 3, we mainly study the dynamic behavior of the reaction-diffusion system presented as System (1.3), including the existence, stability and instability of the positive steady-state solution and the existence of Hopf bifurcating periodic solutions and non-constant, positive, steady-state bifurcating periodic solutions. In the Appendix section, we list some of our proofs conducted in Section 3.

2.   Dynamics of the ODE system

The corresponding ODE system of System (1.3) takes the following form:

System (2.1) has a unique positive equilibrium solution $(u_ \lambda, v_ \lambda)$, where

where it is assumed that $\lambda > \rho_2\mu^2/\rho_1^2$.

In what follows, we shall always assume that $\lambda > \rho_2\mu^2/\rho_1^2$ so that $(u_ \lambda, v_ \lambda)$ is the positive equilibrium solution of (2.1). We shall fix $\rho_1, \rho_2$ and $k$ and choose $\lambda$ and $\mu$ as two main bifurcation parameters.

We have the following results on the stability and instability of $(u_ \lambda, v_ \lambda)$:

Theorem 2.1. Let $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$ so that $(u_ \lambda, v_ \lambda)$ is the unique positive equilibrium solution of (2.1). Define

Then, $0 < \lambda_0 < \dfrac{2\rho_2\mu^2}{\rho_1^2}$. In particular, the following results hold true:

1. Suppose that $\mu\geq\dfrac{\rho_1^2}{2\rho_2}$ or, equivalently that $\lambda_0\leq \dfrac{\rho_2\mu^2}{\rho_1^2}$ holds. Then, for all $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, +\infty\bigg)$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to System (2.1);

2. Suppose that $\mu < \dfrac{\rho_1^2}{2\rho_2}$ or, equivalently that $0 < \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_0\bigg(<\dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$ holds. Then, the following conclusions hold:

(a) for any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_0\bigg)$, $(u_ \lambda, v_ \lambda)$ is unstable;

(b) for any $\lambda\in\big(\lambda_0, +\infty\big)$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to System (2.1).

(c) At $\lambda = \lambda_0$, System (2.1) undergoes a Hopf bifurcation at $\lambda = \lambda_0$. That is, there exists an $s^* > 0$ such that, for $s\in(0, s^*)$, there exists $(\lambda(s), Z(s), u(\cdot, s), v(\cdot, s))$ so that $(u(\cdot, s), v(\cdot, s))$ is a periodic solution of (2.1) with the minimum period $Z(s)\rightarrow 2\pi/\sqrt{D(\lambda_0)}$ and $(\lambda(s), u(\cdot, s), v(\cdot, s))\rightarrow (\lambda_0, u_{ \lambda_0}, v_{ \lambda_0})$ as $s\rightarrow0$. Moreover, the bifurcation is subcritical (i.e., the bifurcating periodic solution is unstable) and forward.

Proof. It is trivial to check that $0 < \lambda_0 < \dfrac{2\rho_2\mu^2}{\rho_1^2}$ under the assumption of $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$.

The linearized operator of System (2.1) evaluated at $(u_ \lambda, v_ \lambda)$ is given by

where we use the fact that $\lambda = \rho_2\mu^2(1+ku_ \lambda^2)/\rho_1^2$.

The eigenvalue problem of $J(\lambda)$ is governed by $\gamma^2-T(\lambda)\gamma+D(\lambda) = 0$, where

Clearly, $D(\lambda) > 0$ holds for all $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$. In particular, $T(\lambda) < 0$ for all $\lambda > \lambda_0$, $T(\lambda) > 0$ for all $0 < \lambda < \lambda_0$ and $T(\lambda_0) = 0$. Since $T\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg) < 0$, we have that $\lambda_0 < \dfrac{2\rho_2\mu^2}{\rho_1^2}$.

Part 1. Clearly, $\lambda_0\leq \dfrac{\rho_2\mu^2}{\rho_1^2}$ is equivalent to $\mu\geq\dfrac{\rho_1^2}{2\rho_2}$. Since $\lambda_0\leq \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda$, we have that $T(\lambda) < 0$ for all $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$. Since $D(\lambda) > 0$ holds for all $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$, we complete the proof.

Part 2. Case (a). Since $\dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda < \lambda_0$, we have that $T(\lambda) > 0$. Thus, $(u_ \lambda, v_ \lambda)$ is unstable; Case (b) is similar to Part 1.

We now prove Case (c). At $\lambda = \lambda_0$, we have that $T(\lambda_0) = 0$ and $D(\lambda_0) > 0$. Then, the eigenvalue problem has a pair of purely imaginary eigenvalues $\gamma = \pm{i}\sqrt{D(\lambda_0)}$. Let $\beta(\lambda)\pm{i}\omega(\lambda)$ be the eigenvalues of the eigenvalue problem satisfying $\beta(\lambda_0) = 0$ and $\omega(\lambda_0) = \sqrt{D(\lambda_0)}$. Then, we have

which shows that the transversality condition holds. Thus, at $\lambda = \lambda_0$, System (2.1) undergoes a Hopf bifurcation at $\lambda = \lambda_0$. That is, there exists an $s^* > 0$ such that, for $s\in(0, s^*)$, there exists $(\lambda(s), Z(s), u(\cdot, s), v(\cdot, s))$ so that $(u(\cdot, s), v(\cdot, s))$ is a periodic solution of (2.1) with the minimum period $Z(s)\rightarrow2\pi/\sqrt{D(\lambda_0)}$ and $(\lambda(s), u(\cdot, s), v(\cdot, s))\rightarrow (\lambda_0, u_{ \lambda_0}, v_{ \lambda_0})$ as $s\rightarrow0$.

Next, we shall study the bifurcation direction and the stability of the bifurcating periodic solutions. Rewrite (2.1) in the following form:

where

where all of the partial derivatives are evaluated at $(\lambda, u_ \lambda, v_ \lambda)$.

For $\lambda$ close to $\lambda_0$, we define a real $2$-by-$2$ matrix

where

and, for $\lambda$ close to $\lambda_0$,

Clearly, the matrix $Y$ is well defined since $J_{12}(\lambda_0) = - \frac{\rho_1}{\rho_2} \lambda_0 < 0$ implies that, for $\lambda$ close to $\lambda_0$, $J_{12}(\lambda)\neq0$.

By using the linear transformation $(u, v)^T = Y(x, y)^T$, we can reduce (2.7) to the following system:

where

With the Taylor expansion of $F(\lambda, u, v)$, we have

where

where

With the Taylor expansion of $G(\lambda, x, y)$, we have

where

with

where

Indeed, we have

where

Following page 90 of ^[12], we define

where

where all quantities are to be evaluated at $(\lambda_0, u_{ \lambda_0}, v_{ \lambda_0})$.

Substituting (2.15) into (2.14), we have

In Part 2, we assume that $\mu < \dfrac{\rho_1^2}{2\rho_2}$. Thus, we have that $\mu < \dfrac{3\rho_1^2}{2\rho_2}$, which is equivalent to $\lambda_0 < \mu$. Thus, we have ${\rm Re(c_1(\lambda_0))} > 0$. Together with (2.6), and in consideration of ^[13], it follows that the bifurcation is subcritical (i.e., the bifurcating periodic solution is unstable) and the bifurcation direction is forward.

Remark 2.2. Biological meaning: Remember that $\mu$ and $\lambda$ are the decay rates of the activator $u$ and the inhibitor $v$, respectively. Choosing $\lambda$ as the main bifurcation parameter (by fixing $\mu$) means that we want to know how the dynamics of System (2.1) changes as the decay rate of $v$ changes. Indeed, one can also choose $\mu$ as the main bifurcation parameter by fixing $\lambda$. (We did not consider this case in the paper, but we expect that System (2.1) can also exhibit similar results.) Our results indicate that, as the decay rate of $v$ changes, System (2.1) might have periodic solutions. That is, the density functions of $u$ and $v$ may undergo temporal oscillations with the evolution of time.

3.   Dynamics of the diffusive system

For convenience, we copy System (1.3) here:

We shall study the dynamics of System (3.1). Without loss of generality, we assume that $\Omega = (0, \ell\pi)$, with $\ell > 0$.

3.1.   Stability, instability and Turing instability of $(u_ \lambda, v_ \lambda)$

In this subsection, we shall consider the stability and Turing instability of the unique positive constant equilibrium solution $(u_ \lambda, v_ \lambda)$ with respect to the reaction-diffusion system presented as System (3.1).

We have the following results on the stability and Turing instability of $(u_ \lambda, v_ \lambda)$.

Theorem 3.1. Let $\mu\geq\dfrac{\rho_1^2}{2\rho_2}$ or, equivalently, $0 < \lambda_0\leq \dfrac{\rho_2\mu^2}{\rho_1^2}$ hold. Then, the following conclusions hold true:

1.For any $\lambda\in\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}, +\infty\bigg)$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1); in this case, the Turing instability of $(u_ \lambda, v_ \lambda)$ cannot feasibly occur;

2. For any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$,

(a) if, in addition, $d_1\geq\dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1); in this case, the Turing instability of $(u_ \lambda, v_ \lambda)$ cannot feasibly occur;

(b) if, in addition, $d_1 < \dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then we define $n_{d_1}$, with $1\leq n_{d_1} < +\infty$, to be the largest positive integer such that, for any integer $n\in[1, n_{d_1}]$,

Define $\widehat d_2: = \min\limits_{1\leq n\leq n_{d_1}}d_2^{(n)}$, where

Then, for any $d_2 < \widehat d_2$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). If $d_2 > \widehat d_2$ holds, then $(u_ \lambda, v_ \lambda)$ undergoes Turing instability in the reaction-diffusion system given by System (3.1).

Proof. The proof is moved to the Appendix.

Similarly, we have the following results:

Theorem 3.2. Let $\mu < \dfrac{\rho_1^2}{2\rho_2}$ or, equivalently, $\dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_0 < \dfrac{2\rho_2\mu^2}{\rho_1^2}$ hold. Then, the following conclusions hold true:

1. For any $\lambda\in\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}, +\infty\bigg)$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1); in this case, the Turing instability of $(u_ \lambda, v_ \lambda)$ cannot feasibly occur;

2. For any $\lambda\in\bigg(\lambda_0, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$,

(a) if, in addition, $d_1\geq\dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1).

(b) if, in addition, $d_1 < \dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then, for any $d_2 < \widehat d_2$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). If $d_2 > \widehat d_2$ holds, then $(u_ \lambda, v_ \lambda)$ undergoes Turing instability in the reaction-diffusion system given by System (3.1).

3. For any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_0\bigg)$, $(u_ \lambda, v_ \lambda)$ is unstable with respect to the reaction-diffusion system given by System (3.1); in this case, the Turing instability of $(u_ \lambda, v_ \lambda)$ cannot feasibly occur.

Remark 3.3. If one of the following cases holds, i.e., (1) of Theorem 3.1, 2(a) of Theorem 3.1, (1) of Theorem 3.2 or 2(a) of Theorem 3.2, then $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the reaction-diffusion system given by System (3.1). In this case, the Turing instability of $(u_ \lambda, v_ \lambda)$ will never occur. This means that, when the decay rate of $v$ and the diffusion rates are chosen in certain ranges, Turing patterns are unlikely to occur; If one of the following cases holds, i.e., 2(b) of Theorem 3.1 or 2(b) of Theorem 3.2, then $(u_ \lambda, v_ \lambda)$ will undergo Turing instability in the reaction-diffusion system given by System (3.1). This means that, when the decay rate of $v$ and the diffusion rates are chosen in certain ranges, Turing patterns can be expected.

3.2.   Hopf bifurcation analysis: Spatially non-homogeneous periodic solutions

In this subsection, we shall consider the occurrence of spatially non-homogeneous periodic solutions bifurcating from Hopf bifurcations.

When $\lambda_0\leq\dfrac{\rho_2\mu^2}{\rho_1^2}$, we have that $T(\lambda)\leq0$ (and, hence, $T_n(\lambda) \leq 0, T_n(\lambda)$ is defined in (A.1)) for $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$. Thus, to expect the spatially non-homogeneous periodic solution bifurcating from Hopf bifurcations, we need to concentrate on the case when $T(\lambda)\geq0$, or, equivalently,

In this case, under (3.2), $(u_ \lambda, v_ \lambda)$ is unstable with respect to the reaction-diffusion system given by System (3.1).

Following ^[14], we shall identify the Hopf bifurcation values, denoted by $\lambda_H$, which satisfy the following: there exists $n\in \mathbb{N}_0$ such that

and, for the unique pair of complex eigenvalues near the imaginary axis $\alpha(\lambda)\pm i \omega(\lambda)$,

Since $\lambda_0 > \dfrac{\rho_2\mu^2}{\rho_1^2}$ and $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_0\bigg]$, we have that $T(\lambda)\geq0$. One can check that $T'(\lambda)$ is decreasing in $(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_0)$, given that $T(\lambda_0) = 0$ and $\lim\limits_{ \lambda\rightarrow0^+}T(\lambda) = +\infty$. Define

Then, for any $\ell > \ell_1$, there exists an $n\in\mathbb N_0\backslash\{0\}$ such that $\ell_n < \ell < \ell_{n+1}$. Then, for any $1\leq j\leq n$, $T(\lambda) = \dfrac{(d_1+d_2)j^2}{\ell^2}$ has a unique root, denoted by $\lambda_j^H$, satisfying

and if $n\neq m$, $\lambda_n^H\neq \lambda_m^H$.

It remains necessary to derive precise conditions so that, for all $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_0\bigg)$, $D_n(\lambda) > 0$ for all $n\in\mathbb N_0$. In fact, since

we can regard $D_n(\lambda)$ as the quadratic function of $n^2/\ell^2$. Define

Then, we have that $D_n(\lambda) = \mathcal F(n^2/\ell^2)$. The discriminant of $\mathcal F(x)$ is given by

By (3.7), to let $\Delta_{\mathcal F} < 0$, we need first to assume that $d_1 \lambda-d_2\mu < 0$ or, equivalently that $\dfrac{d_2}{d_1} > \dfrac{ \lambda}{\mu}$. Then, solving $\Delta_{\mathcal F} < 0$ from (3.7), we have

If $8-4\sqrt{3} < \dfrac{\rho_1^2 \lambda}{\rho_2\mu^2} < 8+4\sqrt{3}$ holds, then we have

Then, to let $\Delta_{\mathcal F} < 0$, we need

If $\dfrac{\rho_1^2 \lambda}{\rho_2\mu^2} > 8+4\sqrt{3}$ or $1 < \dfrac{\rho_1^2 \lambda}{\rho_2\mu^2} < 8-4\sqrt{3}$ holds, then we have

Then, to $\Delta_{\mathcal F} < 0$, we need

So far, by Theorem 2.1 of ^[14], we are in the position to state the following results on Hopf bifurcations:

Theorem 3.4. Suppose that (3.2) holds and that

Then, for any $\ell > \ell_1$ ($\ell_n$ is defined in (3.5)), there exists $n$ points $\lambda_j^H(\ell)$, $1\leq j\leq n$, satisfying

such that the reaction-diffusion system undergoes a Hopf bifurcation at $\lambda = \lambda_j^H$ or $\lambda = \lambda_0$ and the following holds true:

1. The bifurcating periodic solutions from $\lambda = \lambda_0^H$ are spatially homogeneous and unstable, which coincides with the periodic solution of the corresponding ODE system;

2. The bifurcating periodic solutions from $\lambda = \lambda_{j}^H$ are spatially non-homogeneous and unstable.

Remark 3.5. 1. For the bifurcating spatially homogeneous periodic solutions, it is always unstable in the reaction-diffusion system given by System (3.1) since it is unstable in the corresponding ODE system given by System (2.1). On the other hand, the bifurcating spatially non-homogeneous periodic solutions from $\lambda = \lambda_{j}^H$ are always unstable, since, under (3.2), $(u_ \lambda, v_ \lambda)$ is unstable with respect to the reaction-diffusion system given by System (3.1).

2. Biological meaning: Our results showed that, with the inclusion of the spatial diffusion, as the decay rate of $v$ changes, the reaction-diffusion system given by System (3.1) might not only undergo temporal oscillations, but also spatiotemporal oscillations with the evolution of time. That is, the densities of the activator $u$ and the inhibitor $v$ will oscillate with respect to $t$, but they will also depend on the spatial variable $x$ (in the case of spatially non-homogeneous periodic solutions).

3.3.   Steady-state bifurcation analysis: Non-constant positive equilibrium solutions

In this subsection, we shall use the steady-state bifurcation theory to show the existence of the bifurcating non-constant positive equilibrium solutions.

Following ^[14], we identify steady-state bifurcation values, denoted by $\lambda_S$, which satisfy the following: there exists $n\in \mathbb{N}_0$ such that

and that

By Theorems 3.1 and 3.2, for any $\lambda\in\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}, +\infty\bigg)$, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). Hence, any potential bifurcation point $\lambda_S$ must be in the interval $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$. We shall study the steady-state bifurcation points in this interval.

We rewrite $D_n(\lambda)$ in the form of

where $p: = n^2/\ell^2$ and

Solving $p$ from $D_n(\lambda) = 0$, we have

where

By direct calculation, we have

where

Clearly, $g_1(\lambda) = 0$ has a unique positive root, denoted by $\lambda_{g_1}$. In particular, $g_1(\lambda) < 0$ for $\lambda\in(0, \lambda_{g_1})$, while $g_1(\lambda) > 0$ for $\lambda\in(\lambda_{g_1}, +\infty)$. On the other hand, since $g_2(0) = -4d_2\rho_2^2\mu^5 < 0$ and $g_2(+\infty) = +\infty$, $g_2(\lambda) = 0$ has at least one positive root in $(0, +\infty)$. In particular, $g_2(\lambda) = 0$ may have one unique positive root or three positive roots in $(0, +\infty)$.

In what follows, for our particular interests, we assume that $g_2(\lambda) = 0$ has a unique positive root in $(0, +\infty)$. (The case of $g_2(\lambda) = 0$ having more than one positive root in $(0, +\infty)$ is much more complicated and will be considered in our future investigations.)

We have the following results on the properties of $h(\lambda)$ in the interval $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$:

Lemma 3.6. Let $\lambda_{g_1}$ be the unique positive root of $g_1(\lambda) = 0$. Assume, in addition that, $g_2(\lambda) = 0$ also has a unique root, denoted by $\lambda_{g_2}$ in $(0, +\infty)$. Then, the following conclusions hold true:

1. Suppose that $0 < \dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ holds. Then,

(a) $\lambda_{g_1} < \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_{g_2} < \dfrac{2\rho_2\mu^2}{\rho_1^2}$;

(b) for any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $h(\lambda) > 0$ and

(c) $h(\lambda)$ is decreasing in $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{g_2}\bigg)$, while it is increasing in $\lambda\in\bigg(\lambda_{g_2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$; at $\lambda = \lambda_{g_2}$, $h(\lambda)$ attains its positive minimum value $h(\lambda_{g_2})$.

2. Suppose that $\dfrac{d_2}{d_1} > \dfrac{2\rho_2\mu}{\rho_1^2}$ holds. Then,

(a) $\lambda_{g_2} < \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_{g_1} < \dfrac{2\rho_2\mu^2}{\rho_1^2}$;

(b) for any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $h(\lambda) > 0$ and

(c) $h(\lambda)$ is decreasing in $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{g_1}\bigg)$, while it is increasing in $\lambda\in\bigg(\lambda_{g_1}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$; at $\lambda = \lambda_{g_1}$, $h(\lambda)$ attains its positive minimum value $h(\lambda_{g_1})$.

Proof. The proof is moved to the Appendix.

Define

We have the following results on the sign of $\Sigma_{\mathcal D}(\lambda)$ in the interval $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$:

Lemma 3.7. Let $\lambda_{g_1}$ be the unique positive root of $g_1(\lambda) = 0$. Assume, in addition that, $g_2(\lambda) = 0$ also has a unique root, denoted by $\lambda_{g_2}$ in $(0, +\infty)$. Then, the following conclusions hold true:

1. Suppose that $0 < \dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ holds.

(a) if $\dfrac{4\mu}{h(\lambda_{g2})} < \dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ holds, then, for all $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $\Sigma_{\mathcal D}(\lambda) > 0$;

(b) if $\dfrac{d_2}{d_1} < \min\bigg(\dfrac{4\mu}{h(\lambda_{g2})}, \dfrac{2\rho_2\mu}{\rho_1^2}\bigg)$ holds, then there exist $\underline{ \lambda}, \overline{ \lambda}\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, with $\underline{ \lambda} < \overline{ \lambda}$, such that

and, for any $\lambda\in(\underline{ \lambda}, \overline{ \lambda})$, $\Sigma_{\mathcal D}(\lambda) < 0$, while, for $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \underline{ \lambda}\bigg]\cup \bigg[\overline{ \lambda}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $\Sigma_{\mathcal D}(\lambda) > 0$.

2. Suppose that $\dfrac{d_2}{d_1} > \dfrac{2\rho_2\mu}{\rho_1^2}$ holds. Then, there exists a unique point $\lambda_{*}\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$ such that $h(\lambda_{*}) = 4\mu\bigg(\dfrac{d_1}{d_2}\bigg)$. In particular, $\Sigma_{\mathcal D}(\lambda) < 0$ for any $\lambda\in\bigg(\lambda_{*}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, while $\Sigma_{\mathcal D}(\lambda) > 0$ for any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$.

Proof. The proof is moved to the Appendix.

For clarity of our later discussions, we divide our discussions into the following cases:

Case 1: $\dfrac{4\mu}{h(\lambda_{g2})} < \dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ and $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$;

Case 2: $\dfrac{d_2}{d_1} < \min\bigg(\dfrac{4\mu}{h(\lambda_{g2})}, \dfrac{2\rho_2\mu}{\rho_1^2}\bigg)$ and $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \underline{ \lambda}\bigg]\cup \bigg[\overline{ \lambda}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$;

Case 3: $\dfrac{d_2}{d_1} > \dfrac{2\rho_2\mu}{\rho_1^2}$ and $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$, where $\lambda_*$ is defined in Lemma 3.7.

Theorem 3.8. Suppose that either Case 1 or 2 holds. Then, no steady-state bifurcation around $(u_ \lambda, v_ \lambda)$ occurs.

Proof. Suppose that either Case 1 or 2 holds. Then, $p_{\pm}(\lambda)$ is well defined. Clearly, $\dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$. Then, by (a) of Lemma 3.6, we have that $\lambda_{g_1} < \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_{g_2} < \dfrac{2\rho_2\mu^2}{\rho_1^2}$. One can check that $A(\lambda)$ is decreasing since

Since $A(\lambda_{g_1}) = 0$, it follows that, for all $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $A(\lambda) < 0$. Thus, $p_{\pm}(\lambda) < 0$ whenever they are well defined. In this case, no steady-state bifurcation around $(u_ \lambda, v_ \lambda)$ occurs.

Lemma 3.9. Suppose that Case 3 holds. Then, there exists $\lambda_{*}\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, with $\lambda_{*} < \lambda_{g_1}$, such that $p_{+}(\lambda)$ is decreasing while $p_{-}(\lambda)$ is increasing in $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$. In particular,

and $\lim\limits_{ \lambda\rightarrow \lambda{*}}p'_{-}(\lambda) = +\infty$ and $\lim\limits_{ \lambda\rightarrow \lambda{*}}p'_{+}(\lambda) = -\infty$.

Proof. The proof is moved to the Appendix.

Therefore, for any $n > 0$, if $p_{-}\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}\bigg) < n^2/\ell^2 < p_+\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}\bigg)$, then there exists $\lambda_{n}^S\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_*\bigg)$ such that $p_{-}(\lambda_n^S) = 0$ or $p_+(\lambda_n^S) = 0$.

Define

Then, for any $\ell\in(\ell_n^+, \ell_n^-)$, there exists $\lambda_n^S$ such that $D_n(\lambda_n^S) = 0$. These points $\lambda_n^S$ are potential steady-state points. As remarked in ^[14], however, it is possible that, for some $i < j$, $p_{-}(\lambda_i^S) = p_+(\lambda_j^S)$. In this case, for $\lambda = \lambda_i^S = \lambda_j^S$, $0$ is not a simple eigenvalue of $L(\lambda)$, and we shall consider bifurcations at such points. Let $L^E$ be the set of such points.

Summarizing the aforementioned observations, by Theorem 3.10 of ^[14], we have the following results on the existence of steady-state bifurcations:

Theorem 3.10. Suppose that Case 3 holds. Let $\ell_n^+$ and $\ell_n^-$ be defined in (3.20). If, for some $n\in N$, $\ell\in(\ell_n^+, \ell_n^-)\backslash\{L^E\}$ and there exists a point $\lambda_{n}^S\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_*\bigg)$ such that $p_+(\lambda_{n}^S) = 0$ or $p_{-}(\lambda_{n}^S) = 0$. Then, there is a smooth curve of positive solutions of the reaction-diffusion system bifurcating from $(\lambda, u, v) = (\lambda_n^S, u_{ \lambda_n^S}, v_{ \lambda_n^S})$.

Remark 3.11. Biological meaning: Our results showed that, for any $\ell > 0$ (the length of the spatial domain), there exists an $n > 0$ such that, for a suitable $\lambda$ (the decay rate of the inhibitor $v$) and suitable diffusion rates $d_1$ and $d_2$, the reaction-diffusion system given by System (3.1) might have positive non-constant steady-state solutions with the eigen-mode $\cos(nx/\ell)$; That is, the densities of the activator $u$ and the inhibitor $v$ have a non-uniform spatial distribution in $\Omega$. From the viewpoint of pattern dynamics, in this case, System (3.1) will undergo spatial patterns which are different from Turing patterns.

4.   Conclusions

In this study, we were mainly concerned with the spatiotemporal patterns and multiple bifurcations of a reaction-diffusion model for hair follicle spacing.

First, we consider the stability and instability of the equilibrium solution of the ODE system. In particular, by using the center manifold theory, normal form methods, as well as the standard Hopf bifurcation theory, we were able to prove the existence of the Hopf bifurcating periodic solutions bifurcating from the equilibrium solution. By calculating the first Lyapunov coefficient, we found that the Hopf bifurcating periodic solutions are always unstable. This is one of the novel points of this paper. Since the bifurcating periodic solutions are unstable in ODEs, they will never undergo Turing instability (see ^[9,15,16] for more details on the Turing instability of periodic solutions) in the reaction-diffusion equations given by System (3.1). This is different from the equilibrium solutions, which may experience Turing instability under certain suitable conditions on the diffusion rates ($d_1$ and $d_2$) and the decay rate of $v$. This was shown in the analysis of the reaction-diffusion system.

Second, we studied the stability and instability of the constant equilibrium solution in the reaction-diffusion system. In particular, by using Hopf bifurcation theory and steady-state bifurcation theory, we were capable of showing the existence of spatially non-homogeneous Hopf bifurcating periodic solutions, as well as the non-constant bifurcating steady-state solutions for the reaction-diffusion equations. Moreover, Turing instability of the constant equilibrium solution was investigated in detail. If Turing instability of the constant equilibrium solutions occurs, then Turing patterns emerge.

However, in application, for many reaction-diffusion equations, Turing-Hopf bifurcations can also be observed. In this study, we did not consider the existence of Turing-Hopf bifurcations. Indeed, the analytical analysis of Turing-Hopf bifurcation will be much more difficult than that of either Turing bifurcation or Hopf bifurcation. We will consider Turing-Hopf bifurcations of this particular reaction-diffusion system in our future investigations.

Regarding the numerical simulations, we would like to mention that we did not include numerical simulations in the paper. The reasons are as follows: from our analytical analysis, it was found that the Hopf bifurating periodic solutions are always unstable, so it is hard to simulate numerically; on the other hand, for the steady-state bifurcations, under our conditions, the constant equilibrium solutions and bifurcating steady-state solutions are unstable. Again, it is hard to simulate numerically.

Finally, we ought to remark that our results in the current paper tend to be much more analytical. Although it will make contributions to the field of bifurcation theory with applications, as well as to the field of mathematical biology, further work needs to be done to use our analytical results to understand the factors that affect the hair follicle spacing. We shall study it in our future investigations.

Acknowledgments

Zhili Zhang, Aying Wan and Hongyan Lin were partially supported by the National Natural Science Foundation of China (12061033).

Conflict of interest

The authors declare that there is no conflict of interest.

A.   Appendix

1. Proof of Theorem 3.1. To understand the stability and instability of $(u_ \lambda, v_ \lambda)$, we need to know the linearized operator of System (3.1) evaluated at $(u_{ \lambda}, v_{ \lambda})$, which is given by

It is well known that (see ^{[14,17,18,19,20]}) the eigenvalue problem

has eigenvalues $\widehat\mu_n = \frac{n^2}{\ell^2}$ ($n = 0, 1, 2, \cdots$), with corresponding eigenfunctions $\varphi_n(x) = \cos \frac{n}{\ell}x$.

Let

be an eigenfunction for $L(\lambda)$ with the eigenvalue $\gamma(\lambda)$. Then, we have

where

Then, the eigenvalues of $L(\lambda)$ are given by the eigenvalues of $L_n(\lambda)$ for $n = 0, 1, 2, \cdots$.

Indeed, the characteristic equation of $L_n(\lambda)$ is

where

Note that $\mu\geq\dfrac{\rho_1^2}{2\rho_2}$ or, equivalently, $\lambda_0\leq \dfrac{\rho_2\mu^2}{\rho_1^2}$ holds; we have that $T_n(\lambda) < 0$ for any $\lambda > \dfrac{\rho_2\mu^2}{\rho_1^2}$ and $n\in\mathbb N_0$ since $T(\lambda) < 0$. We now consider the sign of $D_n(\lambda)$, with $n\in\mathbb N_0\backslash\{0\}$:

Suppose that $\lambda\in\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}, +\infty\bigg)$ holds. Then, for any $n\in\mathbb N_0$, we have that $D_n(\lambda) > 0$. Thus, in this case, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). By Theorem 2.1, $(u_ \lambda, v_ \lambda)$ is also locally asymptotically stable with respect to the ODE system. Thus, Turing instability of $(u_ \lambda, v_ \lambda)$ does not occur.

Suppose that $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$ holds.

If, in addition, $d_1\geq\dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then, for any $n\in\mathbb N_0$, we have

In this case, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). Again, Turing instability of $(u_ \lambda, v_ \lambda)$ does not occur.

If, in addition $d_1 < \dfrac{(2\rho_2\mu^2-\rho_1^2 \lambda)\mu\ell^2}{\rho_1^2 \lambda}$ holds, then, by the definition of $n_{d_1}$, for any integer $n > n_{d_1}$, we have that $D_n(\lambda) > 0$. Moreover, for any $d_2 < \widehat d_2$, $D_n(\lambda) > 0$ with $n\in\mathbb N_0\backslash\{0\}$. In this case, $(u_ \lambda, v_ \lambda)$ is locally asymptotically stable with respect to the diffusive system given by System (3.1). If $d_2 > \widehat d_2$ holds, then there exists at least an $n\in[1, n_{d_1}]$ such that $D_n(\lambda) < 0$. In this case, $(u_ \lambda, v_ \lambda)$ is unstable with respect to the diffusive system given by System (3.1). Thus, $(u_ \lambda, v_ \lambda)$ undergoes Turing instability in the reaction-diffusion system given by System (3.1). We thus complete the proof.

2. Proof of Lemma 3.6. We only prove the first part, since the proof of the second part is similar. Clearly, $g_1(\lambda)$ is increasing in $(0, +\infty)$. By direct calculation, we have

If $\dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ holds, then, by (A.3), $g_1\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}\bigg) > 0$ and $g_2\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}\bigg) < 0$. This implies that $\lambda_{g_1} < \dfrac{\rho_2\mu^2}{\rho_1^2} < \lambda_{g_2} < \dfrac{2\rho_2\mu^2}{\rho_1^2}$. In this case, the monotonic properties of $h(\lambda)$ can be obtained. We thus complete the proof.

3. Proof of Lemma 3.7. We only prove Part 1, since the other parts can be proved similarly. Assume that $0 < \dfrac{d_2}{d_1} < \dfrac{2\rho_2\mu}{\rho_1^2}$ holds. Then, we have that $4\mu\bigg(\dfrac{d_1}{d_2}\bigg) < h\bigg(\dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$.

If, in addition, $4\mu\bigg(\dfrac{d_1}{d_2}\bigg) < h(\lambda_{g_2})$ or, equivalently that $\dfrac{d_2}{d_1} > \dfrac{4\mu}{h(\lambda_{g_2})}$ holds, then, for any $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, $\Sigma_{\mathcal D}(\lambda) > 0$.

If, in addition, $\dfrac{d_2}{d_1} < \min\bigg(\dfrac{4\mu}{h(\lambda_{g2})}, \dfrac{2\rho_2\mu}{\rho_1^2}\bigg)$ holds, then there exist $\underline{ \lambda}, \overline{ \lambda}\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$, with $\underline{ \lambda} < \overline{ \lambda}$, such that

and $\Sigma_{\mathcal D}(\lambda) < 0$ for any $\lambda\in(\underline{ \lambda}, \overline{ \lambda})$, while $\Sigma_{\mathcal D}(\lambda) > 0$ for $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \underline{ \lambda}\bigg]\cup \bigg[\overline{ \lambda}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$. We thus complete the proof.

4. Proof of Lemma 3.9. Suppose that Case 3 holds. By Lemma 3.7, $p_{\pm}(\lambda)$ is well defined in $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$, where $\lambda_{*}\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{g_1}\bigg)$. Thus, in this interval, $p_{\pm}(\lambda) > 0$.

Since $A(\lambda)$ is decreasing, $C(\lambda)$ is increasing in $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \dfrac{2\rho_2\mu^2}{\rho_1^2}\bigg)$ and $\dfrac{d_2}{d_1} > \dfrac{2\rho_2\mu}{\rho_1^2}$; by

we know that $p_{+}(\lambda)$ is decreasing in $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$.

On the other hand, by $p_{+}(\lambda)p_{-}(\lambda) = \dfrac{\mu{C(\lambda)}}{4d_1d_2}$, we know that $p_{-}(\lambda)$ is increasing in $\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$, since $p_{+}(\lambda)$ is decreasing and $C(\lambda)$ is increasing. By direct calculation, it follows that

Since $A'(\lambda) < 0$, we can obtain that $\lim\limits_{ \lambda\rightarrow \lambda{*}}p'_{-}(\lambda) = +\infty$ and $\lim\limits_{ \lambda\rightarrow \lambda{*}}p'_{+}(\lambda) = -\infty$.

Noticing that $h(\lambda_{*}): = \frac{A^2(\lambda_{*})}{C(\lambda_{*})} = 4\mu\dfrac{d_1}{d_2}$, we have

Since $p_{-}(\lambda)$ is increasing and $p_{+}(\lambda)$ is decreasing, we have that $p_{-}(\lambda) < p_{-}(\lambda_{*}) = p_{+}(\lambda_{*}) < p_{+}(\lambda)$ for all $\lambda\in\bigg(\dfrac{\rho_2\mu^2}{\rho_1^2}, \lambda_{*}\bigg]$. In summary, we have (3.19). We thus complete the proof.