Gonché Danesh
email: gonche.danesh@gmail.com

Emma Saulnier

Olivier Gascuel

Marc Choisy

Samuel Alizon

TiPS: rapidly simulating trajectories and phylogenies from compartmental models

Keywords: R package, stochastic simulations, phylogenies, compartmental models, population dynamics

We introduce TiPS, an R package to generate trajectories and phylogenetic trees associated with a compartmental model. Trajectories are simulated using Gillespie's exact or approximate stochastic simulation algorithm, or a newly proposed mixed version of the two. Phylogenetic trees are simulated from a trajectory under a backward-in-time approach (i.e. coalescent). TiPS is based on the Rcpp package and therefore combines the flexibility of R for model definition and the speed of C++ for simulation execution. The model is defined in R with a set of reactions that allow us to capture heterogeneity in life cycles or any kind of population structure. TiPS converts the model into C++ code and compiles it into a simulator, which is interfaced in R via a function. Furthermore, the package allows one to define time periods in which the model parameters can take different values. This package, available on the CRAN at https://cran.r-project.org/ package=TiPS, is particularly well suited for population genetics and phylodynamic studies that require generating a large number of phylogenies used for population dynamics studies.

Introduction

Stochastic simulations of population dynamics are routinely used in ecology and epidemiology to generate trajectories (i.e. time series of population sizes) and genealogies that capture the relatedness between individuals [START_REF] Otto | A biologist's guide to mathematical modeling in ecology and evolution[END_REF]Keeling and Rohani, 2008;[START_REF] Lenormand | Stochasticity in evolution[END_REF]. The increasing amount of genetic data is fuelling interest in linking population dynamics and genealogies because the former can leave footprints in individuals' genomes [START_REF] Grenfell | Unifying the Epidemiological and Evolutionary Dynamics of Pathogens[END_REF][START_REF] Volz | Viral phylodynamics[END_REF][START_REF] Frost | Eight challenges in phylodynamic inference[END_REF]. Such phylodynamics studies involve computer-intensive methods that can require the simulation of many trajectories and genealogies [START_REF] Ratmann | Phylodynamic inference and model assessment with approximate bayesian computation: influenza as a case study[END_REF][START_REF] Gascuel | How Ecology and Landscape Dynamics Shape Phylogenetic Trees[END_REF][START_REF] Saulnier | Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study[END_REF].

A common method to simulate population dynamics trajectories is Gillespie's exact stochastic simulation algorithm (SSA) [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], which is rooted in probability theory [START_REF] Kurtz | Solutions of ordinary differential equations as limits of pure jump markov processes[END_REF]. In the R software environment, it is implemented in packages such as GillespieSSA [START_REF] Pineda-Krch | GillespieSSA: implementing the stochastic simulation algorithm in R[END_REF], adaptivetau [START_REF] Johnson | adaptivetau: efficient stochastic simulations in R. R Package Version[END_REF], or epimdr [START_REF] Bjørnstad | Epidemics. Models and Data using R[END_REF]. The computational speed of these software packages is facilitated by the fact that they do not keep track of the history of the process, that is the trajectory. Conversely, in this same environment, geiger [START_REF] Pennell | geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees[END_REF], phytools [START_REF] Revell | phytools: An R package for phylogenetic comparative biology (and other things)[END_REF], ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF], and TreeSim can simulate phylogenies. However, the underlying model is often very simple, e.g. birth-death model. Also in R, some packages such as nosoi [START_REF] Lequime | nosoi: A stochastic agent-based transmission chain simulation framework in r[END_REF] allow the user to implement detailed agent-based models but their computational time is slow and the outputs are difficult to compare to classical compartmental models based on differential equations.

However, some software packages can simulate both trajectories and genealogies. In R, rcolgem, which was updated to phydynR [START_REF] Volz | Complex population dynamics and the coalescent under neutrality[END_REF], combines an Euler-Maruyama integration and the structured coalescent to allow the user to rapidly simulate phylogenies from any compartmental model. This will be our main reference in the following in terms of accuracy and computational speed. Another exception is the software package MASTER [START_REF] Vaughan | A stochastic simulator of birth-death master equations with application to phylodynamics[END_REF] in the BEAST2 platform [START_REF] Bouckaert | BEAST 2: A Software Platform for Bayesian Evolutionary Analysis[END_REF]. Although MASTER is a useful tool to simulate both time series and genealogies, the specification of the model of interest in the XML language is not as intuitive as the packages of the R environment. Furthermore, although MASTER is one of the fastest options to simulate a few phylogenies (because it does not need to compile the code), its execution time quickly becomes limiting when simulating thousands (or millions) of phylogenies.

We introduce TiPS, a flexible and easy-to-use R package to rapidly simulate population trajectories and phylogenies using a backwards-in-time, i.e. coalescent, process with either pre-defined sampling dates or a stochastic sampling scheme. We also introduce a new approximate version of the Gillespie algorithm to increase the calculation speed. A brief benchmarking analysis shows that TiPS is faster than adaptivetau to simulate trajectories, especially for large populations (Fig. 3a). It is also at least one order of magnitude faster than phydynR to simulate phylogenies (Fig. 3b).

Methods

Package overview

TiPS generates two types of stochastic simulation output: population dynamics trajectories and phylogenies. These are obtained using a continuous-time individual-based model defined in R as a system of reactions. The model is first transcribed in C++ and then compiled, before being linked back to a simulation function in R thanks to the Rcpp package [START_REF] Eddelbuettel | Rcpp: Seamless R and C++ Integration[END_REF]. The general structure of the pipeline is illustrated in Figure 1.

The TiPS package is available on the CRAN at https://cran.r-project.org/web/packages/ TiPS/index.html and is maintained on GitLab at https://gitlab.in2p3.fr/ete/tips/. (Keeling and Rohani, 2008). The functions of the R package are in blue. The simulator of trajectories, which is built as a function, is in orange. The variable traj, in red, is the output trajectory of class simutraj. It can be plotted using our plot method. TiPS used the simulated trajectory and 100 sampling dates that we generated (variable sampling_dates in green) to simulate a sampled phylogeny. The output simulated phylogeny is a phylo class object from the ape R package [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF], which can be used for plotting.

Model description

TiPS simulates trajectories from a user-specified compartmental model. These models divide the population (animals, cells, etc) into distinct categories (geographic, clinical state, etc) or so-called compartments in which the sub-population behaves uniformly. In these models, the population may progress between the different compartments.

Here, we illustrate the use of TiPS with the SIR epidemiological model, where individuals can have three clinical states: susceptible (S), infected (I), and removed (R) (Keeling and Rohani, 2008). The model can be captured with a system of two individual-based reactions:

S βSI --→ I (1a) I γI -→ R (1b)
where β and γ are the transmission and recovery rates. The rate of occurrence of each reaction is indicated above the transition arrow and the corresponding population-based system of ordinary differential equations (ODE) is shown in the top-right panel of Figure 1.

A simulating function is generated from the individual-based reactions of the model of interest. These are entered by the user as a string vector (top left box of Figure 1).

Simulating trajectories

The simulation function takes as arguments a named numeric vector that contains the initial number of individuals in each compartment, a named list with the parameter values, a vector of the time limits of the simulation, and the type of algorithm to use for the simulation. Users can also enter a vector of breakpoints, which allows parameter values to vary over time. These breakpoints are indicated in the time limits vector, and the corresponding parameter values are ordered chronologically in the named list of parameter values.

Three simulation algorithms are implemented in the Rcpp simulating function:

