**Hopf bifurcation analysis**

   The differential system:

    f1[x_,y_]:=a x (1 - x/k) - b x y; 
    f2[x_,y_]:=-c y + d x y; 
    F[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@{f1[x,y],f2[x,y]};
    X={x, y};
    μ={a,b,c,d,k};
$$
\begin{align}
 &\dot{x}=a x\left(1-\frac{x}{k} \right)- bxy\\
 &\dot{y}= dxy - cy
\end{align}
$$

The Jacobian matrix:

    J[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@D[F[X,μ],{X}]

   The non-trivial equilibrium point:

    X0=Normal[Simplify[SolveValues[F[X,μ]==0&&Variables[F[X,μ]]>0,X]]][[1]];
    MatrixForm@X0
$$
\begin{align}
 P_{0}(x,y)=\left(\frac{c}{d} ,\frac{a}{b}\left(1-\frac{c}{dk} \right)\right)
\end{align}
$$
   The linear approximation at $P_{0}$ (coexistence equilibrium point):

    J0=Simplify@J[X0,μ];
    MatrixForm@J0
$$
\begin{align}
 J(P_{0})=\left(
\begin{array}{cc}
\hspace{-0.25cm}-\displaystyle\frac{a c}{d k} & -\displaystyle\frac{b c}{d}\hspace{0.3cm} \\\\
 \hspace{0.2cm}\displaystyle\frac{a (d k-c)}{b k} &\hspace{0.2cm} 0 \\
\end{array}
\right)
\end{align}
$$
Under the Hopf bifurcation conditions, $\text{tr}(J(P_{0},\mu_{0}))=0$ and $\text{det}(J(P_{0},\mu_{0}))>0$, where $\mu_{0}$ is the critical bifurcation value for some parameter of our system. In our case, the parameters are strictly positive and $\text{tr}(J(P_{0}))$ cannot be zero. Therefore, Hopf bifurcation not take place at $P_{0}$. The non-trivial equilibrium $P_{0}$ is always locally stable and the only condition that must be fulfilled is given by the following inequality
$$
\frac{c}{d k}<1
$$
Code for time series and phase portrait:

Time series

    s = ParametricNDSolve[{x'[t] == a x[t] (1 - x[t]/k) - bx[t]*y[t], 
    y'[t] == -c y[t] + d x[t]*y[t], x[0] == 11/5, y[0] == 4/5}, {x, y}, {t, 0, 1000}, {a, b, c, d, k}];
    Plot[Evaluate[x[1/4, 1/3, 1/2, 1/4, 10][t] /. s], {t, 0, 300}, PlotRange -> All, PlotPoints -> 500, 
    PlotStyle -> {Blue, Thickness[0.003]},AxesStyle -> Directive[Black, Small], Background -> Lighter[Gray, 0.95]]

[![Time series][1]][1]

Phase portrait:

    ParametricPlot[Evaluate[{x[1/4, 1/3, 1/2, 1/4, 10][t], y[1/4, 1/3, 1/2, 1/4, 10][t]} /. s], 
    {t, 0, 300}, PlotRange -> All, PlotPoints -> 500, PlotStyle -> {Blue, Thickness[0.003]}, 
    AxesStyle -> Directive[Black, Small],Background -> Lighter[Gray, 0.95]]

[![Phase portrait][2]][2]

**Example: Hopf bifurcation in the Brusselator system**

`Calculation of the first Lyapunov coefficient`

The Brusselator system is given by:
$$
\begin{align}
 &\dot{x}=\alpha-(\beta+1)x + x^2 y\\
 &\dot{y}= \beta x - x^2 y
\end{align}
$$

Assuming $\alpha> 0$ fixed and taking $\beta$ as a bifurcation parameter, we show that at $\beta = 1 + \alpha^2$ the system exhibits a supercritical Hopf bifurcation.

The Brusselator system code:

    f1[x_, y_] := α - (β + 1) x + x^2 y;
    f2[x_, y_] := β x - x^2 y;
    F[{x_, y_}, {α_, β_}] := Evaluate@{f1[x, y], f2[x, y]};
    X = {x, y};
    μ = {α, β};

The Jacobian matrix and its transpose:

    J[{x_, y_}, {α_, β_}] = D[F[X, μ], {X}];
    Jt[{x_, y_}, {α_, β_}] = Transpose[J[X, μ]];
    MatrixForm[J[X, μ]]
    MatrixForm[Jt[X, μ]]

Stability analysis (Routh-Hurwitz criterion):

    X0[{α_, β_}] = SolveValues[F[X, μ] == 0, X][[1]]
    polJX0 = Collect[CharacteristicPolynomial[J[X0[μ], μ], λ], λ,Simplify];
    a0 = CoefficientList[polJX0, λ][[3]];
    a1 = CoefficientList[polJX0, λ][[2]];
    a2 = CoefficientList[polJX0, λ][[1]];
    Reduce[a1 > 0 && a2 > 0 && α > 0 && β > 0, β]
    (*α > 0 && 0 < β < 1 + α^2*)

Note that $a_{1}=0$ if and only if  $β=β_{0}$, where $β_{0}=1+α^2$. Then, the Brusselator is locally asymptotically stable at $X_{0}(\mu)$ for $β<β_{0}$ and locally asymptotically unstable for $β>β_{0}$ (appears a stable limit cycle surrounded the unstable equilibrium point). We verify the previous conclusion (transversality condition) with the sign of the following derivative:

    D[-a1, β]
    (*1*)

The analysis at the critical bifurcation value $β_{0}$:

Solve the following system of equations:
$$
\left\{\begin{align}
F\left((x,y),(\alpha,\beta)\right) &=0, \\
\operatorname{tr}(J((x,y),(\alpha,\beta))) &=0,
\end{align}\right.
$$
for $(x,y,\beta)$ and we must check that det $\operatorname{det}J((x,y),(\alpha,\beta))>0$ when $\beta = \beta_{0}$ for the solution found, where $\beta_{0}$ is the Hopf critical bifurcation value.

The code for the above system of equations:

    X0μ0 = Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, β}], 1], {3}]
    μ0 = Prepend[Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, β}], 1], {{1}, {2}}], α]
    Det[J[X0μ0, μ0]]
