Comments (4)
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.
from modsimpy.
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
from modsimpy.
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()
}
from modsimpy.
Interesting! I will check it out.
(Sorry for the slow reply -- I didn't get a notification.)
from modsimpy.
Related Issues (20)
- update_func() with UNITS (Pint) an Frame HOT 2
- Issue: Multiple assignment of State in code HOT 1
- Installing Pint with Python 3.7 HOT 3
- modsim install instructions incorrect in version 3 book HOT 11
- Errors on importing library HOT 4
- Decode Error, con't decode with "gbk"
- Cannot run dir() on ModSimSeries objects HOT 1
- Can't pip install modsimpy - Error Message HOT 3
- Convert notebooks to LaTeX
- Open in Colab links HOT 1
- Installing pint with Python 3.8 HOT 2
- chapter 2 'plot' is not defined HOT 3
- Pandas dataframe doesn't display on website HOT 2
- No figure of spiderman on website HOT 1
- Math rendering error on website HOT 1
- Math dollar signs rendered on yoyo example on website HOT 1
- Ch 01 velocity magnitude HOT 1
- Deprecated use of `pd.Series.iteritems` HOT 1
- Issue with Pint - Chapter 1 HOT 3
- Is `modsim/test_modsim.py` still a valid unit test HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from modsimpy.