1. Gillespie's Direct Algorithm (GDA, the default option) [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF]) is an exact algorithm that simulates the time until the next event δ t by assuming that waiting times are exponentially distributed. A limitation is that its computational complexity scales linearly with the number of simulated events.

2. Gillespie's Tau-Leap Algorithm (GTA) [START_REF] Gillespie | Approximate accelerated stochastic simulation of chemically reacting systems[END_REF]) is an approximate algorithm that introduces a fixed time step τ during which the number of events of each type is assumed to be Poisson distributed. This algorithm is limited in terms of computation time if the time step is small compared to the rate at which events occur.

3. The Mixed Simulation Algorithm (MSA) is a new algorithm that switches from GDA to GTA depending on the respective values of δ t and τ . The algorithm switches from GDA to GTA if 10 successive estimations of δ t are shorter than τ /10. The algorithm switches back to GDA if the total number of realised events is less than the number of possible events. The MSA shares similarities with the slow-scale stochastic simulation algorithm [START_REF] Cao | The slow-scale stochastic simulation algorithm[END_REF] or the adaptive explicit-implicit tau-leaping method for an optimised tau-leap selection [START_REF] Cao | Adaptive explicit-implicit tau-leaping method with automatic tau selection[END_REF].

The output of a trajectory simulation is a named list containing the simulated trajectory and some details about the model and the simulation, such as the reactions of the model, the parameter values, the time range, the algorithm used to perform simulations, the time step in case the algorithm is the GTA or the MSA, and, for reproducibility purposes, the random seed used to initialise a pseudorandom number generator. If specified by the user, instead of storing the trajectory as a Rcpp data frame, TiPS can write the simulated trajectory (i.e. at each time step, the time, the size of each compartment, the reaction and the number of reactions simulated) directly into an tab delimited output file.

Simulating phylogenies

A phylogeny is the representation of the evolutionary history and relationships between genes, organisms, or groups of organisms. The root of the tree represents the ancestor of all lineages and the leaves represent the most recent descendants of that ancestor. TiPS simulates binary phylogenies rooted in time.

Infectious disease transmission models

In the context of infectious diseases, under models such as the SIR model, when simulating phylogenies, TiPS traces back the epidemiological history of a sampled pathogen infection, where a forward-in-time transmission event is represented as a coalescence between two lineages under a backward-in-time process, and an end of infection (e.g. caused by death, treatment, or sampling) is represented by a leaf in the phylogeny.

In the SIR model, the pathogen is only present in infected individuals that belong to the I compartment, referred to as 'deme' individuals in the deme compartment. The other compartments (i.e. S and R) are referred to as 'non deme' compartments.

There are four types of reactions in epidemiological models in TiPS:

-transmission: this reaction corresponds to the generation of a new deme individual (Equation 1a in the SIR model);

-removal: this corresponds to the removal of a deme individual from its deme compartment or the displacement of a deme individual from its deme compartment to a non-deme compartment caused, for example due to treatment or death (Equation 1b the SIR model);

-migration: this corresponds to the displacement of a deme individual from its deme compartment to another deme compartment (e.g. from exposed to infectious following the reaction E → I in a SEIR model);

-sampling: this corresponds to the sampling of deme individuals (I in the SIR model) and leads to the end of the infectious period (e.g. through quarantine or change of behaviour).

Note that not all deme individuals can be sampled. For example, in an SEIR model where the E and I individuals are demes, one can consider that only infectious I individuals can be sampled. In other models, the sampled individuals can be only hospitalised individuals, or those in a chronic phase of a disease.

TiPS uses a coalescent approach [START_REF] Kingman | The coalescent[END_REF] to simulate phylogenetic trees based on trajectories, which correspond to a list of dated events (or 'reactions') and sampling dates (e.g. based on observed data). Each node in the simulated phylogenetic tree is associated with a state, i.e. a compartment name, and its height, i.e. its distance to the root.

TiPS can simulate the phylogeny of the entire trajectory or that of a sampled phylogeny if sampling dates are provided. These dates can either be provided by the user as a vector (bottom left panel of Fig. 1) or generated at random during trajectory simulation by adding a sampling reaction in the model. Since parameter values can vary over time, TiPS can reproduce temporal variations in sampling intensity. If sampled individuals can belong to more than one deme compartment, the user can choose between defining the state associated with each sampling event or defining a proportion of sampling events associated with each state. In the last case, TiPS randomly associates a state with each sampling date.

After this preprocessing, the sampling dates are organised as a named list containing a vector of decimal dates (with a column named 'Date') and a vector containing the reactions indicating the state of individuals to sample (a column named 'Reaction'). TiPS then incorporates these pieces of information into the recorded trajectory (which also contains dates and reactions) in chronological order.

The simulation of a phylogeny (sampled or not) starts from the most recent sampling date (or the most recent death event) and progresses through the simulated trajectory backward-in-time.

Each of the four types of reactions previously mentioned can result in a modification in the simulated tree. A forward-in-time transmission event can be represented by a coalescence between two lineages under a backwards-in-time process (or a branching under a forward-in-time process). A sampling event interrupting the transmission of the pathogen is represented as an external node (or leaf) in the phylogeny, and a death event, if observed, is also represented as a leaf. A migration event of an individual from one deme compartment to another will change and update the state of its corresponding lineage.

In some cases, especially with the tau-leap method, more than one event may occur on the same date (e.g. three new transmission events). To determine the number of events that lead to a change in the phylogeny (e.g. a coalescence), we draw a number from a hypergeometric distribution, which is appropriate since it describes the number of successful events (k) when drawing n times (without replacement) from a total population of size N that contains K elements corresponding to a 'success'. For example, when a transmission event in the trajectory is encountered in the phylogeny simulation, the algorithm draws, using the hypergeometric distribution, the number of sampled child lineages and the number of sampled parent lineages to coalesce into one sampled parent lineage. If, at time t, the number of sampled individuals (i.e., sampled nodes) is smaller than the number of simulated transmission events (n) in the trajectory, then the number of coalescence events between two sampled lineages is smaller than n. Therefore, all transmission events at time t may not lead to an observed coalescence in the tree. Similarly, when a migration event is encountered, the algorithm draws the number of sampled lineages associated with a change in state. Fig. 2 illustrates how the tree is simulated from the trajectory and further details can be found in the Appendix.

The output simulated phylogeny is an R object of class phylo as defined in the ape package [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF]. The simulated phylogenetic tree can be written in an output file in Newick (by default) or Nexus format if asked and specified by the user.

An illustration of how to simulate a phylogeny can be found in the 'Simulating phylogenies' box in Fig. 1 and in the Appendix for a logistic growth ecological model.

Ecological population dynamics models

The use of TiPS to simulate phylogenetic trees can also be applied to ecological population dynamics models. As in epidemiological models, there are 4 types of reactions: generation of new demes, migrations, removals, and samplings.

For example, TiPS can simulate the demographic history and the underlying phylogeny of two populations (N 1 and N 2 ) belonging to the same species living in different patches under a logistic growth assumption. In such a case, both compartments are deme compartments. The generation of a new deme individual (or 'birth') in a compartment depends on the growth rate and on the population size, and can be represented as coalescence between two lineages in a backward-in-time process. Migration corresponds to the actual migration of a deme individual from one patch to the other and leads to an update of the state of the corresponding lineage in the tree. A removal event corresponds to the death and, if observed, is represented as a leaf. Finally, sampling events correspond to an observation process that interrupts the biological process (e.g. birth) and are also represented as leaves in the tree. These events can correspond to the sterilisation or capture of the sampled individuals.

N1, N2 - ∅ r1N1(1-N1/K1) -----------→ N1 N1 → N1 + N1 new deme generation
Two population from the same species

