QUALITATIVE ANALYSIS ON POSITIVE STEADY-STATES FOR AN AUTOCATALYTIC REACTION MODEL IN THERMODYNAMICS

. In this paper, a reaction-diﬀusion system known as an autocatalytic reaction model is considered. The model is characterized by a system of two diﬀerential equations which describe a type of complex biochemical reaction. Firstly, some basic characterizations of steady-state solutions of the model are presented. And then, the stability of positive constant steady-state solution and the non-existence, existence of non-constant positive steady-state solutions are discussed. Meanwhile, the bifurcation solution which emanates from positive constant steady-state is investigated, and the global analysis to the system is given in one dimensional case. Finally, a few numerical examples are provided to illustrate some corresponding analytic results.


1.
Introduction. In recent decades, researchers paid much attention to the resemble behaviors between dynamical instability and the equilibrium phase changing on the steady-states of certain systems when they were far away from the thermodynamic equilibrium. For example, they emerge tremendous fluctuation near the critical points, their dynamical behaviors may occur the phenomenon of critical slowing down, and then the system may give rise to critical transition to pattern formation when the phases are changed [3,24]. These phenomena enlighten people to give some specific classifications on critical behaviors of thermodynamic systems by some mathematical techniques, such as bifurcation theory, degree theory and perturbation theory, etc.
In general, most dynamical behaviors in chemical reactions are nonlinear. However, when systems approach to the thermodynamic equilibrium, their kinetic behaviors may be studied by linear non-equilibrium thermodynamics. When systems are far away from the thermodynamic equilibrium, sometimes, the effects of the nonlinear terms may become the principal factors. These nonlinear phenomena can couple with linear diffusive behaviors, and then the spontaneous appearances of the order and chaotic pattern, that is, the chemical oscillation may occur. For these reasons, chemical reactions which display oscillatory solutions and steady-state solutions have been of great interests for both theoreticians and experimentalists for over many years. Some experimental examples of oscillatory behavior in chemical systems including Bray-Liebhafsky, Belousov-Zhabotinsky and Briggs-Rausher reactions, for which periodic variations in concentrations can be visualized via changes in color, see [5] for a review of these reactions and other oscillatory phenomena. To study the chemical oscillation in chemical reaction-diffusion systems, it is important to master the kinetic generalities of these systems near the critical points. These are the main concerns of autocatalytic reaction models, which have been applied to various problems in chemistry and biology, see [10,11,13,15,21,22,29,30] and the references therein for example.
The autocatalytic reaction is a large class of nonlinear reactions most closely associated with order creation. These are reactions in which one or more of the products are the same as one or more of the reactants, that is, the presence of one or more chemical substances stimulate the production of more of them. In the process of reaction, it is reasonable to assume that the chemical substance diffuses with some diffusivity, and the reaction rate is affected not only by reactant concentration, but also by reaction product concentration. It is characterized by: (i) the autocatalytic reaction is initiated by adding a tiny amounts of products, (ii) the reaction is very slow at the beginning, and then the reaction rate is accelerated as the accumulation of reaction products, after that, the reaction rate will go down as the consumption of reaction products. For example, the process of fermentation of some substances is a typical autocatalytic reaction, NO 2 produced in thermal decomposition of gunpowder is able to take catalytic effect to the reaction. Some other examples, such as in the process of DNA replication, Haloform reaction, permanganate with oxalic acid reaction and in the process of α-bromination of acetophenone with bromine, there arise quantities of autocatalytic reaction.
Simple autocatalytic reactions are known to oscillate in time, in which the temporal order can be created. Perhaps, the simplest autocatalytic reaction can be written as A + B → 2B, where A and B represent the reactant and autocatalyst, respectively. This reaction is one in which a molecule of species A interacts with a molecule of species B. An A molecule is converted into a B molecule. The final product consists of the original B molecule plusing the B molecule created in the reaction. Other simple reactions can generate spatial separation of chemical species generating spatial order.
A slightly more complicated autocatalytic reaction is given by the following where P is the density of the reaction precursor, A and B also represent the reactant and autocatalyst, C is some certain inert product. The supply of reactant is proceeded by the gradually degenerative process of the precursor. The reaction can proceed well at normal temperature. It can be assumed that the reactant is effectively immobilised within the reactor, and the autocatalyst is made to flow through the reactor with a constant velocity as well as being able to diffuse, see [2] for details. Much of the previous work has been concerned with the iodate-arsenous acid systems, for which the autocatalysis is a good approximation in the arsenous acid excess case. For iodate-arsenous acid system, the reactant A and the autocatalyst B represent, respectively, IO − 3 and I − [8,20]. The basic idea of this type reaction is that the diffusion of the autocatalyst has a stabilizing effect on the planar waves, whereas the diffusion of the substrate has a destabilizing effect. Thus if the latter effect is sufficiently strong relative to the first, then the wave will become transversely unstable. This leads to a critical value of diffusion rate at which the stability of a planar wave changes [14,18]. More complex autocatalytic reactions are involved in metabolic pathways and metabolic networks in biological systems.
In the present paper, we consider a cubic autocatalytic reaction-diffusion model from the mathematical point of view. The reaction is a type of biochemical reaction in which the reaction products can increase the reaction rate. The simplified and dimensionless form of autocatalytic reaction biochemical model based on cubic reaction term is as the following where Ω ⊂ R n is bounded with smooth boundary ∂Ω, the variables u and v represent the dimensionless concentrations or densities of the reactant and autocatalyst respectively and so are usually assumed to be non-negative. a is a parameter representing the initial concentration of the reaction precursor, d 1 and d 2 are the diffusion rates of u and v respectively, a, d 1 and d 2 are positive constants. ∂ ∂ν denotes the differential in the direction of the outer normal to ∂Ω. The chemical reaction is driven within a closed container.
The key feature of equations in (1) is that they are nonlinear, the term on the right varies as the square of the concentration of v. This feature can lead to multiple fixed points of the system, much like a quadratic equation can have multiple roots. Physically, multiple fixed points allow for multiple states of the system. A system existing in multiple macroscopic states is more orderly than a system in a single state. From this perspective, system (1) seems complicated and particular in some extent.
Our main aim in this paper is to investigate the steady-state solutions of system (1) mathematically, that is, the positive classical solutions of the following elliptic system The idea in constructing this paper is partly due to the techniques developed in [23]. An outline of the paper is as follows: In Section 2, we give some basic properties of non-homogeneous steady-state solutions of (1), including a priori estimate (since the boundedness of positive solutions seems important in latter discussions).
In Section 3, we discuss the stability of the unique constant steady-state solution (u * , v * ) = 1 a , a in detail. In Section 4, we establish non-existence results of nonconstant positive solutions for large or small diffusion rate. In Section 5, by using the fixed-point index theory in Banach space and linear operator theory in functional analysis, we investigate the existence of non-constant positive solutions. For fixed d 1 , d 2 , in Section 6, by taking a as a parameter, we analyze the bifurcation solution which emanates from the constant solution 1 a , a . Meanwhile, in one dimensional case, we give some global analysis to the system. Section 7 aims at some numerical examples in verifying our analytic results.
2. Some characterizations of positive solutions. This section is devoted to some basic properties of non-homogeneous steady-state solutions of (1), namely, positive solutions of (2)-(4). We first set out to seek for the a priori estimate of positive solutions of (2)-(4).
The following two lemmas derived from [16] (Lemma 2.1 and Lemma 2.3) are the main tools in obtaining our estimate. Where the proof of first lemma is based on the result which is well known as a local result for weak super-solutions of linear elliptic equations (see, for example, [9], Theorem 8.18), while the second one is an analogue of Proposition 2.2 in [17]. On the methods in proving the boundedness of positive solutions of a model, the readers can also refer to [7,27], where the authors proved the boundedness of positive solutions of a so-called Sel'kov model in detail, they gave their estimates for the positive solutions by restricting the dimension of Ω. In the present paper, mainly using the techniques introduced in [16] and Proposition 2.2 in [17], we give the a priori estimate of positive solutions of (2)-(4).
Lemma 2.1. [16] Let Ω be a bounded Lipschitz domain in R n . Let Λ be a nonnegative constant and suppose that v ∈ W 1,2 (Ω) is a nonnegative weak solution of inequalities Then, for any q ∈ [1, n n−2 ), there is a constant C 0 , determined only by q, Λ and Ω, Lemma 2.2. [16] Let Ω be a bounded Lipschitz domain in R n and g ∈ C(Ω × R).
Suppose that there is a constant M such that g(x, z) < 0 for z > M . If w ∈ W 1,2 (Ω) is a weak solution of inequalities We now consider nonnegative solutions of system (2)-(4). Our a priori estimate has the following form.
Remark 1. We note that the average value of v on Ω is a, this leads to inf Ω v ≤ a. Therefore a C0 |Ω| ≤ a in view of v ≥ a C0 |Ω|. This shows that the estimate inequality on v is valid, and therefore the inequality on u also holds deservedly.

