An Introduction to the`TiPS`package

Gonché Danesh, Emma Saulnier, Olivier Gascuel, MarcChoisy, Samuel Alizon

library(TiPS)library(ape)

Simulating trajectories

Building the simulator

We use the classical SIR epidemiological model to illustrate thefunctionning ofTiPS.

This SIR model can be described by a system of reactions such as

\[ S+I \xrightarrow[]{\beta S I} I+I\]

and

\[ I \xrightarrow[]{\gamma I} R\]

or by a system of differential equations such as

\[\frac{dS}{dt} = -\beta S I\]

and

\[\frac{dI}{dt} = \beta S I - \gamma I\]

In R, withTiPS, this model can either be written as

reactions <- c("S [beta*S*I] -> I",               "I [gamma*I] -> R")

Let’s now build the simulator:

We then build the simulator that will allow us to run multipletrajectories:

sir_simu <- build_simulator(reactions)

This typically takes 10-15’’ as it involves compilation.

Defining simulations parameters

To run numerical simulations, we define the initial values of thestate variables,

initialStates <- c(I = 1, S = 9999, R = 0)

the time range of the simulations,

time <- c(0, 20)

and the parameters values

theta <- list(gamma = 1, beta = 2e-4)

For the$\tau$-leap and mixedalgorithms, a time step is also required:

dT <- 0.001

Running simulations

In some simulations, the population size of a deme compartment may bezero before the upper time limit is reached, because of stochasticity orparameter values. In this case, the simulation is considered to havefailed and is halted.

To bypass these failures, we can define the following wrapper:

safe_run <- function(f, ...) {  out <- list()  while(! length(out)) {out <- f(...)}  out}

A safe version of our simulatorsir_simu() is then:

safe_sir_simu <- function(...) safe_run(sir_simu, ...)

Direct method

A trajectory using Gillespie’s direct method is obtained by

traj_dm <- safe_sir_simu(  paramValues = theta,  initialStates = initialStates,  times = time,  method = "exact")

The output consists of a named list containing the reactions of themodel (with$reactions), the parameter values (with$values), the time range (with$times), thealgorithm used to simulate (with$method), the time-step incase the algorithm is$\tau$-leap orthe mixed algorithm (with$dT), the random seed that wasgenerated and was used to simulate trajectories (with$seed) and finally the simulated trajectory (with$traj) :

names(traj_dm)#> [1] "reactions" "values"    "times"     "method"    "tau"       "seed"     #> [7] "traj"

The simulated trajectory is also a named list, where each simulatedreaction event is recorded$Reaction, along with the timeat which it occured$Time, the number of times it occured$Nrep (if$\tau$-leap ormixed algorithm chosen), and the size of each compartment througt time,here$I$R$S.

head(traj_dm$traj)#>        Time          Reaction Nrep    S I R#> 1 0.0000000              init    1 9999 1 0#> 2 0.5957257 S [beta*S*I] -> I    1 9998 2 0#> 3 0.7271658 S [beta*S*I] -> I    1 9997 3 0#> 4 0.8464166 S [beta*S*I] -> I    1 9996 4 0#> 5 0.8805381  I [gamma*I] -> R    1 9996 3 1#> 6 0.9309763  I [gamma*I] -> R    1 9996 2 2

The trajectory can readily be plotted using theplot()function:

plot(traj_dm)

You can also specify the seed to generate reproducible results withthe parameterseed. By default, its value is null (set to0), in which case the seed is randomly generated. Let’s fix the seed to166 :

traj_dm <- sir_simu(  paramValues = theta,  initialStates = initialStates,  times = time,  method = "exact",  seed = 166)

When running this multiple times, you will always get the exact sameresults :

head(traj_dm$traj)#>         Time          Reaction Nrep    S I R#> 1 0.00000000              init    1 9999 1 0#> 2 0.06279117 S [beta*S*I] -> I    1 9998 2 0#> 3 0.10048223 S [beta*S*I] -> I    1 9997 3 0#> 4 0.13123844  I [gamma*I] -> R    1 9997 2 1#> 5 0.22783501  I [gamma*I] -> R    1 9997 1 2#> 6 0.35976195 S [beta*S*I] -> I    1 9996 2 2

$\tau$-leap method