dN1 dt = r1N1(1 -N1 K1 ) + µ(N2 -N1) ∅ r2N2(1-N2/K2) -----------→ N2 N2 → N2 + N2 new deme generation
that live in different patches (1 and 2)

dN2 dt = r2N2(1 -N2 K2 ) + µ(N1 -N2) N1 µN1 ---→ N2 N1 → N2 migration that are linked by migration events. N2 µN2 ---→ N1 N2 → N1 migration N1 dN1 ---→ ∅ N1 → ∅ removal N2 dN2 ---→ ∅ N2 → ∅ removal SEIR Epidemiological model: E, I S, R S βSI --→ E I → I + E new deme generation S: susceptible dS dt = -βSI E σE --→ I E → I migration E: exposed dE dt = βSI -σE I γI -→ R I → ∅ removal I: infectious dI dt = σE -γI R: removed dR dt = γI
Table 1 illustrates demes and reaction types for specific ecological and epidemiological models.

Benchmarking

To evaluate TiPS's performances, we designed a benchmarking analysis on both modules of the software package, i.e. the trajectory and the phylogeny simulators, comparing with existing R packages (adaptivetau and phydynR). Table S1 summarises the main features of the approaches used.

We first evaluated the computational speed and the accuracy of the trajectories (i.e. populating dynamics). For five initial population sizes and three different R packages, we simulated 10, 000 trajectories of the epidemiological Susceptible-Infected-Recovered (SIR) model and measured the execution time.

We then compared the time to simulate phylogenies under an SIR model. We varied the sampling proportions to obtain target phylogenies of various sizes (10, 100, 500, 1000, and 1500 leaves). We simulated 1,000 phylogenies under each sampling scheme using phydynR and our package.

Furthermore, assuming a more detailed epidemiological model with two host types, as described in [START_REF] Danesh | Quantifying transmission dynamics of acute hepatitis C virus infections in a heterogeneous population using sequence data[END_REF], we simulated an epidemiological trajectory using the tau-leap algorithm and a complete phylogeny, such that each end of infection event corresponded to a leaf in the tree. The simulation generated a full transmission chain corresponding to a phylogeny with 154, 507 leaves. From this complete phylogeny, we generated 10 subtrees by randomly sampling and keeping 1, 000 leaves for each subtree. We then used TiPS and phydynR to simulate 1, 000 phylogenies with each package under the same epidemiological model with the same parameter values. Note that the 1, 000 dates of each target subtree were imposed when simulating these phylogenies using a backwards-in-time approach. To compare the 2, 000 simulated phylogenies with the target one, we computed 60 summary statistics for each of them using the methods described in [START_REF] Saulnier | Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study[END_REF].

Results

In Figure 3(a), we show the median execution time for one trajectory simulation and the 50% interquartile envelope. TiPS and adaptivetau rely on Gillespie's Direct method (GDA), whereas phydynR uses the Euler-Maruyama integration (EMI). As expected, the population size, and hence the number of events per unit of time, increases the execution time for GDA-based packages but not for the EMI-based package. However, TiPS remains faster than the other two software packages for large populations (10 7 individuals). In Figure 3(b), we perform the same simulations using approximations of the Gillespie algorithm with fixed time steps. The computational speed of this GTA implemented in TiPS and adaptivetau is comparable and much faster than the GDA. Our new MSA algorithm outperforms existing methods and improves computational speed compared to our GTA, especially for small population sizes.

In Figure 3(c), we show the deterministic trajectory and, for each algorithm used, the mean simulated trajectory and its 90% confidence envelope. Among the packages studied, TiPS (in blue) is the one that yields the trajectories that are the closest to the deterministic prediction (in black). Note that there is a temporal shift for all the stochastic simulations, with a more rapid increase in population size compared to the deterministic model. This comes from the fact that stochasticity tends to favours trajectories that spread faster than average because they are less likely to go extinct. This known effect in population genetics models has also been described in epidemiology [START_REF] Hartfield | Epidemiological feedbacks affect evolutionary emergence of pathogens[END_REF].

We then analyse the median execution time to simulate phylogenies for each tool and each sampling In b, the cross shows the target phylogeny. The three colors represent three analyses using different target sub-trees. The phylogenies summary statistics used for the PCA are the same as in [START_REF] Saulnier | Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study[END_REF].

scheme (Figure 4(a)). TiPS's median simulation time is several orders of magnitude faster than that of the other method, with a more pronounced advantage for large phylogenies.

