treesimulator

treesimulator

Simulation of rooted phylogenetic trees under a given Multi-Type Birth–Death (MTBD) model, with or without contact tracing (CT), and with or without Skyline.

Preprint

treesimulator is described in the S1 Appendix of the following article:

Anna Zhukova, Olivier Gascuel. Accounting for contact tracing 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 (In practice, we ask the user to provide an m x m matrix for MTBD transition rates. μ_ii must be 0 for all i)

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

Using these probabilities, one can calculate the equilibrium frequencies $π_i$ of the model's states, where 0 ≤ π_i ≤ 1 and π₁ + ... + π_m = 1.

The MTBD model has the following epidemiological parameters:

• R_i = (Σ_1≤j≤m λ_ij)/ψ_i -- reproduction number of state i
• d_i = 1/(Σ_{1≤j≤m;i≠j}μ_ij + ψ_i) -- exit time from state i

Contact Tracing (CT)

Contact tracing extension was introduced by Zhukova & Gascuel [medRxiv 2024]. It adds two parameters to the initial MTBD model:

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

and a meta-parameter κ, which defines how many most recent contacts can be notified by each index case. Each is notified independently, with the probability υ.

CT extension adds m notified contact states (1_C, ..., m_C) to its original MTBD model, where i_C is a notified version of the state i. Transition rates for a notified state i_C are analogous to those of i (μ_iCjC=μ_ij), while transitions from non-notified states to notified ones and vice versa are not allowed (μ_iCj=μ_jiC=0). Transmission rates for a notified state i_C are the same as those of i (λ_iCj=λ_ij), where the recipients are always in a non-notified state (λ_iCjC=λ_ijC=0). The removal rate for a notified state i_C is φ.

For CT models exit times from notified contact states i_C is calculated as:

• d_iC = 1/(Σ_{1≤j≤m;i≠j}μ_iCjC + φ)

Skyline

Skyline was introduced by Stadler et al. [PNAS 2013] and extended to MTBD by Kühnert et al. [MBE 2016]. It enables piece-wise constant parameter value changes. To use a skyline with k models, one needs to specify k-1 model change times t₁, ..., _k-1, and k sets of model parameters (see above). At time 0 the simulation starts with model 1, it switches to models 2 at time t₁, etc. All the models in the Skyline must have the same states. CT-related parameters can also change at skyline changing times, in that case if some skyline intervals do not have CT, υ=0 and any value for φ must be specified for them. The same κ value is shared among all the skyline -CT models.

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

BD

1 state:

• I (infectious)

3 parameters:

• λ = λ_I -- transmission rate
• ψ = ψ_I -- removal rate
• p = p_I -- sampling probability upon removal

Epidemiological parameters:

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

BD-CT(κ)

2 states:

• I, infectious
• I_C, notified infectious contact

5 parameters:

• λ = λ_I = λ_IC -- transmission rate
• ψ = ψ_I -- removal rate of I
• p = p_I -- sampling probability upon removal of I
• υ -- probability to notify contacts upon sampling
• φ -- notified contact removal and sampling rate of I_C: φ >> ψ

BDEI

2 states:

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

4 parameters:

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

BDEI-specific epidemiological parameter:

• d_E = 1/μ -- incubation period

BDEI-CT(κ)

4 states:

• E, exposed, i.e. infected but not yet infectious
• I, infectious
• E_C, notified exposed contact
• I_C, notified infectious contact

6 parameters:

• μ = μ_EI = μ_ECIC-- transition rate from an exposed state (notified or not) to the corresponding infectious state (becoming infectious)
• λ = λ_IE = λ_ICE -- transmission rate from an infectious state (notified or not) to E
• ψ = ψ_I -- removal rate of I
• p = p_I -- sampling probability upon removal
• υ -- probability to notify contacts upon sampling
• φ -- notified contact removal and sampling rate: φ >> ψ

BDSS

2 states:

• I, standard infectious individual (a.k.a. normal spreader)
• S, superspreader

