Finding ODE solutions

Finding ODE solutions#

PyGOM allows the user to evaluate both the deterministic and stochastic time evolution of their ODE system using the class methods solve_determ() and solve_stochast() respectively. These methods work with both fixed and random parameters as introduced in the previous section.

We begin by defining the series of ODEs and parameters which define our SIR system. This we do from scratch rather than loading in a previously defined model in order to present a more comprehensive example of the workflow.

from pygom import SimulateOde, Transition, TransitionType
import matplotlib.pyplot as plt
import numpy as np
import random

###################
# ODE specification
###################

# Define SIR model
stateList = ['S', 'I', 'R']
paramList = ['beta', 'gamma', 'N']
transitionList = [Transition(origin='S', destination='I', equation='beta*S*I/N', transition_type=TransitionType.T),
                  Transition(origin='I', destination='R', equation='gamma*I', transition_type=TransitionType.T)]

n_pop=1e4  # Total population is fixed

############
# Parameters
############

beta_mn=0.35   # Infectivity, beta. Gives the actual value for fixed params and mean for random distribution.
gamma_mn=0.25  # Recovery rate, gamma.

#######
# Fixed
#######
fixed_param_set=[('beta', beta_mn), ('gamma', gamma_mn), ('N', n_pop)]

########
# Random
########

# Recovery rate, gamma
gamma_var=(gamma_mn/10)**2  # Set the standard deviation to be 1/10th of the mean value
gamma_shape=(gamma_mn**2)/gamma_var
gamma_rate=gamma_mn/gamma_var

# Infectivity parameter, beta
beta_var=(beta_mn/10)**2  # Set the standard deviation to be 1/10th of the mean value
beta_shape=(beta_mn**2)/beta_var
beta_rate=beta_mn/beta_var

from pygom.utilR import rgamma
random_param_set = dict()  # container for random param set
random_param_set['gamma'] = (rgamma,{'shape':gamma_shape, 'rate':gamma_rate})
random_param_set['beta'] = (rgamma,{'shape':beta_shape, 'rate':beta_rate})
random_param_set['N'] = n_pop

Since this notebook will involve stochastic processes, we set the random number generator seed to make outputs reproducible.

np.random.seed(1)

In order to determine the time evolution of the system, we must supply initial conditions as well as the desired time points for the numerical solver. Timesteps should be sufficiently short to reduce numerical integration errors, but not too short such that the computational time costs become too large.

import math

# Initial conditions
i0=10
x0 = [n_pop-i0, i0, 0]

# Time range and increments
tmax=200                              # maximum time over which to run solver
dt=0.1                                # timestep
n_timestep=math.ceil(tmax/dt)         # number of iterations
t = np.linspace(0, tmax, n_timestep)  # times at which solution will be evaluated

Deterministic evolution#

To solve for the deterministic time evolution of the system, PyGOM uses scipy.integrate.odeint() which is wrapped by the member function solve_determ(). We begin with the simple (and likely familiar) case of fixed parameters.

Fixed parameters#

First, we initialise a SimulateOde object with our fixed parameters, fixed_param_set:

# Set up pygom object (D_F suffix implies Deterministic_Fixed)
ode_D_F = SimulateOde(stateList, paramList, transition=transitionList)
ode_D_F.initial_values = (x0, t[0])  # (initial state conditions, initial timepoint)
ode_D_F.parameters=fixed_param_set

The solution is then found via solve_determ, specifying the required time steps (not including the initial one).

solution_D_F = ode_D_F.solve_determ(t[1::])

Plotting the output yields the familiar result, where infecteds initially increase in number exponentially until critical depletion of susceptibles results in epidemic decline.

Hide code cell source 
f, axarr = plt.subplots(1,3, layout='constrained', figsize=(10, 2.5))

# Plot colours
colours=["C1", "C0", "C2"]

for i in range(0, 3):
    axarr[i].plot(t, solution_D_F[:,i], color=colours[i])

for idx, state in enumerate(stateList):
    axarr[idx].set_ylabel(state, rotation=0)
    axarr[idx].set_xlabel('Time')

axarr[1].set_title("Deterministic simulation with fixed parameters")

plt.show()
../_images/00c46ab2830b10180e9d261072fdf576ce9729a160a06520f45813448dec77de.png

Random parameters#

We now solve the same system, but for 1000 repetitions using randomly drawn parameters for each simulation. This time we initialise the parameters withrandom_param_set, but still use solve_determ() to find solutions as before

