Case Study: Susceptible-Infectious Model - Uncertainty

We model a simple susceptible-infectious epidemic in which susceptible persons can get infected from infectious person and thereafter become infectious themselves. We use the model to investigate stochasticity and uncertainty.

Topics:

Uncertainty Quantification
Linkage Between Microscopic and Macroscopic Models

Objectives:

How do agent-based and differential equation models compare?
What is the impact of locality in agent-based models?
How do we quantify stochasticity in stochastic models?

import os
import random

import holoviews as hv
import numpy as np

# suppress warnings
if "LD_PRELOAD" in os.environ:
    del os.environ["LD_PRELOAD"]

from itertools import product
from typing import Self

from joblib import Parallel, delayed
from scipy.stats import norm

from modelling_and_simulation_cookbook import Clock

# default settings for the plots
hv.extension("bokeh")
hv.opts.defaults(
    hv.opts.Curve(
        responsive=True,
        height=400,
        xlabel="timesteps",
        ylabel="number of infected agents",
        tools=["hover"],
        active_tools=["pan"],
        backend_opts={"plot.toolbar.autohide": True},
    )
)
hv.output(widget_location="bottom")

We are considering a disease that is transmitted through contact and renders an individual permanently infectious once infected. This simplified picture provides us with a so-called SI (susceptible-infectious) model and can be implemented at both the microscopic (e.g. agent-based) and macroscopic (e.g. system-dynamics / ordinary differential equation) scope. At first, we want to implement the model as an agent-based one.

Conceptual Modelling - Agent-Based Approach¶

The $N$ agents of the model represent persons and every agent $i=1\dots N$ has one single binary state infected defining whether it is susceptible ( $a_i(t)=S$ ) or infectious ( $a_i(t)=I$ ) at time / model step $t$ . Furthermore, the agent has a fixed list of contacts (i.e. other agents) with whom he is in regular contact causing the spread of the pathogen. We want to leave the number and sampling strategy of the contacts open for the parametrisation. However, a fraction of $x_0$ randomly drawn agents will be initialised infectious whereas the others are created susceptible.

The model is updated in discrete time steps of length $\Delta t$ in which every agent selects one randomly drawn neighbor and has a contact with it. If the agent is infectious, there is probability $\alpha$ that the contact partner will be rendered infectious after the time step. The time-dynamics of the total number of infectious agents is the observable of the model.

class Agent:
    """
    Agent class for the agent-based SI simulation
    :param infected: whether the agent is infected on creation
    """

    def __init__(self, infected: bool) -> None:
        """Initialises an agent"""
        self.infected = infected
        self.neighbours: list[Self] = list()

    def infect(self) -> None:
        """Renders the agent infected"""
        self.infected = True

    def is_infected(self) -> bool:
        """:returns: True if the agent is infected"""
        return self.infected

    def is_susceptible(self) -> bool:
        """:returns: True if the agent is susceptible"""
        return not self.infected

    def get_neighbours(self) -> list[Self]:
        """:returns: list of neighbour agents"""
        return self.neighbours

    def add_neighbor(self, other: Self) -> None:
        """
        Adds another agent as neighbour
        :param other: other agent
        """
        self.neighbours.append(other)


class SISimulationABM:
    """
    Class that handles the agent-based SI Simulation.
    :param N: number of agents
    :param alpha: infection probability
    :param x0: initial fraction of infected agents
    :param neighbours: number of neighbours per agent
    :param seed: Random number seed. If None, a random seed is used
    """

    def __init__(
        self,
        N: int,
        alpha: float,
        x0: float,
        neighbours: int,
        seed: int = None,
        **kwargs,
    ):
        """Initializes an agent-based SI simulation"""
        self.rng = random.Random(seed)
        self.agents = list()
        self.alpha = alpha
        for i in range(N):
            a = Agent(True) if i < (N * x0) else Agent(False)
            self.agents.append(a)
        for i in range(N):
            for _ in range(neighbours):
                i2 = self.rng.randint(0, N - 1)
                self.agents[i].add_neighbor(self.agents[i2])

    def run(self, tend: int, **kwargs) -> np.ndarray:
        """
        Runs the agent-based model for "tend" time-steps
        :param tend: number of time-steps
        :returns: list with total number of infected agents per step
        """
        X = [sum([1 for x in self.agents if x.is_infected()])]
        for _ in range(tend):
            infectees = list()  # make sure to update all agents simultaneously
            for a in self.agents:
                if a.is_infected():
                    nbs = a.get_neighbours()
                    nb = self.rng.choice(nbs)
                    if nb.is_susceptible() and self.rng.random() < self.alpha:
                        infectees.append(nb)
            for a in infectees:
                a.infect()
            X.append(sum([1 for x in self.agents if x.is_infected()]))
        return np.array(X)


