# Hopf Bifurcation: The Emergence of Limit-cycle Dynamics

*Cem Özen*, May 2017.

A *Hopf bifurcation* is a critical point in which a periodic orbit appears or disappears through a local change in the stability of a fixed point in a dynamical system as one of the system parameters is varied. Hopf bifurcations occur in many of the well-known dynamical systems such as the Lotka-Volterra model, the Lorenz model, the Selkov model of glycolysis, the Belousov-Zhabotinsky reaction model, and the Hodgkin-Huxley model for nerve membrane.

In this notebook, I will consider a system of chemical reactions known by the name *Brusselator* in literature (see: https://en.wikipedia.org/wiki/Brusselator for more information) as a model for Hopf bifurcations. The Brusselator reactions are given by

$A \longrightarrow X$ <br>
$2X + Y\longrightarrow 3X$ <br>
$B + X \longrightarrow Y + D$ <br>
$X \longrightarrow E$ <br>

For the sake of simplicity, we will assume that the reaction constants of all these reactions are unity (i.e. in all the reactions, $k=1$ ). Furthermore let's assume that the reactant concentrations $A$ and $B$ are so large that  they remain constant. Therefore, only $X$ and $Y$ concentrations will be dynamical. 

The rate equations for $X$ and $Y$ are then given by <br>


$$
\begin{eqnarray}
\dot{X} & = & A + X^2Y - BX - X, \\
\dot{Y} & = & BX - X^2Y
\end{eqnarray}
$$

The X-nullcline and the Y-nullcline are given by the conditions of  $0 =  A + X^2Y - BX - X$ and $0 =  BX - X^2Y$ respectively. From these equations, we obtain:

$$
\begin{eqnarray}
Y(X) & = & \frac{-A + X(B+1)}{X^2},  & \quad \textrm{(X-nullcline)}  \\
Y(X) & = & \frac{B}{X}, & \quad   \textrm{(Y-nullcline)} 
\end{eqnarray}
$$

In this notebook, I will also demonstrate how one can perform symbolical computations using Python's `SymPy` library. We also need extra Jupyter Notebook functionality to render nice display of the resulting equations. (Notice that we are using LaTex in typesetting this document particularly for the purpose of producing nice looking equations). 


```python
import numpy as np 
from numpy.linalg import eig
from scipy import integrate
import sympy
from IPython.display import display, Math, Latex
import matplotlib.pyplot as plt
sympy.init_printing(use_latex='mathjax')
%matplotlib inline
```

Let's obtain the nullcline equations using `SymPy`: 


```python
X, Y, A, B = sympy.symbols('X Y A B')   # you need to introduce the sysmbols first

# let's get the X-nullcline as a function of X:
sympy.solve(sympy.Eq(A + X**2 * Y - B * X - X, 0), Y)
```




$$\left [ \frac{1}{X^{2}} \left(- A + X \left(B + 1\right)\right)\right ]$$




```python
# let's get the Y-nullcline as a function of X:
sympy.solve(sympy.Eq(B * X - X**2 * Y, 0), Y)
```




$$\left [ \frac{B}{X}\right ]$$



Now let's find the fixed points ($X^*, Y^*$) of this 2-D system (there is only one, actually). The fixed point is given by the simultaneous solution of the X-nullcline and Y-nullcline equations, therefore 

$$ (X^*, Y^*)  =  (A, \frac{B}{A}) $$

For the sake of using `SymPy`, let's obtain this solution once again:


```python
# Solve the system of equations defined by the X-nullcline and Y-nullcline with respect to X and Y: 
sympy.solve([A + X**2 * Y - B * X - X,  B * X - X**2 * Y], [X, Y])
```




$$\left [ \left ( A, \quad \frac{B}{A}\right )\right ]$$



Now, a bifurcation analysis of the Brusselator model requires us to keep track of the local stability of its fixed point. This can be done, according to the *linearized stability analysis* by evaluating the Jacobian  matrix at the fixed point.  <br>


The Jacobian matrix at the fixed point is given by:

$$
\begin{eqnarray}
J & = &  \left\vert\matrix{{\partial f \over \partial x} & {\partial f\over \partial y}  \cr 
                           {\partial g \over \partial x} & {\partial g\over \partial y}
                           }\right\vert_{(X^*, Y^*)} \\
  & = & \left\vert\matrix{ -B + 2XY - 1 & X^2  \cr 
                            B - 2XY & -X^2
                          }\right\vert_{(X^*, Y^*)}  \\
  & = & \left\vert\matrix{ B - 1 & A^2  \cr 
                           -B & -A^2
                          }\right\vert
\end{eqnarray}
$$

This result can also be obtained very easily using `SymPy`:


```python
# define the Brusselator dynamical system as a SymPy matrix
brusselator = sympy.Matrix([A + X**2 * Y - B * X - X, B * X - X**2 * Y])
```


```python
# Jacobian matrix with respect to X and Y
J = brusselator.jacobian([X, Y])
J
```




$$\left[\begin{matrix}- B + 2 X Y - 1 & X^{2}\\B - 2 X Y & - X^{2}\end{matrix}\right]$$




```python
# Jacobian matrix evaluated at the coordinates of the fixed point
J_at_fp = J.subs({X:A, Y:B/A})   # substitute X with A and Y with B/A
J_at_fp
```




$$\left[\begin{matrix}B - 1 & A^{2}\\- B & - A^{2}\end{matrix}\right]$$



A limit-cycle can emerge in a 2-dimensional, attractive dynamical system if the fixed point of the system  goes unstable. In this case, trajectories must be pulled by a limit cycle. (According to the Poincare-Bendixson theorem, a 2-dimensional system cannot have strange attractors). In this case, the Hopf bifurcation is called a *supercritical Hopf bifurcation*, because limit cycle is stable. 

In the following, we will see how the stable fixed point (spiral) of the Brusselator goes unstable, giving rise to a limit cycle in turn. Conditions for the stability are determined by the trace and the determinant of the Jacobian. So let's evaluate them:


```python
Delta = J_at_fp.det()  # determinant of the Jacobian
Delta.simplify()
```




$$A^{2}$$




```python
tau = J_at_fp.trace() # trace of the Jacobian
tau
```




$$- A^{2} + B - 1$$



To have an unstable spiral we need:

$$
\begin{eqnarray}
\tau & > & 0  \quad \Rightarrow \quad & B > A^2 + 1  \quad  \textrm{required} \\
\Delta & > & 0  \quad {} \quad & \textrm{automatically satisfied} \\
\tau^2 & < & 4 \Delta \quad {} \quad & \textrm{automatically satisfied}
\end{eqnarray}
$$

The second and third conditions were satisfied because of the first condition, automatically. Therefore we need to have:

$$ B > A^2 + 1 $$ for limit cycles.

## Birth of A Limit Cycle: Hopf Bifurcation
### Numerical Simulation of the Brusselator System
In the following, I perform a numerical simulation of the (supercritical) Hopf bifurcation in the Brusselator system by varying the parameter $B$ while the value of $A=1$. 


```python
# Brusselator System:
def dX_dt(A, B, X, t):
    x, y = X[0], X[1]
    return np.array([A + x**2 * y - B * x -x, 
                     B * x - x**2 * y])

T = 50 * np.pi                      # simulation time
dt = 0.01                           # integration time step

# time steps to be used in integration of the Brusselator system
t=np.arange(0, T, dt)    

# create a canvas and 3 subplots..we will use each one for different choice of A and B paraeters
fig, ax = plt.subplots(1, 3)
fig.set_size_inches(13,5)

def plotter(A, B, ax):
    """
    This function draws a phase portrait by assigning a vector characterizing how the concentrations
    change at a given value of X and Y. It also draws a couple of example trajectories.
    """
    # Drow direction fields using matplotlib 's quiver function, similar to what we did
    # in class but qualitatively
    xmin, xmax  = 0,  5                           # min and max values of x axis in the plot
    ymin, ymax  = 0,  5                           # min and max values of y axis in the plot
    x = np.linspace(xmin, xmax, 10)               # divide x axis to intervals
    y = np.linspace(ymin, ymax, 10)               # divide y axis to intervals
    X1 , Y1  = np.meshgrid(x, y)            # from these intervals create a grid
    DX1, DY1 = dX_dt(A, B, [X1, Y1], t)     # compute rate of change of the concentrations on grid points
    M = (np.hypot(DX1, DY1))                # norm of the rate of changes 
    M[ M == 0] = 1.                         # prevention against divisions by zero
    DX1 /= M                                # we normalize the direction field vector (each has unit length now)
    DY1 /= M                                # we normalize the direction field vector (each has unit length now)
    Q = ax.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=plt.cm.jet)

    num_traj = 5  # number of trajectories 

    # choose several initial points (x_i, y_j),  for i and j chosen as in linspace calls below
    X0  = np.asarray(list(zip(np.linspace(xmin, xmax, num_traj), np.linspace(ymin, ymax, num_traj))))                         
    vcolors = plt.cm.jet_r(np.linspace(0., 1., num_traj))  # colors for each trajectory

    # integrate the Brusellator ODE's using all initial points to produce corresponding trajectories
    X = np.asarray([integrate.odeint(lambda x, t: dX_dt(A, B, x, t), X0i,t)
                          for X0i in X0])
    # plot the trajectories we obtained above 
    for i in range(num_traj):
        x, y = X[i, :, :].T          # x and y histories for trajectory i
        ax.plot(x, y, '-', c=vcolors[i], lw=2)
    # set limits, put labels etc..    
    ax.set_xlim(xmin=xmin, xmax=xmax)
    ax.set_ylim(ymin=ymin, ymax=ymax)
    ax.set_xlabel("X", fontsize = 20)
    ax.set_ylabel("Y", fontsize = 20)
    ax.annotate("A={}, B={}".format(A, B), xy = (0.4, 0.9), xycoords = 'axes fraction', fontsize = 20, color = "k")

    
# Now let's prepare plots for the following choices of A and B:    
plotter(A=1, B=1, ax=ax[0])
plotter(A=1, B=2, ax=ax[1])
plotter(A=1, B=3, ax=ax[2])
```

Above you see how a limit cycle can be created in a dynamical system, as one of the system parameters is varied.
Here we have kept $A=1$ but varied $B$ from 1 to 2.5. Note that $B=2$ is the borderline value, marking the change in the stability of the fixed point. For $B<1$ the fixed point is stable but as we cross the value $B=2$, it changes its character to unstable and a limit cycle is born. This phenomenon is an example of *Hopf Bifurcation*.

On the leftmost panel we have a stable spiral. Technically, this means that the Jacobian at the fixed point has two complex eigenvalues (a complex conjugate pair). The fact that the eigenvalues are complex is responsible for the spiralling effect. In stable spirals, the real part of the eigenvalues are negative, which is why these spiralling solutions decay, that is, trajectories nearby fall on to the fixed point. As the bifurcation parameter (here $B$) varies, the real part of the complex eigenvalues increase, reach zero at certain value of $B$, and keep growing now on the positive side. If the real part of the eigenvalues are positive, the fixed point is  an unstable spiral; trajectories nearby are pushed out of the fixed point (see the rightmost panel and plot below). Since this 2-D dynamical system is attractive, by the Poincare-Bendixson theorem, the emergence of the unstable spiral accompanies the birth of a limit-cycle. Notice that the panel in the middle is the borderline case between the stable and unstable spirals: in this case the real part of the eigenvalues is exactly zero (see plots below); the linear stabilization theory falsely predicts a neutral oscillation (i.e a center) at $B=2$---due to purely imaginary eigenvalues. However, the fixed point is still a stable spiral then.

### Eigenvalues of the Jacobian 


```python
# Eigenvalues of the Jacobian at A=1, B=1  (fixed point is stable spiral)
J_numeric = np.asarray(J_at_fp.evalf(subs={A:1, B:1})).astype(np.float64)
w, _ = eig(J_numeric)
w
```




    array([-0.5+0.8660254j, -0.5-0.8660254j])




```python
# Eigenvalues of the Jacobian at A=1, B=3  (fixed point is unstable spiral)
J_numeric = np.asarray(J_at_fp.evalf(subs={A:1, B:3})).astype(np.float64)
w, _ = eig(J_numeric)
w
```




    array([ 0.5+0.8660254j,  0.5-0.8660254j])



Let's prepare plots showing how the real and imaginary parts of the eigenvalues change as $B$ is varied.


```python
from numpy.linalg import eig

a = 1
eigen_real, eigen_imag = [], []
B_vals = np.linspace(1, 3, 20)
for b in B_vals: 
    J_numeric = np.asarray(J_at_fp.evalf(subs={A:a, B:b})).astype(np.float64)
    w, _ = eig(J_numeric)
    eigen_real.append(w[0].real)
    eigen_imag.append(abs(w[0].imag))
eigen_real = np.asanyarray(eigen_real)
eigen_imag = np.asarray(eigen_imag)
```


```python
fig, ax = plt.subplots(1, 2)
fig.set_size_inches(10,5)
fig.subplots_adjust(wspace=0.5)
ax[0].axhline(y=0, c="k", ls="dashed")
ax[0].plot(B_vals, eigen_real)
ax[0].set_ylabel(r"$\mathfrak{Re}(\lambda)$", fontsize = 20)
ax[0].set_xlabel(r"$B$", fontsize = 20)
ax[1].set_ylabel(r"$|\mathfrak{Im}(\lambda)|$", fontsize = 20)
ax[1].set_xlabel(r"$B$", fontsize = 20)
ax[1].plot(B_vals, eigen_imag)
```

Hopf bifurcation, is only one type of bifurcation, albeit a very important one for physical and biological systems. There are other types of bifurcation in which one can create or destroy fixed points or alter their properties in different ways than a Hopf bifurcations does. If you are curious, I suggest you to perform your own numerical experiements by playing with the values of $A$, $B$
or both.

### An Animation of the Hopf Bifurcation 


```python
from matplotlib import animation, rc
from IPython.display import HTML

# Brusselator System:
def dX_dt(A, B, X, t):
    x, y = X[0], X[1]
    return np.array([A + x**2 * y - B * x -x, B * x - x**2 * y])

T = 50 * np.pi                      # simulation time
dt = 0.01                           # integration time step

# time steps to be used in integration of the Brusselator system
t=np.arange(0, T, dt)    


num_traj = 5                          # number of trajectories
xmin, xmax  = 0,  5                   # min and max values of x axis in the plot
ymin, ymax  = 0,  5                   # min and max values of y axis in the plot
A = 1.                                # we will keep A parameter constant                    

# vary B parameter 
Bmin, Bmax, numB = 1., 3., 100        # min, max, number of steps for varying B
Bvals = np.linspace(Bmin, Bmax, numB)

# set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
fig.set_size_inches(8,8)
ax = plt.axes(xlim=(xmin, xmax), ylim=(ymin, ymax))
ax.set_ylabel("Y", fontsize = 20)
ax.set_xlabel("X", fontsize = 20)

    
# choose a set of initial points for our trajectories (in each frame we will use the same set)
X0  = list(zip(np.linspace(xmin, xmax, num_traj), np.linspace(ymin, ymax, num_traj)))     

# choose a color set for our trajectories
vcolors = plt.cm.jet_r(np.linspace(0., 1., num_traj))  

# prepare the mesh grid    
x = np.linspace(xmin, xmax, 15)             # divide x axis to intervals
y = np.linspace(ymin, ymax, 15)             # divide y axis to intervals
X1 , Y1  = np.meshgrid(x, y)                # from these intervals create a grid
        
# set up the lines, the quiver and the text object
lines = [ax.plot([], [], [], '-', c=c, lw=2)[0] for c in vcolors]
Q = ax.quiver(X1, Y1, [], [], [], pivot='mid', cmap=plt.cm.jet)
text = ax.text(0.02, 0.95, '', fontsize=20, transform=ax.transAxes)

# initialization function: plot the background of each frame. Needs to return each object to be updated
def init():
    for line in lines:
        line.set_data([], [])
    Q.set_UVC([], [], [])
    text.set_text("")
    return Q, lines, text

# animation function.  This is called sequentially
def animate(i):
    B = Bvals[i]
    DX1, DY1 = dX_dt(A, B, [X1, Y1], t)     # compute rate of change of the concentrations on grid points
    M = (np.hypot(DX1, DY1))                # norm of the rate of changes 
    M[ M == 0] = 1.                         # prevention against divisions by zero
    DX1 /= M                                # we normalize the direction field vector (each has unit length now)
    DY1 /= M                                # we normalize the direction field vector (each has unit length now)
    Q.set_UVC(DX1, DY1, M)
    # integrate the Brusellator ODE's for the set of trajectories, store them in X
    for line, X0i in zip(lines, X0):
        X = integrate.odeint(lambda x, t: dX_dt(A, B, x, t), X0i,t)    
        x, y = X.T        # get x and y for current trajectory
        line.set_data(x, y) 
    text.set_text("A={:.2f}, B={:.2f}".format(A, B))
    return Q, lines, text 
              
# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=100, interval=30, blit=False)
# instantiate the animator.
#anim = animation.FuncAnimation(fig, animate, init_func=init, frames=1000, interval=200, blit=True)
#HTML(anim.to_html5_video())
rc('animation', html='html5')
plt.close()
anim
```








In the animation above, we see how the direction field gets modified as $B$ is varied. Also shown several trajectories that are initalized at various points (I have chosen them on the $Y=X$ line here).

## Notes:

Should you encounter difficulty in running the embedded animation, try to launch Jupter Notebook using the command:<br>
`jupyter notebook --NotebookApp.iopub_data_rate_limit=10000000000`

## References:

1) Strogatz, S.H (2015). *Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Second Edition*, Boulder, USA: Westview Press.  <br>  
2) https://en.wikipedia.org/wiki/Brusselator

ea2d14f88f0858fe87c6488a202633885c6ef02c
ipynb
Jupyter Notebook
.ipynb_checkpoints/brusselator_hopf-checkpoint.ipynb
cemozen/pattern_formation_in_reaction-diffusion_systems
7788c2dec71bcbe47758cabdc9816d99f88df589
MIT
2021-04-04T02:01:50.000Z
2021-04-04T02:01:50.000Z
.ipynb_checkpoints/brusselator_hopf-checkpoint.ipynb
cemozen/pattern_formation_in_reaction-diffusion_systems
7788c2dec71bcbe47758cabdc9816d99f88df589
MIT
.ipynb_checkpoints/brusselator_hopf-checkpoint.ipynb
cemozen/pattern_formation_in_reaction-diffusion_systems
7788c2dec71bcbe47758cabdc9816d99f88df589
MIT
Qwen/Qwen-72B
1. YES
2. YES