5(+1) parameters:

• λ_nn = λ_II -- transmission rate from I to I

• λ_ns = λ_IS -- transmission rate from I to S

• λ_sn = λ_SI -- transmission rate from S to I

• λ_ss = λ_SS -- transmission rate from S to S

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

• ψ = ψ_I = ψ_S -- removal rate

• p = p_I = p_S -- sampling probability upon removal

BDSS-specific epidemiological parameters:

• X_S=λ_ss/λ_ns=λ_sn/λ_nn -- super-spreading transmission ratio
• f_S=λ_ss/(λ_sn + λ_ss) -- super-spreading fraction

BDSS-CT(κ)

4 states:

• I, standard infectious individual (a.k.a. normal spreader)
• S, superspreader
• I_C, notified normal spreader
• S_C, notified superspreader

7(+1) parameters:

• λ_nn = λ_II = λ_ICI -- transmission rate from a normal spreader (notified or not) to I

• λ_ns = λ_IS = λ_ICS -- transmission rate from a normal spreader (notified or not) to S

• λ_sn = λ_SI = λ_SCI -- transmission rate from a superspreader (notified or not) to I

• λ_ss = λ_SS = λ_SCS -- transmission rate from a superspreader (notified or not) to S

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

• ψ = ψ_I = ψ_S -- removal rate of a non-notified state

• p = p_I = p_S -- sampling probability upon removal of a non-notified state

• υ -- probability to notify contacts upon sampling

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

BDEISS

3 states:

• E, exposed, i.e. infected but not yet infectious
• I, standard infectious individual (a.k.a. normal spreader)
• S, infectious superspreader

6 parameters:

• μ_n = μ_EI -- transition rate from E to I (becoming infectious for normal spreaders)
• μ_s = μ_ES -- transition rate from E to S (becoming infectious for superspreaders)
• λ_n = λ_IE -- transmission rate from I to E
• λ_s = λ_SE -- transmission rate from S to E
• ψ = ψ_I = ψ_S -- removal rate of I and of S (the same)
• p = p_I = p_S -- sampling probability upon removal (the same for I and S)

BDEISS-specific epidemiological parameters:

• X_S = λ_s/λ_n -- super-spreading transmission ratio
• f_S = μ_s/(μ_n + μ_s) -- super-spreading fraction (among infectious individuals)
• d_E = 1/(μ_n + μ_s) -- incubation period

BDEISS-CT(κ)

6 states:

• E, exposed, i.e. infected but not yet infectious
• N, standard infectious individual (a.k.a. normal spreader)
• S, superspreader
• E_C, notified exposed individual
• I_C, notified normal spreader
• S_C, notified superspreader

8 parameters:

• μ_n = μ_EI = μ_ECIC -- transition rate from an exposed state (notified or not) to the corresponding normal infectious state (becoming infectious for normal spreaders)
• μ_s = μ_ES = μ_ECSC -- transition rate from an exposed state (notified or not) to the corresponding superspreader infectious state (becoming infectious for superspreaders)
• λ_n = λ_IE = λ_ICE -- transmission rate from I or I_C to E
• λ_s = λ_SE = λ_SCE -- transmission rate from S or S_C to E
• ψ = ψ_I = ψ_S -- removal rate of I and of S (the same)
• p = p_I = p_S -- sampling probability upon removal of I or S
• υ -- probability to notify contacts upon sampling
• φ -- notified contact removal and sampling rate: φ >> ψ

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, BD-CT(κ) and BD-CT(κ)-Skyline

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-CT(1) model, with λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify only the most recent contact 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_contacts 1 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BD-CT(1)-Skyline model with two time intervals, with λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0 between t=0 and t=3, and λ=1, ψ=0.25, p=0.75, φ=2.5, υ=0.2 starting at t=3, and allowing to notify only the most recent contact 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 1 --psi 0.25 0.25 --p 0.5 0.75 \
--phi 2.5 2.5 --upsilon 0 0.2 --max_notified_contacts 1 \
--skyline_times 3 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bd --help

