Gillespie SIR

Simulate an SIR model using the Gillespie method

Parameters

library(ABM)
# the population size
N = 10000
# the initial number of infectious agents
I0 = 10
# the transmission rate
beta = 0.4
# the recovery rate
gamma = 0.2

An exponential random number generator

We need this to handle the special case that the rate is 0, which should return Inf. The rexp function does not cut it.

wait.exp = function(rate) if (rate == 0) Inf else rexp(1, rate)

Create a simulation

An SIR model's state can be characterized by three values, the number of susceptible, infectious and recovered individuals, namely S, I and R.

The state of an agent is an R list. At most one of the value can be unnamed, others must be named. The names are called domains, and their values can be any R value.

sim = Simulation$new()
# the initial state
sim$state = list(S=N-I0, I=I0, R=0)

Schedule events

The Gillespie method works because it assumes exponential waiting times for all events. When the earlier event happens, and the state of the system changes, the other events become invalid. So we need to regenerate them. Because of the memoryless property, this is done by unscheduling all events, then recalculate and reschedule new events. The following function does so.

Note that an event is defined by the event time and an event handler. Different types of events can then be characterized by the event handler.

reschedule = function(time, agent) {
  clearEvents(agent)
  state = getState(agent)
  t.inf = time + wait.exp(beta*state$I*state$S/N)
  schedule(agent, newEvent(t.inf, handler.infect))
  t.rec = time + wait.exp(gamma*state$I)
  schedule(agent, newEvent(t.rec, handler.recover))
}

Events

There are two types of events for an SIR model, infection and recovery. An infection event causes I to increase by 1 and S to decrease by 1. A recovery event causes I to decrease by 1 and R to increase by 1. When either event happen, we need to reschedule the events. These are handled by the following event handlers.

An event handler take 3 arguments

• time is the current simulation time
• sim is an external pointer to the Simulation object.
• agent is the agent that the event is scheduled to

# The infection event handler
handler.infect = function(time, sim, agent) {
  x = getState(agent)
  x$S = x$S - 1
  x$I = x$I + 1
  setState(agent, x)
  reschedule(time, agent)
}
# The recovery event handler
handler.recover = function(time, sim, agent) {
  x = getState(agent)
  x$R = x$R + 1
  x$I = x$I - 1
  setState(agent, x)
  reschedule(time, agent)
}

Event loggers

Without an a logger, no observations will be made. When the state of an agent changes, it is logged by a logger. The StateLogger logs the a specific state of an agent. Creating a state logger needs three arguments

• name: the name of the observation variable
• agent: the agent which state will be logged, NULL meaning the simulation object
• state: the domain (i.e., name) of the state to log

# log the S state of the simulation.
sim$addLogger(newStateLogger(name="S", agent=NULL, state="S"))
sim$addLogger(newStateLogger("I", NULL, "I"))
sim$addLogger(newStateLogger("R", sim$get, "R"))

schedule the first event

# schedule an infection event and a recovery event
reschedule(0, sim$get)

Run the simulation and get the observations

The run method takes a vector of observation times, run the simulation, and return a data.frame of observations.

result = sim$run(0:100)
print(result)