Stability and bifurcation with singularity for a glycolysis model under no-flux boundary condition

In this paper, a glycolysis model subject to no-flux boundary condition is considered. First, by discussing the corresponding characteristic equation, the stability of constant steady state solution is discussed, and the Turing's instability is shown. Next, based on Lyapunov-Schmidt reduction method and singularity theory, the multiple stationary bifurcations with singularity are analyzed. In particular, under no-flux boundary condition we show the existence of nonconstant steady state solution bifurcating from a double zero eigenvalue, which is always excluded in most existing works. Also, the stability, bifurcation direction and multiplicity of the bifurcation steady state solutions are investigated by the singularity theory. Finally, the theoretical results are confirmed by numerical simulations. It is also shown that there is no Hopf bifurcation on basis of the condition \begin{document}$ (C) $\end{document} .

1. Introduction. Glycolysis, which occurs in the cytosol, is thought to be the archetype of a universal metabolic pathway for cellular energy requirement. The wide occurrence of glycolysis indicates that it is one of the most ancient known metabolic pathways and a common way of providing limited energy for the organism in living nature. However, its significance lies in that it can supply the energy not only with a rapid speed, but more importantly under oxygen free conditions such as strenuous exercise and high altitude hypoxia. On the other hand, it is the main way of providing energy for the tumor cells [1,20,35].
Glycolysis model turns out to be a classic and representative system in biochemical reaction. All glycolysis models are based on the same reaction scheme. The different models stem from the difference in the reaction mechanism for key enzyme. The first model of glycolysis was proposed by Higgins [15], whose reaction in the later stage of glycolysis reaction, whose corresponding dimensionless system subject to no-flux boundary condition in the one-dimension spatial domain is x ∈ (0, l), t > 0, u x = v x = 0, x = 0, l, t > 0, u(x, 0) = u 0 (x) ≥ 0, v(x, 0) = v 0 (x) ≥ 0, x ∈ (0, l). (1) Here, u and v represent the chemical concentrations, d 1 and d 2 are the diffusion coefficients, δ is the dimensionless input flux, and k is the dimensionless rate constant for the low activity state. For k = 0, the model is Sel'kov model studied in [8,11,13,19,23,25,33]. Throughout this paper we assume that all constants d 1 , d 2 , δ and k are positive and 0 < k < 1/8 for the range of rational experiment data. The reaction mechanism and the generic class of model (1) are given in [22,28,32]. Concerning this model for a two cell system without spatial diffusion, the spatially homogeneous periodic solution is obtained in [10], which shows the sequential structure of time in the later stage of biochemical reaction. In view of spatial diffusion, there are some existence and stability results of constant steady states (see [3,21]). Under the fixed Dirichlet boundary condition, the authors [36] discuss the stability of constant steady state solution and the existence of non-constant steady state solutions not only from a simple eigenvalue, but more difficultly from a double one, which are confirmed by numerical simulations. Based on the following condition (C), no Hopf bifurcation can occur for the model (1) by taking the bifurcation parameter d 1 . Then, we aim to continue the stationary bifurcation of (1) for no-flux boundary condition. However, the steady state bifurcation mostly focuses on the case of the simple eigenvalue, such as [4,14,16,17,18,34,37,38], but rather less seems to be known about the non-simple case.
The corresponding method for the simple bifurcation is the traditional theory, see [7,26]. Just as [4,14,16,17,18,34,37,38], the bifurcation analysis may be from either a simple or non-simple zero eigenvalues, but the latter case is excluded in the most existing works. Motivated by this, we are interested in the formation of steady state spatially inhomogeneous solutions bifurcating from the double zero eigenvalue. This results in the steady-state solution with mode interaction of cos πj l x and cos πm l x, which is investigated in Section 3 and confirmed by numerical simulations in Section 5. We also show that no steady state solutions bifurcating from the higher eigenvalue can exist for model (1).
The main contribution of this paper is to focus on the singularity bifurcation under the no-flux boundary condition, i.e. the bifurcation where one or more of the hypotheses of the traditional theory fail. In this case, the involved methods are Lyapunov-Schmidt reduction and singularity theory. The singularity theory as developed by Golubitsky and Schaeffer [12] allows one to classify bifurcation into equivalence classes. It offers an extremely useful approach to the bifurcation problem, whose principal advantage is that it adapts well to the singularity bifurcation, such as the double bifurcation. Here, it is worth noting that Lyapunov-Schmidt reduction and singularity theory are also applied to analyze the bifurcation direction, multiplicity and stability of the solutions bifurcating from the simple eigenvalue.
The rest of the present paper is organized as follows. In Section 2, we first briefly analyze the Turing's instability of constant steady state solution, and then study bifurcation direction and multiplicity of the steady state solutions of (1) from simple bifurcation. In Section 3, we mainly discuss the bifurcation from a double eigenvalue for the no-flux boundary condition. In this case, the classical Crandall-Rabinowitz theorem can no longer be applied, and we investigate the conditions in detail for the bifurcation structure by use of Lyapunov-Schmidt procedure and singularity theory. In Section 4, we discuss the stability of the solutions from simple and double bifurcations by the singularity theory. Finally, in Section 5, numerical simulations are performed to confirm the analytic results.

