Simulating infections

Here's a little R simulation of fashion/disease infection/group identity in a population. There are two groups, red and white. Individuals move around at random. Sometimes, they change group. Whether they change depends on the number of their neighbours who are from a different group. Also, individuals who have just changed are "infectious" (or, trendsetters). Individuals who have not changed for a long time are not very infectious, and are also easy to infect.

It's fun to play with the parameters. It would also be fun to implement this with real people, as e.g. a mobile app, using real time movement.

##############################################
N <- 500
maxdist <- .1 # distance within which "neighbours" are counted
maxt <- 240   # number of repetitions
delay <- 0    # time delay between reps
step <- FALSE # ask after each step?
speed <- 0.01 # speed indivs move
decay <- 1    # speed decaycols <- c("red", "white", "pink", "yellow") # last two for changers
##############################################
 prop2changeprob <- function(props) {
  props^2
}
infectivity <- function(lifespan) 1/lifespan
dist <- function(ind, oths) sqrt((ind$x-oths$x)^2+(ind$y-oths$y)^2)
par(bg="gray")

inds <- data.frame(x=runif(N), y=runif(N), col=sample(1:2, N, replace=TRUE))
inds$oldx <- inds$x
inds$oldy <- inds$y
inds$lifespan <- 0 # since last change

for (t in 1:(maxt)) {
  inds$lifespan <- inds$lifespan + 1
  propother <- sapply(1:nrow(inds), function(i) {
    i <- inds[i,]
    close <- dist(i, inds) < maxdist
    if (any(close)) {
      wt <- infectivity(inds$lifespan)
      wt[!close] <- 0
      sum(wt * (inds$col!=i$col)) / sum(wt) # weighted by
    } else 0
  })
  change <- runif(N) < prop2changeprob(propother)
  displaycol <- inds$col
  inds$col <- ifelse(change, 3-inds$col, inds$col)
  inds$lifespan[change] <- 0
  displaycol <- ifelse(displaycol==inds$col, displaycol, inds$col+2)
  inds$oldx <- inds$x
  inds$oldy <- inds$y
  inds$x <- inds$x + runif(N, -speed,speed) + decay*(inds$x-inds$oldx)
  inds$y <- inds$y + runif(N, -speed,speed)+ decay*(inds$y-inds$oldy)
  #inds$x <- pmax(0, pmin(1, inds$x)) # square
  #inds$y <- pmax(0, pmin(1, inds$y))
  inds$x <- inds$x %% 1 # torus
  inds$y <- inds$y %% 1
 
  plot(inds$x, inds$y, col=cols[displaycol], pch=19, cex=.7, bg="grey",
    xlim=0:1, ylim=0:1)
  Sys.sleep(delay)
  if (step) {
    rl <- readline()
    if (rl=="q") stop("Quitting")
  }
}