Documentation

Overview

Origins and design

The PhyloJunction project was born from our need to simulate arbitrarily complex SSE processes, but not being able to do so with many existing applications. A series of related diversification models had been implemented in excellent computational methods, but each of those tools made unique modeling assumptions that we wanted to relax.

Instead of reverse-engineering and modifying code written by others, however, it seemed like the path of least resistance was to write our own simulator. We also realized that our program could potentially come in handy in multiple projects, current and future, though only if written as more than a “one-off” chunk of code.

Ideally, PhyloJunction would have to be “modular” in its architecture, both in terms of its codebase as well as how models should be represented and specified. This modularity would allow future models to be more easily implemented – potentially by other developers – and seamlessly patched into PhyloJunction.

With the above in mind, we reasoned our work would be more likely to bear fruits (and faster) if we used a programming language that:

1. was easy to prototype and debug model code in, and

2. cross-platform,

3. supported object-oriented programming, and for which

4. biology and data-science libraries were available.

Python was our language of choice. As for item (4), for example, we could make heavy use of the great Dendropy library, maintained by Jeet Sukumaran and collaborators.

Modularity was achieved by building PhyloJunction around a graphical modeling architecture. This design (used in the past, see below) would not only make PhyloJunction a more flexible simulator, but also more generally useful in the future, as a tool for method development and inference.

Graphical models

PhyloJunction follows the graphical model paradigm popularized in the last decade by programs like RevBayes, BEAST and BEAST 2. In what follows, we shall offer just a brief refresher on what probabilistic graphical models are, and one simple example of their use in evolutionary modeling. Those interested in a more detailed exposition will find it in Höhna et al. (2014) and references therein (also, take a look here).

Just as with the RevBayes and BEAST programs, specifying a model in PhyloJunction amounts to writing down any number of conditional dependencies among random variables. When looked at collectively, random variables comprise a type of probabilistic graphical model called a Bayesian network. As the name suggest, such models can be described as graphs – directed acyclic graphs (DAG), more specifically – where each node (or vertex) represents a random variable, and each edge a conditional dependency between the two nodes (variables) it connects.

A Bayesian network representing a very simple model can be seen in figure 1a:

Figure 1. A simple probabilistic model represented by (a) a Bayesian network, and (b) a factor graph. From the notation in the main text, \(\theta=\{sd, m\}\). The dashed box is called a “plate”, and for each node it envelops it denotes multiple (in this case, precisely five) i.i.d. random variables. White and gray nodes denote random variables whose values we do (i.e., data) and do not know (i.e., parameters), respectively.

If we were to pen down what this Bayesian network represents in probability terms, we would arrive at a joint density (assuming continuous variables) given by expression:

(1)\[f_{\Theta,D}(\theta,d) = f_{\Theta|D}(\theta|D=d)f_D(d) = f_{D|\Theta}(d|\Theta=\theta)f_\Theta(\theta),\]

where \(\theta=\{sd, m\}\).

By comparing the Bayesian network graph and the expression above, the attentive reader will notice that functions are not represented in the graph. We do not know exactly how random variables are related, other than \(d_1\) through \(d_5\) depending somehow on \(sd\) and \(m\). This is why it can be helpful to adopt a more general graphical representation: a factor graph (Fig. 1b).

The model underlying the factor graph in figure 1b is the same as that shown by the Bayesian network, the main difference being the additional factor nodes (the filled squares). Factor nodes can make the relationships between two or more random variables more explicit, especially if their (factor) functions are annotated.

In the example above we can see a factor whose function \(f_{D|\Theta}(d|\Theta=\theta)\) gives us the probability of \(\boldsymbol{d} = \{d_i: 1 \leq i \leq n\}\) given \(\theta\). It is annotated as “Normal”, so we in fact know a close-form expression for this factor function; it is the probability density function of a normal distribution. Random variables \(sd\) and \(m\) stand for the standard deviation and the mean of that normal distribution. (For those curious to see just how complicated factor graphs can be, an example can be found in Zhang et al., 2023; see their Supplementary Fig. X.)

How to read a DAG?

Most biologists specifying a DAG on a computer do so as the first step of statistical inference. They are interested in estimating unknown quantities about the natural world. These quantities – the random variables! – whose values we do not know (but would like to estimate) are referred to as the parameters of the model. Of course, parameter estimation requires that we know the value(s) of one or more random variables, which we naturally refer to as data.

In the example in figure 1, the parameters are \(sd\) and \(m\), jointly referred to as \(\theta\); the data is node labeled \(d_i\). (Note how they are colored differently, gray for data, white for parameter.) Carrying out statistical inference thus implies reading the DAG from its data nodes, normally at the bottom, towards the parameter nodes above.

The most natural approach to parameter estimation requires that we simply take the middle and right-hand side terms in equation (1), and solve for \(f_{\Theta|D}(\theta|D=d)\), the posterior density function:

(2)\[f_{\Theta|D}(\theta|D=d) = \frac{f_{D|\Theta}(d|\Theta=\theta)f_\Theta(\theta)}{f_D(d)}\]

This expression is known as Bayes theorem. What programs like RevBayes, BEAST and BEAST 2 do is evaluate the posterior density function at several values of \(\theta\), and output the resulting (posterior) distribution.

---

PhyloJunction, however, is among other things a collection of simulation (rather than estimation) engines. Borrowing the jargon used above, we are primarily interested in generating values for random variables, some of which are parameters, some data.

Here, it makes sense to think of factors as the distributions from which values will be sampled (though as we will see below, factors can also represent deterministic functions). As opposed to estimation, simulation starts at the top (or “outer”) layers of the DAG, from parameters whose values we do know, and flows in the direction pointed to by arrows (normally downwards; Fig. 1b).

In figure 1b, for example, we start from known values for the parameters of the distributions at the top. The exponential distribution from which we sample (i.e., simulate) \(sd\) has a rate of 1.0; the standard normal distribution (\(Z\)) from which we sample \(m\), by definition, has a mean of 0.0 and a standard deviation of 1.0. We then define a normal distribution from the sampled values of \(sd\) and \(m\), and in turn sample five times from that normal distribution, obtaining our data values in \(\boldsymbol{d}\).

Specifying a model (an example)

PhyloJunction takes a model specification approach that sits between those adopted by the BEAST and RevBayes communities, and that largely intersects with the LinguaPhylo project. DAG-building instructions in PhyloJunction are written in its own programming language, phylojunction (written in lowercase), whose syntax evokes RevBayes’ Rev language. But unlike Rev, phylojunction is not a fully fledged scripting language; it is lightweight and behaves more like a markup language such as XML, BEAST’s format of choice.

Commands in phylojunction can be read as mathematical statements, and are naturally interpreted as instructions for building a node in a DAG. Here is what a series of those commands would look like if written as a phylojunction script:

Example script 1. Script written in phylojunction specifying a time-homogenous birth-death model.

# hyperprior
m <- 0.0 # mean of log birth rate
sd <- 0.1 # standard deviation of log birth rate

# rate values
d <- 1.0 # death
b ~ lognormal(meanlog=m, sdlog=sd) # birth

# deterministic rate containers
dr := sse_rate(name="death_rate", value=d, event="extinction")
br := sse_rate(name="birth_rate", value=b, event="speciation")

O <- 2.0 # origin age

# deterministic parameter stash
s := sse_stash(flat_rate_mat=[dr, br], n_states=1, n_epochs=1) # parameter stash

# phylogenetic tree
T ~ discrete_sse(stash=s, stop="age", stop_value=O, origin="true")

Each line in the above script is a command instructing the application engine to take some form of user input, produce some value from it, and then store that value permanently in a new variable created on the spot.

Every command string consists of an assignment operator (e.g., <-, :=, ~) placed between the variable being created (on its left side) and some user input (on its right side). We can look at some of these commands individually:

• d <- 1.0 creates a variable d (the death rate) that is then passed and henceforth stores constant value 1.0. We use <- when we know the value of a random variable and want to create a constant node in the DAG;

• b ~ lognormal(meanlog=m, sdlog=sd) creates a variable b (the birth rate) that will store a random value drawn from a user-specified distribution. The logarithm of the birth rate is normal with mean m and standard deviation sd.

• dr := sse_rate(name="death_rate", value=d, event="extinction") calls deterministic function sse_rate, which creates variable dr from the value stored in d and some other user input. There is no stochasticity to deterministic assignments; they help make explicit steps where variables are transformed, annotated or combined with others.

