simulate_deterministic_hbd: Simulate a deterministic homogenous birth-death model.

In castor: Efficient Phylogenetics on Large Trees

Simulate a deterministic homogenous birth-death model.

Description

Given a homogenous birth-death (HBD) model, i.e., with speciation rate \lambda, extinction rate \mu and sampling fraction \rho, calculate various deterministic features of the model backwards in time, such as the total diversity over time. The speciation and extinction rates may be time-dependent. “Homogenous” refers to the assumption that, at any given moment in time, all lineages exhibit the same speciation/extinction rates (in the literature this is sometimes referred to simply as “birth-death model”; Morlon et al. 2011). “Deterministic” refers to the fact that all calculated properties are completely determined by the model's parameters (i.e. non-random), as if an infinitely large tree was generated (aka. “continuum limit”).

Alternatively to \lambda, one may provide the pulled diversification rate (PDR; Louca et al. 2018) and the speciation rate at some fixed age, \lambda(\tau_o). Similarly, alternatively to \mu, one may provide the ratio of extinction over speciation rate, \mu/\lambda. In either case, the time-profiles of \lambda, \mu, \mu/\lambda or the PDR are specified as piecewise polynomial functions (splines), defined on a discrete grid of ages.

Usage

simulate_deterministic_hbd(LTT0,
                           oldest_age,
                           age0           = 0,
                           rho0           = 1,
                           age_grid       = NULL,
                           lambda         = NULL,
                           mu             = NULL,
                           mu_over_lambda = NULL,
                           PDR            = NULL,
                           lambda0        = NULL,
                           splines_degree = 1,
                           relative_dt    = 1e-3,
                           allow_unreal   = FALSE)

Arguments

LTT0	The assumed number of sampled extant lineages at age0, defining the necessary initial condition for the simulation. If the HBD model is supposed to describe a specific timetree, then LTT0 should correspond to the number of lineages in the tree ("lineages through time") at age age0.
oldest_age	Strictly positive numeric, specifying the oldest time before present (“age”) to include in the simulation.
age0	Non-negative numeric, specifying the age at which LTT0, lambda0 and rho are given. Typically this will be 0, i.e., corresponding to the present.
rho0	Numeric between 0 (exclusive) and 1 (inclusive), specifying the sampling fraction of the tree at age0, i.e. the fraction of lineages extant at age0 that are included in the tree (aka. "rarefaction"). Note that if rho0<1, lineages extant at age0 are assumed to have been sampled randomly at equal probabilities. Can also be NULL, in which case rho0=1 is assumed.
age_grid	Numeric vector, listing discrete ages (time before present) on which either \lambda and \mu, or the PDR and \mu, are specified. Listed ages must be strictly increasing, and must cover at least the full considered age interval (from age0 to oldest_age). Can also be NULL or a vector of size 1, in which case the speciation rate, extinction rate and PDR are assumed to be time-independent.
lambda	Numeric vector, of the same size as age_grid (or size 1 if age_grid==NULL), listing speciation rates (\lambda, in units 1/time) at the ages listed in age_grid. Speciation rates should be non-negative, and are assumed to vary polynomially between grid points (see argument splines_degree). If NULL, then PDR and lambda0 must be provided.
mu	Numeric vector, of the same size as age_grid (or size 1 if age_grid==NULL), listing extinction rates (\mu, in units 1/time) at the ages listed in age_grid. Extinction rates should be non-negative, and are assumed to vary polynomially between grid points (see argument splines_degree). Either mu or mu_over_lambda must be provided, but not both.
mu_over_lambda	Numeric vector, of the same size as age_grid (or size 1 if age_grid==NULL), listing the ratio of extinction rates over speciation rates (\mu/\lambda) at the ages listed in age_grid. These ratios should be non-negative, and are assumed to vary polynomially between grid points (see argument splines_degree). Either mu or mu_over_lambda must be provided, but not both.
PDR	Numeric vector, of the same size as age_grid (or size 1 if age_grid==NULL), listing pulled diversification rates (in units 1/time) at the ages listed in age_grid. PDRs can be negative or positive, and are assumed to vary polynomially between grid points (see argument splines_degree). If NULL, then lambda must be provided.
lambda0	Non-negative numeric, specifying the speciation rate (in units 1/time) at age0. Either lambda0 or lambda must be provided, but not both.
splines_degree	Integer, either 0,1,2 or 3, specifying the polynomial degree of the provided lambda, mu and PDR between grid points in age_grid. For example, if splines_degree==1, then the provided lambda, mu and PDR are interpreted as piecewise-linear curves; if splines_degree==2 they are interpreted as quadratic splines; if splines_degree==3 they are interpreted as cubic splines. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on.
relative_dt	Strictly positive numeric (unitless), specifying the maximum relative time step allowed for integration over time. Smaller values increase integration accuracy but increase computation time. Typical values are 0.0001-0.001. The default is usually sufficient.
allow_unreal	Logical, specifying whether HBD models with unrealistic parameters (e.g., negative \mu) should be supported. This may be desired for example when examining model congruence classes with negative \mu.

Details

This function supports the following alternative parameterizations of HBD models:

• Using the speciation rate \lambda and extinction rate \mu.

• Using the speciation rate \lambda and the ratio \mu/\lambda.

• Using the pulled diversification rate (PDR), the extinction rate and the speciation rate given at some fixed age0 (i.e. lambda0).

• Using the PDR, the ratio \mu/\lambda and the speciation rate at some fixed age0.

