Simulating trajectories and phylogenies from a multi-type birth-death model

We sometimes have multiple demes, i.e. different types of individuals that contribute to the phylogeny or that can be sampled (e.g. juveniles vs. adults or acute vs. chronic infections).

We illustrate this example using a birth-death (BD) model with two patches (labeled 1 and 2) and migration between these patches (at a rate μ).

Initializing the system

\begin{aligned}
\frac{dI_1}{dt} = \beta_1 I_1 - \gamma_1 I_1 - \mu_1 I_1 + \mu_2 I_2 \\
\frac{dI_2}{dt} = \beta_2 I_2 - \gamma_2 I_2 - \mu_2 I_2 + \mu_1 I_1 \\
\end{aligned}

The associated reactions are:

reactions <- c("0 [beta1 * I1] -> I1",
               "I1 [gamma1 * I1] -> 0",
               "I1 [mu1 * I1] -> I2",
               "0 [beta2 * I2] -> I2",
               "I2 [gamma2 * I2] -> 0",
               "I2 [mu2 * I2] -> I1")

We then build the simulator:

bd_simu <- build_simulator(reactions)

The initial state variables values are

initialStates <- c(I1 = 0, I2 = 1)

The time range of the simulation is between 1975 and 2018:

time <- c(1975, 1998, 2018)

the parameters values are

theta <- list(gamma1 = c(0.2, 0.09), gamma2 = 0.1, mu1 = 0.025, mu2 = 0.021, beta1 = c(0.26,0.37), beta2 = 0.414)

and the time step (for the τ-leap and mixed algorithms) is:

dT <- 0.01

A safe version of the simulator bd_simu() is:

safe_bd_simu <- function(...) safe_run(bd_simu, ...)

Tau-leap trajectory simulation

We perform the simulations using:

trajbd_tl <- safe_bd_simu(
  paramValues = theta,
  initialStates = initialStates,
  times = time,
  method = "approximate",
  tau = 0.001)

We obtain:

head(trajbd_tl$traj)
#>       Time              Reaction Nrep I1 I2
#> 1 1975.000                  init    1  0  1
#> 2 1975.249  0 [beta2 * I2] -> I2    1  0  2
#> 3 1977.053  0 [beta2 * I2] -> I2    1  0  3
#> 4 1978.664 I2 [gamma2 * I2] -> 0    1  0  2
#> 5 1979.283 I2 [gamma2 * I2] -> 0    1  0  1
#> 6 1981.790  0 [beta2 * I2] -> I2    1  0  2

Graphically, we get:

plot(trajbd_tl)

Phylogeny simulation

With known sampling dates and known proportion of sampling

Instead of loading a vector, we assume we have 100 samples at 100 sampling dates between 2015 and 2018. We can generate the dates vector as:

dates_bd <- seq(from=2015, to=2018, length.out=100)

We then simulate a phylogeny where 20% of the sampling dates correspond to the I1 compartment, and 80% to the I2 compartment:

bd_tree <- simulate_tree(
  simuResults = trajbd_tl,
  dates = dates_bd,
  deme = c("I1", "I2"),
  sampled = c(I1 = 0.2, I2 = 0.8), # the type of individuals that are sampled and their proportion of sampling
  root = "I2", # type of individual at the root of the tree
  isFullTrajectory = FALSE, # deads do not generate leaves
  nTrials = 3,
  addInfos = TRUE) # additional info for each node

This is done using a coaslescence process informed by the trajectory. Therefore, each internal node of the phylogeny corresponds to a coalescence event and is associated with a label stoed in $node.label.

In our two-patches example, there are two types of coalesence: I2 individuals creating a new I2 individual, and I1 individuals creating a new I1 individual.

We can plot the phylogeny and color the internal nodes based on the type of coalescence.

First we generate a vector of colors for the nodes (if we find I2 in the node label we color it in blue, otherwise in red):

inode_cols <- ifelse(grepl(x=bd_tree$node.label,pattern="I2"),"blue","red")

