Modeling the spread of rabies in one dimension¶

In this notebook, we will move from the 2 connected patches described in Notebook 8 to a 1-D grid of connected population patches. This is the scenario explored to describe the spatial spread of rabies in foxes in the work of Kallen, Arcuri, and Murray, Journal of Theoretical Biology (1985) 116, 377-393. The relevant model equations (note that these are now PDEs rather than the ODEs of previous notebooks) are:

$$ \frac{\partial S}{\partial t} = -KIS $$

$$ \frac{\partial I}{\partial t} = D \frac{\partial^2 I}{\partial x^2} + KIS - \mu I $$

and come to the conclusion that this system supports traveling waves with velocity bounded below by

$$ c = 2 \sqrt{D \mu (1/r -1)} $$

where $r = \frac{\mu}{K S_0}$ and $S_0$ is the initial density of susceptibles.

The analogy to the SIR model equations we have been using is clear, though this model conceptually different in specific ways from that model as we now have spatially dimensional quantities. For example, I & S are infective/susceptible population densities rather than population counts, which affects translating between the $K$ of this model and the $\beta$ of the SIR; the derived quantity $r$ is clearly analogous to $\frac{1}{R_0}$ of the SIR model, but requires the presence of the initial susceptible density to be properly dimensionless, and rabies being invariably fatal, $\mu$ represents the mortality of rabies rather than the $\gamma$ we would use for recovery.

With all that said, if we are careful about translation between parameters, and appropriately construct the transmission network to appropriately reflect 1-D diffusion, we ought to be able to reproduce this model in our SIR implementation, and recapitulate the traveling wave behavior of the outbreak. As we do this, we will also run into two well-understood issues in simulating traveling waves of the expected wavespeed on a finite lattice, both biasing us to produce slower waves on average.

The first comes from discretization - we are transitioning from a PDE continuous in space, time, and the relevant variable $I$, to a simulation that is discrete in all 3. This discretization will make it impossible for infection to travel more than 1 lattice spacing per simulation timestep, so in parameter regimes where the expected wavespeed approaches a lattice spacing per simulation timestep, our model will be incapable of demonstrating that.

The second slowdown ill come from discretization of the population itself, and from stochasticity. This phenomenon is sometimes referred to as the "attofox" problem, and is discussed, e.g., in the chapter on Spatial Modeling in Keeling and Rohani, Modeling Infectious Diseases in Humans and Animals. In a model with deterministic dynamics and continuous populations, any fraction of infectivity deposited in a population can generate fractional infectives, and deterministically kick off the start of a local outbreak. In a discrete population, stochastic model, we cannot create new infections in smaller units than a single infected individual. This causes stochastic delays relative to the continuous model - these delays can arise both in exportation - it can take longer to infect an individual in an uninfected location, and it can take longer for the outbreak to then progress from 1 individual to more, and the wavespeed will on average be slower.

In fact, as long as $D>0$ and $1/r >1$, the continuous equation is guaranteed to support traveling waves, but for small D and small 1/r the stochastic version may simply die out in the initial node before successfully propagating to adjacent nodes, or it may die out at some subsequent node along the way.

So as we fit this result, we should expect to see the observed wavespeed consistently be slower than expected at very high propagation speeds. At lower propagation speeds, the match between model and theory will vary with the values of D, $\beta$, $\gamma$ - we should expect better match between theory and observation but with a general bias for observation to be slower.

In [1]:

Copied!





import numpy as np
import pandas as pd
from laser.core.propertyset import PropertySet
import laser.core.distributions as dists
from laser.generic import SIR
from laser.generic import Model
from laser.core.utils import grid
import laser.core
import laser.generic
import matplotlib.pyplot as plt
import os
from scipy.optimize import fsolve

print(f"{np.__version__=}")
print(f"{laser.core.__version__=}")
print(f"{laser.generic.__version__=}")
import numpy as np
import pandas as pd
from laser.core.propertyset import PropertySet
import laser.core.distributions as dists
from laser.generic import SIR
from laser.generic import Model
from laser.core.utils import grid
import laser.core
import laser.generic
import matplotlib.pyplot as plt
import os
from scipy.optimize import fsolve

