2
$\begingroup$

TL:DR

How can one implement a bifurcation diagram of a seasonally forced epidemiological model such as SEIR (susceptible, exposed, infected, recovered) in Python? I already know how to implement the model itself and display a sampled time series (see this stackoverflow question), but I am struggling with reproducing figure from a textbook.

Context and My Attempt

I am trying to reproduce figures from the book "Modeling Infectious Diseases in Humans and Animals" (Keeling 2007) to both validate my implementations of models and to learn/visualize how different model parameters affect the evolution of a dynamical system.

enter image description here

I have found implementations of bifurcation diagrams for examples using the logistic map (see this ipython cookbook this pythonalgos bifurcation, and this stackoverflow question).

For the logistic map ($x_{n+1}=rx_n(1-x_n)$), the bifurcation diagram is plotted by computing a matrix $M \in \mathcal{R}^{T \times P}$ for $T$ iterations and $P$ parameters such that $\{r_j\}_{j=1}^{P}$ and $M_{ij}=LogisticMap(x=M_{(i-1)j}, r=r_j)$. The bifurcation diagram plots $M_{ij}\ \text{vs}\ r_j$. Based on this, I calculate the matrix $M$ for the SEIR model, which is a system of differential equations, by numerically integrating over $T$ timesteps for $P$ parameters such that $\{\beta_{1(j)}\}_{j=1}^{P}$ and the resultant matrix $M \in \mathcal{R}^{P \times T}$. I simply transpose this matrix to $M \in \mathcal{R}^{T \times P}$ and then plot this as $M_{ij}\ \text{vs}\ \beta_{1(j)}$, but with strange results as shown below.

enter image description here

enter image description here

Questions

Why do my panel outputs not match those in the figure? I assume I can't try to reproduce the bifurcation diagram from the textbook without my panels matching the output in the textbook.

I have tried to implement a function to produce a bifurcation diagram, but the output looks really strange. Am I misunderstanding something about the bifurcation diagram?

NOTE: I receive no warnings/errors during code execution.

Code to Reproduce my Figures

from typing import Callable, Dict, List, Optional, Any
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def seasonal_seir(y: List, t: List, params: Dict[str, Any]):
    """Seasonally forced SEIR model.
    
    Function parameters much match with those required
    by `scipy.integrate.odeint`

    Args:
        y: Initial conditions.
        t: Timesteps over which numerical solution will be computed.
        params: Dict with the following key-value pairs:
            beta_zero -- Average transmission rate.
            beta_one  -- Amplitude of seasonal forcing.
            omega     -- Period of forcing.
            mu        -- Natural mortality rate.
            sigma     -- Latent period for infection.
            gamma     -- Recovery from infection term.

    Returns:
        Tuple whose components are the derivatives of the
        susceptible, exposed, and infected state variables
        w.r.t to time.

    References:
        [SEIR Python Program from Textbook](http://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/Program_2.6/Program_2_6.py)
        [Seasonally Forced SIR Program from Textbook](http://homepages.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter5/Program_5.1/Program_5_1.py)
    """
    beta_zero = params['beta_zero']
    beta_one = params['beta_one']
    omega = params['omega']
    mu = params['mu']
    sigma = params['sigma']
    gamma = params['gamma']

    s, e, i = y 
    beta = beta_zero*(1 + beta_one*np.cos(omega*t))
    sdot = mu - (beta * i + mu)*s
    edot = beta*s*i - (mu + sigma)*e
    idot = sigma*e - (mu + gamma)*i
    return sdot, edot, idot

