Changes of n, p, lambda and the Strehler-Mildvan correlation in yeast natural isolates

The red line is the observed lnR-G correlation in yeast natural isolates. Effect on R and G by changing n, p, and lambda are plotted in the background. Each continuous curve is colored by n.



##################################################
#######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')

The above plot is conditioned on lambda=1/300. It makes that large n leads to smaller R and G. The dashed line care caused by change of p.

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.