## CLADE CROWN AGE SIMULATION # Requires R-packages Ape, Geiger and TreeSim # Group birth/death parameters (b/d) determined by stem age and species diversity (n) with given extinction fraction (exf = d/b = 0.33). # r = b-d = 1/t*ln(n(1-exf)+exf) # Output: lists of i simulations crown age per group. # Input birth/death parameters need to be set, raw output file of stem ages needs to be summarized. Pyg <- 0 Car <- 0 Cre <- 0 Dip <- 0 Sph <- 0 Aga <- 0 Bli <- 0 Eug <- 0 Geh <- 0 Var <- 0 Ela <- 0 Pyt <- 0 Ege <- 0 i <- 1 while(i<1001) { p1 <- sim.bdtree(b=0.1879, d=0.062, stop='taxa', n=90, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='AgaIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Aga[i] <- bt[89] p1 <- sim.bdtree(b=0.1769, d=0.0584, stop='taxa', n=50, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='BliIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Bli[i] <- bt[49] p1 <- sim.bdtree(b=0.0683, d=0.0225, stop='taxa', n=30, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='CarIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Car[i] <- bt[29] p1 <- sim.bdtree(b=0.0627, d=0.0207, stop='taxa', n=12, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='CreIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Cre[i] <- bt[11] p1 <- sim.bdtree(b=0.1511, d=0.0499, stop='taxa', n=96, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='DipIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Dip[i] <- bt[95] p1 <- sim.bdtree(b=0.2548, d=0.0841, stop='taxa', n=50, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='EgeIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Ege[i] <- bt[49] p1 <- sim.bdtree(b=0.3786, d=0.1249, stop='taxa', n=97, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='ElaIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Ela[i] <- bt[96] p1 <- sim.bdtree(b=0.2674, d=0.0882, stop='taxa', n=124, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='EugIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Eug[i] <- bt[123] p1 <- sim.bdtree(b=0.1859, d=0.0614, stop='taxa', n=30, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='GehIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Geh[i] <- bt[29] p1 <- sim.bdtree(b=0.0766, d=0.0253, stop='taxa', n=48, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='PygIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Pyg[i] <- bt[47] p1 <- sim.bdtree(b=0.2171, d=0.0716, stop='taxa', n=14, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='PytIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Pyt[i] <- bt[13] p1 <- sim.bdtree(b=0.2195, d=0.0724, stop='taxa', n=276, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='SphIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Sph[i] <- bt[275] p1 <- sim.bdtree(b=0.2234, d=0.0737, stop='taxa', n=24, extinct=FALSE) p2 <- drop.extinct(p1) write.tree(p2, file='VarIGstembd33_Gtrees.phy', append=TRUE, digits=6) bt <- getx(p2) Var[i] <- bt[23] i <- i+1 } write.table(Aga, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=FALSE, col.names=FALSE) write.table(Bli, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Car, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Cre, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Dip, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Ege, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Ela, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Eug, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Geh, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Pyg, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Pyt, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Sph, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE) write.table(Var, file="IGstembd33_Gtimes.txt", row.names=FALSE, append=TRUE, col.names=FALSE)