Frédéric Cérou
email: frederic.cerou@inria.fr

Arnaud Guyader
email: arnaud.guyader@upmc.fr

Mathias Rousset
email: mathias.rousset@inria.fr

Adaptive Multilevel Splitting: Historical Perspective and Recent Results

Keywords: Sequential Monte-Carlo, Interacting particle systems, Rare events 2010 Mathematics Subject Classification: 82C22, 82C80, 65C05, 60J25, 60K35, 60K37

About ten years ago, the Adaptive Multilevel Splitting algorithm (AMS) was proposed to analyse rare events in a dynamical setting. This review paper first presents a short historical perpective of the importance splitting approach to simulate and estimate rare events, with a detailed description of several variants. We then give an account of recent theoretical results on these algorithms, including a central limit theorem for Adaptive Multilevel Splitting. Considering the asymptotic variance in the latter, the choice of the importance function, called the reaction coordinate in molecular dynamics, is also discussed. Finally, we briefly mention some worthwhile applications of AMS in various domains.

Introduction

This paper proposes a review on the Adaptive Multilevel Splitting (AMS) algorithm to simulate rare events associated with a stochastic dynamical system. First, let us explain what we call a rare event: it is an event with non-zero, but very small probability. We assume that the probability is so small that typically we do not have any realization of the event of interest within a reasonable simulation time through a naive Monte-Carlo approach.

To give an idea, this probability will typically be smaller than 10 -10 , so that the number of simulations needed to observe only a handful of realizations is not tractable.

We can give two generic examples where it is required to precisely estimate small probabilities. First, when the rare event is some kind of catastrophe, and it is needed to know exactly how small it is, e.g. for air traffic management, or insurance. Second is when we need to estimate the mean time of return into some set for a stochastic dynamical system. This can be a transition that is not rare at the macroscopic scale, but if the dynamical system can only be simulated with a very small timestep, then that event becomes rare in the simulation timescale. We will have these two situations in mind throughout the present paper.

In the case of such a rare event, it is easy to see mathematically why a naive Monte-Carlo approach, also called crude Monte-Carlo or direct numerical simulation, is not suited. Assume we want to estimate a probability p = P(X ∈ R), with p > 0 but very low, for some random variable or process X, and a measurable set R. The naive approach is to draw an N i.i.d. sample (X 1 , . . . , X N ), with the same distribution as X, and compute the empirical probability

p = 1 N N i=1 1 R (X i ).
Since the variables 1 R (X i ) are i.i.d. Bernoulli with probability of success p, the variance of p is simply V(p) = p(1 -p)/N . If we consider the normalized variance V(p/p) = (1 -p)/(N p) ≈ 1/(N p), we see that the estimator is getting quickly worse when p goes to 0. To keep it within reasonable bounds, we would need to take N of order 1/p, which is intractable for a very low probability p.

Therefore, to construct a good estimate, one has to use some variance reduction technique. For this problem, there are broadly two families of solutions: one is Importance Sampling, and the other one is Importance Splitting. The latter is the main subject of this paper, but let us write a few words on the former.

The idea of Importance Sampling is to sample from an auxiliary distribution in order to make the rare event less rare. Let us denote by (Y 1 , . . . , Y N ) the new i.i.d. sample. One usual requirement is that the distribution of X is absolutely continuous w.r.t. the law of Y , and we denote the corresponding Radon-Nikodym derivative by dP X dP Y . We then estimate the probability p by

p = 1 N N i=1 1 R (Y i ) dP X dP Y (Y i ).
This estimator is clearly unbiased. Nonetheless, it is not always obvious to choose a good importance distribution to sample from, that is one which will give an estimator with a small normalized variance. In fact, if the sampling distribution is badly chosen, things can be very bad, as one is not even guaranteed to have a finite variance (see for example [START_REF] Glasserman | Counterexamples in importance sampling for large deviations probabilities[END_REF]).

Importance Sampling is a very common approach to reduce variance, and there is a huge amount of literature on the subject, that we will not discuss here. Let us just mention the monograph [START_REF] Bucklew | Introduction to rare event simulation[END_REF] in the context of rare events, and [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF] for illustrations in molecular dynamics.

The idea of Importance Splitting, that is the family AMS belongs to, is to simulate according to the original distribution in a sequential way, to discard the trajectories (or samples) going far away from R, and to split/branch/clone those that get closer. This will be made more precise in the sequel. Ä Before proceeding, let us just mention that, in what follows, we will focus our attention on the so-called "dynamic case". By this, we mean that the rare event of interest writes p = P(X τ ∈ R), with X a strong Markov process and τ a stopping time (see Section 2.1). The so-called "static case" corresponds to the situation where, typically, X is a random vector in R d , S is a function from R d to R, and one wants to estimate the probability p = P(S(X) > q) where q is known and such that p is strictly positive but very low. There are of course some strong connections between multilevel splitting methods for the dynamic case and the static case, but the results obtained in both situations are not exactly the same ones. However, even if we have chosen to focus on the dynamic case, we will also mention here and there some results available in the static case.

Finally, note that the key word "multilevel" in Adaptive Multilevel Splitting is completely unrelated to the "Multilevel Monte Carlo" developed by Giles (see, e.g., [START_REF] Giles | Multilevel monte carlo path simulation[END_REF] or [START_REF] Giles | Multilevel monte carlo methods[END_REF]) although there are some attempts in the literature to use both ideas in conjunction (see for example [START_REF] Ullmann | Multilevel estimation of rare events[END_REF]).

A little history 2.1 The origins

Interestingly enough, one of the first presentations of both Importance Sampling and Importance Splitting is generally thought to be [START_REF] Kahn | Estimation of particle transmission by random sampling[END_REF]. Their original problem came from particle physics, where the goal is to simulate and study the particle transmission through an obstacle. Of special interest is the probability that the particle goes through the obstacle without being absorbed. If, for instance, a nuclear device is to be considered safe, that probability should be extremely low.

Let us state the problem in a more abstract form. Let X be a time homogenous strong Markov process in R d starting at t = 0 with a known distribution η 0 . Let F X t = σ(X s , 0 s t) denote its natural filtration that we will assume for simplicity to be right-continuous. These assumptions ensure that the hitting time τ of any (measurable) set is a stopping time, and that the law of X conditional on the past of τ is again the law of X with initial condition X τ .

Note that, in many cases, one may assume that we have preliminarily added the time in the state space (i.e., one coordinate of the state X is the time variable). Then there is no loss of generality in considering time inhomogenous processes and observables that may depend on time.

We suppose that τ is a.s. finite, which includes the case of a deterministic time. The problem is then to construct an estimator of

p = P(X τ ∈ R),
and to give a sample of trajectories such that its empirical measure is an approximation of the distribution of (X s , 0 s τ ) given that X τ ∈ R.

To be more specific, let us give two different instances of the latter. For the original problem in [START_REF] Kahn | Estimation of particle transmission by random sampling[END_REF], we have a random killing time σ a (absorption of the particle), and R is the outside of a confinement tank. If we denote σ R the hitting time of the outside, then τ = σ a ∧ σ R .

Another typical setting arises in molecular dynamics, detailed below in Section 4.4. Here X is a diffusion process, with a drift that derives from a potential, and a constant (full rank) noise intensity. We want to simulate reactive trajectories, which are trajectories that reach a region R (typically another well in an energy landscape), before visiting a recurrent set (typically a neighborhood of the bottom of the current well of the potential). If we denote σ r the latter stopping time, we have in that case τ = σ R ∧ σ r .

The basic idea is to use a real-valued function Φ, also called an importance function or a reaction coordinate in molecular dynamics, to measure, even roughly, how "close" the process may be from the rare event. Ideally, imagine that Φ(x) = Φ (x) = q (x)/p where

q (x) = P(X τ ∈ R|X 0 = x)
is the so-called committor function, which we will extensively discuss later (to write the latter in all generality, we need to assume that the time t is a component of the state vector x, but this is not necessary for the two examples given just above). Then each time a sample trajectory crosses a surface {Φ (x) = 2 k }, for k ∈ {1, . . . , log p log(1/2) }, we split, or clone that trajectory in 2, that is afterwards, we simulate 2 independent trajectories, issued from the same passage point into the level set {Φ (x) = 2 k }. Of course we also divide by 2 the mass of each cloned trajectory. If we do that until step k = log p log(1/2) , the remaining paths have a probability larger than 1/2 of hitting R. The latter probability can therefore be estimated by crude Monte-Carlo. Moreover, the number N of trajectories is random, but its expectation is constant and equal to N 0 , the size of the initial sample.

There is a huge difficulty when using this simple method on a practical example, because the function Φ (x) = P(X τ ∈ R|X 0 = x)/p is unknown, and its computation is intractable, otherwise the rare event problem would be solved. In practice a function Φ which is imprecise is used, often resulting from heuristics proposed by practitioners. By "imprecise", we mean that best asymptotic variance for the AMS estimator of p is reached when taking Φ = Φ , which unfortunately is impossible in most situations of interest. This crucial point will be discussed in Section 3.5.

Note also that the number of clones 2 was suggested by the authors in [START_REF] Kahn | Estimation of particle transmission by random sampling[END_REF], but of course any other value could be chosen: if it is not an integer, we can randomise the number of clones in order to keep the mean value constant. The Multilevel Splitting (MS) algorithm with an integer number of clones (to keep it simple) is detailed in Algorithm 1, and illustrated in Figure 1. We just assume that we have a continuous function Φ and a real number L max such that, defining the stopping time