Remark 2.
Obviously, the estimate inequalities on u and v are valid if D 1 > ε and D 2 < M for some small ε > 0 and large M > 0.
We prove Theorem 2.3 without restricting the dimension of Ω. Thus, in the forthcoming discussion of the paper, if it involves the boundedness of the solution of (2)-(4), we deservedly think that the solution is bounded without considering the dimension of Ω. Now, we study some other basic properties of positive solutions of (2)-(4). The main tools we used here are integration and a few well-known inequalities.
Proof. (1) The result follows immediately from the proof of Theorem 2.3.
In the following, we give a few results on ξ and η (Theorem 2.5). We denote by µ 1 the second eigenvalue of the operator −∆ with homogenous Neumann boundary condition in the sequel of this paper.
Theorem 2.5. Assume that (u, v) is a non-constant positive solution of (2)-(4). Then the following integral inequalities hold.
(1) Ω ξηdx < 0 and Ω ∇ξ∇ηdx < 0. (2) (3) There are positive constants C 1 := C 1 (µ 1 , Ω), (2) adds to (3) which leads to ∆w = v − a. Further, we have Multiplying (5) by w = d 1 u + d 2 v and then integrating the result over Ω, we obtain This shows that Multiplying (5) by η and then integrating over Ω, we get which means that (2) Multiply ∆w = v − a by ξ + η and then integrate the result over Ω, then the left is And the right is (3) It suffices to prove the first inequality, and the second one can be proved similarly.