And this is the DAG such script instantiates:

Figure 2. Factor graph representing the time-homogenous birth-death model specified in example script 1.

Note how this DAG has a factor node characterized by a different type of function: deterministic function sse_rate. Such functions are denoted by filled diamonds instead of filled squares (those represent distributions), and will have their output enclosed in a hollow diamond.

Model replicates and IID vectors

A replicate is one simulated realization of the complete model. Generate ten replicates with pjcli model.pj --runs 10, or set Replicates in GUI Settings before loading a script. For each replicate, the model is evaluated in statement order, drawing fresh values for the stochastic variables.

Within a replicate, iid draws independent observations with shared distribution parameters, as in the plated node in figure 1:

Example script 2. Five IID observations per replicate.

sd ~ exponential(rate=1)
m ~ normal(mean=0.0, sd=1.0)
d ~ iid(5, normal(mean=m, sd=sd))

Ten replicates produce ten pairs of m and sd, each with five observations in d. Assign constants to m and sd to hold the parameters fixed across replicates. Within a script, d[1] selects the first observation; indices start at one. Nested vectors retain their structure: x[i][j] selects element j of vector x[i].

Assignments and named distributions

A distribution can be named and reused: prior := normal(0, 1) constructs it, and x ~ prior draws a value. Construction alone draws nothing. Reassigning a name affects later statements without recomputing earlier values. A complete rerun evaluates the whole script again. Commands entered interactively apply to every existing replicate. If executing a command fails for any replicate, the command leaves the model and its values unchanged.

Vectorization and comprehensions

A function expecting a scalar argument can instead receive a list, producing one result per entry. Scalar arguments are shared; multiple list arguments are paired and must have equal lengths. A one-element list is not a scalar. Functions that accept lists consume them intact: for example, categorical([1,2], [1,1]) constructs one distribution.

mus <- [0,2]
sds <- [1,3]
distributions := normal(mus,sds)
xs ~ distributions
ys ~ [normal(mu,sd) for (mu,sd) in zip_strict(mus,sds)]

Both xs and ys contain two independent normal draws, with means 0 and 2 and standard deviations 1 and 3, respectively. The two statements make separate draws. zip_strict pairs corresponding entries and rejects unequal lengths.

Vectorization handles one level. Use nested comprehensions for nested lists:

means <- [[0],[1,2]]
values ~ [[normal(mu,1) for mu in row] for row in means]

Here values contains a one-element list followed by a two-element list. ~ samples each distribution in a nested list, preserving the nesting. Lists used with ~ must contain only distributions or lists of distributions.

Comprehensions accept list iterables and bind local names, including flat tuple unpacking as above. Use nested comprehensions for multiple levels; filters and multiple for clauses in one comprehension are not supported. Each statement must occupy one line. Ordinary lists must have compatible element types and nesting; integers and reals may mix, and inner lists may differ in length. Empty lists are allowed.

Graphical user interface (GUI)

PhyloJunction comes with a (simple!) graphical user interface (GUI; see Fig. 3) that allows users to incrementally build a model while simultaneously inspecting any simulation output (among other things).

Figure 3. PhyloJunction’s GUI main window.

Like any modern computer application, PhyloJunction’s GUI exposes its features to users via a menu (Fig. 3, number 2). On the main tab (“Model”), one can build a model (DAG) with commands entered in the prompt (Fig. 3, number 1), or by loading a phylojunction script or a model from a previous PhyloJunction section. The latter can be done by clicking the application’s menu (top-left corner of the Desktop on Apple machines), and clicking “File > Read script” or “File > Load model”. Model building commands can be revisited by going into the menu’s “Command log” tab. Saved models from the older value-bearing DAG format cannot be loaded; load their PJ scripts instead.

Warning

If mistakes are made when entering a command, PhyloJunction will flash a “Warning produced” message in red, at the bottom of the main window. The warning message can be seen by visiting the “Warnings” tab.

Given a model, DAG node values can be navigated as a plot, text string, or both (Fig. 3, numbers 3 and 4). Users can select a replicate and a vector element to inspect its value or plot. A selected tree shows its individual statistics, such as root age and tip count. Across replicates, numeric values and tree statistics receive means and standard deviations. Vector summaries compare corresponding elements across replicates; they do not average elements within a vector. Vectors whose shapes vary between replicates have no aggregate summary. The screenshot uses older control labels.

After a model is specified, saving the main plot (from the “Model” tab) to disk can be done by clicking “Save plot as”. Simulated data tables, or the model itself, in turn, can be saved to disk by clicking the application menu (top menu), and clicking “File > Save data as” or “File > Save model as”.

Implementation comparison

The “Compare” tab (under “Model”) streamlines the comparison of PhyloJunction’s simulations with those from an external simulator, with respect to a summary statistic. There are two main moments when carrying out this type of comparison may be useful or required. The first is when a model is being first implemented in PhyloJunction, and there happens to exist an independent implementation elsewhere. This latter implementation can be used to provide a baseline expectation for PhyloJunction.

Beyond model development within the platform, the “Compare” tab was further designed to expedite activities related to teaching. Students in an evolutionary modeling course may be asked to code their own Yule model, for example, in whatever programming language. After Yule trees are obtained, all they need to do is collect summary statistics that are also monitored by PhyloJunction (e.g., the root age), and upload them into the “Compare” tab. PhyloJunction then automatically plots the two distributions (one from the student, one from PhyloJunction itself) side by side, revealing if implementations behave similarly.

In Settings, set Replicates to 10, then use File > Read script to load examples/multiple_scalar_tree_plated.pj. Its rv5 variable holds two normal observations per replicate; trs2 holds two trees. Compare one component across replicates at a time.

In Compare, click Compare to .csv (…) and load examples/compare_files/comparison_multiple_scalar_tree_plated_rv5.csv. Select rv5[1] from Node to compare, then click Draw. The CSV’s value column must match the selected component name, here rv5[1]. It also contains program, sample, and replicate columns identifying the external observations. The plot compares the two sets of normal draws; small samples can differ visibly.

Figure 4. Comparing a scalar variable in the GUI. This screenshot uses older node labels.

To compare trees, select trs2[1] and choose root_age under Summary statistics. Load examples/compare_files/comparison_multiple_scalar_tree_plated_trs2.csv and click Draw. For trees, the CSV column must match the selected statistic rather than the node name.

Coverage validation

Coverage measures how often posterior intervals contain the true parameter used to simulate each dataset. Across repeated draws from the prior and model, a 95% posterior interval should contain the true value about 95% of the time. Finite samples vary; coverage alone does not prove an implementation correct.

This example uses exact normal posterior intervals to demonstrate the GUI without an external inference program. From the repository root, run:

python examples/coverage_files/prepare_normal.py

The preparation script executes this model 100 times with seed 42 and writes normal.csv in the same directory:

# Use 100 replicates and seed 42 to match prepare_normal.py.
mu ~ normal(0, 1)
y ~ normal(mu, 1, clamp=true)

For this model, the posterior of mu has mean y / 2 and variance 1 / 2. The CSV contains its mean and 95% HPD limits for each replicate, in replicate order.

1. In Settings, set Replicates to 100 and Random seed to 42.

2. Use File > Read script to load examples/coverage_files/normal.pj.

3. In Coverage, select mu under Non-det. nodes.

4. Click Read HPDs .csv (…), load examples/coverage_files/normal.csv, then click Draw.

Figure 5. Exact normal posterior intervals compared with the simulated values of mu.

The horizontal axis shows true parameter values; the vertical axis shows posterior means. Blue intervals contain the true value and red intervals miss it. The top-right panel gives the fraction containing the truth.

For your own analysis, supply one CSV row per replicate, in the same order, with columns posterior_mean, lower_95hpd, and higher_95hpd. Select a scalar parameter or a fixed vector component such as mu[1]. Each component contributes one value per replicate. Changing the model, seed, or replicate count requires matching posterior summaries.

Alternatively, use Directory to .log’s (…) to compute summaries from posterior samples. Supply one tab-separated .log file per replicate, named in replicate order (for example, run1.log, run2.log, …). Set Parameter name in .log if its column name differs from the selected parameter. The GUI uses all rows; remove burn-in samples before loading the files.

Command-line interface (CLI)

PhyloJunction has a command-line interface (CLI) as an entry point for some of its functionalities – it is particularly handy when running PhyloJunction is a step in a longer pipeline (perhaps being run on a cluster). The easiest way to see what is possible with the CLI is to call it on the Terminal with the -h flag, like so:

pjcli -h

The help message should list a series of flags, which we expand on below:

--runs: Number of model replicates (default: 1).

--output-format: separate-csv (default), combined-csv, or json; see Exporting data.

-d: Simulated data tables are written to disk.

-f: Plots for the selected DAG nodes (whose names and replicate ranges are passed as arguments) are written to disk. The argument for -f should be passed within single-quotes, e.g., ‘birth_rate,trs;1-2,1-10’.

-o: Directory where to place simulated data tables and figures.

-p: Prefix to prepend output files with (e.g., “simulation1”)

-nex: Certain simulated data is written to disk in NEXUS format.

-r: Random seed so that simulations always produce the same output.

Note that if no argument is passed to -o, data tables written to disk (prompted by -d) will be placed in the directory from which PhyloJunction is called. Figures (prompted by -f), in turn, will be placed in a “figures/” sub-directory, next to the data tables.

Exporting data

Use pjcli model.pj --runs 10 -d -o output to save simulated values, or select an output format in the GUI before saving data. The CLI --output-format choices are:

• separate-csv (default): one CSV per variable, with @replicate and @element identifying the replicate and vector element. Both use one-based indices; nested element paths look like [2][1]. Scalars have an empty element path.

• combined-csv: one replicate per row, with columns such as x[1] and x[2][1]. Each variable must retain its full vector shape across replicates. Empty vectors appear as [].

• json: a samples.mcon file with an MCON 0.2 header followed by one JSON object per replicate, in replicate order. This format preserves nested vectors of changing lengths.

The GUI labels these formats Separate files (CSV), Combined (CSV), and JSON (MCON). Exports include simulated and assigned values; distributions and other model-construction objects are omitted. Trees appear as annotated Newick strings, with additional files for tree statistics, traits, and fossils. Stochastic maps can be saved separately.

RevBayes inference export

Export a simulated model to estimate its parameters in RevBayes. Put clamp=true on the outer call of a sampling statement to treat its simulated values as observed data. Here mu is inferred from the observations in each replicate:

mu ~ normal(0,1)
observations ~ iid(5, normal(mu,1), clamp=true)

Save this as model.pj. The following Python code writes one Rev script per replicate:

from pathlib import Path
from phylojunction.interface.cmdbox.cmd_parse import compile_model
from phylojunction.inference.revbayes.rb_inference import dag_obj_to_rev_inference_spec

results = compile_model(Path("model.pj").read_text()).execute(runs=10, random_seed=42)
models, analyses, directories = dag_obj_to_rev_inference_spec(
    results, str(Path("inference").resolve()), mcmc_chain_length=1000)
for directory in directories:
    Path(directory).mkdir(parents=True, exist_ok=True)
for i, (model, analysis) in enumerate(zip(models, analyses), 1):
    Path(directories[1], f"run{i}.rev").write_text(model + "\n" + analysis + "\n")

Run a script with RevBayes, for example rb inference/rev_scripts/run1.rev. Its posterior samples are written to inference/inference_results/run1.log. The example uses 1,000 MCMC iterations; choose a suitable run length and assess convergence for an actual analysis.

Export supports Normal, Lognormal, Uniform, Exponential, Gamma, Beta and Cauchy distributions, including named distributions, vectorization, comprehensions, nested lists, and iid. At least one numeric stochastic value must remain unclamped. Tree distributions and other unsupported model functions raise export errors.

Place clamp=true only on the outer sampling call; nested calls, named aliases, and comprehension right-hand sides do not accept it. Collection sizes and distribution choices must stay fixed during inference. Empty collections without a known Rev element type, including literal [], cannot be exported; zero-length iid draws are supported.

Sandbox: Bypassing the interfaces

Python callers can execute scripts and inspect their results directly. From the repository root, this example draws two Yule trees in each of two replicates:

from phylojunction.interface.cmdbox.cmd_parse import compile_model

source = """birth := sse_rate(name="lambda", value=1.0, event="speciation")
stash := sse_stash(flat_rate_mat=[birth])
trees ~ iid(2, discrete_sse(stash=stash, stop="size", stop_value=10))
"""
results = compile_model(source).execute(runs=2, random_seed=42)
trees = results.value("trees", 0)
first_tree = trees[0]
print(first_tree.get_stats_dict())

results.value(name, execution_index) returns a variable’s value in one replicate; results.values(name) returns its values across replicates. Python indices start at zero, including indices into the returned lists. PJ script indices start at one.

For an existing script, pass Path("model.pj").read_text() to compile_model after importing Path from pathlib. Further runnable examples are in src/phylojunction/interface/pj_sandbox.py.

To sample a Python distribution directly, call distribution.draw() for one value or distribution.draw_many(count) for repeated draws with the same parameters.

Annotated tree traits in Python

Python callers describe a tree’s primary trait when constructing an AnnotatedTree:

from phylojunction.data.tree import AnnotatedTree
from phylojunction.data.trait import DiscreteTrait, ContinuousTrait

# Given DendroPy trees with the corresponding node attributes and status flags:
discrete = AnnotatedTree(discrete_tree, DiscreteTrait(states=2))
continuous = AnnotatedTree(continuous_tree, ContinuousTrait())

The default node attribute names are state and trait, respectively. Both descriptions accept an explicit name, such as DiscreteTrait(states=2, name="habitat") or ContinuousTrait(name="body_size"). The named attribute must already exist on the nodes; a description does not rename attributes or their Newick annotations.

Both descriptions provide tree.trait.name; discrete descriptions also provide tree.trait.states. Plotting uses the described node attribute by default. Pass an explicit plotting attribute to display a different node attribute.

Lexicon

Parametric distributions

PhyloJunction provides the following parametric distributions for use in model scripts.

Distribution objects and IID vectors

Construct a distribution with := and draw from it with ~. Use iid for repeated independent draws with the same parameters:

d := normal(0, 1)
x ~ d
y ~ iid(3, d)
z ~ iid(2, iid(3, d))
selected <- z[2][1]

Here x is a scalar, y has three elements, and z contains two three-element lists. iid(n, distribution) requires a nonnegative integer length; zero produces an empty list. Indices start at one and must be integer-valued and within bounds. Use model replicates to redraw the whole model, or vectorization and comprehensions to vary parameters between elements.

Uniform

The function for assigning a uniform (unif) distribution to a random variable has two parameters:

min (real number, required): Lower bound (inclusive).

max (real number, required): Upper bound (inclusive).

Example command. Drawing two IID values of a uniformly distributed random variable \(x\), from -1.0 and 1.0.

x ~ iid(2, unif(min=-1.0, max=1.0))

Exponential

The function for assigning a exponential (exponential) distribution to a random variable has two parameters:

rate (positive real, required): Rate (or scale) of distribution.

rate_parameterization (string boolean, optional): Whether the value passed in rate is the distribution’s rate. Defaults to “true”.

Example command. Drawing two IID values of a exponentially (scale \(\frac{1}{\lambda}=1.0\)) distributed random variable \(x\).

x ~ iid(2, exponential(rate=1.0, rate_parameterization="false"))

Gamma

The function for assigning a gamma (gamma) distribution to a random variable has three parameters:

shape (positive real, required): Shape of distribution.

scale (positive real, required): Scale (or rate) of distribution.

rate_parameterization (string boolean, optional): Whether the value passed in scale is the distribution’s rate. Defaults to “false”.

Example command. Drawing two IID values of a gamma-distributed (shape \(k=10.0\), scale \(\theta=1.0\)) random variable \(x\).

x ~ iid(2, gamma(shape=10.0, scale=1.0))

Cauchy

cauchy(location, scale) describes a real-valued distribution with heavy tails. Both parameters are required: location is a finite median and mode, and scale is finite and strictly positive. Scale is not a standard deviation; the Cauchy has no finite mean or variance. Sample averages need not stabilize with more draws.

initial := cauchy(location=0, scale=1)
x ~ initial
xs ~ iid(10, initial)

Use it as start_distribution=initial in quasse to redraw the origin trait for each biological attempt. examples/quasse_cauchy_origin.pj demonstrates this. Cauchy draws also export to RevBayes as dnCauchy(location=..., scale=...).

Normal

The function for assigning a normal (normal) distribution to a random variable has two parameters:

mean (real, required): Mean (location) of distribution.

sd (positive real, required): Standard deviation (scale) of distribution.