Regarding the accuracy of the simulations, we see in Figure 4(b) that the cross, which indicates the projection of the summary statistics from the target phylogeny, is contained in the envelope containing 90% of the phylogenies simulated using TiPS but not of those simulated using phydynR, and that for three different target phylogenies used represented by distinct colours. We observed the same results for seven other target phylogenies used. This means that the target phylogeny cannot be distinguished from phylogenies simulated using the same model with the same parameters using our package. The discrepancy observed for phydynR could originate from the trajectory, which strongly differs from the deterministic prediction (Figure 3c). Furthermore, TiPS is favoured in this analysis because the same algorithm (Gillespie's Tau-Leap) was used to simulate the target phylogeny and the 1, 000 sub-trees. Supplementary Figure S7 shows the distributions of summary statistics computed from the phylogenies simulated using TiPS and phydynR for each analysis using a different target tree.

Discussion

We developed a flexible R package to rapidly simulate trajectories and phylogenies from compartmental models. Its structure allows the user to include several sources of heterogeneity between individuals or between populations, e.g. different life stages or metapopulation structures. The simulation of the phylogeny relies on the trajectory and involves a coalescent approach.

Our benchmarking analyses show that TiPS is comparable to or outperforms existing R packages in terms of speed when generating numerous trajectories or phylogenies. The same is true for the accuracy of the simulation outputs.

The first asset of this software package is its flexibility since it can readily be used for any compartmental model. Another asset is its computation speed as it can simulate trajectories in a matter of milliseconds on a regular desktop computer. These properties have already been used for infection phylodynamics studies involving Approximate Bayesian Computing (ABC) [START_REF] Danesh | Quantifying transmission dynamics of acute hepatitis C virus infections in a heterogeneous population using sequence data[END_REF] or to illustrate the effect of superspreading events [START_REF] Alizon | Superspreading genomes[END_REF].

In the context of infectious diseases, when simulating phylogenetic trees with TiPS, we assume that the time of a transmission event is the same as the time of coalescence in the tree. However, this is a simplification since the time of coalescence in the phylogenetic tree should take place before transmission, during the infection of the 'donnor' host [START_REF] Ypma | Relating Phylogenetic Trees to Transmission Trees of Infectious Disease Outbreaks[END_REF]. A variety of withinhost evolutionary processes may further weaken this assumption. This limitation is common to virus phylodynamics studies and an active line of research [START_REF] Volz | Phylodynamic Inference across Epidemic Scales[END_REF].

Beyond epidemiology, this software package can be used more broadly to simulate population dynamics and associated genealogies. Some future extensions of TiPS will consist in introducing non-Markovian dynamics, simulating multifurcating phylogenies, and implementing other optimised stochastic simulation algorithms [START_REF] Cao | The slow-scale stochastic simulation algorithm[END_REF][START_REF] Cao | Adaptive explicit-implicit tau-leaping method with automatic tau selection[END_REF].

reaction can also be written as Z + Xi -→ Xj + Xi if it includes a non-deme individual. In this section, we will refer to a donor lineage as the parent lineage and a recipient lineage as the child lineage. A birth reaction leads to one of two types of modification in the tree: a coalescence or an 'invisible' coalescence. A coalescence corresponds to the coalescence of two sampled lineages from two individuals, the recipient (or child) and the donor (or parent), into one sampled individual lineage (the donor). An 'invisible' coalescence is defined as the coalescence of the sampled lineage of the recipient (here X j ) and an unsampled lineage of the donor (X j ) into one unsampled lineage (the donor, X i ).

First, we need to determine the number of recipient lineages nrecipient. Again, assuming a hypergeometric distribution, we compute the number of recipient lineages nrecipient if we sample nrep times without replacement in a sample of size |Xj| containing |X j | individuals (see equation S2a). If there are no recipient lineages (nrecipient = 0), there is no coalescence or invisible coalescence and, therefore, no change in the tree. Otherwise, we need to determine the number of donor lineages n donor . Since the sampling is performed without replacement, the total sample size is updated to |Xj| -nrep and the number of daughter lineages Y represented by nodes that are available for the coalescence becomes |X j | -nrecipient. Thus, we compute the number of donor lineages n donor if we sample nrep times without replacement in a sample of size |Xj| -nrep containing |X j | -nrecipient individuals with i = j S2b, or in a sample of size |Xj| containing |X j | individuals with i = j S2c. Mathematically, we can write:

nrecipient ∼ HyperGeom(nrep, |X j |, |Xj|) (S2a) n donor ∼ HyperGeom(nrep, |X i | -nrecipient, |Xi| -nrep) if i = j (S2b) n donor ∼ HyperGeom(nrep, |X i |, |Xi|) if i = j (S2c)
Knowing the number of recipient lineages, we need to determine the number of coalescences (nC ), i.e. the number of lineages among the n donor sampled donor lineages that will coalesce with the recipient ones. The number of visible coalescences is drawn from a hypergeometric law where we sample nrecipient times in a sample of total size nrep containing n donor sampled donor lineages S3a. The number of invisible coalescences, nIC , is the number of remaining recipient lineages that have not coalesced with sampled donor lineages S3b.

nC ∼ HyperGeom(nrecipient, n donor , nrep) (S3a) nIC = nrecipient -nC (S3b)
Upon a coalescence event, |X j | is decreased by one, a node is randomly picked in X i , its height is updated to t, and it is linked to another node that is removed from X j . For an invisible coalescence event, |X j | is decreased by one and both |X i | and |X i | are increased by one. In this case, a new node from X is created with height t and linked to a node randomly picked in X j . In both cases, |Xi| is decreased by one (which is already known from the trajectory).

Removal event. By default, a removal reaction does not require any tree modification. However, if we want to simulate a full tree instead of a sampled one, the sampling events will correspond to the death events, which will therefore lead to the addition of a node. In this case, the number of sampled removal events is nSD = nrep.

Migration event.

Only migration events involving deme individuals can lead to a modification in the tree.

These migration reactions can be written as Xi -→ Xj, with j = i. We assume that the number of migrations that lead to a tree modification is given by the following hypergeometric distribution:

nM ∼ HyperGeom(nrep, |X j |, |Xj|) (S4a)
A migration increases |X i | and decrease |X j | by one. A new node is created in X i with height t and linked to a node randomly picked in X j , which is then removed from X j . Furthermore, Xi is incremented by one and Xj is decremented by one.

S2 Application to an epidemiological SI a I c R model

We illustrate the functioning of TiPS using an SIaIcR epidemiological compartmental model, where individuals can be susceptible (with density S), infected in acute phase (Ia), infected in chronic phase (Ic), and removed (R). The corresponding ODE system is

dS(t) dt = -β S(t) Ia(t) -β S(t) Ic(t) (S5a) dIa(t) dt = β S(t) Ia(t) + β S(t) Ic(t) -σ Ia(t) (S5b) dIc(t) dt = σ Ia(t) -γ Ic(t) -α Ic(t) (S5c) dR(t) dt = γ Ic(t) (S5d)
where β is the infectious contact rate, γ the recovery rate, α the virulence, and 1/σ the expected duration of the acute phase. The model can be described as an individual-based model using a system of reactions:

S βSIa ---→ Ia (S6a) S βSIc ---→ Ia (S6b) Ia σIa --→ Ic (S6c) Ic γIc --→ R (S6d) Ic αIc --→ ∅ (S6e)
where the rate of occurrence of each reaction is indicated above the transition arrow. Reactions S6a and S6b are transmission (i.e. new deme individual creation) reactions, S6c and S6d are migration reactions, and finally S6e is a removal reaction event.

TiPS will first build and generate a function to simulate trajectories using this system of reactions. Population dynamics are simulated by providing parameter values using one of the algorithm described in the main text. In the following, we focus on the model captured by system of equations S6.

In this approach, under this epidemiological model, we trace back the epidemiological history of the sampled virus. Hence, the deme compartments, i.e. the ones sampled, where the virus is present and then contributing to the phylogeny, are Ia and Ic.

Once the trajectory is simulated, to simulate a sampled phylogeny, TiPS requires the sampling dates and the proportion of the sampling dates to be associated with each type of deme. Let us assume, for example, that 15% of the sampling dates are associated with the Ia deme and 85% with the Ic deme compartment. TiPS randomly assigns each sampling date to a deme compartment based on these ratios and adds the dates to the list of events in the simulated trajectory (see Supplementary Figure S5). The R code to build the simulator, simulate a trajectory and a phylogeny are shown in Supplementary Figure S1.

Deme compartments (Ia and Ic) can be composed of individuals represented by nodes in the simulated tree (belonging to the sub-compartments I a and I c ). These individuals represented by nodes can also be still unsampled (belonging to I a and I c ).

TiPS starts the simulation of the phylogeny from the most recent sampling event and follows the trajectory backwards in time.

If the backward step in the trajectory leads to one or multiple migration events, i.e., in this model, from the acute phase compartment (Ia) to the chronic infectious phase compartment (Ic) as in reaction S6c, the number of tree modifications is given by the following hypergeometric distribution:

nM ∼ HyperGeom(nrep, |I c |, |Ic|) (S7)
An illustration of the tree update is shown in Supplementary Figure S2.

If the backward step in the trajectory leads to a transmission reaction, two tree modifications are possible: a coalescence and an invisible coalescence. Here, we consider the donor as the host transmitting the pathogen, and the recipient the one that has been infected. In this model, there are two different demes and two different transmission reactions (see reactions S6a and S6b). Given the transmission reaction S6a, the number of coalescences nC and invisible coalescences nIC leading to tree modifications are governed by the following

I a →I c Figure S2: event in the SI a I c R model.
The evolution of a pathogen lineage in an individual during the acute phase of its infection is represented in red. The individual then enters the chronic phase of his infection and the pathogen lineage continues to evolve (represented in blue branches). The solid branch represents the evolution of the sampled pathogen lineage and the dashed branch represents the evolution of an unsampled pathogen lineage. The unsampled lineage is eventually removed and not represented in the final simulated phylogeny.

n RS(Ia) ∼ HyperGeom n S(Ia) rep , |I a |, |Ia| -(|Ia| -|I a |) (S10a) n S(Ia) = n S(Ia) rep -n RS(Ia) (S10b) n RS(Ic) ∼ HyperGeom n S(Ic) rep , |I c |, |Ic| -(|Ic| -|I c |) (S10c) n S(Ic) = n S(Ic) rep -n RS(Ic) (

S4.2 Benchmarking: phylogenies

To evaluate the accuracy of phylogeny simulations, we generated 10 different sub-trees of 1, 000 leaves, using detailed epidemiological model with two host types, as described in [START_REF] Danesh | Quantifying transmission dynamics of acute hepatitis C virus infections in a heterogeneous population using sequence data[END_REF]. We then used TiPS and phydynR to simulate 1, 000 phylogenies with each package under the same epidemiological model with the same parameter values, imposing the 1, 000 dates of each target sub-tree under a backwards-in-time approach.

To compare the simulated phylogenies with the target one, we computed summary statistics for each of them using the methods described in [START_REF] Saulnier | Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study[END_REF]. Supplementary Figure S7 shows the distributions of summary statistics computed from the phylogenies simulated using TiPS (in red) and using phydynR (in orange) for each analysis using a different target tree (in black).

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q maxHeight meanA meanE meanI_1 medianA medianE minHeight Target 10 TiPS 10 phydynR 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q meanI_2 meanIE_1 meanIE_2 medianI_1medianI_2medianIE_1 tpsT_1 Target 10 TiPS 10 phydynR 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q meanI_3 meanIE_3 medianI_3medianIE_2 medianIE_3 tpsT_2 Target 10 TiPS 10 phydynR 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q dRate gdRatio gRate maxLin timeMaxLin LTT Target 10 TiPS 10 phydynR 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q colless ILnodes maxDiffW to the summary statistics computed from the phylogenies simulated using TiPS, orange distributions correspond to the summary statistics computed from the phylogenies simulated using phydynR, for each target tree analysis. Black dots represent the values of summary statistics computed from each target tree.

Figure 1 :

 1 Figure 1: General structure of the TiPS pipeline. The equations and outputs correspond to the SIR epidemiological model(Keeling and Rohani, 2008). The functions of the R package are in blue. The simulator of trajectories, which is built as a function, is in orange. The variable traj, in red, is the output trajectory of class simutraj. It can be plotted using our plot method. TiPS used the simulated trajectory and 100 sampling dates that we generated (variable sampling_dates in green) to simulate a sampled phylogeny. The output simulated phylogeny is a phylo class object from the ape R package[START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF], which can be used for plotting.

Figure 2 :

 2 Figure2: Tree simulation for a SIR model. The trajectory shows the series of epidemiological events and sampling events. The phylogeny represents the epidemiological history of the individuals carrying sampled viruses in solid lines. The rest of the history is shown in dashed lines. In this representation, at each transmission event (branching), the donor's pathogen lineage is deviated to the right side and the recipient's (i.e. newly infected individual) pathogen lineage to the left side.

Figure 3 :

 3 Figure 3: Benchmark analyses of trajectories simulations. a) Median computation speed and 50% confidence intervals (CI) for one simulation using GDA or EMI methods, b) median computation speed and 50% CI for one simulation using approximating GDA methods, and c) mean resulting trajectories for prevalence time series and the 90% CI. GDA stands for Gillespie's Direct Method and EMI for Euler-Maruyama Integration. 10, 000 trajectories simulations were performed under an epidemiological SIR model, for five initial population sizes N varying from 10 3 to 10 7 and with parameter values R 0 = 2, γ = 1/3 and β = R0γ N .

