Lecture slides for UCLA LS 30B, Spring 2020

Public worksheets/Lecture slides - Spring 2020 / Mini-Lecture 7-10 - Hopf bifurcation in the HPG model.ipynb

¹⁷⁹⁴⁴ views
License: GPL3
ubuntu2004

Kernel: SageMath 9.2

In [1]:

import numpy as np
import plotly.graph_objects
import plotly.subplots

In [2]:

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"

In [3]:

def plotly_text(text, location, **options):
    useroptions = options
    options = dict(showarrow=False)
    if useroptions.get("paper", False):
        options.update(xanchor="left", yanchor="bottom")
    options.update(useroptions)
    options.setdefault("font_color", options.pop("color", "black"))
    if options.pop("paper", False):
        options.update(xref="paper", yref="paper")
    x, y = np.array(location, dtype=float)
    return plotly.graph_objects.layout.Annotation(text=text, x=x, y=y, **options)

In [4]:

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)

def latex_eigenvalues(M, digits=4, show_abs=False):
    eigenvalues = []
    for l in sorted(M.eigenvalues(), key=lambda l: -abs(l)):
        if l.imag_part() < 0:
            continue
        if l.imag_part():
            text = fr"{float(l.real()):.{digits}f} \pm {float(l.imag()):.{digits}f} i"
            if show_abs:
                text += fr" \text{{ (abs. value {round(abs(l), digits)}) }}"
            eigenvalues.append(text)
        else:
            eigenvalues.append(f"{float(l):.{digits}f}")
    return r", \ ".join(eigenvalues)

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 terms of eigenvalues.

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

In [5]:

def hpg_hopf():
    fig = plotly.graph_objects.FigureWidget()
    fig.layout.margin.update(t=20, b=0, l=0, r=0)
    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()

Out[5]:

Interactive function <function hpg_hopf.<locals>.update at 0x7f61611ee040> with 1 widget
  n: TransformIntSlid…