print(f"{np.__version__=}")
print(f"{laser.core.__version__=}")
print(f"{laser.generic.__version__=}")

np.__version__='2.2.6'
laser.core.__version__='0.9.1'
laser.generic.__version__='0.0.0'

In [2]:

Copied!





#Set up the network and initial conditions
npatches = 201
pop = 1000
scenario = grid(M=1, N=npatches, node_size_degs=0.05, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 3
scenario["I"] = 0
scenario["R"] = 0
scenario["S"] = scenario.population
scenario.loc[scenario.nodeid==101,"I"] = initial_infected
scenario.loc[scenario.nodeid==101,"S"] = scenario.loc[scenario.nodeid==101,"population"]-initial_infected
#Set up the network and initial conditions
npatches = 201
pop = 1000
scenario = grid(M=1, N=npatches, node_size_degs=0.05, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 3
scenario["I"] = 0
scenario["R"] = 0
scenario["S"] = scenario.population
scenario.loc[scenario.nodeid==101,"I"] = initial_infected
scenario.loc[scenario.nodeid==101,"S"] = scenario.loc[scenario.nodeid==101,"population"]-initial_infected

The transmission network¶

Spatial transmission in this rabies model is contained in the diffusive term $D \frac{\partial^2 I}{\partial x^2}$. We are operating on a discretized 1-D grid, and will take advantage of the usual way to discretize a second derivative

$$ D \frac{\partial^2 I}{\partial x^2} \rightarrow D \frac {I_{x-h} - 2 I_{x} + I_{x+h}} {h^2} $$

To translate this into a network structure for the model, we can arbitrarily set the length scale $h$ to be 1 in units of the grid spacing. Looking at the above, from the perspective of an arbitrary node $i$, the transmission matrix must "send" $D I_i$ units of infectivity to the adjacent nodes on either side, and remove that infectivity from the local pool. Looking at the spatial transmission code (line 1006-8 in src/laser_generic/models/components.py), we see that spatial transmission is handled by multiplying the transmission matrix and the vector of infectivities by node, and that it internally handles subtracting the "outbound infectivity" from the local source. Thus, the above discretization can be accomplished by filling the value $D$ on the first superdiagonal and subdiagonal of the transmission matrix, and 0 everywhere else.

In [3]:

Copied!





infmean = 10.0
R0 = 10.0
nticks = 2 * 365
D = .01
parameters = PropertySet(
    {
        "seed": np.random.randint(0, 1000000),
        "nticks": nticks,
        "verbose": True,
        "beta": R0 / infmean,
        "inf_mean": infmean,
    }
)

infdurdist = dists.exponential(scale=parameters.inf_mean)
model = Model(scenario, parameters)
model.network = 0*model.network
np.fill_diagonal(model.network[1:], D)
np.fill_diagonal(model.network[:,1:], D)

model.components = [
    SIR.Susceptible(model),
    SIR.Recovered(model),
    SIR.Infectious(model, infdurdist),
    SIR.Transmission(model, infdurdist),
]
model.run()
infmean = 10.0
R0 = 10.0
nticks = 2 * 365
D = .01
parameters = PropertySet(
    {
        "seed": np.random.randint(0, 1000000),
        "nticks": nticks,
        "verbose": True,
        "beta": R0 / infmean,
        "inf_mean": infmean,
    }
)

infdurdist = dists.exponential(scale=parameters.inf_mean)
model = Model(scenario, parameters)
model.network = 0*model.network
np.fill_diagonal(model.network[1:], D)
np.fill_diagonal(model.network[:,1:], D)

model.components = [
    SIR.Susceptible(model),
    SIR.Recovered(model),
    SIR.Infectious(model, infdurdist),
    SIR.Transmission(model, infdurdist),
]
model.run()

201,000 agents in 201 node(s): 100%|██████████| 730/730 [00:00<00:00, 1314.70it/s]

In [4]:

Copied!





plt.figure(figsize=(9, 4))
plt.imshow(model.nodes.I.T, aspect='auto', origin='lower')
plt.xlabel('Time (tick)')
plt.ylabel('Patch')
plt.title('Infectious (I) over time and space')
plt.colorbar(label='Infectious')
plt.show()
plt.figure(figsize=(9, 4))
plt.imshow(model.nodes.I.T, aspect='auto', origin='lower')
plt.xlabel('Time (tick)')
plt.ylabel('Patch')
plt.title('Infectious (I) over time and space')
plt.colorbar(label='Infectious')
plt.show()

No description has been provided for this image

In [5]:

Copied!





from scipy.stats import linregress

# Find the time (tick) of peak I for each node
peak_times = model.nodes.I.argmax(axis=0)  # shape: (npatches,)

# Get the node indices for right-going (101 to 200) and left-going (100 to 0) waves
right_nodes = np.arange(101, npatches)
left_nodes = np.arange(100, -1, -1)

# Fit a line to peak_times vs node index for each direction to estimate velocity (slope = 1/velocity)

# Right-going wave
right_peak_times = peak_times[right_nodes]
slope_r, intercept_r, r_value_r, p_value_r, std_err_r = linregress(right_nodes[right_peak_times>0], right_peak_times[right_peak_times > 0])
velocity_right = 1 / slope_r if slope_r != 0 else np.nan

# Left-going wave
left_peak_times = peak_times[left_nodes]
slope_l, intercept_l, r_value_l, p_value_l, std_err_l = linregress(left_nodes[left_peak_times>0], left_peak_times[left_peak_times > 0])
velocity_left = -1 / slope_l if slope_l != 0 else np.nan  # negative because left

print(f"Estimated right-going wavespeed: {velocity_right:.3f} patches/tick")
print(f"Estimated left-going wavespeed: {velocity_left:.3f} patches/tick")
print(f"Expected wavespeed: {2*np.sqrt(D/infmean*(R0-1)):.3f} patches/tick")
from scipy.stats import linregress

# Find the time (tick) of peak I for each node
peak_times = model.nodes.I.argmax(axis=0)  # shape: (npatches,)

# Get the node indices for right-going (101 to 200) and left-going (100 to 0) waves
right_nodes = np.arange(101, npatches)
left_nodes = np.arange(100, -1, -1)

# Fit a line to peak_times vs node index for each direction to estimate velocity (slope = 1/velocity)

# Right-going wave
right_peak_times = peak_times[right_nodes]
slope_r, intercept_r, r_value_r, p_value_r, std_err_r = linregress(right_nodes[right_peak_times>0], right_peak_times[right_peak_times > 0])
velocity_right = 1 / slope_r if slope_r != 0 else np.nan

# Left-going wave
left_peak_times = peak_times[left_nodes]
slope_l, intercept_l, r_value_l, p_value_l, std_err_l = linregress(left_nodes[left_peak_times>0], left_peak_times[left_peak_times > 0])
velocity_left = -1 / slope_l if slope_l != 0 else np.nan  # negative because left

print(f"Estimated right-going wavespeed: {velocity_right:.3f} patches/tick")
print(f"Estimated left-going wavespeed: {velocity_left:.3f} patches/tick")
print(f"Expected wavespeed: {2*np.sqrt(D/infmean*(R0-1)):.3f} patches/tick")

Estimated right-going wavespeed: 0.179 patches/tick
Estimated left-going wavespeed: 0.180 patches/tick
Expected wavespeed: 0.190 patches/tick

In [6]:

Copied!





# Set up parameter ranges
n_samples = 100
D_samples = np.logspace(np.log10(3e-4), np.log10(3e-1), n_samples)
R0_samples = np.random.uniform(1.5, 10, n_samples)
infmean_samples = np.random.uniform(7, 50, n_samples)

results = []

for i, (D, R0, infmean) in enumerate(zip(D_samples, R0_samples, infmean_samples)):
    # Set up parameters for this run
    parameters = PropertySet(
        {
            "seed": np.random.randint(0, 1000000),
            "nticks": nticks,
            "verbose": False,
            "beta": R0 / infmean,
            "inf_mean": infmean,
        }
    )
    infdurdist = dists.exponential(scale=parameters.inf_mean)
    model = Model(scenario.copy(), parameters)
    model.network = 0 * model.network
    np.fill_diagonal(model.network[1:], D)
    np.fill_diagonal(model.network[:, 1:], D)
    model.components = [
        SIR.Susceptible(model),
        SIR.Recovered(model),
        SIR.Infectious(model, infdurdist),
        SIR.Transmission(model, infdurdist),
    ]
    model.run(f"Sim {i+1}/{n_samples}")
    # Estimate wavespeed as before
    peak_times = model.nodes.I.argmax(axis=0)
    right_nodes = np.arange(101, npatches)
    left_nodes = np.arange(100, -1, -1)
    right_peak_times = peak_times[right_nodes]
    left_peak_times = peak_times[left_nodes]
    # Right
    mask_r = right_peak_times > 0
    if np.any(mask_r):
        slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
        velocity_right = 1 / slope_r if slope_r != 0 else np.nan
    else:
        velocity_right = np.nan
    # Left
    mask_l = left_peak_times > 0
    if np.any(mask_l):
        slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
        velocity_left = -1 / slope_l if slope_l != 0 else np.nan
    else:
        velocity_left = np.nan
    # Expected wavespeed (from theory)
    expected_speed = 2 * np.sqrt(D / infmean * (R0 - 1))
    results.append({
        "D": D,
        "R0": R0,
        "infmean": infmean,
        "velocity_right": velocity_right,
        "velocity_left": velocity_left,
        "expected_speed": expected_speed,
    })

results_df = pd.DataFrame(results)
# Set up parameter ranges
n_samples = 100
D_samples = np.logspace(np.log10(3e-4), np.log10(3e-1), n_samples)
R0_samples = np.random.uniform(1.5, 10, n_samples)
infmean_samples = np.random.uniform(7, 50, n_samples)

results = []

for i, (D, R0, infmean) in enumerate(zip(D_samples, R0_samples, infmean_samples)):
    # Set up parameters for this run
    parameters = PropertySet(
        {
            "seed": np.random.randint(0, 1000000),
            "nticks": nticks,
            "verbose": False,
            "beta": R0 / infmean,
            "inf_mean": infmean,
        }
    )
    infdurdist = dists.exponential(scale=parameters.inf_mean)
    model = Model(scenario.copy(), parameters)
    model.network = 0 * model.network
    np.fill_diagonal(model.network[1:], D)
    np.fill_diagonal(model.network[:, 1:], D)
    model.components = [
        SIR.Susceptible(model),
        SIR.Recovered(model),
        SIR.Infectious(model, infdurdist),
        SIR.Transmission(model, infdurdist),
    ]
    model.run(f"Sim {i+1}/{n_samples}")
    # Estimate wavespeed as before
    peak_times = model.nodes.I.argmax(axis=0)
    right_nodes = np.arange(101, npatches)
    left_nodes = np.arange(100, -1, -1)
    right_peak_times = peak_times[right_nodes]
    left_peak_times = peak_times[left_nodes]
    # Right
    mask_r = right_peak_times > 0
    if np.any(mask_r):
        slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
        velocity_right = 1 / slope_r if slope_r != 0 else np.nan
    else:
        velocity_right = np.nan
    # Left
    mask_l = left_peak_times > 0
    if np.any(mask_l):
        slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
        velocity_left = -1 / slope_l if slope_l != 0 else np.nan
    else:
        velocity_left = np.nan
    # Expected wavespeed (from theory)
    expected_speed = 2 * np.sqrt(D / infmean * (R0 - 1))
    results.append({
        "D": D,
        "R0": R0,
        "infmean": infmean,
        "velocity_right": velocity_right,
        "velocity_left": velocity_left,
        "expected_speed": expected_speed,
    })

results_df = pd.DataFrame(results)

Sim 1/100: 100%|██████████| 730/730 [00:00<00:00, 1589.45it/s]
Sim 2/100: 100%|██████████| 730/730 [00:00<00:00, 1795.47it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:41: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
Sim 3/100: 100%|██████████| 730/730 [00:00<00:00, 1771.04it/s]
Sim 4/100: 100%|██████████| 730/730 [00:00<00:00, 1731.69it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:41: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
Sim 5/100: 100%|██████████| 730/730 [00:00<00:00, 1782.85it/s]
Sim 6/100: 100%|██████████| 730/730 [00:00<00:00, 1615.69it/s]
Sim 7/100: 100%|██████████| 730/730 [00:00<00:00, 1776.41it/s]
Sim 8/100: 100%|██████████| 730/730 [00:00<00:00, 1766.32it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:41: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
Sim 9/100: 100%|██████████| 730/730 [00:00<00:00, 1773.55it/s]
Sim 10/100: 100%|██████████| 730/730 [00:00<00:00, 1693.85it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:41: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
Sim 11/100: 100%|██████████| 730/730 [00:00<00:00, 1800.68it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:48: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
Sim 12/100: 100%|██████████| 730/730 [00:00<00:00, 1757.24it/s]
Sim 13/100: 100%|██████████| 730/730 [00:00<00:00, 1759.43it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:48: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
Sim 14/100: 100%|██████████| 730/730 [00:00<00:00, 1822.74it/s]
Sim 15/100: 100%|██████████| 730/730 [00:00<00:00, 1821.18it/s]
Sim 16/100: 100%|██████████| 730/730 [00:00<00:00, 1800.06it/s]
Sim 17/100: 100%|██████████| 730/730 [00:00<00:00, 1774.61it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:48: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
Sim 18/100: 100%|██████████| 730/730 [00:00<00:00, 1656.45it/s]
Sim 19/100: 100%|██████████| 730/730 [00:00<00:00, 1760.82it/s]
Sim 20/100: 100%|██████████| 730/730 [00:00<00:00, 1786.95it/s]
Sim 21/100: 100%|██████████| 730/730 [00:00<00:00, 1570.53it/s]
Sim 22/100: 100%|██████████| 730/730 [00:00<00:00, 1749.97it/s]
Sim 23/100: 100%|██████████| 730/730 [00:00<00:00, 1779.38it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:48: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
Sim 24/100: 100%|██████████| 730/730 [00:00<00:00, 1822.08it/s]
Sim 25/100: 100%|██████████| 730/730 [00:00<00:00, 1812.67it/s]
Sim 26/100: 100%|██████████| 730/730 [00:00<00:00, 1564.90it/s]
Sim 27/100: 100%|██████████| 730/730 [00:00<00:00, 1789.10it/s]
Sim 28/100: 100%|██████████| 730/730 [00:00<00:00, 1745.11it/s]
Sim 29/100: 100%|██████████| 730/730 [00:00<00:00, 1771.04it/s]
Sim 30/100: 100%|██████████| 730/730 [00:00<00:00, 1794.92it/s]
Sim 31/100: 100%|██████████| 730/730 [00:00<00:00, 1812.69it/s]
Sim 32/100: 100%|██████████| 730/730 [00:00<00:00, 1807.62it/s]
Sim 33/100: 100%|██████████| 730/730 [00:00<00:00, 1793.64it/s]
Sim 34/100: 100%|██████████| 730/730 [00:00<00:00, 1796.68it/s]
Sim 35/100: 100%|██████████| 730/730 [00:00<00:00, 1781.41it/s]
Sim 36/100: 100%|██████████| 730/730 [00:00<00:00, 1822.00it/s]
Sim 37/100: 100%|██████████| 730/730 [00:00<00:00, 1536.01it/s]
Sim 38/100: 100%|██████████| 730/730 [00:00<00:00, 1744.07it/s]
Sim 39/100: 100%|██████████| 730/730 [00:00<00:00, 1728.67it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:48: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_l, *_ = linregress(left_nodes[mask_l], left_peak_times[mask_l])
Sim 40/100: 100%|██████████| 730/730 [00:00<00:00, 1770.27it/s]
Sim 41/100: 100%|██████████| 730/730 [00:00<00:00, 1822.54it/s]
Sim 42/100: 100%|██████████| 730/730 [00:00<00:00, 1835.71it/s]
Sim 43/100: 100%|██████████| 730/730 [00:00<00:00, 1475.71it/s]
Sim 44/100: 100%|██████████| 730/730 [00:00<00:00, 1697.52it/s]
Sim 45/100: 100%|██████████| 730/730 [00:00<00:00, 1812.07it/s]
Sim 46/100: 100%|██████████| 730/730 [00:00<00:00, 1791.60it/s]
Sim 47/100: 100%|██████████| 730/730 [00:00<00:00, 1788.70it/s]
Sim 48/100: 100%|██████████| 730/730 [00:00<00:00, 1791.85it/s]
Sim 49/100: 100%|██████████| 730/730 [00:00<00:00, 1820.97it/s]
Sim 50/100: 100%|██████████| 730/730 [00:00<00:00, 1784.19it/s]
Sim 51/100: 100%|██████████| 730/730 [00:00<00:00, 1816.44it/s]
Sim 52/100: 100%|██████████| 730/730 [00:00<00:00, 1782.10it/s]
Sim 53/100: 100%|██████████| 730/730 [00:00<00:00, 1783.28it/s]
Sim 54/100: 100%|██████████| 730/730 [00:00<00:00, 1767.16it/s]
Sim 55/100: 100%|██████████| 730/730 [00:00<00:00, 1713.80it/s]
Sim 56/100: 100%|██████████| 730/730 [00:00<00:00, 1463.66it/s]
Sim 57/100: 100%|██████████| 730/730 [00:00<00:00, 1727.47it/s]
Sim 58/100: 100%|██████████| 730/730 [00:00<00:00, 1701.79it/s]
Sim 59/100: 100%|██████████| 730/730 [00:00<00:00, 1769.79it/s]
Sim 60/100: 100%|██████████| 730/730 [00:00<00:00, 1781.24it/s]
Sim 61/100: 100%|██████████| 730/730 [00:00<00:00, 1754.85it/s]
Sim 62/100: 100%|██████████| 730/730 [00:00<00:00, 1754.89it/s]
Sim 63/100: 100%|██████████| 730/730 [00:00<00:00, 1396.86it/s]
Sim 64/100: 100%|██████████| 730/730 [00:00<00:00, 1819.93it/s]
Sim 65/100: 100%|██████████| 730/730 [00:00<00:00, 1716.20it/s]
Sim 66/100: 100%|██████████| 730/730 [00:00<00:00, 1759.22it/s]
Sim 67/100: 100%|██████████| 730/730 [00:00<00:00, 1738.89it/s]
Sim 68/100: 100%|██████████| 730/730 [00:00<00:00, 1787.02it/s]
Sim 69/100: 100%|██████████| 730/730 [00:00<00:00, 1787.82it/s]
Sim 70/100: 100%|██████████| 730/730 [00:00<00:00, 1777.94it/s]
Sim 71/100: 100%|██████████| 730/730 [00:00<00:00, 1750.98it/s]
Sim 72/100: 100%|██████████| 730/730 [00:00<00:00, 1899.78it/s]
Sim 73/100: 100%|██████████| 730/730 [00:00<00:00, 1703.63it/s]
Sim 74/100: 100%|██████████| 730/730 [00:00<00:00, 1741.50it/s]
Sim 75/100: 100%|██████████| 730/730 [00:00<00:00, 1784.48it/s]
Sim 76/100: 100%|██████████| 730/730 [00:00<00:00, 1761.36it/s]
/var/folders/fl/_ns4br_j2qxd92vjrgl3t4n00000gn/T/ipykernel_96092/1872132150.py:41: SmallSampleWarning: One or more sample arguments is too small; all returned values will be NaN. See documentation for sample size requirements.
  slope_r, *_ = linregress(right_nodes[mask_r], right_peak_times[mask_r])
Sim 77/100: 100%|██████████| 730/730 [00:00<00:00, 1782.17it/s]
Sim 78/100: 100%|██████████| 730/730 [00:00<00:00, 1777.98it/s]
Sim 79/100: 100%|██████████| 730/730 [00:00<00:00, 1755.14it/s]
Sim 80/100: 100%|██████████| 730/730 [00:00<00:00, 1716.72it/s]
Sim 81/100: 100%|██████████| 730/730 [00:00<00:00, 1742.29it/s]
Sim 82/100: 100%|██████████| 730/730 [00:00<00:00, 1763.68it/s]
Sim 83/100: 100%|██████████| 730/730 [00:00<00:00, 1750.55it/s]
Sim 84/100: 100%|██████████| 730/730 [00:00<00:00, 1479.91it/s]
Sim 85/100: 100%|██████████| 730/730 [00:00<00:00, 1773.55it/s]
Sim 86/100: 100%|██████████| 730/730 [00:00<00:00, 1758.66it/s]
Sim 87/100: 100%|██████████| 730/730 [00:00<00:00, 1770.89it/s]
Sim 88/100: 100%|██████████| 730/730 [00:00<00:00, 1924.05it/s]
Sim 89/100: 100%|██████████| 730/730 [00:00<00:00, 1745.05it/s]
Sim 90/100: 100%|██████████| 730/730 [00:00<00:00, 1754.55it/s]
Sim 91/100: 100%|██████████| 730/730 [00:00<00:00, 1756.36it/s]
Sim 92/100: 100%|██████████| 730/730 [00:00<00:00, 1730.06it/s]
Sim 93/100: 100%|██████████| 730/730 [00:00<00:00, 1761.04it/s]
Sim 94/100: 100%|██████████| 730/730 [00:00<00:00, 1821.00it/s]
Sim 95/100: 100%|██████████| 730/730 [00:00<00:00, 1725.48it/s]
Sim 96/100: 100%|██████████| 730/730 [00:00<00:00, 1860.89it/s]
Sim 97/100: 100%|██████████| 730/730 [00:00<00:00, 1877.07it/s]
Sim 98/100: 100%|██████████| 730/730 [00:00<00:00, 1737.56it/s]
Sim 99/100: 100%|██████████| 730/730 [00:00<00:00, 1377.37it/s]
Sim 100/100: 100%|██████████| 730/730 [00:00<00:00, 1863.82it/s]

In [7]:

Copied!





plt.figure(figsize=(8, 6))
plt.plot(results_df["expected_speed"], results_df["velocity_left"], 'r.', label="Velocity left", markersize=12)
plt.plot(results_df["expected_speed"], results_df["velocity_right"], 'b.', label="Velocity right", markersize=12)
plt.plot([1e-3, 1], [1e-3, 1], 'k--', label="y=x")
plt.xscale('log')
plt.yscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Measured wavespeed (patches/tick)")
plt.title("Comparison of expected and observed wavespeeds")
min_val = np.nanmin(np.abs([results_df["expected_speed"].min(), results_df["velocity_right"].min(), results_df["velocity_left"].min()]))
max_val = np.nanmax([results_df["expected_speed"].max(), results_df["velocity_right"].max(), results_df["velocity_left"].max()])
plt.xlim(min_val, max_val)
plt.ylim(min_val, max_val)
plt.legend()
plt.tight_layout()
plt.show()
plt.figure(figsize=(8, 6))
plt.plot(results_df["expected_speed"], results_df["velocity_left"], 'r.', label="Velocity left", markersize=12)
plt.plot(results_df["expected_speed"], results_df["velocity_right"], 'b.', label="Velocity right", markersize=12)
plt.plot([1e-3, 1], [1e-3, 1], 'k--', label="y=x")
plt.xscale('log')
plt.yscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Measured wavespeed (patches/tick)")
plt.title("Comparison of expected and observed wavespeeds")
min_val = np.nanmin(np.abs([results_df["expected_speed"].min(), results_df["velocity_right"].min(), results_df["velocity_left"].min()]))
max_val = np.nanmax([results_df["expected_speed"].max(), results_df["velocity_right"].max(), results_df["velocity_left"].max()])
plt.xlim(min_val, max_val)
plt.ylim(min_val, max_val)
plt.legend()
plt.tight_layout()
plt.show()

Results¶

As expected, at high speeds the observed wavespeed is consistently slower than theoretical expetactation. At slower wavespeeds, the model is more able to reproduce waves traveling around the correct speed, but there is still a bias towards the model being slower than expected. Below, we will take a look at how the errors behave with model parameters.

In [8]:

Copied!





# Compute relative errors
rel_err_left = (results_df["velocity_left"] - results_df["expected_speed"]) / results_df["expected_speed"]
rel_err_right = (results_df["velocity_right"] - results_df["expected_speed"]) / results_df["expected_speed"]

# Plot: colored by R0
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=results_df["R0"]/results_df['infmean'], cmap='viridis', label="Right", marker='o')
plt.colorbar(sc, label="beta")
plt.scatter(results_df["expected_speed"], rel_err_left, c=results_df["R0"]/results_df['infmean'], cmap='viridis', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by beta)")
plt.legend()
plt.tight_layout()
plt.show()

# Plot: colored by D
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=np.log10(results_df["D"]), cmap='plasma', label="Right", marker='o')
plt.colorbar(sc, label="log10(D)")
plt.scatter(results_df["expected_speed"], rel_err_left, c=np.log10(results_df["D"]), cmap='plasma', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by log10(D))")
plt.legend()
plt.tight_layout()
plt.show()

# Plot: colored by infmean
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=results_df["infmean"], cmap='cividis', label="Right", marker='o')
plt.colorbar(sc, label="inf_mean")
plt.scatter(results_df["expected_speed"], rel_err_left, c=results_df["infmean"], cmap='cividis', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by inf_mean)")
plt.legend()
plt.tight_layout()
plt.show()
# Compute relative errors
rel_err_left = (results_df["velocity_left"] - results_df["expected_speed"]) / results_df["expected_speed"]
rel_err_right = (results_df["velocity_right"] - results_df["expected_speed"]) / results_df["expected_speed"]

# Plot: colored by R0
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=results_df["R0"]/results_df['infmean'], cmap='viridis', label="Right", marker='o')
plt.colorbar(sc, label="beta")
plt.scatter(results_df["expected_speed"], rel_err_left, c=results_df["R0"]/results_df['infmean'], cmap='viridis', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by beta)")
plt.legend()
plt.tight_layout()
plt.show()

# Plot: colored by D
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=np.log10(results_df["D"]), cmap='plasma', label="Right", marker='o')
plt.colorbar(sc, label="log10(D)")
plt.scatter(results_df["expected_speed"], rel_err_left, c=np.log10(results_df["D"]), cmap='plasma', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by log10(D))")
plt.legend()
plt.tight_layout()
plt.show()

# Plot: colored by infmean
plt.figure(figsize=(8, 6))
sc = plt.scatter(results_df["expected_speed"], rel_err_right, c=results_df["infmean"], cmap='cividis', label="Right", marker='o')
plt.colorbar(sc, label="inf_mean")
plt.scatter(results_df["expected_speed"], rel_err_left, c=results_df["infmean"], cmap='cividis', marker='x', label="Left")
plt.xscale('log')
plt.xlabel("Expected wavespeed (patches/tick)")
plt.ylabel("Relative error (measured - expected) / expected")
plt.title("Relative error vs. expected speed (colored by inf_mean)")
plt.legend()
plt.tight_layout()
plt.show()

Discussion¶

We can see clear vertical gradients in the error plots colored by $\beta$ and infectious duration - lower values of $beta$, and longer infectious durations are both clearly correlated with producing observed waves that are slower than theoretical expectations. This makes some sense - lower values of both parameters translate into lower per-infective, per-timestep force of infection, which can cause delays in both the important processes of spatial spread of infections and the growth of local outbreaks in each node. The dependence on D is present but less obvious, as this parameter can only affect the delays related to spatial spread.