# Set up pygom object
ode_D_R = SimulateOde(stateList, paramList, transition=transitionList)
ode_D_R.initial_values = (x0, t[0])
ode_D_R.parameters=random_param_set

n_param_draws=1000  # number of parameters to draw
Ymean, solution_D_R = ode_D_R.solve_determ(t[1::], n_param_draws, full_output=True)

Note

A message may be printed above where PyGOM is trying to connect to an mpi backend, as our module has the capability to compute in parallel using the IPython.

Here we visualise the output in 2 ways, first by viewing 50 randomly selected trajectories and secondly by viewing the confidence intervals (here 95% and 50%) and median calculated over the full 1000 solutions.

Hide code cell source 
y_D_R=np.dstack(solution_D_R)  # unpack the data

#########################
# Individual trajectories
#########################

# Select 50 simulations to plot
i_rand=random.sample(range(n_param_draws), 50)

######################
# Confidence intervals
######################

# Calculate 95%, 50% CIs and median.
y_D_R_lolo=np.percentile(y_D_R, 2.5, axis=2)
y_D_R_lo=np.percentile(y_D_R, 25, axis=2)
y_D_R_hi=np.percentile(y_D_R, 75, axis=2)
y_D_R_hihi=np.percentile(y_D_R, 97.5, axis=2)
y_D_R_md=np.percentile(y_D_R, 50, axis=2)

f, axarr = plt.subplots(2,3, layout='constrained', figsize=(10, 5))

# Plot colours
colours=["C1", "C0", "C2"]

for i in range(0,3):
    # Plot individual trajectories
    for j in i_rand:
        axarr[0][i].plot(t, y_D_R[:,i,j], color=colours[i], alpha=0.2)

    # Plot CI's
    axarr[1][i].fill_between(t, y_D_R_lolo[:,i], y_D_R_hihi[:,i], alpha=0.2, facecolor=colours[i])
    axarr[1][i].fill_between(t, y_D_R_lo[:,i], y_D_R_hi[:,i], alpha=0.4, facecolor=colours[i])
    axarr[1][i].plot(t, y_D_R_md[:,i], color=colours[i])

# Add titles
for idx, state in enumerate(stateList):
    axarr[0][idx].set_ylabel(state, rotation=0)
    axarr[1][idx].set_ylabel(state, rotation=0)
    axarr[0][idx].set_xlabel('Time')
    axarr[1][idx].set_xlabel('Time')

axarr[0][1].set_title("50 deterministic simulations, each with randomly drawn parameters")
axarr[1][1].set_title("Median (line), 50% CI (dark shaded) and 95% CI (light shaded) over 1000 simulations")

plt.show()
../_images/3129dbc62512581475b52e955ea4a59bf4ae67887a16f501648e1713d14c38d7.png

Stochastic evolution#

The approximation that numbers of individuals in each state may be treated as a continuum break down when their sizes are small. In this regime, transitions between states do not represent continuous flows but are instead stochastic events that occur at rates governed by the current state of the system. The simplifying assumption that waiting times for these events to occur are exponentially distributed (memoryless), allows for quicker evaluation of the dynamics.

Two common algorithms have been implemented for use during simulation; the reaction method 5 and the \(\tau\)-Leap method 2. The two change interactively depending on the size of the states.

Fixed parameters#

As previously, we define a model and pass our fixed parameters fixed_param_set. However, this time we employ the function solve_stochast to allow for stochastic time evolution:

# Set up pygom object
ode_S_F = SimulateOde(stateList, paramList, transition=transitionList)
ode_S_F.initial_values = (x0, t[0])

n_sim=1000  # number of simulations
ode_S_F.parameters = fixed_param_set
solution_S_F, simJump, simT = ode_S_F.solve_stochast(t, n_sim, full_output=True)

Illegal jump, x: [4.769e+03 1.000e+00 5.230e+03], new x: [ 4.769e+03 -1.000e+00  5.232e+03]

Illegal jump, x: [5.032e+03 1.000e+00 4.967e+03], new x: [ 5.032e+03 -1.000e+00  4.969e+03]

Before we inspect the epidemic time series, we plot the distribution of final epidemic sizes.

Hide code cell source 
y_S_F=np.dstack(solution_S_F)  # unpack the data
final_size_S_F=y_S_F[-1, 2, :]

