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.
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)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.
Strictly positive numeric, specifying the oldest time before present (``age'') to include in the simulation.
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.
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.
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.
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.
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.
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.
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.
Non-negative numeric, specifying the speciation rate (in units 1/time) at age0. Either lambda0 or lambda must be provided, but not both.
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.
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.
A named list with the following elements:
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.
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).
Numerical vector of size NG, listing the predicted (deterministic) total diversity (number of extant species, denoted \(N\)) at the ages given in ages[].
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.
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.
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.
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.
Numeric vector of size NG, listing the speciation rate (in units 1/time) at the ages given in ages[].
Numeric vector of size NG, listing the extinction rate (in units 1/time) at the ages given in ages[].
Numeric vector of size NG, listing the net diversification rate (\(\lambda-\mu\)) at the ages given in ages[].
Numeric vector of size NG, listing the pulled diversification rate (PDR, in units 1/time) at the ages given in ages[].
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\).
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\).
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).
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)\).
Non-negative numeric, specifying the product of the sampling fraction and the speciation rate at age0, \(\rho\cdot\lambda(\tau_o)\).
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.
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.
# NOT RUN {
# 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')
# }
Run the code above in your browser using DataLab