S Lmax = inf {s 0, Φ(X s ) > L max } , we have the condition X τ ∈ R ⇒ τ S Lmax ,
which precisely ensures that the trajectories contributing to the target event {X τ ∈ R} have all reached the open set {Φ(x) > L max }. The choice of an open set {Φ(x) > L max } is merely conventional, and the case {Φ(x) L max } could be treated similarly. In the same way, the choice {Φ(X i s ) > L j } instead of, say {Φ(X i s ) L j }, in the definition of the thresholds is conventional. The only important condition is that the successive thresholds constitute a decreasing sequence of sets for inclusion.

Algorithm 1 Multilevel Splitting (MS)

Require: Initial distribution η 0 , Importance Function Φ and levels

-∞ = L 0 < L 1 < • • • < L J = L max
for a given number of levels J, cloning rates r 1 , . . . , r J , initial sample size N 0 Initialization: X 1 0 , . . . , X N 0 0 i.i.d. from η 0 for j = 1 to J do for i = 1 to N j-1 do Run trajectory i until next level {Φ(X i s ) > L j } or final time τ i (the first reached) end for Discard trajectories that did not reach L j Clone r j times those which did Denote N j the total number of resulting trajectories Reorder the trajectories with number from 1 to N j end for for i = 1 to N J do Run trajectory i until final time τ i end for Estimate the probability of the rare event by

pms = 1 N 0 1 J j=1 r j N J i=1 1 R (X i τ i )
A crucial choice in Algorithm 1 is the branching rates r j . Let us denote by S L j the hitting time of the threshold {Φ > L j }:

S L j = inf{s, Φ(X s ) > L j },
and by p L j = P(S L j τ ) the associated probability. We will also denote

θ j = p L j /p L j-1 = P(S L j τ |S L j-1 τ ). L max L 2 L 1 X 0 r 1 = 3 r 2 = 2 0 N 0 = 2 Figure 1: Algorithm 1 (MS)
, where X is to reach L max before 0.

To have a useful algorithm, we need that the products r j θ j are all close to 1.

If they are much smaller, there is a high probability that no trajectory hits R (like in naive Monte-Carlo); if they are much larger, the algorithm simulates a very large number of correlated trajectories. An analysis of the cost/variance trade-off of this approach can be found in [START_REF] Lagnoux-Renaudie | A two-step branching splitting model under cost constraint for rare event analysis[END_REF][START_REF] Lagnoux | Multilevel branching and splitting algorithm for estimating rare event probabilities[END_REF] for discrete dynamical models. In a nutshell, to get the maximum efficiency, the probabilities from one level to the next should all be equal, and the r j θ j should be equal to 1.

Also note that a version of the algorithm suitable for computing µ ∞ (R), where µ ∞ is the invariant probability measure of the process X, is proposed in [START_REF] Villén-Altamirano | RESTART : a method for accelerating rare event simulation[END_REF][START_REF] Villén-Altamirano | Enhancement of the accelerated simulation method RESTART by considering multiple thresholds[END_REF] and subsequent papers by the authors of [START_REF] Villén-Altamirano | RESTART : a method for accelerating rare event simulation[END_REF]. Added to the splitting mechanism, it also provides a pruning rule for the trajectories that go downwards in terms of the importance function Φ, keeping only one trajectory all along the time axis. The use of pruning mechanisms to speed up the algorithm even when not considering the invariant probability measure is discussed in [START_REF] L'ecuyer | Splitting for rare-event simulation[END_REF], which also provides a discussion on the early variants of importance splitting.

As we will see below, we can modify the algorithm so that we do not need to choose the branching rates a priori, and get a tractable algorithm without any real drawback.

Further developments

The fact that the number of active trajectories in Algorithm 1 is random makes it difficult to use. Moreover, if the branching rates are not well chosen, it becomes also inefficient. This is why it was proposed in [START_REF] Garvels | The splitting method in rare event simulation[END_REF][START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF][START_REF] Cérou | Genetic genealogical models in rare event analysis[END_REF] a variant that keeps the number of active trajectories (or replicas) constant, which is detailed here as Algorithm 2 and Figure 2. In this algorithm, one may have at step j the equality N j = N , in which case no splitting occurs (see also the discussion in Section 3.4). Moreover, in the latter again, the choice {Φ(X i s ) > L j } instead of, say {Φ(X i s ) L j }, in the definition of the thresholds is just conventional.

Algorithm 2 Sequential Monte-Carlo (SMC) Require: Initial distribution η 0 , Importance Function Φ and levels -∞ = L 0 < L 1 < • • • < L J = L max for a given number of levels J, sample size N Initialization: X 1 0 , . . . , X N 0 i.i.d. from η 0 for j = 1 to J do for i = 1 to N do Run trajectory i until next level {Φ(X i s ) > L j } or final time τ i (the first reached)
end for Discard trajectories that did not reach L j Set N j the number of remaining trajectories, and I j the set of their indices for i ∈ {1, . . . , N } \ I j do Choose uniformly at random an index in I j , clone it, and replace former trajectory i by it end for end for Estimate the probability of the rare event by

psmc = J j=1 N j N × 1 N N i=1 1 R (X i τ i )
This version is sometimes referred to as Fixed Effort (in opposition to Fixed Splitting, see [START_REF] Garvels | The splitting method in rare event simulation[END_REF]) but we propose to call it the Sequential Monte-Carlo (SMC) variant of multilevel splitting. Indeed, it is in fact a special case of a very general particle method usually named Sequential Monte-Carlo (see for example [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF]). It was originally designed for non-linear filtering problems arising in Signal Processing, but has been reshaped as a very powerful abstract framework related to discrete time Feynman-Kac formulae in the monographs [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF][START_REF] Moral | Mean field simulation for Monte Carlo integration[END_REF]. We will see in Section 3 how these general results can be used to deduce mathematical properties of Algorithm 2.

L max L 2 L 1 X 0 0 N = 3 Figure 2: Algorithm 2 (SMC)
, where X is to reach L max before 0.

Recall that, in the static case, X is a random vector in R d , S is a function from R d to R, and one wants to estimate the probability p = P(S(X) > q) where q is known and such that p is strictly positive but very low. The corresponding algorithm in this context was proposed and studied by several others, see for example [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] and [START_REF] Botev | An efficient algorithm for rare-event probability estimation, combinatorial optimization, and counting[END_REF]. More recently, the authors of [START_REF] Gobet | Rare event simulation using reversible shaking transformations[END_REF] propose a framework that encompasses both the static and the dynamic case.

Adaptive Multilevel Splitting

To make the algorithm even more easy to use, the next idea is to set the levels adaptively. Instead of choosing the levels to be reached, one chooses the number K of trajectories to discard out of N , which is equivalent to choosing the proportion (N -K)/N to be kept. The level used at each iteration is then an empirical quantile of the maximal levels reached by the set of trajectories. The resulting algorithm, known as Adaptive Multilevel Splitting, is given as Algorithm 3 and illustrated in Figure 3. It was first mentioned in [START_REF] Garvels | The splitting method in rare event simulation[END_REF], and then formalized and studied in dimension 1 in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF].

An argument to go from Algorithms 1 and 2 to Algorithm 3 is to consider the variance of a standard splitting approach, and to remark that the optimal way to choose the intermediate levels is to fix them such that the successive conditional probabilities are constant. This point is detailed in Section 3.5.

Note that depending on the choice of K, which can depend on N , we can have two interesting regimes. First in the case where (N -K)/N = θ ∈ (0, 1) (see [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]), and another regime when

K = 1, meaning that (N -K)/N = 1 -1/N L max X 0 0 N = 3 1st branching 2nd branching K = 1 Figure 3: Algorithm 3 (AMS)
, where X is to reach L max before 0.

(see [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF]). The case K = N is not interesting because it is equivalent to a naive Monte-Carlo method. Also notice that if, at some stage m, one has

I m = {1, . . . , N }, meaning that M (K) = . . . = M (N )
, then the algorithm is stuck and there is a failure. However, there are some ways to avoid this pathological situation (see Remark 2.1 below).

It turns out that one can adapt easily the AMS method to the static case. The case (N -K)/N = θ ∈ (0, 1) is then often referred to as Subset Simulation (see for example [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF][START_REF] Au | Subset simulation and its application to seismic risk based on dynamic analysis[END_REF] and [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF][START_REF] Cérou | Fluctuation analysis of adaptive multilevel splitting[END_REF]). The case (N -K)/N = 1 -1/N was introduced and studied in [START_REF] Guyader | Simulation and estimation of extreme quantiles and extreme probabilities[END_REF] and shares some connections with Nested Sampling (see [START_REF] Skilling | Nested sampling for general Bayesian computation[END_REF][START_REF] Skilling | Nested sampling for Bayesian computations[END_REF] and [START_REF] Chopin | Properties of nested sampling[END_REF]).

Coming back to Algorithm 3, a first remark is that the number of iterations is random, but that does not mean the computing time is out of control. Indeed, if (N -K)/N → θ ∈ (0, 1), then it is easily checked formally that the number of iterations converges in probability when N → ∞ to

log p Lmax log θ
. This was proved rigorously in dimension 1 in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF], and in the static case in [START_REF] Cérou | Fluctuation analysis of adaptive multilevel splitting[END_REF], Theorem 3.1, whose proof can easily be adapted to the present framework.

When (N -K)/N → 1, some work has been done when K = 1, the so-called "last particle" case, which corresponds to the maximum number of iterations. Then it is shown in [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF], for X a uniformly elliptic diffusion, that the number of iterations is of order O P (-N log p). Hence it grows linearly with N (as usual) but only logarithmically with p, which is a nice property when p is very low.

A second remark is that, in Algorithm 3, the number K m of discarded trajectories becomes random when more than one trajectory have the same current Algorithm 3 Adaptive Multilevel Splitting (AMS) Require: Initial distribution η 0 , Importance Function Φ, sample size N , minimal number K of trajectories to discard at each step, final level L max Initialization:

X 1 0 , . . . , X N 0 i.i.d. from η 0 Set m ← 0 (iteration index) for i = 1 to N do
Run each trajectory to its end

τ i Set M i ← max 0 s τ i Φ(X i s ) end for Sort the M i 's from low to high, so that M (1) . . . M (N ) Set current level L ← M (K) while L < L max do m ← m + 1
Discard all the trajectories for which M i L Let K m be the number of such trajectories (hence K m K) Define I m as the set of indices of the remaining trajectories for i ∈ {1, . . . , N } \ I m do Choose uniformly at random an index in I m Clone the corresponding trajectory until the first time it enters {Φ > L} From that time, simulate the cloned trajectory up to its end time τ i Replace the trajectory with index i by that new trajectory Set M i ← max 0 s τ i Φ(X i s ) end for Sort the M i 's from low to high, so that M (1) . . . M (N ) Set current level L ← M (K) end while Set M = m the total number of iterations Estimate the probability of the rare event by

pams = M m=1 N -K m N × 1 N N i=1 1 R (X i τ i )
score L. If we simplify and discard exactly K at each iteration, the algorithm is still usable, but unsuited for parallelization due to the presence of a (small) bias. This will be discussed in more details in Section 3.2.

Remark 2.1. [A resampling variant]

There is a trick which enables to treat differently the equality case in the order statistics M (1) . . . M (N ) in Algorithm 3 which may result in K m > K. The variant is as follows:

• First impose arbitrarily a total order among the particles' scores M [START_REF] Amrein | A variant of importance splitting for rare event estimation: Fixed number of successes[END_REF] . . . M (N ) by choosing it uniformly at random among the compatible ones.

• Second, discard the K first trajectories with respect to that order, instead of the K m K trajectories with a lower or equal score.

• Third and last, the duplicated remaining particles trajectories are cloned until either the first time they enter {Φ L}, or the first time they enter {Φ > L} (if the latter exists). Note that a decision rule has to be chosen here.

The probability of the rare event is now estimated by

M m=1 N -K N 1 N N i=1 1 R (X i τ i ).
This variant, with an appropriate decision rule in the third step, is expected to obey similar theoretical results (unbiasedness, consistency, CLT), although no rigorous proof is currently available. One advantage of this variant is that it usually prevents extinction, the latter being replaced by more trials attempting to strictly increase the current level.

One should also remark that this variant has been numerically tested in Section 5.1 of [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF] (the so-called Version 2), in the case where the remaining particles trajectories are cloned until the first time they enter {Φ L}. The authors observe that the associated estimator of the rare event probability is then biased. For other, unbiased variants which prevent extinctions, we refer the interested reader to Section 3.5.1 in [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF].

3 Mathematical results

Some notation

We will denote S the hitting time of the open threshold , that is

S = inf{s, Φ(X s ) > }.
For simplicity in the presentation, from now on we will assume that the rare event is exactly

X τ ∈ R ⇐⇒ S Lmax τ, (3.1)
so that the trajectories contributing to the rare event are exactly those who reach the threshold {Φ > L max }. In particular, this implies that the rare event probability becomes p = p Lmax = P(S Lmax τ ).

We will also need the following "committor-like" function q(ϕ), defined for any test function ϕ by

q(ϕ)(x) := E[ϕ(X τ )1 R (X τ )|X 0 = x]. (3.2)
Note that, by construction, the committor function mentioned earlier is just q = q(1).

Let us then define the probability measure η as the distribution of X τ , given that X τ ∈ R. That is, for every test function ϕ,

η(ϕ) = 1 p E[ϕ(X τ )1 R (X τ )] = 1 p E[ϕ(X τ )1 S Lmax τ ] = 1 p η 0 (q(ϕ)).
We also introduce the unnormalized measure

γ(ϕ) = p η(φ) = E[ϕ(X τ )1 R (X τ )] = E[ϕ(X τ )1 S Lmax τ ] = η 0 (q(ϕ)).
Let us denote pms , psmc and pams the estimates of p = P(X τ ∈ R) for the previous algorithms, meaning respectively Multilevel Splitting (Algorithm 1), Sequential Monte-Carlo (Algorithm 2) and Adaptive Multilevel Splitting (Algorithm 3). Note that condition (3.1) simplifies formulas since 1 R (X i τ i ) = 1 almost surely in each estimator. Let us also denote (X 1 ms , . . . , X N J ms ) and ηms = 1

N J N J i=1 δ X i ms
the set of particles at the end of Algorithm 1 and the corresponding empirical measure. Accordingly, γms = pms ηms stands for the unnormalized empirical measure. Note that γ(1) = p and γms (1) = pms . We define ηsmc , γsmc , ηams , γams in the same manner.

A few words about the bias

It turns out that the estimates pk of the probability p for k ∈ {ms, smc, ams}, and more generally the unnormalized measures γk (ϕ) for any test function ϕ, are all unbiased. For Algorithm 1 it is quite straightforward, see [START_REF] Amrein | A variant of importance splitting for rare event estimation: Fixed number of successes[END_REF].

For Algorithm 2, it is a general property of Sequential Monte-Carlo methods, as expected from a Monte-Carlo method used to compute an expectation, here in the form of an unnormalized Feynman-Kac measure (see for example [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF] page 112). The key point is to remark that γsmc (ϕ) can be expressed as the terminal value of a martingale with initial value given by the empirical average of q(ϕ) obtained after the initialization, that is 1

N N i=1 q(ϕ)(X i 0 ) which is obviously unbiased.
For Algorithm 3, the property is a consequence of a very general result given in [START_REF] Bréhier | Unbiasedness of some generalized adaptive multilevel splitting algorithms[END_REF]. The discussion in Section 3.4 of the present paper explains (without rigorous proof) that Algorithm 3 can be interpreted as a limit of Algorithm 2 when the number J of levels goes to infinity. This suggests that the unbiasedness property is inherited by Algorithm 3. Alternatively, a specific martingale analysis could be performed as done in a diffusive case in [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF]. In Algorithm 3, it is important to take care of the multiple values for the maximum of Φ along the set of trajectories. If no care is taken, then the resulting algorithm might be biased. Note that in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF] no such result was obtained precisely because no care was taken of these multiple values.

The variant in Remark 2.1 is also expected to yield unbiased estimators of unnormalized quantities, although up to now there exists no rigorous proof for this variant. The equality case requires to handle in some way a martingale indexed by two indices, endowed with lexicographic order.

It is known that the normalized estimates ηk (ϕ) are all biased with a bias of order 1/N for k ∈ {smc, ams}. The estimate ηms (ϕ) is expected to be biased as well but has still not been carefully studied.

If we want to parallelize the computation, and we are only interested in the estimate of the probability p, the unbiasedness property allows us to run several independent versions of the algorithm in parallel, and then take the empirical mean of all the resulting estimates. If we want a conditional mean of an observable, then we have to parallelize within each algorithm. In Algorithm 3, we can use for example K as the number of trajectories we can resimulate in parallel. Note that they will not have the same length, so they all need to wait for the longest one to finish.

Sequential Monte-Carlo (discrete levels)

This section details theoretical results for Algorithm 2. Let us recall that S L j stands for the hitting time of L j , i.e. S L j = inf{s, Φ(X s ) > L j } where the successive levels -∞ = L 0 < L 1 < • • • < L J = L max are given a priori. Recall that p L j = P S L j τ and

θ j = p L j /p L j-1 = P(S L j τ |S L j-1 τ ).
First we need to mention that there is a small but non zero probability of extinction, that is all the particles fail to reach some level L j . In that case, we consider that we estimate the probability p as 0, just as we would do in naive Monte-Carlo. Fortunately, the probability of such a failure is soon very small, meaning with a reasonable number N of particles. A simple version of Theorem 7.4.1 in [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF] is given in the following proposition. Recall that, in all the paper, we assume that p > 0.

Proposition 1. Define the extinction event:

E := {Algorithm 2 fails with an extinction}.

There exist two constants (depending on the problem being solved) A > 0 and B > 0 such that for all N 1, P(E) Ae -BN .

Note that there are versions of the algorithm without extinction (see [START_REF] Gland | A sequential particle algorithm that keeps the particle system alive[END_REF] for a biased version, [START_REF] Amrein | A variant of importance splitting for rare event estimation: Fixed number of successes[END_REF] for an unbiased version, and [START_REF] Del Moral | The Alive Particle Filter and Its Use in Particle Markov Chain Monte Carlo[END_REF] for a generalization to non-negative potentials), but then the computation time is not bounded.

We also have a law of large numbers which ensures the convergence of the algorithm, as well as a Central Limit Theorem. The asymptotic variance in the latter is of special interest. Its definition requires to introduce the following conditional distributions:

η L j (ϕ) := E[ϕ(X S L j )|S L j τ ].
Recall that the committor function is defined by q (x) = P(X τ ∈ R|X 0 = x).

Theorem 3.1. The estimator psmc satisfies

1 E c √ N (p smc -p) Law ---→ N →∞ N (0, σ 2 )
where

σ 2 = J-1 j=1 θ j p 2 L j-1 -p 2 L j V η L j (q ) + p 2 J j=1 1 θ j -1 . (3.3)
A similar CLT can be obtained for the final empirical distribution

ηsmc = 1 N N i=1 δ X i τ i
which estimates the conditional distribution η = L(X τ |X τ ∈ R), as obtained at the end of Algorithm 2 (see Section 9.4.2 of [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF]):

1 E c √ N (η smc (ϕ) -η(ϕ)) Law ---→ N →∞ N (0, σ 2 (ϕ -η(ϕ)))
with an asymptotic variance given for any bounded observable ϕ by