X = SISimulationABM(N=1000, alpha=0.1, x0=0.05, neighbours=50, seed=12345).run(100)


def create_frame(X, t):
    """creates a single frame for the visualization"""
    current_data = X[: t + 1]
    time = list(range(t + 1))
    frame = hv.Curve((time, current_data)).opts(
        title=f"Single Simulation Run (t = {t})",
    )
    return frame


frames = {t: create_frame(X, t) for t in range(101)}
animated_curve = hv.HoloMap(frames, kdims=["Time"])
hv.output(animated_curve, holomap="scrubber")  # remove holomap option to get a slider

Monte-Carlo Simulation¶

Clearly one simulation run is not representative. In order to get an image of the underlying distribution we need to apply Monte-Carlo (MC) simulation and run the simulation several times.

params = {"N": 2_000, "alpha": 0.1, "x0": 0.05, "neighbours": 50}
mc_iters = 16


def run_mc_single(
    model_params: dict, iters: int = None, tend: int = 100, root_seed: int = None
) -> list[np.ndarray]:
    """
    Performs single threaded Monte Carlo simulation
    :param model_params: model parameters
    :param iters: number of Monte Carlo iterations
    :param tend: number of timesteps for each simulation
    :param root_seed: seed from which all further seeds are created
    :returns: list of all simulation results
    """
    xs = list()
    seeds = np.random.SeedSequence(root_seed).spawn(iters)
    for s in seeds:
        seed = int(s.generate_state(1)[0])
        X = SISimulationABM(**model_params, seed=seed).run(tend)
        xs.append(X)
    return xs


clock = Clock()
clock.start()
xs = run_mc_single(model_params=params, iters=mc_iters, tend=100, root_seed=12345)
st_time = clock.stop(f"computation time for {mc_iters} simulations (single threaded)")
st_time = st_time.total_seconds()

line_dict = dict()
for i, x in enumerate(xs):
    line_dict[f"simulation {i}"] = hv.Curve(x)

plot = hv.NdOverlay(line_dict)
plot.opts(show_legend=False, title="Simulation Runs for different Seeds")

computation time for 16 simulations (single threaded): 2.777491 seconds

Since MC runs are (need to be) independent, we may use multi-threading/processing to do this efficiently. In the way this is done below, this only improves runtime if the computation time of one single simulation is long enough. Otherwise, multithreading creates too much overhead.

Note that we used joblib.Parallel and joblib.delayed because of the memory structure of Jupyter notebooks. Also note, that the results are only equivalent if the seeds are assigned to the simulation instances by the main file (as done here).

def _mc_worker_task(params: dict) -> np.ndarray:
    tend = params.pop("tend")
    return SISimulationABM(**params).run(tend)


def run_mc_multi(
    model_params: dict | list[dict],
    iters: int = 8,
    tend: int = 100,
    root_seed: int = None,
    n_threads: int = 8,
) -> list[np.ndarray]:
    """
    Performs multithreaded Monte Carlo simulation
    :param model_params: model parameters
    :param iters: number of Monte Carlo iterations
    :param tend: number of timesteps for each simulation
    :param root_seed: seed from which all further seeds are created
    :param n_threads: maximum number of threads to use
    :returns: list of all simulation results
    """
    assert isinstance(model_params, dict) or len(model_params) <= iters, (
        f"List of parameters must have a length of at least {iters}!"
    )

    xs, params_list = list(), list()
    seeds = np.random.SeedSequence(root_seed).spawn(iters)
    for i, s in enumerate(seeds):
        if isinstance(model_params, list):
            params = model_params[i]
        else:
            params = model_params.copy()
        params["seed"] = int(s.generate_state(1)[0])
        params["tend"] = tend
        params_list.append(params)
    xs = Parallel(n_jobs=min(len(seeds), n_threads))(
        delayed(_mc_worker_task)(params) for params in params_list
    )
    return xs