Example command. Drawing two IID values of a normally (mean \(\mu=0.0\), standard deviation \(\sigma=1.0\)) distributed random variable \(x\).

x ~ iid(2, normal(mean=0.0, sd=1.0))

Log-normal

The function for assigning a log-normal (lognormal) distribution to a random variable has three parameters:

meanlog (real, required): Mean (location) of distribution over the logarithm of the random variable.

sdlog (positive real, required): Standard deviation (scale) of distribution over the logarithm of the random variable.

log_space (string boolean, optional): Whether the value of meanlog is in log-space. Defaults to “true”.

Example command. Drawing two IID values of a log-normally (mean log(\(\mu)=-3.25\), standard deviation \(\sigma=0.25\)) distributed random variable \(x\).

x ~ iid(2, lognormal(meanlog=-3.25, sdlog=0.25))

Categorical

The function for assigning a categorical (categorical) distribution to a random variable has two parameters:

cats (list, required): Category values.

probs (list, required): Nonnegative weights, one for each category; total weight must be positive.

Example command. Drawing two IID values of a random variable representing to two equiprobable categories,

x ~ iid(2, categorical(cats=[0,1], probs=[0.5, 0.5]))

Beta

beta(alpha, beta) constructs a distribution on values between zero and one. Both shape parameters must be finite and strictly positive. The mean is alpha / (alpha + beta).

prior := beta(alpha=2, beta=5)
probability ~ prior

Tree distributions

Phylogenetic tree stochastic nodes can be added to a model’s DAG using tree distributions. PhyloJunction provides discrete SSE and continuous-trait QuaSSE tree distributions.

Discrete SSE

Discrete SSE (discrete_sse) stands for discrete state-dependent speciation and extinction process, a general model that includes many popular special cases documented below.

Rate types

The implementation and usage of discrete_sse in PhyloJunction is very flexible and takes inspiration from Tim Vaughan’s great ReMASTER simulator program: a model with any number of rates of whatever types may be specified. Of course, this flexibility comes at the price of user responsibility!

SSE rates are added to a DAG by passing deterministic function sse_rate a stochastic or constant node containing a rate’s value (the value can also be passed directly):

Script excerpt 1. Command in phylojunction to add a birth rate (deterministic node) to the DAG.

birth_rate ~ exponential(rate=1.0)

det_birth_rate := sse_rate(name="lambda", value=birth_rate, event="speciation")

# alternatively, one could pass values directly to sse_rate()
# det_birth_rate := sse_rate(name="lambda", value=[0.95, 1.05], event="speciation")

Parameters for sse_rate include:

name (string, required): The name of the rate parameter. Not particularly critical other than for printing sensical messages to the screen, when using the GUI. The argument to name does not have to match anything else in the model, and can in practice be any string enclosed in double quotes.

value (positive real, required): The rate value. This parameter can also take the name of other nodes in the DAG (in which case the value of those nodes are extracted and used).

event (string, required): A string specifying the rate type (see below).

epoch (positive integer, optional): Integer specifying the epoch to which this rate applies. See more in the “Time-heterogeneity” section below. Defaults to 1.

Arguments for the event parameter should be enclosed in double quotes and can be:

1. "speciation" or "w_speciation" for within-state birth rates;

2. "b_speciation" for between-state birth rates;

3. "asym_speciation" for asymmetric birth rates;

4. "extinction" for extinction rates;

5. "transition" for anagenetic transition rates;

6. "anc_sampling" for ancestor sampling (e.g., “fossilization”) rates.

Sampling probabilities

The discrete SSE model in PhyloJunction supports not only state-dependent rates, but also state-dependent sampling probabilities. Specifying SSE probabilities can be done like so:

Script excerpt 2. Command in the phylojunction language to add a sampling probability parameter (deterministic node) to the DAG.

det_sampling_prob := sse_prob(name="rho", value=1.0, state=0)
# by setting the value to 1.0, we assume complete sampling!

The parameters for sse_prob are:

name (string, required): The name of the probability parameter. Behaves the same as the name of SSE rate parameters.

value (real between 0.0 and 1.0, required): The probability value. As with SSE rates, this parameter can also take the name of DAG nodes that hold probability values.

state (positive integer starting at 0, optional): This integer represents the state to which the probability is associated. Defaults to 0.

epoch (positive integer, optional): Integer specifying the epoch to which this probability applies. See more in the “Time-heterogeneity” section below. Defaults to 1.

Note that the state parameter defaults to 0, and is the state value used in non-state-dependent tree models, like the Yule or birth-death models.

The SSE stash (and time-heterogeneity)

Before discrete_sse can be assigned as the tree distribution to a DAG stochastic node, state-dependent rate and probability parameters must be stashed together and combined with some more information. This is accomplished with deterministic function sse_stash, whose parameters are:

flat_rate_mat (string vector, required): A vector (deterministic) DAG node names holding the objects created with sse_rate. DAG node names do not have to appear in any specific order.

flat_prob_mat (string vector, required): A vector (deterministic) DAG node names holding the objects created with sse_prob. DAG node names do not have to appear in any specific order.

n_states (positive integer, optional): How many states the SSE process has. Defaults to 1.

n_epochs (positive integer, optional): How many epochs characterize the skyline of the SSE process. Defaults to 1.

seed_age (positive real, optional): The age of the node at which the SSE process starts (either the origin or root node). This parameter must receive an argument if parameters are time-heterogeneous, and it must match the value given to stop_value when specifying discrete_sse with stop="age" (see more below).

epoch_age_ends (positive real vector, optional): The age of the ends of each epoch (i.e., the “younger” limit of each epoch, measured as an age), going from the oldest to the youngest epoch. The end of the youngest epoch (which contains the present) should be ommitted as it is added by default. The number of provided age ends should thus be the number of epochs (in n_epochs) minus one. Defaults to [0.0].

The discrete SSE model in PhyloJunction is a “skyline” model, that is, it allows both state-dependent rates and parameters to vary over time (i.e., to be time-heterogeneous). In the presence of time-heterogeneity, these parameters are assumed to vary according to a user-defined piecewise constant function (this is coded into the model via the value and epoch parameters of sse_rate and sse_prob). The limits of each “piece” of that function are given by the argument to epoch_age_ends (see above).

Script excerpt 3. Command in the phylojunction language to organize rate and probability parameters into a stash (a deterministic node in the DAG).

det_birth_rate1 := sse_rate(name="lambda1", value=1.0, event="speciation", epoch=1)
det_birth_rate2 := sse_rate(name="lambda2", value=2.0, event="speciation", epoch=2)

stash := sse_stash(flat_rate_mat=[det_birth_rate1, det_birth_rate2], n_epochs=2, seed_age=2.0, epoch_age_ends=[1.0])
# there are two epochs, the oldest one ends at age 1.0, the youngest one ends at the present
# the birth-rate of the oldest epoch is 1.0, but that rate doubles in the youngest epoch

The discrete SSE distribution

We have seen how to add to the DAG a series of ingredients required by the discrete_sse distribution. The parameters of this function are:

stash (string, required): Name of DAG node holding the object created by sse_stash.

start_state (nonnegative integer, optional): State at which each tree starts. Defaults to 0.

stop (string, required): String specifying the simulation stopping condition. If set to "age", the simulation will stop when the origin or root’s age (see parameter origin) reaches the value provided in stop_value. If set to "size", the simulation will stop when the tree has the number of sampled nodes provided in stop_value.

stop_value (positive numeric, required): The stopping condition threshold. If stop="age", this can be any positive real number. If stop="size", this can be any positive integer.

origin (string boolean, required): Flag specifying if process starts at origin node (as opposed to the root node). Defaults to “true”.

cond_spn (string boolean, optional): Condition the process on at least one speciation event happening. Note that this speciation event may or not be the event represented by the root of the reconstructed tree. Defaults to “false”.

cond_surv (string boolean, optional): Condition the process on having at least one surviving lineage at the present moment. This conditioning only makes sense when stop="age". Defaults to “true”.

cond_obs_both_sides (string boolean, optional): Condition the process on having sampled nodes on both sides of the complete tree’s root node. Defaults to “false”.

min_rec_taxa (positive integer, optional): Minimum number of sampled (extant and sampled ancestors) nodes the tree should have. Defaults to 0.

max_rec_taxa (positive integer, optional): Maximum number of sampled (extant and sampled ancestors) nodes the tree should have. Defaults to 1e12.