BDEI, BDEI-CT(κ) and BDEI-CT(κ)-Skyline

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-CT(2) model, with μ=1, λ=0.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify last two contacts 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_contacts 2 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BDEI-CT(2)-Skyline model with three time intervals, with μ=1, λ=0.5, ψ=0.25, p=0.2, φ=2.5, υ=0.2, between t=0 and t=2, with μ=1, λ=0.5, ψ=0.3, p=0.3, φ=2.5, υ=0.3, between t=2 and t=3, and μ=1, λ=0.5, ψ=0.5, p=0.5, φ=5, υ=0.3 starting at t=3, and allowing to notify last two contacts 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 1 1 --la 0.5 0.5 0.5 --psi 0.25 0.3 0.5 --p 0.2 0.3 0.5 \
--phi 2.5 2.5 5 --upsilon 0.2 0.3 0.3 --max_notified_contacts 2 \
--skyline_times 2 3 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bdei --help

BDSS, BDSS-CT(κ) and BDSS-CT(κ)-Skyline

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-CT(3) 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 last three contacts 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_contacts 3 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BDSS-CT(3)-Skyline model with two time intervals, with λ_nn=0.1, λ_ns=0.3, λ_sn=0.5, λ_ss=1.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2 between t=0 and t=2, and λ_nn=0.1, λ_ns=0.3, λ_sn=1, λ_ss=3, ψ=0.25, p=0.5, φ=5, υ=0.5 starting at t=2, and allowing to notify last three contacts 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 0.1 --la_ns 0.3 0.3 --la_sn 0.5 1 --la_ss 1.5 3 --psi 0.25 0.25 --p 0.5 0.5 \
--phi 2.5 5 --upsilon 0.2 0.5 --max_notified_contacts 3 \
--skyline_times 2 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bdss --help

BDEISS, BDEISS-CT(κ) and BDEISS-CT(κ)-Skyline

The following command simulates a tree with 200-500 tips under the BDEISS model, with μ_n=0.1, μ_s=0.3, λ_n=0.5, λ_s=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_bdeiss --min_tips 200 --max_tips 500 \
--mu_n 0.1 --mu_s 0.3 --la_n 0.5 --la_s 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 BDEISS-CT(1) model, with μ_n=0.1, μ_s=0.3, λ_n=0.5, λ_s=1.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2, and allowing to notify the last contact 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_bdeiss --min_tips 200 --max_tips 500 \
--mu_n 0.1 --mu_s 0.3 --la_n 0.5 --la_s 1.5 --psi 0.25 --p 0.5 \
--phi 2.5 --upsilon 0.2 --max_notified_contacts 1 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under the BDSS-CT(1)-Skyline model with two time intervals, with μ_n=0.1, μ_s=0.3, λ_n=0.5, λ_s=1.5, ψ=0.25, p=0.5, φ=2.5, υ=0.2 between t=0 and t=2, and μ_n=0.1, μ_s=0.3, λ_n=0.5, λ_s=1.5, ψ=0.25, p=0.2, φ=2.5, υ=0.5 starting at t=2, and allowing to notify the last contact 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_bdeiss --min_tips 200 --max_tips 500 \
--mu_n 0.1 0.1 --mu_s 0.3 0.3 --la_n 0.5 0.5 --la_s 1.5 1.5 --psi 0.25 0.25 --p 0.5 0.2 \
--phi 2.5 2.5 --upsilon 0.2 0.5 --max_notified_contacts 1 \
--skyline_times 2 \
--nwk tree.nwk --log params.csv

To see detailed options, run:

generate_bdss --help

User-defined MTBD, MTBD-CT(κ) and MTBD-CT(κ)-Skyline models

