For this tutorial, we'll be using the example influenza data which is installed along with BDMM-Prime. This data, assembled from publicly-available data downloaded from NCBI GenBank, is the same set used in this MultiTypeTree tutorial. It consists of an aligned set of 60 H3N2 HA sequences sampled from New Zealand and Hong Kong. To make finding this data easy, we first select the package working directory.

Select "Set working dir" from the File menu, then choose "BDMM-Prime" from the sub-menu.

To load the data, simply choose "Import alignment" from the File menu and select one or more FASTA files.

Select "Import alignment" from the File menu to display a file selection dialog box. Using this, open the examples subdirectory and select the file h3n2_2deme.fna.

The BEAUti window should now look something like the following:

The sequences in this data set are sampled at different times, so use the following instructions to import this information into the analysis. (This step is unnecessary if sequences were collected at the same time, or near enough given the anticipated scale of the tree.)

Select the Tip Dates panel and check the "Use tip dates" option. Then press the "Auto-configure" button to set the times based on values extracted from the sequence headers. Finally, use the radio button and dropdown menu to interpret everything after the last "_" character as a numerical tip time.

For this example, we will use the basic HKY substitution model, with equilibrium nucleotide frequencies fixed to the empirical frequencies of characters in our alignment.

Select the "Site model" panel. Select the HKY model with "empirical" nucleotide frequencies.

This choice to fix the equilibrium frequencies is made here only to reduce the computational complexity of the tutorial analysis. For production analyses, these should probably be estimated, and Gamma-distributed site-to-site rate heterogeneity should also be used.

Since we have serially-sampled data, our analysis requires some kind of molecular clock. For the sake of simplicity, we use a strict clock, setting the initial value of the clock rate to $5\times 10^{-3}$ substitutions per site per year, which is close to the true value for influenza.

Select the "Clock model" panel. Set the "mean clock rate" value to 0.005.

The tree prior (also known as the phylodynamic likelihood) is the component of the analysis which BDMM-Prime provides. Here we configure the particular birth-death model we will use to relate the tree to the various population-level parameters we'd like to learn about.

Select the "Priors" panel. To choose the BDMM-Prime tree prior, find the drop-down menu next to Tree.t:h3n2_2deme and select "BDMMPrime".

The next thing we want to do is to specify the locations associated with each of the samples included in our analysis.

Scroll down to end of the tree prior section and find the table associating tree leaves with locations. Notice that the sample names have the form ID_Location_Date. To extract the locations from the sample names, click the Auto-configure button, select split on character, ensure "_" is specified as the delimiter, and select "2" from the take group(s) drop-down menu. Once this is complete, press "Ok".

The table should now be populated with types extracted from the sample names.

Now we get to the main part of the analysis setup, where we decide exactly how to model the generation of our data. In our case, we want to model a small part of the global transmission dynamics of H3N2.

The first thing to do is to is to select our desired parameterization of the birth-death model. In our case, as we're analysing pathogen sequences in an attempt to reconstruct transmission dynamics, we'll use the epidemiological parameterization.

Expand the tree prior section by clicking the arrow next to Tree.t:h3n2 on the left-hand side of of the screen. Then select "Epi Parameterization" from the drop-down menu at the top of the expanded section.

With this done, the next task is to define the dimensionality and initial values of the multi-type birth-death parameters. In BDMM-Prime, each of these parameters is considered a "skyline parameter" which can change in a piecewise-constant fashion through time.

The first skyline parameter to consider is the Re parameter, representing the effective reproductive number. By default parameter values are assumed to be the same across all types,

For the Re skyline parameter, uncheck the "Scalar values" check box to allow this parameter to take type-dependent values.

The next parameter is the "become uninfectious rate" parameter, representing the rate at which infected individuals become uninfectious. As this parameter is predominantly a function of the pathogen rather than location, we will assume it takes a single value which is shared among types. Furthermore, as influenza infections often take a relatively short time to pass, we will initialise this value to correspond to an infectious period of 1 week.

