R code for a dynamic swarm simulation and exporting frames to make an animation.
As a part of our BIO 479: Introduction to Complexity Science course this semester, we have been learning about networks, collective behavior, swarm intelligence, and programming using R and the tidyverse.
This simulation started out to address a more complicated question about how honeybee colonies use quorum sensing during collective decision making processes, but the basics of figuring out how to model a group of animals moving in space turned into a neat visualization on its own.
Moving past an earlier draft based on brownian motion, there were a lot of decisions to make so that the animal trajectories were more biologically realistic. In this draft of the code, animals move following changes in direction and speed randomly sampled from an exponential distribution. The animal movements are restricted to a bounded area and at this point there are no interactions.
I shared the animation on twitter and in addition to a lot of positive feedback, many folks were interested in the code. Unfortunately I haven’t learned github or other code sharing tools yet, so I’m just sharing it here. I tried to post the source.R script file but couldn’t make it available through this wix site.
The code was drafted relatively quickly and isn’t meant to be highly efficient or elegant in any way, but feedback is certainly welcome.
# Draft simulation of quorum-sensing collective behavior
# Introduction to Complexity Science
# April 14, 2017
# Goal:
# Simulate a group of bees, moving randomly in two dimensions.
# If and when a quorum is reached, the bees should aggregate.
# Density-dependent positive and negative feedback regulate speed.
# Libraries needed
library(tidyverse)
library(gridExtra)
library(RColorBrewer)
# number of bees
n <- 50
# some starting coordinates and angles to sample the initial conditions from
coords <- 475:525
angles <- 0:360
# making a data frame to store the tracking information
tracks <- data.frame(bee.id=NA, time=NA, bee.x=NA, bee.y=NA, bee.theta=NA, bee.speed=NA)
# seeding the first time point with random sampled values
tracks[1:n,] <- c(1:n, rep(1, n), sample(coords, n), sample(coords, n), sample(angles, n), rep(10, n))
tracks
# Exponential distribution
# These values will be used to determine how much the
# direction and speed change each time step
dist.vals <- rexp(1000, .1)
dist.vals <- dist.vals/max(dist.vals)
hist(dist.vals, breaks=100)
# Making the bees fly
# In each run of the loop the next time step is calculated
# Based on a speed and a direction, each of which can
# change from one time step to the next, a new coordinate
# is calculated and saved as a new row.
# How long should the simulation run?
time.steps <- 400
# How to account for all the bees?
# Trick: two for-loops
# the first loop (i) goes through each time step
# the second loop (j) goes through each bee
# Also: edited the foraging range to [0, 1000] in x- and y-coordinates
for(i in 2:time.steps){
for(j in 1:n){
last.track <- tracks[tracks$time == i-1 & tracks$bee.id == j,]
time <- i
delta.theta <- 180*sample(dist.vals, 1)*sample(c(-1,1),1)
delta.speed <- 2*sample(dist.vals, 1)*sample(c(-1,1),1)
bee.theta <- last.track$bee.theta + delta.theta
bee.speed <- last.track$bee.speed + delta.speed
bee.x.new <- last.track$bee.x + bee.speed*cos(bee.theta*pi/180)
bee.y.new <- last.track$bee.y + bee.speed*sin(bee.theta*pi/180)
# Check to see if bees go out of bounds
# bounds: from 0:1000 on x and y axes
# If new coordinate falls outside of bounds, then
# the bee stays in its current spot until next timestep
bee.x <- ifelse(bee.x.new > 0 & bee.x.new < 1000, bee.x.new, last.track$bee.x)
bee.y <- ifelse(bee.y.new > 0 & bee.y.new < 1000, bee.y.new, last.track$bee.y)
new.row <- c(j, time, bee.x, bee.y, bee.theta, bee.speed)
tracks <- rbind(tracks, new.row)
}
}
head(tracks)
# View trajectories
ggplot(tracks, aes(x=bee.x, y=bee.y, color=time)) + geom_path(aes(group=bee.id))
ggplot(tracks, aes(x=bee.x, y=bee.y, color=time, alpha=time)) + geom_path(aes(group=bee.id))
# Notes:
# Can play with the parameters above to change trajectory characteristics
# e.g. could make delta.theta always positive
# e.g. could change the curviness of the paths by changing scale (180)
# Looking at how speed and direction change over time for each bee
ggplot(tracks, aes(x=time, y=bee.speed, color=bee.id)) + geom_line(aes(group=bee.id))
ggplot(tracks, aes(x=time, y=bee.theta, color=bee.id)) + geom_line(aes(group=bee.id))
# Color housekeeping
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
r <- rf(32)
# Animating the path (and saving, if you want to make a movie)
# Set your working directory (different on your computer)
# setwd("")
# Save the frames of the movie as an image sequence
# to build animation in Quicktime
for(i in 1:time.steps){
p.map <- ggplot(tracks[tracks$time <= i,], aes(x=bee.x, y=bee.y, color=bee.id, alpha=time)) + geom_path(aes(group=bee.id)) + annotate("point", x=tracks$bee.x[tracks$time == i], y=tracks$bee.y[tracks$time == i]) + guides(color=F, alpha=F)
p.density <- ggplot(tracks[tracks$time == i,], aes(x=bee.x, y=bee.y)) + stat_bin2d(bins=10, aes(alpha=log(..count..))) + scale_fill_gradientn(colours=r) + guides(fill=F, alpha=F) + xlim(0, 1000) + ylim(0, 1000) + coord_fixed() + theme_void() # + ggtitle("Density heatmap")
p.hist <- ggplot(tracks[tracks$time <= i,], aes(bee.speed)) + geom_histogram() + theme_bw()
p.speed <- ggplot(tracks[tracks$time <= i,], aes(x=time, y=bee.speed)) + geom_line(aes(color=as.factor(bee.id), alpha=time)) + guides(color=F, alpha=F) + theme_bw()
bot.grid <- grid.arrange(p.hist, p.speed, ncol=2)
side.grid <- grid.arrange(p.density, bot.grid, ncol=1)
p <- grid.arrange(p.map, side.grid, ncol=2)
png(filename=paste("plot", i, ".png", sep=""), width=825, height=400)
grid.arrange(p.map, side.grid, ncol=2)
dev.off()
}