abort_at_alive_count (positive integer, optional): The count of a tree’s living nodes (irrespective of node sampling) at which point the tree is considered too large for simulation to continue. Trees can grow out of control when stop="age" and one of the rates allows it to grow too big before the stopping condition is met. This parameter can be used to throw away such trees and preventing PhyloJunction from crashing. Defaults to 1e12.

runtime_limit (positive integer, optional): Maximum number of seconds to obtain one accepted tree. Sampling is aborted at this point. Defaults to 300.

max_n_attempts (positive integer, optional): Maximum number of (failed) attempts before sampling is aborted. Defaults to 200.

Below, an example of assigning a discrete SSE to a stochastic DAG node representing phylogenetic tree samples:

Script excerpt 4. Commands in phylojunction to build a discrete SSE distribution and assign it to a stochastic DAG node.

det_birth_rate1 := sse_rate(name="lambda1", value=1.0, event="speciation", epoch=1)
det_birth_rate2 := sse_rate(name="lambda2", value=2.0, event="speciation", epoch=2)

stash := sse_stash(flat_rate_mat=[det_birth_rate1, det_birth_rate2], n_epochs=2, seed_age=2.0, epoch_age_ends=[1.0])

trs ~ discrete_sse(stash=stash, start_state=0, stop="age", stop_value=2.0, origin="true")

Each distribution draws one tree. Use iid(2, discrete_sse(...)) for two trees with the same parameters, and Replicates or --runs to repeat the complete model. Vectorization and comprehensions allow different parameters for each tree. Each tree has its own runtime and failed-attempt allowance, including retries.

Special cases of discrete SSE tree models

1. Yule (pure-birth) model;

2. Birth-death model;

3. Birth-death with incomplete sampling model;

4. Fossilized birth-death (FBD) model;

5. Binary state-dependent speciation and extinction (BiSSE) model;

6. Time-heterogeneous binary state-dependent speciation and extinction (skyline BiSSE) model;

7. Geographic state-dependent speciation and extinction (GeoSSE) model;

8. Time-heterogeneous geographic state-dependent speciation and extinction (skyline GeoSSE) model;

Reading trees and selecting elements

read_tree always returns a list, even when the input contains one tree. Select a scalar tree explicitly with [1]:

trees <- read_tree(string="((A:1,B:1):1,C:2);")
first <- trees[1]

Reading and applying stochastic histories

A history records the state changes along a tree’s branches for one stochastic-map iteration. Read the trees and histories, then attach each history to a copy of the corresponding tree:

trees := read_tree(file_path="examples/trees_maps_files/geosse_dummy_tree1.tre", node_name_attr="index")
histories := read_maps(file_path="examples/trees_maps_files/geosse_dummy_tree1_maps.tsv", tip_attr_file_path="examples/trees_maps_files/geosse_dummy_tree1_tip_states.tsv", n_regions=2, geo=true)
mapped := map_attr(tree=trees[1], history=histories, attr_name="state")

read_maps returns one history per iteration identifier in first-appearance order. map_attr copies its tree before attaching the selected history. A scalar tree is shared across a vector of histories; equally sized tree and history lists pair elementwise. Singleton lists are not expanded: trees[1] selects the one tree explicitly. Both operations also work inside comprehensions. History node identifiers must match the tree’s index attributes; node_name_attr="index" reads those identifiers from the Newick annotations.

QuaSSE

QuaSSE couples a real-valued Brownian trait to trait-dependent speciation, extinction and optional fossil sampling. Two simulation methods are available: the default diversitree approximation and local_thinning, which evaluates rates at Brownian-bridge candidate traits and has an explicit probability-error budget under the assumptions below.

With PhyloJunction installed, run these commands from the repository root (where examples/ is located):

pjcli --runs 2 examples/quasse.pj -d -f 'trees;1-2' -r 42 -o /tmp/quasse-output
pjcli --runs 2 examples/quasse_size.pj -d -r 42 -o /tmp/quasse-size
pjcli --runs 2 examples/quasse_conditioned.pj -d -r 42 -o /tmp/quasse-conditioned
pjcli --runs 2 examples/quasse_local_thinning.pj -d -r 42 -o /tmp/quasse-local

Each .pj statement must occupy one line. -d writes data; -f selects figures. The GUI’s normal tree display also colors continuous terminal traits.

Model

# Run the complete model with pjcli --runs 2 (or set Replicates in Settings).
# Brownian diffusion is variance per unit time. All commands occupy one line.
birth := quasse_logistic(y0=0.4, y1=0.8, midpoint=0.0, slope=2.0)
death := quasse_logistic(y0=0.02, y1=0.1, midpoint=1.0, slope=2.0)
trees ~ iid(2, quasse(birth_rate=birth, death_rate=death, start_trait=0.0, drift=0.0, diffusion=0.1, stop='age', stop_value=3.0, method='diversitree', k=500))

Download quasse.pj.

quasse_logistic matches Java SSE’s link function:

f(x) = y0 + (y1 − y0) / (1 + exp(slope × (midpoint − x)))

The plateau values are finite and nonnegative. The midpoint and slope are finite. Either slope sign and either plateau ordering are supported. Equal plateaus produce a rate independent of the trait. Rate functions have no explicit time dependence; rates change as traits evolve.

Other rate families are available for speciation, extinction or fossil sampling:

Constructor and required arguments	Rate at trait x
quasse_constant(rate)	rate
quasse_gaussian(baseline, center_rate, center, width)	baseline + (center_rate − baseline) exp(−z²/2), z = (x − center)/width
quasse_step(left, right, threshold)	left when x < threshold; right otherwise, including equality
quasse_linear(intercept, slope)	max(0, intercept + slope × x)
quasse_skew_gaussian(baseline, amplitude, location, width, skew)	baseline + 2 amplitude exp(−z²/2) Φ(skew × z), z = (x − location)/width
quasse_quadratic(baseline, strength, center)	baseline + strength × (x − center)²
quasse_absolute(baseline, strength, center)	baseline + strength × abs(x − center)

Supply arguments by name in scripts. All parameters must be finite; rate levels, amplitudes and strengths must be nonnegative, and widths positive. Locations, slopes, skew and the linear intercept can have either sign. Gaussian curves include troughs when center_rate < baseline. Φ is the standard normal cumulative distribution function. The asymmetric peak’s location and amplitude are not generally its peak position and height; skew=0 gives a symmetric Gaussian peak. All constructor parameters are constant through time within a simulation.

Every constructor accepts optional cap: the evaluated rate is min(cap, f(x)). The cap must be finite and nonnegative, may be below the baseline, and may be zero. Without it, linear, quadratic and absolute-value rates remain unbounded when their slope or strength is nonzero. There is no hidden upper truncation.

An explicit cap changes the biological model; computational resource limits do not cap rates. For example, pure branching Brownian motion with quadratic speciation can have an infinite expected population after a finite time even though each finite-time population is finite almost surely. A large realization or exceeded node budget is not proof of infinitely many nodes. General birth–death models should not be assigned an explosion classification from these runtime checks.

The trait increment over duration δ has mean drift × δ and variance diffusion × δ. diffusion is variance per unit time, not a standard deviation. Drift and diffusion are independent of the trait and constant through time within each simulation. Zero diffusion is supported. Each tree starts with one lineage at the origin. Daughters inherit the parent’s trait and then evolve independently. Optional fossil_rate supplies a rate object for non-destructive sampling along each lineage, in observations per lineage per unit time. It defaults to zero (None in Python) and accepts every rate family. Fossils record the exact trait at sampling; they do not remove the lineage or reset its trait. Sampling at death, removal probabilities and measurement error are not supported.

Rate-constructor and distribution parameters are scalar within a draw. Use iid to draw several trees with the same parameters and model replicates to redraw the whole model. In PJ expressions, list arguments to scalar parameters construct lists of distributions through one-level vectorization; comprehensions give explicit control over nested collections. Rejected trees do not resample DAG parents. A start_distribution is instead drawn inside the tree sampler for each biological attempt; a numerical restart retains that starting trait. For example, start_distribution=normal(0, 1) passes a distribution object directly. For a heavy-tailed starting trait, use start_distribution=cauchy(0, 1); the second argument is a scale, not a standard deviation.

Diversitree method and stopping

At the current living traits, compute total birth–death–sampling rate R. A full step has duration 1/(kR), and one event occurs with probability 1/k. Choose the event type and lineage using the current rates. Apply the event first, then extend and evolve the resulting living lineages. A newborn daughter therefore evolves during that step, while a dead lineage does not. A fossil records the start-of-step trait and its continuing lineage evolves for the whole step.