clock.start()
xs = run_mc_multi(model_params=params, iters=mc_iters, tend=100, root_seed=12345)
mt_time = clock.stop(f"computation time for {mc_iters} simulations (multithreaded)")
mt_time = mt_time.total_seconds()
print(f"  \u2192 multithreading {st_time / mt_time:.2f} times faster!")

line_dict = dict()
for i, x in enumerate(xs):
    line_dict[f"simulation {i}"] = hv.Curve(x)

plot = hv.NdOverlay(line_dict)
plot.opts(show_legend=False, title="Simulation Runs for different Seeds")

computation time for 16 simulations (multithreaded): 3.021276 seconds
  → multithreading 0.92 times faster!

Uncertainty¶

The model is ineretly stochastic, meaning that it incorporates a native way to depict uncertainty. However, in case parameters or inputs are also uncertain, it is possible to add additional randomess.

Below we investigate the impact of additional uncertainty w.r. to the initial value and the parameter alpha in form of a triangular distributed random factor $X\sim Tri(0.7,1.0,1.3)$ .

params = {"N": 100, "alpha": 0.1, "x0": 0.05, "neighbours": 50}
mc_iters = 16

params_list = list()
for _ in range(mc_iters):
    params_temp = params.copy()
    params_temp["x0"] = params["x0"] * random.triangular(0.7, 1.0, 1.3)
    params_list.append(params_temp)
xs = run_mc_multi(model_params=params_list, iters=mc_iters, tend=100, root_seed=12345)

line_dict = dict()
for i, x in enumerate(xs):
    line_dict[f"(run {i}) x0 = {params_list[i]['x0']:.3f}"] = hv.Curve(x)

plot = hv.NdOverlay(line_dict)
plot.opts(show_legend=False, title="Additional Input-Uncertainty in x0")

params = {"N": 100, "alpha": 0.1, "x0": 0.05, "neighbours": 50}
mc_iters = 16

params_list = list()
for _ in range(mc_iters):
    params_temp = params.copy()
    params_temp["alpha"] = params["alpha"] * random.triangular(0.7, 1.0, 1.3)
    params_list.append(params_temp)
xs = run_mc_multi(model_params=params_list, iters=mc_iters, tend=100, root_seed=12345)

line_dict = dict()
for i, x in enumerate(xs):
    line_dict[f"(run {i}) \u03b1 = {params_list[i]['alpha']:.3f}"] = hv.Curve(x)

plot = hv.NdOverlay(line_dict)
plot.opts(show_legend=False, title="Additional Parameter-Uncertainty in \u03b1")

Prediction Intervals¶

In the following code cell there are two different approaches to calculate the prediction intervals of future simulation runs.

First, we look at the gaussian approach. For example, roughly 68% of the simulation results should fall in the interval $[\mu-\sigma, \mu+\sigma]$ for every time step. By performing Monte-Carlo simulations, we can calculate the mean $\hat\mu$ and the standard deviation $\hat\sigma$ of the simulation results. However, this approach assumes that our simulation results are normally distributed. Without further statistical tests, this assumption cannot be made.

If we have no knowledge of the underlying distribution, we can still make predictions by looking at quantiles of our simulation results. Assuming we performed enough Monte-Carlo simulations, we can presume that the distribution of the results is representative of the real distribution. Therefore, if 95% of our simulation results lie in a specific interval, it is probable that this will be the case for future simulation runs as well.

Keep in mind that for both of these approaches the reliability of the prediction intervals depends on the number of Monte-Carlo simulations. The more simulations we perform, the more reliable the prediction intervals will be.

params = {"N": 100, "alpha": 0.1, "x0": 0.05, "neighbours": 50, "tend": 100}
mc_iters = 100
seeds = np.random.SeedSequence(1234).spawn(mc_iters)

X = []  # list where we store our individual simulation results
Xm = np.zeros(shape=(params["tend"] + 1,))  # mean of the results for every timestep
Xm2 = np.zeros(
    shape=(params["tend"] + 1,)
)  # mean of the squared results for every timestep