The following command simulates a tree with 200-500 tips under a generic MTBD model, with two states A and B, with μ_ab=0.6, μ_ba=0.7 (note that μ_aa=μ_bb=0 as only transitions between different states are possible!), λ_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 0.6 0.7 0 \
--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-CT(1) model, with two states A and B, with μ_ab=0.6, μ_ba=0.7, λ_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 contact 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 0.6 0.7 0 \
--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_contacts 1 \
--nwk tree.nwk --log params.csv

The following command simulates a tree with 200-500 tips under a generic MTBD-CT(1)-Skyline model, with two states A and B, with μ_ab=0.6, μ_ba=0.7, λ_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 between t=0 and t=8, and μ_ab=1.6, μ_ba=1.7, λ_aa=1.1, λ_ab=1.2, λ_ba=1.3, λ_bb=1.4, ψ_a=1.05, ψ_b=1.08, p=_a0.1, p=_b0.6, φ=3.5, υ=0.4 starting at t=8, and allowing to notify only the most recent contact 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 0.6 0.7 0 0 1.6 1.7 0 \
--transmission_rates 0.1 0.2 0.3 0.4 1.1 1.2 1.3 1.4 \
--removal_rates 0.05 0.08 1.05 1.08 \
--sampling_probabilities 0.15 0.65 0.1 0.6 \
--phi 2.5 3.5 --upsilon 0.2 0.4 --max_notified_contacts 1 \
--skyline_times 8 \
--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, BirthDeathExposedInfectiousWithSuperSpreadingModel, CTModel

# 1. BD, BD-CT(1) and BD-CT(1)-Skyline
## BD model
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')
## Adding -CT to the model above
bdct_model = CTModel(model=bd_model, upsilon=0.2, phi=2.5)
[bdct_tree], _, _ = generate([bdct_model], min_tips=200, max_tips=500, max_notified_contacts=1)
save_forest([bdct_tree], 'BDCT_tree.nwk')
## BD-CT(1)-Skyline models
bdct_model_1 = CTModel(BirthDeathModel(p=0.5, la=0.5, psi=0.25),
                       upsilon=0, phi=2.5)
bdct_model_2 = CTModel(BirthDeathModel(p=0.75, la=1, psi=0.25),
                       upsilon=0.2, phi=2.5)
[bdct_skyline_tree], _, _ = generate([bdct_model_1, bdct_model_2], skyline_times=[3],
                                     min_tips=200, max_tips=500, max_notified_contacts=1)
save_forest([bdct_skyline_tree], 'BDCTSkyline_tree.nwk')

# BDEI, BDEI-CT(2) and BDEI-CT(2)-Skyline
## BDEI model
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')
## Adding -CT to the model above
bdeict_model = CTModel(model=bdei_model, upsilon=0.2, phi=2.5)
[bdeict_tree], _, _ = generate([bdeict_model], min_tips=200, max_tips=500, max_notified_contacts=2)
save_forest([bdeict_tree], 'BDEICT_tree.nwk')
## BDEI-CT(2)-Skyline with three time intervals
bdeict_model_1 = CTModel(model=BirthDeathExposedInfectiousModel(p=0.2, mu=1, la=0.5, psi=0.25), upsilon=0.2,
                         phi=2.5)
bdeict_model_2 = CTModel(model=BirthDeathExposedInfectiousModel(p=0.3, mu=1, la=0.5, psi=0.3), upsilon=0.3, phi=2.5)
bdeict_model_3 = CTModel(model=BirthDeathExposedInfectiousModel(p=0.5, mu=1, la=0.5, psi=0.5), upsilon=0.3, phi=5)
[bdeict_skyline_tree], _, _ = generate([bdeict_model_1, bdeict_model_2, bdeict_model_3], skyline_times=[2, 3],
                                       min_tips=200, max_tips=500, max_notified_contacts=2)
save_forest([bdeict_skyline_tree], 'BDEICTSkyline_tree.nwk')