Figure 4 :

 4 Figure 4: Benchmarking analysis of phylogenies simulations. a) Median computation speed for simulating one phylogeny of a given size and b) First two axes of a Principal Component Analysis (PCA) of phylogenies summary statistics. In a, the shaded area shows the 95% confidence intervals.In b, the cross shows the target phylogeny. The three colors represent three analyses using different target sub-trees. The phylogenies summary statistics used for the PCA are the same as in[START_REF] Saulnier | Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study[END_REF].

  S+I a →I a +I a (b) Births without tree modifications S+I c →I a +I c (d) Births without tree modifications S+I a →I a +I a (a) Births inducing tree modifications S+I c →I a +I c (c) Births inducing tree modifications

Figure S3 :

 S3 Figure S3: Transmission events in the SI a I c R model. Each transmission event from a living infectious individual to a susceptible individual can be represented by a branching. We show the donnor lineage on the right side of the branching and the recipient pathogen lineage (i.e. the newlyinfected individual) on the left side. Dots correspond to nodes in the resulting phylogeny. Color and branch line codes are identical to Figure S2.

Figure S6 :

 S6 Figure S6: Simulating a trajectory and a phylogeny using TiPS given a logistic growth model. The equations and outputs correspond to the logistic growth model. The functions of the R package are in blue. The simulator of trajectories built as a function is in orange. The variable traj in red is the output trajectory. The phylogeny is plotted using a function from the ape R package. Paradis and Schliep (2019).

Figure S7 :

 S7 Figure S7: Distributions of summary statistics. The summary statistics are grouped into three families : branch lengths (denoted BL), tree topology (TOPO) and Lineage-Through-Time plot (LTT). The dots represent the median and the horizontal lines represent the 95% HPD. Red distributions correspond

Table 1 :

 1 Description of demes and reactions for an ecological and an epidemiological model. 'ODEs' stands for ordinary differential equations.

	Model Description	ODEs	Deme	Non-deme	Individual-based	Tree deme individual-based	Tree reaction
				compartments compartments	reactions	reactions	type
	N1N2	Logistic growth in two patches model :

  Samplings in the SI a I c R model. Colors, branches, and dots code are identical to FigureS3. Each sampling event leads to the addition of a node. Re-sampling events (left side of the figure) occur when the pathogen lineage has already been sampled but the individual currently carrying it has never been sampled (individuals can only be sampled once but can transmit after sampling). Otherwise, we have a classical sampling event (right side of the figure), i.e. the pathogen has not been sampled yet. S+I a →I a +I a S+I a →I a +I a S+I a →I a +I a I a →I c S+I c →I a +I c I c S+I c →I a +I c I a →I c I c → S+I a →I a +I a I c I a →I c I c S+I a →I a +I a I a →I c I a

	S+I a →I a +I a S+I c →I a +I c trajectory Figure S4: Epidemiological simulated events sampling events →I a	I a I c
		S10d)

Table S1 :

 S1 Main algorithms of the R packages compared. GDA stands for Gillespie's Direct Algorithm, GTA for Gillespie's Tau-Leap Algorithm, and MSA for Mixed simulation algorithm.

	Package	Methods	Features
	TiPS	GDA (exact)
		GTA (approximate) Requires a fixed time-step τ
		MSA (mixed)	Requires a fixed time-step τ
	adaptivetau GDA (exact)
		GTA (approximate) Requires a fixed time-step τ
	phydynR	Euler-Maruyama	Requires a time-step τ

Supplementary Information S1 Phylogenies simulation algorithm

TiPS simulates rooted and binary phylogenies, which means that every internal node has exactly two daughter nodes. The root of the tree represents the ancestor of all lineages, and the tips of the branches (or leaves) represent the most recent descendants of that ancestor. A coalescence between two lineages using a backwardin-time process is equivalent to a branching event in the phylogeny in the forward-in-time process. Under the hypothesis of neutral evolution that phylodynamics relies on, a branching event can represent, for example, a viral transmission event, and a leaf can represent the end of viral infection. Note that we will refer to the height of a node as its distance to the root.

In this section, we distinguish two types of compartments: the deme compartments, where individuals contribute directly or indirectly to the phylogeny, and the non-deme compartments that cannot be placed in the genealogical process. For example, in a SIR epidemiological model, the deme and non-deme compartments would be respectively the I compartment and the S and R compartments. Each deme compartment is denoted Xi, with i ranging from 1 to the number of deme compartments in the model. The sampled individuals belong to the sub-compartment X i (X i ⊂ Xi) and are all associated with a leaf in the tree. We also introduce X i , the sub-compartment of X i (X i ⊂ X i ) corresponding to individuals in X i that have not yet been but may be sampled in a backwards-in-time process. The discrete size of compartments Xi, X i and X i at time t are denoted as |Xi|, |X i | and |X i |, respectively. A non-deme compartment is denoted Z and its discrete size |Z|.