for i in range(mc_iters):
    params["seed"] = int(seeds[i].generate_state(1)[0])
    X.append(SISimulationABM(**params).run(**params))
    Xm = (i * Xm + X[-1]) / (i + 1)
    Xm2 = (i * Xm2 + X[-1] ** 2) / (i + 1)
    if i > 1:
        Sm2 = (
            i / (i - 1) * (Xm2 - Xm**2)
        )  # i / (i - 1) corrector to become free of bias
Sm = np.array([x**0.5 for x in Sm2])

# plotting
t = np.arange(params["tend"] + 1)
sim_data = [np.column_stack((t, x)) for x in X]
sim_paths = hv.Path(sim_data, label="Simulation Runs").opts(
    color="deepskyblue", alpha=0.3, show_legend=True
)
mean_curve = hv.Curve((t, Xm), label="Mean").opts(color="crimson", line_width=2)
pred_68 = hv.Area(
    (t, Xm - 1.00 * Sm, Xm + 1.00 * Sm), vdims=["y1", "y2"], label="68% PI"
).opts(color="crimson", fill_alpha=0.4, line_alpha=0, show_legend=True)
pred_95 = hv.Area(
    (t, Xm - 1.96 * Sm, Xm + 1.96 * Sm), vdims=["y1", "y2"], label="95% PI"
).opts(color="crimson", fill_alpha=0.2, line_alpha=0, show_legend=True)

plot_gauss = sim_paths * pred_95 * pred_68 * mean_curve
plot_gauss.opts(legend_position="bottom_right", title="Gaussian predictions intervals")

pc1, pc2, pc3, pc4 = np.percentile(np.array(X), [2.5, 16, 84, 97.5], axis=0)
p3 = hv.Area((t, pc2, pc3), vdims=["y1", "y2"], label="68% PI").opts(
    color="crimson", fill_alpha=0.4, line_alpha=0, show_legend=True
)
p4 = hv.Area((t, pc1, pc4), vdims=["y1", "y2"], label="95% PI").opts(
    color="crimson", fill_alpha=0.2, line_alpha=0, show_legend=True
)

plot_empiric = sim_paths * p4 * p3 * mean_curve
plot_empiric.opts(
    legend_position="bottom_right", title="Empirical (percentile) predictions intervals"
)

layout = hv.Layout([plot_gauss, plot_empiric]).cols(2)
layout

Early Stopping¶

Depending on the complexity of the model, Monte-Carlo simulation can be very computationally expensive. It is therefore useful to stop the simulation as soon as the results are stable enough. For further details on the theory underlying this section, see Monte Carlo Simulation.

In the cell below, we start out with a warm-up period where 6 simulations runs are performed regardless of the stopping rule. Afterwards we calculate the necessary statistics after every simulation run and check if the stopping rule is satisfied. If it is, the Monte-Carlo simulation is stopped.

params = {"alpha": 0.1, "x0": 0.05, "neighbours": 50, "tend": 100}
N_list = [1_000, 100, 50]
p_list = [0.95, 0.9]
M_limit = 10_000
combinations = list(product(N_list, p_list))
seeds = np.random.SeedSequence(1234).spawn(len(combinations))


def stopping_rule(sigma2, M, p, delta):
    return (1 - p) - 2 * norm.cdf(-((M / sigma2[1:]) ** 0.5) * delta)


plots, t = list(), np.arange(params["tend"] + 1)