σ 2 (ϕ) = V η (ϕ) + J j=1 θ j p 2 L j -p 2 L j-1 p 2 V η L j (q(ϕ -η(ϕ))). (3.4)
Using the same technical apparatus (see Section 9.4.2 of [START_REF] Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF]), the latter CLT can be extended to path observables of the empirical distribution

ηpath = 1 N N i=1 δ X i 0 s τ i (3.5)
which estimates the law η path of the full trajectory conditional on the rare event. The generalized pathwise variance is the same as in (3.4), except that the observable dependent committor function

V η L j (q(ϕ)) = V E ϕ(X τ )1 S Lmax τ |X S L j .
is now replaced by the variance of the analogous path functional

V E ϕ(X 0 s τ )1 S Lmax τ |X 0 s S L j . (3.6)

Estimate of the variance

Recently, [START_REF] Lee | Variance estimation in the particle filter[END_REF] proposed and showed the convergence of an estimator of the variance applicable to a range of SMC algorithms, including particle filter.

The theory applies to a modification of Algorithm 2, in which the resampling of particle is instead a multinomial resampling. By that we mean that we do not keep necessarily all the successful particles at the next level, but make N i.i.d. uniform draws among them. Of course, this will lead to more dependence than Algorithm 2, and thus more variance. We nevertheless think that even if the theory does not apply strictly to Algorithm 2, this variance estimate can in practice give a useful insight. But clearly, more investigations are needed here.

If we restrict it to our rare event setting, with Algorithm 2, we consider

V N J (ϕ) = η N J (ϕ) 2 - N J-1 (N -1) J+1 A i 0,J =A k 0,J ϕ(X i J )ϕ(X k J ),
where A i 0,J is the ancestor at time 0 of X i J . Then N V N J (ϕ) is a consistent estimator (as N → ∞) of σ 2 (ϕ), and N p2 smc V N J (1) a consistent estimator of σ 2 . This estimator is actually an unbiased modification of the one proposed earlier in [START_REF] Chan | A general theory of particle filters in hidden Markov models and some applications[END_REF].

These estimators are almost free to compute because we only need to store the ancestors of the particles all along the algorithm. The drawback is that if N is not large enough, all the particle will share a very small set of ancestors (or even the same one), and the estimators will not be useful.

Adaptive Multilevel Splitting (continuous levels)

Let us now discuss some theoretical results related to Algorithm 3.

Well-posedness and CLT for the case K = 1

For the case K = 1 ("last particle" version), some results have been obtained under additional assumptions. For simplicity, we will only explain here the diffusive case. We refer the interested reader to [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF] for details and complements.

Suppose that the process (X s ) s 0 is a strong solution of a Stochastic Differential Equation

dX s = b(X s ) + σ(X s )dW s ,
where b and σ are smooth. Moreover, we assume that the level function Φ is smooth and there exists δ > 0 such that (∇Φ) T a∇Φ δ.

In this context, Algorithm 3 is well-defined in the sense that, almost surely, only one particle is discarded at each step (no equality in the scores) and the algorithm stops after a finite number of steps, meaning that a.s. M < ∞.

We then have a law of large numbers which ensures the convergence of the algorithm when N → ∞, as well as a Central Limit Theorem. Interestingly, as will be explained later, the asymptotic variance in Theorem 3.2 can be seen as the "continuous levels" limit of the Sequential Monte-Carlo variance ("discrete levels" case) as given in (3.3). More precisely, let us write the conditional probability

η (ϕ) := E[ϕ(X S )|S τ ].
In the same way we denote

p = P(S τ ).
We stress that, as before, is not an integer, but a real number corresponding to a level. In this respect, we recall that p = p Lmax is assumed strictly greater than 0.

Theorem 3.2. The unbiased estimator pams satisfies the CLT

√ N (p ams -p) Law ---→ N →∞ N (0, σ 2 )
where

σ 2 = -p 2 ln p- Lmax -∞ V η (q )d p 2 = -p 2 ln p-2 Lmax -∞ V η (q )p dp . (3.7)
In these formulas, the integration is with respect to the level , meaning that goes from -∞ to L max . Also note that in the expression of the asymptotic variance σ 2 , both terms are positive since 0 < p < 1, and → p is decreasing, making dp negative. As in the discrete case (SMC), a CLT can be obtained for the final empirical distribution at the end of Algorithm 3, which we also denote

ηams = 1 N N i=1 δ X i τ i .
The latter estimates the conditional distribution η = L(X τ |X τ ∈ R). We prove in [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF] the following CLT:

Theorem 3.3. If ϕ is bounded and continuous, then √ N (η ams (ϕ) -η(ϕ)) Law ---→ N →∞ N (0, σ 2 (ϕ -η(ϕ))) with σ 2 (ϕ) = V η (ϕ) - Lmax -∞ V η (q(ϕ))d p 2 /p 2 . (3.8)
Here again, the latter CLT can been extended to pathwise bounded and continuous -with respect to uniform convergence on compact time intervals -observables of the empirical distribution

ηpath = 1 N N i=1 δ X i 0 s τ i ,
which estimates the law η path of the full trajectory conditional on the rare event. The generalized pathwise variance is the same as (3.5), except that the variance of the observable dependent committor function V η l (q(ϕ)) is again replaced by (3.6). The extension of the CLT to non continuous observables remains open.

Interpretation as a limit of the discrete levels case when K = 1

Now we would like to discuss a topic that has not been investigated in the literature so far, namely the fact that Algorithm 3 (AMS), which has been introduced as an adaptive version of the discrete levels Algorithm 2 (SMC), can in fact be understood as a continuous levels limit of the latter.

For this purpose, let us assume that the discrete levels are chosen of the form

L j = L 1 + (j -1) L max -L 1 J -1 and Φ(X 0 ) L 1 a.s.
Then, modify slightly the formulation of Algorithm 2 by simulating in the main loop the whole path -up to the final stopping time τ i -of each newborn particle i, instead of simulating it level brackets by level brackets. Note that this modification does not change the probability distribution of the whole particle system.

It is then easy to check (e.g., with a drawing) that Algorithm 3 is exactly Algorithm 2 in the case where the function Φ takes its values in a finite set, for instance using Φ ε approximated by

Φ ε = ε Φ/ε . (3.9)
For J large compared to N , most iterations in Algorithm 2 become useless (nothing happens), except when the level value L j coincides with the value of the smallest particle score within the current particle system. Note that iterations in Algorithm 3 exclusively select the latter events: this explains why the iteration index was denoted with a different letter, namely m instead of j.

Finally, for instance when X is a pure jump process which takes only a finite number of values, it is clear that the latter algorithm will be the same if performed with any Φ ε , provided that ε is small enough. This implies that, at least formally, Algorithm 3 is the "limit" when J → ∞ of Algorithm 2.

It is then possible to compute formally the J → ∞ limit in the asymptotic variances obtained in the previous CLTs. Indeed, suppose for simplicity that p L 1 = 1, then one has

J j=1 1 θ j -1 = J j=1 - p L j -p L j-1 p L j → - Lmax L 1 dp p = -ln p Lmax = -ln p.
In the same way, for any continuous function f ,

J j=1 f (L j )θ j p 2 L j-1 -p 2 L j → - Lmax L 1 f ( ) d(p 2 ).
We thus recover the continuous levels variances (3.7) and (3.8) from the limits, when J → ∞, of the discrete levels variances (3.3) and (3.4). In turn, this suggests that the CLTs of Theorems 3.2 and 3.3 are in fact valid in a completely general setting, as soon as the extinction probability becomes small enough when N → ∞.

Finally, note that a similar interpretation of Algorithm 3 as a limit of Algorithm 2 in the case where K = 1 is possible, but requires a substantial modification of Algorithm 2. In such a modified algorithm, the number of particles a each iteration is non constant and decreasing -particles are not duplicated -until at least K particles are killed, which triggers the creation of new particles to recover the maximal population size N .

About the asymptotic behavior for the case K/N → θ > p Interestingly enough, a kind of converse remark can be made: in the case where 1 -K/N → θ ∈ (p, 1) when N → ∞, Algorithm 3 is very similar to the discrete levels case of Algorithm 2, for the specific choice of equiprobable levels, meaning that the levels L j , j ∈ {1, . . . , J -1}, are such that

     J = log p/ log θ + 1 θ j = θ for j = 1, . . . , J -1 θ J = p/θ J-1 ∈ (θ, 1]. (3.10)
This point view has been systematically studied in [START_REF] Cérou | Fluctuation analysis of adaptive multilevel splitting[END_REF] where a CLT is obtained in the static case. However, the assumptions therein, in particular the use of multinomial resampling, do not apply to Algorithm 3.

Nonetheless, the relationship between the two algorithms can be described assuming at least that for each j ∈ {1, . . . , J -1}, there exists a unique level L j such that P(τ S L j ) = p L j = θ j . In that case, the random level reached at iteration m J-1 of Algorithm 3, denoted Lm , should converge in probability towards the deterministic level L m associated with the probability θ m . As a consequence, at each iteration in both algorithms, a fraction K/N (or an empirical estimation of it for Algorithm 2) of particles fail to reach the next level in the sequence L 1 < . . . < L J (or in an empirical estimation of it for Algorithm 3), both algorithms being thus formally similar when N → +∞.

Finally, the results in [START_REF] Cérou | Fluctuation analysis of adaptive multilevel splitting[END_REF] suggest that a CLT may be obtained for pams and ηams with exactly the same variances as in (3.3) and (3.4) in the case (3.10). Although these results are strongly believed to be true, they still require a rigorous proof.

On the importance function

In practice, the main source of variance in Algorithms 2 and 3 comes from a bad choice of the importance function Φ (reaction coordinate). For this reason, it is of crucial interest to try to minimize the asymptotic variances (3.3) or (3.7) with respect to the choice of Φ.