# BDSS, BDSS-CT(3) and BDSS-CT(3)-Skyline
## BDSS model
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')
## Adding -CT to the model above
bdssct_model = CTModel(model=bdss_model, upsilon=0.2, phi=2.5)
[bdssct_tree], _, _ = generate([bdssct_model], min_tips=200, max_tips=500, max_notified_contacts=3)
save_forest([bdssct_tree], 'BDSSCT_tree.nwk')
## BDSS-CT(3)-Skyline with two time intervals, using the model above for the first interval
bdssct_model_2 = CTModel(
  model=BirthDeathWithSuperSpreadingModel(p=0.5, la_nn=0.1, la_ns=0.3, la_sn=1, la_ss=3, psi=0.25),
  upsilon=0.5, phi=5)
[bdssct_skyline_tree], _, _ = generate([bdssct_model, bdssct_model_2], skyline_times=[2], min_tips=200, max_tips=500,
                                       max_notified_contacts=3)
save_forest([bdssct_skyline_tree], 'BDSSCTSkyline_tree.nwk')

# BDEISS, BDEISS-CT(1) and BDEISS-CT(1)-Skyline
## BDEISS model
bdeiss_model = BirthDeathExposedInfectiousWithSuperSpreadingModel(p=0.5, mu_n=0.1, mu_s=0.3, la_n=0.5, la_s=1.5,
                                                                  psi=0.25)
print(bdeiss_model.get_epidemiological_parameters())
[bdeiss_tree], _, _ = generate([bdeiss_model], min_tips=200, max_tips=500)
save_forest([bdeiss_tree], 'BDEISS_tree.nwk')
## Adding -CT to the model above
bdeissct_model = CTModel(model=bdeiss_model, upsilon=0.2, phi=2.5)
[bdeissct_tree], _, _ = generate([bdeissct_model], min_tips=200, max_tips=500, max_notified_contacts=1)
save_forest([bdeissct_tree], 'BDEISSCT_tree.nwk')
## BDEISS-CT(1)-Skyline with two time intervals, using the model above for the first interval
bdeissct_model_2 = CTModel(
  model=BirthDeathExposedInfectiousWithSuperSpreadingModel(p=0.2, mu_n=0.1, mu_s=0.3, la_n=0.5, la_s=1.5, psi=0.25),
  upsilon=0.5, phi=5)
[bdeissct_skyline_tree], _, _ = generate([bdeissct_model, bdeissct_model_2], skyline_times=[2], min_tips=200,
                                         max_tips=500,
                                         max_notified_contacts=1)
save_forest([bdeissct_skyline_tree], 'BDEISSCTSkyline_tree.nwk')

# MTBD, MTBD-CT(1) and MTBD-CT(1)-Skyline
## MTBD model with two states: A and B
mtbd_model = Model(states=['A', 'B'], transition_rates=[[0, 0.6], [0.7, 0]],
                   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')
## Adding -CT to the model above
mtbdct_model = CTModel(model=mtbd_model, upsilon=0.2, phi=2.5)
[mtbdct_tree], _, _ = generate([mtbdct_model], min_tips=200, max_tips=500, max_notified_contacts=1)
save_forest([mtbdct_tree], 'MTBDCT_tree.nwk')
## MTBD-CT(1)-Skyline with two time intervals, using the model above for the first interval
mtbdct_model_2 = CTModel(model=Model(states=['A', 'B'], transition_rates=[[0, 1.6], [1.7, 0]],
                                     transmission_rates=[[1.1, 1.2], [1.3, 1.4]],
                                     removal_rates=[1.05, 1.08], ps=[0.1, 0.6]),
                         upsilon=0.4, phi=3.5)
[mtbdct_skyline_tree], _, _ = generate([mtbdct_model, mtbdct_model_2], skyline_times=[8],
                                       min_tips=200, max_tips=500, max_notified_contacts=1)
save_forest([mtbdct_skyline_tree], 'MTBDCTSkyline_tree.nwk')

Run with apptainer

Once apptainer is installed, run the following command:

apptainer run docker://evolbioinfo/treesimulator

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.