minextant<-10 #condition of the model (Set to 0 for No Condition, 1 for Extant, or minimum species diversity for Min or Min-Max) maxextant<-41 #Set the maximum number of species. This is only important if the Min-Max model is being used tmax<-235 #maximum number of time steps, or the age of the clade in 100,000 year time steps pcaxes<-4 #Set the number of PC axes to be used lambda<-c( 81.5723374, 37.7899324, 19.7300884, 13.6455636) #First four eigenvalues of empirical covariance matrix r0<-0.005 #Set starting morphological change rate p<-.01 #Set starting value of p q<-(p*.9) #Set starting value of q qprecent<-0.9 rparameterincrementnumber=20 #Number of r values to be explored rparameterincrement=.0075 #Size of increase of r values for each additional set of simulations pqparameterincrementnumber<-10 #Number of values of p and q to be explored pqparameterincrement<-.02 #Size of increase of p values for each additional set of simulations nrep<-1000 #number of replicates n0<-1 #starting diversity discrete.birthdeath<-function(p,q,tmax,n0,pcaxes,lambda,r) #Create a function to simulate speciation-dependent morphological evolution in discrete time { extant<-n0 #Set the number of extant species to starting diversity lineagecount<-1 #Total lineages generated for this clade (including extinct lineages) morphmatrix<-matrix(data=rep(0,4),ncol=4) #Set aside a matrix for the morphology of each lineage finalmorphmatrix<-matrix(data=rep(0,4),ncol=4) #Set aside a matrix for the morphology of each lineage extant at the end of the trial counter<-1 #Reset time step counter extinct<-numeric() #Set aside a vector to record which species are extinct extinct[1]<-0 #Set species 1 to extant (entry in the extinct vector = 0 for extant species, and the time step of extinction for extinct species) n<-vector(length=tmax) #set aside vector for standing diversity of each time step while(counter<(tmax+1)&&extant>0) #Loop to end trial when all lineages are extinct or after specified number of time steps { for(z in 1:lineagecount) #In each time step, give each species a chance to go extinct, branch, or remain the same { change<-runif(1, 0, 1) #Random number to determine what specified lineage will do in this time step if(change=minextant&&Results$extant<=maxextant) #Only record the results if they meet the conditions. #Note that for the Min-Max version of the model, this line must be if (Results$extant>=minextant&&Results$extant<=maxextant) { irep<-(irep+1) #If conditions are met, move to the next trial G[irep]<-Results$variance #Record the final variance of the clade } } G<-G[-1] #delete first simulation from data tables write(t(N),file="ectodinistandingdiversitymatrix.xls",ncolumns=tmax,append=TRUE,sep = "\t") #Write the standing diversities of all 1000 clades to file G[is.na(G)]<-0 #Set all zero variances (those of clades with only 1 or 0 extant species) to 0 VarianceResults[rparameter,]<-G #Add all final variances to the variance matrix } write(t(VarianceResults),file="EctodiniECVariancresults.xls",ncolumns=nrep,append=TRUE,sep = "\t") #Write the final variance for all clades for the set of 20 r values to file }