Multiplying (3) by η and then integrating over Ω, we obtain
Use the Hölder inequality, it then follows that Combining the Poincaré inequality with (8) we have Further, it follows from the above that Ω |∇η| 2 dx ≤ C 2 |Ω| By (7) we have Thus, it follows that Therefore we have

QUALITATIVE ANALYSIS FOR AN AUTOCATALYTIC REACTION MODEL 4793
On the other hand, for any b, c ∈ R and s > 0, using the Cauchy inequality we get Especially, if we take s = 1 µ1 , then bc ≤ µ1 It follows that 3 4 It then follows by (4) that .
Moreover, if d 2 ≤ 2, then (9) implies that Therefore, we get The proof is completed.
3. Stability of positive constant solution. In this section, we discuss the stability of constant positive solution of (2)-(4). As with the boundedness of positive solutions established in Theorem 2.3, the stability result also takes an important role in Section 5 in obtaining the existence theorem.
Then system (1) transforms into the following form Considering the Turing instability [26] of the constant solution, we take into account the following ordinary differential system Clearly, (u * , v * ) is also the unique constant solution of (13). The Jacobian matrix of (13) Since the Jacobian determinant Det(J) = a 2 > 0, the eigenvalues of J have negative real parts if and only if the trace of J satisfies Tr(J) = 1 − a 2 < 0. In this case, (u * , v * ) is the stable solution of (13). Now, we denote by 0 = µ 0 < µ 1 ≤ µ 2 ≤ µ 3 ≤ · · · all eigenvalues of the operator −∆ subject to homogenous Neumann boundary condition and τ i ≥ 1 the multiplicity of each µ i , i = 0, 1, 2, · · · . Let ϕ ij , 1 ≤ j ≤ τ i be the normalized eigenfunction corresponding to µ i . Then the sequence {ϕ ij }, i ≥ 0, 1 ≤ j ≤ τ i constitutes a complete orthonormal basis in L 2 (Ω). In the following, we investigate the stability of (u * , v * ).
Then we have Clearly, λ is an eigenvalue of L d1d2 if and only if for some i ≥ 0, the determinant of the matrix L d1d2 − λI is equal to 0. Simple calculation gives for all i ≥ 0. Clearly, Q 0 > 0. Thus, if a > 1 and µ 1 ≥ 1 d2 , then we must have P i , Q i > 0 for all i ≥ 0. Therefore, the real part Reλ of λ satisfies Reλ < 0, and in this case, the constant steady-state solution (u * , v * ) is asymptotically stable.
The stability of positive constant solution reads as follows.
The statement is obvious. Indeed, when a < 1, we have at least P 0 < 0, so L d1d2 has eigenvalues whose real parts are positive.
1+µd2 . Then It is easy to check that f (µ) has a unique stationary point µ = is increasing in (0, ) and decreasing in ( , ∞). f (µ) has two zero points 0 and 1 d2 in (0, ∞). Then the stable results in this part follow. As to the instability results, we have Q 1 (resp. Q i0 ) < 0 under our assumptions. Therefore L d1d2 has eigenvalues whose real parts are positive.
for some i. So when a > 1, it can be ensured that (u * , v * ) is asymptotically stable if we take d 2 suitably large or d 1 suitably small.