Find the Become Uninfectious Rate parameter, double-click the starting value and change it to 52 (per year). Be sure to press "Enter" to cause BEAUti to accept the new value.

Now let's consider the Sampling Proportion parameter. This parameter requires special care, as the sampling assumptions of the model can dramatically influence inference results. In our case we want (a) to allow for type-dependent sampling proportions and (b) to fix the sampling proportion to zero for all times earlier than the first sample.

In BDMM-Prime, this configuration is reasonably simple to configure.

Find the Sampling Proportion parameter, then follow these steps:

1. Uncheck the "Scalar value" checkbox to allow for type-dependent values.
2. Check "Display visualization" to see the distribution of sample times relative to the most recent sample.
3. Set the "Number of change times" value to 1. (Notice the appearance of the epoch boundary marker on the visualisation.)
4. Set the change time for the boundary between epoch 1 and 2 to 5.7. (Notice that this value ensures the oldest sample remains in epoch 1.)
5. In the "Values" table, double-click the entries corresponding to the sampling proportion values for Epoch 2 and change each to 0.0, remembering to press "enter" after each change.

Note skyline parameter values of 0 are always fixed in the analysis, regardless of whether or not the "estimate values" option is checked.

Finally, we need to define how lineages change type in our multi-type model. Like BDMM, BDMM-Prime allows both direct "migration" as well as "birth to a new type" styles of multi-type models. Given we are modelling the movement of infected influenza hosts, migration suits our situation best. We also acknowledge that migration rates between different pairs of types/locations may be different from one another.

To incorporate these decisions, follow these instructions:

Find the Migration Rate parameter, then:

1. Double-click on the value and change it from it's default to the smaller rate of 0.1 (any infected individual has a 10% chance of moving between the two locations in any given year).
2. Uncheck the "Scalar values" checkbox to allow rates to vary between pairs of locations.
3. Check the "Estimate values" checkbox to ensure these migration rates are estimated.

With the multi-type model is configured, what remains is to define priors for each of the parameters we're estimating. As in any ther BEAST analysis, this includes parameters each of the various site, clock and phylodynamic models we've configured.

