SIR Births Deaths

An SIR Model with Births and Deaths

Since Version 0.3, the framework allows agents leaving a population. This allows two new modeling schemes:

1. Births and deaths: if an agent leaves a population without entering another population, then the agent leaves the simulation is can thus be considered dead. Note that adding an agent was already possible in v0.2 (except for networks, which is implemented in v0.3).
2. Patch models and migration: an agent can leave a population and enter another, which models migration. Note that the framework already support multiple populations in v0.2, but without migration it was not very useful, because contacts are restricted to agents within the same population.

In this example, I demonstrate how to incorporate births and deaths in the simulation. We will simulate an SIR model with natural births and deaths. The births happen at a constant rate lambda, and the each individual dies at a rate mu. The equilibrium population size is thus lambda/mu.

lambda = 100
mu = 0.01

Deaths

There are two methods to simulate the death of na agent:

1. Calling the leave() function causes the agent to the population. Te agent is not in the simulation anymore, until it enters another population, and thus can be considered dead.
2. The setDeathTime() function can be used to set the death time of an agent. When the time arrives, the agent leaves the population.

• Note that calling this function multiple times causes the agent to die at the earliest time. Thus this can be used to simulate multiple causes, m such as disease induced deaths, or natural deaths.

In this example, we will use the second method. The death times of newborn agents are set when the agents are created. See the section below. For the initial agents that are created with the simulation (see, e.g., the SEIR model example), we then assign the death time of each agent.

  for (i in 1:n) {
    a = sim$agent(i)
    setDeathTime(a, rexp(1, mu))
  }

where n is. the population size, and sim is the created population. Please see the full program below.

Births

Create a new agent using the newAgent() function, then add it to a population using the addAgent(population, agent) function.

Because the agents are born at a constant rate, we will create a birth event and a schedule it to the simulation object. The event has an exponentially distributed waiting time with the rate lambda. In the handle, we create an agent, add it to the simulation, then schedule a new birth event.

# Here time is the current time int eh simulation
birth.event = function(time) {
newEvent(time + rexp(1, lambda), function(time, sim, agent) {
  a = newAgent("S", death_time = time + rexp(1, mu))
  addAgent(sim, a)
  schedule(sim, birth.event(time))
})
}
sim$schedule(birth.event(0))

The whole program

# Parameters
gamma = 0.2 # the recovery rate
beta = 0.4 # the transmission rate
n = 10000 # the initial population size
mu = 0.01 # the natural death rate
lambda = 100 # the birth rate
I0 = 5 # the initial number of infectious agents
# A function to create a birth event, the argument "time" is the current simulation time
birth.event = function(time) {
  newEvent(time + rexp(1, lambda), function(time, sim, agent) {
    a = newAgent("S", death_time = time + rexp(1, mu))
    addAgent(sim, a)
    schedule(sim, birth.event(time))
  })
}
# Initialize simulation, create the agents, and set the death times for the agents
sim = Simulation$new(as.list(c(rep("I", I0), rep("S", n-I0))))
for (i in 1:n) {
  a = sim$agent(i)
  setDeathTime(a, rexp(1, mu))
}
# schedule a birth event
sim$schedule(birth.event(0))
# the remaining is a standard SIR model
sim$addLogger(newCounter("S", "S"))
sim$addLogger(newCounter("I", "I"))
sim$addLogger(newCounter("R", "R"))
m = newRandomMixing()
sim$addContact(m)
sim$addTransition("I"->"R", gamma)
sim$addTransition("I" + "S" -> "I" + "I" ~ m, beta)
x = sim$run(0:200)
print(x)