For stop="age", shorten the final interval to the remaining duration δ and use event probability Rδ. Surviving branches end exactly at the observation age. Completely extinct trees retain the actual extinction times, with node ages measured relative to observation.

For stop="size", N counts sampled living tips. The initial lineage receives a Bernoulli(sampling_prob) mark. At birth, one randomly chosen daughter inherits its parent’s mark and the other receives an independent fresh mark. Thus a birth adds zero or one to the sampled count. A sampled death subtracts one; fossils preserve the continuing mark.

Stop immediately before a proposed birth would increase that count from N to N+1, omitting that birth and its following trait update. Merely reaching N does not stop the simulation: zero-increment births can occur at N and sampled deaths can lower the count again. Do not resample marks at the end. 0 < sampling_prob <= 1 is required for size stopping; zero is permitted for age stopping. All living lineages, including unsampled ones, evolve and count toward max_alive. Zero sampled tips with living lineages is not extinction. If survival conditioning is disabled, extinct outcomes can also be returned.

On a fixed biological tree independent of the marks, these leaf marks are independent Bernoulli draws. The flag-dependent stopping time conditions their distribution; the returned size-stopped tree is defined by this marked process, not independent post-sampling of a tree stopped by its total living count. At complete sampling the old behavior is preserved.

The interior algorithm follows diversitree/R/simulate-quasse.R. Intentional differences are:

• The final age-limited interval and event probability are shortened, preventing overshoot.

• Size stopping omits diffusion after the terminating birth, matching PhyloJunction’s convention.

• No artificial initial branch length of 1e-8 is added; genuine zero-length branches are valid.

• Tree storage, observation sampling, and conditioning use PhyloJunction’s conventions.

The method is approximate. Increasing k generally improves resolution at greater cost; k does not specify an absolute time-step tolerance. It defaults to 500 and must be an integer at least one. The default method is "diversitree".

Optional dt_max is a positive finite duration. When supplied, each interval is at most min(1/(kR), dt_max, remaining age), with event probability R times that interval. This can resolve trait changes where present rates are small, but is not an error tolerance or an exact simulation method. Step discontinuities and rapidly varying or unbounded rates require resolution checks using smaller dt_max and larger k. Omitting it preserves the original step selection for positive total rates.

Brownian bridge local thinning

Select method="local_thinning" and supply a positive finite block_duration in the model’s time units. Each living lineage samples a Brownian endpoint up to that duration away. The endpoints and a Brownian-bridge tail bound define a likely containing trait interval. The sum of bounds on birth, death and fossil rates over that interval bounds a Poisson proposal process. At each candidate time the simulator samples the actual intermediate trait from the stored bridge and selects birth, death, fossil sampling or rejection using actual rates. Rejected candidates and fossil events retain the same endpoint and proposal bound. Sampling only replaces the living node by its continuing child; it does not allocate another bridge or error allowance. At birth the parent ends at that time and its daughters evolve independently. There is no one-event-per-step restriction and no requirement that rates change slowly.

Longer blocks save endpoint work but can enlarge proposal bounds and increase rejected candidates. Thus block_duration primarily controls cost. It need not be extremely small for accuracy. Local thinning does not use dt_max or k: a supplied dt_max or nondefault k is rejected; the inherited default k=500 is ignored.

bridge_error defaults to 1e-8 and must be finite and strictly between zero and one. It budgets excursions outside the likely trait regions for one accepted-tree draw(), including its biological rejections and numerical restarts. Every tree in an IID vector gets its own runtime limit, failed-attempt limit and bridge-error allowance. Block j receives probability allowance δ/[j(j+1)], whose infinite sum is δ. A direct simulate() call gets its own budget. Separate calls start new budgets. All existing rate families are supported, including uncapped linear, quadratic and absolute-value functions: a finite local bound suffices, without a global rate cap.

A finite evaluated total rate exceeding its proposal bound automatically restarts the current tree with the same model parameters. RNG progression, error budget, deadline and remaining work allowance persist. Completed earlier trees are retained. This does not consume a biological rejection attempt. A loose valid bound merely wastes candidate work; an invalid bound can miss events, so restarting does not make this an exact sampler.

With exact arithmetic and sampling, valid interval bounds, and a terminating, nonexplosive reference process, the total variation distance from the exact complete output law is at most bridge_error. The coupling agrees until the first region excursion, which also covers detected violations and their restarts. This is an absolute probability bound, not a relative bound on population means or other unbounded statistics. It is not a claim of second-order convergence in block duration. Floating-point error and computational interruptions are outside this guarantee. Biological rejection-limit failure is included as an output; discarding failures and conditioning on success needs a separate bound (at most 2δ/p, capped at one, where p is exact success probability and both are positive).

Age stopping ends living branches exactly at the requested age. Size stopping retains the same N-to-N+1 convention as diversitree, but uses the continuous candidate birth time and advances every unfinished bridge to that time. The stored future endpoints are retained when sampling final traits. Both methods use the same inherited marks for size stopping and independent end-of-simulation sampling for age stopping. The probability bound includes the marked stopping law: before a region excursion, shared mark draws give identical stopping decisions. Adding fossil intensity to the same region’s bound needs no extra error allowance.

The derivation and assumptions are in TODO/simulator.tex in the source repository.

Observation and conditions

Argument	Default	Meaning
sampling_prob	1	Independent terminal sampling for age stopping; inherited Bernoulli marks for size stopping
cond_surv	"true"	Require at least one living lineage
cond_spn	"false"	Require a retained speciation event
cond_obs_both_sides	"false"	Require living sampled descendants on both sides of the original root
min_rec_taxa	0	Minimum living sampled-tip count, for age stopping
max_rec_taxa	unbounded	Maximum living sampled-tip count, for age stopping

Equal count bounds request an exact observed count. Sampling precedes these checks. Survival does not imply observation: with sampling_prob=0, a surviving tree has no sampled living tips, but may retain fossil observations. Root-side checks use the original root, before pruning can replace it. For unconditional age-based simulation set cond_surv="false" and leave other conditions unset. Explicit reconstructed-taxon count bounds with size stopping are errors. Fossils do not satisfy living-tip count, root-side or survival conditions; sampling alone is not speciation.

Rejected attempts use fresh simulation randomness with unchanged model parameters. Rare or impossible conditions can exhaust the rejection budget. Limits are:

Argument	Default	Scope
max_n_attempts	200	Failed attempts allowed while obtaining one tree with draw()
runtime_limit	300 seconds	One draw() call, including retries
max_steps	1,000,000	Diversitree steps or local queue actions per attempt, shared across numerical restarts
max_alive	100,000	Living lineages per attempt
max_nodes	unset	All nodes per attempt, including extinct nodes, fossils, sampling dummy nodes, origin and initial lineage

All limits are positive; count limits are integers. Reaching the rejection budget raises an error. Except for the finite local bound-exceedance recovery described above, numerical failures and exceeded resource limits raise errors without retry. In particular, transiently exceeding max_rec_taxa does not reject an attempt.

max_nodes is checked before creating the initial two nodes or the two new nodes at a birth or fossil event. The omitted size-stopping birth does not consume nodes. This is a memory safeguard, not an infinite-population detector; hitting it raises an error without returning a truncated tree.

When all three rate functions are identically zero, age stopping uses a direct Brownian transition, even with dt_max set. Identically zero birth and death rates make size stopping unreachable even if fossil sampling is positive. For diversitree, if the total rate is zero only at the current traits (including underflow), supply dt_max so traits can advance and rates be reevaluated; without it, simulation raises an error. Local thinning instead renews Brownian blocks even when a local proposal bound is zero.

Trees, plots, and output

Complete trees retain extinct and unsampled tips. Reconstructed trees retain paths to sampled living tips and fossils, including fossils with extinct or unobserved descendants and trees with fossils alone. A lone fossil retains its original age and origin-to-observation stem. All complete-tree nodes carry the trait at that node’s time; no full trait trajectory is stored. Black branches and colored terminal markers show traits using a continuous colorbar. Ancestral fossils appear as trait-colored markers along their continuing branches; terminal fossils have their own rows. sa_along_branches=False gives every fossil its own row. The complete tree’s terminal range, including fossils, determines colors in both views. In complete-tree plots, paths from retained observations back to the full-tree root are drawn normally; other edges are lighter (25% opacity) and thinner (65% linewidth). This includes any displayed ancestral stem above the reconstructed root. Fossil markers, tip colors and labels remain unchanged, and reconstructed plots retain their usual styling.

