Bifurcation Analysis and Its Applications

1. Introduction

Continuous dynamical systems that involve differential equations mostly contain parameters. It can happen that a slight variation in a parameter can have significant impact on the solution. The main questions of interest in this chapter are: How to continue equilibria and periodic orbits of dynamical systems with respect to a parameter? How to compute stability boundaries of equilibria and limit cycles in the parameter space? How to predict qualitative changes in system’s behavior (bifurcations) occurring at these equilibrium points? This chapter will also cover the classification of bifurcations in terms of equilibria and periodic orbits. Especially it will present the specific bifurcation called ” Hopf bifurcation” which refers to the development of periodic orbits from stable equilibrium point, as a bifurcation parameter crosses a critical value. Since the theory of bifurcation from equilibria based on center manifold reduction and Poincaré-Normal forms, the direction of bifurcations for the mathematical models will also be explained using this theory. Finally, by introducing several software packages and numerical methods this chapter will also cover the techniques to determine and continue in some control parameters all local bifurcations of periodic orbits of dynamical systems and relevant normal form computations combined with the center manifold theorem, including periodic normal forms for periodic orbits.

In general, in a dynamical system, a parameter is allowed to vary, then the differential system may change. An equilibrium can become unstable and a periodic solution may appear or a new stable equilibrium may appear making the previous equilibrium unstable. The value of parameter at which these changes occur is known as ”bifurcation value” and the parameter that is varied is known as the ”bifurcation parameter”.

In this chapter, we also discuss several types of bifurcations, saddle node, transcritical, pitchfork and Hopf bifurcation. Among these types, we especially focus on Hopf bifurcation. The first three types of bifurcation occur in scalar and in systems of differential equations. The fourth type called Hopf bifurcation does not occur in scalar differential equations because this type of bifurcation involves a change to a periodic solution. Scalar autonomous differential equations can not have periodic solutions. Hopf bifurcation occurs in systems of differential equations consisting of two or more equations. This type is also referred to as a ”Poincare-Andronov-Hopf bifurcation”.

For a given system of differential equations first we shall consider the stability and the local Hopf bifurcation. By using the Hopf bifurcation theorem we prove the occurrence of the Hopf bifurcation. And then, based on the normal form method and the center manifold reduction introduced by Hassard et al.,[10], we derive the formulae determining the direction, stability and the period of the bifurcating periodic solution at the critical value of the bifurcation parameter. To verify the theoretical analysis, numerical simulations for bifurcation analysis are given in this chapter. For references see [1]-[22].

We also introduce the Hopf bifurcation for continuous dynamical systems and state the Hopf bifurcation theorem for these models. As it is well known, Hopf bifurcations occur when a conjugated complex pair of eigenvalues crosses the boundary of stability. In the time-continuous case, a limit cycle bifurcates. It has an angular frequency which is given by the imaginary part of the crossing pair. In the discrete case, the bifurcating orbit is generally quasi-periodic, except that the argument of the crossing pair times an integer gives just $2Ϟ$. If we consider an ordinary differential equation (ODE) that depends on one or more parameters $Ϗ$

