treesimulator

Simulation of rooted phylogenetic trees under a given Multitype Birth–Death (MTBD) model, with or without partner notification (PN).

Preprint

Anna Zhukova, Olivier Gascuel. Accounting for partner notification in epidemiological birth-death-models. medRxiv 2024.09.09.24313296; doi:10.1101/2024.09.09.24313296

MTBD

The MTBD models were introduced by Stadler & Bonhoeffer [Philos. Trans. R. Soc. B 2013].

An MTBD model with m states has

m(m-1) state transition rate parameters:

• μ_ij -- transition rate from state i to state j (1 ≤ i, j ≤ m; i ≠ j), where μ_ij ≥ 0

m² transmission rate parameters:

• λ_ij -- transmission rate from state i (donor) to state j (recipient) (1 ≤ i, j ≤ m), where λ_ij ≥ 0

m removal (becoming non-infectious) rate parameters:

• ψ_i -- removal rate of state i (1 ≤ i ≤ m), where ψ_i ≥ 0

m sampling probability upon removal parameters:

• p_i -- probability to sample the pathogen of an individual in state i upon removal (1 ≤ i ≤ m), where 0 < p_i ≤ 1

Partner Notification (PN)

Partner notification adds two parameters to the initial MTBD model:

• υ -- probability to notify partners upon sampling
• φ -- notified partner removal and sampling rate: φ >> ψ_i ∀i (1 ≤ i ≤ m). The pathogen of a notified partner is sampled automatically (with a probability of 1) upon removal.

We pay particular interest to the classical BD model, the BD Exposed-Infectious (BDEI) model, and BD with superspreading (BDSS), as they are described in [Voznica et al. 2021], and to their -PN versions [Zhukova et al. 2024].

BD

3 parameters:

• λ -- transmission rate
• ψ -- removal rate
• p -- sampling probability upon removal

Epidemiological parameters:

• R₀=λ/ψ -- reproduction number
• 1/ψ -- infectious time

BD-PN

5 parameters:

• λ -- transmission rate
• ψ -- removal rate
• p -- sampling probability upon removal
• υ -- probability to notify partners upon sampling
• φ -- notified partner removal and sampling rate: φ >> ψ

Epidemiological parameters:

• R₀=λ/ψ -- reproduction number
• 1/ψ -- infectious time
• 1/φ -- notified partner removal time

BDEI

2 states:

• E, exposed, i.e. infected but not yet infectious
• I, infectious

4 parameters:

• μ -- transition rate from E to I (becoming infectious)
• λ -- transmission rate from I to E
• ψ -- removal rate of I
• p -- sampling probability upon removal

Epidemiological parameters:

• R₀=λ/ψ -- reproduction number
• 1/ψ -- infectious time
• 1/μ -- incubation period

BDEI-PN

2 states:

• E, exposed, i.e. infected but not yet infectious
• I, infectious

6 parameters:

• μ -- transition rate from E to I (becoming infectious)
• λ -- transmission rate from I to E
• ψ -- removal rate of I
• p -- sampling probability upon removal
• υ -- probability to notify partners upon sampling
• φ -- notified partner removal and sampling rate: φ >> ψ

Epidemiological parameters:

• R₀=λ/ψ -- reproduction number
• 1/ψ -- infectious time
• 1/μ -- incubation period
• 1/φ -- notified partner removal time

BDSS

2 compartments:

• N, standard infectious individual
• S, superspreader

6 parameters:

• λ_nn -- transmission rate from N to N

• λ_ns -- transmission rate from N to S

• λ_sn -- transmission rate from S to N

• λ_ss -- transmission rate from S to S

  (with a constraint that λ_ss/λ_ns=λ_sn/λ_nn)

• ψ -- removal rate of S and of N (the same)

• p -- sampling probability upon removal (the same for N and S)

Epidemiological parameters:

• R₀=(λ_nn + λ_ss)/ψ -- reproduction number
• 1/ψ -- infectious time
• X=λ_ss/λ_ns=λ_sn/λ_nn -- super-spreading transmission ratio
• f=λ_ss/(λ_sn + λ_ss) -- super-spreading fraction

BDSS-PN

2 states:

• N, standard infectious individual
• S, superspreader

8 parameters:

• λ_nn -- transmission rate from N to N