def plot_panels(
    model: Callable,
    model_params: Dict, 
    panel_param_space: List, 
    panel_param_name: str,
    initial_conditions: List,
    timesteps: List,
    odeint_kwargs: Optional[Dict] = dict(),
    x_ticks: Optional[List] = None,
    time_slice: Optional[slice] = None,
    state_var_ix: Optional[int] = None,
    log_scale: bool = False):
    """Plot panels that are samples of the parameter space for bifurcation.

    Args:
        model: Function that models dynamical system. Returns dydt.
        model_params: Dict whose key-value pairs are the names
            of parameters in a given model and the values of those parameters.
        bifurcation_parameter_space: List of varied bifurcation parameters.
        bifuraction_parameter_name: The name o the bifurcation parameter.
        initial_conditions: Initial conditions for numerical integration.
        timesteps: Timesteps for numerical integration.
        odeint_kwargs: Key word args for numerical integration.
        state_var_ix: State variable in solutions to use for plot.
        time_slice: Restrict the bifurcation plot to a subset
            of the all solutions for numerical integration timestep space.
    
    Returns:
        Figure and axes tuple.
    """

    # Set default ticks
    if x_ticks is None:
        x_ticks = timesteps

    # Create figure
    fig, axs = plt.subplots(ncols=len(panel_param_space))

    # For each parameter that is varied for a given panel
    # compute numerical solutions and plot
    for ix, panel_param in enumerate(panel_param_space):

        # update model parameters with the varied parameter
        model_params[panel_param_name] = panel_param

        # Compute solutions
        solutions = odeint(
            model,
            initial_conditions,
            timesteps,
            args=(model_params,),
            **odeint_kwargs)

        # If there is a particular solution of interst, index it
        # otherwise squeeze last dimension so that [T, 1] --> [T]
        # where T is the max number of timesteps
        if state_var_ix is not None:
            solutions = solutions[:, state_var_ix]

        elif state_var_ix is None and solutions.shape[-1] == 1:
            solutions = np.squeeze(solutions)
        
        else:
            raise ValueError(
                f'solutions to model are rank-2 tensor of shape {solutions.shape}'
                ' with the second dimension greater than 1. You must pass'
                ' a value to :param state_var_ix:')

        # Slice the solutions based on the desired time range
        if time_slice is not None:
            solutions = solutions[time_slice]

        # Natural log scale the results
        if log_scale:
            solutions = np.log(solutions)

        # Plot the results
        axs[ix].plot(x_ticks, solutions)

    return fig, axs

def plot_bifurcation(
    model: Callable,
    model_params: Dict,
    bifurcation_parameter_space: List,
    bifurcation_param_name: str,
    initial_conditions: List,
    timesteps: List,
    odeint_kwargs: Optional[Dict] = dict(),
    state_var_ix: Optional[int] = None,
    time_slice: Optional[slice] = None,
    log_scale: bool = False):
    """Plot a bifurcation diagram of state variable from dynamical system.

    Args:
        model: Function that models system. Returns dydt.
        model_params: Dict whose key-value pairs are the names
            of parameters in a given model and the values of those parameters.
        bifurcation_parameter_space: List of varied bifurcation parameters.
        bifuraction_parameter_name: The name o the bifurcation parameter.
        initial_conditions: Initial conditions for numerical integration.
        timesteps: Timesteps for numerical integration.
        odeint_kwargs: Key word args for numerical integration.
        state_var_ix: State variable in solutions to use for plot.
        time_slice: Restrict the bifurcation plot to a subset
            of the all solutions for numerical integration timestep space.
        log_scale: Flag to natural log scale solutions.

    Returns:
        Figure and axes tuple.
    """
    
    # Track the solutions for each parameter
    parameter_x_time_matrix = []

    # Iterate through parameters
    for param in bifurcation_parameter_space:

        # Update the parameter dictionary for the model
        model_params[bifurcation_param_name] = param

        # Compute the solutions to the model using 
        # dictionary of parameters (including the bifurcation parameter)
        solutions = odeint(
            model,
            initial_conditions,
            timesteps,
            args=(model_params, ),
            **odeint_kwargs)

        # If there is a particular solution of interst, index it
        # otherwise squeeze last dimension so that [T, 1] --> [T]
        # where T is the max number of timesteps
        if state_var_ix is not None:
            solutions = solutions[:, state_var_ix]

        elif state_var_ix is None and solutions.shape[-1] == 1:
            solutions = np.squeeze(solutions)
        
        else:
            raise ValueError(
                f'solutions to model are rank-2 tensor of shape {solutions.shape}'
                ' with the second dimension greater than 1. You must pass'
                ' a value to :param state_var_ix:')

        # Update the parent list of solutions for this particular 
        # bifurcation parameter
        parameter_x_time_matrix.append(solutions)

    # Cast to numpy array
    parameter_x_time_matrix  = np.array(parameter_x_time_matrix)

    # Transpose: Bifurcation plots Function Output vs. Parameter
    # This line ensures that each row in the matrix is the solution
    # to a particular state variable in the system of ODEs
    # a timestep t
    # and each column is that solution for a particular value of 
    # the (varied) bifurcation parameter of interest
    time_x_parameter_matrix = np.transpose(parameter_x_time_matrix)

    # Slice the iterations to display to a smaller range
    if time_slice is not None:
        time_x_parameter_matrix = time_x_parameter_matrix[time_slice]

    # Make bifurcation plot
    fig, ax = plt.subplots()

    # For the solutions vector at timestep plot the bifurcation
    # NOTE: The elements of the solutions vector represent the 
    # numerical solutions at timestep t for all varied parameters
    # in the parameter space
    # e.g.,
    # t  beta1=0.025     beta1=0.030   ....   beta1=0.30
    # 0  solution00      solution01    ....   solution0P
    for sol_at_time_t_for_all_params in time_x_parameter_matrix:

        if log_scale:
            sol_at_time_t_for_all_params = np.log(sol_at_time_t_for_all_params)

        ax.plot(
            bifurcation_parameter_space, 
            sol_at_time_t_for_all_params,
            ',k',
            alpha=0.25)

    return fig, ax