for (N, p), seed in zip(combinations, seeds, strict=True):
    delta = 0.01 * N  # absolute difference allowed
    params["N"] = N
    sample_seeds = iter(seed.spawn(M_limit))

    clock.start()
    X = []  # list where we store our individual simulation results
    Xm = np.zeros(shape=(101,))  # mean of the results for every timestep
    Xm2 = np.zeros(shape=(101,))  # mean of the squared results for every timestep

    M = 0
    while M <= 5:
        params["seed"] = int(next(sample_seeds).generate_state(1)[0])
        X.append(SISimulationABM(**params).run(**params))
        Xm = (M * Xm + X[-1]) / (M + 1)
        Xm2 = (M * Xm2 + X[-1] ** 2) / (M + 1)
        M += 1

    Sm2 = np.maximum(M / (M - 1) * (Xm2 - Xm**2), 1e-10)
    out = stopping_rule(Sm2, M, p, delta)

    while not all(out >= 0) and M_limit > M:
        params["seed"] = int(next(sample_seeds).generate_state(1)[0])
        X.append(SISimulationABM(**params).run(**params))
        Xm = (M * Xm + X[-1]) / (M + 1)
        Xm2 = (M * Xm2 + X[-1] ** 2) / (M + 1)
        M += 1
        Sm2 = np.maximum(M / (M - 1) * (Xm2 - Xm**2), 1e-10)  # avoid division by zero
        out = stopping_rule(Sm2, M, p, delta)

    if M_limit <= M:
        print(
            f"Monte Carlo simulation stopped after {M} iterations "
            f"- stopping rule not achieved."
        )

    clock.stop(f"Elapsed time for N = {N:4d} and p = {p:.2f} ({M:4d} runs)")

    # create the plots
    path_data = [np.column_stack((t, x)) for x in X]
    sims = hv.Path(path_data, label="Simulation Runs").opts(
        color="deepskyblue", alpha=0.2
    )
    mean = hv.Curve((t, Xm), label="Mean").opts(color="crimson", line_width=2)
    plot = (sims * mean).opts(
        title=f"{M} Runs with N = {N} and p = {p}", legend_position="bottom_right"
    )
    plots.append(plot)

layout = hv.Layout(plots).cols(max(2, int(len(plots) / 3)))
layout.opts(shared_axes=False)
layout

Elapsed time for N = 1000 and p = 0.95 (  63 runs): 5.723105 seconds

Elapsed time for N = 1000 and p = 0.90 (  36 runs): 2.882337 seconds

Elapsed time for N =  100 and p = 0.95 ( 521 runs): 4.237106 seconds

Elapsed time for N =  100 and p = 0.90 ( 351 runs): 2.753829 seconds

Elapsed time for N =   50 and p = 0.95 ( 949 runs): 3.713852 seconds

Elapsed time for N =   50 and p = 0.90 ( 640 runs): 2.036707 seconds

Influence of Agent Count¶

Before we have a look at the influence varying agent-counts have on the simulations results we establish the mean-field model for the ABM approach using ordinary differential equations:

Let $I(t)=|\{i:a_i(t)=I\}|$ denote the total number of susceptible agents and $S(t)=|\{i:a_i(t)=S\}|=N-I(t)$ the total number of infectious agents at time $t$ , then we may apply mean-field theory to derive the differential-equation-model pendent of the agent-based model by looking into the corresponding state-change probabilities of an agent per time. Let $Nbh(i)$ refer to all agents in the beighborhood of agent $i$ , then we find from the agent-based model that