• λ_ns -- transmission rate from N to S

• λ_sn -- transmission rate from S to N

• λ_ss -- transmission rate from S to S

  (with a constraint that λ_ss/λ_ns=λ_sn/λ_nn)

• ψ -- removal rate of S and of N (the same)

• p -- sampling probability upon removal (the same for N and S)

• υ -- probability to notify partners upon sampling

• φ -- notified partner removal and sampling rate: φ >> ψ

Epidemiological parameters:

• R₀=(λ_nn + λ_ss)/ψ -- reproduction number
• 1/ψ -- infectious time
• X=λ_ss/λ_ns=λ_sn/λ_nn -- super-spreading transmission ratio
• f=λ_ss/(λ_sn + λ_ss) -- super-spreading fraction
• 1/φ -- notified partner removal time

Installation

There are 4 alternative ways to run treesimulator on your computer: with docker, apptainer, in Python3, or via command line (requires installation with Python3).

Installation in python3 or command-line

You could either install python (version 3.6 or higher) system-wide and then install treesimulator via pip:

sudo apt install -y python3 python3-pip python3-setuptools python3-distutils
pip3 install treesimulator

or alternatively, you could install python (version 3.6 or higher) and treesimulator via conda (make sure that conda is installed first).

(Optional) to install treesimulator in a new conda environment (e.g., called phyloenv below), first create and activate the environment:

conda create --name phyloenv python=3.6
conda activate phyloenv

Install treesimulator with conda

conda install treesimulator

Basic usage in a command line

If you installed treesimulator in a conda environment (here named phyloenv), do not forget to first activate it, e.g.

conda activate phyloenv

BD and BD-PN

The following command simulates a tree with 200-500 tips under the BD model, with λ=0.5, ψ=0.25, p=0.5, and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv:

generate_bd --min_tips 200 --max_tips 500 \
--la 0.5 --psi 0.25 --p 0.5 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BD-PN model, with λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify only the most recent partner of each sampled index case. The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv:

generate_bd --min_tips 200 --max_tips 500 \
--la 0.5 --psi 0.25 --p 0.5 \
--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bd --help

BDEI and BDEI-PN

The following command simulates a tree with 200-500 tips under the BDEI model, with μ=1, λ=0.5, ψ=0.25, p=0.5, and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv:

generate_bdei --min_tips 200 --max_tips 500 \
--mu 1 --la 0.5 --psi 0.25 --p 0.5 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BDEI-PN model, with μ=1, λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify only the most recent partner of each sampled index case. The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv:

generate_bdei --min_tips 200 --max_tips 500 \
--mu 1 --la 0.5 --psi 0.25 --p 0.5 \
--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bdei --help

BDSS and BDSS-PN

The following command simulates a tree with 200-500 tips under the BDSS model, with λ_nn=0.1, λ_ns=0.3, λ_sn=0.5, λ_ss=1.5, ψ=0.25, p=0.5, and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv:

generate_bdss --min_tips 200 --max_tips 500 \
--la_nn 0.1 --la_ns 0.3 --la_sn 0.5 --la_ss 1.5 --psi 0.25 --p 0.5 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BDSS-PN model, with λ_nn=0.1, λ_ns=0.3, λ_sn=0.5, λ_ss=1.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify only the most recent partner of each sampled index case. The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv:

generate_bdss --min_tips 200 --max_tips 500 \
--la_nn 0.1 --la_ns 0.3 --la_sn 0.5 --la_ss 1.5 --psi 0.25 --p 0.5 \
--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bdss --help

User-defined MTBD and MTBD-PN models

The following command simulates a tree with 200-500 tips under a generic MTBD model, with two states A and B, with μ_aa=0.5, μ_ab=0.6, μ_ba=0.7, μ_bb=0.8, λ_aa=0.1, λ_ab=0.2, λ_ba=0.3, λ_bb=0.4, ψ_a=0.05, ψ_b=0.08, p=_a0.15, p=_b0.65, and saves it to the file tree.nwk, while saving the parameters to the comma-separated file params.csv:

generate_mtbd --min_tips 200 --max_tips 500 \
--states A B \
--transition_rates 0.5 0.6 0.7 0.8 \
--transmission_rates 0.1 0.2 0.3 0.4 \
--removal_rates 0.05 0.08 \
--sampling_probabilities 0.15 0.65 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under a generic MTBD model, with two states A and B, with μ_aa=0.5, μ_ab=0.6, μ_ba=0.7, μ_bb=0.8, λ_aa=0.1, λ_ab=0.2, λ_ba=0.3, λ_bb=0.4, ψ_a=0.05, ψ_b=0.08, p=_a0.15, p=_b0.65, φ=2.5, υ=0.2, and allowing to notify only the most recent partner of each sampled index case. The simulated tree is saved to the file tree.nwk, while the model parameters are saved to the comma-separated file params.csv:

generate_mtbd --min_tips 200 --max_tips 500 \
--states A B \
--transition_rates 0.5 0.6 0.7 0.8 \
--transmission_rates 0.1 0.2 0.3 0.4 \
--removal_rates 0.05 0.08 \
--sampling_probabilities 0.15 0.65 \
--phi 2.5 --upsilon 0.2 --max_notified_partners 1 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_mtbd --help

Basic usage in Python3

To simulate trees with 200-500 tips under the above models and settings:

from treesimulator.generator import generate
from treesimulator import save_forest
from treesimulator.mtbd_models import Model, BirthDeathModel, BirthDeathExposedInfectiousModel, \
  BirthDeathWithSuperSpreadingModel, PNModel

# BD and BD-PN
bd_model = BirthDeathModel(p=0.5, la=0.5, psi=0.25)
print(bd_model.get_epidemiological_parameters())
[bd_tree], _, _ = generate(bd_model, min_tips=200, max_tips=500)
save_forest([bd_tree], 'BD_tree.nwk')
bdpn_model = PNModel(model=bd_model, upsilon=0.2, partner_removal_rate=2.5)
[bdpn_tree], _, _ = generate(bdpn_model, min_tips=200, max_tips=500)
save_forest([bdpn_tree], 'BDPN_tree.nwk')

# BDEI and BDEI-PN
bdei_model = BirthDeathExposedInfectiousModel(p=0.5, mu=1, la=0.5, psi=0.25)
print(bdei_model.get_epidemiological_parameters())
[bdei_tree], _, _ = generate(bdei_model, min_tips=200, max_tips=500)
save_forest([bdei_tree], 'BDEI_tree.nwk')
bdeipn_model = PNModel(model=bdei_model, upsilon=0.2, partner_removal_rate=2.5)
[bdeipn_tree], _, _ = generate(bdeipn_model, min_tips=200, max_tips=500)
save_forest([bdeipn_tree], 'BDEIPN_tree.nwk')

# BDSS and BDSS-PN
bdss_model = BirthDeathWithSuperSpreadingModel(p=0.5, la_nn=0.1, la_ns=0.3, la_sn=0.5, la_ss=1.5, psi=0.25)
print(bdss_model.get_epidemiological_parameters())
[bdss_tree], _, _ = generate(bdss_model, min_tips=200, max_tips=500)
save_forest([bdss_tree], 'BDSS_tree.nwk')
bdsspn_model = PNModel(model=bdss_model, upsilon=0.2, partner_removal_rate=2.5)
[bdsspn_tree], _, _ = generate(bdsspn_model, min_tips=200, max_tips=500)
save_forest([bdsspn_tree], 'BDSSPN_tree.nwk')

# MTBD and MTBD-PN
mtbd_model = Model(states=['A', 'B'], transition_rates=[[0.5, 0.6], [0.7, 0.8]], 
                   transmission_rates=[[0.1, 0.2], [0.3, 0.4]],
                   removal_rates=[0.05, 0.08], ps=[0.15, 0.65])
[mtbd_tree], _, _ = generate(mtbd_model, min_tips=200, max_tips=500)
save_forest([mtbd_tree], 'MTBD_tree.nwk')
mtbdpn_model = PNModel(model=mtbd_model, upsilon=0.2, partner_removal_rate=2.5)
[mtbdpn_tree], _, _ = generate(mtbdpn_model, min_tips=200, max_tips=500)
save_forest([mtbdpn_tree], 'MTBDPN_tree.nwk')

Run with apptainer

Once apptainer is installed, run the following command:

apptainer run docker://evolbioinfo/treesimlator

This will launch a terminal session within the container, in which you can run treesimulator following the instructions for the command line ("Basic usage in a command line") above.