Remark 4.
We see that the constant steady-state solution (u * , v * ) is unstable for both systems (11)-(12) and (13) if a < 1. For a > 1, (u * , v * ) is stable for (13), but may be unstable for (11)- (12). Therefore, applying the implication of Turing instability, for the corresponding parabolic system (1), the Turing instability of the constant solution (u * , v * ) may occur if a > 1. This implies that it is just diffusion that causes the occurrence of instability.
Remark 5. In view of the unstability of (u * , v * ) for system (13) for a < 1, here, by the way, we intend to consider the Hopf bifurcation which would bifurcate from (u * , v * ) for system (13) by using the Hopf bifurcation theory [12,19]. We restrict 0 < a ≤ 1.
, then the two eigenvalues of linearized operator of (13) at (u * , v * ) are where i is the imaginary unit. For a = 1, the two eigenvalues are a pair of pure imaginaries We denote by Re(a) the real part of the eigenvalue of linearized operator of (13) at (u * , v * ). Though the derivative Re (a)| a=1 = −1, it is easy to see that the real parts of λ + , λ − are positive for √ 2 − 1 < a < 1. Therefore, there is no Hopf bifurcation in this case either.
The discussion above shows that the local system (13) has no periodic solution arising from Hopf bifurcation. Therefore, system (13) does not permit any branch of periodic solution bifurcating from (u * , v * ) for 0 < a ≤ 1 (However, we will see in Section 7 that the periodic solutions may emerge for system (1)).
4. Non-existence of non-constant positive solution. This section is concerned with the non-existence of non-constant positive solution of (2)-(4). We prove that (2)-(4) has no non-constant positive solution if the effective diffusion rates are suitably chosen. Here we note that the non-existence results obtained in this section play a critical role in Section 5 in obtaining the existence theorem. Considering Theorem 2.3, for simplicity and convenience in use, we denote by  Proof. Suppose on the contrary that (u, v) is a non-constant positive solution of (2)-(4). Multiplying (3) by η and then integrating over Ω, by Theorem 2.3, (1) of Theorem 2.5, the Poincaré inequality and the fact that Ω ηdx = 0, we obtain which leads to a contradiction on the condition that d 2 ≥ d 3 . Therefore, (2)-(4) does not admit a non-constant positive solution when d 2 ≥ d 3 .
Using the similarly techniques as above, we can prove following results. Proof. Let (u, v) be a non-constant positive solution of (2)-(4). Multiplying (2) by ξ and then integrating the result over Ω, by Theorem 2.3, the Cauchy inequality and the fact that Ω ξdx = 0, we obtain (Note that here we take s = a|Ω| C0 2 in applying the Cauchy inequality (10)).
However, by (4) of Theorem 2.5, we get Ω |∇η| 2 dx <   Proof. Let (u, v) be a non-constant positive solution of (2)-(4). Combining the equality with the proofs of Theorem 4.1 and Theorem 4.2, by Theorem 2.3, (1) of Theorem 2.5, the Cauchy inequality and the Poincaré inequality, for some small s > 0, we which leads to a contradiction provided that d 1 ≥ d 5 , d 2 ≥ d 6 . Therefore, (2)-(4) does not admit a non-constant positive solution under our assumptions.

