simulate_sbm {castor}

Simulate Spherical Brownian Motion on a tree.

Description

Given a rooted phylogenetic tree and a Spherical Brownian Motion (SBM) model for the evolution of the geographical location of a lineage on a sphere, simulate random outcomes of the model on all nodes and/or tips of the tree. The function traverses nodes from root to tips and randomly assigns a geographical location to each node or tip based on its parent's previously assigned location and the specified model parameters. The generated states have joint distributions consistent with the SBM model (Perrin 1928; Brillinger 2012). This function generalizes the simple SBM model to support time-dependent diffusivities.

Usage

simulate_sbm(tree, 
             radius, 
             diffusivity,
             time_grid      = NULL,
             splines_degree = 1,
             root_latitude  = NULL, 
             root_longitude = NULL)

Arguments

tree

A rooted tree of class "phylo". The root is assumed to be the unique node with no incoming edge. Edge lengths are assumed to represent time intervals or a similarly interpretable phylogenetic distance.

radius

Strictly positive numeric, specifying the radius of the sphere. For Earth, the mean radius is 6371 km.

diffusivity

Either a single numeric, or a numeric vector of length equal to that of time_grid. Diffusivity ("D") of the SBM model (in units distance^2/time). If time_grid is NULL, then diffusivity should be a single number specifying the time-independent diffusivity. Otherwise diffusivity specifies the diffusivity at each time point listed in time_grid.

Under a planar approximation the squared geographical distance of a node from the root will have expectation 4LD, where L is the node's phylogenetic distance from the root. Note that distance is measured in the same units as the radius (e.g., km if the radius is given in km), and time is measured in the same units as the tree's edge lengths (e.g., Myr if edge lengths are given in Myr).

time_grid

Numeric vector of the same length as diffusivity and listing times since the root in ascending order, or NULL. This can be used to specify a time-variable diffusivity (see details below). If NULL, the diffusivity is assumed to be constant over time and equal to diffusivity (which should be a single numeric). Time is measured in the same units as edge lengths, with root having time 0.

splines_degree

Integer, either 0,1,2 or 3, specifying the polynomial degree of the provided diffusivity between grid points in time_grid. For example, if splines_degree==1, then the provided diffusivity is interpreted as a piecewise-linear curve; if splines_degree==2 it is interpreted as a quadratic spline; if splines_degree==3 it is interpreted as a cubic spline. 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.

root_latitude

The latitude of the tree's root, in decimal degrees, between -90 and 90. If NULL, the root latitude is chosen randomly according to the stationary probability distribution of the SBM.

root_longitude

The longitude of the tree's root, in decimal degrees, between -180 and 180. If NULL, the root longitude is chosen randomly according to the stationary probability distribution of the SBM.

Details

For short expected transition distances this function uses the approximation formula by Ghosh et al. (2012). For longer expected transition distances the function uses a truncated approximation of the series representation of SBM transition densities (Perrin 1928).

The pair time_grid and diffusivity can be used to define a time-dependent diffusivity, with time counted from the root to the tips (i.e. root has time 0) in the same units as edge lengths. For example, to define a diffusivity that varies linearly with time, you only need to specify the diffusivity at two time points (one at 0, and one at the time of the youngest tip), i.e. time_grid and diffusivity would each have length 2. Note that time_grid should cover the full time range of the tree; otherwise, diffusivity will be extrapolated as a constant when needed.

If tree$edge.length is missing, each edge in the tree is assumed to have length 1. The tree may include multifurcations as well as monofurcations.

Value

A list with the following elements:

success

Logical, specifying whether the simulation was successful. If FALSE, then an additional return variable error will contain a brief description of the error that occurred, and all other return variables may be undefined.

tip_latitudes

Numeric vector of length Ntips, listing simulated decimal latitudes for each tip in the tree.

tip_longitudes

Numeric vector of length Ntips, listing simulated decimal longitudes for each tip in the tree.

node_latitudes

Numeric vector of length Nnodes, listing simulated decimal latitudes for each internal node in the tree.

node_longitudes

Numeric vector of length Nnodes, listing simulated decimal longitudes for each internal node in the tree.

Author(s)

Stilianos Louca

References

F. Perrin (1928). Etude mathematique du mouvement Brownien de rotation. 45:1-51.

D. R. Brillinger (2012). A particle migrating randomly on a sphere. in Selected Works of David Brillinger. Springer.

A. Ghosh, J. Samuel, S. Sinha (2012). A Gaussian for diffusion on the sphere. Europhysics Letters. 98:30003.

S. Louca (2021). Phylogeographic estimation and simulation of global diffusive dispersal. Systematic Biology. 70:340-359.

See Also

simulate_ou_model, simulate_rou_model, simulate_bm_model, fit_sbm_const

Examples

## Not run: 
# generate a random tree
tree = generate_random_tree(list(birth_rate_intercept=1),max_tips=100)$tree

# simulate SBM on the tree
simulation = simulate_sbm(tree, radius=6371, diffusivity=1e4,
                          root_latitude=0, root_longitude=0)

# plot latitudes and longitudes of the tips
plot(simulation$tip_latitudes,simulation$tip_longitudes)

## End(Not run)