A light dotted vertical line marks the present, behind the branches and node labels. It uses the simulation’s observation horizon, so fossil-only or extinct trees can end before the line. The same reference appears in complete and reconstructed views and in saved plots. Plotting uses no randomness. An empty reconstructed tree displays No sampled tips.

Data output includes complete/reconstructed Newick tables, annotated variants, generic tree summaries, and <node>_traits.tsv. The trait table has columns @replicate, @element, node, trait, alive, and sampled; replicate indices and element-path components start at one. It includes internal nodes and the origin as well as living and extinct tips. Discrete state tables and discrete NEXUS character matrices are not generated for these trees.

<node>_fossils.tsv contains @replicate, @element, node, time, age, trait. Replicate indices and element-path components are one-based. time is forward from the origin, and age is relative to the simulation horizon. There is one row per fossil in the complete tree, or a header-only table when none occur. Annotated Newick identifies fossils with is_sa and sampling nodes with is_sa_dummy_parent. The generic Direct ancestor count counts all fossil observations, including those made terminal by reconstruction. Total taxon count retains its existing biological-terminal meaning.

QuaSSE RevBayes inference export is not supported and fails before script generation.

Python interface

import numpy as np
from phylojunction.calculation.continuous_sse import LogisticRate
from phylojunction.distribution.dn_quasse import DnQuaSSE

np.random.seed(42)
simulator = DnQuaSSE(
    birth_rate=LogisticRate(0.4, 0.8, 0.0, 2.0),
    death_rate=LogisticRate(0.02, 0.1, 1.0, 2.0),
    stop="size", stop_value=10, diffusion=0.1,
)
tree = simulator.draw()
trees = simulator.draw_many(2)

draw() applies conditions and returns one tree. simulate() returns one attempt after observation sampling, without rejection conditioning. draw_many(2) draws two trees with independent resource allowances. The random_seed argument to execute() and CLI -r seed both Python and NumPy’s global random streams; distribution objects do not own RNG seeds.

The Python rate classes are ConstantRate, LogisticRate, GaussianRate, StepRate, LinearRate, SkewGaussianRate, QuadraticRate and AbsoluteRate, all in phylojunction.calculation.continuous_sse. They accept the same named arguments as the script constructors and evaluate NumPy arrays or scalars. Their shared RateFunction interface exposes is_identically_zero and upper_bound (a conservative finite bound, or None when unavailable). A zero rate at one trait is not an identically zero function. bound_on_interval(left, right) returns a finite conservative bound on a closed finite interval, including zero-width intervals, or raises ValueError if unavailable. Positive bounds have an outward floating-point guard; this is not certified interval arithmetic.

Additional scripts

Size stopping (quasse_size.pj):

# Run the complete model with pjcli --runs 2 (or set Replicates in Settings).
# Return ten living tips, immediately before the birth that would produce eleven.
birth := quasse_logistic(y0=0.4, y1=0.8, midpoint=0.0, slope=2.0)
death := quasse_logistic(y0=0.02, y1=0.1, midpoint=1.0, slope=2.0)
trees ~ quasse(birth_rate=birth, death_rate=death, diffusion=0.1, stop='size', stop_value=10, cond_surv='true', method='diversitree', k=500)

Conditioning on observed tip count (quasse_conditioned.pj):

# Run the complete model with pjcli --runs 2 (or set Replicates in Settings).
# Fix age, sample each survivor independently, and retain exactly three observed tips.
birth := quasse_logistic(y0=0.4, y1=0.8, midpoint=0.0, slope=2.0)
death := quasse_logistic(y0=0.02, y1=0.1, midpoint=1.0, slope=2.0)
trees ~ quasse(birth_rate=birth, death_rate=death, diffusion=0.1, stop='age', stop_value=3.0, sampling_prob=0.8, min_rec_taxa=3, max_rec_taxa=3, method='diversitree', k=500)

Bounded and explicitly unbounded rates (quasse_rates.pj):

# Caps change rates; max_nodes only aborts excessive allocations. dt_max is a duration.
bounded := quasse_gaussian(baseline=0.1, center_rate=0.8, center=0, width=1, cap=0.6)
unbounded := quasse_quadratic(baseline=0.1, strength=0.05, center=0)
death := quasse_constant(rate=0.02)
bounded_trees ~ quasse(birth_rate=bounded, death_rate=death, diffusion=0.1, stop="age", stop_value=2, dt_max=0.01, max_nodes=10000)
unbounded_trees ~ quasse(birth_rate=unbounded, death_rate=death, diffusion=0.1, stop="age", stop_value=2, dt_max=0.01, max_nodes=10000)

Local thinning with uncapped quadratic speciation (quasse_local_thinning.pj):

# Run the complete model with pjcli --runs 2 (or set Replicates in Settings).
# Uncapped quadratic speciation with Brownian bridge local thinning.
b := quasse_quadratic(baseline=0.6,strength=0.05,center=0)
d := quasse_constant(rate=0.2)
trees ~ iid(2, quasse(birth_rate=b, death_rate=d, stop='age', stop_value=2, start_trait=0, drift=0, diffusion=1, method='local_thinning', block_duration=0.5, bridge_error=1e-08, cond_surv='false', max_nodes=10000))

Fossil sampling and incomplete size sampling (quasse_fossils.pj):

# Run the complete model with pjcli --runs 2 (or set Replicates in Settings).
birth := quasse_constant(rate=1.5)
death := quasse_constant(rate=0.5)
fossil := quasse_constant(rate=0.5)
mixed ~ iid(2, quasse(birth_rate=birth, death_rate=death, fossil_rate=fossil, stop='age', stop_value=4, diffusion=0.5, sampling_prob=0.5, cond_surv='false', method='local_thinning', block_duration=0.5, max_nodes=10000))
fossil_only ~ iid(2, quasse(birth_rate=birth, death_rate=death, fossil_rate=fossil, stop='age', stop_value=4, diffusion=0.5, sampling_prob=0, cond_surv='false', method='local_thinning', block_duration=0.5, max_nodes=10000))
sampled_size ~ iid(2, quasse(birth_rate=birth, death_rate=death, fossil_rate=fossil, stop='size', stop_value=15, diffusion=0.5, sampling_prob=0.5, method='local_thinning', block_duration=0.5, max_nodes=10000))

The three calls produce mixed observations, fossil-only observations, and size-stopped trees with exactly five sampled living tips. The last call can retain additional unsampled living tips and any number of fossils; neither contributes to its target count.

Reading complete trees

read_tree(..., is_complete=true) restores biological status from Newick annotations. Each ordinary tip requires alive=true or alive=false; dead tips additionally require sampled=true or sampled=false. Sampled fossil tips survive reconstruction. Explicit status resolves near-present fossils even when branch lengths are rounded; lengths are unchanged. keep_unsampled_alive=true treats explicitly unsampled living tips as sampled for reconstruction.

The default, is_complete=false, reads an already reconstructed tree: every tip is sampled, and living status is inferred from depth. Explicit living-status annotations or unsampled tips are rejected in this mode. Both modes return a list of trees; select one with [1].

n_rec_branches(tree) counts reconstructed edges with a parent and duration greater than that tree’s numerical tolerance. Sampled-ancestor marker edges are excluded. Edges on either side of an observation remain separate. Branch-specific vectors follow preorder traversal. Hidden-speciation counts (n_stub) belong to each reconstructed edge separately; observations do not transfer counts from their incoming edge to their continuing edge.

Gamma-spike branch rates

gamma_spike_clock(tree, attr_name, clock_rate, spike_mean, spike_shape, br_specific_rel_rate, strict_clock_when_no_stubs=false, apply_to_complete=false, br_specific_rel_rate_complete=none) describes a distribution on a copy of one tree with rates on its reconstructed edges. Draw with ~. The input’s times, traits and stochastic maps are preserved. br_specific_rel_rate is a vector in the preorder used by n_rec_branches; generate it separately with iid.

The gradual component is clock_rate * relative_rate. Each hidden speciation on an edge and each genuine parent speciation contributes a Gamma spike with mean spike_mean and shape spike_shape. Sampling observations contribute no spike. Dividing the total spike size by edge duration gives its additional rate. spike_mean=0 removes the spike effect. All parameters must be finite; relative rates and spike_mean may be zero, whereas clock_rate and spike_shape must be positive. The optional strict_clock_when_no_stubs forces stub-free edges to the constant clock_rate; it is a diagnostic option, not the usual statistical model.

