1. Introduction

Many differential equations have been proposed (see [8,11,13], [17]-[19], [21]-[22], [24,27] and references therein) to model the dynamic changes of substrate concentration and product one in enzyme-catalyzed reactions. Among those models, a typical form ([7]) is the following skeletal system

$ \left\{ $$\begin{array}{ll} \dot{x}=v-V_1(x, y)-V_{3}(x), \\ \dot{y}=q(V_{1}(x, y)-V_{2}(y)), \end{array}$$ \right. \label{sp0} $

(1)

where $x$ and $y$ denote the concentrations of the substrate and the product respectively, $v$ and $q$ are both positive constants, $V_1(x, y)$ and $V_2(y)$ denote the enzyme reaction rate and the output rate of the product respectively and satisfy that

$ V_{1}(0, y)=0, ~~ \partial V_{1}/\partial x>0, ~~ \partial V_{1}/\partial y>0, ~~ V_2(y)\geq0, ~~ \forall x, y>0, $

and $V_3(x)$ denotes the branched-enzyme reaction rate. Figure 1 shows the scheme of the enzyme-catalyzed reaction which comprises a branched network from the substrate. In Figure 1, $S$ and $P$ represent the substrate and product, respectively, and $E_1, E_2$ and $E_3$ are the three enzymes.

The case that $V_3(x)\equiv0$ in system (1), which represents an unbranched reaction, has been discussed extensively in [1,6,7,9,20]. Recently, more efforts were made to the case that $V_3(x)\not\equiv0$. One of the efforts ([12,23]) is made for $V_1(x, y)=x^my^n, V_2(y)=y$ and $V_3(x)=lx$ and $v=1$, with which system (1) reduces to

$ \left\{ $$\begin{array}{ll} \dot{x}=1-x^my^n-lx, \\ \dot{y}=q(x^my^n-y), \end{array}$$ \right. $

called the multi-molecular reaction model sometimes, where $m, n\geq1$ are integers and $l\geq0$ is real. All local bifurcations of this system such as saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation were discussed in [12] and [23]. Reference [15] is concerned with the case that $V_1(x, y)=\gamma x^my^n$, $V_3(x)=\beta x$, $q=1$ and $V_2(y)$ is a saturated reaction rate, i.e., $V_2(y)=v_2y/(\mu_2+y)$, with which (1) reduces to

$ \left\{ $$\begin{array}{ll} \dot{x}=v-\gamma x^my^n-\beta x, \\ \dot{y}=\gamma x^my^n-\frac{v_2y}{\mu_2+y}, \end{array}$$ \right. $

where $v, \gamma>0, {\mu}_{2}, v_2$ and $\beta\geq0$. Results on existence and nonexistence of periodic solutions on Hopf bifurcation were obtained in [15] with $n=1$ and $\beta=0$. When $V_2(y)$ and $V_3(x)$ are both saturated reaction rates, system (1) was considered in [16] as

$ \left\{ $$\begin{array}{ll} \dot{x}=v-V_1(x, y)-\frac{v_3 x}{u_3+x}, \\ \dot{y}=q(V_1(x, y)-\frac{v_2 y}{u_2+y}) \end{array}$$ \right. $

with $V_1(x, y)=v_1x(1+x)(1+y)^2/[L+(1+x)^2(1+y)^2]$, where $L$ is the allosteric constant of $E_1$. Varying the parameter $v_2$ but fixing the other parameters, Liu ([16]) investigated numerically how the enzyme saturation affects the emergence of dynamical behaviors such as the change from a stable oscillatory state to a divergent state. Later, Davidson and Liu ([3]) discussed the saddle-node bifurcation, Hopf bifurcation and the global bifurcation corresponding to the appearance of homoclinic orbit. When $V_2(y)$ and $V_3(x)$ are both saturated reaction rates, system (1) was also considered in [4] as

$ \label{sp1} \left\{ $$\begin{array}{ll} \dot{x}=v-v_{1}x y-\frac{v_{3}x}{u_{3}+x}, \\ \dot{y}=q(v_{1}x y-\frac{v_{2}y}{u_{2}+y}) \end{array}$$ \right. $

(2)

with $V_1(x, y)=v_{1}xy$. With a change of variables $x=u_{3} \tilde{x}, \ y=u_{2}\tilde{y}$ and the time rescaling $t\to v_{1}^{-1}\mu_{2}^{-1}t$, system (2) can be written as

$ \left\{ $$\begin{array}{ll} \dot{x}=a-x y-\frac{bx}{1+x}, \\ \dot{y}=\kappa y(x -\frac{c}{1+y}), \end{array}$$ \right. \label{xy0} $

(3)

where we still use $x$, $y$ to present $\tilde{x}$, $\tilde{y}$ and take notations $a:=v_{1}^{-1}u_{3}^{-1}u_{2}^{-1}v$, $b:=v_{1}^{-1}u_{3}^{-1}u_{2}^{-1}v_{3}$, $c:=v_{1}^{-1}u_{3}^{-1}u_{2}^{-1}v_{2}$ and $\kappa:=u_{2}^{-1}q u_{3}$ for positive constants. Actually, system (3) is orbitally equivalent to the following quartic polynomial differential system

$ \left\{ $$\begin{array}{ll} \dot{x}=(1+y)\{(1+x)(a-x y)-bx\}, \\ \dot{y}=\kappa(1+x)y\{(1+y)x -c\}, \end{array}$$ \right. \label{xy1} $

(4)

in the first quadrant ${\mathcal Q}_+:=\{(x, y):x\ge0, y\ge0\}$ by a time scaling $d\tau=(x+1)(y+1)dt$. In [4] Davidson, Xu and Liu discussed the case that $k=1$ and $a<c$, where the system has at most two equilibria, giving the existence of limit cycles (by the Poincaré-Bendixson Theorem seen in [10] or [26]) and the non-existence of periodic orbits (by the Dulac Criterion seen in [10] or [26]), proving the uniqueness of limit cycles (by reducing to the form of Liénard system) with some restrictions, and illustrating with the software AUTO saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation for fixed $\kappa=1, b=1.5$ and $c=5$. Recently, the general case that $\kappa, a, b, c>0$ was discussed in [27], where all codimension-one bifurcations such as saddle-node, transcritical and pitchfork bifurcations were investigated and the weak focus was proved to be of at most order $2$.

In this paper we continue the work of [27] to give conditions for the existence of a cusp and compute the parameter curves for the Bogdanov-Takens bifurcation, which induces the appearance of homoclinic orbits and periodic orbits, indicating the tendency to steady-states or a rise of periodic oscillations for the concentrations of the substrate and product.

---

2. Condition for cusp

It is proved in [27] that system (4) has at most 3 equilibria, i.e., $E_0: (a/(b-a), 0)$, $E_1:(p_1, c/p_1-1)$ and $E_2:(p_2, c/p_2-1)$, where

$ $$\begin{array}{c} p_1:=-\frac{1}{2}\big\{(a-b-c+1)-[(a-b-c+1)^2-4(a-c)]^{1/2}\big\}, \\ p_2:=-\frac{1}{2}\big\{(a-b-c+1)+[(a-b-c+1)^2-4(a-c)]^{1/2}\big\}. \end{array}$$ \label{p12} $

(5)

Moreover, if $a=a_*:=c+(b^{1/2}-1)^2$, then $E_1$ and $E_2$ coincide into one, i.e., the equilibrium $E_*: (b^{1/2}-1, c(b^{1/2}+1)/(b-1)-1)$. There are found in [27] totally 6 bifurcation surfaces

$ $$\begin{eqnarray*} {\mathcal T}_{E_0} &:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), b\ne(c+1)^2\}:=\bigcup\limits_{i=1}^{4}{\mathcal T}_{E_0}^{(i)}, \\ {\mathcal P}_{E_0} &:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), b=(c+1)^2\}, \\ {\mathcal H}_{E_1} &:=& \{(a, b, c, \kappa)\in\mathbb{R}_+^4|\kappa=\kappa_1, bc/(1+c)<a<c, 0<b\leq 1\} \nonumber\\ && \cup \{(a, b, c, \kappa)\in\mathbb{R}_+^4| \kappa=\kappa_1, bc/(1+c)<a<c + (b^{1/2} - 1)^2, 1 < b < (c + 1)^2\}, \\ {\mathcal SN}_{E_*} &:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~1<b<(c+1)^2, \kappa\neq \kappa_*\}:=\bigcup\limits_{i=1}^{4}{\mathcal SN}_{E_*}^{(i)}, \\ {\mathcal B}_1 &:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=c\}, \nonumber\\ {\mathcal B}_2 &:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=b\}, \end{eqnarray*}$$ $

which divide $\mathbb{R}_+^4:=\{(a, b, c, \kappa): a>0, b>0, c>0, \kappa>0\}$ into 8 subregions

