##################################################
#######second simulation to improve the range
lambda_v = 1 / seq(200,600, by=100)
m = 1000 ; #1000 essential genes
n_v = seq(3,10,by=1); #each essential gene interacts with n non-essential genes
p_v = seq(0.5, 0.95, by=0.005) ; #the chance that each gene interaction is active at t=0
#c = 1 / (1 - (1-p)^n ) ; #for un-normalized simulation
#t0 = (1-p)/(lambda * p)
sim_names = c( "m", "n","p", "c","lambda", "R","G","t0", "t_median")
sim2 = t(rep(NA, length(sim_names)))
sim2 = data.frame(sim2)
names(sim2) = sim_names
for( lambda in lambda_v) {
for( n in n_v ){
for( p in p_v ) {
c = 1 / (1 - (1-p)^n )
t0 = (1-p)/(lambda * p)
G =lambda * (n-1) * p/(1-p)
R = c*m*n*(p*lambda)*(1-p)^(n-1)
t_median = log(1 + log(2)*G/R)/G
ret = c(m, n, p, c, lambda, R, G, t0, t_median)
sim2 = rbind(sim2, ret)
}
}
}
sim2.good = sim2[sim2$G<0.2 & sim2$G>0.05
& sim2$R<0.01 & sim2$R>1E-4, ]
sim2.good = sim2.good[sim2.good$t_median>20 & sim2.good$t_media < 50, ]
summary(sim2.good)
my.colors = c("red","green","blue","purple","brown","black","darkgreen")
model.sim2 = lm( log10(sim2.good$R) ~ sim2.good$G )
summary(model.sim2)
plot( log10(sim2.good$R) ~ sim2.good$G, col=my.colors[sim2.good$n - 3], pch=3 )
points( log10(tb$R0) ~ tb$G, pch=19, xlim=c(0.05, 0.19), ylim=c(-3, -2.1), col="red" )
title("p=0.007, Strehler-Mildvan, natural isolates")
text( tb$G*1.02, log10(tb$R0)*1.025, tb$strain)
model = lm(log10(tb$R0) ~ tb$G )
summary( model ) #p=0.007, strehler-mildvan correlation!
abline(model, col='red')
library(RColorBrewer);
#hmcol = colorRampPalette(brewer.pal(5,"RdBu"))(8);
hmcol = colorRampPalette(brewer.pal(3,"Blues"))(4);
#sim2.good = sim2[sim2$G<0.2 & sim2$G>0.05 & sim2$R<0.01 & sim2$R>1E-4, ]
sg2.lambda = sim2.good[sim2.good$lambda == 1/300 & sim2.good$n<=7, ]
#plot( log10(tb$R0) ~ tb$G, pch=19, xlim=c(0.05, 0.19), ylim=c(-3, -2.5), col="red" )
plot( log10(tb$R0) ~ tb$G, pch=19, ylim=c(-4,-2), xlim=c(0.05,0.2), col="red" )
#lines( log10(sg2.lambda$R) ~ sg2.lambda$G, col=my.colors[sg2.lambda$n - 3] )
for( n in unique(sg2.lambda$n)) {
#tmp = sg2.lambda[ sg2.lambda$n==n & sg2.lambda$t_median>20 & sg2.lambda$t_media < 50, , ]
tmp = sg2.lambda[ sg2.lambda$n==n , ]
lines( log10(tmp$R) ~ tmp$G, col=hmcol[n-3], lty=2 )
text(max(tmp$G, na.rm=T), log10(min(tmp$R, na.rm=T)), paste('n=',n))
}
title(paste("lambda=", "1/300"))
text( tb$G*1.02, log10(tb$R0)*1.025, tb$strain)
abline(model, col='red')
lambda = 1/350 is better than 1/300.
