Mini-Lecture 7 - How a limit cycle attractor forms.ipynb

Kernel: SageMath 9.0

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 how a limit cycle attractor forms in the state space.

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

def holling_tanner_trajectories():
    var("r_1, r_2, k, j, w, d")
    field(N, P) = ( r_1*N*(1 - N/k) - w*P*N/(d + N), r_2*P*(1 - j*P/N) )
    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)
    }
    field = field.subs(params)
    b = (d - k + k*w/j/r_1)/2
    N_eq = sqrt(b^2 + d*k) - b
    P_eq = N_eq/j
    N0, P0 = (N_eq.subs(params), P_eq.subs(params))

    t_range = srange(0, 2000, 0.1)
    initial_states = {
        (220, 0.8): "red", 
        (1500, 3): "orange", 
        (N0 + 20, P0 + 0.03): "fuchsia", 
    }

    trajectories = []
    attractor = None
    for state, color in initial_states.items():
        solution = desolve_odeint(field, state, t_range, [N, P])
        if attractor is None:
            attractor = list_plot(solution[-2000:], plotjoined=True, color="purple", thickness=2)
        trajectories.append(list_plot(solution[:-2000], plotjoined=True, color=color))

    @interact(which=slider(srange(1, len(trajectories) + 1, 1) + ["All"], default=1), 
              show_eqpt=checkbox(False, label="Show equilibrium point"))
    def update(which, show_eqpt):
        p = plot_vector_field(field, (N, 0, 2000), (P, 0, 4), color="limegreen", frame=False)
        if which == "All":
            p += sum(trajectories)
        else:
            p += trajectories[which - 1]
        p += attractor
        if show_eqpt:
            p += point([N0, P0], size=40, color="purple")
        p.show(axes_labels=("$N$ (Prey)", "$P$ (Predators)"), figsize=6)
holling_tanner_trajectories()

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

def hpg_trajectories():
    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)"
    xmin, ymin, zmin = (0, 0, 0)
    xmax, ymax, zmax = (1, 2, 6)
    fig.add_scatter3d(x=(xmin, xmax), y=(ymin, ymax), z=(zmin, zmax), mode="markers", marker_size=0.1)

    n = 12
    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
    t_range = srange(0, 400, 0.05)
    initial_states = {
        (0.5, 1, 1): "red", 
        (0.05, 0.8, 4): "orange", 
        (H0 + 0.005, P0 + 0.02, G0 + 0.01): "fuchsia", 
    }

    trajectories = []
    attractor = None
    for state, color in initial_states.items():
        solution = desolve_odeint(field, state, t_range, state_vars)
        if attractor is None:
            solution, attractor = solution[:-2000], solution[-2000:]
            fig.add_scatter3d(x=attractor[:,0], y=attractor[:,1], z=attractor[:,2], 
                              mode="lines", line_color="purple", line_width=10)
            attractor = fig.data[-1]
        fig.add_scatter3d(x=solution[:,0], y=solution[:,1], z=solution[:,2], 
                          mode="lines", line_color=color)
        trajectories.append(fig.data[-1])
    fig.add_scatter3d(x=[H0], y=[P0], z=[G0], mode="markers", marker_color="purple", marker_size=3)
    eqpt = fig.data[-1]

    @interact(which=slider(srange(1, len(trajectories)+1, 1) + ["All"], default=1), 
              show_eqpt=checkbox(False, label="Show equilibrium point"))
    def update(which, show_eqpt):
        with fig.batch_update():
            for i, trajectory in enumerate(trajectories):
                trajectory.visible = (which == i + 1 or which == "All")
            eqpt.visible = show_eqpt
    return fig
hpg_trajectories()

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

(from the textbook, section 4.1)

def rayleigh_trajectories():
    vectorfield(X, V) = ( V, -X - (V^3 - V) )
    t_range = srange(0, 100, 0.1)
    initial_states = {
        (0, -2): "red", 
        (1.6, 1.6): "orange", 
        (0, 2): "blue", 
        (0.01, 0): "fuchsia", 
    }
    X0, V0 = (0, 0)

    trajectories = []
    attractor = None
    for state, color in initial_states.items():
        solution = desolve_odeint(vectorfield, state, t_range, [X, V])
        if attractor is None:
            solution, attractor = solution[:-200], solution[-200:]
            attractor = list_plot(attractor, plotjoined=True, color="purple", thickness=2)
        trajectories.append(list_plot(solution, plotjoined=True, color=color))

    @interact(which=slider(srange(1, len(trajectories) + 1, 1) + ["All"], default=1), 
              show_eqpt=checkbox(False, label="Show equilibrium point"))
    def update(which, show_eqpt):
        p = plot_vector_field(vectorfield, (X, -2.2, 2.2), (V, -2.2, 2.2), color="limegreen", frame=False)
        if which == "All":
            p += sum(trajectories)
        else:
            p += trajectories[which - 1]
        p += attractor
        if show_eqpt:
            p += point([(X0, V0)], size=40, color="purple")
        p.show(axes_labels=("$X$ (position)", "$V$ (velocity)"), figsize=6)
rayleigh_trajectories()

Conclusion:

A limit cycle attractor forms when the following two conditions are met:

1. There is an equilibrium point that is an unstable spiral.

2. In some part of the state space around that equilibrium point, other trajectories spiral inward.

Notes:

• We will learn (later this quarter) how to detect condition (1).

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

Pitfalls:

• The “other trajectories spiral inward” does not mean that they are “stable spirals”. That term only applies to an equilibrium point. In general, you can't say that a trajectory is stable. Just say that those trajectories spiral inward.