$ $$\begin{eqnarray*} {\mathcal R}_1 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|c < a< a_*, 1<b<c, c>1, \\ && ~~{\rm or} ~~b < a< a_*, c < b < (c + 1)^2/4, c>1\}, \\ {\mathcal R}_2 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|b < a< c, 0<b<c\} \\ {\mathcal R}_3 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|bc/(1+c) < a< b, 0 < b < c ~~{\rm or} ~~bc/(1+c) < a< c, c < b < c + 1\}, \\ {\mathcal R}_4 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|0 < a < bc/(1+c), 0 < b < c + 1~~{\rm or} ~~0 < a < c, b > c + 1\}, \\ {\mathcal R}_5 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|c < a< bc/(1+c), b > c + 1\}, \\ {\mathcal R}_6 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4|c < a< b, c< b < (c + 1), c > 3 \\&& ~~{\rm or} ~~bc/(1 + c) < a< b, c + 1< b < (c + 1)^2/4, c > 3 \\&& ~~{\rm or} ~~bc/(1 + c) < a< a_*, (c + 1)^2/4 < b < (c + 1)^2, c > 3 \\&& ~~{\rm or} ~~c < a< b, c < b < (c + 1)^2/4, 1 < c \leq 3 \\&& ~~{\rm or} ~~c < a< a_*, (c + 1)^2/4 < b < c + 1, 1 < c \leq 3 \\&& ~~{\rm or} ~~bc/(1+c) < a< c+ (b^{1/2}-1)^2, (c + 1) < b < (c + 1)^2, c \leq 3 \\&& ~~{\rm or} ~~c < a< a_*, 1 < b < c + 1, c \leq 1 \}, \\ {\mathcal R}_7 &:=& \{(a, b, c, \kappa) \in \mathbb{R}^4_+|c+ (b^{1/2}-1)^2 < a < b, (c + 1)^2/4 < b < (c + 1)^2, c > 1 \\&& ~~{\rm or} ~~bc/(1 + c) < a< b, b > (c + 1)^2 ~~{\rm or} ~~ c < a< b, c < b < 1, c \leq 1 \\&& ~~{\rm or} ~~c+ (b^{1/2}-1)^2 < a< b, 1 < b < (c + 1)^2, c \leq 1\}, \\ {\mathcal R}_0 &:=& \mathbb{R}^4_+\backslash\big\{ {\mathcal P}_{E_0}\cup{\mathcal SN}_{E_*}\cup{\mathcal T}_{E_0}\cup \big(\bigcup\limits_{i=1}^2{\mathcal B}_i\big)\cup{\mathcal B}\cup\big(\bigcup\limits_{i=1}^7{\mathcal R}_i\big)\big\}, \end{eqnarray*}$$ $

where

$ $$\begin{array}{ll} {\mathcal T}_{E_0}^{(1)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), 0<b<c+1\}, \\ {\mathcal T}_{E_0}^{(2)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), c+1<b<(c+1)^2\}, \\ {\mathcal T}_{E_0}^{(3)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), b>(c+1)^2\}, \\ {\mathcal T}_{E_0}^{(4)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=bc/(1+c), b=c+1\}, \\ {\mathcal SN}_{E_*}^{(1)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~1<b<(c+1)^2/4, c>1, \kappa\neq \kappa_*\}, \\ {\mathcal SN}_{E_*}^{(2)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~b=(c+1)^2/4, c>1, \kappa\neq \kappa_*\}, \\ {\mathcal SN}_{E_*}^{(3)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~(c+1)^2/4<b<(c+1)^2, c>1, \kappa\neq \kappa_*\}, \\ {\mathcal SN}_{E_*}^{(4)}:=& \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~1<b<(c+1)^2, c\le1, \kappa\neq \kappa_*\}, \\ \kappa_1:=& p_1^{-2}\{(p_1+1)(c-p_1)\}^{-1}c\{p_1(c-p_1)+a\}, \\ \kappa_*:=& (c-b^{1/2}+1)^{-1}(b^{1/2}-1)^{-2}c^2. \end{array}$$ \label{pa0k0} $

(6)

The following lemma is a summary of Theorems 1, 2 and 3 of [27].

Lemma 2.1. (ⅰ) System (4) has a saddle-node $E_0$ if $(a, b, c, \kappa)\in {\mathcal T}_{E_0}\cup{\mathcal P}_{E_0}$. Moreover, as $(a, b, c, \kappa)$ crosses either ${\mathcal T}_{E_0}^{(1)}$ from ${\mathcal R}_3$ to ${\mathcal R}_4$, ${\mathcal T}_{E_0}^{(2)}$ from ${\mathcal R}_6$ to ${\mathcal R}_5$, or ${\mathcal T}_{E_0}^{(4)}$ from ${\mathcal R}_6$ to ${\mathcal R}_4$, a saddle $E_0$ and a stable (resp. unstable) node $E_1$ merge into a stable node $E_0$ on the boundary of the first quadrant ${\mathcal Q}_+$ for $\kappa<\kappa_1$(resp. $\kappa>\kappa_1$) through a transcritical bifurcation; as $(a, b, c, \kappa)$ crosses ${\mathcal T}_{E_0}^{(3)}$ from ${\mathcal R}_5$ to ${\mathcal R}_7$, a stable node $E_0$ and a saddle $E_2$ merge into a saddle $E_0$ on the boundary of ${\mathcal Q}_+$ through a transcritical bifurcation; as $(a, b, c, \kappa)$ crosses ${\mathcal P}_{E_0}$ from ${\mathcal R}_7$ to ${\mathcal R}_5$, a saddle $E_0$ changes into a stable node $E_0$, a saddle $E_2$ through a pitchfork bifurcation at $E_0$ on the boundary of ${\mathcal Q}_+$.

(ⅱ) System (4) has a weak focus $E_1$ of at most order 2 for $(a, b, c, \kappa)\in {\mathcal H}_{E_1}$, which is of order $\ell$ exactly and produces at most $\ell$ limit cycles through Hopf bifurcations as $(a, b, c, \kappa)\in {\mathcal H}_{E_1}^{(\ell)}$, $\ell=1, 2$, where ${\mathcal H}_{E_1}^{(1)}:={\mathcal H}_{E_1}\backslash{\mathcal H}_{E_1}^{(2)}$ and

$ $$\begin{eqnarray*} {\mathcal H}_{E_1}^{(2)} &:=& \{(a, b, c, \kappa)\in {\mathcal H}_{E_1}: 2p_1(p_1+1)a^3+\{(p_1^2+p_1+1)c^2+p_1(2p_1^2+p_1-2)c \\ && -3p_1^3(p_1+1)\}a^2-(c-p_1)\{(p_1^3+3p_1^2+p_1+1)c^2+2p_1^2(p_1^2+3p_1+3)c \\ && +3p_1^4(p_1+1)\}a+p_1^2\{(p_1+2)c+p_1^2\}\{c-p_1(p_1+1)\}(c-p_1)^2=0\}. \end{eqnarray*}$$ $

(ⅲ) System (4) has a saddle-node $E_*$ if $(a, b, c, \kappa)\in {\mathcal SN}_{E_*}$. Moreover, as $(a, b, c, \kappa)$ crosses either ${\mathcal SN}_{E_*}^{(1)}$ from ${\mathcal R}_0$ to ${\mathcal R}_1$, ${\mathcal SN}_{E_*}^{(2)}$ from ${\mathcal R}_0$ to ${\mathcal R}_6$, or ${\mathcal SN}_{E_*}^{(3)}\cup{\mathcal SN}_{E_*}^{(4)}$ from ${\mathcal R}_7$ to ${\mathcal R}_6$, a stable (resp. unstable) node $E_1$ and a saddle $E_2$ arise through a saddle-node bifurcation for $\kappa<\kappa_1$(resp. $\kappa>\kappa_1$).

The above Lemma 2.1 does not consider parameters in the set

$ \mathcal{B} := \{(a, b, c, \kappa)\in \mathbb{R}_+^4| a=a_*, ~1<b<(c+1)^2, \kappa=\kappa_*\}, \label{BT} $

(7)

where $a_*$ is given below (5) and $\kappa_*$ is given in (6). $\mathcal{B}$ is actually the intersection of the saddle-node bifurcation surface ${\mathcal SN}_{E_*}$ and the Hopf bifurcation surface ${\mathcal H}_{E_1}$, which are described by the curves $\widehat{{\mathcal SN}_{E_*}}$ and $\widehat{{\mathcal H}_{E_1}}$ respectively on the section $\{(a, b, c, \kappa)\in \mathbb{R}_+^4 | b=2, ~c=1\}$ in Figure 2. The intersection of $\widehat{{\mathcal SN}_{E_*}}$ and $\widehat{{\mathcal H}_{E_1}}$ indicates ${\mathcal B}$.

This paper is devoted to bifurcations in $\mathcal{B}$. For $(a, b, c, \kappa)\in \mathcal{B}$, equilibrium $E_*$ is degenerate with two zero eigenvalues. In the following lemma we prove that $E_*$ is a cusp.

Lemma 2.2. If $(a, b, c, \kappa)\in \mathcal{B}\backslash\mathcal{C}$, where

$ \mathcal{C}:=\big\{(a, b, c, \kappa)\in \mathcal{B} | c=\varsigma(b):=\frac{1}{4b^{1/2}}(b^{1/2}-1)\{b^{1/2}+2+(17b-12b^{1/2}+4)^{1/2}\}\big\}, $

then equilibrium $E_*$ is a cusp in system (4).

Proof. For simplicity in statements, we use the notation

$ p:=b^{1/2}-1. \label{p} $

(8)

For $(a, b, c, \kappa)\in{\mathcal B}$, system (4) can be transformed into the form