For simplicity, we assume that the initial condition is deterministic X 0 = x 0 , and that Φ is at least continuous. We can then remark that, on the one hand, the target open set {x, Φ(x) > L max } depends on Φ only through its boundary {x, Φ(x) = L max } ⊂ R d , while, on the other hand, Algorithms 2 and 3 depend on all the intermediate level sets {x, Φ(x) = } from = L min = Φ(x 0 ) to = L max .

This implies that the latter algorithms are unchanged when Φ is multiplied by a constant, so we can assume without loss of generality that L max = 1. More importantly, the target set is unchanged if Φ is modified while keeping the set {x, Φ(x) = L max } fixed. As a consequence, we will now try to optimize σ 2 in the set of continuous Φ 1 such that Φ = 1 on the latter set.

We first start with the AMS variance as given by (3.7) (continuous levels case), and immediately remark that the term -p 2 ln p cannot be modified while the other term vanishes as soon as the following condition holds true: .11) In the diffusive case -and, interestingly, only in the diffusive case -, that is when (X s ) s 0 has continuous trajectories, the condition is satisfied for an explicit choice of Φ which turns out to be (any function of) the committor function Φ = f (q ) ⇒ V η (q ) = 0, where f is any continuous strictly increasing function. Indeed, by continuity of the trajectories, the conditional distribution η has necessarily for each its support in the associated level set

V η (q ) = 0 ∀ ∈ [L min , L max ]. ( 3
η (1 {Φ= } ) := P(Φ(X S ) = |S τ ) = 1, so that V η (f -1 (Φ)) = 0.
Next, for the SMC asymptotic variance (3.3) (discrete levels case) and the choice Φ = f (q ) for the importance function, the remaining variance term reduces to

p 2 J j=1 1 θ j -1 ,
which, in turn, can be optimized for each fixed J. Indeed, the only constraint is θ 1 × . . . θ J = p, so that a standard convex optimization implies that the optimum is reached for θ j = p 1/J for all j. As a consequence, the minimal possible variance is

σ 2 = p 2 J p -1/J -1 ,
which is decreasing towards -p 2 ln p when J → ∞, i.e. when θ j → 1. This accounts for the choice K = 1 in order to minimize the variance.

This simple remark also explains why in the AMS Algorithm 3 (with no equality in the scores for simplicity), the number K of particles which is discarded is the same at each step. In full generality, one could decide to discard K 1 particles at the first step, then K 2 at the second step, etc., run the same algorithm and consider the estimator

pams = M m=1 N -K m N × 1 N N i=1 1 R (X i τ i ).
But the constrained optimization problem above proves that, in the ideal situation where one would have the committor function at disposal, the best thing to do is to discard the same number of trajectories at each step. In other words, the idea is to minimize in (θ 1 , . . . , θ J ) the lower bound for σ 2 given by (3.3).

One can also notice that the asymptotic variance (3.8) for the conditional empirical distribution can be bounded as follows (see Section 2.4 in [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF])

V η (ϕ) σ 2 (ϕ -η(ϕ)) V η (ϕ) + ϕ -η(ϕ) 2 ∞ σ 2 p 2 -ln p ,
with σ 2 as in (3.7). The lower bound is the variance we would get with an i.i.d. sample from η. As noticed above, at best the second term in the r.h.s. reduces to -2 ϕ -η(ϕ) 2 ∞ ln p. Concerning the asymptotic variance σ 2 of pams in Theorem 3.2, we can also show (see Corollary 2.10 in [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF]) that one always has -p 2 log p σ 2 2p(1 -p).

In other words, it is impossible with AMS to do better than -p 2 log p, even with the optimal importance function at disposal, but it is also impossible to do worse than 2p(1 -p), which is twice the variance of a naive Monte-Carlo method (see Section 1). In comparison, it is well-known that for Importance Sampling, there always exists an optimal (usually out of reach) sampling distribution such that the resulting variance is equal to 0 (see for example Section 6.2 in [START_REF] Lelièvre | Partial differential equations and stochastic methods in molecular dynamics[END_REF] for the connection with the committor function). But on the opposite, as mentioned in the introduction, a bad choice for the sampling distribution may lead to an infinite variance. Hence, the take-home message is that in the best case, Importance Sampling is much better than AMS, but in the worst case, Importance Sampling is much worse than AMS.

To conclude this section, let us say a few words about efficiency. To take into account both computational complexity and variance, Hammersley and Handscomb [START_REF] Hammersley | Monte Carlo methods[END_REF] have proposed to define the efficiency of a Monte-Carlo method as "inversely proportional to the product of the sampling variance and the amount of labour expended in obtaining this estimate." If we consider that the cost of the simulation of a single trajectory is one, than the inverse of the efficiency of a naive Monte-Carlo method to estimate p is just p(1 -p). For AMS with K = 1, it is shown in Corollary 2.9 of [START_REF] Cérou | On the Asymptotic Normality of Adaptive Multilevel Splitting[END_REF] that the number of iterations is -N log p + O P ( √ N ). Taking into account the (quick)sorting of the particles and the fact that an iteration amounts to simulate one new trajectory, the complexity scales like -N log N log p. As explained above, the asymptotic variance satisfies -p 2 log p σ 2 2p(1 -p), hence the following bounds for the inverse of the efficiency of AMS:

(p log p) 2 log N E -1 ams -2p(1 -p) log p log N.
Therefore, in the best case, AMS is much more efficient than naive Monte-Carlo, whereas in the worst case, it is less efficient: larger variance (by a factor 2) and larger complexity (by a factor -log p log N ). However, notice that the latter happens only for a very bad choice of the importance function.

Examples and applications 4.1 Return times

A generic application, proposed in [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF], to rare event probability estimation is to compute return times to some unfrequent event for a stationary ergodic process X with a.s. continuous trajectories. In this work, the return time is defined as the average of the hitting time of a given event with an initial condition X 0 distributed according to the stationary distribution. The considered event is often in the form {Φ(x) L max } for some function Φ.

Let us denote r(L max ) the expectation of the corresponding return time and let us assume that this return time is much larger than the mixing time of X, which happens for L max large enough. In that asymptotics, the authors argue that the sequence of the hitting times of the considered event is distributed according to a Poisson process with parameter 1/r(L max ), so that the return time is explicitly related to the probability of reaching the considered event on a unit time interval. Thus a standard approach to estimate is r(L max ) is to use a block estimator: one simulates X for a long time M * T . On each time block [(m -1)T, mT ], let s m (L max ) = 1 if the level L max is reached by X, and zero otherwise. One can then use the following block estimator, valid when T r(L max ):

r1 B (L max ) = M T M m=1 s m (L max )
.

A better estimator, based on the exponential distribution, and this time valid when T and r(L max ) are comparable, is the modified block estimator r2

B (L max ) = T log(1 -1 M M m=1 s m (L max ))
.

Using AMS Algorithm 3, one can adapt the above estimator while using AMS to estimate the probability of reaching level L max on [0, T ]. In AMS, we have a set of M trajectories that evolve on [0, T ], and we denote by pams the associated estimator of the probability of reaching L max defined in Algorithm 3. One can then adapt the above estimator as follows r3 B (L max ) = T log(1 -pams ) .

Note that, although the initial distribution in the latter AMS variant is supposed to be the stationary distribution of the process, this assumption could be relaxed if we still consider the case where the mixing time of the process is very small as compared to T , so that the influence of a non stationary initial condition is very limited.

Finally, because AMS is used with a deterministic final time, it is worth noted that the best importance function should depend on time (see [START_REF] Lestang | Computing return times or return periods with rare event algorithms[END_REF] section IIIC, see also [START_REF] Bréhier | On a new class of score functions to estimate tail probabilities of some stochastic processes with adaptive multilevel splitting[END_REF]). Also note that using all these trajectories with intermediate levels , it is easy to generate an approximation of the curve → r( ).

Neutronics

The AMS Algorithm 3 was successfully used for problems of neutral particle transport in [START_REF] Louvin | Adaptive Multilevel Splitting for Monte Carlo particle transport[END_REF][START_REF] Louvin | Development of an adaptive variance reduction technique for Monte Carlo particle transport[END_REF]. The adaptation of AMS is quite straightforward, with care to be taken on how to efficiently store the particles trajectories. As we have seen above, the choice of the importance function/reaction coordinate is crucial. Here the authors were able to use some code already developed for other variance reduction techniques as a good importance function, leading to a high figure of merit. The authors also used the fact that AMS does not only give an approximation of a rare event probability, but also a set of empirical trajectories reaching it, which can in turn be used to estimate any observable, given the rare event.

Air Traffic Management

The Multilevel Splitting approach to rare events simulation is now quite widely used in Air Traffic Management (ATM). It is an important application because it has a real life impact, it is not just a simulation to illustrate a theoretical result. Problems faced by people working in ATM include checking that a proposed new regulation in air traffic will not increase the risk of accident (e.g., two planes crashing on each other) above the required safety level. Typically, these safety levels are very low, and therefore one clearly needs to use some variance reduction techniques.

The splitting approach is well suited to this application because it allows one to use a large simulation code already developed separately to simulate air traffic scenarios. The interesting point here is that the splitting paradigm can be used in conjunction with a complex simulation code, that was not initially designed for this purpose. A detailed account of this application can be found in [START_REF] Blom | Rare event estimation for a large-scale stochastic hybrid system with air traffic application[END_REF] and references therein. Note that due to the presence of a discrete component in the state space, the SMC Algorithm 2 had to be adapted to be really efficient to make sure that all the discrete modes are represented.

Molecular dynamics

Recently, the AMS Algorithm 3 (case K = 1) was used in a real Molecular Dynamics problem in [START_REF] Teo | Adaptive multilevel splitting method for molecular dynamics calculation of benzamidinetrypsin dissociation time[END_REF]. The problem is to estimate the expected dissociation time between a protein and a ligand. An abstract view of the problem can be given as follow. The state of the system (configuration of the molecules) can be modeled by a stochastic differential equation in R d (d can be large):

dX t = -∇V (X t )dt + √ 2εdW t ,
where V is a (sufficiently smooth) potential, ε a temperature parameter, and W a standard Brownian motion.

We consider two sets of configurations A and B. A is the set where the protein and its ligand are tightly associated, and B the set where they are completely dissociated. Starting in A, we want to estimate the expected time until the diffusion hits the set B. A is actually a neighborhood of a local minimum of V , and it is assumed to be recurrent. A trajectory leaving a neighborhood A δ of A and hitting B before going back to A is called a reactive trajectory. Algorithm 3 is designed to sample these reactive trajectories efficiently. It represents a rare event because the timescale at which the diffusion can be simulated is orders of magnitudes lower than the timescale of the transitions from A to B.

The solution is not completely straightforward here because AMS simulates reactive trajectories, and estimates their duration, but not the time spent leaving A δ and going back to A again without reaching B. We are also not in the situation considered in Section 4.1 because the duration of a reactive trajectory is much smaller than the time needed to see one happen. The approach in [START_REF] Teo | Adaptive multilevel splitting method for molecular dynamics calculation of benzamidinetrypsin dissociation time[END_REF], also found in [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF], is to write the expectation of the transition time T from A to B as

E[T ] = 1 -p p E[T 1 ] + E[T 2 ],
where p is the (very low) probability of a reactive trajectory, T 1 is the time taken starting from A (at equilibrium) to leave A δ , and go back directly to A, and T 2 the time taken starting from A to leave A δ and go to B without going back to A (reactive trajectory).

Note that this decomposition is exact if the boundaries of A and A δ are level sets of the committor q (see [START_REF] Cérou | A multiple replica approach to simulate reactive trajectories[END_REF] section 4.2). We can nevertheless reasonably conjecture that in most practical cases, it will provide a sharp approximation. We can use Algorithm 3 to estimate p and E[T 2 ], and a direct numerical simulation to estimate E[T 1 ] (no rare events here), and thus get an estimate of E[T ].

Finally, let us mention another recent application of AMS to isomerization when using the NAMD simulation program [START_REF] Lopes | Adaptive multilevel splitting method: Isomerization of the alanine dipeptide[END_REF].

Acknowledgements

This work was partially supported by the French Agence Nationale de la Recherche, under grant ANR-14-CE23-0012, and by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492.


Rafael I Ponce-Toledo¹

David Moreira¹
email: david.moreira@u-psud.fr.

Purificación López-García¹

Philippe Deschamps¹
email: philippe.deschamps@u-psud.fr

Secondary Plastids of Euglenids and Chlorarachniophytes Function with a Mix of Genes of Red and Green Algal Ancestry

Keywords: Chlorarachniophyta, Euglenida, endosymbiotic gene transfer, phylogenomics, plastids

Endosymbiosis has been common all along eukaryotic evolution, providing opportunities for genomic and organellar innovation. Plastids are a prominent example. After the primary endosymbiosis of the cyanobacterial plastid ancestor, photosynthesis spread in many eukaryotic lineages via secondary endosymbioses involving red or green algal endosymbionts and diverse heterotrophic hosts. However, the number of secondary endosymbioses and how they occurred remain poorly understood. In particular, contrasting patterns of endosymbiotic gene transfer (EGT) have been detected and subjected to various interpretations. In this context, accurate detection of EGTs is essential to avoid wrong evolutionary conclusions. We have assembled a strictly selected set of markers that provides robust phylogenomic evidence suggesting that nuclear genes involved in the function and maintenance of green secondary plastids in chlorarachniophytes and euglenids have unexpected mixed red and green algal origins. This mixed ancestry contrasts with the clear red algal origin of most nuclear genes carrying similar functions in secondary algae with red plastids.

Photosynthesis in eukaryotes takes place in a specialized compartment: the plastid. This organelle first evolved in a common ancestor of Archaeplastida (i.e., Viridiplantae + Rhodophyta + Glaucophyta) through the endosymbiosis of a cyanobacterium inside a eukaryotic host [START_REF] Moreira | Sure facts and open questions about the origin and evolution of photosynthetic plastids[END_REF][START_REF] Archibald | The puzzle of plastid evolution[END_REF][START_REF] Keeling | The number, speed, and impact of plastid endosymbioses in eukaryotic evolution[END_REF]. This primary endosymbiotic event entailed massive endosymbiotic gene transfer (EGT) from the cyanobacterial genome to the host nucleus [START_REF] Weeden | Genetic and biochemical implications of the endosymbiotic origin of the chloroplast[END_REF]; [START_REF] Kleine | DNA transfer from organelles to the nucleus: the idiosyncratic genetics of endosymbiosis[END_REF]. Consequently, most proteins required for the proper functioning of primary plastids are encoded in the nuclear genome and addressed to the plastid lumen via specialized signal sequences and a translocation apparatus [START_REF] Gutensohn | structure and function of protein transport machineries in chloroplasts[END_REF]. Other photosynthetic eukaryotic phyla obtained their plastids through secondary endosymbiosis, i.e. the symbiosis of either green or red algae within another eukaryotic cell, or even through tertiary endosymbiosis (symbiosis of secondary photosynthetic eukaryotes within eukaryotic hosts) [START_REF] Delwiche | Tracing the thread of plastid diversity through the tapestry of life[END_REF][START_REF] Archibald | The puzzle of plastid evolution[END_REF][START_REF] Keeling | The number, speed, and impact of plastid endosymbioses in eukaryotic evolution[END_REF]. Euglenida (Excavata) and Chlorarachniophyta (Rhizaria) carry green algal secondary plastids ('green plastids') acquired through two independent endosymbioses involving Prasinophyceae and Ulvophyceae green algae, respectively [START_REF] Rogers | The complete chloroplast genome of the chlorarachniophyte Bigelowiella natans: evidence for independent origins of chlorarachniophyte and euglenid secondary endosymbionts[END_REF][START_REF] Hrdá | The plastid genome of Eutreptiella provides a window into the process of secondary endosymbiosis of plastid in euglenids[END_REF][START_REF] Suzuki | Plastid genome sequences of Gymnochlora stellata, Lotharella vacuolata, and Partenskyella glossopodia reveal remarkable structural conservation among chlorarachniophyte species[END_REF]. Photosynthetic species in the Cryptophyta, Alveolata, Stramenopiles and Haptophyta (CASH) lineages have plastids derived from red algae ('red plastids') but so far it has been impossible to retrace a consensual evolutionary history [START_REF] Lane | The eukaryotic tree of life: endosymbiosis takes its TOL[END_REF][START_REF] Archibald | The puzzle of plastid evolution[END_REF][START_REF] Keeling | The number, speed, and impact of plastid endosymbioses in eukaryotic evolution[END_REF]. Whereas phylogenomic analyses of plastid-encoded genes support the monophyly of all CASH plastids, arguing for a single red algal secondary endosymbiosis [START_REF] Yoon | The single, ancient origin of chromist plastids[END_REF][START_REF] Muñoz-Gómez | The new red algal subphylum Proteorhodophytina comprises the largest and most divergent plastid genomes known[END_REF], most of the phylogenies based on host nuclear genes do not retrieve their monophyly [START_REF] Baurain | Phylogenomic evidence for separate acquisition of plastids in cryptophytes, haptophytes, and stramenopiles[END_REF][START_REF] Burki | Untangling the early diversification of eukaryotes: a phylogenomic study of the evolutionary origins of Centrohelida, Haptophyta and Cryptista[END_REF]). To reconcile these incongruent results, some authors have proposed the hypothesis that a unique phylum (which may have gone extinct or evolved into one of the extant CASH phyla) acquired a red alga through secondary endosymbiosis and originated the first lineage of red secondary algae.

