A population in ReMASTER is simply a RealParameter (or other Function) with a specific ID. The value of the RealParameter is used as its initial population size in the simulation.

For example, the following defines a population named "X" with an initial size of 10:

<population id="X" spec="RealParameter" value="10"/>

If the RealParameter is a vector, its elements represent subpopulations. For example, the following population "Y" has two subpopulations initially occupied by 1 and 10 individuals, respectively:

<population id="Y" spec="RealParameter" value="1 10"/>

Populations may also be marked as "sample" populations by using the appropriate tag:

<samplePopulation id="sample" spec="RealParameter" value="0"/>

For trajectories, sample populations behave just like any other population, besides being forbidden from appearing on the reactant side of a reaction. individuals of these populations take on a special meaning when trees are simulated, however, as the production of a sample corresponds to the creation of a node (internal or external) in the resulting tree.

Reactions connect zero or more reactant individuals with zero or more product individuals. The dynamics that ReMASTER simulates are simply the result of applying reactions at random or pre-determined times.

There are two distinct kinds of reactions which can be specified: "continuous" reactions and "punctual" reactions.

By composing populations and reactions as defined above you can construct birth-death (and later coalescent) models of very large complexity. However, for models having many distinct types, explicitly describing every possible reaction between these types can become extremely tedious.

For instance, consider a model composed of 5 populations X,Y,Z,W, and V, in which every population has initially 1 individual and each is subject to a birth reaction.

To set this up directly would involve something like the following:

<population id="X" spec="RealParameter" value="1"/>
<population id="Y" spec="RealParameter" value="1"/>
<population id="Z" spec="RealParameter" value="1"/>
<population id="W" spec="RealParameter" value="1"/>
<population id="V" spec="RealParameter" value="1"/>

<reaction spec="Reaction" rate="1.0"> X -> 2X </reaction>
<reaction spec="Reaction" rate="1.0"> Y -> 2Y </reaction>
<reaction spec="Reaction" rate="1.0"> Z -> 2Z </reaction>
<reaction spec="Reaction" rate="1.0"> W -> 2W </reaction>
<reaction spec="Reaction" rate="1.0"> V -> 2V </reaction>

If we also wanted to include reactions involving pairs of different populations, the model specification could easily run to many pages.

Happily, BEAST 2 provides a compact way to express repeated a repeated fragment of XML: the <plate> element. To use this, we surround a fragment of XML to repeat as follows:

<plate var="i" range="A,B,C,...">
  <!-- fragment to repeat -->
</plate>

This causes the body of this element to be repeated a certain number of times, defined by the number of comma-delimited values in the string passed to the range input. The var input defines the name of a variable which will be substituted for each of these values when it appears in the fragment.

Using <plate> notation, our model specification can be written much more simply as

<plate var="pop" range="X,Y,Z,W,V">
  <population id="$(pop)" spec="RealParameter" value="1"/>
  <reaction spec="Reaction" rate="1.0"> $(pop) -> 2$(pop) </reaction>
</plate>

We can also nest <plate> elements to create reactions involving multiple populations. For instance, the following defines migration between each possible pair of populations:

<plate var="source" range="X,Y,Z,W,V">
  <plate var="dest" range="X,Y,Z,W,V">
    <reaction spec="Reaction" rate="0.5"> $(source) -> $(dest) </reaction>
  </plate>
</plate>

(Note that this will produce some reactions that do nothing, but this is not an issue.)

Besides a comma-delimited list of symbols, range can also take a sequence notation "start:end" to represent a sequence of integers between start and end inclusive. This can be helpful in creating reactions involving many subpopulations:

<population id="X" spec="RealParameter" dimension="100" value="1"/>
<plate var="i" range="0:99">
  <reaction spec="Reaction" rate="1.0"> X[$(i)] -> 2X[$(i)] </reaction>
</plate>

Once the population and reaction elements are specified, defining a stochastic trajectory in ReMASTER's XML is very easy.

For example, the following represents a stochastic birth-death-sampling process trajectory with constant rates:

<trajectory spec="StochasticTrajectory"
            maxTime="10">
  <population id="X" spec="RealParameter" value="1"/>

  <reaction spec="Reaction" rate="1.2"> X -> 2X </reaction>
  <reaction spec="Reaction" rate="1"> X -> 0 </reaction>
</trajectory>

The maxTime attribute specifies the maximum end time of the simulation. For stochastic trajectories this attribute is actually optional, since there are many birth-death processes (although not this one) are guaranteed to terminate in finite amount of time. (ReMASTER always ceases the simulation once no reactions are able to fire.)