5.
Existence of non-constant positive solutions. In this and the subsequent sections, we investigate the existence of non-constant positive solutions of (2)-(4). Typically, there are two methods in establishing the existence of non-trivial solutions for elliptic systems, one is a bifurcation technique that we will use in the sequel, and the other is a singular perturbation. On bifurcation technique, a variation procedure is effective through a so-called fixed-point index theory. In this section, by using the fixed-point index theory developed in [1] for compact operators in Banach space, we consider the existence of non-constant positive solutions of (2)-(4). Before establishing the main result of this section, we first introduce some notations and preliminaries.
Then (2) and (3) shift into the forms where f 3 and f 4 are higher order terms of u and v with Thus, the constant steady-state solution (u * , v * ) transforms into (0, 0). For w(x) ∈ C α (Ω), α ∈ (0, 1), assume that are solutions of problems respectively. Then (u, v) ∈ C 2+α (Ω) × C 2+α (Ω) is unique, and the operators G d1 , G d2 are continuous and compact. Denote Then (14)- (15) can be interpreted as the equation where and H are compact on E with For some small ε > 0, define Here C 3 is the same as that of in the proof of (3) of Theorem 2.5. Our existence result reads as follows.
Theorem 5.1. Suppose that µ 1 < 1 d2 < µ 2 and the multiplicity of µ 1 is odd. Then Proof. By Theorem 2.3 and the definition of ( u, v), we see that the equation (16) has no solution on the boundary ∂S of S. Using the homotopy invariance on degree [1], we know that deg( where σ is the sum of the algebraic multiplicities of all positive eigenvalues of Suppose that ϕ ψ is an eigenfunction of K(d 1 , d 2 ) − I corresponding to eigenvalue µ. Then we have Thus, all eigenvalues of K(d 1 , d 2 )−I exactly consist of all roots of the characteristic equation The constant term of (19) can be reformed as By the nonexistence results in Section 4, we can choose d 1 , d 2 suitably such that (2)-(3) has no non-constant positive solution, further, (16) has no nonnegative solution in S except (0, 0). And by Theorem 3.1, we can also take d 1 , d 2 suitably such that K(d 1 , d 2 ) − I has no nonnegative eigenvalue, that is, K(d 1 , d 2 ) − I is invertible. Therefore Since , then it is easy to check that So the characteristic equation (19) has a negative root and a positive root, say, λ 1 for i = 1. On the other hand, by 1 d2 < µ 2 , we know that Q i > 0 for all i = 1. Therefore, when i = 1, (19) has no non-negative real part root, and K(d 1 , d 2 ) − I has exactly one positive eigenvalue λ 1 , while other eigenvalues of K(d 1 , d 2 ) − I have negative real parts. Thus, the fixed-point index index(I − K(d 1 , d 2 ), (0, 0)) can be calculated by (18) when a 2 < µ1d1(1−µ1d2)
(1) a i = a j for any integer i = j.
By the above discussion, under assumptions of Theorem 6.1, we see that system (2)-(4) has a bifurcation solution emanating from positive constant solution. Using the Crandall-Rabinowitz bifurcation theorem, we know that there exist δ > 0 and a one-parameter family of smooth functions β : and denote U (s) = u(s) v(s) . Then (U (s); a i0 (s)) is the unique positive solution of (2)-(4) near U * i0 ; a i0 . Remark 11. If Ω is one dimensional, then all µ i are simple. Therefore, for each i, (U * i ; a i ) is always a bifurcation point to (2)-(4) provided that a = a i and a i = a j for i = j.
Remark 12. In Theorem 6.1, the assumption a i = a j for any positive integer i = j is only for the technical reason. In fact, if a i = a j for some positive integer i = j, then it is easy to see that dimN (G U (U * i ; a i )) > 1. The above local bifurcation result gives a precise description for the structure of positive solutions near the bifurcation points, it provides no information on the bifurcation curve far from the equilibrium. Therefore, a further study is necessary in order to understand the bifurcation curve when it is far away from the bifurcation point. In the following, we investigate the coexistence of the steady-state solution of (1) by considering the global bifurcation. The global bifurcation result shows that the bifurcation curve reaches to infinity. For simplicity, we suppose that Ω is one dimensional, say Ω = (0, 1). Note that the following discussion is mainly based on previous sections, especially Section 5.
In one dimensional case, a steady-state of (1) is a positive solution to the ordinary differential equations system Theorem 6.2. Suppose that Ω is one dimensional. If a = a i for some i and a i = a j for i = j, where a i is given as the same as in Theorem 6.1, then the global bifurcation curve Γ i of positive solution of system (27)- (29) which occurs at (U * i ; a i ) tends to infinity.
Proof. Consider the eigenvalue problem It is well known that all eigenvalues {µ i } ∞ i=0 of (30) are simple and µ i = i 2 π 2 , i = 0, 1, 2, · · · , whose corresponding normalized eigenfunctions are given by and all these eigenfunctions {ϕ i } ∞ i=0 form an orthonormal basis in L 2 ((0, 1)). As in Section 5, we rewrite the equations (27) and (28) in new forms such that the standard global bifurcation theory [25,28] can be conveniently applied. Let Then (27) and (28) where f and g are higher order terms of u and v (In fact, f and g have the same expressions as f 3 and f 4 defined in Section 5). Thus, the positive constant solution (27) and (28) transforms into U * = (0,0) of systems (32)-(33).
For w(x) ∈ C([0, 1]), assume that 1]) is unique, and the operators G d1 , G d2 are continuous and compact. Denote Then (32)-(33) can be interpreted as the equation where and H are compact on E, In order to apply the global bifurcation theorem, we first need to show that 1 is an eigenvalue of K(a) with odd algebraic multiplicity for a = a i .
In Section 5, we define the operator K(d 1 , d 2 ), it is easy to see that K(d 1 , d 2 ) = K(a) in one-dimensional case. By the proofs of Theorem 5.1 and Theorem 6.1, we then have where ϕ i is an eigenfunction of eigenvalue problem (30). This shows that 1 is an eigenvalue of K(a i ), and dimN (K(a i ) − I) = 1. Meanwhile, by the proof of Theorem 5.1 again, for n = 1, we also know that Therefore, 1 is indeed an eigenvalue of K(a i ) with algebraic multiplicity one. For small δ > 0 and each a i , denote by U (a i , δ) the δ-neighborhood of a i . Take 0 < a ∈ U (a i , δ), a = a i . Then Now we need to verify that the above index is different as a crosses a i , that is, ∀ε > 0, satisfying ε < δ, the following inequality holds.
Let µ be an eigenvalue of K(a) corresponding to eigenfunction ϕ ψ . Then we are the eigenfunctions of eigenvalue problem (30). Then we obtain Thus, all eigenvalues of K(a) exactly consist of all roots of the characteristic equation It is easy to check that if µ = 1 satisfies (36), then for a = a j , we must have a i = a j by the definition of a i , a j , and further, i = j by the assumption. This means that for a = a j , if i = j, then µ = 1. So, for i = j in (36), K(a) has the same number and multiplicities of eigenvalues larger than 1 for all a tending to a i if we do not count those eigenvalues corresponding to i = j in (36).
For i = j and a = a i , it can be easily seen that the two roots of equation (36) denoted by µ(a) and µ(a), respectively, have positive real parts, and For all a ∈ U (a i , δ), a = a i , we always have µ(a) < 1 since µ is continuous in a. Furthermore, we note that the constant term 2a 2 a 2 +d1µi of (36) is increasing in a, and then µ(a) is a decreasing function of a. This implies that Therefore, for each small ε > 0, K(a i − ε) has exactly one more eigenvalue that is larger than 1 than K(a i + ε) does. By a similar argument as in Section 5, we can still prove that such eigenvalue has algebraic multiplicity one, and (35) holds.
With the help of [25], we conclude that the bifurcation curve Γ i which emanates from (U * i ; a i ) either reaches to infinity or meets some other bifurcation point (U * j ; a j ), where i = j, a j > 0 (In the latter case, we may think that Γ i is bounded). In the following, we claim that the first alternative must occur.
Assume that Γ i dose not reach to infinity. Then Γ i must meet another bifurcation point, say (U * k ; a k ), and can not meet other bifurcation point (U * l ; a l ), where l > k.

