Skip to content

Agent SEIR

Jump to bottom
junlingm edited this page Mar 8, 2023 · 5 revisions

An Agent Based SEIR Model

Parameter values

  gamma = newExpWaitingTime(0.2) # the recovery rate
  beta = 0.4 # the transmission rate
  sigma=0.5 # the rate leaving the latent stage and becoming infectious
  n = 10000 # the population size

Create a Simulation Object and Initialize the Agents

A Simulation is represented by the R6 class Simulation. The constructor takes a single argument: the number of agents. If it is not given, then it starts with no agents (yet a simulation itself is an agent).

If a state is specified as an R value that is not a list, it is equivalent to a list with the value as the unnamed value.

The following create a simulation with n agents, and initialize the states of the agent (the first 5 as list("I"), and the following are list("S")).

Since version 0.2, the constructor of Simulation accepts an initializer, which is a function taking an agent index (starting from 1) and return the initial state of the agent.

  sim = Simulation$new(n, function(i) if (i <= 5) "I" else "S")

Counters

Without a logger, there will be no observations. A counter is a special type of logger that counts state changes. It take three arguments

  • name: the name of the counter, which is the name of the observation variable
  • from: its meaning depends on the following argument.
  • to: optional. If to is NULL, then counter counts the number of agents in the state that matches from. Otherwise, it logs the transitions that agent jumped from a state that matched "from" to one matched "to".

A state matches if the the corresponding domains in "from" (or "to") have identical values as the agent's state.

  sim$addLogger(newCounter("S", "S"))
  sim$addLogger(newCounter("I", "I"))
  sim$addLogger(newCounter("R", "R"))

Contacts

Contacts of agents are managed by the Contact objects. These objects defined the contacts of an agent. We create a RandomMixing object which is a subclass of Contact. When an agent initialize the contact, it randomly select an agent in the population as the contact.

The contact objects must be added to the simulation.

  m = newRandomMixing()
  sim$addContact(m)

Transitions

There are two types of transitions.

  • Spontaneous transitions, e.g., showing system or recovery.
  • Transitions caused by contacts.

The transition rules are specified by a formula:

  • "->" denotes the direction of transition.
  • A contact transition needs to have two agents state specified, the agent who initiates the contact, and the contacting agent. There states are specified by the "+" operator. The first is the state of the initiator, the second is the state of the contact.
  • A contact transition is always associated with a contact object, which determines the contacts of an agent. The contact is specified after "~".

Each transition has to have a waiting time, which must be exponentially distributed. It can be a number specifying the rate. Though not used in the example, a transition can have a to_change_callback and a changed_callback. The first callback is a function that return a logical value. If TRUE, the transition will occur. Otherwise the transition will be prevented. The latter callback is executed after the state change happened.

The following defines two spontaneous state transitions

  • becoming infectious: "I" -> "R"
  • recovery: "I" -> "R" and a contact transition
  • "I" +"S" -> "I" + "E" ~ m, which means an agent with a state "I" contacts an agent with state "S", the first agent's state remains "I" while the contact's state becomes "E". The contact rule is specified by "m".
  sim$addTransition("E"->"I", sigma)
  sim$addTransition("I"->"R", gamma)
  sim$addTransition("I" + "S" -> "I" + "E" ~ m, beta) #, changed_callback=changed)

Disease controls

The transmission rate may be a time dependent function when disease control measures are considers. This is why the waiting times take the current time as input. The following example defines a waiting time which corresponds to a step function in the transmission rate beta at time T, and the reduction factor is p.

  T = 15 # time when the control measure is implemented
  p = 0.5 # the reduction in transmission rate due to control measure
  beta.t = function(time) {
    if (time>=T) return(rexp(1, beta * p))
    t = rexp(1, beta1)
    if (time + t <= T) t else T - time + rexp(1, beta * p)
  }

Then the last sim$addTransition call above is replaced by

  sim$addTransition("I" + "S" -> "I" + "E" ~ m, beta.t) #, changed_callback=changed)

Run the simulation

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

  result = sim$run(0:200)
  print(result)
Clone this wiki locally