FigureWidget({
    'data': [{'marker': {'color': 'purple', 'size': 3},
              'mode': 'markers',
      …

The Jacobian of the HPG model:

$\qquad \begin{cases} H' = \frac{1}{1 + G^n} - 0.2H \\ P' = H - 0.2P \\ G' = P - 0.2G \end{cases}$

In [6]:

state_vars = list(var("H, P, G"))
n = var("n")
system = (
    1/(1 + G^n) - 1/5*H,
    H - 1/5*P,
    P - 1/5*G,
)
field(H, P, G) = system
J = jacobian(field, (H, P, G))
show(r"J = " + latex(J(H, P, G)))

Out[6]:

In [7]:

def hpg_hopf_eigendiagram():
    state_vars = list(var("H, P, G"))
    n = var("n")
    vectorfield(H, P, G) = (1/(1 + G^n) - 0.2*H, H - 0.2*P, P - 0.2*G)

    fig = plotly.subplots.make_subplots(rows=int(1), cols=int(2), specs=[[{"type": "scene"}, {}]])
    fig = plotly.graph_objects.FigureWidget(fig)
    fig.layout.margin.update(t=40, b=0, l=0, r=0)
    fig.layout.scene.domain.x = [0.00, 0.70]
    fig.layout.scene.domain.y = [0.00, 0.95]
    fig.layout.xaxis.domain   = [0.75, 1.00]
    fig.layout.yaxis.domain   = [0.18, 0.51]
    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.layout.xaxis.title.text = "Real"
    fig.layout.yaxis.title.text = "Imaginary"
    fig.layout.xaxis.title.standoff = 5
    fig.layout.yaxis.title.standoff = 0
    fig.layout.xaxis.range = [-0.9, 0.9]
    fig.layout.yaxis.range = [-0.5, 0.5]

    fig.add_annotation(plotly_text("Trajectory", (0.35, 0.97), paper=True, font_size=18, xanchor="center"))
    fig.add_annotation(plotly_text("Equilibrium point:", (0.72, 0.97), paper=True, font_size=18))
    fig.add_annotation(plotly_text("", (0.76, 0.92), paper=True, font_size=14))
    fig.add_annotation(plotly_text("Jacobian:", (0.72, 0.82), paper=True, font_size=18))
    fig.add_annotation(plotly_text("", (0.76, 0.65), paper=True, font_size=14))
    fig.add_annotation(plotly_text("Eigenvalues:", (0.72, 0.53), paper=True, font_size=18))
    fig.add_annotation(plotly_text("", (0.74, 0.00), paper=True, font_size=14))
    l1, l2, eqpt_label, l3, jacobian_label, l4, eigenvalues_label = fig.layout.annotations
    fig.add_scatter3d(row=int(1), col=int(1), mode="markers", marker_color="purple", marker_size=3)
    fig.add_scatter3d(row=int(1), col=int(1), mode="lines", line_color="red")
    fig.add_scatter3d(row=int(1), col=int(1), mode="lines", line_color="fuchsia")
    fig.add_scatter3d(row=int(1), col=int(1), mode="lines", line_color="purple", line_width=10)
    fig.add_scatter(row=int(1), col=int(2), mode="markers", marker_color="red")
    eqpt, trajectory1, trajectory2, attractor, eigenvalues = fig.data

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

        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_label.text = fr"$({float(H0):.04f}, {float(P0):.04f}, {float(G0):.04f})$"
            jacobian_label.text = fr"${latex(J.round(3))}$"
            eigenvalues_label.text = fr"${latex_eigenvalues(J)}$"
            eigenvalues.x = [a.real() for a in (e1, e2, e3)]
            eigenvalues.y = [a.imag() for a in (e1, e2, e3)]
            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_eigendiagram()

Out[7]:

Interactive function <function hpg_hopf_eigendiagram.<locals>.update at 0x7fa2245096a8> with 1 widget
  n: Tra…

FigureWidget({
    'data': [{'marker': {'color': 'purple', 'size': 3},
              'mode': 'markers',
      …

A bifurcation diagram for the HPG model:

In [8]:

def hpg_bifurcation_diagram():
    n = var("n")
    vectorfield(H, P, G) = (1/(1 + G^n) - 0.2*H, H - 0.2*P, P - 0.2*G)
    J = jacobian(vectorfield, (H, P, G))

    def stability(n):
        G0 = find_root(1/(1 + G^n) - 0.008*G, 0, 20)
        J0 = matrix(RDF, J.subs(n=n)(0.04*G0, 0.2*G0, G0))
        e1, e2, e3 = J0.eigenvalues()
        return e1.real() if abs(e1.imag()) > 1e-10 else e2.real()

    bif_point = find_root(stability, 1, 12)
    base = plot(stability, (n, 1, 14), axes_labels=("$n$", "Real part"))
    @interact(b=("Show bifurcation point", False), l=("Show labels", False))
    def update(b, l):
        p = base
        if b:
            p += line(((bif_point, p.ymin()), (bif_point, p.ymax())),
                color="red", linestyle="dashed")
            p += text("n = {}".format(bif_point), (bif_point, p.ymax()),
                color="black", vertical_alignment="bottom")
        if l:
            p += text("Stable spiral\n(oscillations don't last)", (bif_point/2, p.ymin()/2),
                color="red", fontsize=12, horizontal_alignment="left")
            p += text("Limit cycle\nattractor", ((bif_point + p.xmax())/2, p.ymax()/2),
                color="red", fontsize=12, horizontal_alignment="right", vertical_alignment="bottom")
        show(p, xmin=0, figsize=5)

hpg_bifurcation_diagram()

Out[8]:

Interactive function <function hpg_bifurcation_diagram.<locals>.update at 0x7fa224509d08> with 2 widgets
  b: …

Conclusion:

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

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

Learning goals:

Product

Resources

Company