search

AllenDowney / ModSimPy

Text and supporting code for Modeling and Simulation in Python
https://allendowney.github.io/ModSimPy/
MIT License
822 stars 1.76k forks source link

Suggestion for a biological modelling exercise #1

Closed tobsecret closed 5 years ago

tobsecret commented 7 years ago

Hi, Saw your lightning talk at SciPy, thought I would make a suggestion, as you had invited.

I had to reproduce a bacterial simulation for my bioinformatics class and I think it is a pretty doable example.

The paper is the following: Counteraction of antibiotic production and degradation stabilizes microbial communities. Kelsic ED, Zhao J, Vetsigian K, Kishony R. Nature. 521(7553):516-9. https://www.nature.com/nature/journal/v521/n7553/full/nature14485.html

The actual panel depicting a visualization of the simulation is Figure 1D.

I could supply the simulation in python, at the moment I only have a pretty version of it in R.

ofenerci commented 7 years ago

Could you provide R version of bacterial simulation. I could work on it and try to convert it into Python.

My another suggestion for the book is to provide some simple examples in neuroscience or cardiac dynamics.

tobsecret commented 7 years ago

Fig 1D left (credit @jichen1010)

ptm <- proc.time()
#install.packages('d3heatmap')
#install.packages('installr')
#require(installr)
#install.ImageMagick()
#install.packages("animation")
#library(d3heatmap)
#library(animation)
library(gplots)

size = 100
grid = matrix(0, nrow = size, ncol = size)
anti.a = grid
anti.b = grid
anti.c = grid

wrap = function(coords){
  # coords is a matrix of coordinates of two columns
  coords = coords %% 100
  coords[which(coords == 0, arr.ind = T)] = coords[which(coords == 0, arr.ind = T)] + 100
  return(coords)
}
# A kills B, B kills C, C kills A
# Strain A = 1;
# Strain B = 2;
# Strain C = 3;
rdis = 3
rp = 5
rdeg = 2