apply_to_complete=true additionally annotates every branch of the complete tree, not just the reconstructed tree’s – this only has an effect on a simulated tree, since a tree read in from Newick has no complete topology behind it to annotate. Each reconstructed branch’s own, already-drawn spike total is split across the complete-tree branches hidden along it (a Dirichlet draw, following the Gamma-to-Dirichlet connection: i.i.d. Gamma draws summing to a known total have their normalized shares distributed Dirichlet, independently of that total), while branches entirely absent from the reconstructed tree – extinct lineages, or ones still alive today but never sampled – each receive a fresh, independent draw instead. br_specific_rel_rate_complete supplies that fresh draw’s relative rates, one per branch counted by n_complete_only_branches(tree), in the same preorder n_complete_only_branches and complete_only_branch_nodes agree on; it is required when apply_to_complete=true and ignored otherwise.

n_rec_branches/n_complete_only_branches can both be zero for a given tree (e.g. an unusually small simulated replicate, or no extinct/unsampled lineages at all) – an empty br_specific_rel_rate/br_specific_rel_rate_complete vector is valid in that case, not an error.

Displaying rates

A tree with rate annotations initially displays the reconstructed tree in evolutionary-distance units (time × rate). Uncheck Evolutionary distance to use time units and enable the Reconstructed checkbox. Switching back to time restores its previous reconstruction choice. The present-time reference line appears only in time views.

GUI summaries describe the selected complete or reconstructed topology. Ages remain in time units in both views. Exported statistics and comparison/coverage calculations are independent of the display selection. Python plotting also accepts branch_rate_attr to color edges by a continuous rate attribute without replacing biological trait annotations.

The complete tree can also be plotted in evolutionary-distance units from Python, provided apply_to_complete=true was used – call tree.plot_node(axes, evolutionary_distance=true, draw_reconstructed=false) directly. Without apply_to_complete, the complete tree has no rates to plot against, and this raises rather than falling back to the reconstructed tree or producing an incomplete plot. The GUI itself currently only ever displays the reconstructed tree in evolutionary-distance units.

Exporting stub counts

Add --stub-count-output (or -nstub) alongside -d to write one <variable>_n_stub.tsv table per tree-valued variable. Rows identify the model replicate and nested vector element; columns give reconstructed-edge stub counts by node label. A missing node is an empty cell, not a zero count. In Python, use write_results(..., write_stub_counts=True). This option does not depend on GUI views.

Reading and mapping external trees

read_tree(file_path=...) reads one Newick tree per line; string=... reads a literal Newick tree. Specify exactly one. node_name_attr="index" uses an annotation to name nodes, which lets external stochastic histories refer to the same branches. epoch_age_ends retains the existing time-slice interface. Use [1] to select one tree from the returned list. For example, a complete tree can distinguish sampled fossils from unsampled extinction:

trees <- read_tree(string="(living[&alive=true]:2,fossil[&alive=false,sampled=true]:1,extinct[&alive=false,sampled=false]:0.5);", is_complete=true)

read_maps(file_path, tip_attr_file_path, n_regions=1, geo=false) reads the external histories. The maps TSV identifies anagenetic and cladogenetic changes by node index; the tip TSV supplies tip states. geo=true interprets states as combinations of areas across n_regions. map_attr(tree, history, attr_name="state") returns a copy annotated with one history, leaving topology and times unchanged. A vector of histories maps the same tree once per history. These discrete annotations are retained when rates are added later.

Runnable examples

The first example assigns rates to two FBD trees read in from Newick – already reconstructed, so there is no complete-only branch to annotate. The second retains every external history as a separate mapped tree before adding rates. The third simulates (rather than reads in) an FBD tree, so it has a real complete topology behind it, and uses apply_to_complete to also rate the branches absent from the reconstructed tree. Use pjcli --runs to repeat the whole model.

# Two reconstructed FBD trees; model replicates are selected with pjcli --runs.
trs <- read_tree(file_path="examples/trees/two_FBD_trees.tre")
# Mean-one relative rates: meanlog = -sdlog^2/2, with sdlog=0.5.
n_rec_br := n_rec_branches(trs)
br_specific_rel_rate ~ [iid(n, lognormal(-0.125, 0.5)) for n in n_rec_br]
# Each pair supplies one tree and its own vector of rates in reconstructed preorder.
phylo ~ [gamma_spike_clock(t, "mutation_rate", 0.001, 0.01, 2.0, r) for (t, r) in zip_strict(trs, br_specific_rel_rate)]

# A GeoSSE tree with every supplied stochastic history retained as a mapped tree.
trs <- read_tree(file_path="examples/trees_maps_files/geosse_dummy_tree2.tre", node_name_attr="index")
histories := read_maps(file_path="examples/trees_maps_files/geosse_dummy_tree2_maps.tsv", tip_attr_file_path="examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv", n_regions=2, geo=true)
mapped_trs := map_attr(tree=trs[1], history=histories, attr_name="state")
# Clock draws preserve the discrete maps; each mapped tree receives its own rate vector.
n_rec_br := n_rec_branches(mapped_trs)
br_specific_rel_rate ~ [iid(n, lognormal(-0.125, 0.5)) for n in n_rec_br]
phylo ~ [gamma_spike_clock(t, "mutation_rate", 0.001, 0.01, 2.0, r) for (t, r) in zip_strict(mapped_trs, br_specific_rel_rate)]

# FBD tree, simulated (speciation + extinction + ancestor sampling,
# i.e. serial/fossil sampling through time) rather than read in from
# Newick -- see 'examples/fbd.pj' for the simulation alone, and
# 'examples/fbd_read_gamma_spike.pj' for the read-in-tree version of
# this same clock-model pipeline. Run the complete model with
# 'pjcli --runs 2' (or set Replicates in Settings).
birth_rate <- 1.0
death_rate <- 0.2
fossil_rate <- 0.5
det_birth_rate := sse_rate(name="lambda", value=birth_rate, states=[0,0,0], event="w_speciation")
det_death_rate := sse_rate(name="mu", value=death_rate, states=[0], event="extinction")
det_fossil_rate := sse_rate(name="psi", value=fossil_rate, states=[0], event="anc_sampling")
stash := sse_stash(flat_rate_mat=[det_birth_rate, det_death_rate, det_fossil_rate], n_states=1)
ori <- 3.0
# 'min_rec_taxa'/'max_rec_taxa' reject-and-resimulate until the
# reconstructed tree has between 5 and 15 extant, sampled tips -- a
# plain age-stopped birth-death process can otherwise occasionally
# produce a tree with very few (or very many) tips, depending on how
# the simulation happens to play out
trs ~ discrete_sse(stash=stash, start_state=0, stop="age", stop_value=ori, origin="true", cond_spn="true", cond_surv="true", min_rec_taxa=5, max_rec_taxa=15)

# per-branch relaxed-clock rates for the reconstructed tree's own
# branches, drawn independently of the clock model below; mean 1 in
# real space (meanlog = -sdlog^2 / 2), sdlog=0.5 as in Douglas et al.
# (2025)'s own empirical analyses ('clockSD')
n_rec_br <- n_rec_branches(trs)
br_specific_rel_rate ~ iid(n_rec_br, lognormal(-0.125, 0.5))

# unlike a tree read in from Newick, a simulated tree has a real,
# known complete topology behind it -- so it also has branches absent
# from the reconstructed tree (going extinct, or surviving to the
# present without being sampled), each needing its own, independently
# drawn relative rate
n_absent_br <- n_complete_only_branches(trs)
br_specific_rel_rate_complete ~ iid(n_absent_br, lognormal(-0.125, 0.5))

# Gamma-spike punctuated relaxed clock (Douglas et al. 2025), with
# clock_rate/spike_mean/spike_shape as in their empirical analyses --
# 'apply_to_complete="true"' additionally annotates the complete
# tree's own branches, not just the reconstructed tree's
phylo ~ gamma_spike_clock(trs, "mutation_rate", 0.001, 0.01, 2.0, br_specific_rel_rate, apply_to_complete="true", br_specific_rel_rate_complete=br_specific_rel_rate_complete)

The Gamma-spike model follows Douglas, Bouckaert, Harris, Carter, and Wills (2025), Evolution is coupled with branching across many granularities of life, Proceedings of the Royal Society B 292:20250182, https://doi.org/10.1098/rspb.2025.0182.