The tree simulation starts from the last (i.e. most recent) sampling date and progresses through the simulated trajectory backward-in-time. The number of events that lead to a change in the phylogeny is drawn from a hypergeometric distribution. The hypergeometric distribution is appropriate as it describes the number of events k from a sample n drawn from a total population N without replacement. Each of the four types of reactions (sampling, new deme generation (or birth), removal, and migration) can result in a modification in the simulated tree : a new external node (or leaf) or the coalescence of two lineages.

Sampling event.

We define a sampling event as an event that interrupts the biological process of an individual in the population of interest and a re-sampling event as an observation event. In an ecological context, the sterilisation of an animal would be a sampling event whereas marking the animal would be a re-sampling event. In an epidemiological context, that is, when tracing back the history of pathogen spread, a sampling event usually corresponds to a host individual being detected. Usually, this is assumed to coincide with the end of the infectious period, because the individual will be isolated or use adequate protection to prevent further transmission. However, if this assumption does not apply, it is possible to treat the event as a re-sampling event, such that the host individual will continue to transmit the pathogen to other individuals.

When we know the sampling dates, we know the number of sampling events occurring at time t (n S rep ) and, therefore, the number of tips to create with height t. However, since we allow sampled lineages to continue to have an offspring, we need to determine if the nodes we sampled at time t are associated with re-sampling events or not (i.e. if they should be linked to a node that has already been sampled after time t in our coalescent approach). The number of re-sampling (nRS) and sampling events (nS) at time t is governed by the following relationships:

where n S rep is the total number of phylogeny nodes generated at t. Using the hypergeometric distribution, we compute the number of re-sampling events nres if we sample n S rep times without replacement in a sample of size |Xi| -

The number of 'classical' samplings, nS, is the difference between nRS and the number of tips to be generated (n S rep ). In the case of a re-sampling, we randomly pick a node from X i (the pool of individuals who have not yet been sampled), update its height to t, and link it to a node from X i (the pool that has been sampled). Each re-sampling event decreases |X i | by one. In the case of a 'classical' sampling event, a new node with height t is created in X i and |X i | increases by one.

Birth event.

A birth reaction corresponds to the generation of a new individual of the deme and can be written as follows: Xi -→ Xj + Xi where i and j range from 1 to the number of deme compartments. The birth

Given birth reaction S6b, the number of coalescences and invisible coalescences leading to tree modifications are governed by the following relationships:

Supplementary Figure S3 illustrates possible tree updates.

Tree

S3 Application to a logistic growth model

In this section we use the African Savannah elephant (loxodonta africana) population in the Kruger National Park as an example that has demonstrated an important increase in the last two decades. A study has shown that a high number of the elephants can damage the Park's ecosystem [START_REF] Cumming | Elephants, woodlands and biodiversity in southern Africa[END_REF]. To prevent that, until the early 90's, a solution was the culling of elephants. Due to ethical issues, another solution has been proposed where female elephants would be on contraceptives [START_REF] Whyte | Managing the elephants of Kruger National Park[END_REF].

Here we use a logistic growth model to study the population dynamics. The ODE system is:

where N is the elephant population density, K is the carrying capacity of the environment, and r is the intrinsic growth rate of the population where r = bd with b the birth rate and d the death rate. The model can be described as an individual-based model using a system of reactions:

where the rate of occurrence of each reaction is indicated above the transition arrow. Reaction S12a is a birth reaction and S12b a removal reaction. TiPS allows to simulate a phylogeny without sampling dates, where events interrupting the biological process, such as removal events or here sterilisation events, simulated in the trajectory are represented as leaves in the phylogeny. To illustrate this module of the tool, we add a sterilisation rate to the model to simulate the events in the trajectory. The system of reactions becomes:

where reaction S13c is the sterilisation reaction with a rate of its occurrence indicated above the transition error.

Using this system of reactions, TiPS will first build and generate a function to simulate trajectories. Population dynamics are simulated by providing parameter values using one of the algorithm described in the main text. The R code to build the simulator, simulate a trajectory and a phylogeny are shown in Supplementary Figure S6.

We assume that a sampling event is an event that interrupts a biological process. In this example, where no sampling dates are required, the removal events and the elephant sterilisation events will be considered as the sampling events and will be represented as leaves in the simulated phylogeny.

We introduce here the compartments N , N and N where N ⊆ N ⊆ N . All the elephants are in compartment N , the sampled elephants are sub-compartment N and the elephants that in N that have not yet been but may be sampled are in sub-compartment N .

TiPS starts the simulation of the phylogeny from the most recent sampling event follows the trajectory backwards-in-time.

When the backward step in the trajectory at time t leads to n death or sterilisation events, n new nodes are created with height t in N and |N | increases by one. Note that, each node is labelled with the corresponding reaction, so we can distinguish and visualise them when plotting the phylogeny.

If the backward step in the trajectory leads to a birth reaction, two tree modifications are possible: a coalescence and an invisible coalescence. A coalescence in the phylogeny corresponds to the coalescence of two sampled lineages each representing an elephant, into one sampled individual (the donor). An 'invisible' coalescence is the coalescence of the sampled lineage of the recipient individual (N ) and an unsampled lineage of the donor individual (N ) into one unsampled lineage (the donor, N ). Given the birth reaction S13a, the number of coalescences nC and invisible coalescences nIC leading to tree modifications are governed by the following relationships:

where nrep is the number of birth events in the trajectory at time t of the trajectory, nrecipient is the number of recipient lineages and n donor the number of donor lineages. Upon a coalescence, |N | is decreased by one, a node is randomly picked from N with height t and is linked to another node that is removed from N . For an invisible coalescence, |N | is decreased by one both and N and N are increased by one. A new node from N is created with height t ans is linked to a node randomly picked in N . In both cases, |N | is decreased by one as recorded already in the simulated trajectory.

S4 Benchmarking

S4.1 Benchmarking: trajectories

To evaluate our simulator we performed a benchmarking analysis on both modules of the tool, i.e. the trajectory and the phylogeny simulators, using two existing R packages. adpativetau performs simulations of trajectories of continuous-time Markov processes by using Gillespie's stochastic simulation algorithms. phydynR is a package for phylodynamic inference using population genetic models and performs simulations of trajectories as well using the Euler-Maruyama integration method, and provides methods for simulating trees conditional on a demographic process. The different algorithms of each package are presented in Table S1.


Abdelmounaim Bouroudi
email: abdelmounaim.bouroudi@nokia.com

Abdelkader Outtagarts
email: abdelkader.outtagarts@nokia-bell-labs.com

Yassine Hadjadj-Aoul
email: yassine.hadjadj-aoul@irisa.fr

Robust Deep Reinforcement Learning Algorithm for VNF-FG Embedding

Keywords: Virtual Network Embedding, Deep Reinforcement Learning, Algorithm Selection, Stationary bandit, B5G/6G

Network slicing, also known as the virtual network embedding (VNE) problem, is an NP-hard optimization problem. Compared to traditional approaches, the methods relying on deep reinforcement learning yield better performance without exhibiting issues such as stacking at local minima and/or solutions' space exploration limits. These algorithms present, however, different performances according to the employed approach, and the problem to be treated, resulting in robustness problems. To overcome these limits, we propose the adoption of the best algorithm, from a selection of learning strategies, in terms of reward and sample efficiency at each time step. The proposed strategy acts as a meta-algorithm that brings more robustness to the network by dynamically selecting the best solution for a specific scenario. Our solution proved its efficiency and managed to dynamically select the best algorithm in terms of the best acceptance ratio of the deployed services and outperform all the stand-alone algorithms.

