Converting a math problem into a discrete-event simulation

Question

I am trying to implement the SIR epidemic model. Basically, the way I understand the model is that at each time step, it tells us how many nodes get infected and how many nodes are recovered. I am now trying to convert this into an discrete-event simulation and am confused at a point.

The model is generally solved using Euler's method. Now, when I am converting it into a discrete event simulation, I am doing something like this (the numbers are used for clarity):

Initialize 100 members
At every time step t,
  //Determine how many get infected
  for i = 1 to 100
     let i pass a message to its neighbors
     When the neighbor receives the message from an infected member, it generates a random number and if it is less than beta*(infected/total), where beta is the infection rate, then the member gets infected
     Update the count for infected, recovered, susceptible

  //Determine how many are recovered
  for i = 1 to 100
      Generate a random number from a uniform distribution and check if it is less than gamma*infected. If it is, then this member is recovered.
      Update the count for infected, recovered, susceptible

I was essentially wondering if the above approach is right. Any suggestions?

Answer 1

I see at least two problems here, both stemming from the same root cause.

1) If you introduce some kind of a neighborhood structure, then beta*(infected/total) should be just beta.

2) The line that has gamma*infected apples only to infected individuals, and you don't need gamma*infected -- it should be just gamma.

Maybe there are more problems -- it's difficult to say from the pseudo-code as it is. You are welcome to refine it a bit, and I'll have another look.

Answer 2

Looks pretty good as a start, except that for the first loop you need to remember that only susceptible individuals can become infected, and for the second one, that only infected individuals can become recovered. I also believe the probabilities of transition on each event (susceptible receiving message from infected neighbor, infected possibly recovering) are not a function of the current number of infected individuals -- they're constants (I think you're misleading the "mass effect" probabilities as applying to each individual episode at a time step -- they don't).

A bit subtler is how you do the first loop (not obvious to me from the SIR model): I think you want to determine all "messages" first, then which ones cause transitions susceptible -> infected -- i.e., two loops rather than one -- because an individual that's just gotten infected at this time step cannot infect yet others in this same time step but only in the future; also, the transition infected -> recovered is not possible for individual who just became infected at this very time step, so you'd have to arrange your loops a bit differently!

Consider modeling each individual with two "state" attributes:

-- nummsgs, number of "messages" received this time step
-- compartment (susceptible, infected or recovered)

as well as a fixed set of neighbors. Then:

for each individual:
    if individual.compartment != infected:
        continue
    for each neighbor of the individual:
        neighbor.nummsgs += 1
    if (random number says so):
        individual.compartment = recovered

for each individual:
    if individual.compartment != susceptible:
        continue
    maybe (depending on random number & nummsgs):
        individual.compartment = infected

for each individual:
    individual.nummsgs = 0

this seems to better capture the overall flow (net of collecting and logging overall counts, which you can conceptually do as part of the last loop).

Answer 3

First, this approach is not usually called discrete-event simulation, but time-step simulation. It's not incorrect, strictly-speaking, but in practice the 'discrete-event' term is used for simulations where time between events is variable. For example, if each infected node would generate "infect" message to a random neighbor once in a random time.

Second, your code misses priority or weights of messages. What if node receives both infect and recover messages? Should probability of recovery really be the same for isolated node and node in the middle of infected blob? It would lead to a kind of jitter of internal nodes between infected and healthy states.

Actually, I don't see a reason to use messages here. Probability of infection of each node is simply a function of number of it's infected neighbors - why not use just that? Depending on your goals, it might be better to replace probabilities with mathematical expectations.

Finally, the problem is closely related to the famous 'Game of Life'. Take a look at various implementations.