Mini-Lecture 8 - How a Hopf bifurcation occurs.ipynb

Kernel: SageMath 9.3

import numpy as np
import plotly.graph_objects

mydefault = plotly.graph_objects.layout.Template()
mydefault.layout.hovermode = False
mydefault.layout.scene.hovermode = False
mydefault.layout.xaxis.showspikes = False
mydefault.layout.yaxis.showspikes = False
mydefault.layout.scene.xaxis.showspikes = False
mydefault.layout.scene.yaxis.showspikes = False
mydefault.layout.scene.zaxis.showspikes = False
plotly.io.templates["mydefault"] = mydefault
plotly.io.templates.default = "mydefault"

Learning goal:

Be able to explain the two conditions needed for a Hopf bifurcation to occur.

Recall the Holling–Tanner model: $\qquad \begin{cases} N' = r_1 N \cdot \left( 1 - \frac{N}{k} \right) - w \cdot \frac{N}{d + N} \cdot P \\ P' = r_2 P \cdot \left( 1 - \frac{jP}{N} \right) \end{cases}$

Parameters: $$\begin{aligned} r_1 &= \text{natural per-capita growth rate of the prey population} \\ r_2 &= \text{natural per-capita growth rate of the predator population} \\ j &= \text{size of prey population needed to support each predator} \\ w &= \text{maximum rate at which one predator consumes prey (per year)} \\ d &= \text{prey population at which the predation rate will be half of that maximum} \\ k &= \text{maximum prey population supported by the environment (carrying capacity)} \end{aligned}$$

Here I will be using the following values of these parameters:

$$r_1 = 0.4, \qquad r_2 = 0.03, \qquad j = 150, \qquad w = 300, \qquad d = 1000, \qquad k = 3000$$

var("r_1, r_2, k, j, w, d")
def holling_tanner_hopf(vary):
    field(N, P) = ( r_1*N*(1 - N/k) - w*P*N/(d + N), r_2*P*(1 - j*P/N) )
    J = jacobian(field, (N, P))
    b = (d - k + k*w/j/r_1)/2
    N_eq = sqrt(b^2 + d*k) - b
    P_eq = N_eq/j
    J_eq = J(N_eq, P_eq)

    params = {
        r_1: 0.4,      # Natural per-capita growth rate of the prey population
        r_2: 0.03,     # Natural per-capita growth rate of the predator population
        j:  150,       # Size of prey population needed to support each predator
        w:  300,       # Maximum rate at which one predator consumes prey (per year)
        d:  1000,      # Prey population at which the predation rate will be half of that maximum
        k:  3000,      # Maximum prey population supported by the environment (carrying capacity)
    }
    param_ranges = {
        r_1: (0.2, 0.8, 0.01),
        r_2: (0.01, 0.10, 0.001),
        j:   (50, 300, 10),
        w:   (200, 800, 10),
        d:   (300, 2000, 10),
        k:   (1800, 5000, 10),
    }
    t_range = srange(0, 8000, 0.1)
    initial_state = (220, 0.8)

    @interact(value=slider(*param_ranges[vary], default=params[vary], label="${}$".format(vary)))
    def update(value):
        params[vary] = value
        myfield = field.subs(params)
        N0, P0 = (N_eq.subs(params), P_eq.subs(params))
        realpart = J_eq.subs(params).trace()/2
        p = plot_vector_field(myfield, (N, 0, 2000), (P, 0, 4), color="limegreen", frame=False)
        p += point([(N0, P0)], size=40, color="purple")
        spiralin = desolve_odeint(myfield, initial_state, t_range, [N, P])
        if realpart > 0:
            p += list_plot(spiralin[:-4000], plotjoined=True, color="red")
            near_eqpt = (N0, P0 + 0.1)
            spiralout = desolve_odeint(myfield, near_eqpt, t_range, [N, P])
            p += list_plot(spiralout, plotjoined=True, color="fuchsia")
            p += text("Unstable spiral eq pt", (1250, 3.7), color="black", fontsize=18)
            p += text("(and LCA)", (1250, 3.2), color="purple", fontsize=18)
            p += list_plot(spiralin[-4000:], plotjoined=True, color="purple", thickness=2)
        else:
            p += list_plot(spiralin, plotjoined=True, color="red")
            p += text("Stable spiral eq pt", (1250, 3.7), color="black", fontsize=18)
        p.show(xmax=2000, ymax=4, axes_labels=("$N$ (Prey)", "$P$ (Predators)"), figsize=6)

Hopf bifurcation in the Holling–Tanner model, as we vary $r_1$:

holling_tanner_hopf(r_1)

Hopf bifurcation in the Holling–Tanner model, as we vary $d$:

holling_tanner_hopf(d)

Hopf bifurcation in the Holling–Tanner model, as we vary $j$:

holling_tanner_hopf(j)

Hopf bifurcation in the HPG model: $\qquad \begin{cases} H' = \frac{1}{1 + G^n} - 0.2H \\ P' = H - 0.2P \\ G' = P - 0.2G \end{cases}$

def hpg_hopf():
    fig = plotly.graph_objects.FigureWidget()
    fig.layout.showlegend = False
    fig.layout.scene.xaxis.title.text = "H (GnRH)"
    fig.layout.scene.yaxis.title.text = "P (FSH and LH)"
    fig.layout.scene.zaxis.title.text = "G (estrogen)"
    fig.add_scatter3d(mode="markers", marker_color="purple", marker_size=3)
    eqpt = fig.data[-1]
    fig.add_scatter3d(mode="lines", line_color="red")
    trajectory1 = fig.data[-1]
    fig.add_scatter3d(mode="lines", line_color="fuchsia")
    trajectory2 = fig.data[-1]
    fig.add_scatter3d(mode="lines", line_color="purple", line_width=10)
    attractor = fig.data[-1]

    @interact(n=slider(1, 20, 1, default=4))
    def update(n):
        state_vars = list(var("H, P, G"))
        system = (
            1/(1 + G^n) - 0.2*H,
            H - 0.2*P,
            P - 0.2*G,
        )
        field(H, P, G) = system
        G0 = find_root(1/(1 + G^n) - 0.008*G, 0, 100)
        P0 = 0.2*G0
        H0 = 0.2*P0
        J = jacobian(field, state_vars)(H0, P0, G0)
        e1, e2, e3 = J.eigenvalues()
        realpart = e1.real().n() if abs(e1.imag()) > 1e-10 else e2.real().n()

        initial_state = (0.5, 1, 1)
        t_range = srange(0, 400, 0.05)
        spiralin = desolve_odeint(field, initial_state, t_range, state_vars)
        if realpart > 0:
            spiralout = desolve_odeint(field, (H0 + 0.005, P0 + 0.02, G0 + 0.01), t_range, state_vars)

        with fig.batch_update():
            eqpt.x, eqpt.y, eqpt.z = [H0], [P0], [G0]
            end = len(spiralin) - (1000 if realpart > 0 else 0)
            trajectory1.x = spiralin[:end,0]
            trajectory1.y = spiralin[:end,1]
            trajectory1.z = spiralin[:end,2]
            if realpart > 0:
                attractor.x = spiralin[end:,0]
                attractor.y = spiralin[end:,1]
                attractor.z = spiralin[end:,2]
                trajectory2.x = spiralout[:,0]
                trajectory2.y = spiralout[:,1]
                trajectory2.z = spiralout[:,2]
            trajectory2.visible = (realpart > 0)
            attractor.visible = (realpart > 0)
    return fig
hpg_hopf()

Hopf bifurcation in the Rayleigh clarinet model: $\qquad \begin{cases} X' = V \\ V' = -X - (V^3 - aV) \end{cases}$

def rayleigh_hopf():
    t_range = srange(0, 100, 0.1)
    initial_state = (1.6, 1.6)
    X0, V0 = (0, 0)

    @interact(a=slider(-1, 1, 0.1, default=-1))
    def update(a):
        vectorfield(X, V) = ( V, -X - (V^3 - a*V) )
        realpart = jacobian(vectorfield, (X, V))(X0, V0).trace().n() / 2
        p = plot_vector_field(vectorfield, (X, -2.2, 2.2), (V, -2.2, 2.2), color="limegreen", frame=False)
        p += point([(X0, V0)], size=40, color="purple")
        spiralin = desolve_odeint(vectorfield, initial_state, t_range, [X, V])
        if realpart > 0:
            spiralout = desolve_odeint(vectorfield, (X0 + 0.01, V0), t_range, [X, V])
            p += list_plot(spiralout, plotjoined=True, color="fuchsia")
            p += list_plot(spiralin[:-200], plotjoined=True, color="red")
            p += list_plot(spiralin[-200:], plotjoined=True, color="purple", thickness=2)
        else:
            p += list_plot(spiralin, plotjoined=True, color="red")
        p.show(axes_labels=("$X$ (position)", "$V$ (velocity)"), figsize=6)
rayleigh_hopf()

Conclusion:

A Hopf bifurcation occurs when, as we change a parameter, the following happens:

1. An equilibrium point changes from a stable spiral to an unstable spiral (or vice-versa).

2. Meanwhile, in some part of the state space around that equilibrium point, other trajectories continue to spiral inward.

Notes:

• Just as before, condition (2) is hard to check. We'll just rely on simulation for that.

• By the end of this quarter, you will learn how to actually detect condition (1) mathematically. As I have said many times now, this is a major unsolved problem at this point in LS 30, but one that we will solve eventually. This is A Big Deal.