Steady state solutions when
m for any m = j. Obviously, (u * , v * ) = ( δ k + δ 2 , δ) is the unique constant steady state solution of the PDE system (1) and the corresponding ODE system of (1). For the bifurcation from the constant solution, it is essential to analyze the Turing instability of the constant solution.
It is well known that the eigenvalue problem has eigenvalues λ i = (πi/l) 2 , i = 0, 1, 2, · · · with corresponding normalized eigenfunctions i > 0. Throughout this paper, we always assume where is defined to be the largest positive integer such that λ i < g 1 d 2 , and then denote d 1min = min{d Here, the linear part is and the nonlinear part is where Then the stationary solutions of (1) are corresponding to the solutions of the elliptic problem F (w, λ) = 0, w = (u, v), x ∈ (0, l) (4) with the boundary condition Thus we have F (0, λ) = 0 and take λ instead of d 1 as bifurcation parameter for the later discussions.
Based on the condition (C), it follows that (i) no Hopf bifurcation occurs for taking λ (or d 1 ) as the bifurcation parameter.
(ii) the model (1) is an activator-substrate system where the activator u consume the substrate v.
The following lemma shows that the Turing instability is found for the glycolysis model (1).
According to Eq. (2), it is readily found that at most two elements in {d 1 for j ∈ [1, Λ] to be the main discussed mode and divide our discussions into two cases: (i) d for some m = j in the next section.
Denote the inner product of Y by and the linearized operator by and X 1 = X ∩ R(L 0 ). We further define the operator P on Y by Then R(P ) = N (L 0 ), and it is easy to verify that P 2 = P which implies that P is the projection onto N (L 0 ). Thus, Q = I − P is a projection onto R(L 0 ) in Y . Hence, the system (4) can be transformed into According to the decomposition of X, we rewrite w ∈ X as w = sΦ j + W , where s ∈ R and W ∈ X 1 . It follows from the implicit theorem that the second equation of (5) is uniquely solvable for W in the neighborhood of zero. Let W (s, λ) := W (sΦ j , λ) denotes this solution for all small s and λ. It is obvious that W s (0, 0) = 0, W (0, λ) ≡ 0 for small λ, and then W λ (0, 0) = W λλ (0, 0) = · · · = 0. Substituting W (s, λ) into the first equation of (5), we can get Rewrite where H(w, λ) := λM w + N (w). (8) Therefore, the solutions of (4) are further determined by the zeros of the reduced equation (7).
In order that the given reduction problem (that is, the algebraic equation which results from the original partial differential equation by the Lyapunov-Schmidt reduction mehtod) is equivalent to such a normal form, certain conditions have to be met by the singularity theory.
It is obvious that h s (0, 0) = 0 from W s (0, 0) = 0 and h λ (0, 0) = h λλ (0, 0) = · · · = 0 from H(0, λ) = 0. To begin the following discussions, by straightforward calculations, we further show that the second and third derivatives of H at the origin are given as From (3) we further get the second and third Fréchet derivatives of N are separately given by and According to the second equation of (5), the second order derivative of W with respect to s at the origin denoted by W ss meets where Thus we obtain Hence, from (9)- (11) and (13), it is easy to get .

Steady state solutions when
for some m = j. From (2), one may verify that d which is the basic assumption in this section. It leads to d . Without loss of generality, we assume j < m.
In {d 1 (m = j) has been discussed in the previous section, and the case d for some m = j, then 0 is not a simple eigenvalue of L 0 and the Crandall-Rabinowitz theorem is no longer applied. In this case, 0 is a double eigenvalue of L 0 , and the Lyapunov-Schmidt reduction technique and singularity theory will be used to investigate whether or not the bifurcation occurs.
For d where (17) and X 1 = X ∩ R(L 0 ). By letting P be the projection of Y onto N (L 0 ), we must take the form of P as (4) is determined by a pair of equations P F (w, λ) = 0, QF (w, λ) = 0.
The second derivatives of H at the origin are given as with s j = s, s m = τ . From (10) and (21), we have The following vectors defined at the origin play a central role in the computations of Taylor coefficients and are exhibited here.
∂s i ∂τ j ∂λ k , and the second and third Fréchet derivatives of N are given by (10) and (11), respectively.
In this case, we have In the same way, the derivative W ss meets When m = 2j, we have When m = 2j, we know According to But it is worth noting that |B 2j | = 0, so we cannot solve c 2j d 2j as above. Based on the equation we solve algebraically for c 2j d 2j to get infinitely many solutions for any number k. Then it follows that

MEIHUA WEI, YANLING LI AND XI WEI
However, we note that W ss ∈ X 1 is an important point to be considered here, and W ss subtracts off an appropriate multiple of Φ j and/or Φ m to be the solution to Then we can find the unique solution W ss in X 1 . It is obvious that the first term of W ss is in X 1 , and then we only need the second term, denoted by W 2 , belongs to X 1 . We know that Thus, when m = 2j, we get Therefore, by the calculations, we obtain (l 2j,j + 2), m = 2j, . Then we have . Similarly, the second order derivative of W with respect to s and τ at the origin denoted by W sτ meets with e = e j + e m 2 . In the same way as above, for m = 2j, we have we can obtain to get infinitely many solutions for any number k. Since W sτ ∈ X 1 , for m = 2j, by the same method way as above we get the only solution in X 1 as Hence, we have , a 2 010 = Φ * j , , ν m+j , whereê = eδ =ê j +ê m 2 . Then we have a 010 = a 1 010 + a 2 010 + a 3 Similarly, we conclude Thus, we have Theorem 3.2. If A 1 C 0 p 1 p 2 p 3 = 0 for m = 2j andÃ 1 C 0p1p2p3 = 0 for m = 2j , then the reduced problem (20) is equivalent to the normal form s(ε 1 s 2 + ρτ 2 + λ + µ 1 s m−2 τ j ) τ (ε 2 τ 2 + κs 2 + λ + µ 2 s m τ j−2 ) with the ε i and ρ, κ given by The parameters µ 1 , µ 2 are determined as follows: By the condition (C), we have P i (d (29), and then L 0 has a positive eigenvalue. Thus, according to the perturbation theory of linearized operator, the bifurcation solutions of simple and double bifurcation described by (15) and (27) are unstable for j = m 1 . The proof is completed. Proof. For j = m 1 , it follows that by the Fredholm alternative, and then 0 a simple eigenvalue of L 0 . From (28) and (29), we have that for all i, ) > 0, i = 0, 1, 2, · · · , m 1 − 1, m 1 + 1, · · · .
Hence, 0 is a simple eigenvalue of L 0 with the largest real part, and all the other eigenvalues of L 0 lie in the left half complex plane. Thus we complete the proof. Therefore, from (14), we have the following stability result for the solution (15) of simple bifurcation for j = m 1 on the basis of Lemma 4.2 and [12], where the stability theorem in [31] is not valid. for any integer m = j. If h sss (0, 0) < 0 (> 0), then the bifurcation solution (15) is stable (unstable) for both s < 0 and s > 0.
5. Numerical results. The goal of this section is to present the numerical simulations for complementing the analytic results in the previous sections. We perform the initial-boundary-value problem (1) numerically by use of a standard implicit method, that is, the Crank-Nicholson scheme. Here, we transform the spatial domain from 0 < x < l to 0 <x < 1 by puttingx = x/l, and still denotex by x. The numerical results presented below are plotted at sufficiently large times so that the system reaches a steady state.