Rather than describe in detail how to set these priors, here we will simply list a set of simple priors which can be used for the purpose of demonstrating the model. (Detailed information on setting up priors in BEAUti can be found elsewhere, see for instance the Prior Selection tutorial on the Taming the BEAST website.

Set the following parameter priors:

Parameter	Prior
ReEpi.t:h3n2_2deme	LogNormal(0,0.5)
becomeUninfectiousRateEpi.t:h3n2_2deme	1/X
migrationRateEpi.t:h3n2_2deme	Exp(0.1)
samplingProportionEpi.t:h3n2_2deme	Unif(0,1)

Finally, we need to configure the details of the MCMC chain to run and the trees and parameters to be logged.

Switch to the MCMC panel and change the Chain Length to 1000000 ($10^6$).

For this demonstration we will leave everything else as-is. Notice however that beyond the regular trace, screen and tree logs, there are three additional loggers. We'll discuss these in more detail

typedTreeLogger

A logger which uses stochastic mapping to generate samples from the marginal posterior over edge-typed trees, which amount to hypotheses about when and where each type change occurred on the tree.

nodeTypedTreeLogger

Similar to typedTreeLogger, but the sampled trees includes type information only at internal nodes. Used to build summary trees with BEAST 2's TreeAnnotator utility.

trajLogger

A logger which uses particle filtering to log samples from the marginal posterior of multi-type population trajectories.

Each of these loggers is optional, and can be enabled/disabled by clicking the arrow buttons to their left and checking/unchecking the "Enable logger" box. By default edge-typed tree logging is enabled, but trajectory logging is disabled.

Firstly, let's open the parameter log file using the Tracer app (distributed with BEAST 2).

Open Tracer, select "Import trace file…" from the File menu. Use the interface to explore the marginal distributions as estimated from the MCMC samples.

We see immediately that the ESS for many of the parameters is low, which is unsurprising for such a short run. However, the focus here is on understanding how BDMM-Prime analysis results are presented in the log file rather than trying to seriously interpret the results.

With this in mind, see firstly that each of the estimated parameters including Re, becomeUninfectiousRate and migrationRate are represented in the log file by ReSPEpi, becomeUninfectiousRateSPEpi and migrationRateSPEpi. The "SP" indicates that the parameter is a "Skyline parameter" (as opposed to BEAST 2's regular RealParameter).

Also notice that when we have allowed type-specific (not scalar) values, as in the case of Re, the parameter name appears multiple times with a suffix indicating the type associated with the specific value.

Finally, notice that when we have allowed the parameter to change at a particular time, each epoch/interval receives its own set of parameter values. In this case an additional ".endtime" value is also recorded, which identifies the end time of a specific interval, measured relative to the start time of the process.

The untyped tree log file can be dealt with in the same way as all BEAST 2 tree files. However we also have tree log files containing sampled ancestral type information which we can use to learn about the type change process as it occurred on the tree. In this influenza case, this contains potentially interesting information about global influenza spatial dynamics.

While in principle the .typed.trees file contains the most detailed information, it is easier to jointly summarize these ancestral dynamics with the tree itself by applying the TreeAnnotator program to the .typed.node.trees file.

1. Open TreeAnnotator (distributed with BEAST 2)
2. For the node heights option, select "mean heights". (This produces less biased estimates of divergence times than the default.)
3. Select h3n2.h3n2_2deme.typed.node.trees as the input file.
4. Choose an appropriate output file name (e.g. summary.tree)
5. Click "Run".

Once complete, you can visualise the summary tree using any BEAST-compatible tree viewer, such as FigTree or IcyTree.

1. Navigate to https://icytree.org in a javascript-enabled web browser.
2. Select "Load from file" from the File menu and select the summary.tree file produced by TreeAnnotator.
3. From the Style menu choose "colour nodes by" and then "type".
4. Also from the Style menu choose "Node height error bars" and then "height_95%_HPD".

This should produce something quite like the figure below:

Notice that we've chosen to colour the nodes themselves rather than the edges. This is done precisely because in the case of these node-typed symmary trees, it is only the types at the nodes that are being represented. Please resist the urge to colour the edges based on these values! While tree viewers make this easy, and it is frequently done, it is a misleading representation.

The "canonical" parameterization is in a sense the native mathematical description of the birth-death model; the its parameters reflect the rates of the basic reactions which make up the process.

Birth Rate ($\lambda_{ii}$)

This is the (per individual per unit time) rate at which each individual in the population $i$ produce new individuals of their own type.

Death Rate ($\mu_i$)

Similarly, this is the rate at which individuals in population $i$ are removed from the population.

Sampling Rate ($\psi_i$)

This is the rate at which individuals in population $i$ are "sampled". A non-zero rate here corresponds to a process that will produce samples that are distributed through time. (Sometimes referred to as $ψ$-sampling.)

Removal probability ($r_i$)

This is the probability with which sampled individuals in population $i$ are also removed from the population.

Rho sampling ($\rho_i$)

This parameterizes binomial sampling at a particular time (usually the present). If enabled, each individual in population $i$ is considered in turn, and sampled with this probability.

In addition, the following rates parameterize reactions which allow type changes to occur. In general these are specified as matrices, allowing for a unique rate for every ordered pair of distinct types.

Migration Rate ($m_{ij}$)

This is the rate at which an individual with some type becomes an individual of another type.

Birth Among Demes ($\lambda_{ij}$)

This is the rate at which an individual with some type produces an individual of another type.

Finally, the following parameter expresses the total duration of the birth-death-sampling process.

ProcessLength

the time between the start of the process and the most recent sample. (Here equivalent to the "origin" of the process.)

The assumption behind the "epidemiological" parameterization is that individuals in the model represent infected individuals in a population. In place of the canonical birth, death and sampling parameters we have

Re ($R_{ii}$)

This is the type-specific effective reproductive number, defined by $R_i=\lambda_{ii}/(\mu_i+r_i\psi_i)$. Values above 1 correspond to growing (in expectation) numbers of infected individuals.

Become uninfectious rate ($\delta_i$)

This is the rate at which infected individuals are removed from the infectious pool, either with or without sampling. In terms of the canonical parameters, this is $\mu_i+r_i\psi_i$.

Sampling proportion ($s_i$)

This is the probability with which infected individuals are sampled through time, defined as $\psi_i/(\mu_i+r_i\psi_i)$.

Besides the migration matrix, epidemiological compartment changes are parameterized by:

Re Among Demes ($R_{ij}$)

This is a generalisation of the effective reproductive number to account for transmissions between epidemiological compartments. It is defined as $R_{ij}=\lambda_{ij}/(\mu_i+r_i\psi_i)$.

(The rho sampling, removal probability, migration and process length parameters are identical to those described for the canonical parameterization.)

The FBD parameterization is most appropriate when individuals in the model are intended to represent whole species, and birth and death events correspond to speciation and extinction. In place of the canonical birth, death and sampling rate parameters we have

Diversification ($d_{ii}=\lambda_{ii}-\mu_i$)

This is the rate at which species richness increases through time.

Turnover ($t_{i}=\mu_i/\lambda_{ii}$)

For $d_{ii}>0$, this parameter indicates the relative frequency with which species go extinct before speciation, and satisfies $0\leq t_i < 1$.

Sampling proportion ($s_i=\psi_i/(\mu_i+r_i\psi_i)$)

Analogous to the epdemiological parameterization, this indicates the probability with which species leave fossils before going extinct.

Besides the migration rate, species type changes are parameterized by:

Diversification between types ($d_{ij}=\lambda_{ij}-\mu_i$)

This is a generalisation of diversification rate, but allowing for one of the daughter species to carry a new type.

(The rho sampling, removal probability, migration and process length parameters are identical to those described for the canonical parameterization.)

Note: Unlike the Canonical and Epidemiological parameterizations, the FBD parameterization on its own does not permit exploration of the whole parameter space. In particular, it makes the assumption that diversification is always increasing ($d_{ii}>0$ and $t_{i}<1$). If you want to also allow negative diversification, use the Canonical parameterization instead.

In this section we delve into the details of the BDMM-Prime skyline parameter user interface. Although focus here on vector skyline parameters, everything here is also applicable to vector skyline parameters.

The input editor BDMM-Prime uses for specifying skyline parameters is quite complex, reflecting the large number of ways in which these parameters can be configured. The general layout is as follows:

Here you can see that the configuration editor consists of three main parts: the epoch configuration, the value configuration and an epoch visualiser for helping with the configuration. The following sections detail each of these parts in turn.

For our example here, we use the trajectory file generated by the analysis from the tutorial section with trajectory logging enabled as above.

Trajectory data generated by BDMM-Prime is stored as a simple table in the plain text file with the extension ".traj". The format of this file is described in detail in the appendix.

The easiest way to plot trajectories sampled using BDMM-Prime is using the tidyverse suite of R packages. To install these in R (if you haven't already), simply issue the following command from the R command line:

install.packages("tidyverse")

(You'll only have to do that once.) With the packages installed you'll need make them available in your current R session by issuing the command:

library(tidyverse)

The first step is to load the trajectory file using the read_tsv command.

traj <- read_tsv("h3n2.h3n2_2deme.traj")

When doing this it is a good idea to omit some of the early trajectories to account for burn-in. (Examining the parameter log file is one way to get an idea of how much to remove.) As described in the trajectory file format section, the table has several columns, one of which is "Sample", which indicates the MCMC iteration a particular trajectory belongs. We can use this to remove trajectories

For example, the following removes trajectories belonging to the first 10% of samples, assuming $10^6$ samples overall:

traj <- traj %>% filter(Sample>=100000)

Another column in the table is named "variable", and is used to indicate the type of information a given row provides about the trajectory. To plot population sizes, we are only interested in rows where the variable has the value "N", which means that the row contains information about a population size. Thus, we can further filter our table to include only these rows:

traj <- traj %>% filter(variable=="N")

To actually plot the distribution, we use the fact that each row now contains an "age" relative to the end of the process, a "type" indicating the type of the population, and a "value" indicating the size of the population immediately before age.

ggplot(traj, aes(age, value, col=factor(type), group_by=factor(Sample))) +
  geom_line(alpha=0.2) + xlim(0,10) +
  ylab("Population size")

Here we've chosen to simply plot all of the sampled trajectories together, with each trajectory coloured by its type. Note the group_by option, which is used to tell ggplot that rows with the same "Sample" value belong to the same trajectory. The use of factor() is important, as without this ggplot will treat the type and Sample indices as continuous variables rather than integers. The xlim(0,10) command is included to focus on the time domain covered by the majority of the trajectories.

Notice that the types are indicated using an integer rather than the original label. As described in the section on start type prior probabilities, these integers are simply the indices of the type names when sorted lexicographically (generalized alphabetically). Thus in this case, 0 stands for Hong Kong and 1 stands for New Zealand.

While plotting all of the trajectories in this way can be interesting, it's often more useful to display precise credible intervals for the population sizes through time, similar to the approach often used to visualize the results of Bayesian skyline plot analyses.

To do this, we first define a vector of ages which will span the time domain of our plot. In our case, we'll use 1001 values spread evenly between ages 0 and 10:

ages <- seq(0,10,length.out=1001)

Then, we need to determine the population size in each trajectory at every one of these ages. We can do this by making use of the approx() function, which is a base R procedure for interpolating function sampled at discrete time points.

gridded_traj <- traj %>%
  group_by(Sample, type) %>%
  reframe(N=approx(age, value, ages, method="constant", f=0, yright=0)$y,
          age=ages)

We can then summarize the distribution at each time point by computing the 50%, 2.5%, and 97.5% quantiles, corresponding respectively to the median and the lower and upper bounds of the 95% central posterior density (CPD).

gridded_traj <- gridded_traj %>%
  group_by(type,age) %>%
  summarize(low=quantile(N,0.025),
            high=quantile(N,0.975),
            med=quantile(N,0.5))

These values can be nicely plotted using geom_ribbon() to display the CPD boundaries and a geom_line() to display the medians.

ggplot(gridded_traj,
       aes(age, N, col=factor(type), fill=factor(type), y=med, ymin=low, ymax=high)) +
    geom_ribbon(alpha=0.5) +
  geom_line() +
  ylab("Population size")

This yields the following graphic:

This can be further improved by including actual calendar times on the time axis. To do this we use our knowledge that the final included sample is from September 1, 2005, meaning that this is the date relative to which our ages are computed. (Our ages are the number of years prior to this date.)

ggplot(gridded_traj %>%
       mutate(Date=ymd("2005/09/01")-age*365) %>%
       mutate(Location=recode(factor(type),
                         "0"="Hong Kong",
                         "1"="New Zealand")),
       aes(date, N, col=Location, fill=Location, y=med, ymin=low, ymax=high)) +
    geom_ribbon(alpha=0.5) +
    geom_line() + ylab("Population size") +
    scale_x_date(date_breaks="1 year", date_labels="%Y")

As an aside, you may have noticed that the population sizes inferred in these example analyses are extremely small, given that these are intended to reflect global Influenza H3N2 dynamics with sampled sequences drawn from different countries across several years. Recall however that both $R_e$ and the become uninfectious rate $\delta$ were assumed constant in the tutorial. By fitting such a simple constant-rate birth-death model to this dataset, our parameters are forced to represent "effective" rates and "effective" population sizes which best fit the birth-death dynamics on this scale.

While these effects are made very obvious by the population trajectory posteriors, they are no less present when performing birth-death analyses without trajectory sampling. It is therefore always important to keep the possibility of these kinds of model-fitting effects in mind when interpreting birth-death results and phylodynamic inference results in general, even when applied at smaller scales.

<log spec="bdmmprime.util.StartTypePosteriorProbsLogger">
  <!-- Inputs -->
</log>

This logger logs the posterior distribution of the start type conditional on the current parameter values. This posterior is computed by the BirthDeathMigrationDistribution while evaluating the tree prior. When averaged over the entire MCMC chain, these probabilities yield the complete marginal posterior distribution over the starting state conditional on the data.

This object takes the following input:

bdmmTreePrior (BirthDeathMigrationDistribution)

Reference to the BirthDeathMigrationDistribution used in the analysis.

<log spec="bdmmprime.mapping.TypeMappedTree">
  <!-- Inputs -->
</log>

Stochastic mapping is handled in BDMM-Prime by a TypeMappedTree object, which is intended to be logged. The stochastic mapping calculation is performed on the input state only when the output has to be appended to the log file.

The important inputs are:

parameterization (Parameterization)

bdmmDistrib (BirthDeathMigrationDistribution)

One of these two inputs must be provided, as these provide the parameters needed to sample the ancestral type distribution.

untypedTree (Tree, required)

The tree on which the type changes are to be stochastically mapped.

typeTraitSet (TraitSet, required)

Trait set specifying the taxon types.

startTypePriorProbs (RealParameter, required)

The prior probabilities for the initial individual type.

typeLabel (String, required)

The label to be used for traits in the generated metadata. In other words, this is the key used in the key-value pairs which will be used to annotate the tree with the ancestral type information.

finalSampleOffset (Function, default 0)

If provided, the difference in time between the final sample and the end of the BD process.

Because the TypeMappedTree class extends the Tree class it accepts many inputs that are not particularly relevant. These have been omitted.

<logger fileName="output.traj" logEvery="10000">
  <log spec="bdmmprime.trajectories.SampledTrajctory">
  </log>
</logger>

Trajectory sampling is handled in BDMM-Prime by a SampledTrajectory object, which is intended to be logged. This implements a particle filtering calculation which is performed only when the output has to be appended to the log file.

Important: The SampledTrajectory must be the only input passed to its parent logger. Broken output files will result if this rule is violated!

Objects of type SampledTrajectory accept the following inputs:

parameterization (Parameterization)

bdmmDistrib (BirthDeathMigrationDistribution)

One of these two inputs must be provided, as these provide the parameters needed to sample the population dynamics.

typeMappedTree (Tree, required)

An edge-typed tree. Ideally a reference to a TypeMappedTree.

typeLabel (String, required)

The type label which was used to store type information on the tree.

finalSampleOffset (Function, default 0)

If provided, the difference in time between the final sample and the end of the BD process.

nParticles (Integer, default 1000)

Size of particle ensemble used in particle filtering. (Larger values give more accurate results, but make the calculation more involved.)

resampThresh (Double, default 0.5)

Relative effective size of the weighted particle ensemble at which to trigger a resampling step.

useTauLeaping (Boolean, default false)

If true, use the faster but approximate tau-leaping algorithm to simulate population histories. Otherwise, use the slower but statistically correct Gillespie algorithm.

minLeapCount (Integer, default 100)

When tau leaping is enabled, at least this many steps are used to simulate each population trajectory.

epsilon (Double, default 0.03)

Controls the adaptive time step used in the tau leaping algorithm, with lower values causing the step size propsals to be more conservative (smaller steps). For details, refer to Cao et al., 2006.