Here, the Hopf critical bifurcation value is $\beta_{0}=1+\alpha^2$ and $\operatorname{det}J((x,y),(\alpha,\beta_{0}))=\alpha^2>0$. Thus, the Brusselator at $\beta_{0}=1+\alpha^2$ has the equilibrium
$$
\begin{align}
X_{0}(\mu_{0})=\left(\alpha, \displaystyle\frac{1+\alpha^2}{\alpha} \right)
\end{align}
$$
and the linear approximation at $X_{0}(\mu_{0})$ has purely imaginary eigenvalues $\lambda_{1,2}=\pm \omega i$, $\omega=\alpha$.

The code for the linear approximation and its transpose at $X_{0}(\mu_{0})$:

    α= ω;
    JX0μ0 = Simplify@J[X0μ0, μ0];
    JtX0μ0 = Simplify@Transpose@JX0μ0;
    MatrixForm@JX0μ0
    MatrixForm@JtX0μ0
    Eigenvalues[JX0μ0]

The next step is to translate the equilibrium $X_{0}(\mu_{0})$ to the origin of coordinates:

    bb = {0, 0};
    F0[{x_, y_}, {α_, β_}] = Collect[Expand@F[X + X0μ0, μ0], {x, x^2, y, y^2, x y, x^2 y},Factor]
    MatrixForm@F0[bb,μ]

Now, to  obtain the normal form of the Hopf bifurcation, we need the Taylor expansion of the third order for $F_{0}((x,y),(\alpha,\beta))$:

    (*Rank 3 tensor*)
    D2[{x_, y_}, {α_, β_}] = Simplify@D[F0[X, μ], {X, 2}]
    D2X0μ0 = Simplify@D2[bb, μ0]
    (*Rank 4 tensor*)
    D3[{x_, y_}, {α_, β_}] = Simplify@D[F0[X, μ], {X, 3}]
    D3X0μ0 = Simplify@D3[bb, μ0]

