Generate a random phylogenetic tree via simulation of a Poissonian speciation/extinction (birth/death) process. New species are added (born) by splitting of a randomly chosen extant tip. Per-capita birth and death rates evolve under a Brownian motion model constrained to a finite interval (via reflection). Thus, the probability rate of a tip splitting or going extinct depends on the tip, with closely related tips having more similar per-capita birth and death rates.
generate_tree_with_evolving_rates(parameters = list(),
max_tips = NULL,
max_time = NULL,
max_time_eq = NULL,
coalescent = TRUE,
as_generations = FALSE,
Nsplits = 2,
tip_basename = "",
node_basename = NULL,
include_birth_times = FALSE,
include_death_times = FALSE,
include_rates = FALSE)A named list specifying the birth-death model parameters, with one or more of the following entries:
birth_rate_diffusivity: Non-negative number. Diffusivity constant for the Brownian motion model of the evolving per-capita birth rate. In units 1/time^3. See simulate_bm_model for an explanation of the diffusivity parameter.
min_birth_rate_pc:
Non-negative number. The minimum allowed per-capita birth rate of a clade. In units 1/time. By default this is 0.
max_birth_rate_pc:
Non-negative number. The maximum allowed per-capita birth rate of a clade. In units 1/time. By default this is 1.
death_rate_diffusivity: Non-negative number. Diffusivity constant for the Brownian motion model of the evolving per-capita death rate. In units 1/time^3. See simulate_bm_model for an explanation of the diffusivity parameter.
min_death_rate_pc:
Non-negative number. The minimum allowed per-capita death rate of a clade. In units 1/time. By default this is 0.
max_death_rate_pc:
Non-negative number. The maximum allowed per-capita death rate of a clade. In units 1/time. By default this is 1.
rarefaction:
Numeric between 0 and 1. Rarefaction to be applied at the end of the simulation (fraction of random tips kept in the tree).
Note that if coalescent==FALSE, rarefaction may remove both extant as well as extinct clades. Set rarefaction=1 to not perform any rarefaction.
Maximum number of tips of the tree to be generated. If coalescent=TRUE, this refers to the number of extant tips. Otherwise, it refers to the number of extinct + extant tips. If NULL or <=0, the number of tips is unlimited (so be careful).
Maximum duration of the simulation. If NULL or <=0, this constraint is ignored.
Maximum duration of the simulation, counting from the first point at which speciation/extinction equilibrium is reached, i.e. when (birth rate - death rate) changed sign for the first time. If NULL or <0, this constraint is ignored.
Logical, specifying whether only the coalescent tree (i.e. the tree spanning the extant tips) should be returned. If coalescent==FALSE and the death rate is non-zero, then the tree may include non-extant tips (i.e. tips whose distance from the root is less than the total time of evolution). In that case, the tree will not be ultrametric.
Logical, specifying whether edge lengths should correspond to generations. If FALSE, then edge lengths correspond to time.
Integer greater than 1. Number of child-tips to generate at each diversification event. If set to 2, the generated tree will be bifurcating. If >2, the tree will be multifurcating.
Character. Prefix to be used for tip labels (e.g. "tip."). If empty (""), then tip labels will be integers "1", "2" and so on.
Character. Prefix to be used for node labels (e.g. "node."). If NULL, no node labels will be included in the tree.
Logical. If TRUE, then the times of speciation events (in order of occurrence) will also be returned.
Logical. If TRUE, then the times of extinction events (in order of occurrence) will also be returned.
Logical. If TRUE, then the bper-capita birth & death rates of all tips and nodes will also be returned.
A named list with the following elements:
Logical, indicating whether the simulation was successful. If FALSE, an additional element error (of type character) is included containing an explanation of the error; in that case the value of any of the other elements is undetermined.
A rooted bifurcating (if Nsplits==2) or multifurcating (if Nsplits>2) tree of class "phylo", generated according to the specified birth/death model.
If coalescent==TRUE or if all death rates are zero, and only if as_generations==FALSE, then the tree will be ultrametric. If as_generations==TRUE and coalescent==FALSE, all edges will have unit length.
Numeric, giving the time at which the tree's root was first split during the simulation.
Note that if coalescent==TRUE, this may be later than the first speciation event during the simulation.
Numeric, giving the final time at the end of the simulation. If coalescent==TRUE, then this may be greater than the total time span of the tree (since the root of the coalescent tree need not correspond to the first speciation event).
Numeric, giving the first time where the sign of (death rate - birth rate) changed from the beginning of the simulation, i.e. when speciation/extinction equilibrium was reached. May be infinite if the simulation stoped before reaching this point.
Total number of birth events (speciations) that occurred during tree growth. This may be lower than the total number of tips in the tree if death rates were non-zero and coalescent==TRUE, or if Nsplits>2.
Total number of deaths (extinctions) that occurred during tree growth.
Numeric vector, listing the times of speciation events during tree growth, in order of occurrence. Note that if coalescent==TRUE, then speciation_times may be greater than the phylogenetic distance to the coalescent root.
Numeric vector, listing the times of extinction events during tree growth, in order of occurrence. Note that if coalescent==TRUE, then speciation_times may be greater than the phylogenetic distance to the coalescent root.
Numeric vector, listing the per-capita birth rate of each tip and node in the tree. The length of an edge in the tree was thus drawn from an exponential distribution with rate equal to the per-capita birth rate of the child tip or node.
Numeric vector, listing the per-capita death rate of each tip and node in the tree.
If max_time==NULL, then the returned tree will always contain max_tips tips. In particular, if at any moment during the simulation the tree only includes a single extant tip, the death rate is temporarily set to zero to prevent the complete extinction of the tree. If max_tips==NULL, then the simulation is ran as long as specified by max_time. If neither max_time nor max_tips is NULL, then the simulation halts as soon as the time exceeds max_time or the number of tips (extant tips if coalescent is TRUE) exceeds max_tips. If max_tips!=NULL and Nsplits>2, then the last diversification even may generate fewer than Nsplits children, in order to keep the total number of tips within the specified limit.
The per-capita birth and death rates of the root are chosen randomly, independently and uniformly within their allowed intervals.
D. J. Aldous (2001). Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Statistical Science. 16:23-34.
# NOT RUN {
# Generate tree
parameters = list(birth_rate_diffusivity = 1,
min_birth_rate_pc = 1,
max_birth_rate_pc = 2,
death_rate_diffusivity = 0.5,
min_death_rate_pc = 0,
max_death_rate_pc = 1)
simulation = generate_tree_with_evolving_rates(parameters,max_tips=1000,include_rates=TRUE)
tree = simulation$tree
Ntips = length(tree$tip.label)
# plot per-capita birth & death rates of tips
plot( x=simulation$birth_rates_pc[1:Ntips],
y=simulation$death_rates_pc[1:Ntips],
type='p',
xlab="pc birth rate",
ylab="pc death rate",
main="Per-capita birth & death rates across tips",
las=1)
# }
Run the code above in your browser using DataLab