plt.figure(figsize=(8, 4))
plt.title("Epidemic final size distribution")
plt.xlabel("Final size")
plt.ylabel("Count")
plt.hist(final_size_S_F, bins=100, color="C0")
plt.show()
../_images/2a3a06043a37cce1023acbde4ac475a1b2cea0c440f7c837d33a86dc8f67b30a.png

We see that 28 of the 1000 simulations (2.8%) have final sizes very close to zero. In these cases, the chain of infections reaches stochastic extinction before an epidemic is able to take off, as can be seen when plotting a selection of these failed epidemics.

Hide code cell source 
# Plot failed epidemics
i_success_S_F=np.where(final_size_S_F>2000)
y_fail_S_F=np.delete(y_S_F, obj=i_success_S_F, axis=2) 
f, ax = plt.subplots(figsize=(8, 4))

ax.set_title("Stochastically extinct epidemic trajectories")
ax.set_xlabel("Time")
ax.set_ylabel("Infected count")
for i in range(min(y_fail_S_F.shape[2], 3)):
    ax.plot(t[0:300], y_fail_S_F[0:300,1,i], color="C0", alpha=0.5)
../_images/3ab873e59da4579d18d73e5ee674421f00cdb916582603532f4e3f8602d93d57.png

We proceed to eliminate these when plotting our simulation outputs, since we’d like to compare the output here with our deterministic results (which, by definition, do not exhibit stochastic extinction):

Hide code cell source 
# Delete failed epidemics
i_fail_S_F=np.where(final_size_S_F<2000)
y_success_S_F=np.delete(y_S_F, obj=i_fail_S_F, axis=2) 

#########################
# Individual trajectories
#########################

# Select 50 simulations to plot
i_rand=random.sample(range(y_success_S_F.shape[2]), 50)

######################
# Confidence intervals
######################

# Calculate 95%, 50% CIs and median.
y_S_F_lolo=np.percentile(y_success_S_F, 2.5, axis=2)
y_S_F_lo=np.percentile(y_success_S_F, 25, axis=2)
y_S_F_hi=np.percentile(y_success_S_F, 75, axis=2)
y_S_F_hihi=np.percentile(y_success_S_F, 97.5, axis=2)
y_S_F_md=np.percentile(y_success_S_F, 50, axis=2)


f, axarr = plt.subplots(2,3, layout='constrained', figsize=(10, 5))

# Plot colours
colours=["C1", "C0", "C2"]

for i in range(0,3):
    # Plot individual trajectories
    for j in i_rand:
        axarr[0][i].plot(t, y_success_S_F[:,i,j], color=colours[i], alpha=0.2)

    # Plot CI's
    axarr[1][i].fill_between(t, y_S_F_lolo[:,i], y_S_F_hihi[:,i], alpha=0.2, facecolor=colours[i])
    axarr[1][i].fill_between(t, y_S_F_lo[:,i], y_S_F_hi[:,i], alpha=0.4, facecolor=colours[i])
    axarr[1][i].plot(t, y_S_F_md[:,i], color=colours[i])

# Add titles
for idx, state in enumerate(stateList):
    axarr[0][idx].set_ylabel(state, rotation=0)
    axarr[1][idx].set_ylabel(state, rotation=0)
    axarr[0][idx].set_xlabel('Time')
    axarr[1][idx].set_xlabel('Time')

axarr[0][1].set_title("50 stochastic simulations with shared parameters")
axarr[1][1].set_title("Median (line), 50% CI (dark shaded) and 95% CI (light shaded) over 1000 simulations")

plt.show()
../_images/c5fe55afc9c8201692791a430e3b954d9a25e60c734f5bb107039ceed1e64b3d.png

Random parameters#

Finally, we explore a scenario where we wish to solve for a stochastic system with parameters drawn randomly from some distributions. Using the random parameter set, random_param_set, defined previously

# Set up pygom object
ode_S_R = SimulateOde(stateList, paramList, transition=transitionList)
ode_S_R.initial_values = (x0, t[0])

# In this case n_sim and n_param_draw mean the same thing since we have new parameter set per simulation
ode_S_R.parameters = random_param_set
solution_S_R, simJump, simT = ode_S_R.solve_stochast(t, n_sim, full_output=True)

Again, we begin our analysis by inspecting the final size distribution.

Hide code cell source 
y_S_R=np.dstack(solution_S_R)