Multilinear forms:

    (*Bilinear form*)
    BB[{x_, y_}, {u_, v_}] = Collect[Expand[D2X0μ0.X.U], {u x, v y, v x, v y}, FullSimplify];
    MatrixForm[BB[{x, y}, {u, v}]]
    (*Trilinear form*)
    CC[{x_, y_}, {u_, v_}, {r_, s_}] = Collect[Expand[D3X0μ0.X.U.R], {r u x, r u y, r v x, r v y, s u x, s v x, s u y, s v y}, FullSimplify]
    MatrixForm[CC[{x, y}, {u, v}, {r, s}]]

We verify that the first three terms of the Taylor series expansion of $F_{0}((x,y),(\alpha,\beta))$ are correct:

    MatrixForm@FullSimplify[F0[X, μ] - (JX0μ0.X + 1/2! BB[X, X] + 1/3! CC[X, X, X]) /. {x -> t x, y -> t y}]
    (*{0,0}*)

Now, we compute the critical eigenvectors of $J((0,0),\mu_{0})$ and its transpose:

    (*Eigenvectors for J[X0,μ0]*)
    vp = ComplexExpand@Eigenvectors[JX0μ0]
    q = vp[[2]];
    qc = vp[[1]];
    MatrixForm@q
    MatrixForm@qc
    MatrixForm@Simplify[JX0μ0.q - I ω q]
    (*Eigenvectors for Transpose[J[X0,μ0]]*)
    vpt = ComplexExpand[Eigenvectors[JtX0μ0]]
    (*Normalization constant*)
    cn = ComplexExpand[Conjugate[vp[[1]] . vpt[[1]]]]
    p = Expand@Simplify[vpt[[2]]/cn];
    pc = ComplexExpand[Conjugate[p]];
    MatrixForm@p
    MatrixForm@pc
    Simplify[JtX0μ0.p-I ω p]

We verify the normalization condition $\langle p,q\rangle=1$

    Simplify@(p.q)
    (*1*)
Finally, we compute the first Lyapunov coefficient:
$$
\begin{align}
l_1(0,\mu_{0})= &\frac{1}{2\omega_0} {\rm Re}\left[\langle p,C(q,q,\bar{q}) \rangle - 2 \langle p, B(q,A_0^{-1}B(q,\bar{q}))\rangle 
+\\\hspace{0.5cm} \langle p, B(\bar{q},(2i\omega_0 I_n-A_0)^{-1}B(q,q))\rangle \right]
\end{align}
$$

Before to calculate $l_1(0,\mu_{0})$, we clean $\alpha$

    Clear[α]
    ω0=α;

The code for the first Lyapunov coefficient:

    Factor@ComplexExpand[Re[1/(2 ω) (p.D3X0μ0.q.q.qc- 2(p.D2X0μ0.q.Inverse[JX0μ0].D2X0μ0.q. qc) +
    p.D2X0μ0.qc. Inverse[2 I ω*IdentityMatrix[2]-JX0μ0].D2X0μ0.q. q)] /. ω -> ω0]
    (*-((2 + α^2)/(2 α (1 + α^2)))*)

$$
\begin{align}
l_1(0,\mu_{0})=-\frac{\alpha^2+2}{2 \alpha \left(\alpha^2+1\right)}
\end{align}
$$
The first Lyapunov coefficient is clearly negative for all positive $\alpha$. Thus, the Hopf bifurcation is nondegenerate and always supercritical.

The above expression is the result that Kuznetsov arrives at on page 105 of his book (see [Elements of Applied Bifurcation Theory][3]).

Limit cycle:

[![Limit cycle][4]][4]

For more details see: [Andronov-Hopf bifurcation][5].

###### In case anyone wishes to collaborate on these topics, please, communicate.


  [1]: https://i.sstatic.net/pRAqT.png
  [2]: https://i.sstatic.net/rQ6Nk.png
  [3]: https://doi.org/10.1007/978-1-4757-3978-7
  [4]: https://i.sstatic.net/SCEwC.png
  [5]: http://www.scholarpedia.org/article/Andronov-Hopf_bifurcation