getradius = function(radius){
  dispersal = matrix(0, nrow = 2, ncol = 1)
  for(j in 0:radius){
    for(k in 1:radius){
      if(sqrt((j-0)^2+(k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(j,k))
      }
      if(sqrt((j-0)^2+(-k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(j,-k))
      }
      if(j != 0 & sqrt((-j-0)^2+(k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(-j,k))
      }
      if(j != 0 & sqrt((-j-0)^2+(-k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(-j,-k))
      }
    }
    if(j>0 & sqrt((j-0)^2) < radius+0.5){
      dispersal = cbind(dispersal, c(j,0))
      dispersal = cbind(dispersal, c(-j,0))
    }
  }
  return(dispersal)
}

adis = getradius(rdis)
ap = getradius(rp)
adeg = getradius(rdeg)

# randomly localize three strains
sampling = sample(1:(size^2), 30)
for(i in 1:10){
  grid[sampling[3*i-2]] = 1
  grid[sampling[3*i-1]] = 2
  grid[sampling[3*i]] = 3
}

# 100 iterations in total
#windows()
for(i in 1:200){
  A = which(grid == 1, arr.ind = T)
  B = which(grid == 2, arr.ind = T)
  C = which(grid == 3, arr.ind = T)

  for(j in 1:nrow(A)){
    k.area = t(wrap(ap+A[j,]))
    anti.a[k.area] = 1
  }
  for(j in 1:nrow(B)){
    k.area = t(wrap(ap+B[j,]))
    anti.b[k.area] = 1
  }
  for(j in 1:nrow(C)){
    k.area = t(wrap(ap+C[j,]))
    anti.c[k.area] = 1
  }
  #for(k in 1:nrow(degraders)){}
  kill = function(Grid, kill.area, sensitive.number){
    anti.idx = which(kill.area == 1, arr.ind = T)
    killed.idx = anti.idx[which(Grid[anti.idx] == sensitive.number),]
    if(class(killed.idx) == 'matrix'){
      Grid[killed.idx] = 0
    } else if(class(killed.idx) == 'integer'){
      Grid[killed.idx[1], killed.idx[2]] = 0
    }
    return(Grid)
  }
  grid = kill(grid, anti.a, 2)
  grid = kill(grid, anti.b, 3)
  grid = kill(grid, anti.c, 1)
  # remove antibiotics
  anti.a = matrix(0, nrow = size, ncol = size)
  anti.b = matrix(0, nrow = size, ncol = size)
  anti.c = matrix(0, nrow = size, ncol = size)

  # colonization; probability based on abundance
  tmp = grid
  empty = which(tmp == 0, arr.ind = T)
  for(j in 1:nrow(empty)){
    dispersal = t(wrap(adis+empty[j,]))
    abund.a = length(which(tmp[dispersal] == 1))
    abund.b = length(which(tmp[dispersal] == 2))
    abund.c = length(which(tmp[dispersal] == 3))
    t.abund = abund.a+abund.b+abund.c
    if(t.abund >= 1){
      colon = sample(1:3, 1, prob = c(abund.a/t.abund, abund.b/t.abund, abund.c/t.abund))
      grid[empty[j,1], empty[j,2]] = colon
    }
  }
  if(min(grid) == 0){
    col.vector = c('black', 'red', 'yellow', 'green')
  } else if(min(grid) == 1){
    col.vector = c('red', 'yellow', 'green')
  }
  #pic.handle = paste('r_dispersal_3_nondegrading_', i, '.png', sep = "")
  #png(pic.handle)
  heatmap.2(grid, dendrogram="none", Colv = FALSE, 
            Rowv = FALSE, trace="none",
            col = col.vector,
            #col = colorRampPalette(c("black", "red","yellow", "green"))(50),
            labRow=FALSE, labCol = FALSE, key=FALSE)
  #dev.off()
}
proc.time() - ptm
tobsecret commented 7 years ago

Fig 1D right (credit @jichen1010):

library(gplots)

size = 100
grid = matrix(0, nrow = size, ncol = size)
anti.a = grid
anti.b = grid
anti.c = grid

wrap = function(coords){
  # coords is a matrix of coordinates of two columns
  coords = coords %% 100
  coords[which(coords == 0, arr.ind = T)] = coords[which(coords == 0, arr.ind = T)] + 100
  return(coords)
}
# A kills B, B kills C, C kills A
# Strain A = 1; produce = 8
# Strain B = 2; produce = 7
# Strain C = 3; produce = 9
rdis = 3
rp = 5
rdeg = 2

getradius = function(radius){
  dispersal = matrix(0, nrow = 2, ncol = 1)
  for(j in 0:radius){
    for(k in 1:radius){
      if(sqrt((j-0)^2+(k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(j,k))
      }
      if(sqrt((j-0)^2+(-k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(j,-k))
      }
      if(j != 0 & sqrt((-j-0)^2+(k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(-j,k))
      }
      if(j != 0 & sqrt((-j-0)^2+(-k-0)^2) < radius+0.5){
        dispersal = cbind(dispersal, c(-j,-k))
      }
    }
    if(j>0 & sqrt((j-0)^2) < radius+0.5){
      dispersal = cbind(dispersal, c(j,0))
      dispersal = cbind(dispersal, c(-j,0))
    }
  }
  return(dispersal)
}

adis = getradius(rdis)
ap = getradius(rp)
adeg = getradius(rdeg)

# randomly localize three strains
sampling = sample(1:(size^2), 30)
for(i in 1:10){
  grid[sampling[3*i-2]] = 1
  grid[sampling[3*i-1]] = 2
  grid[sampling[3*i]] = 3
}

# 100 iterations in total
#windows()
for(i in 1:200){
  A = which(grid == 1, arr.ind = T)
  B = which(grid == 2, arr.ind = T)
  C = which(grid == 3, arr.ind = T)

  # produce antibiotics
  for(j in 1:nrow(A)){
    k.area = t(wrap(ap+A[j,]))
    anti.a[k.area] = 1
  }
  for(j in 1:nrow(B)){
    k.area = t(wrap(ap+B[j,]))
    anti.b[k.area] = 1
  }
  for(j in 1:nrow(C)){
    k.area = t(wrap(ap+C[j,]))
    anti.c[k.area] = 1
  }

  # degrade antibiotics: C degrades A, B degrades C, A degrades B;
  for(j in 1:nrow(A)){
    d.area = t(wrap(adeg+A[j,]))
    anti.b[d.area] = 0
  }
  for(j in 1:nrow(B)){
    d.area = t(wrap(adeg+B[j,]))
    anti.c[d.area] = 0
  }
  for(j in 1:nrow(C)){
    d.area = t(wrap(adeg+C[j,]))
    anti.a[d.area] = 0
  }

  kill = function(Grid, kill.area, sensitive.number){
    anti.idx = which(kill.area == 1, arr.ind = T)
    killed.idx = anti.idx[which(Grid[anti.idx] == sensitive.number),]
    if(class(killed.idx) == 'matrix'){
      Grid[killed.idx] = 0
    } else if(class(killed.idx) == 'integer'){
      Grid[killed.idx[1], killed.idx[2]] = 0
    }
    return(Grid)
  }
  grid = kill(grid, anti.a, 2)
  grid = kill(grid, anti.b, 3)
  grid = kill(grid, anti.c, 1)
  # remove antibiotics
  anti.a = matrix(0, nrow = size, ncol = size)
  anti.b = matrix(0, nrow = size, ncol = size)
  anti.c = matrix(0, nrow = size, ncol = size)

  # colonization; probability based on abundance
  tmp = grid
  empty = which(tmp == 0, arr.ind = T)
  for(j in 1:nrow(empty)){
    dispersal = t(wrap(adis+empty[j,]))
    abund.a = length(which(tmp[dispersal] == 1))
    abund.b = length(which(tmp[dispersal] == 2))
    abund.c = length(which(tmp[dispersal] == 3))
    t.abund = abund.a+abund.b+abund.c
    if(t.abund >= 1){
      colon = sample(1:3, 1, prob = c(abund.a/t.abund, abund.b/t.abund, abund.c/t.abund))
      grid[empty[j,1], empty[j,2]] = colon
    }
  }
  if(min(grid) == 0){
    col.vector = c('black', 'red', 'yellow', 'green')
  } else if(min(grid) == 1){
    col.vector = c('red', 'yellow', 'green')
  }
  pic.handle = paste('r_dispersal_3_degrading_', i, '.png', sep = "")
  png(pic.handle)
  heatmap.2(grid, dendrogram="none", Colv = FALSE, 
            Rowv = FALSE, trace="none",
            col = col.vector,
            #col = colorRampPalette(c("black", "red","yellow", "green"))(50),
            labRow=FALSE, labCol = FALSE, key=FALSE)
  dev.off()
}
AllenDowney commented 6 years ago

Interesting! I will check it out.

(Sorry for the slow reply -- I didn't get a notification.)