I. INTRODUCTION

The sixth generation of mobile networks is expected to provide a multi-purpose network with massive machine-tomachine communications and a multitude of services with diverse requirements, thereby requiring a new architecture to better meet such needs [START_REF] Dogra | A survey on beyond 5g network with the advent of 6g: Architecture and emerging technologies[END_REF]. To meet such requirements, network providers are adopting new technologies including network function virtualization (NFV) and software-defined network (SDN). These technologies allow the virtualization of the physical infrastructure and the dynamic creation and sharing of network instances and services which enable the concept of network slicing.

The coexistence of different network functions on the same substrate network has led to the NFV placement problem (also known as Virtual Network Embedding "VNE"), in which we look to efficiently allocate resources for network functions an links connecting them (forming a Forwarding Graph). The VNE is an NP-hard problem with several existing solutions [START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF]. Mathematical optimization is only able to deal with very small network instances [START_REF] Allybokus | Virtual function placement for service chaining with partial orders and antiaffinity rules[END_REF]. Heuristics converge extremely fast, but reach typically only local minima [START_REF] Khebbache | Scalable and cost-efficient algorithms for vnf chaining and placement problem[END_REF]. Meta-heuristics, which are only applicable to medium-sized networks, suffer from a lack of knowledge accumulation [START_REF]A multi-objective non-dominated sorting genetic algorithm for vnf chains placement[END_REF].

Given these limitations, machine learning approaches have been explored, and more particularly Deep Reinforcement Learning (DRL), which are more appropriate to deal with the unpredictable behavior of the virtualized services. The authors, in [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF], [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF], proposed a new DRL strategy based on the DDPG algorithm [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF], with improvements in the exploration capacity combined with a heuristic in order to guarantee in the worst case the performance of the latter. The proposed solution provides significantly better results compared to a vanilla DDPG. However, the main drawback of this solution is their lack of robustness, which is a critical aspect to consider in next generation cellular networks [START_REF] Li | Trustworthy deep learning in 6g-enabled mass autonomy: From concept to quality-oftrust key performance indicators[END_REF].

Since existing DRL-based algorithms are case-specific solutions, their performances can be affected by considering new unseen scenarios. This creates a risk to the network making such standalone solutions unsafe. Moreover, these agents present different learning potentials, where some algorithms allow quick learning but prematurely converge to a local optimum, while others learn slowly in the beginning but yield better performance at the end.

To solve these problems, we propose a new robust resource allocation algorithm based on the Epochal Stochastic Bandit Algorithm Selection (ESBAS) framework [START_REF] Laroche | Reinforcement learning algorithm selection[END_REF], that we adapt to the context of the VNE problem. ESBAS is an off-line Algorithm Selection (AS) framework for off-policy RL algorithms. Its principle is to select among a set of RL algorithms the one maximizing the expected return. The algorithm splits the execution time to epochs of exponential size, in which the set of algorithms is not updated, and the selection is done using a rebooted stochastic bandit with accumulation of knowledge using a buffer to store the executed episodes. This framework allows both a competitive execution between different RL agents by selecting the most promising agent at each iteration, and a cooperative exploration of the action space through the transition sharing mechanism.

The rest of this paper is organized as follows. We start by a background section, where we present the VNE problem and formulate it as a Markov Decision Process (MDP). Then, we present an offline AS algorithm based on ESBAS. In the fifth section, the simulation results of offline selection are presented alongside with a discussion and a comparison with already existing standalone algorithms. Finally, we provide a conclusion in section V.

II. BACKGOUND

In this section, we introduce the modeling of the VNE problem and its formulation using an MDP, which is common to all the considered learning techniques.

A. VNE Problem Formulation

In the VNE problem we aim to deploy a set of V VNF-Forwarding Graphs (VNF-FG) noted G = {G 1 , ..., G V } into a substrate network. On one hand, Each VNF-FG G ζ is described by a directed graph with a set of N Virtual Network Functions (VNFs) where each VNF n requests a set of k VNF type of resources represented by a vector h n = [h n ,0 , .., h n ,KVNF-1 ]. These VNFs are connected by a set of L virtual links and each link l requests k VL types of resources denoted by a vector f l = [l n ,0 , .., l n ,KVL-1 ]. On the other hand, the substrate network is also represented by a directed graph where each node n has a specific amount of each resource k denoted r n,k , and each link has also a specific amount of resource k denoted r l,k .

To deploy a VNF-FG, we need to deploy all its VNFs and its VLs. A VNF can be deployed on a specific node if this last has sufficient resources for each requested resource and Each VNF is deployed on only one substrate node. The same requirements are applied on the VLs with an aditional condition that to deploy a VLs, all VNFs it connects must be deployed.

We define an auxiliary binary variable g ζ indicating if the VNF-FG ζ is deployed. The objective is to maximize the number of deployed VNF-FG called the Acceptance Ratio (AR), which is defined by:

AR = ζ g ζ V (1)

B. Markov Decision Process formulation

Reinforcement Learning consist in learning through trial and errors using an agent that interact with an environment [START_REF] Kaelbling | Reinforcement learning: A survey[END_REF]. Each time steps t, the agent receives the state of the environment s t , and based on it, the agent chooses an action a t to execute on the environment following a policy (π : S → A). The environment return a feedback on the executed action known as reward r t with the new environment state s t+1 . This feedback is used to optimize the policy with the goal of maximizing the expected reward. This process can be modeled as an MDP. In this paper, we use the same MDP formulation presented in our previous paper [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF]: , w l,l t ) ∈ [0, 1] showing the favor of assigning VNF n to node n (resp. virtual link l to substrate link l ) at time-step t. As the action is sparse, we propose using a binary versions of the action vector to convert it to a feasible solution by applying the WF2A heuristic [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF].

• Reward: to assess the performance of the RL algorithm, the Acceptance Ratio defined in ( 1) is adopted as a reward measure.

III. A NEW FRAMEWORK FOR DRL-BASED ALGORITHMS' SELECTION

As DRL based algorithms have different performances depending on the instances of the problem, it is relevant to have a meta-algorithm solution allowing the association of each algorithm to an instance of the problem to obtain better performances. This solution is a selection algorithm framework.

Formally, the algorithms' selection is a sequential decision problem D with a set of instances I, a set of algorithms A = A 1 , A 2 , ..., A n to solve D and a metric m : A × I → R allowing the qualitative measurement of the performance of each algorithm for a problem instance [START_REF] John | Rice. the algorithm selection problem[END_REF]. The goal is, thus, to define a selector S that maps each instance i ∈ I of the problem to an algorithm S(i) ∈ A that allows the optimization of the overall performance according to the metric m.

In our case, D is the VNE problem, the algorithms are the different DRL agents, and the metric is the Acceptance Ratio returned by the agent. Constructing such a selection algorithm allows the combination of the strengths of multiple heterogeneous agents (algorithms). This is achieved through a competition and a cooperation mechanisms [START_REF] Zheng | Competitive and cooperative heterogeneous deep reinforcement learning[END_REF]. The competitive aspect appears in the execution between different RL agents by selecting the most promising agent at each iteration, while the cooperative aspect concerns the exploration of the action space through the transition sharing mechanism. In the following subsections, we first describe the unifying framework for RL agents. Then we introduce the proposed offline algorithm for the selection of the best agent.

A. Unifying view for Reinforcement Learning agents