# Define initial conditions based on figure
s0 = 6e-2
e0 = i0 = 1e-3
initial_conditions = [s0, e0, i0]

# Define model parameters based on figure
# NOTE: omega is not mentioned in the figure, but 
# omega is defined elsewhere as 2pi/365
days_per_year = 365

mu = 0.02/days_per_year
beta_zero = 1250
sigma = 1/8
gamma = 1/5
omega = 2*np.pi / days_per_year

model_params = dict(
        beta_zero=beta_zero,
        omega=omega,
        mu=mu,
        sigma=sigma,
        gamma=gamma)

# Define timesteps
nyears = 200
ndays = nyears * days_per_year
timesteps = np.arange(1, ndays + 1, 1)

# Define different levels of seasonality (from figure)
beta_ones = [0.025, 0.05, 0.25]

# Define the time range to actually show on the plot
min_year = 190
max_year = 200

# Create a slice of the iterations to display on the diagram
time_slice = slice(min_year*days_per_year, max_year*days_per_year)

# Get the xticks to display on the plot based on the time slice
x_ticks = timesteps[time_slice]/days_per_year

# Plot the panels using the infected state variable ix
infection_ix = 2

# Plot the panels 
panel_fig, panel_ax = plot_panels(
    model=seasonal_seir,
    model_params=model_params,
    panel_param_space=beta_ones,
    panel_param_name='beta_one',

    initial_conditions=initial_conditions,
    timesteps=timesteps,
    odeint_kwargs=dict(hmax=5),

    x_ticks=x_ticks,
    time_slice=time_slice,
    state_var_ix=infection_ix,
    log_scale=False)

# Label the panels
panel_fig.suptitle('Attempt to Reproduce Panels from Keeling 2007')
panel_fig.supxlabel('Time (years)')
panel_fig.supylabel('Fraction Infected')
panel_fig.set_size_inches(15, 8)

# Plot bifurcation
bi_fig, bi_ax = plot_bifurcation(
    model=seasonal_seir,
    model_params=model_params,
    bifurcation_parameter_space=np.linspace(0.025, 0.3),
    bifurcation_param_name='beta_one',

    initial_conditions=initial_conditions,
    timesteps=timesteps,
    odeint_kwargs={'hmax':5},

    state_var_ix=infection_ix,
    time_slice=time_slice,
    log_scale=False)