P(S\rightarrow I):=P(a_i(t+1)=I|a_i(t)=S)=\alpha\cdot P(\text{randomly drawn neighbour is infected})=\alpha\cdot P\left(\frac{|\{j:a_j(t)=True,j\in Nbh(i)|}{|\{j:j\in Nbh(i)\}|}\right)

(1)

Under the assumption of homogeneous mixing, the likelihood of a neighboured agent to be infected can be approximated by the infection likelihood of an arbitrary agent in the model, which allows to reduce the probability calculation to the macroscopic states $S,I$ :

\alpha\cdot P\left(\frac{|\{j:a_j(t)=I,j\in Nbh(i)|}{|\{j:j\in Nbh(i)\}|}\right)\approx \alpha\cdot P\left(\frac{|\{j:a_j(t)=I\}|}{|\{j\}|}\right)=\alpha \frac{I(t)}{N}.

(2)

Since the likelihood for the inverse state change $P(I\rightarrow S):=P(a_i(t+1)=S|a_i(t)=I)$ is zero in the agent-based SI model, we get the mean field equations

\frac{dS}{dt}(t)=I(t)P(I\rightarrow S)-S(t)P(S\rightarrow I)=-\alpha S(t)\frac{I(t)}{N}

(3)

\frac{dI}{dt}(t)=S(t)P(S\rightarrow I)-I(t)P(I\rightarrow S)=\alpha S(t)\frac{I(t)}{N}

(4)

which can be collapsed into a single equation with $S(t)=N-I(t)$ to

\frac{dI}{dt}(t)=\alpha (N-I(t))\frac{I(t)}{N}.

(5)

params = {"N": 50, "alpha": 0.1, "x0": 0.05, "neighbours": 50, "tend": 100}


def si_simulation_ode(N: int, alpha: float, x0: float, tend: int = 100, **kwargs):
    x, t = x0 * N, 0
    X, T = [x], [t]

    def rhs(t, x, alpha):
        return alpha * x * (N - x) / N

    h = 0.01
    while t < tend:
        x = x + h * rhs(t, x, alpha)  # solution with explicit euler method
        t += h
        X.append(x)
        T.append(t)
    return T, X


T, X_ode = si_simulation_ode(**params)
hv.Curve((T, X_ode)).opts(show_grid=True, title="SI Simulation results with ODE")

In the following cell we keep the product of MC iterations and agent-count constant and find that increasing the agent-count and increasing the MC iterations have roughly the same effect on how well the mean from the ABM matches the mean-field model results.

params = {"alpha": 0.1, "x0": 0.05, "neighbours": 50, "tend": 100}
N_list = [50, 100, 1_000, 5_000]
seeds = np.random.SeedSequence(1234).spawn(len(N_list))
plots = list()

t = np.arange(params["tend"] + 1)
for N, seed in zip(N_list, seeds, strict=True):
    params["N"] = N
    runs = round(
        20_000 / N
    )  # product of MC iterations and agent-count should be 20_000
    sample_seeds = iter(seed.spawn(runs))

    [T_ode, X_ode] = si_simulation_ode(**params)
    X_abm = list()

    clock.start()
    for _ in range(runs):
        params["seed"] = int(next(sample_seeds).generate_state(1)[0])
        X_abm.append(SISimulationABM(**params).run(**params))
    X_mean = sum(X_abm) / len(X_abm)
    clock.stop(f"Elapsed time for {runs:3d} simulations runs and N = {N:4d} agents")

    path_data = [np.column_stack((t, x)) for x in X_abm]
    plot_abm = hv.Path(path_data, label="ABM results").opts(
        color="deepskyblue", alpha=0.3, show_legend=True
    )
    plot_mean = hv.Curve((t, X_mean), label="ABM Mean").opts(color="blue", line_width=2)
    plot_ode = hv.Curve((T_ode, X_ode), label="ODE result").opts(color="crimson")
    plot = (plot_abm * plot_mean * plot_ode).opts(
        title=f"{runs} Runs with N = {N}", legend_position="bottom_right"
    )
    plots.append(plot)

layout = hv.Layout(plots).cols(max(2, int(len(plots) / 3)))
layout.opts(shared_axes=False)
layout

Elapsed time for 400 simulations runs and N =   50 agents: 1.834374 seconds

Elapsed time for 200 simulations runs and N =  100 agents: 1.25758 seconds

Elapsed time for  20 simulations runs and N = 1000 agents: 1.409306 seconds

Elapsed time for   4 simulations runs and N = 5000 agents: 1.461782 seconds

Influence of Neighbourhood¶

This time we fix the agent-count to 500 and perform 20 MC iterations. We only vary the neighbourhood from 500 down to 1. Looking at the influence of different neighbourhood we find that the more local effects the ABM has, the more the result diverges from the mean-field model.

params = {"N": 500, "alpha": 0.1, "x0": 0.05, "tend": 100}
neigbours_list = [params["N"], 50, 5, 1]
[T_ode, X_ode] = si_simulation_ode(**params)
seeds = np.random.SeedSequence(1234).spawn(len(neigbours_list))
runs = 20
plots = list()

t = np.arange(params["tend"] + 1)
for neighbours, seed in zip(neigbours_list, seeds, strict=True):
    sample_seeds = iter(seed.spawn(runs))
    params["neighbours"] = neighbours
    X_abm = list()

    clock.start()
    for _ in range(runs):
        params["seed"] = int(next(sample_seeds).generate_state(1)[0])
        X_abm.append(SISimulationABM(**params).run(**params))
    X_mean = sum(X_abm) / len(X_abm)
    clock.stop(
        f"Elapsed time for {runs} simulations runs and {neighbours:3d} neighbour(s)"
    )

    path_data = [np.column_stack((t, x)) for x in X_abm]
    plot_abm = hv.Path(path_data, label="ABM results").opts(
        color="deepskyblue", alpha=0.5, show_legend=True
    )
    plot_mean = hv.Curve((t, X_mean), label="ABM Mean").opts(color="blue", line_width=2)
    plot_ode = hv.Curve((T_ode, X_ode), label="ODE result").opts(color="crimson")
    plot = (plot_abm * plot_mean * plot_ode).opts(
        title=f"{runs} Runs with {neighbours} neighbour(s)",
        legend_position="bottom_right",
    )
    plots.append(plot)

layout = hv.Layout(plots).cols(max(2, int(len(plots) / 3)))
layout.opts(shared_axes=False)
layout

Elapsed time for 20 simulations runs and 500 neighbour(s): 3.138742 seconds

Elapsed time for 20 simulations runs and  50 neighbour(s): 0.665406 seconds

Elapsed time for 20 simulations runs and   5 neighbour(s): 0.423097 seconds
Elapsed time for 20 simulations runs and   1 neighbour(s): 0.190452 seconds

This development shows the importance of considering local effects and how relevant microscopic models are for this. The provided mean-field model can only depict a situation in which the population is homogeneously mixed and everyone can infect everyone else. While it is not (with reasonable effort) possible to refine the model to any local neighborhood smaller than $N$ , there is at least a way to construct a (better) mean-field model for the special case with neighborhood 1.

Mean-Field Model for Neighbourhood Size 1¶

We apply the same logic as for the derivation of the other mean-field model and investigate the likelihood for state changes:

P(S\rightarrow I)=P(a_i(t+1)=I|a_i(t)=S)=\alpha\cdot P(\text{neighbour is infected}).

(6)

The key difference to the first mean-field model is, that the probability to the right does not scale with $I$ but actually remains constant with value $x_0$ . The intuition behind this is, that every new infectee can, at most, infect one other agent because the one-dimensional structure of the disease network. The mean-field model

\frac{dI}{dt}(t)=\alpha x_0(N-I(t))

(7)

results. We refer to the final section of this notebook for a formal derivation, or, alternatively to Bicher (2017) Chapter 5.3.3, where the mean-field equation is derived from an inifinite coupled system of differential equations.

As seen below, the model only matches the initial upswing of the curve. The reason for this is, that the contact network which is spanned by the agent-neighbour relations, has a high likelihood of not being fully connected or being a tree-structure rather than a linear-loop. This creates unreachable/uninfectable agents and earlier local saturation effects.

def si_simulation_1d_ode(N: int, alpha: float, x0: float, tend: int = 100, **kwargs):
    x, t = x0 * N, 0
    X, T = [x], [t]

    def rhs(t, x, alpha):
        return alpha * x0 * (N - x)

    h = 0.01
    while t < tend:
        x = x + h * rhs(t, x, alpha)  # solution with explicit euler method
        t += h
        X.append(x)
        T.append(t)
    return T, X


params = {"N": 500, "alpha": 0.1, "x0": 0.05, "tend": 100}
neigbours_list = [1]
[T_ode, X_ode] = si_simulation_ode(**params)
[T_ode_1d, X_ode_1d] = si_simulation_1d_ode(**params)
runs = 20
plots = list()

t = np.arange(params["tend"] + 1)
neighbours = 1
seed = np.random.SeedSequence(1234)
sample_seeds = iter(seed.spawn(runs))
params["neighbours"] = neighbours
X_abm = list()

clock.start()
for _ in range(runs):
    params["seed"] = int(next(sample_seeds).generate_state(1)[0])
    X_abm.append(SISimulationABM(**params).run(**params))
X_mean = sum(X_abm) / len(X_abm)
clock.stop(f"Elapsed time for {runs} simulations runs and {neighbours:3d} neighbour(s)")

path_data = [np.column_stack((t, x)) for x in X_abm]
plot_abm = hv.Path(path_data, label="ABM results").opts(
    color="deepskyblue", alpha=0.5, show_legend=True
)
plot_mean = hv.Curve((t, X_mean), label="ABM Mean").opts(color="blue", line_width=2)
plot_ode = hv.Curve((T_ode, X_ode), label="ODE result").opts(color="crimson")
plot_ode_1d = hv.Curve((T_ode_1d, X_ode_1d), label="ODE 1D result").opts(
    color="crimson", line_dash="dashed"
)
plot = (plot_abm * plot_mean * plot_ode * plot_ode_1d).opts(
    title=f"{runs} Runs with {neighbours} neighbour(s)",
    legend_position="bottom_right",
)
display(plot)

Elapsed time for 20 simulations runs and   1 neighbour(s): 0.192434 seconds

Formal derivation of the Mean-Field Model for Neighbourhood Size 1¶

Key for finding the mean-field model is to quantify $P(\text{neighbour is infected})$ . One way to evaluate this formally is to look at (agent,neighbour) tuples and their states: An agent can only get infectd if the corresponding tuple is in an $(S,I)$ state, if the tuple has state $(S,S)$ no infection is possible. Let furthermore $SI(t),SS(t),IS(t),II(t)$ stand for the total number of (agent,neighbour) tuples with states $(S,I),(S,S),(I,S),(I,I)$ in the simulation at time $t$ , then

P(\text{neighbour is infected})=\frac{SI(t)}{S(t)}.

(8)

We now look at the expected evolution of this ration, when an infection occurs in the model. Clearly, quantity $S(t)$ is reduced by one, however, there are two possibilities on what happens with the total number of $(S,I)$ tuples:

case 1: if the agent in the $(S,I)$ tuple in which the infection occurs is, itself, neighbour in an $(S,S)$ tuple, then the original tuple will become $(I,I)$ and the second one will become $(S,I)$ . As a result, $SI(t)$ will remain unchanged.
case 2: if the agent in the $(S,I)$ tuple in which the infection occurs is, itself, neighbour in an $(I,S)$ tuple, then the original tuple will become $(I,I)$ and the second one will become $(I,I)$ as well. As a result, $SI(t)$ will decrease by one.

To compute the expected change of $SI(t)/S(t)$ , we furthermore need to evaluate the probabilities for cases 1 and 2, which can be done by contrasting the numbers of $(I,S)$ and $(S,S)$ tuples:

P(\text{case 1})=\frac{SS(t)}{SS(t)+IS(t)},\ P(\text{case 2})=\frac{IS(t)}{SS(t)+IS(t)}.

(9)

For symmetry reasons, the number of $(S,I)$ tuples is always equal to the number of $(I,S)$ tuples, and the sum of $(S,I)$ and $(S,S)$ tuples is equal to the total number of agents with state $S$ . Therefore

P(\text{case 1})=\frac{SS(t)}{S(t)},\ P(\text{case 2})=\frac{SI(t)}{S(t)}.

(10)

We put things together and evaluate the expected value $X$ of the ratio $SI(t)/S(t)$ after an infection:

X = P(\text{case 1})\frac{SI(t)}{S(t)-1}+P(\text{case 2})\frac{SI(t)-1}{S(t)-1}

(11)

=\frac{SS(t)}{S(t)}\frac{SI(t)}{S(t)-1}+\frac{SI(t)}{S(t)}\frac{SI(t)-1}{S(t)-1}=\frac{SI(t)}{S(t)}\frac{SS(t)+SI(t)-1}{S(t)-1}=\frac{SI(t)}{S(t)}\frac{S(t)-1}{S(t)-1}=\frac{SI(t)}{S(t)}

(12)

As a result, the ratio remains constant in expectation and, as said, does not scale with the number of infected.

It remains to determine its value. This can be done by looking at the initial configuration of the model. Since $x_0\cdot N$ agents are initially rendered infectious and randomly picked from the population, the expected total number of $(S,I)$ pairs is precisely $N\cdot (1-x_0)x_0$ whereas the expected number of susceptible agents is $N\cdot (1-x_0)$ . As a result

E\left(\frac{SI(t)}{S(t)}\right)=\frac{SI(0)}{S(0)}=x_0.

(13)

We finally put things together and setup the mean-field model:

\frac{dS}{dt}(t)=-S(t)P(S\rightarrow I)= -S(t)\alpha\frac{SI(t)}{S(t)}=-\alpha \cdot x_0\cdot S(t).

(14)

\frac{dI}{dt}(t)=\alpha x_0(N-I(t)).

(15)

References¶

Bicher, M. (2017). Classification of microscopic models with respect to aggregated system behaviour [Phdthesis]. Technische Universität Wien.