Mini-Lecture 7-09 - How a Hopf bifurcation occurs, part 2.ipynb

Kernel: SageMath 9.3

import numpy as np

def round_complex(z, digits):
    if z.imag_part():
        return round(z.real_part(), digits) + round(z.imag_part(), digits) * I
    return round(z, digits)

Learning goals:

• Be able to explain the two conditions needed for a Hopf bifurcation to occur. (review)

• Be able to explain one of those two conditions in a mathematical (numerical) way. Specifically, in terms of eigenvalues.

• Be able to use this to find the exact parameter value at which a Hopf bifurcation occurs.

Recall the definition of a Hopf bifurcation:

Definition: When changing a parameter of the model causes the long-term (steady-state) behavior to change from equilibrium behavior to oscillating, we say that a Hopf bifurcation has occurred.

Or, in other words, when changing a parameter causes a limit cycle attractor to appear where there previously was not one, that's a Hopf bifurcation.

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, 2000, 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)

Review of our conclusion from way back then:

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:

• Condition (2) is hard to check. We'll just rely on simulation for that.

• We now finally know how to actually detect condition (1) mathematically!

What we now know:

• An equilibrium point is a stable spiral if the eigenvalues of the Jacobian are complex, with a negative real part.

• An equilibrium point is an unstable spiral if the eigenvalues of the Jacobian are complex, with a positive real part.

The Holling–Tanner equations, again:

$$\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}$$

Its Jacobian:

$$\left( \begin{array}{cc} -r_{1} {\left(\frac{N}{k} - 1\right)} - \frac{N r_{1}}{k} - \frac{P w}{N + d} + \frac{N P w}{{\left(N + d\right)}^{2}} & -\frac{N w}{N + d} \\ \frac{P^{2} j r_{2}}{N^{2}} & -{\left(\frac{P j}{N} - 1\right)} r_{2} - \frac{P j r_{2}}{N} \end{array} \right)$$

var("r_1, r_2, k, j, w, d")
def holling_tanner_hopf_eigenvalues(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, 2000, 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))
        J0 = matrix(RDF, J_eq.subs(params))
        realpart = J0.trace() / 2
        imagpart = abs(J0.eigenvalues()[0].imag())
        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)
        print(f"Eigenvalues of the Jacobian:  {float(realpart):.5f}  ±  {float(imagpart):.5f}i")

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

holling_tanner_hopf_eigenvalues(r_1)

var("r_1, r_2, k, j, w, d")
def holling_tanner_hopf_eigenvalue_diagram(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, 2000, 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))
        J0 = matrix(RDF, J_eq.subs(params))
        realpart = J0.trace() / 2
        imagpart = abs(J0.eigenvalues()[0].imag())
        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.axes_labels(("$N$ (Prey)", "$P$ (Predators)"))
        p.set_axes_range(xmin=0, xmax=2000, ymin=0, ymax=4)
        p.set_aspect_ratio(0.75*2000/4)
        eigendiagram = point(J0.eigenvalues(), size=40, color="red")
        eigendiagram += text("Eigenvalues:", (0,0.40), color="black", fontsize=16)
        eigendiagram += text(fr"${float(realpart):.5f} \pm {float(imagpart):.5f} i$", (0,0.35), 
                             color="black", fontsize=16)
        eigendiagram.axes_labels(("Real", "Imaginary"))
        eigendiagram.set_axes_range(xmin=-0.1, ymin=-0.2, xmax=0.1, ymax=0.2)
        both = multi_graphics([
            (p,            (0   , 0   , 0.65, 1.00)), 
            (eigendiagram, (0.70, 0.20, 0.30, 0.40)), 
        ])
        both.show(figsize=10)

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

holling_tanner_hopf_eigenvalue_diagram(r_1)

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

holling_tanner_hopf_eigenvalue_diagram(d)

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

holling_tanner_hopf_eigenvalue_diagram(j)

var("r_1, r_2, k, j, w, d")
def holling_tanner_bifurcation_diagram_1D():
    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),
    }

    @interact(vary=("Parameter to vary", params.keys()), 
              show_bifpt=("Show bifurcation point", False),
              show_labels=("Show labels", False))
    def update(vary, show_bifpt, show_labels):
        myparams = params.copy()
        del myparams[vary]
        realpart = J_eq.subs(myparams).trace() / 2
        vmin, vmax, vdef = param_ranges[vary]
        bif_point = find_root(realpart, vmin, vmax)
        p = plot(realpart, (vary, vmin, vmax), axes_labels=(f"${vary}$", "Real part"))

        if show_bifpt:
            p += line(((bif_point, p.ymin()), (bif_point, p.ymax())),
                color="red", linestyle="dashed")
            p += text(f"${vary} = {float(bif_point):.5f}$", (bif_point, p.ymax()),
                color="black", fontsize=16, vertical_alignment="bottom")
        if show_labels:
            for xval in (vmin + bif_point)/2, (bif_point + vmax)/2:
                yval = realpart.subs({vary: xval})
                if yval < 0:
                    p += text("Stable spiral\n(oscillations don't last)",
                        (xval, yval), color="red", fontsize=12)
                else:
                    p += text("Limit cycle\nattractor",
                        (xval, yval), color="red", fontsize=12)
        show(p)

holling_tanner_bifurcation_diagram_1D()

var("r_1, r_2, k, j, w, d")
def holling_tanner_bifurcation_diagram_2D():
    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),
    }

    @interact(vary1=("First parameter to vary (horiz. axis)", params.keys()),
              vary2=("Second parameter to vary (vert. axis)", params.keys()))
    def update(vary1, vary2):
        if vary1 is vary2:
            print("Please select two *different* parameters.")
            return
        myparams = params.copy()
        del myparams[vary1]
        del myparams[vary2]
        realpart = J_eq.subs(myparams).trace()/2
        v1min, v1max, v1def = param_ranges[vary1]
        v2min, v2max, v2def = param_ranges[vary2]

        p = region_plot(realpart > 0, (vary1, v1min, v1max), (vary2, v2min, v2max), 
            incol="lightpink", outcol="lightgreen", bordercol="blue")
        options = {"fontsize": 14, "horizontal_alignment": "left", "axis_coords": True}
        p += text("Legend:",       (1.02,0.90), color="black", **options)
        p += text("Limit cycle",   (1.02,0.84), color="red",   **options)
        p += text("Stable spiral", (1.02,0.78), color="green", **options)
        p.set_axes_range(xmin=v1min, xmax=v1max, ymin=v2min, ymax=v2max)
        p.show(axes_labels=(f"${vary1}$", f"${vary2}$"), aspect_ratio="auto", figsize=9)

holling_tanner_bifurcation_diagram_2D()

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).

   • This means the eigenvalues of the Jacobian at that equilibrium point are complex, and the real part changes sign (from negative to positive or vice-versa).

   • In particular, the Hopf bifurcation point is the parameter value at which the real part of the eigenvalues is exactly $0$.

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

Note: Condition (2) is still hard to check. We'll just rely on simulation for that.