# Label the bifurcation
bi_fig.suptitle('Attempt to Reproduce Bifurcation Diagram from Keeling 2007')
bi_fig.supxlabel(r'$\beta_1$')
bi_fig.supylabel('Fraction Infected')
bi_fig.set_size_inches(15, 8)
Improve this question
$\endgroup$
4
  • $\begingroup$ Can you check the value of $\beta_0=1250$ again in your sources, if it appears in other places than the caption of the text. This value appears to be rather large, is there a decimal dot missing, for a value of $12.50$? $\endgroup$
    – Lutz Lehmann
    Oct 10, 2022 at 10:21
  • $\begingroup$ Changing $\beta_0$ to $1.250$ as well as $12.50$ does not lead to plots that match the figure in the textbook. The textbook only mentions $\beta_0=1250$ for this particular plot and it doesn't look like $\beta_0$ appears in other places. $\endgroup$
    – Jared Frazier
    Oct 10, 2022 at 13:19
  • $\begingroup$ At a second glance, as the component $I$ is rather small, it is not unreasonable that the coefficient can be large. A test would be if for $\beta_1=0$ the then easily computable stationary points are in the range of the plots. $\endgroup$
    – Lutz Lehmann
    Oct 10, 2022 at 13:22
  • $\begingroup$ No, at $β_0=1250$ the diagram is a simple line, the solutions are strictly sinusoid as your upper plots show, the sections at the max of the cosine have all the same value. In our previous discussion, the most "chaotic" $R_0=17$ corresponds to $β_0=R_0·\gamma=3.4$. And indeed, a plot similar to the published one is obtained for $β_0=3.6$, generally $3.5$ to $5.5$, with a rapid reduction in features around that. $\endgroup$
    – Lutz Lehmann
    Oct 10, 2022 at 16:42

1 Answer 1

Reset to default
6
$\begingroup$

Recently in https://math.stackexchange.com/questions/4542008/how-to-loop-parameter-a-in-henon-map I came into contact with the idea of an "adiabatically gliding" parameter where one gets the full bifurcation diagram in one integration.

Applying this to the current problem would take the maximal time $T$ of the time steps and set in the ODE function $\beta_1(t)=\beta_{1,\max}\frac{t}{T}$. The longer the time span the slower the changes in $\beta_1$.

Note that the plot is kind of a Poincaré section, it is a scatter plot of the value at the same day of the year, not a continuous plot.

This gives bifurcating plots of the given type for $\beta_0$ roughly between 3 and 7, with the closest plots around $\beta_0=3.6$

enter image description here

This is not horribly precise, especially where the "equilibrium cycle" changes rapidly or qualitatively the process may not remain long enough in the region to reach a new stable behavior.

I've shortened the code to the essential, removing the interfaces for parametric plot production and documentation.

nyears = 10000
days_per_year = 365
timesteps = np.arange(0, nyears+1, 1) * days_per_year
T = timesteps[-1]

# Model parameters (see Figure 5.6 description)
mu = 0.02 / days_per_year
sigma = 1/8
gamma = 1/5
omega = 2 * np.pi / days_per_year  # frequency of oscillations per year

beta_zero = 3.6
beta_one = lambda t: t/T*0.35


def seir(y, t):
    s, e, i = y
    beta = beta_zero * (1 + beta_one(t) * np.cos(omega * t))
    
    sdot = mu - (beta*i + mu)*s
    edot = beta*s*i - (mu  + sigma)*e
    idot = sigma*e - (mu + gamma)*i
    return sdot, edot, idot

# sigma*beta*s*i = sigma*(mu  + sigma)*e = (mu  + sigma)*(mu + gamma)*i
beta = beta_zero
s0 = (mu  + sigma)*(mu + gamma)/(sigma*beta)
i0 = mu*(1-s0)/(beta*s0)
e0 = (mu + gamma)*i0/sigma
y0 = [s0, e0, i0]
print(y0,seir(y0,0))