$$x^{\mathrm{\text{'}}}=f(x,Ϗ),$$E1

where, for simplicity, we assume $Ϗ$ to be the only parameter. There is the possibility that under variation of $Ϗ$ nothing interesting happens to Equation (1). There is only a quantitatively different behavior. Let us define Equation (1) to be structurally stable in the case there are no qualitative changes occurring. However, the ODE (Ordinary Differential Equation) might change qualitatively. At that point, bifurcations will have occurred.

Many of the basic principles for one dimensional systems apply also for two-dimensional systems. Let us define a two-dimensional system

$$x^{\mathrm{\text{'}}}=f(x,y,Ϗ),$$E2

$y^{\mathrm{\text{'}}}=g(x,y,Ϗ),$

where biologically we mostly interpret x as prey or resource and y as predator or consumer. Equilibria can be found by taking the equations equal to zero, i.e.,

$f(x,y,Ϗ)=0,$$g(x,y,Ϗ)=0\mathrm{.}$

We have three possibilities for the stability of an equilibrium. Next to the stable and unstable equilibrium, there is the saddle equilibrium. A two-dimensional stable equilibrium is attracting in two directions, while a two-dimensional unstable equilibrium is repelling in two directions. A saddle point is attracting in one direction and repelling in the other direction. In the less formal literature saddles are often considered just unstable equilibria. A second remark is that also the dynamics of the system around the equilibria can differ.

The attracting or repelling can occur via straight orbits (a node) or via spiralling orbits (a spiral or focus). Note that, it is not possible to have a saddle focus in two dimensions. It is possible though in three of higher dimensional systems.

To prove the existence of Hopf bifurcation, we first obtain the Hopf bifurcation theorem hypothesis, i.e., the existence of purely imaginary eigenvalues of the corresponding characteristic equation with respect to the parameter $Ϗ$ and also we prove the transversality condition

$$\frac{dϙ}{dϏ}(Ϗ={Ϗ}_0)\stackrel{?}{=}0$$E3

at the critical value ${Ϗ}_0$ where Hopf bifurcation occurs. Then based on the normal form approach and the center manifold theory introduced by Hassard et. al,[10], we derive the formula for determining the properties of Hopf bifurcation of the model.

Finally in this chapter, to support these theoretical results, we illustrate them by numerical simulations. In numerical analysis, generally MATLAB solver packages are used to analyze the dynamics of nonlinear models. In these solvers, differential equation systems are simulated by difference equations. In this chapter, we also give some examples from biology (such as well known predator-prey models with time delay) with numerical simulations and by graphing the solutions in two or three dimensions, we illustrate the occurrence of periodic solutions. For some examples see references by Çelik, [4], [5] and [6].

2. Basic concepts of bifurcation analysis

As it is stated above, in dynamical systems, a bifurcation occurs when a small smooth change made to the parameter values (the bifurcation parameters) of a system causes a sudden ” qualitative" or topological change in its behavior. Generally, at a bifurcation, the local stability properties of equilibria, periodic orbits or other invariant sets changes. It has two types;

Local bifurcations, which can be analyzed entirely through changes in the local stability properties of equilibria, periodic orbits or other invariant sets as parameters cross through critical thresholds; and

Global bifurcations, which often occur when larger invariant sets of the system ”collide” with each other, or with equilibria of the system. They cannot be detected purely by a stability analysis of the equilibria (fixed or equilibrium points, see the next section).

2.1. Equilibrium points

In dynamical systems, only the solutions of linear systems may be found explicitly. The problem is that in general real life problems may only be modeled by nonlinear systems. The main idea is to approximate a nonlinear system by a linear one (around the equilibrium point). Of course, we do hope that the behavior of the solutions of the linear system will be the same as the nonlinear one. But this is not always true. Before the linear stability analysis, we give some basic definitions below.

Definition (Equilibrium Point): Consider a nonlinear differential equation

$x^{\mathrm{\text{'}}}\left(t\right)=f\left(x\right(t),u(t\left)\right),$

where $f$ is a function mapping $R^{n}x{R}^3\to R^n\mathrm{.}$ A point $\stackrel{\text{̣}}{x}$ is called an equilibrium point if there is a specific $\stackrel{\text{̣}}{u}ϓR^m$ such that

$f\left(x\right(t),u(t\left)\right)=0_n\mathrm{.}$

Suppose $\stackrel{\text{̣}}{x}$ is an equilibrium point (with the input $\stackrel{\text{̣}}{u}$). Consider the initial condition $x\left(0\right)=\stackrel{\text{̣}}{x}$, and applying the input $u\left(t\right)=\stackrel{\text{̣}}{u}$ for all $t\le t_0$, then resulting solution $x\left(t\right)$ satisfies

$x\left(t\right)=\stackrel{\text{̣}}{x},$

for all $t\le t_0$. That is why it is called an equilibrium point or solution.

Example: As an example, consider the logistic growth equation (the rate of population density)

$x^{\mathrm{\text{'}}}=\mathit{rx}(1-\frac{x}{K}),$

where $x\left(t\right)$ denotes the population density at time $t$, $r$ and $K$ are positive constants, $K$ is the carrying capacity. Then by setting right hand side function equal to zero,

$f\left(x\right)=\mathit{rx}(1-\frac{x}{K})=0,$

we obtain two equilibrium points $x=0$ and $x=K\mathrm{.}$

2.2. Linear stability analysis

Linear stability of dynamical equations can be analyzed in two parts: one for scalar equations and the other for two dimensional systems;

2.2.1. Linear stability analysis for scalar equations

To analyze the ODE

$x^{\mathrm{\text{'}}}=f\left(x\right)$

locally about the equilibrium point $x=\stackrel{\text{̣}}{x}$, we expand the function $f\left(x\right)$ in a Taylor series about the equilibrium point $\stackrel{\text{̣}}{x}$. To emphasize that we are doing a local analysis, it is customary to make a change of variables from the dependent variable $x$ to a local variable. Now let

$x\left(t\right)=\stackrel{\text{̣}}{x}+ϓ\left(t\right),$

where it is assumed that $ϓ\left(t\right)\gneqq 1$, so that we can justify dropping all terms of order two and higher in the expansion. Substituting $x\left(t\right)=\stackrel{\text{̣}}{x}+ϓ\left(t\right)$ into the RHS of the ODE yields;

$f\left(x\right(t\left)\right)=f(\stackrel{\text{̣}}{x}+ϓ(t\left)\right)=f\left(\stackrel{\text{̣}}{x}\right)+f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)ϓ\left(t\right)+f^{\mathrm{\text{'}\text{'}}}\left(\stackrel{\text{̣}}{x}\right)\frac{{ϓ}^2\left(t\right)}{2}+\mathrm{...}=0+f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)ϓ\left(t\right)+O\left({ϓ}^2\right),$

and dropping higher order terms, we obtain

$f\left(x\right)\ncong f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)ϓ\left(t\right)\mathrm{.}$

Note that dropping these higher order terms is valid since $ϓ\left(t\right)\gneqq 1\mathrm{.}$ Now substituting $x\left(t\right)=\stackrel{\text{̣}}{x}+ϓ\left(t\right)$ into the LHS of the ODE,

${ϓ}^{\mathrm{\text{'}}}\left(t\right)=f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)ϓ\left(t\right)\mathrm{.}$

The goal is to determine if we have growing or decaying solutions. If the solutions grows, then the equilibrium point is unstable. If the solution decays, then the fixed point is stable.

To determine whether or not the solution is stable or unstable we simply solve the ODE and get the solution as

$ϓ\left(t\right)={ϓ}_0\exp \left(f^{\mathrm{\text{'}}}\right(\stackrel{\text{̣}}{x}\left)t\right),$

where ${ϓ}_0$ is a constant. Hence, the solution is growing if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)>0$ and decaying if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)<0.$ As a result, the equilibrium point is stable if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)<0,$ unstable if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)>0$ as it is stated in the following theorem.

Theorem: Suppose for scalar differential equation

$x^{\mathrm{\text{'}}}=f\left(x\right),$

the derivative function $f^{\mathrm{\text{'}}}$ is continuous on an open interval $I$ where the equilibrium point $\stackrel{\text{̣}}{x}ϓI$. Then the equilibrium point $\stackrel{\text{̣}}{x}$ is locally stable if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)<0$ and it is unstable if $f^{\mathrm{\text{'}}}\left(\stackrel{\text{̣}}{x}\right)>0.$

If the equilibrium point is stable and in addition

$\underset{t\to \mathrm{\propto }}{\lim }x\left(t\right)=\stackrel{\text{̣}}{x},$

then it is called asymptotically stable equilibrium point.

Example: Determine the stability of the fixed points to the Logistic growth equation

$N^{\mathrm{\text{'}}}=f\left(N\right)=\mathit{rN}\left(1-\frac{N}{K}\right),\mathrm{where}r>0.$

We first find the equilibrium points by setting

$f\left(\stackrel{\text{̣}}{N}\right)=0,\mathrm{whichyieldstwopoints}\stackrel{\text{̣}}{N}=0,K\mathrm{.}$

Next we compute $f^{{}^{\mathrm{\text{'}}}}\left(N\right)$ and evaluate it at the equilibrium points.

$f^{\mathrm{\text{'}}}\left(N\right)=r(1-\frac{2N}{K})$$\mathrm{at}\stackrel{\text{̣}}{N}=0,f^{\mathrm{\text{'}}}\left(0\right)=r>0,\mathit{soitisunstable}$$\mathrm{at}\stackrel{\text{̣}}{N}=K,f^{\mathrm{\text{'}}}\left(K\right)=-r<0,\mathit{soitislocallystable.}$

2.2.2. Linear stability analysis for systems

Consider the two dimensional nonlinear system

$x^{\mathrm{\text{'}}}=f(x,y),$$y^{\mathrm{\text{'}}}=g(x,y),$

and suppose that $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})$ is a steady state (equilibrium point), i.e.,

$f\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)=0\mathrm{and}g(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})=0\mathrm{.}$

Now let’s consider a small perturbation from the steady state $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})$

$x=\stackrel{\text{̣}}{x}+u,$$y=\stackrel{\text{̣}}{y}+v,$

where $u$ and $v$ are understood to be small as $u\gneqq 1$ and $v\gneqq 1\mathrm{.}$ Is is natural to ask whether $u$ and $v$ are growing or decaying so that $x$ and $y$ will move away form the steady state or move towards the steady states. If it moves away, it is called unstable equilibrium point, if it moves towards the equilibrium point, then it is called stable equilibrium point.As in scalar equations, by expanding the Taylor’s series for $f(x,y)$ and $g(x,y);$

$u^{\mathrm{\text{'}}}=x^{\mathrm{\text{'}}}=f(x,y)$$=f(\stackrel{\text{̣}}{x}+u,\stackrel{\text{̣}}{y}+v)$$=f\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)+f_x\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)u+f_y\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)v+\mathrm{higherorderterms...}$$=f_x\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)u+f_y\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)v+\mathrm{higherorderterms....}$

Similarly,

$v^{\mathrm{\text{'}}}=y^{\mathrm{\text{'}}}=g(x,y)$$=g(\stackrel{\text{̣}}{x}+u,\stackrel{\text{̣}}{y}+v)$$=g\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)+g_x\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)u+g_y\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)v+\mathrm{higherorderterms...}$$=g_x\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)u+g_y\left(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y}\right)v+\mathrm{higherorderterms....}$

Since $u$ and $v$ are assumed to be small, the higher order terms are extremely small, we can neglect the higher order terms and obtain the following linear system of equations governingthe evolution of the perturbations $u$ and $v,$

$\left[\begin{array}{l}{u}^{\mathrm{\text{'}}}\\ v^{\mathrm{\text{'}}}\end{array}\right]=\left[\begin{array}{ll}{f}_x(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})& f_y(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})\\ g_x(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})& g_y(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})\end{array}\right]\left[\begin{array}{l}u\\ v\end{array}\right],$

where the matrix $\left[\begin{array}{ll}{f}_x& f_y\\ g_x& g_y\end{array}\right]$ is called Jacobian matrix $J$ of the nonlinear system.The above linear system for $u$ and $v$ has the trivial steady state $(u,v)=\left(\mathrm{0,0}\right)$, and the stability of this trivial steady state is determined by the eigenvalues of the Jacobian matrix, as follows:

Theorem: An equilibrium point $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})$ of the differential equation is stable if all the eigenvalues of $J$, the Jacobian evaluated at $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})$ have negative real parts. The equilibrium point is unstable if at least one of the eigenvalues has a positive real part.

As a summary,

Asymptotically stable : A critical point is asymptotically stable if all eigenvalues of the jacobian matrix $J$ are negative, or have negative real parts.

Unstable: A critical point is unstable if at least one eigenvalue of the jacobian matrix $J$ is positive, or has positive real part.

Stable (or neutrally stable) : Each trajectory move about the critical point within a finite range of distance.

Definition(Hyperbolic point): The equilibrium is said to be hyperbolic if all eigenvalues of the jacobian matrix have non-zero real parts.

Hyperbolic equilibria are robust(i.e., the system is structurally stable): Small perturbations of order do not change qualitatively the phase portrait near the equilibria. Moreover, local phase portrait of a hyperbolic equilibrium of a nonlinear system is equivalent to that of its linearization. This statement has a mathematically precise form known as the Hartman-Grobman. This theorem guarantees that the stability of the steady state $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})$ of the nonlinear system is the same as the stability of the trivial steady state $\left(\mathrm{0,0}\right)$ of the linearized system.

Definition(Non-Hyperbolic point): If at least one eigenvalue of the Jacobian matrix is zero or has a zero real part, then the equilibrium is said to be non-hyperbolic.

Non-hyperbolic equilibria are not robust (i.e., the system is not structurally stable). Small perturbations can result in a local bifurcation of a non-hyperbolic equilibrium, i.e., it can change stability, disappear, or split into many equilibria. Some refer to such an equilibrium by the name of the bifurcation (See the section below).

Example: Consider the following nonlinear autonomous system

$$\stackrel{\diamond }{x}\left(t\right)=y\left(t\right)\left[x\right(t)-y(t\left)\right],$$E4

$\stackrel{\diamond }{y}\left(t\right)=x\left(t\right)[2-y(t\left)\right]\mathrm{.}$

The equilibria are the points $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})=\left(\mathrm{0,0}\right)$ and $(\stackrel{\text{̣}}{x},\stackrel{\text{̣}}{y})=\left(\mathrm{1,2}\right)$ and the Jacobian matrix is

$J=\left[\begin{array}{ll}y& x-1\\ 2-y& -x\end{array}\right]\mathrm{.}$

We compute the Jacobian at the equilibrium point (0,0) where $J\left(\mathrm{0,0}\right)=\left[\begin{array}{ll}0& -1\\ 2& 0\end{array}\right]$ which implies that the eigenvalues are purely imaginary

${ϙ}_{1.2}=\pm \sqrt{2}i,$

by solving the characteristic equation

$d\mathit{et}(J-ϙI)=0\mathrm{.}$

Since the point is non-hyperbolic, the linearized system can not tell about the stability. Later on we will show that this is a center.

For the equilibrium point $\left(\mathrm{1,2}\right),$ the Jacobian $J\left(\mathrm{1,2}\right)=\left[\begin{array}{ll}2& 0\\ 0& -1\end{array}\right]$ thus the point is locally unstable (as

${ϙ}_1=1\mathrm{and}{ϙ}_2=-1,$

where one of the eigenvalues is strictly positive). Since it is a hyperbolic equilibrium point, the stability of fixed point is the same as in the linearized system. So it is also unstable.

3. Bifurcation analysis

Consider a family of ODE’s that depend on one parameter $ϙ$

$$x^{\mathrm{\text{'}}}=f(x,ϙ),$$E5

where $f:{\mathbb{R}}^{n+1}\to {\mathbb{R}}^n$ is analytic for $\mathit{ϙϓ}\mathbb{R}$, $xϓ{\mathbb{R}}^n$. Let $x=x_0\left(ϙ\right)$ be a family of equilibrium points of equation (5), i.e., $f\left(x_0\right(ϙ),ϙ)=0\mathrm{.}$ Now let’s set

$z=x-x_0\left(ϙ\right)\mathrm{.}$

Then

$z^{\mathrm{\text{'}}}=A\left(ϙ\right)z+O\left(\right|z{|}^2),$

where $A\left(ϙ\right)=\frac{\complement f}{\complement x}\left(x_0\right(ϙ),ϙ)\mathrm{.}$

Let ${ϙ}_1,{ϙ}_2,\mathrm{...}{ϙ}_n\left(ϙ\right)$ be the eigenvalues of $A\left(ϙ\right)\mathrm{.}$ If, for some $i$, $\mathit{Re}{ϙ}_i\left(ϙ\right)$ changes sign at $ϙ={ϙ}_0,$ we say that ${ϙ}_0$ is a bifurcation point of equation (5).

3.1. Bifurcation in one dimension

We may assume that $n=1$ so that $f:{\mathbb{R}}^2\to {\mathbb{R}}^1,$ and $x_0\left(ϙ\right)$ is a real valued analytic function of $ϙ$ provided

${ϙ}_1\left(ϙ\right)=\frac{\complement f}{\complement x}\left(x_0\right(ϙ),ϙ)=A\left(ϙ\right)=0\mathrm{.}$

Therefore, the equilibrium point is asymptotically stable if ${ϙ}_1\left(ϙ\right)<0,$ and unstable if ${ϙ}_1\left(ϙ\right)>0.$ This implies that ${ϙ}_0$ is a bifurcation point if ${ϙ}_1\left({ϙ}_0\right)=0\mathrm{.}$ Hence, bifurcation points $\left(x_0\right(ϙ),ϙ)$are solutions of

$f(x,ϙ)=0\mathrm{and}\frac{\complement f}{\complement x}(x,ϙ)=0\mathrm{.}$

The most common bifurcation types are illustrated by the following examples.

3.1.1. Saddle-node bifurcation

A saddle-node bifurcation or tangent bifurcation is a collision and disappearance of two equilibria in dynamical systems. In autonomous systems, this occurs when the critical equilibrium has one zero eigenvalue. This phenomenon is also called fold or limit point bifurcation.

Now consider the dynamical system defined by

$$x^{\mathrm{\text{'}}}=a-x^2,\mathrm{for}a\mathrm{isreal.}$$E6

An equilibrium solution (where $x\mathrm{\text{'}}=0)$ is simply $x=\pm \sqrt{a}\mathrm{.}$ Therefore,

if $a<0,$ then we have no real solutions,if $a>0,$ then we have two real solutions.

We now consider each of the two solutions for $a>0$, and examine their linear stability in the usual way. First, we add a small perturbation:

$x=\stackrel{\text{̣}}{x}+ϓ\mathrm{.}$

Substituting this into the equation yields

$\frac{dϓ}{\mathit{dt}}=(a-{\stackrel{\text{̣}}{x}}^2)-2\stackrel{\text{̣}}{x}ϓ-{ϓ}^2,$

and since the term in brackets on the RHS is trivially zero, therefore

$\frac{dϓ}{\mathit{dt}}=-2\stackrel{\text{̣}}{x}ϓ,$

which has the solution

$ϓ\left(t\right)=A\exp (-2\stackrel{\text{̣}}{x}t)\mathrm{.}$

From this, we see that for $x=+\sqrt{a}$, $\mathrm{}\left|x\right|\to 0$ as $t\to \mathrm{\propto }$ (linear stability); for $x=-\sqrt{a}$,$\left|x\right|\to 0$ as $t\to \mathrm{\propto }$ (linear stability).

As sketched in the “ bifurcation diagram” below, therefore, the saddle node bifurcation at $a=0$ corresponds to the creation of two new solution branches. One of these is linearly stable, the other is linearly unstable.

Figure 1.

Bifurcation diagram corresponding to the saddle-node bifurcation

3.1.2. Transcritical bifurcation

In a transcritical bifurcation, two families of fixed points collide and exchange their stability properties. The family that was stable before the bifurcation is unstable after it. The other fixed point goes from being unstable to being stable.

Now consider the dynamical system

$\frac{\mathit{dx}}{\mathit{dt}}=\mathit{az}-b{x}^2,\mathrm{for}x,a,b\mathrm{real.}$

Again, $a$ and $b$ are control parameters. We can find two steady states $(x^{\mathrm{\text{'}}}=0)$ to this system

$x={\stackrel{\text{̣}}{x}}_1=0,⇿a,b$$x={\stackrel{\text{̣}}{x}}_2=\frac{a}{b},⇿a,b,b\stackrel{?}{=}0\mathrm{.}$

We now examine the linear stability of each of these states in turn, following the usual procedure.

For the state ${\stackrel{\text{̣}}{x}}_1$, we add a small perturbation

$x={\stackrel{\text{̣}}{x}}_1+ϓ,$

which yields

$\frac{dϓ}{\mathit{dt}}=aϓ-b{ϓ}^2$

with the linearized form

$\frac{dϓ}{\mathit{dt}}=aϓ$

has the solution

$e\left(t\right)=A\exp \left(\mathit{at}\right)\mathrm{.}$

Therefore, perturbations grow for a > 0 and decay for a < 0. So

$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_1\mathrm{isunstableif}a>0,$$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_1\mathrm{isstableif}a<0.$

Now for the state ${\stackrel{\text{̣}}{x}}_2,$ we add a small perturbation

$x={\stackrel{\text{̣}}{x}}_2+ϓ,$

which yields the linearized form

$\frac{dϓ}{\mathit{dt}}=-aϓ$

has the solution

$e\left(t\right)=A\exp (-\mathit{at})\mathrm{.}$

Therefore, perturbations grow for a > 0 and decay for a < 0. So

$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_2=\frac{a}{b}\mathrm{islinearlystableif}a>0,$$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_2=\frac{a}{b}\mathrm{islinearlyunstableif}a<0.$

It can be easily seen that the bifurcation point $a=0$ corresponds to an exchange of stabilities between the two solution branches.

Figure 2.

Bifurcation diagram corresponding to the transcritical bifurcation

3.1.3. The pitchfork bifurcation

In pitchfork bifurcation one family of fixed points transfers its stability properties to two families after or before the bifurcation point. If this occurs after the bifurcation point, then pitchfork bifurcation is called supercritical. Similarly, a pitchfork bifurcation is called subcritical if the nontrivial fixed points occur for values of the parameter lower than the bifurcation value. In other words, the cases in which the emerging nontrivial equilibria are stable are called supercritical whereas the cases in which these equilibria are called subcritical.

Consider the dynamical system

$x^{\mathrm{\text{'}}}=\mathit{ax}-b{x}^3,\mathrm{for}a,b\mathrm{real.}$

As usual, $a$ and $b$ are external control parameters. Steady states, for which $x^{\mathrm{\text{'}}}=0$ are as follows:

$x={\stackrel{\text{̣}}{x}}_1=0,⇿a,b,$$x={\stackrel{\text{̣}}{x}}_2=-\sqrt{a/b},\mathrm{for}a/b>0,$$x={\stackrel{\text{̣}}{x}}_3=\mathrm{}\sqrt{a/b},\mathrm{for}a/b>0.$

Note that the equilibrium points ${\stackrel{\text{̣}}{x}}_2$ and ${\stackrel{\text{̣}}{x}}_3$ only exist when $a>0$ if $b>0$ and for $a<0$ if $b<0$.

As usual, we now examine the linear stability of each of these steady states in turn. (This can be done for a general $b$). First we write the perturbation for ${\stackrel{\text{̣}}{x}}_1=0,$

$x={\stackrel{\text{̣}}{x}}_1+ϓ$

that yields the linearized equation

$\frac{dϓ}{\mathit{dt}}=aϓ,$

with the solution

$e\left(t\right)=A\exp \left(\mathit{at}\right)\mathrm{.}$

So we see that

$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_1=0\mathrm{islinearlyunstableif}a>0,$$\mathrm{thestate}{\stackrel{\text{̣}}{x}}_2=0\mathrm{islinearlystableif}a<0.$

For the states $x={\stackrel{\text{̣}}{x}}_2$ and $x={\stackrel{\text{̣}}{x}}_3$, setting

$\stackrel{\text{̣}}{x}=\pm \sqrt{a/b}+ϓ,$$\frac{dϓ}{\mathit{dt}}=aϓ-3b{\stackrel{\text{̣}}{x}}^2ϓ,$

with the solution

$e\left(t\right)=A\exp \left(\mathit{ct}\right)\mathrm{where}c=-2a\mathrm{.}$

Thus it is obvious that

$\mathrm{thestates}{\stackrel{\text{̣}}{x}}_2\mathrm{and}{\stackrel{\text{̣}}{x}}_3\mathrm{arelinearlystableif}a>0,$$\mathrm{thestates}{\stackrel{\text{̣}}{x}}_2\mathrm{and}{\stackrel{\text{̣}}{x}}_3\mathrm{arelinearlyunstableif}a<0,$

Figure 3.

Bifurcation diagram corresponding to the pitchfork bifurcation

3.2. Bifurcation in two dimension

3.2.1. Hopf bifurcation

Definition: A Hopf or Poincare-Andronov-Hopf bifurcation is a local bifurcation in which a fixed point of a dynamical system loses stability as a pair of complex conjugate eigenvalues of linearization around the fixed point cross the imaginary axis of the complex plane.

Figure 4.

Hopf bifurcation diagram

3.3. Hopf bifurcation theorem

Consider the two dimensional system

$$\frac{\mathit{dx}}{\mathit{dt}}=f\left(x,y,Ϣ\right),$$E7

$\frac{\mathit{dy}}{\mathit{dt}}=g\left(x,y,Ϣ\right),$

where $Ϣ$ is the parameter and suppose that $\left(\overline{x}\left(Ϣ\right),\overline{y}\left(Ϣ\right)\right)$ is the equilibrium point and $Ϗ\left(Ϣ\right)\pm iϐ\left(Ϣ\right)$ are the eigenvalues of the Jacobian matrix which is evaluated at the equilibrium point.

In addition let’s assume that the change in the stability of the equilibrium point occurs at $Ϣ={Ϣ}^{\mathrm{*}}$ where $Ϗ\left({Ϣ}^{\mathrm{*}}\right)=0\mathrm{.}$

First the system is transformed so that the equilibrium is at the origin and the parameter $Ϣ$ at ${Ϣ}^{\mathrm{*}}=0$ gives purely imaginary eigenvalues. System (7) is rewritten as follows;

$$\frac{\mathit{dx}}{\mathit{dt}}=a_{11}\left(Ϣ\right)x+a_{12}\left(Ϣ\right)y+f_1(x,y,Ϣ),$$E8

$\frac{\mathit{dy}}{\mathit{dt}}=a_{21}\left(Ϣ\right)x+a_{22}\left(Ϣ\right)y+g_1\left(x,y,Ϣ\right)\mathrm{.}$

The linearization of the system (7) about the origin is given by $\frac{\mathit{dX}}{\mathit{dt}}=J\left(Ϣ\right)X$, where $X=\left[\begin{array}{l}x\\ y\end{array}\right]$ and

$J\left(Ϣ\right)=\left[\begin{array}{ll}{a}_{11}\left(Ϣ\right)& a_{12}\left(Ϣ\right)\\ a_{21}\left(Ϣ\right)& a_{22}\left(Ϣ\right)\end{array}\right]$

is the Jacobian matrix evaluated at origin.

Theorem (Hopf bifurcation theorem)

Let $f_1$ and $g_1$, in system (8) have continuous third order partial derivatives in $x$ and $y\mathrm{.}$ Suppose that the origin is an equilibrium point of (8) and that the Jacobian matrix $J\left(Ϣ\right)$ as above, is valid for all sufficiently small $\left|Ϣ\right|\mathrm{.}$ Moreover, assume that the eigenvalues of matrix $J\left(Ϣ\right)$ are $Ϗ\left(Ϣ\right)\pm iϐ\left(Ϣ\right)$ where $Ϗ\left(0\right)=0,ϐ\left(0\right)\stackrel{?}{=}0$ such that the eigenvalues cross the imaginary axis with nonzero speed, i.e., $\frac{dϏ}{dϢ}{|}_{Ϣ=0}\stackrel{?}{=}0\mathrm{.}$

Then in any open set $U$ containing the origin in ${\mathbb{R}}^2$ and for any ${Ϣ}_0>0,$ there exists a value $\stackrel{\text{̣}}{Ϣ},\left|\stackrel{\text{̣}}{Ϣ}\right|<{Ϣ}_0$ such that the system of differential equations (8) has a periodic solution for $Ϣ=\stackrel{\text{̣}}{Ϣ}$ in $U$. (Allen, L.J.S).

Note: The Hopf bifurcation requires at least a two dimensional differential equation system to appear.

Definition: The bifurcation stated in the Hopf bifurcation theorem is called ” supercritical” if the equilibrium point $\left(\mathrm{0,0}\right)$ is asymptotically stable when $Ϣ=0$ (at the bifurcation point) and it is called ” subcritical” if the equilibrium point $\left(\mathrm{0,0}\right)$ is negatively asymptotically stable (as $t\to -\mathrm{\propto })$ when $Ϣ=0\mathrm{.}$

In a supercritical Hopf bifurcation, the limit cycle grows out of the equilibrium point. In other words, right at the parameters of the Hopf bifurcation, the limit cycle has zero amplitude, and this amplitude grows as the parameters move further into the limit-cycle. (See the figure below)

Figure 5.

Bifurcation diagram corresponding to Supercritical Hopf bifurcation

Figure 6.

Supercritical Hopf bifurcation

However in a subcritical Hopf bifurcation, there is an unstable limit cycle surrounding the equilibrium point, and a stable limit cycle surrounding that. The unstable limit cycle shrinks down to the equilibrium point, which becomes unstable in the process. For systems started near the equilibrium point, the result is a sudden change in behavior from approach to a stable focus, to large-amplitude oscillations.(See the figure below).

Figure 7.

Diagram for Subcritical Hopf bifurcation

Figure 8.

Subcritical Hopf bifurcation

Example: Consider the two dimensional system

$$\frac{\mathit{dx}}{\mathit{dt}}=x^2-Ϗ,$$E9

$\frac{\mathit{dy}}{\mathit{dt}}=-(x^2+1)y,$

where $Ϗ$ is a parameter. When we compute the equilibrium points depending on parameter $Ϗ;$

If $Ϗ<0,$ then there is no x-nullclines, hence the system has no equilibrium points.If $Ϗ=0,$ then the system has exactly one equilibrium point at $\left(\mathrm{0,0}\right)\mathrm{.}$If $Ϗ>0,$ then the system has two equilibrium points $(-\sqrt{Ϗ},0)$ and $(\sqrt{Ϗ},0)\mathrm{.}$

Then the Jacobian matrix is

$J=\left[\begin{array}{ll}2x& 0\\ -2\mathit{xy}& -(x^2+1)\end{array}\right],$

where at the equilibrium points $J\left(0.0\right)=\left[\begin{array}{ll}0& 0\\ 0& -1\end{array}\right],J(\sqrt{Ϗ},0)=\left[\begin{array}{ll}-2\sqrt{Ϗ}& 0\\ 0& -Ϗ-1\end{array}\right]$ and $J(\sqrt{Ϗ},0)=\left[\begin{array}{ll}2\sqrt{Ϗ}& 0\\ 0& -Ϗ-1\end{array}\right]\mathrm{.}$

For $Ϗ=0$, there will be a line equilibrium (since one of the eigenvalues is zero) and for $Ϗ>0$, the point $(-\sqrt{Ϗ},0)$ is a sink and $(\sqrt{Ϗ},0)$ is a saddle point so that $Ϗ=0$ is the bifurcation point for this differential equation system.

Example: Consider the two dimensional system

$$\frac{\mathit{dx}}{\mathit{dt}}=y+Ϗx,$$E10

$\frac{\mathit{dy}}{\mathit{dt}}=-x+Ϗy,$

where $Ϗ$ is the bifurcation parameter. We can easily show that the conditions of the Hopf Bifurcation theorem hold. In this system $f_1$ and $f_2$ are zero. Then the Jacobian matrix is

$J=\left[\begin{array}{ll}Ϗ& 1\\ -1& Ϗ\end{array}\right]$

for which the eigenvalues are ${ϙ}_{\mathrm{1,2}}=Ϗ\pm i$ where Re$ϙ\left(Ϗ\right)=Ϗ$ and the imaginary part Im$ϙ\left(Ϗ\right)=\pm 1\mathrm{.}$It follows that Re$ϙ\left(0\right)=0$ and Im$ϙ\left(0\right)\stackrel{?}{=}0\mathrm{.}$ and also

$\frac{\mathit{dRe}ϙ\left(Ϗ\right)}{dϏ}{|}_{Ϗ=0}=1\stackrel{?}{=}0\mathrm{.}$

Hence, we conclude that there exists a periodic solution for $Ϗ=0$ in every neighborhood of origin.

4. Center Manifold teorem

Let $f\left(0\right)=0$, for the dynamical system

$$x^{\mathrm{\text{'}}}=f\left(x\right),x\nabla {\mathbb{R}}^n,$$E11

and let the eigenvalues of the Jacobian matrix be ${ϙ}_1,{ϙ}_2,\mathrm{...},{ϙ}_n\mathrm{.}$ Suppose that, the real parts of the eigenvalues are zero and if not, suppose there are $n_{+}$ numbers of eigenvalues with $\mathit{Re}ϙ>0,n_0$ number of eigenvalues with $\mathit{Re}ϙ=0$ and $n_{-}$ number of eigenvalues with $\mathit{Re}ϙ<0.$ Let $T^c$ be the eigenspace on imaginary axis corresponding to $n_0$ eigenvalues. The eigenvalues on the imaginary axis $\left(\mathit{Re}ϙ=0\right)$ are called the critical eigenvalues as on the eigenspace $T^c\mathrm{.}$ And suppose the function ${Ϥ}^t$ denote the flow corresponding to the equation $\left(11\right)$.

With these assumption, we state the Center Manifold theorem as follows;

Theorem: (Center Manifold theorem)

There exists a locally invariant $C^{\mathrm{\propto }}$ center manifold $W_{\mathit{loc}}^c\left(0\right)$ such that

$W_{\mathit{loc}}^c\left(0\right)=\left\{\right(x,y):y=h(x);|x|<\Upsilon ,h(0)=0;DJ(0)=0\}$

such that the dynamics of the system

$x^{\mathrm{\text{'}}}=A^{c}x+r_1(x,y),$$y^{\mathrm{\text{'}}}=A^{s}y+r_2(x,y),$

(where $A^c$and $A^s$ are are the blocks in the canonical form whose diagonals contain the eigenvalues with $\mathit{Re}ϙ=0$ and $\mathit{Re}ϙ<0$; respectively) restricted to the center manifold are given by

$x^{\mathrm{\text{'}}}=A^c\left(x\right)+r_1(x,h(x\left)\right)\mathrm{.}$

And the manifold $W_{\mathit{loc}}^c$ is called center manifold.

Remark: Center manifolds are not unique.

4.1. Center Manifold reduction for two dimensional systems

Example: Consider the two dimensional system of differential equations

$$x^{\mathrm{\text{'}}}=\mathit{xy},$$E12

$y^{\mathrm{\text{'}}}=-y+x^2\mathrm{.}$

The only equilibrium point is $\left(\mathrm{0,0}\right)$, we linearize around that and obtain

$J\left(\mathrm{0,0}\right)=\left[\begin{array}{ll}0& 0\\ 0& -1\end{array}\right]\mathrm{.}$

Now we look for $y=h\left(x\right)=a{x}^2+b{x}^3+c{x}^4+d{x}^5+O\left(x^6\right)\mathrm{.}$ Then,

$y^{\mathrm{\text{'}}}=h^{\mathrm{\text{'}}}\left(x\right)x^{\mathrm{\text{'}}}=x{h}^{\mathrm{\text{'}}}\left(x\right)h\left(x\right)$$=2a^2{x}^4+5\mathit{ab}{x}^5+O\left(x^6\right)$

and on the other hand,

$y^{\mathrm{\text{'}}}=-h\left(x\right)-x^2$$=-(a+1)x^2-b{x}^3-c{x}^4-d{x}^5+O\left(x^6\right)\mathrm{.}$

Comparing the two expressions, we deduce that $a=-1,b=0,c=-2,d=0$ and the center manifold reduction takes the form

$y=h\left(x\right)=-x^2-2x^4+O\left(x^7\right)\mathrm{.}$

Hence for the last equation $x=0$ is asymptotically stable and therefore $\left(\mathrm{0,0}\right)$ is asymptotically stable for the original system.

Example: Consider the two dimensional system

$$x^{\mathrm{\text{'}}}=x^{2}y-x^5,$$E13

$y^{\mathrm{\text{'}}}=-y+x^2\mathrm{.}$

Again $\left(\mathrm{0,0}\right)$ is an equilibrium point and the jacobian matrix for the linearized system is

$J\left(\mathrm{0,0}\right)=\left[\begin{array}{ll}0& 0\\ 0& -1\end{array}\right]\mathrm{.}$

Consider the transformation $y=h\left(x\right)=a{x}^2+b{x}^3+c{x}^4+d{x}^5+O\left(x^6\right),$ which leads

$y^{\mathrm{\text{'}}}=h^{\mathrm{\text{'}}}\left(x\right)x^{\mathrm{\text{'}}}$$=2a^2{x}^5+\left[2a\right(b-1)+3\mathit{ab}]x^6+O\left(x^7\right),$and $y^{\mathrm{\text{'}}}=-h\left(x\right)+x^2$$=-(a+1)x^2-b{x}^3+O\left(x^4\right)$$=2a^2{x}^4+5a-\mathit{bxb}{x}^5+O\left(x^6\right),$

from which we deduce that $a=1$ and $b=0$ and

$x^{\mathrm{\text{'}}}=x^4+O\left(x^5\right)\mathrm{.}$

For this reduced equation, $x=0$ is unstable and hence, $\left(\mathrm{0,0}\right)$ is also unstable for the original system.

4.2. Center Manifold reduction for Hopf bifurcation

The aim of this section is to give a formal framework for the analytical bifurcation analysis of Hopf bifurcations in delay differential equations

$$x^{\mathrm{\text{'}}}=f(x,Ϣ),xϓ{\mathbb{R}}^3$$E14

with a single fixed time delay $Ϣ$ to be chosen as a bifurcation parameter. Characteristic equations of the delay differential equation form (14) are often studied in order to understand changes in the local stability of equilibria of certain delay differential equations. It is therefore important to determine the values of the delay at which there are roots with zero part. We give a general formalization of these calculations and determine closed form algebraic equations where the stability and amplitude of periodic solutions close to bifurcation can be calculated.

We shall determine the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions by applying the normal form theory and the center manifold theorem by Hassard et al.,[10], and throughout this section, we assume that the three dimensional system of delay differential equations (14) undergoes Hopf bifurcations at the positive equilibrium $(N_0^{\mathrm{*}},P_0^{\mathrm{*}},S_0^{\mathrm{*}})$ at $Ϣ={Ϣ}_k,$ and $i{ϧ}_1$ is the corresponding purely imaginary root of the characteristic equation at the positive equilibrium $(N_0^{\mathrm{*}},P_0^{\mathrm{*}},S_0^{\mathrm{*}})$. For the sake of simplicity, we use the notation $iϧ$ for $i{ϧ}_1$.

We first consider the system (14) by the transformation

$x_1=N-N_0^{\mathrm{*}},x_2=P-P_0^{\mathrm{*}},x_3=S-S_0^{\mathrm{*}},t=\frac{t}{Ϣ},Ϣ={Ϣ}_k+Ϛ$

which is equivalent to the following Functional Differential Equation(FDE) system in $C=C\left(\right[-\mathrm{1,0}],{\mathbb{R}}^3)$,

$$\stackrel{\text{̥}}{x}\left(t\right)=L_{Ϛ}\left(x_t\right)+f(Ϛ,x_t),$$E15

where $x\left(t\right)=\left(x_1\right(t),x_2(t),x_3(t){)}^Tϓ{\mathbb{R}}^3,$ and $L_{Ϛ}:C\to {\mathbb{R}}^3,f:\mathbb{R}\times C\to {\mathbb{R}}^3$ are given respectively, by

$L_{Ϛ}\left(ϳ\right)=({Ϣ}_k+Ϛ)\left[\begin{array}{lll}{a}_1& a_2& a_3\\ a_4& a_5& a_6\\ a_7& a_8& a_9\end{array}\right]\left[\begin{array}{l}{ϳ}_1\left(0\right)\\ {ϳ}_2\left(0\right)\\ {ϳ}_3\left(0\right)\end{array}\right]+({Ϣ}_k+Ϛ)\left[\begin{array}{lll}{b}_1& b_2& b_3\\ b_4& b_5& b_6\\ b_7& b_8& b_9\end{array}\right]\left[\begin{array}{l}{ϳ}_1(-1)\\ {ϳ}_2(-1)\\ {ϳ}_3(-1)\end{array}\right],$and$f(Ϛ,ϳ)=({Ϣ}_k+Ϛ)\left[\begin{array}{l}{f}_{11}\\ f_{12}\\ f_{13}\end{array}\right]\mathrm{.}$

By Riesz representation theorem, there exists a function $\varphi (\varpi ,Ϛ)$ of bounded variation for $\mathit{\varpi ϓ}-\mathrm{1,0}],$ such that

$L_{Ϛ}ϳ={\cup }_{-}^1\mathrm{‎}d\varphi (\varpi ,0)ϳ\left(\varpi \right),\mathrm{for}ϳϓC\mathrm{.}$

Indeed we may take

$\varphi (\varpi ,Ϛ)=({Ϣ}_k+Ϛ)\left[\begin{array}{lll}{a}_1& a_2& a_3\\ a_4& a_5& a_6\\ a_7& a_8& a_9\end{array}\right]\Upsilon \left(\varpi \right)-({Ϣ}_k+Ϛ)\left[\begin{array}{lll}{b}_1& b_2& b_3\\ b_4& b_5& b_6\\ b_7& b_8& b_9\end{array}\right]\Upsilon (\varpi +1),$

where $\Upsilon$ is the Dirac delta function. For $ϳϓC^1\left(\right[-\mathrm{1,0}],{\mathbb{R}}^3),$ define

$A\left(Ϛ\right)ϳ=\left\{\begin{array}{l}\frac{\mathit{dϳ}\left(\varpi \right)}{d\varpi }\mathit{\varpi ϓ}-\mathrm{1,0}),\\ {\cup }_{-}^1\mathrm{‎}d\varphi (Ϛ,s)ϳ\left(s\right)\mathrm{}\varpi =0\end{array}\right.$and$R\left(Ϛ\right)ϳ=\left\{\begin{array}{l}0\mathrm{}\mathit{\varpi ϓ}-\mathrm{1,0}),\\ f(Ϛ,ϳ),\varpi =0\mathrm{.}\end{array}\right.$

Then the system (15) is equivalent to

${\stackrel{\text{̥}}{x}}_t=A\left(Ϛ\right)x_t+R\left(Ϛ\right)x_t,$

where $x_t\left(\varpi \right)=x(t+\varpi )$ for $\mathit{\varpi ϓ}-\mathrm{1,0})\mathrm{.}$ For $\mathit{Ϧϓ}{C}^1\left(\right[-\mathrm{1,0}],({\mathbb{R}}^3{)}^{\mathrm{*}}),$ define

$A^{\mathrm{*}}Ϧ\left(s\right)=\left\{\begin{array}{l}-\frac{dϦ\left(s\right)}{\mathit{ds}}sϓ\left(\mathrm{0,1}\right],\\ {\cup }_{-}^1\mathrm{‎}d{\varphi }^T(t,0)Ϧ(-t)\mathrm{}s=0\end{array}\right.$

and a bilinear inner product

$$⟢Ϧ\left(s\right),ϳ\left(\varpi \right)⟣=\stackrel{\text{̣}}{Ϧ}\left(0\right)ϳ\left(0\right)-{\cup }_{-}^1\mathrm{‎}{\cup }_{Ϝ}^{=}\mathrm{‎}\stackrel{\text{̣}}{Ϧ}(Ϝ-\varpi )d\varphi \left(\varpi \right)ϳ\left(Ϝ\right)dϜ,$$E16

where $\varphi \left(\varpi \right)=\varphi (\varpi ,0)\mathrm{.}$ Then $A\left(0\right)$ and $A^{\mathrm{*}}$ are adjoint operators. Suppose that $q\left(\varpi \right)$ and $q^{\mathrm{*}}\left(s\right)$ are eigenvectors of $A$ and $A^{\mathrm{*}}$ corresponding to $iϧ{Ϣ}_k$ and $-iϧ{Ϣ}_k,$ respectively. Then suppose that $q\left(\varpi \right)=(1,ϐ,\vartheta {)}^T{e}^{iϧ{Ϣ}_k\varpi }$ is the eigenvector of $A\left(0\right)$ corresponding to $iϧ{Ϣ}_k,$ then $A\left(0\right)q\left(\varpi \right)=iϧ{Ϣ}_{k}q\left(\varpi \right)\mathrm{.}$ Then in the following, we use the theory by Hassard et al.,[10], to compute the coordinates describing center manifold $C_{0\mathrm{}}$at $Ϛ=0\mathrm{.}$ Define

$$z\left(t\right)=⟢q^{\mathrm{*}},x_t⟣,\mathrm{}W(t,\varpi )=x_t-2\mathit{Rez}\left(t\right)q\left(\varpi \right)\mathrm{.}$$E17

On the center manifold $C_0$, we have

$W(t,\varpi )=W\left(z\right(t),\stackrel{\text{̣}}{z}(t),\varpi )=W_{20}\left(\varpi \right)\frac{z^2}{2}+W_{11}\left(\varpi \right)z\stackrel{\text{̣}}{z}+W_{02}\left(\varpi \right)\frac{{\stackrel{\text{̣}}{z}}^2}{2}+\mathrm{...},$

where $z$ and $\stackrel{\text{̣}}{z}$ are local coordinates for center manifold $C_{0\mathrm{}}$in the direction of $q$ and ${\stackrel{\text{̣}}{q}}^{\mathrm{*}}\mathrm{.}$ Note that $W$ is real if $x_t$ is real. We consider only real solutions. For the solution $x_tϓC_0,$ since $Ϛ=0$ and (15), we have

$\stackrel{\text{̥}}{z}=iϧ{Ϣ}_{k}z+⟢q^{\mathrm{*}}\left(\varpi \right),F(0,W(z,\stackrel{\text{̣}}{z},\varpi )+2\mathit{Rezq}(\varpi \left)\right)⟣$$=iϧ{Ϣ}_{k}z+{\stackrel{\text{̣}}{q}}^{\mathrm{*}}\left(0\right)F(0,W(z,\stackrel{\text{̣}}{z},0)+2\mathit{Rezq}(0\left)\right)$$\stackrel{\mathit{def}}{=}iϧ{Ϣ}_{k}z+{\stackrel{\text{̣}}{q}}^{\mathrm{*}}\left(0\right)F_0(z,\stackrel{\text{̣}}{z})$$=iϧ{Ϣ}_{k}z+g(z,\stackrel{\text{̣}}{z}),$

where

$$g(z,\stackrel{\text{̣}}{z})={\stackrel{\text{̣}}{q}}^{\mathrm{*}}\left(0\right)F_0(z,\stackrel{\text{̣}}{z})=g_{20}\frac{z^2}{2}+g_{11}z\stackrel{\text{̣}}{z}+g_{02}\frac{{\stackrel{\text{̣}}{z}}^2}{2}+g_{21}\frac{z^2\stackrel{\text{̣}}{z}}{2}+\mathrm{....}$$E18

By using (17), we have $x_t\left(x_{1t}\right(\varpi ),x_{2t}(\varpi ),x_{3t}(\varpi \left)\right)=W(t,\varpi )+\mathit{zq}\left(\varpi \right)+\overline{\mathit{zq}}\left(\varpi \right)$, $q\left(\varpi \right)=(1,ϐ,\vartheta {)}^T{e}^{iϧ{Ϣ}_k\varpi },$ and

$x_{1t}\left(0\right)=z+\stackrel{\text{̣}}{z}+W_{20}^{\left(1\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(0\right)z\stackrel{\text{̣}}{z}+W_{02}^{\left(1\right)}\left(0\right)\frac{{\stackrel{\text{̣}}{z}}^2}{2}+O\left({\left|(z,\stackrel{\text{̣}}{z})\right|}^3\right),$$x_{2t}\left(0\right)=ϐz+\stackrel{\text{̣}}{ϐ}\overline{z}+W_{20}^{\left(2\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(2\right)}\left(0\right)z\stackrel{\text{̣}}{z}+W_{02}^{\left(2\right)}\left(0\right)\frac{{\stackrel{\text{̣}}{z}}^2}{2}+O\left({\left|(z,\stackrel{\text{̣}}{z})\right|}^3\right),$$x_{3t}\left(0\right)=\vartheta z+\overline{\vartheta }\overline{Ϗ}+W_{20}^{\left(3\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(3\right)}\left(0\right)z\stackrel{\text{̣}}{z}+W_{02}^{\left(3\right)}\left(0\right)\frac{{\stackrel{\text{̣}}{z}}^2}{2}+O\left({\left|(z,\stackrel{\text{̣}}{z})\right|}^3\right),$$x_{1t}(-1)=z{e}^{-iϧ{Ϣ}_k\varpi }+\stackrel{\text{̣}}{z}{e}^{iϧ{Ϣ}_k\varpi }+W_{20}^{\left(1\right)}(-1)\frac{z^2}{2}+W_{11}^{\left(1\right)}(-1)z\stackrel{\text{̣}}{z}+W_{02}^{\left(1\right)}(-1)\frac{{\stackrel{\text{̣}}{z}}^2}{2}$$+O\left({\left|(z,\stackrel{\text{̣}}{z})\right|}^3\right)\mathrm{.}$

From the definition of $F(Ϛ,x_t)$, we have

$g(z,\stackrel{\text{̣}}{z})={\stackrel{\text{̣}}{q}}^{\mathrm{*}}\left(0\right)f_0(z,\stackrel{\text{̣}}{z})=\stackrel{\text{̣}}{D}{Ϣ}_k(1,{\stackrel{\text{̣}}{ϐ}}^{\mathrm{*}},{\stackrel{\text{̣}}{\vartheta }}^{\mathrm{*}})\left[\begin{array}{l}{f}_{11}^0\\ f_{12}^0\\ f_{13}^0\end{array}\right]$

and we evaluate $g(z,\stackrel{\text{̣}}{z})\mathrm{.}$

To determine $g_{21}$, we need to compute $W_{20}\left(\varpi \right)$ and $W_{11}\left(\varpi \right)$. By (15) and (18), we have

$$\stackrel{\text{̥}}{W}={\stackrel{\text{̥}}{x}}_t-\stackrel{\text{̥}}{z}q+\stackrel{\mathrm{.}}{\overline{\mathit{zq}}}$$E19

$=\{\begin{array}{l}\mathit{AW}-2\mathit{Re}\left\{{\stackrel{\text{̣}}{q}}^{\mathrm{*}}\right(0\left)f_{0}q\right(\varpi \left)\right\},\mathrm{}\mathit{\varpi ϓ}-\mathrm{1,0})\\ \mathit{AW}-2\mathit{Re}\left\{{\stackrel{\text{̣}}{q}}^{\mathrm{*}}\right(0\left)f_{0}q\right(\varpi \left)\right\}+f_0\varpi =0,\mathrm{}\end{array}$$\stackrel{\mathit{def}}{=}\mathit{AW}+H(z\mathrm{.}\stackrel{\text{̣}}{z},\varpi ),$

where

$$H(z\mathrm{.}\stackrel{\text{̣}}{z},\varpi )=H_{20}\left(\varpi \right)\frac{z^2}{2}+H_{11}\left(\varpi \right)z\stackrel{\text{̣}}{z}+H_{02}\left(\varpi \right)\frac{{\stackrel{\text{̣}}{z}}^2}{2}+\mathrm{....}$$E20

Note that on the center manifold $C_0$ near to the origin,

$$\stackrel{\text{̥}}{W}=W_z\stackrel{\text{̥}}{z}+W_{\stackrel{\text{̣}}{z}}\stackrel{\text{̣}}{z}\mathrm{.}$$E21

Thus we obtain,

$$(A-2iϧ{Ϣ}_k)W_{20}\left(\varpi \right)=-H_{20}\left(\varpi \right),\mathrm{}A{W}_{11}\left(\varpi \right)=-H_{11}\left(\varpi \right)\mathrm{.}$$E22

By using (19), for $\mathit{\varpi ϓ}-\mathrm{1,0}),$

$$H(z\mathrm{.}\stackrel{\text{̣}}{z},\varpi )={\stackrel{\text{̣}}{q}}^{\mathrm{*}}\left(0\right)f_{0}q\left(\varpi \right)-q^{\mathrm{*}}\left(0\right)f_0\left(0\right)\stackrel{\text{̣}}{q}\left(\varpi \right)=-\mathit{gq}\left(\varpi \right)-\stackrel{\text{̣}}{g}\stackrel{\text{̣}}{q}\left(\varpi \right)\mathrm{.}$$E23

Comparing the coefficients with (20), we obtain the following

$$H_{20}\left(\varpi \right)=-g_{20}q\left(\varpi \right)-{\stackrel{\text{̣}}{g}}_{02}\stackrel{\text{̣}}{q}\left(\varpi \right),\mathrm{}{H}_{11}\left(\varpi \right)=-g_{11}q\left(\varpi \right)-{\stackrel{\text{̣}}{g}}_{11}\stackrel{\text{̣}}{q}\left(\varpi \right)\mathrm{.}$$E24

From (22) and (24) and the definition of $A,$ we get

${\stackrel{\text{̥}}{W}}_{20}\left(\varpi \right)=2iϧ{Ϣ}_k{W}_{20}\left(\varpi \right)-g_{20}q\left(\varpi \right)-{\stackrel{\text{̣}}{g}}_{02}\stackrel{\text{̣}}{q}\left(\varpi \right)\mathrm{.}$

Noticing $q\left(\varpi \right)=q\left(0\right)e^{iϧ{Ϣ}_k\varpi },$ we evaluate $W_{20}\left(\varpi \right)$ by $E_1=(E_1^{\left(1\right)},E_1^{\left(2\right)},E_1^{\left(3\right)})ϓ{\mathbb{R}}^3$ is a constant vector. From the definition of $A$ and (22), we obtain

$${\cup }_{-}^1\mathrm{‎}d\varphi \left(\varpi \right)W_{20}\left(\varpi \right)=2iϧ{Ϣ}_k{W}_{20}\left(0\right)-H_{20}\left(0\right),$$E25

and

$${\cup }_{-}^1\mathrm{‎}d\varphi \left(\varpi \right)W_{11}\left(\varpi \right)=-H_{11}\left(0\right),$$E26

where $d\varphi \left(\varpi \right)=\varphi (\varpi ,0)\mathrm{.}$ Next we compute $W_{20}\left(\varpi \right)$ and $W_{11}\left(\varpi \right)$ from (25) and (26) and determine the following values to investigate the qualities of bifurcating periodic solution in the center manifold at the critical value ${Ϣ}_k\mathrm{.}$ For this purpose, we express the direction of Hopf bifurcation in terms of $g_{\mathit{ij}}^{\mathrm{\text{'}}}s$ and eigenvalues $ϙ{Ϣ}_k$. And then we can evaluate the following values;

$c_1\left(0\right)=\frac{i}{2ϧ{Ϣ}_k}(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{|g_{02}{|}^2}{3})+\frac{g_{21}}{2},{Ϛ}_2=-\frac{\mathit{Re}\left\{c_1\right(0\left)\right\}}{\mathit{Re}\left\{{ϙ}^{{}^{\mathrm{\text{'}}}}\right({Ϣ}_k\left)\right\}},$

$${ϐ}_2=2\mathit{Re}\left\{c_1\right(0\left)\right\},T_2=-\frac{\mathit{Im}\left\{c_1\right(0\left)\right\}+{Ϛ}_2\mathit{Im}\left\{{ϙ}^{{}^{\mathrm{\text{'}}}}\right({Ϣ}_k\left)\right\}}{ϧ{Ϣ}_k}$$E27

and we state this as in the following theorem.

Theorem: ${Ϛ}_2$ determines the direction of Hopf bifurcation; if ${Ϛ}_2>0,$ then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exist for $Ϣ>{Ϣ}_0$, if ${Ϛ}_2<0,$ then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exist for $Ϣ<{Ϣ}_0$. ${ϐ}_2$ determines the stability of the bifurcating periodic solutions; bifurcating periodic solutions are stable if ${ϐ}_2<0,$ unstable if ${ϐ}_2>0$. $T_2$ determines the period of the bifurcating solution; the period increases if $T_2>0,$ it decreases if $T_2<0.$

In the following section, we shall give a numerical example to verify the theoretical results.

4.2.1. Numerical example of Center Manifold reduction

Consider the following system with discrete time delay $Ϣ$;

$$\begin{array}{l}\frac{\mathit{dN}\left(t\right)}{\mathit{dt}}=r_{1}N\left(t\right)-ϓP\left(t\right)N\left(t\right),\\ \frac{\mathit{dP}\left(t\right)}{\mathit{dt}}=P\left(t\right)(r_2-\varpi \frac{S\left(t\right)}{N(t-Ϣ)}),\\ \frac{\mathit{dS}\left(t\right)}{\mathit{dt}}=ϏP\left(t\right)-ϏS\left(t\right)\mathrm{.}\end{array}$$E28

Now we present some numerical simulations by using MATLAB(7.6.0) programming (Çelik-3). We simulate the predator-prey system (28) by choosing the parameters $r_1=0.45$, $r_2=0.1$, $\varpi =0.05$, $ϓ=0.03$, and $Ϗ=1$, i.e., we consider the following system

$$\begin{array}{l}\frac{\mathit{dN}\left(t\right)}{\mathit{dt}}=0.45N\left(t\right)-0.03P\left(t\right)N\left(t\right),\\ \frac{\mathit{dP}\left(t\right)}{\mathit{dt}}=P\left(t\right)(0.1-0.05\frac{S\left(t\right)}{N(t-Ϣ)}),\\ \frac{\mathit{dS}\left(t\right)}{\mathit{dt}}=P\left(t\right)-S\left(t\right),\end{array}$$E29

which has only one positive equilibrium $E^{\mathrm{*}}=(N_0^{\mathrm{*}},P_0^{\mathrm{*}},S_0^{\mathrm{*}})=\left(\mathrm{7.5,15,15}\right)$. By algorithms in the previous sections, we obtain ${Ϣ}_0=1.5663$, ${ϧ}_1=0.0045$,$z_1={ϧ}_1^2$ and $g^{\mathrm{\text{'}}}\left(z_1\right)=1.1944\times 10^{-4}>0$ which leads to $\frac{\mathit{dRe}ϙ\left({Ϣ}_0\right)}{dϢ}=5.8982>0$. So by the theorem above, the equilibrium point $E^{\mathrm{*}}$ is asymptotically stable when $\mathit{Ϣϓ}[0,{Ϣ}_0)=\left[\mathrm{0,1.5663}\right)$ and unstable when $Ϣ>1.5663$ and also Hopf bifurcation occurs at $Ϣ={Ϣ}_0=1.5663$ as it is illustrated by computer simulations.

By the theory of Hassard et al.,[10], as it is discussed in previous section, we also determine the direction of Hopf bifurcation and the other properties of bifurcating periodic solutions. From the formulae in Section 5.2 we evaluate the values of ${Ϛ}_2$, ${ϐ}_2$ and $T_2$ as

${Ϛ}_2=0.0981>0,\mathit{}{ϐ}_2=-0.5785<0,\mathit{}{T}_2=7.1165>0,$

from which we conclude that Hopf bifurcation of system (29) occurring at ${Ϣ}_0=1.5663$ is supercritical and the bifurcating periodic solution exists when $Ϣ$ crosses ${Ϣ}_0$to the right, and also the bifurcating periodic solution is stable.

In computer simulations, the initial conditions are taken as $(N_0,P_0,S_0)=\left(\mathrm{50,25,25}\right)$ and MATLAB DDE (Delay Differential Equations) solver is used to simulate the system (29). We first take $Ϣ=1.5<{Ϣ}_0$ and plot the density functions $N\left(t\right)$, $P\left(t\right)$ and $S\left(t\right)$ in Figs.9,10,11 respectively which shows the positive equilibrium is asymptotically stable for $Ϣ<{Ϣ}_{0\mathrm{.}}$. Moreover in Fig.12, we illustrate the asymptotic stability in three dimension.

Figure 9.

The trajectory of prey density versus time with the initial condition $N_0=50$. The graph of solutions of the model (29) when $Ϣ=1.5<{Ϣ}_0$, where the equilibrium point $E^{\mathrm{*}}$ is asymptotically stable.

Figure 10.

The trajectory of predator density versus time with the initial condition $P_0=25$. The graph of solutions of the model (29) when $Ϣ=1.5<{Ϣ}_0$, where the equilibrium point $E^{\mathrm{*}}$ is asymptotically stable.

Figure 11.

The trajectories of $S\left(t\right)$ versus time with the initial condition $S_0=25$. The graph of solutions of the model (29) when $Ϣ=1.5<{Ϣ}_0$, where the equilibrium point $E^{\mathrm{*}}$ is asymptotically stable.

Figure 12.

The trajectory of $N,P,S$ in three dimension with the initial condition $\left(\mathrm{50,25,25}\right)$ when $Ϣ=1.5<{Ϣ}_0$ for which the equilibrium point $E^{\mathrm{*}}$ is asymptotically stable.

Figure 13.

The trajectory of $N\left(t\right)$ when $Ϣ=2>{Ϣ}_0$.

Figure 14.

The trajectory of $P\left(t\right)$ when $Ϣ=2>{Ϣ}_0$.

Figure 15.

The trajectory of $S\left(t\right)$ when $Ϣ=2>{Ϣ}_0$.

Figure 16.

When $Ϣ=2>{Ϣ}_0$ a stable periodic orbit bifurcates from the equilibrium point $E^{\mathrm{*}}$.

However in Figs.13,14,15 and 16 below, we take $Ϣ=2>{Ϣ}_0$ sufficiently close to ${Ϣ}_0$ which illustrates the existence of bifurcating periodic solutions from the equilibrium point $E^{\mathrm{*}}$.

4.3. Exercises

1. Consider the Van der Pol equation

$\stackrel{\text{̦}}{x}-(1-x^2)\stackrel{\text{̥}}{x}+x=0$

and convert this into a system and check the stability of fixed points.

2. Let $y_0$ be an equilibrium point of the equation

$\frac{\mathit{dy}}{\mathit{dt}}=f\left(y\right)\mathrm{.}$Let $g\left(y\right)$ be a positive function. Deduce the stability of $y_0$ as an equilibrium solution of the equation$\frac{\mathit{dy}}{\mathit{dt}}=f\left(y\right)g\left(y\right)$

from its stability as an equilibrium solution of the equation $\frac{\mathit{dy}}{\mathit{dt}}=f\left(y\right)\mathrm{.}$ Repeat the same question if $g\left(y\right)$is a negative function.

3. For the following nonlinear system

$x^{\mathrm{\text{'}}}=-x^3-y,$$y^{\mathrm{\text{'}}}=x-y^3,$

determine the stability of fixed points.

4. For the following nonlinear system

$x^{\mathrm{\text{'}}}=y-x{y}^2,$$y^{\mathrm{\text{'}}}=-x+x^{2}y,$

classify the fixed points.

5. Analyze the bifurcation properties of the following problems choosing $r$ as bifurcation parameter,

1. $x^{\mathrm{\text{'}}}=-x+ϐ\tanh x,$
2. $x^{\mathrm{\text{'}}}=\mathit{rx}-4x^3,$
3. $x^{\mathrm{\text{'}}}=\mathit{rx}-\mathit{sin}\left(x\right),$
4. $x^{\mathrm{\text{'}}}=\mathit{rx}+4x^3,$
5. $x^{\mathrm{\text{'}}}=\mathit{rx}-\sin h\left(x\right),$
6. $x^{\mathrm{\text{'}}}=x+\frac{\mathit{rx}}{1+x^2}\mathrm{.}$

6. Find the equilibrium points and identify the bifurcation in the following system, and sketch the appropriate bifurcation diagram and phase portraits:

$\frac{\mathit{dx}}{\mathit{dt}}=(1+Ϣ)y-x-2x^3,$$\frac{\mathit{dy}}{\mathit{dt}}=x-y-y^3\mathrm{.}$

Then compute the extended center manifold near the bifurcation point by choosing $Ϣ$ as bifurcation parameter.

7. Show that the following system is structurally unstable,

$\frac{\mathit{dx}}{\mathit{dt}}=x+y-x^2+y^2,$$\frac{\mathit{dy}}{\mathit{dt}}=-2x-y+\mathit{xy},$

and the following system is structural stable. Explain our reason.

$\frac{\mathit{dx}}{\mathit{dt}}=2x+y-x(x^2+y^2),$$\frac{\mathit{dy}}{\mathit{dt}}=-x+2y-y(x^2+y^2)\mathrm{.}$

8.

$\frac{\mathit{dx}}{\mathit{dt}}=\mathit{xy},$$\frac{\mathit{dy}}{\mathit{dt}}=1-y^2-x^2,$

for the above system, classify the fixed points.

9. For $r\nabla R$, consider the differential equation

$x^{\mathrm{\text{'}}}=\mathit{rx}-2x^2+x^3$

on the real line.

1. Show that $x^{\mathrm{*}}=0\mathrm{}$ is a fixed point for any value of the parameter r, and determine its stability. Hence identify a bifurcation point $r_{1\mathrm{.}}$

2. Show that for certain values of the parameter $r$ there are additional fixed points. For which values of $r$ do these fixed points exist? Determine their stability and identify a further bifurcation points $r_2\mathrm{.}$

3. Using a Taylor expansion of the differential equation above, determine the normal form of the bifurcation at $r_1$. What type of bifurcation takes place.

4. Similarly, determine the normal from of the bifurcation at $r_2\mathrm{.}$What type bifurcation takes place?

5. Sketch the bifurcation diagram for all values of $r$ and $x^{\mathrm{*}}$. (Use a full line to denote a curve of stable fixed points, and a dashed line for a curve of unstable fixed points).

10. Consider the system of differential equations

$x^{\mathrm{\text{'}}}=y-x{y}^2,$$y^{\mathrm{\text{'}}}=-x+y{x}^2,$

in the $(x,y)$ plane.

1. Determine all fixed points of the system.

2. $r^2=x^2+y^2\mathrm{.}$$r^{\mathrm{\text{'}}}\mathrm{}=0$

3. Sketch the phase portrait in the $(x,y)\mathrm{}$plane, including trajectories through $\left(\mathrm{1,0}\right)\mathrm{}$and $\left(\mathrm{2,0}\right)\mathrm{.}$ Whichfixed point does the trajectory through $\left(\mathrm{2,0}\right)$ approach?

11. Consider the Lorenz system (the model of heat convection by Rayleigh-Benard occurring in the earth’s atmosphere) in three dimension as follows;

$x^{\mathrm{\text{'}}}=ϡ(y-x),$$y^{\mathrm{\text{'}}}=\mathit{rx}+-\mathit{xz}-y,$$z^{\mathrm{\text{'}}}=\mathit{xy}-\mathit{bz},$

where $ϡ,r,b$ are constants. P$\mathit{erf}$om the stability analysis of this nonlinear system.

12. For the Holling-Tanner type Predator-Prey model

$x^{\mathrm{\text{'}}}=x(1-x)-\frac{x}{a+x}y,$$y^{\mathrm{\text{'}}}=y(\Upsilon -ϐ\frac{y}{x}),$

where $x\left(0\right)>0,y\left(0\right)>0$ and $a,\Upsilon ,ϐ$ are positive constants. Find the equilibria and classify them.

13. Consider the Delayed Predator-Prey model

$x^{\mathrm{\text{'}}}\left(t\right)=\mathit{rx}\left(t\right)-\mathit{bx}\left(t\right)x(t-Ϣ)-Ϗx\left(t\right)y\left(t\right),$$y^{\mathrm{\text{'}}}\left(t\right)=-\mathit{cy}\left(t\right)+ϐx\left(t\right)y\left(t\right),$

with time delay $Ϣ$ and positive constants $r,b,Ϗ,c$ and $ϐ\mathrm{.}$ By choosing $Ϣ$ as bifurcation parameter, check whether bifurcating periodic solutions occur around the equilibrium points or not.