Defining deterministic trajectories is almost exactly the same as defining a stochastic trajectory. The only difference is that the spec attribute must be set to "DeterministicTrajectory" and the maxTime attribute is mandatory.

For example, the following represents an exponentially growing population with rate $\lambda-\mu=1.2-1.0=0.2$:

<trajectory spec="DeterministicTrajectory"
            maxTime="10">
  <population id="X" spec="RealParameter" value="1"/>

  <reaction spec="Reaction" rate="1.2"> X -> 2X </reaction>
  <reaction spec="Reaction" rate="1"> X -> 0 </reaction>
</trajectory>

Another minor difference is that the endsWhen condition can only include inequalities, since detecting equality in floating point numbers can be problematic.

Finally, DeterministicTrajectory accepts two numerical attributes, forwardRelativeStepSize and backwardRelativeStepSize. These attributes define the size of the steps taken by the forward-time and backward-time numerical integrators, respectively, relative to the total simulation time. Reducing the size of these values increases the accuracy of the deterministic trajectory calculation and the associated tree simulation, but makes the calculation more computationally demanding. The forward-time step size affects the population size trajectory calculation accuracy, while the backward-time step size affects only the subsequent tree simulation accuracy. Leaving these values to their default settings is usually okay, but smaller step sizes may be necessary to avoid numerical errors in certain situations. (These sometimes manifest in the form of "too many lineages remaining" errors.) Models where populations are changing very quickly are particularly susceptible to such problems.

In addition to the rates and initial population sizes, conditions can be placed on trajectory outcomes so that the simulation finishes when the population sizes reach some specified target or, in the case of stochastic trajectories, the resulting trajectory is only "accepted" when a particular condition is met.

To define an early termination point, provide the endsWhen attribute to the trajectory element. The value of the attribute should be an expression predicating which will cause the simulation to end when it becomes true. For instance

<trajectory spec="StochasticTrajectory"
            endsWhen="X==100">
...
</trajectory>

defines a trajectory whose simulation will be terminated the instant that the $X$ population contains exactly 100 individuals.

Similarly, to define a condition on the finished simulation, use the mustHave attribute. For instance,

<trajectory spec="StochasticTrajectory"
            mustHave="sample>=50">
...
</trajectory>

specifies that the "sample" population must contain at least 50 individuals at the end of the simulation.

The complete syntax for these predicates is specified in section 8.1. Beware that XML (unfortunately) forbids use of characters such as < and & in attribute values, so to use them you'll need to use the escaped versions < and &.

While the above XML fragments are sufficient to represent a simulated trajectory, in order to actually be useful it needs to be embedded in a BEAST XML that produces output. The most direct way to do this is via ReMASTER's Simulator object, which can be used as a <run> element as follows:

<beast version="2.0" namespace="beast.base.inference:beast.base.inference.parameter:remaster">
  <run spec="Simulator" nSims="100">
    <simulate id="trajectory" spec="StochasticTrajectory" maxTime="10">
      <population id="X" spec="RealParameter" value="1"/>
      <reaction spec="Reaction" rate="1.2"> X -> 2X </reaction>
      <reaction spec="Reaction" rate="1"> X -> 0 </reaction>
    </simulate>

    <logger fileName="output.traj">
      <log idref="trajectory"/>
    </logger>
  </run>
</beast>

The Simulator takes the following inputs:

nSims (required)

an integer specifying the number of independent simulations to perform,

simulate (required)

an object representing the simulation, in this case a StochasticTrajectory, and

logger

what to record in the output files for each replicate of the simulation. This input can take multiple values, one for each output file. (In this case we only produce one output file, containing a representation of the whole simulated trajectory.)

To run the simulation and produce the output, copy the above XML into a file named something like "BDtraj.xml" and then run it through BEAST, either via the command line:

$ beast BDtraj.xml

or by running BEAST and selecting the "BDtraj.xml" file using the graphical file choose which appears on startup.

This should produce some screen output, together with a file "output.traj" in the same directory containing the simulated trajectories.

The "output.traj" file generated by the previous example is simply a tab-separated-variable (TSV) file containing the size of the population after each event in every replicate trajectory. This file be loaded directly into R for further analysis or for plotting.

From an R session in the same working directory as the simulation results, use the following command to load the trajectories into a data frame:

traj <- read.table("output.traj", header=T)

Alternatively, the read_tsv command from the readr package (Tidyverse) can be used to load the trajectories into a "tibble", which is often faster for large files:

library(readr)
traj <- read_tsv("output.traj")

The resulting data frame has the following layout:

> traj
# A tibble: 7,605 × 5
   Sample     t population index value
    <dbl> <dbl> <chr>      <dbl> <dbl>
 1      0 0     X              0     1
 2      0 0.117 X              0     0
 3      1 0     X              0     1
 4      1 2.60  X              0     2
 5      1 2.99  X              0     1
 6      1 3.17  X              0     0
 7      2 0     X              0     1
 8      2 0.583 X              0     2
 9      2 0.840 X              0     1
10      2 0.871 X              0     2
# ℹ 7,595 more rows

Each row specifies the size of a particular population (or subpopulation) at a particular time. The columns are as follows:

Sample

A number which indicates to which specific simulated trajectory the results belong.

t

The time with which the population size is associated.

population

The name (id) of the population.

index

An index starting at 0 indicating the subpopulation with which the size is associated. (This column will always be 0 if no subpopulations are defined.)

value

The size of the population.

This format is ideal for plotting using ggplot2:

library(ggplot2)
ggplot(traj, aes(t, value, group_by=factor(Sample))) + geom_step(alpha=0.5)

(The group_by=factor(Sample) ensures that population sizes associated with different trajectories are grouped appropriately.)

Replacing StochasticTrajectory with DeterministicTrajectory in the simulation XML (and setting nSims="1") produces an output file which, when plotted, displays the following figure depicting the expected population size dynamics:

(The deterministic model gives the true expected population size dynamics only for linear models.)

In addition to logging full trajectories, you can log basic statistics of simulated trajectories using the StochasticTrajectoryStatsLogger logger, which can be added using an XML fragment like this:

<logger fileName="traj_stats.log">
  <log spec="StochasticTrajectoryStatsLogger" traj="@traj"/>
</logger>

This will write things like the total event count and the final population states to columns of the output file. Although all of this information is contained in the full trajectory log, occasionally it is useful to log only a few basic statistics instead - the full trajectory logs can be very big.

(Additional statistics may be added to this logger in future.)

A self-contained tree simulation XML is very similar to a self-contained trajectory simulation XML. The only differences are that tree simulation XMLs have a SimulatedTree object for the <simulate> element, and generally have a logger which writes the simulated trees to a file.

For example, here is a full XML description of a tree simulator which produces trees sampled from a linear birth-death-sampling model and saves them to a file named "output.trees":

<beast version="2.0" namespace="beast.base.inference:beast.base.inference.parameter:remaster">
  <run spec="Simulator" nSims="1">
    <simulate id="tree" spec="SimulatedTree">
      <trajectory id="trajectory" spec="StochasticTrajectory" maxTime="10">
        <population id="X" spec="RealParameter" value="1"/>
        <samplePopulation id="samp" spec="RealParameter" value="0"/>

        <reaction spec="Reaction" rate="1.4"> X -> 2X </reaction>
        <reaction spec="Reaction" rate="0.5"> X -> 0 </reaction>
        <reaction spec="Reaction" rate="0.5"> X -> samp </reaction>
      </trajectory>
    </simulate>

    <logger fileName="output.traj">
      <log idref="trajectory"/>
    </logger>

    <logger fileName="output.trees">
      <log idref="tree"/>
    </logger>
  </run>
</beast>

Note that in this example we have left the trajectory logger in place, so in fact this XML produces both output.traj and output.tree files.

Loading the tree file in a phylogenetic tree viewer such as icytree, FigTree or into R using ggtree allows you to view the result:

The example model in the previous section contained only one main population type (besides the sample population). However, many interesting models involve multiple types and result in trees whose edges are associated with different populations. The standard BEAST tree logger used in the previous example will not record this information in its output file.

In order to include this information so that it can be shown on the tree, you'll need to use the TypedTreeLogger in place of the standard logger:

<beast version="2.0" namespace="beast.base.inference:beast.base.inference.parameter:remaster">
  <run spec="Simulator" nSims="1">
    <simulate id="tree" spec="SimulatedTree">
      <trajectory id="trajectory" spec="StochasticTrajectory" maxTime="10">
        <population id="X" spec="RealParameter" value="1"/>
        <population id="Y" spec="RealParameter" value="0"/>
        <samplePopulation id="samp" spec="RealParameter" value="0"/>

        <reaction spec="Reaction" rate="1.4"> X -> 2X </reaction>
        <reaction spec="Reaction" rate="0.1"> X -> Y </reaction>
        <reaction spec="Reaction" rate="0.5"> X -> 0 </reaction>
        <reaction spec="Reaction" rate="0.5"> Y -> samp </reaction>
      </trajectory>
    </simulate>

    <logger fileName="output.trees" mode="tree">
      <log spec="TypedTreeLogger" tree="@tree"/>
    </logger>
  </run>