Subsequently, this lineage would have transmitted the secondary red plastid to other CASH phyla via serial tertiary endosymbioses involving different hosts [START_REF] Larkum | Shopping for plastids[END_REF][START_REF] Sanchez-Puerta | A hypothesis for plastid evolution in chromalveolates[END_REF][START_REF] Bodył | Chromalveolate plastids: direct descent or multiple endosymbioses?[END_REF][START_REF] Baurain | Phylogenomic evidence for separate acquisition of plastids in cryptophytes, haptophytes, and stramenopiles[END_REF][START_REF] Petersen | Chromera velia, endosymbioses and the rhodoplex hypothesis--plastid evolution in cryptophytes, alveolates, stramenopiles, and haptophytes (CASH lineages)[END_REF].

As for the primary endosymbiosis, each secondary or tertiary endosymbiosis was accompanied by numerous EGTs from the nucleus of the endosymbiotic red or green alga to the host nucleus.

Consequently, secondary photosynthetic eukaryotes possess two types of genes that can inform about the phylogenetic identity of their plastids: plastid-encoded genes and nucleus-encoded genes acquired via EGT. Genes encoded in primary plastid genomes and the EGTs found in the genomes of Archaeplastida are related to cyanobacteria and have helped to identify the cyanobacterial lineage at the origin of the first plastid [START_REF] Ponce-Toledo | An earlybranching freshwater cyanobacterium at the origin of plastids[END_REF]. Similarly, plastid-encoded genes and EGTs found in nuclear genomes of secondary photosynthetic eukaryotes are expected to be useful to determine the red or green algal origin of their plastids. Compared to plastid-encoded genes, EGTs have the additional advantage that they can inform about the presence and identity of past plastids in lineages where plastids have been lost or replaced (cryptic plastid endosymbioses).