The PDR is defined as PDR = \lambda-\mu+\lambda^{-1}d\lambda/d\tau, where \tau is age (time before present). To avoid ambiguities, only one of the above parameterizations is accepted at a time. The sampling fraction at age0 (i.e., rho0) should always be provided; setting it to NULL is equivalent to setting it to 1.

Note that in the literature the sampling fraction usually refers to the fraction of lineages extant at present-day that have been sampled (included in the tree); this present-day sampling fraction is then used to parameterize birth-death cladogenic models. The sampling fraction can however be generalized to past times, by defining it as the probability that a lineage extant at any given time is included in the tree. The simulation function presented here allows for specifying this generalized sampling fraction at any age of choice, not just present-day.

The simulated LTT refers to a hypothetical tree sampled at age age_grid[1], i.e. LTT(t) will be the number of lineages extant at age t that survived until age age_grid[1] and have been sampled, given that the fraction of sampled extant lineages at age0 is rho0. Similarly, the returned Pextinct(t) (see below) is the probability that a lineage extant at age t would not survive until age_grid[1]. The same convention is used for the returned variables Pmissing, shadow_diversity, PER, PSR, SER and PND.

Value

A named list with the following elements:

success	Logical, indicating whether the calculation was successful. If FALSE, then the returned list includes an additional 'error' element (character) providing a description of the error; all other return variables may be undefined.
ages	Numerical vector of size NG, listing discrete ages (time before present) on which all returned time-curves are specified. Listed ages will be in ascending order, will cover exactly the range age_grid[1] - oldest_age, may be irregularly spaced, and may be finer than the original provided age_grid. Note that ages[1] corresponds to the latest time point (closer to the tips), while ages[NG] corresponds to the oldest time point (oldest_age).
total_diversity	Numerical vector of size NG, listing the predicted (deterministic) total diversity (number of extant species, denoted N) at the ages given in ages[].
shadow_diversity	Numerical vector of size NG, listing the predicted (deterministic) “shadow diversity” at the ages given in ages[]. The shadow diversity is defined as N_s=N\cdot \rho(\tau_o)\lambda(\tau_o)/\lambda, where \tau_o is age0.
Pmissing	Numeric vector of size NG, listing the probability that a lineage, extant at a given age, will be absent from the tree either due to extinction or due to incomplete sampling.
Pextinct	Numeric vector of size NG, listing the probability that a lineage, extant at a given age, will be fully extinct at present. Note that always Pextinct<=Pmissing.
LTT	Numeric vector of size NG, listing the number of lineages represented in the tree at any given age, also known as “lineages-through-time” (LTT) curve. Note that LTT at age0 will be equal to LTT, and that values in LTT will be non-increasing with age.
lambda	Numeric vector of size NG, listing the speciation rate (in units 1/time) at the ages given in ages[].
mu	Numeric vector of size NG, listing the extinction rate (in units 1/time) at the ages given in ages[].
diversification_rate	Numeric vector of size NG, listing the net diversification rate (\lambda-\mu) at the ages given in ages[].
PDR	Numeric vector of size NG, listing the pulled diversification rate (PDR, in units 1/time) at the ages given in ages[].
PND	Numeric vector of size NG, listing the pulled normalized diversity (PND, in units 1/time) at the ages given in ages[]. The PND is defined as PND=(N/N(\tau_o))\cdot\lambda(\tau_o)/\lambda.
SER	Numeric vector of size NG, listing the “shadow extinction rate” (SER, in units 1/time) at the ages given in ages[]. The SER is defined as SER=\rho(\tau_o)\lambda(\tau_o)-PDR.
PER	Numeric vector of size NG, listing the “pulled extinction rate” (PER, in units 1/time) at the ages given in ages[]. The PER is defined as SER=\lambda(\tau_o)-PDR (Louca et al. 2018).
PSR	Numeric vector of size NG, listing the “pulled speciation rate” (PSR, in units 1/time) at the ages given in ages[]. The PSR is defined as PSR=\lambda\cdot(1-Pmissing).
rholambda0	Non-negative numeric, specifying the product of the sampling fraction and the speciation rate at age0, \rho\cdot\lambda(\tau_o).

Author(s)

Stilianos Louca

References

H. Morlon, T. L. Parsons, J. B. Plotkin (2011). Reconciling molecular phylogenies with the fossil record. Proceedings of the National Academy of Sciences. 108:16327-16332.

S. Louca et al. (2018). Bacterial diversification through geological time. Nature Ecology & Evolution. 2:1458-1467.

See Also

loglikelihood_hbd

Examples

# define an HBD model with exponentially decreasing speciation/extinction rates
Ntips     = 1000
beta      = 0.01 # exponential decay rate of lambda over time
oldest_age= 10
age_grid  = seq(from=0,to=oldest_age,by=0.1) # choose a sufficiently fine age grid
lambda    = 1*exp(beta*age_grid) # define lambda on the age grid
mu        = 0.2*lambda # assume similarly shaped but smaller mu

# simulate deterministic HBD model
simulation = simulate_deterministic_hbd(LTT0       = Ntips, 
                                        oldest_age = oldest_age, 
                                        rho0       = 0.5,
                                        age_grid   = age_grid,
                                        lambda     = lambda,
                                        mu         = mu)

# plot deterministic LTT
plot( x = simulation$ages, y = simulation$LTT, type='l',
      main='dLTT', xlab='age', ylab='lineages')

castor documentation built on June 29, 2024, 9:08 a.m.