To apply the Algorithm Selection framework for RL agents, we propose a unified view for the different agents independently of their learning mechanism. This unified view allows the Algorithm Selection solution to deal with these agents as black boxes and select the best one, relying only on their returns. For this purpose, we assume that we have a set of offpolicy RL algorithms called portfolio that learn from the same shared data-set (fairness assumption). This framework was first introduced in [START_REF] Laroche | Reinforcement learning algorithm selection[END_REF] for general purpose RL algorithms. In this paper, we adapt this framework to the VNE problem and we propose the workflow presented in Fig 1 . The process starts with a selection of an algorithm A i from the portfolio (1). In (2), the selected algorithm is used in the execution environment to generate a trajectory τ (3) and returns a cumulative reward µ( τ ) (4.2). On one hand, the generated trajectory (4.1) is stored in a shared memory (buffer) that will be used to update all the algorithms in the portfolio (5.2). On the other hand, the cumulative reward is used to update the Algorithm Selection policy corresponding to the selected algorithm (5.1). The final global goal is to optimize the expected cumulative reward defined as:

E S T τ =1 µ(ε τ ) (2)
Based on these assumptions, and using these formulations, we propose in the next section two meta-algorithms which allow the selection of DRL agents to solve the VNE problem in two manners, offline selection and online selection with attention.

B. Offline Algorithm Selection

Based on [START_REF] Laroche | Reinforcement learning algorithm selection[END_REF], we propose an adapted version of Epochal Stochastic Bandit Algorithm Selection (ESBAS) framework which is a meta-algorithm used to solve the off-policy RL Algorithm Selection problem by considering it as an Exploration/Exploitation problem. In their work, Authors relied on the Upper Confidence Bound (UCB1) algorithm [START_REF] Katehakis | The multi-armed bandit problem: decomposition and computation[END_REF] to assure a good trade-off between exploration and exploitation of the different agents' performances and exploration abilities. The main idea behind ESBAS is to block the update for all agents during exponential sized epochs in order to create a stationary system allowing the selection using UCB1.

In Alg. 1 we present the ESBAS algorithm used for VNE problem. In this algorithm, we have two main parts: The first part, from line 2 to 5, the algorithm starts by updating all agent's policies using the past recorded trajectories only if we have enough records in the reply buffer (threshold G). The second part starts from the 6th line in which we initialize the UCB1 parameters. From line 7 to 17 we execute UCB1 on exponential sized epochs. During the epoch, the meta-algorithm computes the optimistic guess using UCB1 formula in line 8 and generates a trajectory using the selected algorithm from line 10 to 14 (in VNE problem, a trajectory is a number T of placements). In each trajectory we calculate the accumulated reward by the executed agent, and we store the new record to the buffer δ. The last step, from line 15 to 17, allows the update of the different parameters of the UCB1 algorithm.

IV. SIMULATION RESULTS

To simulate the different algorithms, we had first to setup the DRL parameters as mentioned in [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF]. For the substrate

Algorithm 1: ESBAS for VNE problem Input: δ 0 , A, G, T 1 g ← 0 2 for β ← 0, ∞ do 3 for A k ∈ A do 4 if g ≥ G then 5 π k δ 2 β -1 : learn π k on trajectories of epoch 2 β 6 n ← 0, ∀A k ∈ A, n k ← 0, x k ← 0 7 for τ ← 2 β , 2 β+1 do 8 A kmax = argmax A k x k + ξ log(n) n k 9 µ ← 0 10 for t ∈ 0, T do 11 Execute a k t generated by A kmax → r t , s t+1 12 (s t , a k t , r t , s t+1 ) → δ 13 µ ← µ + r t 14 g ← g + 1 15 x kmax ← n kmax x kmax +µ n kmax 16 n kmax ← n kmax + 1 17 n ← n + 1
network, we use the BtEurope network topology [START_REF] Knight | The internet topology zoo[END_REF] with 24 nodes and 37 full-duplex links. For the VNF-FG request we use 10 connected graphs generated using the Erdos-Renyi model [START_REF] Erdős | On the evolution of random graphs[END_REF]. In table I, we present the amount of resources available on substrate network, as well as the resource values requested randomly by the VNF-FGs. The portfolio consists of three off-policy DRL agents, we used three Off-policy DRL agents already used for VNE problem. Each agents have Actor and Critic Networks. Each network is a Dense layer with 300 hidden units each and a ReLU activation function. For the output we use a Sigmoid activation function for the Actor Network and a linear activation function for Critic Networks.The description of the different agents is as follow:

• W2FA DDPG agent: which is the ordinary DDPG algorithm the W2FA heuristic [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF]. • EEDDPG agent: which is the implement the EEDDPG in [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF] without the evolutionary operation. To test the performances of our solution, we used the OMNeT++ simulation tool to deploy the topology based on the same platform used by authors in [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF]. ESBAS algorithm is run for 10 exponential-sized epochs with β = [[0, 9]]. Fig. 2 shows the Acceptance Ratio of the ESBAS algorithm compared with the three agents. All algorithms are executed on epochs of exponential size, where the updates are run only at the beginning of the epoch. Taking the three agents without ESBAS, the figure shows that the best agent in the large epochs is EEDDPG Evolution, and this is coherent with what authors mentioned in [START_REF]Evolutionary actor-multi-critic model for vnf-fg embedding[END_REF], since EEDDPG with Evolution has the possibility to learn early and quickly reaches good Acceptance Ratio. In exponential epochs, updates are less frequent, and the other agents don't perform as well as EEDDPG with Evolution. Compared with ESBAS, we can easily notice that ESBAS outperforms all the tree agents in all the big epochs starting from the third one. In Fig. 3 we present the selection decisions of ESBAS described by the algorithm usage percentage. We can notice that ESBAS makes bad selection decisions in the first small epochs, and afterward starting from the 5th epoch, it starts selecting EEDDPG and EEDDPG with Evolution. With larger epochs ESBAS makes good selection choices, and we notice that both W2FA DDPG and EEDDPG are vanishing while EEDDPG with evolution replace them more than 80% in the last epoch.

V. CONCLUSION

The problem of virtual network embedding is one of the most important challenges in B5G networks. Due to the heterogeneous nature of future networks and the unpredictable behavior of its services, this problem has been addressed using DRL algorithms. In this paper, we adopt the ESBAS framework, which takes advantage of different DRL algorithms, to solve the VNE problem. We propose an adapted version of the offline algorithm selection framework in which the selection is performed in epochs of exponential size during which the agents' updates are blocked. This technique allows the creation of stationary system to ensure a good algorithm selection. We showed the superiority of our approach over standalone DRL algorithms, and highlighted the robustness of the solution and it's ability to select the best algorithm dynamically depending on the considered scenario.

Fig. 1 .

 1 Fig. 1. Algorithm Selection Flow Diagram.

Fig. 3 .

 3 Fig. 3. Selection Ratio of ESBAS.

•

  State representation: the environment state is represented by the formulation of the VNF-FGs which are described by a vector of size |N | × K VNF + |L | × K VL . This vector describes the requests of VNFs and VLs of all VNF-FGs and is fed into the DRL agent as the environment state.

• Action: each action a t is a vector of size |N | × |N | + |L| × |L | and a t = [a n,n t , w l,l t ] where a n,n t

  • EEDDPG Evolution agent: which is the full version of EEDDPG with evolutionary operation.

	Resource types	Available Resources Requested Resources
	CPU	12-79.8	1, 2, 4, 8
	RAM	128-1024	2, 4, 8, 16, 32
	Storage	1000-10000	1, 5, 20, 100
	Bandwidth	1000-10000	10, 20, 50, 100