Consider the problem
If U is a solution of (37), then using U , we can construct a solution of (27)-(29) by a reflective and periodic extension. For example, let x n = n k , n = 0, 1, 2, · · · , k and Then x ∈ [0, 1] and U (x) is a solution of (27)- (29). Clearly, (U * k ; a k ) is also a bifurcation point of (37). Denote Γ k the bifurcation curve which emanates from (U * k ; a k ). Then by the same argument as the above, it is easy to show that Γ k either reaches to infinity or meets some other bifurcation point (U * k ; a k ), k > k. If the latter alternative occurs, then it shows that Γ i meets (U * k ; a k ) too. This is an obvious contradiction. Therefore, Γ k reaches to infinity, and then by the extension again we know that Γ i reaches to infinity too, which also means that the projection of Γ i on the a-axis is unbounded.
Remark 13. Our above arguments only rule out the possibility that Γ i meets some other bifurcation point without finally tending to infinity. As to if it is possible that Γ i can meet some other bifurcation point and then reaches to infinity, it is unknown. By the uniqueness of bifurcation curve which occurs at each bifurcation point (U * i ; a i ), we know that if Γ i meets some other bifurcation point (U * j ; a j ) and then reaches to infinity, then Γ i coincides with Γ j (which emanates from (U * j ; a j )) when a ≥ a j , that is, Γ j is a part of Γ i . Whereas, if Γ i meets (U * j ; a j ) but does not reach to infinity, then Γ i only exists when a ∈ [a i , a j ].
So far, we deal with the bifurcation solution of (2)-(4) by taking a, the initial concentration of the reaction precursor, as a bifurcation parameter. For fixed a, if we take d 1 or d 2 , the diffusion rate of the reactant or autocatalyst, as a parameter, then by the fixed point index theory, using the techniques in obtaining Theorem 3.1 in [27], we can get other bifurcation results (Theorems 6.3 and 6.4), the proofs of which will not be given, the readers can refer to [27]. Denote Given a and d 2 , let Then we have following bifurcation result. Theorem 6.3. Let d 1 > 0 be any number. Consider the equilibrium ((u, v); d 1 ) = ((u * , v * ); d 1 ) of (2)-(4).
is not a bifurcation point to (2)-(4), where ∅ denotes the empty set.
(ii) Suppose Similarly, for fixed a and d 1 , let Then the next result can be obtained.
is the same as that of in Theorem 6.3.