The default mode of the simulator is the$\tau$-leap method. Themethodargument must be specified and fixed toapproximate. Thetime-steptau is by default set to 0.05. Let’s fix itsvalue to 0.009:

traj_tl <- safe_sir_simu(  paramValues = theta,    initialStates = initialStates,    times = time,  method = "approximate",    tau = 0.009)

We obtain the same type of output as with the direct method:

head(traj_tl$traj)#>    Time          Reaction Nrep    S I R#> 1 0.000              init    1 9999 1 0#> 2 0.603 S [beta*S*I] -> I    1 9998 2 0#> 3 0.720 S [beta*S*I] -> I    1 9997 3 0#> 4 0.981 S [beta*S*I] -> I    1 9996 4 0#> 5 1.008 S [beta*S*I] -> I    1 9995 5 0#> 6 1.125  I [gamma*I] -> R    1 9995 4 1

The trajectory can also be plotted:

plot(traj_tl)

Mixed method

To run simulations with the Mixed Simulation Algorithm (MSA)(basically switching between the direct method and the$\tau$-leap method depending on the numberof reactions occurring per unit of time), themethodargument must be specified and fixed tomixed. Thetime-steptau is by default set to 0.05. Let’s fix itsvalue to 0.009 :

traj_mm <- safe_sir_simu(    paramValues = theta,    initialStates = initialStates,    times = time,    method = "mixed",    tau = 0.009)

Outputs are similar to the other methods and the trajectory can alsobe plotted :

names(traj_mm)#> [1] "reactions" "values"    "times"     "method"    "tau"       "seed"     #> [7] "traj"head(traj_mm$traj)#>        Time          Reaction Nrep    S I R#> 1 0.0000000              init    1 9999 1 0#> 2 0.4163948 S [beta*S*I] -> I    1 9998 2 0#> 3 0.4818173 S [beta*S*I] -> I    1 9997 3 0#> 4 0.4982856 S [beta*S*I] -> I    1 9996 4 0#> 5 0.5009927 S [beta*S*I] -> I    1 9995 5 0#> 6 0.5110239 S [beta*S*I] -> I    1 9994 6 0plot(traj_mm)

The MSA is a new algorithm we introduce that switches from the directmethod (GDA) to the tau-leap algorithm (GTA) if overniterations (optionmsaIt, by defaultmsaIt = 10) the time until the next event$\delta_t$ is below a user-definedthreshold (optionmsaTau, by defaultmsaTau = tau/10), and from GTA to GDA if the total numberof realised events is lower than the number of possible events. Whensimulating trajectories with the mixed algorithm, it is possible tomodify these parameters such as :

traj_mm <- safe_sir_simu(    paramValues = theta,    initialStates = initialStates,    times = time,    method = "mixed",    tau = 0.009,  msaTau = 0.0001,  msaIt = 20)

Trajectory in output file

It is possible to specify to write the trajectory directly in a tabdelimited output file with the optionoutFile. By default,the trajectory is as presented previously andoutFile isempty.

traj_dm <- sir_simu(  paramValues = theta,  initialStates = initialStates,  times = time,  method = "exact",  seed = 166,  outFile = "sir_traj.txt")

In that case,traj_dm$traj isNULL and theoutput file name is indicated intraj_dm$outFile. It ispossible to visualize the simulated trajectory by reading with theread.table R function :

trajectory <- read.table(file = "sir_traj.txt", header = TRUE, sep = "\t", stringsAsFactors = F)

Simulating phylogenies

A great advantage ofTiPS, besides its compulationalefficiency, is that it can generate phylogenies from the populationdynamics trajectories using a coalescent approach.

For this, we may need a vector of sampling dates.

For the SIR example, these are stored here:

dates <- system.file("extdata", "SIR-dates.txt", package = "TiPS")

Thesimulate_tree function simulates a phylogeny from atrajectory object and using a set of sampling dates:

sir_tree <- simulate_tree(  simuResults = traj_dm,  dates = dates,  deme = c("I"), # the type of individuals that contribute to the phylogeny  sampled = c(I = 1), # the type of individuals that are sampled and their proportion of sampling  root = "I", # type of individual at the root of the tree  isFullTrajectory = FALSE, # deads do not generate leaves  nTrials = 5,  addInfos = FALSE, # additional info for each node  seed = 54673) ## to reproduce results (optionnal)

Thesampled option can be used for labelled phylogenies(i.e. with multiple host types) but it requires specifying theproportion of each label type. Theroot option indicatesthe state of the individual initiating the dynnamics. See the nextsection for details.

The full phylogeny can be obtained (therefore neglecting the samplingdates) with the optionisFullTrajectory.

Finally, some runs may fail to simulate a phylogeny from a trajectoryfor stochastic reasons. ThenTrials parameter indicates thenumber of unsuccessful trials allowed before giving up.

The simulated phylogeny can be visualised using:

ape::plot.phylo(sir_tree, cex = .5)

It is possible to write the tree in a output file with the optionoutFile. If a file name is given as input, by default thetree is written in a Newick format (optionformat = nexus).To write the simulated tree in a Nexus format, the optionformat must be specified and fixed tonexus.

sir_tree <- simulate_tree(  simuResults = traj_dm,  dates = dates,  deme = c("I"),  sampled = c(I = 1),  root = "I",  isFullTrajectory = FALSE,  nTrials = 5,  addInfos = FALSE,  outFile = "sir_tree.nexus",  format = "nexus")

Multiple demes

We sometimes have multiple demes, i.e. different types of individualsthat contribute to the pylogeny or that can be sampled (e.g. juvenilesvs. adults or acute vs. chronic infections).

We illustrate this example using an SIR model with two patches(labelled 1 and 2) and migration between these patches (at a rate$\mu$).

Initialising the system

\[\frac{dI_1}{dt} = \beta_1 I_1 - \gamma_1 I_1 - \mu_1 I_1 + \mu_2 I_2 \\\frac{dI_2}{dt} = \beta_2 I_2 - \gamma_2 I_2 - \mu_2 I_2 + \mu_1 I_1 \\\]

The associated reactions are:

reactions <- c("0 [beta1 * I1] -> I1",               "I1 [gamma1 * I1] -> 0",               "I1 [mu1 * I1] -> I2",               "0 [beta2 * I2] -> I2",               "I2 [gamma2 * I2] -> 0",               "I2 [mu2 * I2] -> I1")

We then build the simulator:

bd_simu <- build_simulator(reactions)

The initial state variables values are

initialStates <- c(I1 = 0, I2 = 1)

The time range of the simulation is between 1975 and 2018:

time <- c(1975, 1998, 2018)

the parameters values are

theta <- list(gamma1 = c(0.2, 0.09), gamma2 = 0.1, mu1 = 0.025, mu2 = 0.021, beta1 = c(0.26,0.37), beta2 = 0.414)

and the time step (for the$\tau$-leap and mixed algorithms) is:

dT <- 0.01

A safe version of the simulatorbd_simu() is:

safe_bd_simu <- function(...) safe_run(bd_simu, ...)

Tau-leap trajectory simulation

We perform the simulations using:

trajbd_tl <- safe_bd_simu(  paramValues = theta,  initialStates = initialStates,  times = time,  method = "approximate",  tau = 0.001)

We obtain:

head(trajbd_tl$traj)#>       Time              Reaction Nrep I1 I2#> 1 1975.000                  init    1  0  1#> 2 1977.020  0 [beta2 * I2] -> I2    1  0  2#> 3 1977.547  0 [beta2 * I2] -> I2    1  0  3#> 4 1977.997  0 [beta2 * I2] -> I2    1  0  4#> 5 1978.297   I2 [mu2 * I2] -> I1    1  1  3#> 6 1978.938 I1 [gamma1 * I1] -> 0    1  0  3

Graphically, we get:

plot(trajbd_tl)

Phylogeny simulation

With known sampling dates and known proportion of sampling

Instead of loading a vector, we assume we have 100 samples at 100sampling dates between 2015 and 2018. We can generate the dates vectoras:

dates_bd <- seq(from=2015, to=2018, length.out=100)

We then simulate a phylogeny where 20% of the sampling datescorrespond to the I1 compartment, and 80% to the I2 compartment:

bd_tree <- simulate_tree(  simuResults = trajbd_tl,  dates = dates_bd,  deme = c("I1", "I2"),  sampled = c(I1 = 0.2, I2 = 0.8), # the type of individuals that are sampled and their proportion of sampling  root = "I2", # type of individual at the root of the tree  isFullTrajectory = FALSE, # deads do not generate leaves  nTrials = 3,  addInfos = TRUE) # additional info for each node

This is done using a coaslescence process informed by the trajectory.Therefore, each internal node of the phylogeny corresponds to acoalescence event and is associated with a label stoed in$node.label.

In our two-patches example, there are two types of coalesence: I2individuals creating a new I2 individual, and I1 individuals creating anew I1 individual.

We can plot the phylogeny and color the internal nodes based on thetype of coalescence.

First we generate a vector of colors for the nodes (if we findI2 in the node label we color it in blue, otherwise inred):

inode_cols <- ifelse(grepl(x=bd_tree$node.label,pattern="I2"),"blue","red")

Then we plot the phylogeny:

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, align.tip.label = T)nodelabels(pch=20,col=inode_cols)

With known sampling dates, each assigned to a compartment by theuser

One can give as input, sampling dates assigned to a compartment, inwhich case the optionsampled is not required.

dates_bd <- seq(from=2015, to=2018, length.out=100)dates_bd <- data.frame(Date=sample(dates_bd),Comp=c(rep("I1",20),rep("I2",80)))head(dates_bd)#>       Date Comp#> 1 2015.091   I1#> 2 2017.970   I1#> 3 2016.061   I1#> 4 2015.303   I1#> 5 2015.455   I1#> 6 2016.697   I1

Now let’s simulate a phylogeny with sampling dates assigned to acompartment by the user.

bd_tree <- simulate_tree(  simuResults = trajbd_tl,  dates = dates_bd,  deme = c("I1", "I2"),  root = "I2", # type of individual at the root of the tree  nTrials = 3,  addInfos = TRUE) # additional info for each node

We can plot the phylogeny and color the external nodes given thecompartment.

tips_cols <- ifelse(grepl(x=bd_tree$tip.label,pattern="I2"),"blue","red")

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)tiplabels(pch=20,col=tips_cols)

With only known sampling dates

In the case where the user has no information on the samplingproportions or the assignment of sampling dates on any compartment, thealgorithm will randomly assign each sampling date to a compartment. Theuser gives as input sampling dates:

dates_bd <- seq(from=2015, to=2018, length.out=100)

Now let’s simulate a phylogeny with sampling dates and no informationabout the sampling schemes :

bd_tree <- simulate_tree(  simuResults = trajbd_tl,  dates = dates_bd,  sampled = c("I1","I2"),  deme = c("I1", "I2"),  root = "I2", # type of individual at the root of the tree  nTrials = 10,  addInfos = TRUE) # additional info for each node

ape::plot.phylo(bd_tree, root.edge = T, no.margin = F, show.tip.label = F)

References

For further details see Danesh G, Saulnier E, Gascuel O, ChoisyM, Alizon S. 2020. Simulating trajectories and phylogenies frompopulation dynamics models with TiPS.bioRxiv,2020.11.09.373795. DOI:10.1101/2020.11.09.373795.
This work was supported by a doctoral grant from la Fondationpour la Recherche Medicale (FRM grant number ECO20170637560) to GoncheDanesh.

Movatterモバイル変換

An Introduction to the`TiPS`package

Gonché Danesh, Emma Saulnier, Olivier Gascuel, MarcChoisy, Samuel Alizon

Simulating trajectories

Building the simulator

Defining simulations parameters

Running simulations

Direct method

\(\tau\)-leap method

Mixed method

Trajectory in output file

Simulating phylogenies

Multiple demes

Initialising the system

Tau-leap trajectory simulation

Phylogeny simulation

With known sampling dates and known proportion of sampling

With known sampling dates, each assigned to a compartment by theuser

With only known sampling dates

References

Movatterモバイル変換

An Introduction to theTiPSpackage

Gonché Danesh, Emma Saulnier, Olivier Gascuel, MarcChoisy, Samuel Alizon

Simulating trajectories

Building the simulator

Defining simulations parameters

Running simulations

Direct method

\(\tau\)-leap method

Mixed method

Trajectory in output file

Simulating phylogenies

Multiple demes

Initialising the system

Tau-leap trajectory simulation

Phylogeny simulation

With known sampling dates and known proportion of sampling

With known sampling dates, each assigned to a compartment by theuser

With only known sampling dates

References

An Introduction to the`TiPS`package