y = odeint(seir, y0, timesteps, atol=1e-9, rtol=1e-9, hmax = 5)
s,e,i = y.T

%matplotlib inline
plt.figure(figsize=(12,5))
plt.semilogy(beta_one(timesteps),i,'o', ms=1)
plt.ylim(1e-18,1e-1)
plt.title(f"beta_zero = {beta_zero}")
plt.xlabel("beta_one")
plt.ylabel("infected")
plt.show()

A proper diagram

In the end, it is not too much difference to make a more correct plot. Instead of the slow-moving $\beta_1$, keep it constant, do some iterations to stabilize, then about the same amount to plot. Switch to the next $\beta_1$, use the last state vector as initial point.

The resulting plot looks much more clear

enter image description here

The change in the code is mostly the loop over the single integration runs.

nyears = 200
days_per_year = 365
timesteps = np.arange(0, nyears+1, 1) * days_per_year

# Model parameters (see Figure 5.6 description)
mu = 0.02 / days_per_year
sigma = 1/8
gamma = 1/5
omega = 2 * np.pi / days_per_year  # frequency of oscillations per year

beta_zero = 3.6
beta_ones = np.linspace(0,0.35,600)


def seir(y, t, beta_one):
    s, e, i = y
    beta = beta_zero * (1 + beta_one * np.cos(omega * t))
    
    sdot = mu - (beta*i + mu)*s
    edot = beta*s*i - (mu  + sigma)*e
    idot = sigma*e - (mu + gamma)*i
    return sdot, edot, idot

def seir_log(y,t,*args):
    y = np.exp(y)
    dy = np.array(seir(y,t,*args))
    return dy/y

# sigma*beta*s*i = sigma*(mu  + sigma)*e = (mu  + sigma)*(mu + gamma)*i
beta = beta_zero
s0 = (mu  + sigma)*(mu + gamma)/(sigma*beta)
i0 = mu*(1-s0)/(beta*s0)
e0 = (mu + gamma)*i0/sigma
y0 = [s0, e0, i0]
print(y0,seir(y0,0,0))

b1s, inf = [],[]
for beta_one in beta_ones:
    y = odeint(seir_log, np.log(y0), timesteps, args = (beta_one,),
               atol=1e-9, rtol=1e-9, hmax = 5)
    b1s.append(100*[beta_one])
    
    s,e,i = np.exp(y.T)
    y0 = np.exp(y[-1])
    inf.append(i[-100:])
    
b1s = np.concatenate(b1s)
inf = np.concatenate(inf)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Really cool solution! A few things: (1) Why do you slice the infected solutions with i[-100:]? (2) Why are you using the last state vector computed in the previous iteration (i.e., y0=np.exp(y[-1])) for the initial conditions of numerical integration? (3) It is really frustrating that the author clearly is using parameters/techniques that are not mentioned in the figure. It's discouraging as a learner since it shifts the focus from trying to reproduce results using the parameters and description to guessing what the parameters are or math tricks that the author uses. Any advice? $\endgroup$
    – Jared Frazier
    Oct 10, 2022 at 20:32
  • $\begingroup$ Incidentally, I asked this question on stack overflow also, so you should go to this question and just answer with a link to your answer here on Computational Science. That way I can credit you properly and close my question on stack overflow. $\endgroup$
    – Jared Frazier
    Oct 10, 2022 at 20:35
  • 1
    $\begingroup$ (1) The first 100 years serve to stabilize the cycle, if any, then the last 100 years are used for the plot. (2) In most cases, the last value will be very close to the cycle in the new parameter. Thus faster stabilization. (3) They might be using a home-brew integrator, such as fixed-step RK4. There could have been a chain of errors, the author asks the computer guru to implement the system, the guru delegates the documentation of the methods used to some student, details get lost, ... $\endgroup$
    – Lutz Lehmann
    Oct 10, 2022 at 20:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.