However, if EGTs are valuable to track contemporary and cryptic endosymbioses, their detection within whole nuclear genome sequences remains a complex task [START_REF] Stiller | Experimental design and statistical rigor in phylogenomics of horizontal and endosymbiotic gene transfer[END_REF]). In the case of primary endosymbiosis, EGT detection is rather straightforward because cyanobacterial-type genes are easily distinguishable from typical eukaryotic nuclear genes. The situation is more difficult in the case of secondary endosymbioses. Indeed, detection of EGT genes transferred from the nucleus of green or red algal endosymbionts can be ambiguous due to the poor resolution often found in single gene phylogenies that hampers distinguishing EGTs from vertically inherited nuclear genes, especially considering the short phylogenetic distance between Archaeplastida and several groups of secondary algae. Two studies on red-plastid-bearing algae, the chromerids (Alveolata) and the diatoms (Stramenopiles), illustrate this issue. Both reported an unexpected high number of genes phylogenetically related to green algal homologs. Whereas in the case of the chromerids the green signal was attributed to probable phylogenetic artifacts and the reduced sampling of red algal genome sequences [START_REF] Woehle | Red and problematic green phylogenetic signals among thousands of nuclear genes from the photosynthetic and apicomplexa-related Chromera velia[END_REF], it was interpreted in diatoms as evidence for a cryptic green algal endosymbiont [START_REF] Moustafa | Genomic footprints of a cryptic plastid endosymbiosis in diatoms[END_REF]). However, the subsequent reanalyses of the same genes using richer taxonomic sampling and more robust phylogenetic methods largely erased the evidence for cryptic green endosymbioses in these CASH phyla [START_REF] Burki | Re-evaluating the green versus red signal in eukaryotes with secondary plastid of red algal origin[END_REF][START_REF] Deschamps | Reevaluating the green contribution to diatom genomes[END_REF][START_REF] Moreira | What was the real contribution of endosymbionts to the eukaryotic nucleus? Insights from photosynthetic eukaryotes[END_REF].

The extent and impact of horizontal gene transfer (HGT) on eukaryotic evolution remain controversial topics [START_REF] Leger | Demystifying eukaryote lateral gene transfer[END_REF]. HGTs might be valuable to infer the history of genomes and lineages [START_REF] Abby | Lateral gene transfer as a support for the tree of life[END_REF] but they can also introduce inconvenient noise in phylogenomic analyses, in particular for the study of EGTs [START_REF] Stiller | Experimental design and statistical rigor in phylogenomics of horizontal and endosymbiotic gene transfer[END_REF]. Through time, secondary photosynthetic eukaryotes may have accumulated HGTs in their nuclear genomes from various sources, perhaps even including non-endosymbiotic red or green algae. Unfortunately, gene phylogenies of such HGTs may display topologies comparable to those of EGTs, making them difficult to set apart. In this context, anomalous phylogenetic signal in certain secondary photosynthetic groups has been interpreted as HGT rather than EGT from cryptic endosymbionts. This is the case of the nuclear genome sequence of the green-plastid-containing chlorarachniophyte alga Bigellowiella natans, in which 22% of the genes potentially acquired via HGT appeared to have a red algal origin [START_REF] Curtis | Algal genomes reveal evolutionary mosaicism and the fate of nucleomorphs[END_REF]. Because of the phagotrophic ability of chlorarachniophytes, the presence of these genes was considered to be the result of progressive accumulation of HGTs from red algae or from redplastid-containing CASH lineages, some eventually substituting original 'green' EGTs [START_REF] Archibald | Lateral gene transfer and the evolution of plastid-targeted proteins in the secondary plastid-containing alga Bigelowiella natans[END_REF][START_REF] Yang | An extended phylogenetic analysis reveals ancient origin of "non-green" phosphoribulokinase genes from two lineages of "green" secondary photosynthetic eukaryotes: Euglenophyta and Chlorarachniophyta[END_REF][START_REF] Yang | Phylogenomic analysis of "red" genes from two divergent species of the "green" secondary phototrophs, the chlorarachniophytes, suggests multiple horizontal gene transfers from the red lineage before the divergence of extant chlorarachniophytes[END_REF]. Analogous studies on euglenid species suggested a similar trend for several genes involved in central metabolic pathways [START_REF] Maruyama | Eukaryote-to-eukaryote gene transfer gives rise to genome mosaicism in euglenids[END_REF][START_REF] Yang | An extended phylogenetic analysis reveals ancient origin of "non-green" phosphoribulokinase genes from two lineages of "green" secondary photosynthetic eukaryotes: Euglenophyta and Chlorarachniophyta[END_REF][START_REF] Markunas | Evolutionary history of the enzymes involved in the Calvin-Benson cycle in euglenids[END_REF]. The unexpected presence of those 'red' genes in chlorarachniophytes and euglenids was first considered as the result of multiple HGTs (e.g., [START_REF] Archibald | Lateral gene transfer and the evolution of plastid-targeted proteins in the secondary plastid-containing alga Bigelowiella natans[END_REF][START_REF] Maruyama | Eukaryote-to-eukaryote gene transfer gives rise to genome mosaicism in euglenids[END_REF] but the increasing number of reported cases has prompted some authors to speculate on putative cryptic red algal endosymbioses in both lineages [START_REF] Maruyama | Eukaryote-to-eukaryote gene transfer gives rise to genome mosaicism in euglenids[END_REF][START_REF] Markunas | Evolutionary history of the enzymes involved in the Calvin-Benson cycle in euglenids[END_REF]. A systematic investigation of HGT/EGT is still missing in euglenids and chlorarachniophytes but, as mentioned above, in the context of secondary endosymbioses it can be difficult to distinguish among HGT, EGT, and just unresolved trees on the basis of single-gene phylogenies [START_REF] Deschamps | Reevaluating the green contribution to diatom genomes[END_REF].

In this work, we have focused on a particular group of genes to reduce this uncertainty: genes transferred from the original cyanobacterial plastid endosymbiont into the nuclear genome of Archaeplastida and subsequently transferred from Archaeplastida into the genomes of complex secondary algae. In Archaeplastida, these genes are known to be involved in essential plastid functions and tend to be highly conserved [START_REF] Reyes-Prieto | Cyanobacterial contribution to algal nuclear genomes is primarily limited to plastid functions[END_REF][START_REF] Deschamps | Signal conflicts in the phylogeny of the primary photosynthetic eukaryotes[END_REF], so we expected that they can provide strong phylogenetic signal. To identify them, we queried by BLAST the whole predicted proteomes of Guillardia theta and Bigelowiella natans against a local genome database containing representatives of the three domains of life, in particular a comprehensive collection of genomes and transcriptomes of photosynthetic protists (supplementary table S1, Supplementary Material online). Guillardia and Bigelowiella proteins with hits in other photosynthetic eukaryotes and in cyanobacteria were selected for phylogenetic analysis. Maximum likelihood (ML) phylogenetic trees for these proteins were constructed and manually filtered to retain those fulfilling two criteria: i) trees have to support a clear separation of Viridiplantae and Rhodophyta (with secondary lineages branching within them), and ii) proteins have to be shared by at least three secondary photosynthetic lineages. We identified in this way 82 genes most likely acquired by secondary photosynthetic eukaryotes from Archaeplastida. 70 were cyanobacterial genes likely transferred sequentially through primary and secondary endosymbioses, and 12 were derived from diverse bacterial groups likely transferred to a common ancestor of Archaeplastida and subsequently transferred to secondary photosynthetic groups (supplementary table S3 and figs. S1-S82, Supplementary Material online). Interestingly, most of these genes were absent in nonphotosynthetic eukaryotes, supporting that they were not misinterpreted vertically-inherited ones.

Most of the 82 ML phylogenies were well resolved and enabled us to unambiguously determine, for each secondary lineage, whether the source of the gene was a green or a red alga. As expected, in the great majority of our trees (between 84 and 90%, fig. 1A) the genes of red-plastid-endowed CASH lineages derived from red algae (e.g., fig. 2A and2B). Because of their secondary green plastids, we expected the opposite situation in chlorarachniophytes and euglenids, namely a majority of 'green' genes. However, 42 of the 78 trees where chlorarachniophytes were present (54%, fig. 1A) supported a 'red' origin of the corresponding genes (e.g., fig. 2A). Similarly, 22 of the 61 trees containing euglenids (36%, fig. 1A) also supported a 'red' ancestry (e.g., fig. 2B). These surprisingly high values were in sharp contrast with the small number of trees (<10%, fig. 1A)

showing CASH phyla embedded within green algae. Interestingly, the CASH phyla were monophyletic in 7 of these trees, arguing for a common evolutionary origin of the corresponding 'green' genes. Almost all of the 82 genes identified here encode plastid-targeted proteins involved in essential plastid functions (fig. 1B and supplementary table S4, Supplementary Material online). For instance, in both chlorarachniophytes and euglenids, these nuclear-encoded 'red' genes participate in plastid genome expression (e.g., elongation factors and aminoacyl-tRNA synthetases), light harvesting, chlorophyll biosynthesis, and photosystem II assembly. Keeping these important genes implies a plastid-related selective pressure, which excludes that they could have accumulated in the heterotrophic ancestors of green secondary photosynthetic eukaryotes prior to plastid acquisition.

