Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Lecture slides for UCLA LS 30B, Spring 2020

Views: 14462
License: GPL3
Image: ubuntu2004
Kernel: SageMath 9.2
import numpy as np
import plotly.graph_objects
import plotly.subplots


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"


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)


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: {H′=11+Gn−0.2HP′=H−0.2PG′=P−0.2G\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.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()

The Jacobian of the HPG model:

{H′=11+Gn−0.2HP′=H−0.2PG′=P−0.2G\qquad \begin{cases} H' = \frac{1}{1 + G^n} - 0.2H \\ P' = H - 0.2P \\ G' = P - 0.2G \end{cases}

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

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

A bifurcation diagram for the HPG model: 

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

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

  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.