Then we plot the phylogeny:

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, align.tip.label = T)
nodelabels(pch=20,col=inode_cols)

With known sampling dates, each assigned to a compartment by the user

One can give as input, sampling dates assigned to a compartment, in which case the option sampled is not required.

dates_bd <- seq(from=2015, to=2018, length.out=100)
dates_bd <- data.frame(Date=sample(dates_bd),Comp=c(rep("I1",20),rep("I2",80)))
head(dates_bd)
#>       Date Comp
#> 1 2017.212   I1
#> 2 2015.515   I1
#> 3 2017.606   I1
#> 4 2015.576   I1
#> 5 2017.515   I1
#> 6 2015.030   I1

Now let's simulate a phylogeny with sampling dates assigned to a compartment by the user.

bd_tree <- simulate_tree(
  simuResults = trajbd_tl,
  dates = dates_bd,
  deme = c("I1", "I2"),
  root = "I2", # type of individual at the root of the tree
  nTrials = 3,
  addInfos = TRUE) # additional info for each node

We can plot the phylogeny and color the external nodes given the compartment.

tips_cols <- ifelse(grepl(x=bd_tree$tip.label,pattern="I2"),"blue","red")

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)
tiplabels(pch=20,col=tips_cols)

With known sampling dates, each assigned to a compartment by the user

One can give as input, sampling dates assigned to a compartment, in which case the option sampled is not required.

dates_bd <- seq(from=2015, to=2018, length.out=100)
dates_bd <- data.frame(Date=sample(dates_bd),Comp=c(rep("I1",20),rep("I2",80)))
head(dates_bd)
#>       Date Comp
#> 1 2017.212   I1
#> 2 2015.515   I1
#> 3 2017.606   I1
#> 4 2015.576   I1
#> 5 2017.515   I1
#> 6 2015.030   I1

Now let's simulate a phylogeny with sampling dates assigned to a compartment by the user.

bd_tree <- simulate_tree(
  simuResults = trajbd_tl,
  dates = dates_bd,
  deme = c("I1", "I2"),
  root = "I2", # type of individual at the root of the tree
  nTrials = 3,
  addInfos = TRUE) # additional info for each node

We can plot the phylogeny and color the external nodes given the compartment.

tips_cols <- ifelse(grepl(x=bd_tree$tip.label,pattern="I2"),"blue","red")

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)
tiplabels(pch=20,col=tips_cols)

With only known sampling dates

In the case where the user has no information on the sampling proportions or the assignment of sampling dates on any compartment, the algorithm will randomly assign each sampling date to a compartment. The user gives as input sampling dates:

dates_bd <- seq(from=2015, to=2018, length.out=100)

Now let's simulate a phylogeny with sampling dates and no information about the sampling schemes :

bd_tree <- simulate_tree(
  simuResults = trajbd_tl,
  dates = dates_bd,
  sampled = c("I1","I2"),
  deme = c("I1", "I2"),
  root = "I2", # type of individual at the root of the tree
  nTrials = 10,
  addInfos = TRUE) # additional info for each node

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)

Without given sampling dates

Let's simuate a phylogeny where simulated deaths in the trajectory generate leaves. This can be done with the isFullTrajectory option.

bd_tree <- simulate_tree(
  simuResults = trajbd_tl,
  deme = c("I1", "I2"),
  root = "I2", # type of individual at the root of the tree
  nTrials = 10,
  isFullTrajectory = TRUE, # deads generate leaves
  addInfos = TRUE) # additional info for each node

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)

Comments

Please register or sign in to add a comment.

Tutorial 3

Simulating trajectories and phylogenies from a multi-type birth-death model

Initializing the system

Tau-leap trajectory simulation

Phylogeny simulation

With known sampling dates and known proportion of sampling

With known sampling dates, each assigned to a compartment by the user

With known sampling dates, each assigned to a compartment by the user

With only known sampling dates

Without given sampling dates

Comments