7.
Numerical simulations. In the previous sections, we show the existence of spatial patterns for our model (1) by using the stability analysis, topological degree theory and bifurcation theory, which partially indicates the richness of the dynamics for the model. In this section, we give some numerical examples to illustrate the dynamical behaviors and the spatiotemporal pattern formation for model (1) corresponding to some analytical results, especially Theorem 5.1 and Theorem 6.1. We consider Ω = (0, 1). In this case, µ 1 = π 2 , µ 2 = 4π 2 . It follows from Theorem 5.1 that system (2)-(4) has at least a non-constant positive solution if µ 1 < 1 d2 < µ 2 and a 2 < µ1d1(1−µ1d2)

1+µ1d2
. We choose parameters d 1 = 1, d 2 = 0.05, a = 1.3 and d 1 = 1, d 2 = 0.06, a = 1.1, respectively. Then (1) has a non-constant positive solution, these facts are shown in Figure 1 and Figure  2, respectively, where the initial conditions are small random perturbations of the homogeneous steady state. The last two numerical examples show that system (1) admits periodic solution when a ≈ 1. That is, the Hopf bifurcation may occur for system (1) near positive constant solution, which is consistent with Remark 4.
All these numerical examples demonstrate the richness of spatiotemporal patterns of (1) for certain parameters sets (d 1 , d 2 , a). The numerical simulations show the subtle dependence of the spatiotemporal patterns on the diffusion coefficients and the reaction precursor. Chemically, this can be interpreted as: if the reactant diffuses relatively faster than autocatalyst, then the spatiotemporal patterns are more likely to appear.
8. Discussion and conclusion. In this paper, we investigate an autocatalytic reaction-diffusion system which has been identified as one of typical nonlinear mechanisms in biochemical procedures. The main motivation is to propose the effects of the autocatalyst concentration and the reaction process on the stability and the u(x, t) v(x, t) coexistence of the reactant and autocatalyst for autocatalytic reaction-diffusion systems. We find that the concentration of the reaction precursor and the chemical process both have significant influence on the dynamics of the model. For example, the process of the reactions drastically changes the stability properties of constant solution with regard to the model, this may be reflected by the changes of the reactant concentration and the diffusion process of the reactant and autocatalyst. Roughly speaking, the reaction will fully succeed if the reactant concentration and the diffusion process of the reactant and autocatalyst inside the reactor are controlled reasonably, and the conversion from the reactant to autocatalyst will occur until the concentration of autocatalyst in the reactor reaches a high level, that is, the positive steady-state exists; Conversely, if the catalyst diffuses too fast or reactant diffuses too slow, then the concentration of autocatalyst may eventually drop to zero, that is, the positive steady-state disappears.
For system (2)-(4), further studies may include a number of extensions in future work such as more numerical simulations, bifurcation diagram, existence of periodic solutions, which would also be of great interests for this system. Here, for periodic solutions, we would like to say something a bit more. A periodic solution of an ODE system is also a (spatially homogeneous) solution of the corresponding PDE system. On periodic solutions, Hopf bifurcation should be concerned. A Hopf bifurcation analysis can be performed for both ODE and PDE systems at the same bifurcation point, but from the local uniqueness of periodic solutions near Hopf bifurcation point, only spatial homogeneous periodic solutions exist near bifurcation point. However, the stability of these periodic solutions with respect to PDE systems could be different from that of for ODE systems. First, an unstable periodic solution of an ODE system is also unstable for the corresponding PDE system; Second, an unstable equilibrium solution of a PDE system may be stable for the corresponding ODE system, and the nearby bifurcating periodic solutions through Hopf bifurcation are also unstable. The latter case illustrates the interaction of Hopf instability and Turing instability. In Section 3, we show that the unique positive constant solution is a stable equilibrium point for the ODE system (13) for high concentration of reactant, but it may become unstable for PDE system (2)-(4) through Turing instability. Moreover, for system (13), there does not exist Hopf bifurcation solution for low concentration of reactant even though the constant solution is unstable; However, for system (1), the numerical simulations show that there exists Hopf bifurcation solution indeed, which tells us that the only possible reason which causes the Hopf bifurcation to occur is diffusion.