The marked disproportion of unexpected gene sources in green versus red secondary photosynthetic lineages is intriguing and may be interpreted in different ways. First, the green algal ancestors of chlorarachniophyte and euglenid plastids may have had a high proportion of red algal HGT genes in their genomes. However, such a high HGT proportion involving essential genes has not been reported so far in any green alga. Second, these 'red' genes may have accumulated in chlorarachniophyte and euglenid nuclear genomes by numerous HGTs, for example from food sources. This would imply that, for unknown reasons, HGT is much more frequent in secondary green lineages than in red ones, as well as a long-lasting feeding preference towards red prey in both secondary green lineages. Moreover, the 'red' genes are shared by all the species of the relatively rich taxon sampling available for chlorarachniophytes (fig. 2A), indicating that their acquisition predated the diversification of this group and stopped afterwards (we did not retrieve any tree supporting a recent HGT involving only a subgroup of chlorarachniophytes). Our data therefore argue for an ancient timing of 'red' gene acquisition. These observations may support a third interpretation: the 'red' genes are shared by all SAR lineages (Stramenopiles, Alveolata, and Rhizaria) because they were acquired from a single common secondary red algal endosymbiosis ancestral to the whole SAR supergroup. This original red plastid would have been lost in many phyla and replaced by a green alga in chlorarachniophytes. However, this scenario poses several problems. On the one hand, traces of past presence of red algal plastids, in the form of EGTs, in non-photosynthetic SAR lineages are very often controversial [START_REF] Elias | Sizing up the genomic footprint of endosymbiosis[END_REF]Stiller et al. 2009;[START_REF] Stiller | Experimental design and statistical rigor in phylogenomics of horizontal and endosymbiotic gene transfer[END_REF]). On the other hand, plastid-bearing chlorarachniophytes constitute a relatively late-emerging branch within SAR [START_REF] Sierra | Evolutionary origins of rhizarian parasites[END_REF], implying that if their present-day green plastid replaced a former red one, this red plastid would have had to be present until recently and been lost in all other rhizarian lineages, which may seem unparsimonious. The case of euglenids is even more difficult to interpret as this group of excavates has no close phylogenetic relationship with any other photosynthetic lineage. In addition, massive sequence data remain much more limited for euglenids than for chlorarachniophytes (only a few transcriptomes available, see supplementary table S1, Supplementary Material online), making it difficult to infer the relative age of possible gene transfers. Nonetheless, 'red' genes were often shared by several euglenids in our trees, suggesting a similar pattern of ancient acquisition as in chlorarachniophytes (supplementary figs. S1-S82, Supplementary Material online).

Our results show the presence of an unexpectedly high number of genes of red algal affinity in the two groups of eukaryotic algae with secondary green plastids, the euglenids and chlorarachniophytes, which is significantly higher than the frequency of 'green' genes in algae with secondary red plastids, the CASH lineages. To address this question, we have focused on a subset of genes selected because of their strong phylogenetic signal and their implication in plastid-related activities. It is therefore uncertain whether this conclusion can be applied to the rest of HGTs/EGTs potentially present in the genomes of all these algae. In fact, in addition to the problems inherent to the accurate detection of EGTs, our focus on these specific genes may explain, at least partly, the different results obtained in recent analyses of all potential EGTs in some CASH lineages, not only those of ultimate cyanobacterial origin (e.g., [START_REF] Dorrell | Chimeric origins of ochrophytes and haptophytes revealed through an ancient plastid proteome[END_REF]). However, we could not identify any particular bias in our gene selection process that could have artificially enriched the observed 'red' gene frequency in euglenids and chlorarachniophytes.

Despite the methodological problems inherent to global genome analyses cited above, including a highly unbalanced representation of red and green algal genomes in sequence databases [START_REF] Deschamps | Reevaluating the green contribution to diatom genomes[END_REF], the study of the chlorarachniophyte B. natans genome already pointed in that direction, with 22% of EGT genes of apparent red algal ancestry [START_REF] Curtis | Algal genomes reveal evolutionary mosaicism and the fate of nucleomorphs[END_REF].

The origin of the 'red' genes in euglenids and chlorarachniophytes, either by cumulative HGT or by EGT from cryptic red algal endosymbionts, remains mysterious but our work indicates that they were acquired early in both groups and that they fulfill essential functions for plastid activity and maintenance. Interestingly, indisputable evidence supports that in a third group of complex algae with green plastids, the dinoflagellate genus Lepidodinium, a former red plastid was replaced by the current green one, leading to a mosaic plastid proteome encoded by a mix of red and green algal genes [START_REF] Minge | A phylogenetic mosaic plastid proteome and unusual plastid-targeting signals in the green-colored dinoflagellate Lepidodinium chlorophorum[END_REF], reminiscent of those found in euglenids and chlorarachniophytes. It has been proposed that, since they retain more gene-rich genomes than green ones, red plastids have increased capacity for autonomous metabolism that could explain why they are more widespread across the diversity of eukaryotes as secondary plastids (the "portable plastid" hypothesis [START_REF] Grzebyk | The mesozoic radiation of eukaryotic algae: the portable plastid hypothesis[END_REF]). It is thus tempting to speculate for euglenids and chlorarachniophytes a similar case as for Lepidodinium, with initial red plastids subsequently replaced by green ones. Even if this hypothesis turns out to be wrong and these cryptic red endosymbioses did not exist, the ancient acquisition by another mechanism of a significant number of red algal genes in both groups before their diversification and, especially, their maintenance in the contemporary species through millions of years of evolution, suggest that the 'red' genes were instrumental in the establishment and maintenance of the secondary green plastids. Sequencing and analysis of additional genomes of euglenids, chlorarachniophytes, and their non-photosynthetic relatives will help to refine the inventory of 'red' genes in these lineages and their timing and, eventually, mechanism of acquisition.

Materials and Methods

Sequence Analysis

A local database was constructed to host the predicted proteomes from various nuclear genomes and transcriptomes as well as plastid genomes (for the complete list, see supplementary table S1, Supplementary Material online). All proteins of the Bigelowiella natans (Chlorarachniophyta) and Guillardia theta (Cryptophyta) predicted proteomes were used as queries for BLASTp sequence similarity searches [START_REF] Camacho | BLAST+: architecture and applications[END_REF]) against the local database. We retained up to 350 top hits with an e-value threshold of 1e-05. BLASTp outputs were parsed with a custom Python script to identify the proteins having hits in diverse photosynthetic eukaryotes and that were more similar to cyanobacteria or other bacteria than to non-photosynthetic eukaryotes.

For these proteins, reciprocal BLASTp searches were done against the database to collect up to 600 similar sequences. We then aligned each set using Mafft v7.123b [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] with default parameters. Non-conserved alignment regions were trimmed with BMGE v1.0 [START_REF] Criscuolo | BMGE (Block Mapping and Gathering with Entropy): a new software for selection of phylogenetic informative regions from multiple sequence alignments[END_REF] with the BLOSUM62 matrix and allowing less than 50% gaps per position. Preliminary phylogenetic trees were inferred from trimmed alignments using FastTree v2.1.7 [START_REF] Price | FastTree 2--approximately maximum-likelihood trees for large alignments[END_REF]) with default parameters. These trees were then manually inspected to identify those compatible with an EGT/HGT scenario. For all positive cases, only the sequences corresponding to the portion of interest of each phylogenetic tree (the part showing the photosynthetic eukaryotes and the closest outgroup) were retained for the remaining steps. We then removed very short partial sequences and, to speed up subsequent calculations, several outgroup sequences from all alignments (see supplementary table S2, Supplementary Material online). The final sequence datasets were realigned and trimmed using TrimAL v1.4.rev15 with "gappy-out" parameter [START_REF] Capella-Gutierrez | trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses[END_REF]. ML phylogenetic trees were inferred using IQtree v1.5.1 with the PMSF model of sequence evolution [START_REF] Wang | Modeling site heterogeneity with posterior mean site frequency profiles accelerates accurate phylogenomic estimation[END_REF]) parameterized using guided trees constructed with the LG+G+I model. Statistical support was calculated with 1000 ultrafast bootstrap replicates [START_REF] Minh | Ultrafast approximation for phylogenetic bootstrap[END_REF][START_REF] Nguyen | IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies[END_REF][START_REF] Hoang | UFBoot2: Improving the ultrafast bootstrap approximation[END_REF].

Final selection of trees was done by manual inspection to keep those fulfilling the following two requirements: i) the protein had to be shared by Cyanobacteria (or other bacteria), Archaeplastida and at least three secondary photosynthetic lineages, and ii) the corresponding phylogenetic trees had to support the clear separation of Viridiplantae and Rhodophyta (plus the lineages with secondary green and red plastids nested within them). Finally, the 82 trees passing this final filter (supplementary figs. S1-S82, Supplementary Material online) were inspected to infer the phylogenetic origin of the corresponding genes in the secondary photosynthetic lineages (supplementary table S3, Supplementary Material online).

Gene Functional Annotation

We annotated the functions of the 82 proteins from the final selection (see above) through the EggNOG 4.5 [START_REF] Huerta-Cepas | eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences[END_REF] web portal (http://eggnogdb.embl.de). For each protein we used as queries the ortholog sequences of Guillardia theta and Bigelowiella natans. Functional annotations are shown in supplementary table S4, Supplementary Material online.

Figure captions FIG 1 .

 captions1 Figure captions

Acknowledgments

This study was supported by European Research Council grant ProtistWorld (P.L.-G., agreement no. 322669), the Université Paris-Sud program "Attractivité" (P.D.) and the Agence Nationale de la Recherche (D.M., project ANR-15-CE32-0003 "ANCESSTRAM"). We thank the Associated Editor and two anonymous reviewers for constructive comments.

Data Availability

Protein sequence datasets used in this work are available for download at http://www.ese.upsud.fr/article950.html?lang=en. They include nonaligned sequences and trimmed alignments.

Supplementary Material

Supplementary figures S1-S82 and tables S1-S4 are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).