Mini-Lecture 6 - Behavior of the HPG model.ipynb

Kernel: SageMath 9.4

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"

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

@interact(n=slider(1, 20, 1, default=12))
def steepness(n):
    sigmoid(G) = 1/(1 + G^n)
    p = plot(sigmoid(G), (G, 0, 3), thickness=2)
    p += text(r"$f(G) = \dfrac{1}{1 + G^{" + str(n) + "}}$", (2, 1), color="black", fontsize=20)
    p.show(ymin=0, ymax=1.3, axes_labels=("$G$", "Production rate of $H$"), figsize=5)

Time series of the HPG model:

@interact(n=slider(1, 20, 1, default=12, label="$n$"), 
          tmax=slider([200, 400, 800], default=200, label="$t_{max}$"))
def hpg_timeseries(n, tmax):
    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

    initial_state = (0.5, 1, 1)
    t_range = srange(0, tmax, 0.05)
    solution = desolve_odeint(field, initial_state, t_range, state_vars)
    solution = np.insert(solution, 0, t_range, axis=1)
    p  = list_plot(solution[:,(0,1)], plotjoined=True, color="blue", legend_label="$H$ (GnRH)")
    p += list_plot(solution[:,(0,2)], plotjoined=True, color="red",  legend_label="$P$ (FSH and LH)")
    p += list_plot(solution[:,(0,3)], plotjoined=True, color="gold", legend_label="$G$ (estrogen)")
    p += text(r"$n = {}$".format(n), (tmax/2,4), color="black", fontsize=16)
    p.show(ymin=0, ymax=5, axes_labels=("$t$", "Hormone levels"))

A new type of 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.

Of course, this can also happen in the other direction: changing a parameter can cause the system to change from oscillating behavior to equilibrium behavior. That's still called a Hopf bifurcation.

Searching for the Hopf bifurcation:

@interact(n=slider(7, 9, RDF(0.1), default=8, label="$n$"), 
          tmax=slider([200, 400, 800, 2000], default=200, label="$t_{max}$"))
def hpg_timeseries(n, tmax):
    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

    initial_state = (0.5, 1, 1)
    t_range = srange(0, tmax, 0.05)
    solution = desolve_odeint(field, initial_state, t_range, state_vars)
    solution = np.insert(solution, 0, t_range, axis=1)
    p  = list_plot(solution[:,(0,1)], plotjoined=True, color="blue", legend_label="$H$ (GnRH)")
    p += list_plot(solution[:,(0,2)], plotjoined=True, color="red",  legend_label="$P$ (FSH and LH)")
    p += list_plot(solution[:,(0,3)], plotjoined=True, color="gold", legend_label="$G$ (estrogen)")
    p += text(r"$n = {}$".format(n), (tmax/2,4), color="black", fontsize=16)
    p.show(ymin=0, ymax=5, axes_labels=("$t$", "Hormone levels"))

Biologically, what does the Hopf bifurcation represent in the HPG model?

That is, in the real biological HPG system, what does it mean for the system to change from not oscillating to oscillating?

Recall that we said that oscillations in this system give a very basic mathematical model of how the menstrual cycle occurs in females.

So the Hopf bifurcation corresponds to the onset of the menstrual cycle.

In short: Puberty is a Hopf bifurcation.