</beast>

Loading the resulting tree file into icytree and colouring the edges using the "type" attribute produces a visualisation in which the tree edges are labeled by their corresponding population types:

Occasionally you really do want to log trees without any additional type annotation. With single-population models this can be achieved as shown in 4.3.1. However, multi-population models (or models with sub-populations) often produce trees with internal nodes representing type changes. These nodes will appear in the output trees, even when the edge types are not recorded. This can also produce problems when using ReMASTER SimulatedTree objects to initialise standard BEAST MCMC analyses, as the extra nodes will not be properly understood.

To log output trees with these nodes removed, use TypedTreeLogger as in the previous section, but add the attribute removeSingletonNodes="true". This will produce a tree in which the internal degree 3 nodes (one parent two children) are still annotated, but in which degree 2 nodes (one parent one child) are removed. Such trees are fully compatible with the rest of BEAST.

Alternatively, as mentioned at the start of this section, one can add the same attribute removeSingletonNodes="true" to the SimulatedTree element, causing these nodes to be stripped from the simulated tree itself. One can then use the standard tree logger as in 4.3.1. This allows one to use the simulated tree directly to initialise an MCMC analysis or to simulate sequences, but has the drawback that the original multi-type tree is lost.

Many stochastic models produce trees containing a random number of leaves. Ensembles of trees containing different numbers of leaves are written to BEAST tree files can cause downstream problems, as such files explicitly assume the number of leaves is the same for all trees they contain.

To circumvent this problem, one can use TypedTreeLogger as above but with the attribute noLabels="true" added.

Occasionally it is useful to record the complete tree together with a sampled tree in which lineages not ancestral to samples have been pruned.

This can be accomplished through the use of the PrunedTree class. This takes a SimulatedTree together with one or more sample populations as input, and produces a representation of the input tree on which the lineages not ancestral to the specified sample populations have been removed (pruned).

For example, the following produces two output files, output.full.trees and output.sampled.trees, which contain the full and sampled trees respectively. (Note the use of an additional "removed" sample population.)

<beast version="2.0" namespace="beast.base.inference:beast.base.inference.parameter:remaster">
  <run spec="Simulator" nSims="1">
    <simulate id="tree" spec="SimulatedTree">
      <trajectory id="trajectory" spec="StochasticTrajectory" maxTime="10">
        <population id="X" spec="RealParameter" value="1"/>
        <population id="Y" spec="RealParameter" value="0"/>
        <samplePopulation id="removed" spec="RealParameter" value="0"/>
        <samplePopulation id="sampled" spec="RealParameter" value="0"/>

        <reaction spec="Reaction" rate="1.4"> X -> 2X </reaction>
        <reaction spec="Reaction" rate="0.1"> X -> Y </reaction>
        <reaction spec="Reaction" rate="0.5"> X -> removed </reaction>
        <reaction spec="Reaction" rate="0.5"> Y -> removed + sampled </reaction>
      </trajectory>
    </simulate>

    <logger fileName="output.full.trees" mode="tree">
      <log spec="TypedTreeLogger">
        <tree spec="PrunedTree" samplePops="removed" simulatedTree="@tree"/>
      </log>
    </logger>
    <logger fileName="output.sampled.trees" mode="tree">
      <log spec="TypedTreeLogger">
        <tree spec="PrunedTree" samplePops="sampled" simulatedTree="@tree"/>
      </log>
    </logger>
  </run>
</beast>

As with all reactions, the n, times, rate and changeTimes inputs can be given multiple values. In the context of coalescent models, this is useful for specifying large numbers of leaves at different times, or for allowing lineage migration rates to change in a piecewise-constant way.

For example,

<reaction spec="PunctualReaction" n="5" times="0 1 2 3"> 0 -> pop </reaction>

produces 5 tips at each of the times 0, 1, 2, and 3; i.e. 20 leaf nodes in total.

Since coalescent reactions are applied to lineages backward in time, the meaning of parents and children is reversed.

For instance, in reaction

<reaction spec="Reaction" rate="1"> S -> L </reaction>

(which might appear in a structured coalescent model) applies to all lineages currently associated with the S population. The "parent" in this reaction is the product L, while the child is the reactant S.

Note that, despite the switch in direction, coalescent reactions follow the same rules as birth-death reactions for assigning parents to children. The ":id" suffix can be applied to both reactant and product to force a parent-child relationship where the default assignment is not correct.