Contact Tracing SIR

An SIR Model for Contact Tracing

Contact Tracing

To incorporate contact tracing in an SIR model, we need to track the contacts that each patient made during their infectious period. So we need an additional domain in the state of agent to track the contacts. In addition, a patient may contact a susceptible individual and not transmit the disease. In the standard SIR model, such a contact does not change the disease dynamics, and thus can be ignored by tracking only infectious contacts. That is, the transmission rate in the standard SIR model is a product of the contact rate (here denoted beta) and the transmission probability per contact. Here we need to split these two quantities, because contacted susceptible agents may be contact traced and put into quarantine, potentially (albeit slightly) affecting disease dynamics by reducing the susceptible population. We thus need a new parameter sigma to denote the susceptibility. We assume a randomly mixed population.

To incorporate the susceptibility into the simulation, we will use the to_change_callback argument when specifying a state transition caused by a contact. The callback can be an R function that takes three arguments:

• time: the current time
• agent: the agent initiating the contact
• contact the agent being contacted This callback function should return a logical value. TRUE means that the state transition will occur, and FALSE prevents it.

Thus, our callback should generate a Bernoulli random variable with the success probability alpha to determine if the transmission will occur. No matter occurring or not, we will add the contact to the contact list of the agent, and vice versa.

We assume that each patient recover with a rate gamma, and are diagnosed with a rate tau (with a state "T"). Each diagnosed infectious agent (with a state "T") triggers contact tracing after an exponentially distributed waiting time (with a rate theta), for simplicity. Because the contact tracing is only triggered once, after triggering, we change the state of the patient to "X", so that tracing will not be triggered again.

When a diagnosed patient triggers contact tracing, their contact list is examined, and each contact is traced with a probability p (called the contact tracing coverage probability). If the contact is susceptible (with a state "S"), they are put into quarantine (with a state "Q"). If they have a state "I", it is immediately diagnosed (and their state changes to "T"). Q quarantined susceptible agent leaves quarantine after exactly Tq days.

Note that this mode is a extension of the one introduced in [1], with the quarantine of susceptible contacts.

Parameter values

  N=10000 # the population size
  p = 0.2 # the contact tracing coverage probability
  beta = 0.8 # the contact rate
  sigma = 0.5 # the susceptibility, i.e., the probability of transmission per contact
  gamma = 0.1 # the recovery rate
  tau = 0.15 # the diagnose rate
  I0 = 20 # the initial number of infections
  theta = 1 # the contact tracing rate for diagnosed patients
  Tq = 14 # the quarantine period

Initialize the simulation

Here we use an initializer function to initialize each agent. The first I0 agents are initially infectious, and the remaining are susceptible. Each agent initially has an empty contract list, as disease transmission at not started yet.

We count the agent sin each of the "S", "I", "T", "X", "Q", and "S" state.

  sim = Simulation$new(N, function(i) {
    list(if (i<I0) "I" else "S", contacts=list())
  })
  sim$addLogger(newCounter("S", "S"))
  sim$addLogger(newCounter("I", "I"))
  sim$addLogger(newCounter("T", "T"))
  sim$addLogger(newCounter("X", "X"))
  sim$addLogger(newCounter("Q", "Q"))
  sim$addLogger(newCounter("R", "R"))
  m = newRandomMixing()
  sim$addContact(m)

Transitions and callbacks

Note that we use a R function to generate the quarantine period. This function takes the current time as the input, and returns the waiting time for the transition to occur.

  quarantine.period = function(time) Tq
  sim$addTransition("Q"->"S", quarantine.period)
  sim$addTransition("I"->"R", gamma)
  sim$addTransition("I"->"T", tau)
  sim$addTransition(
    "I" + "S" -> "I" + "I" ~ m, beta,
    to_change_callback = function(time, agent, contact) { 
      sa = getState(agent)
      sa$contacts = c(sa$contacts, contact)
      sc = getState(contact)
      sc$contacts = c(sc$contacts, agent)
      setState(agent, sa)
      setState(contact, sc)
      # this callback returns a boolean value indicating whether this
      # event should happen.
      runif(1) < sigma
    })
  sim$addTransition(
    "T"->"X", theta, 
    changed_callback = function(time, agent) {
      for (c in getState(agent)$contacts) {
        if (runif(1) < p) {
          if (matchState(c, "S")) {
            setState(c, "Q")
          } else if (matchState(c, "I")) {
            setState(c, "T")
          }
        }
      }
    })

Run the simulation

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

Reference

[1] Bednarski S, Cowen LLE, Ma J, Philippsen T, van den Driessche P, Wang M. A contact tracing SIR model for randomly mixed populations. J Biol Dyn. 2022 Dec;16(1):859-879. doi: 10.1080/17513758.2022.2153938. PMID: 36522826.