final_size_S_R=y_S_R[-1, 2, :]
plt.figure(figsize=(8, 4))
plt.title("Epidemic final size distribution")
plt.xlabel("Final size")
plt.ylabel("Count")
plt.hist(final_size_S_R, bins=100)
plt.show()
../_images/18e35efafc47712a765168c6a80e322f19387ae25ffccc9ec8e2b3ea071aecd0.png

We see that the final size distribution is much more spread out, owing to increased variation in the parameters. The threshold value separating failed and successful epidemics is less clear cut in this case and we judge this by eye to be 2000. Again, eliminating the failed epidemics we plot the model output

Hide code cell source 
# Delete failed epidemics
i_fail_S_R=np.where(final_size_S_R<2000)
y_success_S_R=np.delete(y_S_R, obj=i_fail_S_R, axis=2) 

#########################
# Individual trajectories
#########################

# Select 50 simulations to plot
i_rand=random.sample(range(y_success_S_R.shape[2]), 50)

######################
# Confidence intervals
######################

# Calculate 95%, 50% CIs and median.
y_S_R_lolo=np.percentile(y_success_S_R, 2.5, axis=2)
y_S_R_lo=np.percentile(y_success_S_R, 25, axis=2)
y_S_R_hi=np.percentile(y_success_S_R, 75, axis=2)
y_S_R_hihi=np.percentile(y_success_S_R, 97.5, axis=2)
y_S_R_md=np.percentile(y_success_S_R, 50, axis=2)

f, axarr = plt.subplots(2,3, layout='constrained', figsize=(10, 5))

# Plot colours
colours=["C1", "C0", "C2"]

for i in range(0,3):
    # Plot individual trajectories
    for j in i_rand:
        axarr[0][i].plot(t, y_success_S_R[:,i,j], color=colours[i], alpha=0.2)

    # Plot CI's
    axarr[1][i].fill_between(t, y_S_R_lolo[:,i], y_S_R_hihi[:,i], alpha=0.2, facecolor=colours[i])
    axarr[1][i].fill_between(t, y_S_R_lo[:,i], y_S_R_hi[:,i], alpha=0.4, facecolor=colours[i])
    axarr[1][i].plot(t, y_S_R_md[:,i], color=colours[i])

# Add titles
for idx, state in enumerate(stateList):
    axarr[0][idx].set_ylabel(state, rotation=0)
    axarr[1][idx].set_ylabel(state, rotation=0)
    axarr[0][idx].set_xlabel('Time')
    axarr[1][idx].set_xlabel('Time')

axarr[0][1].set_title("50 stochastic simulations, each with randomly drawn parameters")
axarr[1][1].set_title("Median (line), 50% CI (dark shaded) and 95% CI (light shaded) over 1000 simulations")

plt.show()
../_images/39d3ea9b3f55681166595b1cdcd4fb8760af2d8b6abcc4eb4243856eeef1a6bb.png

Summary#

With two possible simulation methods and two possible parameter types, we have explored the four total model configurations. We plot these results side by side for the infected compartment to provide an illustration of their key qualitative differences.

Hide code cell source 
ymax=1250

fig, axarr2 = plt.subplots(2,2, layout='constrained', figsize=(10, 6))

axarr2[0][0].plot(t, solution_D_F[:,1], color="C0")
axarr2[0][0].set_ylim(bottom=0, top=ymax)

i_rand=random.sample(range(n_param_draws), 50)

for i in i_rand:
    axarr2[0][1].plot(t, y_S_F[:,1,i], color="C0", alpha=0.3)
    axarr2[1][0].plot(t, y_D_R[:,1,i], color="C0", alpha=0.3)
    axarr2[1][1].plot(t, y_S_R[:,1,i], color="C0", alpha=0.3)

axarr2[0][0].set_title("Deterministic with fixed params")
axarr2[0][1].set_title("Stochastic with fixed params")
axarr2[1][0].set_title("Deterministic with random params")
axarr2[1][1].set_title("Stochastic with random params")

axarr2[0][0].set_xlabel("Time")
axarr2[0][1].set_xlabel("Time")
axarr2[1][0].set_xlabel("Time")
axarr2[1][1].set_xlabel("Time")

axarr2[0][1].set_ylim(bottom=0, top=ymax)
axarr2[1][0].set_ylim(bottom=0, top=ymax)
axarr2[1][1].set_ylim(bottom=0, top=ymax)
plt.show()
../_images/05ae23eb5feab2c1e23864f56bf6a26d9960a8622aa20e923e61da80abc8b708.png
Contents