$ \left\{ $$\begin{array}{ll} \dot{x}=&y+\frac{c(p^2+cp+c)}{p^3}x^2+\frac{1}{p+1}xy-\frac{p}{c^2(p+1)}y^2 -\frac{c(p^2+c)}{p^4}x^3-\frac{p^2+2pc+2c}{p^2c(p+1)}x^2y \\& -\frac{2p+1}{c^2(p+1)^2}xy^2-\frac{c^2}{p^4}x^4-\frac{2}{p^2(p+1)}x^3y-\frac{1}{c^2(p+1)^2}x^2y^2, \\ \dot{y}=&-\frac{c^3(p+1)}{p^3}x^2-\frac{c^2(p+1)}{p^2(c-p)}xy-\frac{1}{c-p}y^2-\frac{(p+1)(p^2+c)}{p^5(c-p)}x^3 -\frac{c(p^2+2pc+2c)}{p^3(c-p)}x^2y \\& -\frac{2p+1}{p(p+1)(c-p)}xy^2-\frac{c^4(p+1)}{p^5(c-p)}x^4-\frac{2c^2}{p^3(c-p)}x^3y -\frac{1}{p(p+1)(c-p)}x^2y^2, \end{array}$$ \right. \label{bt0} $

(9)

by translating $E_*$ to the origin $O$ and Jordanizing the linear part of system (4). For convenience, introducing new variables $(x, y)\mapsto (u, v)$, where $u=x$ and $v$ denotes the right-hand side of the first equation in (9), we change (9) into the Kukles form

$ \left\{ $$\begin{array}{ll} \dot{u}=&v, \\ \dot{v}=&-\frac{c^3(p+1)}{p^3}u^2+\frac{c\{(2p+2)c^2-(p^2+3p)c-2p^3\}}{p^3(c-p)}uv+\frac{c-2p-1}{(p+1)(c-p)}v^2 +\frac{c^3(p^2+c)}{p^4(c-p)}u^3 \\& -\frac{c\{(p+1)(p+3)c^2+p(p^2-3p-3)c-p^3(3p+2)\}}{p^4(p+1)(c-p)}u^2v -\frac{(5p^2+8p+4)c+2p^2(p+1)}{cp^2(p+1)^2}uv^2 \\&-\frac{1}{c^2(p+1)}v^3 -\frac{c^2(c^2+2p^2c-p^3)}{p^5(c-p)}u^4 +\frac{1}{p^5(p+1)^2(c-p)}\{(p+4)(p+1)^2c^3 \\& +p(7p^3+7p^2-3p-4)c^2-p^3(8p^2+15p+8)c-2p^5(p+1)\}u^3v \\& +\frac{(3p^3+6p^2+6p+2)c^2+p(2p+1)(2p^2+2p-1)c-p^3(p+1)(7p+4)}{cp^3(p+1)^3(c-p)}u^2v^2 \\& -\frac{(3p+4)c^2-3p(p+2)c-2p^3}{c^3p(p+1)^2(c-p)}uv^3 -\frac{2c-3p}{c^4(p+1)^2(c-p)}v^4+O(|u, v|^5). \end{array}$$ \right. \label{bt1} $

(10)

Since the linear part is nilpotent, by Theorem 8.4 in [14] system (10) is conjugated to the Bogdanov-Takens normal form, i.e., the right-hand side of the second equation is a sum of terms of the form $au^k+bu^{k-1}v$. Hence, one can use the transformation $u\rightarrow u$, $v\rightarrow v-\frac{c-2p-1}{(p+1)(c-p)} uv$ together with the time-rescaling $dt=(1-\frac{c-2p-1}{(p+1)(c-p)}u)d\tau$ to change system (10) into the following

$ \left\{ $$\begin{array}{ll} \dot{u}=&v, \\ \dot{v}=&-\frac{c^3(p+1)}{p^3}u^2+\frac{c\{(2p+2)c^2-(p^2+3p)c-2p^3\}}{p^3(c-p)}uv +O(|u, v|^3), \end{array}$$ \right. \label{bt2} $

(11)

where the term of $v^2$ is eliminated and terms of degree 2 are normalized. The term of $u^2$ exists since $-{c^3(p+1)}/{p^3}\neq 0$. For the existence of the term of $uv$, we need to discuss on the quadratic equation

$ c^2-\frac{p^2+3p}{2(p+1)}c-\frac{p^3}{p+1}=0, \label{eq2} $

(12)

which comes from the numerator of the coefficient of $uv$. Since the constant term is negative for $p>0$, the quadratic equation (12) has exactly one positive root

$ c=\frac{1}{4}(p+1)^{-1}p\{p+3+(17p^2+22p+9)^{1/2}\}, $

which defines the function $\varsigma(b)$ as shown in Lemma 2.2 with the replacement (8). It implies by Theorem 8.4 of [14] that for $c\neq\varsigma(b)$, i.e., $(a, b, c, \kappa)\in {\mathcal S}\backslash{\mathcal C}$, $O$ is a cusp of system (11). The proof of this lemma is completed.

---

3. Bogdanov-Takens bifurcation

In this section we discuss in the case that $(a, b, c, \kappa)\in {\mathcal B}\backslash{\mathcal C}$, in which system (4) is of codimension 2. We choose $a, \kappa$ as the bifurcation parameters and unfold the Bogdanov-Takens normal forms of codimensions $2$ when $(a, \kappa)$ is perturbed near the point $(a_*, \kappa_*)$, where $a_*$ is given below (5) and $\kappa_*$ is given in (6).

Theorem 3.1. If $(a, b, c, \kappa)\in \mathcal{B}\backslash{\mathcal C}$, where ${\mathcal B}$ is defined in (7) and ${\mathcal C}$ is defined as in Lemma 2.2, then there are a neighborhood $U$ of the point $(a_*, \kappa_{*})$ in the $(a, \kappa)$-parameter space and four curves

$ $$\begin{eqnarray*} \mathcal{SN}^{+} &:=& \Big\{(a, \kappa)\in U|a=a_*, ~\kappa>\kappa_*, 0<c<\varsigma(b)\Big\} \cup \Big\{(a, \kappa)\in U|a=a_*, ~\kappa<\kappa_*, c>\varsigma(b)\Big\}, \\ \mathcal{SN}^{-} &:=& \Big\{(a, \kappa)\in U|a=a_*, ~\kappa<\kappa_*, 0<c<\varsigma(b)\Big\} \cup \Big\{(a, \kappa)\in U|a=a_*, ~\kappa>\kappa_*, c>\varsigma(b)\Big\}, \\ \mathcal{H} &:=& \Big\{(a, \kappa)\in U|a=a_*- \Big((2b^{1/2}+1)c^2-((b^{1/2}-1)^2+3(b^{1/2}-1))c \\ &&~~~~~~~~~~~~~ -2(b^{1/2}-1)^3\Big)^{-2} b^{1/2}(b^{1/2}-1)^6(c-b^{1/2}+1)^4 (\kappa-\kappa_*)^2 +O(|\kappa-\kappa_*|^{3}), \\ && ~~~~~~~~~~~~~\kappa>\kappa_*, 0<c<\varsigma(b)\Big\} \\&&\cup \Big\{(a, \kappa)\in U|a=a_*- \Big((2b^{1/2}+1)c^2-((b^{1/2}-1)^2+3(b^{1/2}-1))c \\ &&~~~~~~~~~~~~~ -2(b^{1/2}-1)^3\Big)^{-2} b^{1/2}(b^{1/2}-1)^6(c-b^{1/2}+1)^4 (\kappa-\kappa_*)^2 +O(|\kappa-\kappa_*|^{3}), \\&& ~~~~~~~~~~~~~\kappa<\kappa_*, c>\varsigma(b)\Big\}, \\ \mathcal{L} &:=& \Big\{(a, \kappa)\in U|a=a_*-49/25 \Big((2b^{1/2}+1)c^2-((b^{1/2}-1)^2+3(b^{1/2}-1))c \\ &&~~~~~~~~~~~~~ -2(b^{1/2}-1)^3\Big)^{-2} b^{1/2}(b^{1/2}-1)^6(c-b^{1/2}+1)^4 (\kappa-\kappa_*)^2 +O(|\kappa-\kappa_*|^{3}), \\&& ~~~~~~~~~~~~~\kappa>\kappa_*, 0<c<\varsigma(b)\Big\} \\&&\cup \Big\{(a, \kappa)\in U|a=a_*-49/25 \Big((2b^{1/2}+1)c^2-((b^{1/2}-1)^2+3(b^{1/2}-1))c \\ &&~~~~~~~~~~~~~ -2(b^{1/2}-1)^3\Big)^{-2} b^{1/2}(b^{1/2}-1)^6(c-b^{1/2}+1)^4 (\kappa-\kappa_*)^2 +O(|\kappa-\kappa_*|^{3}), \\&& ~~~~~~~~~~~~~\kappa<\kappa_*, c>\varsigma(b)\Big\}, \end{eqnarray*}$$ $

such that system (4) produces a saddle-node bifurcation near $E_{*}$ as $(a, c)$ acrosses $\mathcal{SN}^{+}\cup \mathcal{SN}^{-}$, a Hopf bifurcation near $E_{*}$ as $(a, \kappa)$ acrosses $\mathcal{H}$, and a homoclinic bifurcation near $E_{*}$ as $(a, \kappa)$ acrosses $\mathcal{L}$, where $\kappa_*$ and $\varsigma(b)$ are given in (6) and Lemma 2.2 respectively.

The above bifurcation curve ${\mathcal H}$ is exactly the same as ${\mathcal H}_{E_1}$ given in Lemma 2.1, and the union $\mathcal{SN}^{+}\bigcup \mathcal{SN}^{+}$ is exactly the bifurcation curves $\mathcal{SN}_{E_*}$ given in Lemma 2.1.

Proof. Let $p=b^{1/2}-1$ and

$ \varepsilon_{1}:=a-a_*, \varepsilon_{2}:=\kappa-\kappa_*, \label{eee} $

(13)

and consider $|\varepsilon_{1}|$ and $|\varepsilon_{2}|$ both to be sufficiently small. Expanding system (4) at $E_{*}$, we get

$ \label{btb0} \left\{ $$\begin{array}{ll} \dot{x}=&\frac{c(p+1)}{p}\varepsilon_1+(\frac{-c^{2}(p+1)}{p^{2}}+\frac{c}{p}\varepsilon_1)x +(-c(p+1)+(p+1)\varepsilon_1)y\\ &-\frac{c(c-p)}{p^{2}}x^{2} +(-\frac{c(2+3p)}{p}+\varepsilon_1)xy-p(p+1)y^{2}+O(\|(x, y)\|^3), \\ \dot{y}=&(\frac{c^3(p+1)}{p^{4}}+\frac{c(p+1)(c-p)}{p^2}\varepsilon_2)x +(\frac{c^2(p+1)}{p^2}+(p+1)(c-p)\varepsilon_2)y\\ &+(\frac{c^3}{p^{4}}+\frac{c(c-p)}{p^2}\varepsilon_2)x^{2} +(\frac{c^3(2+3p)-c^2p(2p+1)}{(c-p)p^3}+\frac{c(3p+2)-p(2p+1)}{p}\varepsilon_2)xy\\ &+(\frac{c^2(p+1)}{(c-p)p}+p(p+1)\varepsilon_2)y^{2}+O(\|(x, y)\|^3). \end{array}$$ \right. $

(14)

Introducing new variables $(x, y)\mapsto (\xi_1, \eta_1)$, where $\xi_1=x$ and $\eta_1$ denotes the right-hand side of the first equation in (14), we change (14) into the Kukles form, whose second order truncation is the following

$ \left\{ $$\begin{array}{ll} \dot{\xi_{1}}=&\eta_{1}, \\ \dot{\eta_{1}}=&E_{00}(\varepsilon_{1}, \varepsilon_{2})+E_{10}(\varepsilon_{1}, \varepsilon_{2})\xi_{1} +E_{20}(\varepsilon_{1}, \varepsilon_{2})\xi_{1}^{2}\\ &+ F(\xi_{1}, \varepsilon_{1}, \varepsilon_{2}) \eta_{1} +E_{02}(\varepsilon_{1}, \varepsilon_{2})\eta_{1}^{2}, \end{array}$$ \right. $

(15)

where $ F(\xi_{1}, \varepsilon_{1}, \varepsilon_{2}):=E_{01}(\varepsilon_{1}, \varepsilon_{2}) +E_{11}(\varepsilon_{1}, \varepsilon_{2})\xi_{1} $ and $E_{ij}$ s ($i, j=0, 1, 2$) are given in Appendix. Notice that $(a, b, c, \kappa)\in {\mathcal B}\backslash{\mathcal C}$ implies that $c\neq\varsigma(b)$. From (12) we see that the quadratic equation has exactly one positive root $c=\varsigma(b)$. Thus, for $c\neq\varsigma(b)$ we can check that

$ F(0, 0, 0)=0, ~~~ \frac{\partial F}{\partial\xi_{1}}(0, 0, 0)=E_{11}(0, 0)=(2p+2)(c^2-\frac{p^2+3p}{2(p+1)}c-\frac{p^3}{p+1})\ne 0. $

By the Implicit Function Theorem, there exists a function $\xi_{1}=\xi_{1}(\varepsilon_{1}, \varepsilon_{2})$ defined in a small neighborhood of $(\varepsilon_{1}, \varepsilon_{2})=(0, 0)$ such that $\xi_{1}(0, 0)=0$ and $F(\xi_{1}(\varepsilon_{1}, \varepsilon_{2}), \varepsilon_{1}, \varepsilon_{2})=0$. Thus, from the definition of $F$ we obtain $\xi_{1}(\varepsilon_{1}, \varepsilon_{2})= -{E_{01}(\varepsilon_{1}, \varepsilon_{2})}/{E_{11}(\varepsilon_{1}, \varepsilon_{2})}$ near $(0, 0)$. Then, we use a parameter-dependent shift

$ \xi_{2}=\xi_{1}-\xi_{1}(\varepsilon_{1}, \varepsilon_{2}), \eta_{2}=\eta_{1} $

to vanish the term proportional to $\eta_2$ in the equation for $\eta_2$ from system (15), which leads to the following system

$ \label{btphi} \left\{ $$\begin{array}{ll} \dot{\xi_{2}}=&\eta_{2}, \\ \dot{\eta_{2}}=&\psi_{1}(\varepsilon_{1}, \varepsilon_{2})+\psi_{2}(\varepsilon_{1}, \varepsilon_{2})\xi_{2}+ E_{20}(\varepsilon_{1}, \varepsilon_{2})\xi_{2}^{2}+E_{11}(\varepsilon_{1}, \varepsilon_{2})\xi_{2}\eta_{2} +E_{02}(\varepsilon_{1}, \varepsilon_{2})\eta_{2}^{2}, \end{array}$$ \right. $

(16)

where

$ $$\begin{eqnarray*} \psi_{1}(\varepsilon_{1}, \varepsilon_{2})&:=&E_{00}(\varepsilon_{1}, \varepsilon_{2}) +E_{10}(\varepsilon_{1}, \varepsilon_{2})\xi_{1}(\varepsilon_{1}, \varepsilon_{2}) +E_{20}(\varepsilon_{1}, \varepsilon_{2})\xi_{1}^{2}(\varepsilon_{1}, \varepsilon_{2}), \\ \psi_{2}(\varepsilon_{1}, \varepsilon_{2})&:=&E_{10}(\varepsilon_{1}, \varepsilon_{2}) +2\xi_{1}(\varepsilon_{1}, \varepsilon_{2})E_{20}(\varepsilon_{1}, \varepsilon_{2}). \end{eqnarray*}$$ $

In order to eliminate the $\eta_2^2$ term, one can use the transformation

$ \xi_3= \xi_2, ~~~ \eta_3= \eta_2-E_{02}(\varepsilon_{1}, \varepsilon_{2})\xi_2\eta_2 $

together with the time-rescaling $dt=(1-E_{02}(\varepsilon_{1}, \varepsilon_{2})\xi_2)d\tau$ to change system (16) into the following

$ \label{btphi2} \left\{ $$\begin{array}{ll} \dot{\xi_{3}}=&\eta_{3}, \\ \dot{\eta_{3}}=&\zeta_{1}(\varepsilon_{1}, \varepsilon_{2}) +\zeta_{2}(\varepsilon_{1}, \varepsilon_{2})\xi_{3} +\tilde{E}_{20}(\varepsilon_{1}, \varepsilon_{2})\xi_{3}^{2} +E_{11}(\varepsilon_{1}, \varepsilon_{2})\xi_{3}\eta_{3}, \end{array}$$ \right. $

(17)

where

$ $$\begin{eqnarray*} \label{phi12} \zeta_{1}(\varepsilon_{1}, \varepsilon_{2})&:=& \psi_{1}(\varepsilon_{1}, \varepsilon_{2}), ~~~ \zeta_{2}(\varepsilon_{1}, \varepsilon_{2}):= \psi_{2}(\varepsilon_{1}, \varepsilon_{2})-\psi_{1}(\varepsilon_{1}, \varepsilon_{2})E_{02}(\varepsilon_{1}, \varepsilon_{2}), \\ \tilde{E}_{20}(\varepsilon_{1}, \varepsilon_{2})&:=&E_{20}(\varepsilon_{1}, \varepsilon_{2})-E_{10}(\varepsilon_{1}, \varepsilon_{2})E_{02}(\varepsilon_{1}, \varepsilon_{2}). \end{eqnarray*}$$ $

Further, in order to reduce coefficient of $\xi_3^2$ to $1$, we apply the transformation

$ u=\frac{\tilde{E}_{20}(\varepsilon_{1}, \varepsilon_{2})}{E_{11}^2(\varepsilon_{1}, \varepsilon_{2})}\xi_{3}, v={\rm sign}\Big(\frac{E_{11}(\varepsilon_{1}, \varepsilon_{2})}{\tilde{E}_{20}(\varepsilon_{1}, \varepsilon_{2})}\Big) \frac{\tilde{E}_{20}^2(\varepsilon_{1}, \varepsilon_{2})}{E_{11}^3(\varepsilon_{1}, \varepsilon_{2})}, $

where $ \tilde{E}_{20}(0, 0)=-\frac{c^3(p+1)}{p^3}<0, $ and the time-scaling $dt=|\frac{E_{11}(\varepsilon_{1}, \varepsilon_{2})}{\tilde{E}_{20}(\varepsilon_{1}, \varepsilon_{2})}|d\tau$ to system (17) and obtain

$ \label{btf} \left\{ $$\begin{array}{ll} \dot{u}=&v, \\ \dot{v}=&\phi_{1}(\varepsilon_{1}, \varepsilon_{2}) +\phi_{2}(\varepsilon_{1}, \varepsilon_{2})u+u^{2}+\vartheta uv, \end{array}$$ \right. $

(18)

where $\vartheta={\rm sign}\Big(\frac{E_{11}(0, 0)}{\tilde{E}_{20}(0, 0)}\Big)$,

$ $$\begin{eqnarray*} \phi_{1}(\varepsilon_{1}, \varepsilon_{2}) &:=& \frac{E_{11}^4(\varepsilon_{1}, \varepsilon_{2})}{\tilde{E}_{20}^3(\varepsilon_{1}, \varepsilon_{2})} \zeta_{1}(\varepsilon_{1}, \varepsilon_{2}) \nonumber\\ &=& \frac{\{(2p+2)c^2-(p^2+3p)c-2p^3\}^4\varepsilon_1\phi_{11}(\varepsilon_1, \varepsilon_2)}{p^4(c-p)^4\phi_{12}^2(\varepsilon_1, \varepsilon_2)}, ~ \\ \phi_{2}(\varepsilon_{1}, \varepsilon_{2}) &:=& \frac{E_{11}^2(\varepsilon_{1}, \varepsilon_{2})}{\tilde{E}_{20}^2(\varepsilon_{1}, \varepsilon_{2})} \zeta_{2}(\varepsilon_{1}, \varepsilon_{2}) \nonumber\\ &=& \frac{\sqrt{2}\{(2p+2)c^2-(p^2+3p)c-2p^3\}\phi_{21}(\varepsilon_1, \varepsilon_2)}{c^{3/2}(c-p)^2(p+1)^{1/2}p\phi_{12}^{3/2}(\varepsilon_1, \varepsilon_2)}, \nonumber \end{eqnarray*}$$ $

and polynomials $\phi_{ij}$ s are given in the Appendix.

Let

$ \label{muphi} \mu_{1}=\phi_{1}(\varepsilon_{1}, \varepsilon_{2}), \mu_{2}=\phi_{2}(\varepsilon_{1}, \varepsilon_{2}), $

(19)

where $\phi_{1}$ and $\phi_{2}$ are defined just below (18). Clearly, $\phi_{1}(0, 0)=\phi_{2}(0, 0)=0$. Compute the Jacobian determinant of (19) at the point $(0, 0)$

$ \left| $$\begin{array}{cc} \frac{\partial \phi_1(\varepsilon_{1}, \varepsilon_{2})}{\partial \varepsilon_1}& \frac{\partial \phi_1(\varepsilon_{1}, \varepsilon_{2})}{\partial \varepsilon_2} \\ \frac{\partial \phi_2(\varepsilon_{1}, \varepsilon_{2})}{\partial \varepsilon_1}&\frac{\partial \phi_2(\varepsilon_{1}, \varepsilon_{2})}{\partial \varepsilon_2} \end{array}$$ \right|_{(\varepsilon_{1}, \varepsilon_{2})=(0, 0)} = -\frac{\{(2p+2)c^2-(p^2+3p)c-2p^3\}^5}{p^6c^4(c-p)^4(p+1)}\neq0, \label{Jacob} $

(20)

implying that (19) is a locally invertible transformation of parameters. This transformation makes a local equivalence between system (18) and the versal unfolding system

$ \label{btw} \left\{ $$\begin{array}{ll} \dot{\widetilde{u}}=&\widetilde{v}, \\ \dot{\widetilde{v}}=&\mu_1+\mu_2\widetilde{u}+{\widetilde{u}}^{2}+ \vartheta\widetilde{u}\widetilde{v}, \end{array}$$ \right. $

(21)

where $\vartheta$ is given in (18). As indicated in Section 7.3 of [10], system (21) has the following bifurcation curves

$ \label{btbifur} $$\begin{array}{rl} \mathcal{SN}^{+}& :=\{(\mu_{1}, \mu_{2})\in V_{0}\ |\ \mu_{1}=0, \ \mu_{2}>0\}, \\ \mathcal{SN}^{-}& :=\{(\mu_{1}, \mu_{2})\in V_{0}\ |\ \mu_{1}=0, \ \mu_{2}<0\}, \\ \mathcal{H}& :=\{(\mu_{1}, \mu_{2})\in V_{0}\ |\ \mu_{1}=-\mu_{2}^{2}, \ \mu_{2}>0\}, \\ \mathcal{L}& :=\{(\mu_{1}, \mu_{2})\in V_{0}\ |\ \mu_{1}=-\frac{49}{25}\mu_{2}^{2} +o(|\mu_{2}|^{2}), \ \mu_{2}>0\}, \end{array}$$ $

(22)

where $V_{0}$ is a small neighborhood of $(0, 0)$ in $\mathbb{R}^{2}$.

In what follows, we present above bifurcation curves in parameters $\varepsilon_1$ and $\varepsilon_2$ in explicit forms. For this purpose, we need the relation between $(\varepsilon_1, \varepsilon_2)$ and $(\mu_1, \mu_2)$. Note that $\phi_1$ and $\phi_2$ defined just below (18) are $C^{k}$ near the origin (0, 0)($k$ is an arbitrary integer). By condition (20), the well-known Implicit Function Theorem implies that there are two $C^k$ functions

$ \varepsilon_1=\omega_1(\mu_1, \mu_2), ~~~~\varepsilon_2=\omega_2(\mu_1, \mu_2) \label{eps12} $

(23)

in a small neighborhood of $(0, 0, 0, 0)$ such that $\omega_1(0, 0)=\omega_2(0, 0)=0$ and

$ \mu_1=\phi_1(\omega_1(\mu_1, \mu_2), \omega_2(\mu_1, \mu_2)), ~~~ \mu_2=\phi_2(\omega_1(\mu_1, \mu_2), \omega_2(\mu_1, \mu_2)). \label{compare-coeff} $

(24)

Substitute the second order formal Taylor expansions of $\omega_1$ and $\omega_2$ in (24) while expand $\phi_1$ and $\phi_2$ in (24) to the second order

$ $$\begin{eqnarray*} \phi_1(\varepsilon_1, \varepsilon_2) &=& \{(2p+2)c^2-(p^2+3p)c-2p^3\}^4\varepsilon_1/ \{p^6c^2(c-p)^4(p+1)\} -\{(2p+2)c^2 \nonumber\\&& -(p^2+3p)c -2p^3\}^4(24p^2c^4+42c^4p+21c^4-8p^3c^3-54c^3p^2-44c^3p \nonumber\\&& -36c^2p^4-12p^3c^2 +27p^2c^2+8p^5c+32cp^4+16p^6)\varepsilon_1^2 /\{2c^4p^8(c-p)^6 \nonumber\\&& (p+1)^2\}-\{(2p+2)c^2 -(p^2+3p)c-2p^3\}^4\varepsilon_1\varepsilon_2/\{(c^4p^4(c-p)^3(p+1)\} \nonumber\\&& +o(|\varepsilon_1, \varepsilon_2|^2), \end{eqnarray*}$$ $

(25)

$ $$\begin{eqnarray*} \phi_2(\varepsilon_1, \varepsilon_2) &=& \{(2p+2)c^2-(p^2+3p)c-2p^3\}\varepsilon_1/\{2c^2(p^3-2cp+p^2+c^2p+c^2-2cp^2)p^4\} \nonumber\\&& -\{(2p +2)c^2 -(p^2+3p)c-2p^3\}\varepsilon_2/c^2-\{(2p+2)c^2-(p^2+3p)c-2p^3\} \nonumber\\&& (-243p^3c^3 +832p^3c^4+513p^2c^4+455p^4c^3-594p^5c^2-1347p^3c^5 -1209p^2c^5 \nonumber\\&& +165p^4c^4 +1138p^5c^3-324p^6c^2-424p^7c -200p^5c^4+382p^6c^3+512p^7c^2 \nonumber\\&& -520cp^8-396c^5p -48p^9+108c^6-48p^{10}+384c^6p^3+414c^6p-104cp^9 \nonumber\\&& +264c^2p^8+594c^6p^2-672c^5p^4 +96c^6p^4-136c^5p^5-44c^4p^6 -76c^3p^7)\varepsilon_1^2 \nonumber\\&& /\{4c^3(p+1)^2(c-p)^4p^6\}-\{(2p+2)c^2 -(p^2+3p)c-2p^3\}(8p^2c^4+23c^4p \nonumber\\&& +12c^4+30p^3c^3 +8c^3p^2-22c^3p-58c^2p^4 -85p^3c^2+6p^2c^2-8p^5c+46cp^4 \nonumber\\&& +24p^6)\varepsilon_1\varepsilon_2 /\{4c^4p^2(p+1)(c-p)^2\} +(c-p)p^2 \{(2p+2)c^2-(p^2+3p)c \nonumber\\&& -2p^3\}\varepsilon_2^2/c^4+o(|\varepsilon_1, \varepsilon_2|^2). \end{eqnarray*}$$ $

(26)

Then, comparing the coefficients of terms of the same degree in (24), we obtain the second order approximations

$ $$\begin{eqnarray*} \varepsilon_1 &=& c^2p^6(c-p)^4(p+1)\mu_1/\{(2p+2)c^2-(p^2+3p)c-2p^3\}^4 + c^2p^{10}(c-p)^6(p+1)(32p^2c^4 \nonumber\\&& +56c^4p +27c^4-16p^3c^3-79c^3p^2 -59c^3p-48c^2p^4-19p^3c^2+36p^2c^2+12p^5c \nonumber\\&& +50cp^4+24p^6)\mu_1^2 /\{2\{(2p+2)c^2-(p^2+3p)c-2p^3\}^8\} + c^2p^8(c-p)^5(p+1) \nonumber\\&& \mu_1\mu_2 /\{(2p+2)c^2 -(p^2+3p)c-2p^3\}^5 +o(|\mu_1, \mu_2|^2), \end{eqnarray*}$$ $

(27)

$ $$\begin{eqnarray*} \varepsilon_2 &=& c^2p^2(c-p)^2(-8p^5-12cp^4-18cp^3+8c^3p^2-11p^2c^2 -9c^2p +14c^3p+6c^3)\mu_1 \nonumber\\ && / \{2\{(2p+2)c^2 -(p^2+3p)c-2p^3\}^4\} -c^2\mu_2/\{(2p+2)c^2-(p^2+3p)c-2p^3\} \nonumber\\ && +c^2p^6(c-p)^4(1314c^7p^2 +630p c^7-270p^3c^4+2068p^3c^5 +612p^2c^5 +677p^4c^4 \nonumber\\ && -1134p^5c^3+4387p^5c^4 -1056p^6c^3 -1804p^7c^2-3741c^6p^3+756c^5p^4 +1160c^3p^8 \nonumber\\ && -2268c^6p^4+1176c^7p^3-1272c^5p^6 -352c^6p^5 +384c^7p^4-320p^{11}+108c^7-704cp^{10} \nonumber\\ && +224c^2p^9-2046c^5p^5+4258c^4p^6+832p^7c^4 -1464p^8c^2-2289c^6p^2+1544p^7c^3 \nonumber\\ && -450c^6p-1344cp^9)\mu_1^2/ \{8\{(2p+2)c^2-(p^2+3p)c -2p^3\}^8\} +c^2p^4(c-p)^2(40p^2c^4 \nonumber\\ && +61c^4p+24c^4-78p^3c^3-158c^3p^2 -68c^3p-14c^2p^4 +43p^3c^2+48p^2c^2+32p^5c \nonumber\\ && +62cp^4 +24p^6)\mu_1\mu_2 /\{4\{(2p+2)c^2-(p^2+3p)c-2p^3\}^5\} +c^2p^2(c-p)\mu_2^2 \nonumber\\ && /\{(2p+2)c^2-(p^2+3p)c-2p^3\}^2 +o(|\mu_1, \mu_2|^2). \end{eqnarray*}$$ $

(28)

Then we are ready to express those bifurcation curves in parameters $\varepsilon_1$ and $\varepsilon_2$.

For curves ${\mathcal SN}^\pm$, we need to consider $\mu_1=0$. From the first equality of (19) we see that $\mu_1=0$ if and only if $\varepsilon_1=0$ because in the expression of $\phi_1(\varepsilon_1, \varepsilon_2)$ we have $ {\phi_{11}(0, 0)}/{\phi_{12}^2(0, 0)}={1}/{p^2c^2(p+1)}\neq0. $ Thus, for $\mu_1=0$ we obtain from (28) that

$ \varepsilon_2=-\frac{c^2}{(2p+2)\Psi(c)}\mu_2+O(|\mu_2|^2), \label{Psi(c)} $

(29)

where $\Psi(c)$ is the same quadratic polynomial as given in (12). It follows that the inequality $\mu_2>0$ (or $<0$) together with the sign of $\Psi(c)$ determines the sign of $\varepsilon_2$. From the analysis of the quadratic equation (12) we see that $\Psi(c)<0$ (or $>0$) if $0<c<\varsigma(b)$ (or $c>\varsigma(b)$), where $\varsigma(b)$ is defined in Lemma 2.2. Hence from (22) we obtain that

$ $$\begin{eqnarray*} \mathcal{SN}^{+}:&=&\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=0, \varepsilon_{2}>0, 0<c<\varsigma(b)\} \cup \{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=0, \varepsilon_{2}<0, c>\varsigma(b)\} , \nonumber\\ \mathcal{SN}^{-}:&=&\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=0, \varepsilon_{2}<0, 0<c<\varsigma(b)\} \cup \{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=0, \varepsilon_{2}>0, c>\varsigma(b)\} . \end{eqnarray*}$$ $

For curve $\mathcal{H}$, we need to consider $\mu_{1}=-\mu_{2}^{2}$, which is equivalent to $\Upsilon(\varepsilon_1, \varepsilon_2):= \phi_{1}(\varepsilon_{1}, \varepsilon_{2})+\phi_{2}^{2}(\varepsilon_{1}, \varepsilon_{2})=0$ by (19). Clearly, $\Upsilon(0, 0)=0$ and

$ \frac{\partial\Upsilon}{\partial \varepsilon_1}\Big|_{(\varepsilon_{1}, \varepsilon_{2})=(0, 0)} = \{ (2p+2)\Psi(c) \}^4 / \{p^6c^2(c-p)^4(p+1)\}\ne 0. $

By the Implicit Function Theorem, there exists a unique $C^k$ function $\varepsilon_1=\epsilon_1(\varepsilon_2)$ such that $\epsilon_1(0)=0$ and $\Upsilon(\epsilon_1(\varepsilon_2), \varepsilon_2)=0$. Similarly to (27) and (28), expanding $\Upsilon$ at $(\varepsilon_1, \varepsilon_2)=(0, 0)$ and substituting with a formal expansion of $\epsilon_1(\varepsilon_2)$ of order 2, we obtain by comparison of coefficients that

$ \varepsilon_1=\epsilon_1(\varepsilon_2) = -\frac{p^6(c-p)^4}{4(p+1)\Psi^2(c)}\varepsilon_2^2+o(|\varepsilon_2|^2). \label{Hp} $

(30)

Further, replacing $\mu_1$ with $\mu_{1}=-\mu_{2}^{2}$ in (28), we get

$ \varepsilon_2=-\frac{c^2}{(2p+2)\Psi(c)}\mu_2+o(|\mu_2|). $

Similarly to (29), from (22) we obtain that

$ $$\begin{eqnarray*} \mathcal{H} &:=& \big\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=-\frac{p^6(c-p)^4}{4(p+1)\Psi^2(c)}\varepsilon_{2}^2+o(|\varepsilon_2|^2) , ~~\varepsilon_{2}>0, 0<c<\varsigma(b) \big\} \\ &&\cup \big\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=-\frac{p^6(c-p)^4}{4(p+1)\Psi^2(c)}\varepsilon_{2}^2+o(|\varepsilon_2|^2) , ~~\varepsilon_{2}<0, c>\varsigma(b) \big\}. \end{eqnarray*}$$ $

For curve $\mathcal{L}$, we need to consider $\mu_{1}=-\frac{49}{25}\mu_{2}^{2} +o(|\mu_{2}|^{2})$, i.e., $\phi_{1}(\varepsilon_{1}, \varepsilon_{2}) =-\frac{49}{25}\phi_{2}^{2}(\varepsilon_{1}, \varepsilon_{2}) +o(|\phi_{2}|^{2})$. Similarly to $\mathcal{H}$, we apply the Implicit Function Theorem to obtain

$ \varepsilon_{1}=-\frac{49p^6(c-p)^4}{100(p+1)\Psi^2(c)}\varepsilon_{2}^2+o(|\varepsilon_{2}|^{2}). $

Similarly to (29), from (22) we obtain that

$ $$\begin{eqnarray*} \mathcal{L} &:=& \big\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=-\frac{49p^6(c-p)^4}{100(p+1)\Psi^2(c)}\varepsilon_{2}^2+o(|\varepsilon_{2}|^{2}), ~~\varepsilon_{2}>0, 0<c<\varsigma(b) \big\} \\ &&\cup \big\{(\varepsilon_{1}, \varepsilon_{2})\ |\ \varepsilon_{1}=-\frac{49p^6(c-p)^4}{100(p+1)\Psi^2(c)}\varepsilon_{2}^2+o(|\varepsilon_{2}|^{2}) , ~~\varepsilon_{2}<0, c>\varsigma(b) \big\}. \end{eqnarray*}$$ $

Finally, with the replacement (13) we can rewrite the above bifurcation curves ${\mathcal SN}^\pm, {\mathcal H}$ and ${\mathcal L}$ expressed in parameters $(\varepsilon_{1}, \varepsilon_{2})$ in expressions in the original parameters $(a, b, c, \kappa)$ as shown in Theorem 3.1.

---

4. Conclusions

In this paper we analyzed the dynamics of system (4) near the equilibrium $E_*$ when parameters lie near ${\mathcal B}\backslash{\mathcal C}$. We proved that $E_*$ is a cusp when parameters lie on ${\mathcal B}\backslash{\mathcal C}$. We investigated the Bogdanov-Takens bifurcation near the cusp and compute in Theorem 3.1 the four bifurcation curves $\mathcal{SN}^{+}$, $\mathcal{SN}^{-}$, $\mathcal{H}$ and $\mathcal{L}$ in the practical parameters. Those bifurcation curves can be observed in Figure 3 in the case that $c>1$ and $b=(c+1)^2/4$ (which implies $p=(c-1)/2$). They display the merge of equilibria and the rise of homoclinic orbits and periodic orbits.

More concretely, in this case,

$ a_*=\frac{(c+1)^2}{4}, ~\kappa_*=\frac{8c^2}{(c+1)(c-1)^2}. $

Moreover, the four bifurcation curves divide the neighborhood $U$ of $(a_*, \kappa_*)$ into the following regions:

$ $$\begin{eqnarray*} \mathcal{D}_{I} &:=& \Big\{(a, \kappa)\in U|~a<\frac{(c+1)^2}{4}, ~\kappa\le\frac{8c^2}{(c+1)(c-1)^2}\Big\} \nonumber \\ &&\bigcup \Big\{(a, \kappa)\in U|~a<\frac{(c+1)^2}{4}-\frac{49(c-1)^6(c+1)^3}{3200(2c^2+c+1)^2}\big\{\kappa-\frac{8c^2}{(c+1)(c-1)^2}\big\}^2 \nonumber\\ && +O(|\kappa-\frac{8c^2}{(c+1)(c-1)^2}|^{3}), ~\kappa>\frac{8c^2}{(c+1)(c-1)^2}\Big\}, \nonumber\\ \mathcal{D}_{II} &:=& \Big\{(a, \kappa)\in U|~ \frac{(c+1)^2}{4}-\frac{49(c-1)^6(c+1)^3}{3200(2c^2+c+1)^2}\big\{\kappa-\frac{8c^2}{(c+1)(c-1)^2}\big\}^2 \nonumber\\ && +O(|\kappa-\frac{8c^2}{(c+1)(c-1)^2}|^{3}) <a<\frac{(c+1)^2}{4}-\frac{(c-1)^6(c+1)^3}{128(2c^2+c+1)^2} \nonumber\\ &&\big\{\kappa-\frac{8c^2}{(c+1)(c-1)^2}\big\}^2 +O(|\kappa-\frac{8c^2}{(c+1)(c-1)^2}|^{3}), ~\kappa>\frac{8c^2}{(c+1)(c-1)^2}\Big\}, \nonumber\\ \mathcal{D}_{III} &:=& \Big\{(a, \kappa)\in U|~ \frac{(c+1)^2}{4}-\frac{(c-1)^6(c+1)^3}{128(2c^2+c+1)^2} \big\{\kappa-\frac{8c^2}{(c+1)(c-1)^2}\big\}^2 \\ && +O(|\kappa-\frac{8c^2}{(c+1)(c-1)^2}|^{3}) <a<\frac{(c+1)^2}{4}, ~\kappa>\frac{8c^2}{(c+1)(c-1)^2}\Big\}, \nonumber\\ \mathcal{D}_{IV} &:=& \Big\{(a, \kappa)\in U|~a>\frac{(c+1)^2}{4}\Big\}. \nonumber \end{eqnarray*}$$ $

Theorem 3.1 gives dynamical behaviors of system (4) near $E_*$ in the first quadrant in Table 4. The coordinates of equilibria $E_0:(x_0, 0)$, $E_1:(p_1, q_1)$ and $E_2:(p_2, q_2)$ are given by $x_0:=a/(b-1)$ and

$ p_1:=-\frac{1}{2}\big\{(a-b-c+1)-\{(a-b-c+1)^2-4(a-c)\}^{1/2}\big\}, \\ p_2:=-\frac{1}{2}\big\{(a-b-c+1)+\{(a-b-c+1)^2-4(a-c)\}^{1/2}\big\} $

as in [27]. $E_0$ exists in the first quadrant when $(a, \kappa)\in\mathcal{D}_{I}\cup\mathcal{L}\cup\mathcal{D}_{II}\cup\mathcal{H}\cup\mathcal{D}_{III}$ but disappears when $(a, \kappa)\in\mathcal{D}_{IV}$ (appearing in other quadrants) or $(a, \kappa)\in\mathcal{SN}^+\cup\{(a_*, \kappa_*)\}\cup\mathcal{SN}^-$ (not existing).

The appearance of limit cycle displays a rise of oscillatory phenomenon in system (4). Choosing parameters $a=3.99999, b=4, c=3$ and $\kappa=4.495$ in $\mathcal{D}_{II}$, we used the command ODE45 in the software Matlab Version R2014a to simulate the orbit initiated from $(x_0, y_0)=(1.00432, 1.98662845)$ numerically, which plots an attractive limit cycle in Figure 4 and shows a dynamic balance and permanence of the substrate and the product in the enzyme-catalyzed reaction. The homoclinic loop actually gives a boundary for the break of the dynamic balance and permanence.

In this paper we only considered parameters in ${\mathcal B}\backslash {\mathcal C}$. When parameters lie in ${\mathcal C}$, higher degeneracy may happen at $E_*$. Although efforts have been made for higher degeneracies, for example, versal unfolding was discussed in [5] for a normal form of cusp system of codimension 3, it is still difficult to compute bifurcation curves in original parameters in the case of codimension 3. Such a computation with original parameters is indispensible for practical systems and for system (4) it will be our next work.

---

Appendix. Some coefficients

The functions in system (15) are

$ $$\begin{eqnarray*} E_{00} &:=& \{(2p+2)c^2-(p^2+3p)c-2p^3\}^4\varepsilon_1/\big\{c^2(p+1)p^6(c-p)^4\big\}, \\ E_{10} &:=& -\{(2p+2)c^2-(p^2+3p)c-2p^3\}^2\varepsilon_{1}\big\{(-6c^3p-4c^3p^2-4p^3c^2+3p^2c^2 +4cp^4+4c^4p \\ && +3c^4) -(p^2c^2-3c^3p-3c^2p+cp^2+2cp^3-2p^4)\varepsilon_1 -(p^3c^2-2cp^4+p^5+4c^2p^4 \\ && -5p^5c-p^3c^3+2p^6)\varepsilon_2\big\}/\big\{(p+1)p^4c^3(c-p)^4\big\}, \\ E_{20} &:=& \big\{(-2c^6(p+1)^2(c-p)^2) +(9c^3p^2+4c^2p^4-13c^4p+4p^5c^2+6p^3c^3 +9c^5p-15p^2c^4 \\ && -2p^4c^3+4p^2c^5-4p^3c^4+6c^5)\varepsilon_1 -(2p^7c-6p^7c^2-6p^6c^2-2p^5c^4+6p^6c^3 +2cp^8 \\ && -2p^4c^4+6p^5c^3)\varepsilon_2 +(6p^5c^2-2p^4c^3-6p^6c-6p^7c+6p^6c^2-2p^5c^3+2p^7+2p^8)\varepsilon_1\varepsilon_2 \\ && +(6p^3c^3-4p^2c^4-2c^2p^4-10p^3c^2-9c^4p-2cp^4+17c^3p^2 -2p^5c+13c^3p-9p^2c^2 \\ && -6c^4)\varepsilon_1^2\big\} /\big\{2c^3p^2(c-p)^2(p+1)\big\}, \\ E_{01} &:=& -\{(2p+2)c^2-(p^2+3p)c-2p^3\}\{2c^3\varepsilon_1 +(cp^4-2p^3c^2+c^3p^2)\varepsilon_2 +(2p^4-6cp^3 \\ && +4p^2c^2)\varepsilon_1\varepsilon_2 +(12c^2-6cp)\varepsilon_1^2\}/\big\{p^2(c-p)^2c^3\big\}, \\ E_{11} &:=& \big\{(3c^3p^2-8p^2c^4-p^4c^3+2c^5+2c^2p^4+4c^5p+2p^2c^5-5c^4p+2p^5c^2 +2p^3c^3 \\ && -3p^3c^4) +(3c^2p^4+3c^3p+p^2c^2+2p^5c+3p^2c^4+3p^3c^2+2c^4p+2cp^4 -4p^3c^3 \\ && +c^3p^2)\varepsilon_1 +(5p^6c^2-2p^7c-3p^6c+7p^5c^2 +2c^2p^4 -p^5c-5p^4c^3-p^3c^3+p^3c^4 \\ && +p^4c^4-4p^5c^3)\varepsilon_2 -(5p^6c-4p^5c^2+p^4c^3-5c^2p^4-p^3c^2+7p^5c+p^3c^3 +2cp^4 \\ && -2p^7-p^5-3p^6)\varepsilon_1\varepsilon_2 +(13cp^2-8cp^4+9c^3p^2-38p^2c^2 +5cp^3+10p^4+10p^5 \\ && +19c^3p+10c^3-13p^3c^2 -25c^2p)\varepsilon_1^2\big\} /\big\{c^2(p+1)(-p+c)\big\}, \\ E_{02} &:=& \big\{(c-2p-1)+(5c^3-2c^2p)\varepsilon_1-(3p^3c^2-2c^3p^2-cp^4)\varepsilon_2 +(p^4-cp^3)\varepsilon_1\varepsilon_2-(2cp \\ && -c^2)\varepsilon_1^2\big\} /\big\{(p+1)^2(c-p)^2\big\}. \end{eqnarray*}$$ $

The functions below system (18) are

$ $$\begin{eqnarray*} \phi_{11} &:=& 24c^6p^5+4c^8p^2-16c^7p^4+4c^8p^3-16c^5p^6+4p^7c^4 +24c^6p^4-16c^7p^3-16c^5p^5 \\ && +4c^4p^6 +(9p^4c^4-16p^6c^3+40c^3p^7+68p^5c^4-26p^3c^5+3c^8 -6c^8p +42c^6p^3+36c^6p^4 \\ && -94c^5p^4 +6c^7p^2-4c^4p^6-16c^2p^8 -56c^5p^5-8c^8p^2+8c^7p^3+28c^6p^2-14c^7p)\varepsilon_1 \\ && +(4c^7p^4+40c^5p^7 -4c^2p^9-4c^2p^{10}+20c^3p^8+20c^3p^9 -20c^6p^5+40c^5p^6-20c^6p^6 \\ && -40c^4p^8+4c^7p^5-40p^7c^4)\varepsilon_2 -(40p^2c^5+12p^4c^3+32c^7p^2 +8p^5c^3+12c^7+92p^3c^5 \\ && +8p^6c^2-32p^3c^4-12p^6c^3-28p^7c^2 +4c^5p^4-88c^6p^2 -56p^4c^4+36c^7p+48p^5c^4 \\ && -60c^6p^3+16cp^8-32c^6p)\varepsilon_1^2 +(12cp^9-24p^7c^4 -8c^7p^4-6c^7p^3-88c^5p^5-32c^2p^8 \\ && +20cp^{10} -24c^5p^4+6c^6p^3+2c^3p^7-24p^6c^3+36c^6p^5 +72c^4p^6+96c^3p^8 -76c^2p^9 \\ && +36p^5c^4+40c^6p^4-44c^5p^6+6p^7c^2)\varepsilon_1\varepsilon_2 +(8p^7c-9p^2c^4 -16p^5c^2+6p^3c^3-c^2p^4 \\ && +11p^4c^4 +6p^3c^5 +10p^4c^3-16p^5c^3-18p^2c^5+12p^3c^4 -9c^6p^2+4p^6c-4p^8)\varepsilon_1^3 \\ && +(-34c^4p^6+2c^3p^7+4p^9-28cp^9-16cp^8+6p^4c^4 +8p^{10} -2p^7c+32c^2p^8 +32p^7c^2 \\ && -32p^6c^3-14p^5c^3 +10p^6c^2-6c^6p^4 +26c^5p^5+12p^5c^4)\varepsilon_1^2\varepsilon_2+(4c^3p^7-c^6p^6+44c^3p^9 \\ && -41c^2p^{10} -c^4p^6+4cp^9 +2c^5p^6-4p^{11}+28c^3p^8-32c^2p^9+8c^5p^7-26c^4p^8 +20p^{11}c \\ && -12p^7c^4 -p^{10}-4p^{12}+18cp^{10}-6c^2p^8)\varepsilon_1\varepsilon_2^2, \\ \phi_{12} &:=& (2p^5c^3-4p^3c^4+2p^4c^3+2p^3c^5+2p^2c^5-4p^4c^4) +(9c^3p^2+4c^2p^4-13c^4p+4p^5c^2 \\ && +6p^3c^3 +9c^5p-15p^2c^4-2p^4c^3+4p^2c^5-4p^3c^4+6c^5)\varepsilon_1 +(-2p^7c+6p^7c^2+6p^6c^2 \\ && +2p^5c^4-6p^6c^3 -2cp^8 +2p^4c^4-6p^5c^3)\varepsilon_2 +(6p^3c^3-4p^2c^4-2c^2p^4-10p^3c^2 -9c^4p \\ && -2cp^4+17c^3p^2-2p^5c +13c^3p -9p^2c^2-6c^4)\varepsilon_1^2 +(6p^5c^2 -2p^4c^3-6p^6c-6p^7c \\ && +6p^6c^2-2p^5c^3+2p^7+2p^8)\varepsilon_1\varepsilon_2, \\ \phi_{21} &:=& (6c^{10}+12c^8p^5+69c^8p^4-77c^9p^3+20c^7p^6 +9c^6p^4-33c^7p^3+18c^5p^6-34c^6p^5 \\ && +45c^8p^2 -26c^7p^4+102c^8p^3-27c^9p-80c^9p^2+27c^7p^5 +6c^5p^7+8c^4p^8-12c^5p^8 \\ && -55c^6p^6-12c^6p^7 +8c^4p^9 +22c^{10}p^2-24c^9p^4+8c^{10}p^3+20c^{10}p)\varepsilon_1 +(4p^{10}c^4+20p^9c^6 \\ && -10p^{10}c^5-4p^5c^9 +2p^{11}c^4 -2c^9p^4-2p^6c^9+10p^7c^8+20p^6c^8+10c^8p^5 -20p^9c^5 \\ && +2c^4p^9-40p^7c^7-20c^7p^8 -10c^5p^8 -20c^7p^6 +40c^6p^8+20c^6p^7)\varepsilon_2 +(-12c^9+12c^3p^9 \\ && -47c^8p^4+10c^9p^3-86c^6p^4 -19c^7p^3+102c^5p^6-220c^6p^5+60c^8p^2+159c^7p^4 -40c^8p^3 \\ && +61c^5p^5+2c^4p^6-16p^7c^4 -18c^9p +3c^9p^2+92c^7p^5 +12c^3p^8+26c^5p^7-14c^4p^8 \\ && +53c^8p+35c^6p^3-76c^7p^2-79c^6p^6)\varepsilon_1^2 +(2p^5c^9-34c^8p^5+2c^3p^9+19c^8p^4 -10c^9p^3 \\ && +151c^7p^6-17c^5p^6 +39c^6p^5-45c^7p^4 +26c^8p^3 -2c^3p^{10}+3p^7c^4-6c^9p^2+23c^7p^5 \\ && +77c^5p^7-26c^4p^8+145c^5p^8-85c^6p^6-227c^6p^7-31c^4p^9 -2c^9p^4 -103c^6p^8+83p^7c^7 \\ && +51p^9c^5-2p^{10}c^4-4p^{11}c^3-27p^6c^8)\varepsilon_1\varepsilon_2 +(-4p^7c^8-2p^6c^8 -30p^{10}c^4-60p^{11}c^4 \\ && +40p^9c^5-2p^8c^8+12p^{11}c^3 +24c^7p^8+12c^7p^9-60p^9c^6+12p^{13}c^3 -30p^{10}c^6+40p^{11}c^5 \\ && -30p^{12}c^4 -4p^{13}c^2-2p^{14}c^2-2p^{12}c^2+24p^{12}c^3-30c^6p^8+80p^{10}c^5 +12p^7c^7)\varepsilon_2^2 \\ && +(-30c^8+69p^3c^5-16p^4c^4-212p^5c^4+58p^6c^3+331c^5p^4-232c^6p^4 +79c^7p^3 \\ && +117c^5p^6-65c^6p^5-21c^8p^2-3c^7p^4+5c^8p^3+379c^5p^5 -187c^4p^6+4p^7c^4-94c^3p^8 \\ && +44c^2p^8 +44c^2p^9+91c^7p-53c^8p -263c^6p^3+163c^7p^2-106c^6p^2-38c^3p^7)\varepsilon_1^3 \\ && +(36c^2p^{10} -10c^8p^5 -199c^3p^9-41c^8p^4+9c^7p^6+18c^2p^{11}-166c^6p^4+84c^7p^3 \\ && +165c^5p^6-297c^6p^5-18c^8p^2 +193c^7p^4-48c^8p^3-110c^3p^{10} +164c^5p^5-76c^4p^6 \\ && +47p^7c^4+123c^7p^5-78c^3p^8+18c^2p^9 -208c^5p^7 +351c^4p^8-219c^5p^8-62c^6p^6 \\ && +79c^6p^7+233c^4p^9+12c^3p^7)\varepsilon_1^2\varepsilon_2 +(-2p^{14}c+c^2p^{10} -4c^3p^9-8c^7p^6+21c^2p^{11} \\ && -2p^{12}c-72c^3p^{10}-4p^{13}c -4c^5p^7+6c^4p^8-102c^5p^8+c^6p^6 +45c^6p^7+118c^4p^9 \\ && +58p^9c^6 -114p^{10}c^5+121p^{11}c^4+102c^6p^8-22p^7c^7-212p^9c^5+233p^{10}c^4 -138p^{11}c^3 \\ && +40p^{12}c^2-70p^{12}c^3+20p^{13}c^2-14c^7p^8+p^7c^8+p^6c^8)\varepsilon_1\varepsilon_2^2 (-176p^3c^4+41p^4c^3 \\ && +769p^3c^5+293p^2c^5-388p^4c^4+27p^5c^3 +28p^6c^2-178p^5c^4-58p^6c^3+4p^7c^2+20cp^8 \\ && +20cp^9+603c^5p^4 -192c^6p^4+75c^7p^3+127c^5p^5+34c^4p^6+72c^7-234c^6p-24c^2p^8 \\ && +210c^7p-616c^6p^3+213c^7p^2-658c^6p^2-44c^3p^7)\varepsilon_1^4 +(-286p^5c^4+154p^6c^3-32p^7c^2 \\ && +136c^2p^{10}-56cp^{10}-198c^3p^9 -24cp^9+262c^5p^4-32cp^{11}-330c^6p^4+70c^7p^3 \\ && +438c^5p^6 -284c^6p^5 +68c^7p^4+636c^5p^5-580c^4p^6-210p^7c^4+22c^7p^5-154c^3p^8 \\ && +30c^2p^8 +198c^2p^9 +64c^5p^7+84c^4p^8-122c^6p^3+24c^7p^2-76c^6p^6+198c^3p^7)\varepsilon_1^3\varepsilon_2 \\ && +(4p^{12}+102c^2p^{10}-cp^{10} -158c^3p^9-c^7p^6+198c^2p^{11}-64p^{12}c -33cp^{11}-c^5p^6 \\ && -313c^3p^{10}-32p^13c+4p^7c^4+8p^{13} -6c^3p^8+4c^2p^9 -57c^5p^7+132c^4p^8-126c^5p^8 \\ && +10c^6p^6+26c^6p^7+272c^4p^9+4p^{14} +16c^6p^8 -p^7c^7-70p^9c^5+144p^{10}c^4-161p^{11}c^3 \\ && +100p^{12}c^2)\varepsilon_1^2\varepsilon_2^2. \end{eqnarray*}$$ $

---