Yvon Maday

Olga Mula

An Adaptive Parareal Algorithm

come

Introduction

Solving complex models with high accuracy and within a reasonable computing time has motivated the search for numerical schemes that exploit efficiently parallel computing architectures. In this paper, the model consists of a Partial Differential Equation (PDE) set on a domain D. In this context, one of the main ideas to parallelize a simulation is to break the problem into subproblems defined over subdomains of a partition of D. The domain can potentially have high dimensionality and be composed of different variables like space, time, velocity or even more specific variables for some problems. While there exist algorithms with very good scalability properties for the decomposition of the spatial variable in elliptic and saddle-point problems (see [START_REF] Quarteroni | Domain decomposition methods for partial differential equations[END_REF] or [START_REF] Toselli | Domain decomposition methods: algorithms and theory[END_REF] for an overview), the same cannot be said for the decomposition of time of even simple systems of ODEs. This is despite the fact that research on time domain decomposition is currently very active and has by now a history of at least 50 years (back to at least [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF]) during which several algorithms have been explored (see [START_REF] Gander | 50 years of time parallel time integration[END_REF] for an overview). As a consequence, time domain decomposition is to date only a secondary option when it comes to deciding what algorithm/method distributes the tasks in a parallel cluster.

The main goal of this work is to address this efficiency limitation in the framework of one particular scheme: the parareal in time algorithm. The method was first introduced in [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF] and has been well accepted by the community because it is easily applicable to a relatively large spectrum of problems. (Some specific difficulties are nevertheless encountered on certain types of PDEs as reported in, e.g., [START_REF] Dai | Stable parareal in time method for first-and second-order hyperbolic systems[END_REF][START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF] for hyperbolic systems or [START_REF] Bal | Symplectic parareal[END_REF][START_REF] Dai | Symmetric parareal algorithms for hamiltonian systems[END_REF] for hamiltonian problems). Another ingredient for its success is that, even though its scalability properties are limited, they are in general competitive in comparison with other methods. Without entering into very specific details of the algorithm at this stage, we can summarize the procedure by saying that we build iteratively a sequence to approximate the exact solution of the problem by a predictorcorrector algorithm. At every iteration, predictions are made by a solver which has to be as numerically inexpensive as possible since it is run on the full time interval. It usually involves coarse physics and/or coarse resolution. Corrections involve an expensive solver with highfidelity physics and high resolution which is propagated in parallel over small time subdomains. In the classical version of parareal, the fine solver has a fixed high accuracy across all iterations. It is set to the one that we would use to solve the dynamics at the desired accuracy with a purely sequential solver. It is well-known that this point is the major obstacle to achieve better parallel efficiency. In this paper, we propose an adaptive variant where the accuracy of the fine solver is increased across the iterations. Our main goal is to show that this new point of view overcomes the obstacle of the cost of the fine solver and that the only remaining factor limiting high performance becomes the cost of the coarse solver. We refer to, e.g., [START_REF] Carlberg | Data-driven time parallelism via forecasting[END_REF] for contributions on the lowering of the cost of that coarse solver.

We present in section 2 the new adaptive point of view. This requires to formulate an idealized version of the parareal algorithm in an infinite dimensional function space where the fine propagations are replaced by the exact ones (section 2.2). Since this scheme is obviously not implementable in practice, we formulate a feasible "perturbed" version that involves approximations of the exact propagations at increasing accuracy across the iterations (section 2.3). The accuracies are tightened in such a way that the feasible adaptive algorithm converges at the same rate as the ideal one and with a near-minimal numerical cost. The identified tolerances involve quantities that are difficult to estimate in practice. In addition, they may not be optimal because they are derived from a theoretical convergence analysis based on abstract conditions for the coarse and fine solvers. We bridge this gap between theory and actual implementation by proposing practical guidelines to set these tolerances. We next explain in section 2.4 how the new formulation invites to use adaptive schemes not only in the time variable, but also in other variables that may be involved in the dynamics. The performance of the algorithm could also be enhanced by re-using informations from previous iterations in order to limit the cost of internal solvers. The techniques for this will strongly depend on the nature of the specific problem. We discuss common situations in Appendix A. We close section 2.5 by listing the main advantages of the new framework and how the classical parareal paradigm can be formulated with the optics of the new standpoint.

The parallel performance of the adaptive scheme is difficult to predict a priori but in section 3 we carry a discussion where we show that it will always be superior to the classical approach. In the idealized situation where the cost of the coarse solver and communication delays are negligible, we show that the algorithm would exhibit a very high parallel efficiency.

Finally, in section 4, we illustrate the performance of the algorithm in a numerical example: the Brusselator system. It is a stiff ODE system which is challenging because it does not allow to use a very inexpensive coarse solver. We show that the adaptive parareal algorithm systematically performs around 2.5 times better than the classical one. In addition, we confirm in the example that the only remaining obstacle to achieve a very competitive performance is the cost of the coarse solver. Due to the nature of the algorithm, it is not clear how to overcome this limitation and addressing this issue is deferred to future works.

We conclude this introduction by some bibliographical remarks. To the best of our knowledge, the current abstract and broad formulation of an adaptive version of parareal has never been proposed in the literature. However, previous works have instantiated in a variety of particular applications the idea of re-using information from previous parareal iterations and this is the main point of contact between this work and previous contributions. Among the most relevant ones stand the coupling of the parareal algorithm with spatial domain decomposition (see [START_REF] Maday | The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation[END_REF][START_REF] Guetat | Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine[END_REF][START_REF] Aouadi | Convergence analysis of the coupled parareal-schwarz waveform relaxation method[END_REF]), the combination of the parareal algorithm with iterative high order methods in time like spectral deferred corrections (see [START_REF] Minion | Parareal and spectral deferred corrections[END_REF][START_REF] Minion | A hybrid parareal spectral deferred corrections method[END_REF][START_REF] Minion | Interweaving PFASST and parallel multigrid[END_REF]) and, in a similar spirit, applications of the parareal algorithm to solve optimal control problems (see [START_REF] Maday | The Parareal in Time Iterative Solver: a Further Direction to Parallel Implementation[END_REF][START_REF] Maday | Monotonic parareal control for quantum systems[END_REF]). In appendix A, we briefly explain in what sense the two first applications can be seen as particular instances of the current approach and how our viewpoint could help to give them more solid theoretical foundations. Finally, we identify another relevant scenario where efficiency could be enhanced by reusing information from previous iterations: the solution of time-dependent problems involving internal iterative schemes at every time step. This idea was first proposed in [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] and a more complete analysis will be proposed in a forthcoming work.

An adaptive parareal algorithm

In this section, after introducing some preliminary notations in section 2.1, we formulate an ideal parareal scheme on an infinite dimensional functional setting (section 2.2). We then present feasible realizations involving a fine solver whose accuracy is adaptively increased across the iterations (section 2.3). We prove that the feasible adaptive algorithm converges at the same rate as the ideal one provided that the tolerances of the fine solver are increased at certain rate which will be discussed. Finally, we discuss how the new paradigm can be realized thanks to adaptive schemes and/or the re-use of information from previous steps (section 2.4 and appendix A). In section 2.5, we connect the new adaptive formulation with the classical parareal algorithm and list the main advantages of the new standpoint.

Setting and preliminary notations

Let U be a Banach space of functions defined over a domain

Ω ⊂ R d (d ≥ 1), e.g. U = L 2 (Ω). Let E : [0, T ] × [0, T ] × U → U
be a propagator, that is, an operator such that, for any given time t ∈ [0, T ], s ∈ [0, T -t] and any function w ∈ U, E(t, s, w) takes w as an initial value at time t and propagates it at time t + s. We assume that E satisfies the semi group property E(r, t -r, w) = E(s, t -s, E(r, s -r, w)), ∀w ∈ U, ∀(r, s, t) ∈ [0, T ] 3 , r < s < t.

We further assume that E is implicitly defined through the solution u ∈ C 1 ([0, T ], U) of the time-dependent problem

u (t) + A (t, u(t)) = 0, t ∈ [0, T ], (1)
where A is an operator from [0, T ] × U into U with adequate regularity we shall detail latter. Then, given any w ∈ U, E(t, s, w) denotes the solution to (1) at time t + s with initial condition w at time t ≥ 0. In our problem of interest, we study the evolution given by (1) when the initial condition is u(0) ∈ U. Note that E could also be associated to a discretized version of the evolution equation or be defined through an operator that is not necessary related to an evolution equation (see [START_REF] Gaja | Parallel in time algorithms for nonlinear iterative methods[END_REF]). Since, in general, the problem does not have an explicit solution, we seek to approximate it at a given target accuracy. For any initial value w ∈ U, any t ∈ [0, T [, s ∈ [0, T -t] and any ζ > 0 we denote by [E(t, s, w); ζ] an element of U that approximates E(t, s, w) such that we have

E(t, s, w) -[E(t, s, w); ζ] ≤ ζ s (1 + w ), (2)
where, here and in the following, • denotes the norm in U. Any realization of [E(t, s, w); ζ] involves three main ingredients: i) a numerical scheme to discretize the time dependent problem (1) (e.g. an Euler scheme in time), ii) a certain expected error size associated with the choice of the discretization (e.g. error associated with the time step size of the Euler scheme),

iii) a numerical implementation to solve the resulting discrete systems (e.g. conjugate gradient, Newton method, SSOR, . . . ).

In the following, we will use the term solver to denote a particular choice for i), ii) and iii). Given a solver S, we will use the same notation as for the exact propagator E to express that

S(t,

An idealized version of the parareal algorithm

Let be given a decomposition of the time interval [0, T ] into N subintervals [T N , T N +1 ], N = 0, . . . , N -1. Without loss of generality, we will take them of uniform size ∆T = T /N which means that T N = N ∆T for N = 0, . . . , N . For a given target accuracy η > 0, the primary goal of the parareal in time algorithm is to build an approximation ũ(T N ) of u(T N ) such that max

1≤N ≤N u(T N ) -ũ(T N ) ≤ η. (3)
The classical way to achieve this is to set

ũ(T N ) = S seq (0, T N , u(0)) = [E(0, T N , u(0)); ζ], 1 ≤ N ≤ N ,
where S seq is some sequential solver in [0, T ] with ζ = η/(T (1 + u(0) )) in [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF]. Since this comes at the cost of solving over the whole time interval [0, T ], the main goal of the parareal in time algorithm is to speed up the computing time, while maintaining the same target accuracy η. This is made possible by first decomposing the computations over the time domain. Instead of solving over [0, T ], we perform N parallel solves over each intervals (T N , T N +1 ]of size ∆T . We next introduce an idealized version of it which will not be feasible in practice but will be the starting point of subsequent implementable versions. The algorithm relies on the use of a solver G (known as the coarse solver) with the following properties involving the operator

δG := E -G.

Hypotheses (H):

There exists ε G , C c , C d > 0 such that for any function x, y ∈ U and for any t

∈ [0, T [ and s ∈ [0, T -t], G(t, s, x) = [E(t, s, x), ε G ] ⇔ δG(t, s, x) ≤ s(1 + x )ε G (4a) G(t, s, x) -G(t, s, y) ≤ (1 + C c s) x -y , (4b) δG(t, s, x) -δG(t, s, y) ≤ C d sε G x -y (4c)
Note that these hypothesis are the classical abstract formulations of the properties of numerical schemes related to stability and accuracy. Hypothesis (4b) is a Lipschitz condition and the quantity ε G is a small constant which, in the case of a Euler scheme, would be equal to the time step size.

The idealized version of the algorithm consists in building iteratively a series (y N k ) k of approximations of u(T N ) for 0 ≤ N ≤ N following the recursive formula

             y N +1 0 = G(T N , ∆T, y N 0 ), 0 ≤ N ≤ N -1 y N +1 k+1 = G(T N , ∆T, y N k+1 ) + E(T N , ∆T, y N k ) -G(T N , ∆T, y N k ), 0 ≤ N ≤ N -1, k ≥ 0, y 0 0 = u(0). (5)
At this point, several comments are in order. The first one is that the computation of y N k only requires propagations with E over intervals of size ∆T . As follows from ( 5), for a given iteration k, N propagations of this size are required, each of them over distinct intervals [T N , T N +1 ] of size ∆T , each of them with independent initial conditions. Since they are independent from each other, they can be computed over N parallel processors and the original computation over [0, T ] is decomposed into parallel computations over N subintervals of size ∆T . The second observation is that the algorithm may not be implementable in practice because it involves the exact propagator E. Feasible instantiations consist in replacing

E(T N , ∆T, y N k ) by some approximation [E(T N , ∆T, y N k ), ζ N k ]
with a certain accuracy ζ N k which has to be carefully chosen. We will come to this point in the next section. The third observation is to note that, in the current version of the algorithm, for all N = 0, . . . , N , the exact solution u(T N ) is obtained after exactly k = N parareal iterations. This number can be reduced when we only look for an approximate solution with accuracy η. Depending on the problem, the final number of iterations K(η) can actually be much smaller than N . The convergence result of theorem 2.1 and its proof are helpful to understand the main mechanisms driving the convergence of the algorithm and explaining its behavior. To present it, we introduce the shorthand notation for the error norm

E N k := u(T N ) -y N k , k ≥ 0, 0 ≤ N ≤ N ,
max 0≤N ≤N u(T N ) -y N k ≤ µ τ k+1 (k + 1)! , ∀k ≥ 0. ( 6
)
Proof. The proof is in the spirit of existing results from the literature (see [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF][START_REF] Bal | Parallelization in time of (stochastic) ordinary differential equations[END_REF][START_REF] Maday | The parareal in time algorithm: Basics, stability analysis and more[END_REF][START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF]) but it is instructive to give it for subsequent developments in the paper. We introduce the following quantities

       α := C d ε G ∆T β := 1 + C c ∆T γ := ∆T ε G max 0≤N ≤N (1 + u(T N ) ) (7)
as shorthand notations for the proof.

If k = 0, using definition (5) for y N 0 , we have for 0

≤ N ≤ N -1, E N +1 0 = y N +1 0 -u(T N +1 ) = G(T N , ∆T, y N 0 ) -E(T N , ∆T, u(T N )) ≤ G(T N , ∆T, y N 0 ) -G(T N , ∆T, u(T N )) + G(T N , ∆T, u(T N )) -E(T N , ∆T, u(T N )) ≤ (1 + C c ∆T )E N 0 + ∆T ε G (1 + u(T N ) ) ≤ βE N 0 + γ,
where we have used (4a) and (4b) to derive the second to last inequality. For k ≥ 1, starting from (5), we have

y N +1 k -u(T N +1 ) = G(T N , ∆T, y N k ) + E(T N , ∆T, y N k-1 ) -G(T N , ∆T, y N k-1 ) -E(T N , ∆T, u(T N )) = G(T N , ∆T, y N k ) -G(T N , ∆T, u(T N )) + δG(T N , ∆T, y N k-1 ) -δG(T N , ∆T, u(T N ))
. Taking norms and using (4b), (4c), we derive

E N +1 k ≤ βE N k + αE N k-1
, Following [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF], we consider the sequence (e N k ) N,k≥0 defined recursively as follows. For k = 0,

e N 0 = 0, if N = 0 βe N -1 0 + γ, if N ≥ 1 (8)
and for k ≥ 1,

e N k = 0, if N = 0 αe N -1 k-1 + βe N -1 k , if N ≥ 1. ( 9
)
Since E N k ≤ e N k for k ≥ 0 and N = 0, . . . , N , we analyze the behavior of (e N k ) to derive a bound for E N k . For this, we consider the generating function

ρ k (ξ) = N ≥0 e N k ξ N .
From ( 8) and ( 9) we get

ρ k (ξ) = αξρ k-1 (ξ) + βξρ k (ξ), k ≥ 1 ρ 0 (ξ) = γ ξ 1-ξ + βξρ 0 (ξ), from which we derive ρ k (ξ) = γα k ξ k+1 (1 -ξ) 1 (1 -βξ) k+1 , k ≥ 0.
Since, β ≥ 1, we can bound the term (1 -ξ) in the denominator by (1 -βξ). Next, using the binomial expansion

1 (1 -βξ) k+2 = j≥0 k + 1 + j j β j ξ j ( 10
)
and identifying the term in ξ N in the expansion, we derive the bound

e N k ≤ γα k β N -k-1 N k + 1 .
Hence, using definition [START_REF] Dai | Symmetric parareal algorithms for hamiltonian systems[END_REF] for α, β and γ,

E N k ≤ e N k ≤ (1 + C c ∆T ) N -k-1 max 0≤N ≤N (1 + ||u(T N ) ) C d (k + 1)! [C d ε G e -Cc∆T T N k+1 ,
which ends the proof of the lemma.

Note that at least one step is not sharp in the above proof: it is the step where 1 -ξ is replaced by 1-βξ. Note also that τ is the quantity driving convergence and its speed. It follows that a sufficient condition to converge is to use a coarse solver such that

τ < 1 ⇔ ε G < e Cc∆T C d T .
This gives a certain limitation on the definition of the coarse solver since it imposes some minimal accuracy. In the following, we will work under this assumption for ε G .

Feasible realizations of the parareal algorithm

Feasible versions of algorithm (5) involve approximations of E(T N , ∆T, y N k ) with a certain accuracy ζ N k . This leads to consider algorithms of the form

             y N +1 0 = G(T N , ∆T, y N 0 ), 0 ≤ N ≤ N -1 y N +1 k+1 = G(T N , ∆T, y N k+1 ) + [E(T N , ∆T, y N k ); ζ N k ] -G(T N , ∆T, y N k ), 0 ≤ N ≤ N -1, k ≥ 0, y 0 0 = u(0). ( 11
)
Since no feasible version will converge at a better rate than [START_REF] Carlberg | Data-driven time parallelism via forecasting[END_REF], we analyze here what is the minimal accuracy ζ N k that preserves it. A result in this direction is given in the following theorem. It requires to introduce the quantity

ν p := max 0≤N ≤N (1 + y N p ) max 0≤N ≤N (1 + u(T N ) )
, ∀p ≥ 0.

which tends to 1 as p → ∞.

Theorem 2.2. Let G and δG satisfy Hypothesis (4). Let k ≥ 0 be any given positive integer. If for all 0 ≤ p < k and all 0 ≤ N < N , the approximation

[E(T N , ∆T, ζ N p )] has accuracy ζ N p ≤ ζ p := ε p+2 G (p + 1)! ν p , ( 12
)
then the (y N k ) N of the feasible parareal scheme [START_REF] Emmett | Toward an efficient parallel in time method for partial differential equations[END_REF]

satisfy max 0≤N ≤N u(T N ) -y N k ≤ µ τ k+1 (k + 1)! , ( 13
)
with τ := 1 + C d T e -Cc∆T ε G .
Let us make a couple of remarks before giving the proof of the theorem. First, comparing (6) and ( 13), we see that the feasible parareal algorithm converges at a rate close to the ideal one in the sense that it only deviates by a factor

τ τ = τ + ε G τ = 1 + e Cc∆T C d T .
For a given problem with fixed C c , C d and T , this deviation depends on the size of ∆T , which is itself driven by the number of processors N . The minimal accuracy to update the realization of E(T N , ∆T, y N k ) is given by [START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF]. As a final remark, we note that we can also interpret both convergence bounds ( 6) and ( 13) in terms of ε G . With this viewpoint, the error at iteration k is bounded by C(k)ε k+1 G /(k + 1)! where C(k) is a factor which grows (only) exponentially with k and is thus "swallowed" by the denominator.

Proof. The proof follows the same lines as the one for theorem 2.1 and E N k , α, β, γ are defined exactly as before. In addition, it will be useful to introduce the sequence

g k = ζ k ∆T max 0≤N ≤N (1 + y N k ), ∀k ≥ 0 g -1 = γ
We concentrate on the case k ≥ 1 since the case k = 0 is identical as in theorem 2.1. For k ≥ 1, using [START_REF] Emmett | Toward an efficient parallel in time method for partial differential equations[END_REF], we have

y N +1 k -u(T N +1 ) = G(T N , ∆T, y N k ) -G(T N , ∆T, u(T N )) -G(T N , ∆T, y N k-1 ) + G(T N , ∆T, u(T N )) + [E(T N , ∆T, y N k-1 ); ζ N k-1 ] -E(T N , ∆T, u(T N )) = G(T N , ∆T, y N k ) -G(T N , ∆T, u(T N )) + δG(T N , ∆T, y N k-1 ) -δG(T N , ∆T, u(T N )) + [E(T N , ∆T, y N k-1 ); ζ N k-1 ] -E(T N , ∆T, y N k-1 ).
Taking norms, using (4b), (4c) and the definition (2) applied to [E(T N , ∆T,

y N k-1 ); ζ N k-1 ],
we derive

E N +1 k ≤ [1 + C c ∆T ]E N k + C d ∆T ε G E N k-1 + ζ N k-1 ∆T max 0≤N ≤N (1 + y N k-1 ) ≤ βE N k + αE N k-1 + g k-1 .
Similarly to theorem 2.1, we introduce the sequence (ẽ N k ) N,k≥0 defined for k = 0 as ẽN 0 = e N 0 for all N ≥ 0 and for k ≥ 1,

ẽN k = 0, if N = 0 αẽ N -1 k-1 + βẽ N -1 k + g k-1 , if N ≥ 1
The associated generating function ρk satisfies

   ρk (ξ) = αξ ρk-1 (ξ) + βξ ρk (ξ) + g k-1 ξ 1-ξ , ∀k ≥ 1, ρ0 (ξ) = ρ 0 (ξ) = γξ (1-ξ)(1-βξ) . Hence ρk (ξ) = αξ 1 -βξ ρk-1 (ξ) + ξ (1 -ξ)(1 -βξ) g k-1 = αξ 1 -βξ k ρ0 (ξ) + ξ (1 -ξ)(1 -βξ) k-1 =0 αξ 1 -βξ g k-1-
By replacing again at the denominator the factor (1 -ξ) by (1 -βξ) and using the binomial expansion [START_REF] Dutt | Spectral deferred correction methods for ordinary differential equations[END_REF], we derive the bound

ρk (ξ) ≤ γα k ξ k+1 j≥0 k + 1 + j j β j ξ j + k-1 =0 α ξ +1 g k-1- j≥0 + 1 + j j β j ξ j , = j≥0 k =0 α g k-1-β j + 1 + j j ξ +1+j
where we have used that g -1 = γ. The coefficient associated to the term ξ N above gives the inequality

ẽN k ≤ k =0 α g k-1-β N --1 N + 1 , ∀k ≥ 1,
From the definition of ζ , we have that

g ≤ ∆T ε +2 G ( + 1)! max 0≤N ≤N (1 + u(T N ) ).
Therefore, recalling the definition (7) of α, β and γ, we derive

ẽN k ≤ ε k+1 G max 0≤N ≤N (1 + u(T N ) ) C d k =0 (C d ∆T ) +1 (1 + C c ∆T ) N --1 (k -)! N + 1 ≤ ε k+1 G max 0≤N ≤N (1 + u(T N ) ) C d k =0 C d T e -C C ∆T +1 e C C T ( + 1)! (k -1 -)! ≤ max 0≤N ≤N (1 + u(T N ) )e CcT C d (k + 1)! (1 + C d T e -Cc∆T )ε G k+1 ,
which ends the proof since E N k ≤ ẽN k for N = 0, . . . , N .

Practical realization of

[E(T N , ∆T, y N k ), ζ N k ]
Since the accuracy ζ N k needs to improve with k, the most natural way to build the approxima-

tions [E(T N , ∆T, y N k ), ζ N k ]
is with adaptive techniques and with adaptive refinements at every step k. The implementation ultimately rests on the use of a posteriori error estimators. It opens the door to local time step adaptation in the parareal algorithm as well as spatial coarsening or refinement if the problem involves additional spatial variables.

In principle, as ζ N k decreases with k, the numerical cost increases in terms of degrees of freedom and also in terms of computing time. This actually reveals the key idea of this new approach which is that we would like that only the last fine solver is expensive and the cost of the previous ones is a small fraction of the cost of the last one (we refer to the next sub-section for a more precise statement). By re-using information from previous iterations, we can limit the cost of internal solvers required in [E(T N , ∆T, y N k ), ζ N k ] and enhance the speed-up. This depends of course on the nature of the specific problem. We discuss several common situations in Appendix A.

Connection to the classical formulation of the parareal algorithm and advantages of the current view-point

In the original version of the algorithm, E(T N , ∆T, y N k ) is approximated with an accuracy ζ N k = ζ F which is kept constant in N and across the parareal iterations k. This has usually been done by using a solver F defined in the same spirit as G, but satisfying Hypothesis (4) with a better accuracy ε F < ε G . We have in this case

[E(T N , ∆T, y N k ); ζ F ] = F(T N , ∆T, y N k )
and we recover the classical algorithm (see [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF] and [START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF])

             y N +1 0 = G(T N , ∆T, y N 0 ), 0 ≤ N ≤ N -1 y N +1 k+1 = G(T N , ∆T, y N k+1 ) + F(T N , ∆T, y N k ) -G(T N , ∆T, y N k ), 0 ≤ N ≤ N -1, k ≥ 0, y 0 0 = u(0).
Compared to this classical version of the parareal algorithm, the adaptive approach offers the following important advantages:

1. The algorithm converges to the exact solution u(T N ) and not to the solution achieved by the fixed chosen fine solver F(0, T N , u(0)) (indeed, for any

N , ζ N k -→ 0 as k -→ ∞).
2. For a final target accuracy η, the parallel efficiency will always be superior to the classical approach (see section 3).

3. The numerical cost is designed to be near-minimal because we identify the minimal required accuracies at each iteration (see equation ( 12)). Early iterations use a loose tolerance, thus avoiding unnecessary work due to oversolving, while later iterations use tighter tolerances to deliver accuracy.

4. The dynamical refinements of the fine solver invite to incorporate adaptive solvers with a posteriori error estimators to the parareal scheme.

Parallel efficiency

It is difficult to give accurate a priori estimations for the speed-up and efficiency of the method due to its adaptive nature so the actual performance can only be established through relevant examples. In section 4, we give some results for the case of the Brusselator system. Despite this difficulty in estimation, we make some general remarks in this section, which aim primarily at highlighting the relevance of the cost of the coarse solver. The speed-up is defined as the ratio

speed-up AP/seq (η, [0, T ]) := cost seq (η, [0, T ]) cost AP (η, [0, T ]) (14)
between the cost to run a sequential fine solver achieving a target accuracy η with the cost to run an adaptive parareal algorithm providing at the end the same target accuracy η. The parallel efficiency of the method is then defined as the ratio of the above speed up with the number of processor which gives a target of 1 to any parallel solver:

f N 0 < • • • < f N K(η)
. Neglecting the communication delays, the cost of the adaptive solver is

cost AP (η, [0, T ]) = K(η) k=0 N -1 N =0 g N k + K(η)-1 k=0 N -1 N =0 f N k .
The classical parareal algorithm involves, at every iteration k ∈ {0, ..., K(η)} propagations at the highest accuracy ζ N K(η) . Thus its cost is

cost CP (η, [0, T ]) = K(η) k=0 N -1 N =0 g N k + K(η) N -1 N =0 f N K(η) ,
from which it directly follows that (at least if, with obvious notation, K AP (η) = K CP (η))

eff AP/seq (η, [0, T ]) > eff CP/seq (η, [0, T ]).
Thus the parallel performance of the adaptive algorithm is at least the one of the classical version. Note that this holds even when communications are not negligible since there is the same amount of information exchange in both algorithms.

We next give a more quantitative statement on an admittedly idealized setting. Assume that the cost of the coarse solve is negligible, that there is no communication delay and that the cost to realize [E(T N , ∆T,

y N k ), ζ N k ] is f N k = ∆T (ζ N k ) -1/α . ( 15
)
This assumption for the cost is, for instance, reasonable when we use an explicit time-stepping method of order α > 0. It would also hold for an implicit method where a direct solver can be used. Note that α could actually depend on k but we stick to this simple model for clarity of exposition.

Proposition 3.1. If f N k = ∆T (ζ N k ) -1/α
, for some α > 0 and if the cost of the coarse solver is negligible with respect to f N k for any k ≥ 0, then

eff AP/seq (η, [0, T ]) = 1 -τ 1/α 1 -τ K(η)/α ∼ 1 (1 + ε 1/α G ) . Therefore speed-up AP/seq (η, [0, T ]) ∼ N 1 (1 + ε 1/α G )
.

Proof. The cost of the scalable adaptive parareal scheme after K(η) iterations is

cost AP (η, [0, T ]) = ∆T K(η)-1 k=0 ζ -1/α k ( 16
)
Since we are in a range where the scheme converges, the quantity max 0≤N ≤N y N k is bounded and thus there exists 0 < c ≤ 1 ≤ c such that c ≤ ν k ≤ c for all k ≥ 0. We will account for this with the notation ν k ∼ 1. Note that in fact c and c are close to one. Let us start with the simple case α = 1 and denote K = K(η) -1:

cost AP (η, [0, T ]) = ∆T K k=0 ζ -1 k = ∆T ζ -1 K(η)-1   1 + ε G K + ε 2 G K(K -1) + • • • + ε K-1 G K!   we thus derive cost AP (η, [0, T ]) ≤ ∆T ζ -1 K(η)-1 (1 + ε G ) In the general case (α > 1), the same is true with cost AP (η, [0, T ]) ≤ ∆T ζ -1/α K(η)-1 (1 + ε 1/α G )
Up to this last factor, the current conclusion is that the global cost of the parareal procedure is equal to the last fine solver on each sub-interval with size ∆T (both the coarse and the previous fine propagations are negligible).

Since the accuracy that is obtained at the end of the parareal procedure (see [START_REF] Gaja | Parallel in time algorithms for nonlinear iterative methods[END_REF]) is of the same order as the accuracy provided with classical parareal solver (compare with (6)), it follows that if we now take the last target accuracy ζ -1/α K(η)-1 of the adaptive algorithm as the accuracy of the fine scheme in the classical parareal algorithm, the cost would be

cost CP (η, [0, T ]) = K(η)∆T ζ -1/α K(η)-1 , Therefore, cost AP (η, [0, T ]) cost CP (η, [0, T ]) ∼ 1 K(η) (1 + ε 1/α G ). ( 17
)
In addition, we know that when the cost of the coarse solver is negligible,

speed-up CP/seq (η, [0, T ]) = cost seq (η, [0, T ]) cost CP (η, [0, T ]) = N K(η) , ( 18
)
Dividing ( 18) by ( 17) yields speed-up AP/seq (η, [0, T ]) ∼ N 1

(1 + ε 1/α G ) and eff AP/seq (η, [0, T ]) ∼ 1 (1 + ε 1/α G )
.

In the ideal setting of Proposition 3.1:

• The parallel efficiency of the adaptive parareal algorithm does not depend on the final number of iterations. This is in contrast to the classical version whose efficiency decreases with the final number of iterations K(η) as 1/K(η).

• The efficiency behaves like 1 -o(ε G ) in the adaptive version, and o(ε G ) rapidly goes to zero with ε G . As soon as ε G becomes negligible with respect to 1, we will be in the range of full scalability.

We emphasize that, obviously, the above idealized setting will never hold in practice, but the result is interesting in its own right since it highlights that the cost of the fine solver is no longer the main obstacle for full scalability in the adaptive setting: the cost of the coarse solver becomes now the major obstruction towards full efficiency.

Numerical tests 4.1 Guidelines for a practical implementation

Practical choice of ζ N k

Formula (12) of the convergence analysis of section 2.3 gives an estimate for ζ N k that one could in principle use for the implementation. However, these tolerances may not be optimal because they are derived from a theoretical convergence analysis based on abstract conditions for the coarse and fine solvers. This was confirmed during our numerical tests where we observed that using estimates [START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF] for ζ N k did not deliver satisfactory enough results. This is the reason why it is necessary to devise a practical rule to set ζ N k . We have explored the following choice: if η is the final target accuracy, the classical parareal algorithm is usually run with a solver that delivers a slightly higher accuracy, say η/2. Assume that the classical algorithm converges in K CP (η) = K iterations. We propose to build the tolerances of ζ N k in such a way to target that K AP (η) = K CP (η) and such that the cost of the last fine propagation is of the order of the sum of the previous ones. This motivates to set

ζ N k =    ε 1-k+1 K G η 2 k+1 K , if k < K η/2, if k ≥ K.
The numerical example of the next section uses these tolerances.

Load balancing: For simplicity of exposition, the algorithm has so far been discussed for N subintervals of uniform size ∆T . However, this decomposition may lead to a task imbalance because some time intervals may have more complex dynamics than others, requiring more degrees of freedom, thus more computational time. In order to balance tasks as efficiently as possible, we dynamically adapt the size of the N subintervals in a way to have the fine solver propagations as balanced as possible among processors.

Results for a stiff ODE: the Brusselator system

We consider the brusselator system

x = A + x 2 y -(B + 1)x y = Bx -x 2 y,
with initial condition x(0) = 0 and y(0) = 1. This is a stiff ODE that models a chain of chemical reactions. It was already studied in a previous work on the parareal algorithm (see [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF]). The system has a fixed point at x = A and y = B/A which becomes unstable when B > 1 + A 2 and leads to oscillations. We place ourselves in this oscillatory regime by setting A = 1 and B = 3. The dynamics present large velocity variations in some time subintervals, making the use of adaptive time-stepping schemes particularly desirable for an appropriate treatment of the transient.

For the coarse solver, we set ε G = 0.1, and use an explicit Runge Kutta method of order 5 with an adaptive time-stepping (see [START_REF] Dormand | A family of embedded runge-kutta formulae[END_REF]).

For the fine solver, we use the implicit Runge-Kutta method of the Radau IIA family of order 5 with adaptive time-stepping (see [START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF][START_REF] Hairer | Stiff differential equations solved by radau methods[END_REF]). Both integrators are available in the ODE integration library of Scipy 1 which we have used to produce our results. As already discussed, the target accuracies η N k should be ensured by rigorous a posteriori error estimators. These type of estimators is unfortunately not available in the Scipy library and we are not aware of any mainstream library with this capability, so we have built a priori a "chart" mapping accuracies of the solver against certain tolerance parameters of the library. Once the chart is prepared, we run our parareal algorithm and adapt dynamically the parameters of our library solver to obtain desired accuracy η N k . We use formula [START_REF] Gander | 50 years of time parallel time integration[END_REF] to compare the speed-up of the classical and adaptive parareal algorithm in terms of the number of operations involved in the numerical solution (communication delays have not been taken into account). For the costs g N k and f N k , we take into account:

• the number of time steps (which is adaptively increased as we tightned the accuracy),

• the number of right-hand side evaluations,

• for the fine solver, we additionally count the number of evaluations of the Jacobian matrix and of the number of linear system inversions.

In Figure 1, we plot the obtained speed-up for different configurations:

• the final time T varies from 100 to 900,

• the final target accuracy is η = 10 -6 or η = 10 • the number of processors N varies from 10 to 100.

As anticipated in section 3, the speed-up of the adaptive parareal is always superior to the one of the classical parareal. We observe that the gain is marginal for a moderate accuracy (η = 10 -6 ) but it is about 2.5 times larger for η = 10 -8 . Note that sometimes the speed-up does not increase monotonically as the number of processors N increases. Also, the speed-up generally increases with N but the increase is rather moderate. We believe that this behavior could be improved with a rigorous a posteriori error estimator.

The values significantly differ from the range of full scalability and we next explain why this is mainly due to the cost of the coarse solver. Since the problem is stiff and we consider relatively long time intervals, it has been necessary to use a sufficiently accurate coarse solver. This explains our choice of an explicit Runge-Kutta scheme of order 5. To illustrate the impact of its cost, let us fix T = 500, η = 10 -8 and N = 50 (other parameters would yield similar conclusions). We compare the speed-up and efficiency when we count or do not count the cost of the coarse solver on Table 1. Obviously, when we do not count the cost of the coarse solver, the performance of both algorithms improves but it is particularly increased in the case of the adaptive version. If the cost of G was negligible, it would deliver a very satisfactory efficiency of 75.52%. This is five times larger than what the classical parareal would yield. This analysis illustrates that the major obstacle to achieve competitive scalabilities is no longer the cost of the fine solver like in the classical version, but the cost of the coarse propagator.

We next give some insight on the differences in the convergence behavior of both algorithms. We fix T = 10, η = 10 -8 and N = 10 and plot in Figure 2 the convergence history of the parareal solution in terms of:

Speed-up

Classical parareal

• the errors of the fine solver at every fine time-step

• the maximum error of the parareal solution at the macro-intervals

max N u(T N ) -y N k
Note that the maximum error in the adaptive scheme steadily decreases to the desired accuracy whereas the error in the classical scheme degrades at iteration k = 1 before converging. This type of behavior has been observed for all other configurations and we conjecture that an important difference in accuracy between the coarse and the fine solver at early stages of the algorithm may be the cause. Finally, an inspection of the error of the fine solver shows that the adaptive algorithm succeeds to reduce the error at every time t in a much more uniform way than the classical algorithm.

Conclusions and perspectives

The new adaptive formulation of the parareal algorithm opens the door to improve significantly the parallel efficiency of the method provided that the cost of the coarse solver is moderate. The increasing target tolerances which have to be met at each step allows to use online stopping criteria involving a posteriori estimators. The developed methodology remains theoretical since we have not quantified the impact of communication delays between processors nor potential memory issues (note however that the load balancing is devised with the purpose of equilibrating tasks and memory). In the framework of the ANR project "Ciné-Para (ANR-15-CE23-0019)", we are working on these issues that are of a different level of theory and involve different collaborators. This will be the topic of another paper. In addition to this, several extensions based on the current findings are subject of ongoing works, in particular, the coupling of the adaptive parareal with adaptive space-time schemes and the coupling of parareal with internal iterative solvers like in the discussion of section A.3 of the appendix.

Enriching the input information with previous iterations

Usually, solvers to realize [E(T N , ∆T, y N k ), ζ N k ] are built using only I N k = {T N , ∆T, y N k } as input information. We account for this idea with the notation

S(T N , ∆T, y N k ) -------→ (I N k ,cost N k ) [E(T N , ∆T, y N k ); ζ N k ],
and the numerical cost, denoted cost N k , increase as ζ N k is tightened. In this section, we discuss how to enhance the gain in efficiency of the adaptive algorithm by discussing ways to increase the accuracy of the solver [E(T N , ∆T, y N k ), ζ N k ] across the iterations while maintaining the cost to realize it as independent as possible from ζ N k , N , and k. For this, one possibility is to enrich I N k with data produced during the previous parareal iterations (although it would of course be at the cost of increasing the storage requirements).

Let P N k denote the intermediate information that has been produced at iteration k between [T N , T N +1 ] and by P k := ∪ N -1 N =0 P N k all the information produced at step k. Using

I N k = {I N k , P N -1 k , .
. . , P 0 k , P k-1 , . . . , P 0 }, as input information, the idea is to see whether it could be possible to find a solver S such that

S(T N , ∆T, y N k ) ------→ ( I N k ,cost) [E(T N , ∆T, y N k ); ζ N k ] (19)
with a constant and small complexity cost. Note that using the enriched set of information I N k means that we want to learn from the previous approximations of u(T N +1 ) given by [E(T N , ∆T, y N p ), ζ N p ], 0 ≤ p ≤ k -1, to start the current algorithm closer to u(T N +1 ). After each parareal iteration, we thus improve the accuracy, without increasing the work for solving because we start from a better input, accumulated from the previous parareal iterations.

In the rest of this section, we describe three relevant scenarios where we can approximate E(T N , ∆T, y N k ) by trying to build a scheme in the spirit of [START_REF] Lions | Résolution d'EDP par un schéma en temps pararéel[END_REF]. The first two examples have already been presented in the literature and concern the coupling of parareal with spatial domain decomposition (section A.1) and with iterative high-order time integration schemes (section A.2). In these two cases, there is to date no complete convergence analysis since it remains to show that i) E(T N , ∆T, y N k ) is approximated with accuracy ζ k and ii) the cost of the solver is really constant through the parareal steps. In addition to these two applications, we mention a third scenario where the convergence analysis can be fully proven. It concerns the solution of time-dependent problems involving internal iterative schemes at every time step. This idea was first analyzed in [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] in a restricted setting. In section A.3, we give the main setting and defer the analysis for a forthcoming paper.

A.1 Parareal coupled with spatial domain decomposition

Here, we consider a solver S = DDM which involves spatial domain decomposition over [T N , T N +1 ] [START_REF] Guetat | Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine[END_REF][START_REF] Aouadi | Convergence analysis of the coupled parareal-schwarz waveform relaxation method[END_REF]. We assume DDM involves a time discretization with a small time step δt < ∆T . Let n be the number of time steps on each interval [T N , T N +1 ] so that we have the relations ∆T = nδt and T = N ∆T = N nδt and the total number of time steps in [0, T ] is n := N n. The domain decomposition iterations act on a partition Ω = ∪ L l=1 Ω l of the domain. For 0 ≤ n ≤ n, we denote by u N,n,j k the solution produced by DDM at time t = T N + nδt after j ≥ 0 domain decomposition iterations. The notation J * will denote the last iteration (fixed according to some stopping criterion). At j = 0, these iterations need to be initialized at the interfaces ∂Ω l , 1 ≤ l ≤ L. The idea explored in, e.g., [START_REF] Guetat | Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine[END_REF][START_REF] Aouadi | Convergence analysis of the coupled parareal-schwarz waveform relaxation method[END_REF], is to take the values u N,n,J * k-1 | ∂Ω l at these interfaces as a starting guess for 0 ≤ n ≤ n so that

I N k = {I N k , {u N,n,J * k-1 | ∂Ω l , 0 ≤ n ≤ n}}.
From [START_REF] Guetat | Méthode de parallélisation en temps: Application aux méthodes de décomposition de domaine[END_REF][START_REF] Aouadi | Convergence analysis of the coupled parareal-schwarz waveform relaxation method[END_REF], there is numerical evidence that the computations of DDM(T N , ∆T, y N k ) using

I N k yield [E(T N , ∆T, y N k ); ζ N k ]
after a reduced number of iterations J * which is independent of k. Thus cost would be kept constant and

u N,n,J * k = DDM(T N , ∆T, y N k ) ------→ ( I N k ,cost) [E(T N , ∆T, y N k ); ζ N k ],
where the above "cost" is much small than the cost of the fine solver.

A.2 Parareal coupled with iterative high-order time integration schemes

Spectral Deferred Correction (SDC, [START_REF] Dutt | Spectral deferred correction methods for ordinary differential equations[END_REF]) is an iterative time integration scheme. Starting from an initial guess of u(t) at discrete points, the method adds successive corrections to this guess. The corrections are found by solving an associated evolution equation. Under certain conditions, the correction at every step increases by one the accuracy order of the time discretization.

We carry here a simplified discussion on how to build [E(T N , ∆T, y N k ); ζ N k ] when S = SDC and connect it to the so-called Parallel Full Approximation Scheme in Space-Time (PFASST, [START_REF] Minion | A hybrid parareal spectral deferred corrections method[END_REF][START_REF] Emmett | Toward an efficient parallel in time method for partial differential equations[END_REF]). For a given time interval [T N , T N +1 ], let us consider its n + 1 associated Gauss-Lobatto points {t N,n } n n=0 and quadrature weights {ω n } n n=0 . The Gauss-Lobatto points are such that t N,0 = T N and t N,n = T N +1 . Let us denote u N,n,j k the approximation of u(t N,n ) at parareal iteration k after j ≥ 0 SDC iterations. Assuming that one uses an implicit time-stepping scheme to solve the corrector equations involved in this method, u N,n+1,j k is given by u ]. An additional component of PFASST is that the algorithm also tries to improve the accuracy of the coarse solver G using SDC iterations built with I k N . This has not been taken into account in our adaptive parareal algorithm [START_REF] Emmett | Toward an efficient parallel in time method for partial differential equations[END_REF].

N,n+1,j k = u N,n,j k + (t N,n+1 -t N,n ) A(t N,n+1 , u N,n+1,j k ) -A(t N,

A.3 Coupling with internal iterative schemes

Any implicit discretization of problem (1) leads to discrete linear or nonlinear systems of equations which are often solved with iterative schemes. When they are involved as internal iterations within the parareal algorithm, one could try to speed them up by building good initial guesses based on information from previous parareal iterations. We illustrate this idea in the simple case where:

• A(t, •) is a linear differential operator in U complemented with suitable boundary conditions,

• we use an implicit Euler scheme for the time discretization.

A sequential solution of problem (1) with the implicit Euler scheme goes as follows. At each time t N,n = T N + nδt, the solution u(t N,n ) is approximated by the function u N,n ∈ U which is itself the solution to B(u N,n ) = g N,n , where g N,n = u N,n-1 + δtf (t N,n ) and B(v) := v + δtA(t N,n , v), ∀v ∈ U.

Note that B depends on time but our notation does not account for it in order not to overload the notations. After discretization of U, the problem classically reduces to solving a linear system of the form Bū N,n = ḡN,n , for the unknown ūN,n in some discrete subspace S of U. Usually, the above system is solved either by means of a conjugate gradient method or by a Richardson iteration of the form ūN,n,j = (Id + ωPB)ū N,n,j-1 + ωP b, j ≥ 1 ūN,n,0 ∈ S given.

Here, ω is a suitably chosen relaxation parameter and P can be seen as a pre-conditioner. The internal iterations j are stopped whenever a certain criterion is met (it could be an a posteriori estimator) and we denote by J N,n their final number. Obviously, J N,n depends on the starting guess for which a usual choice is to take the solution at the previous time, that is ūN,n,0 = ūN,n-1,J N,n-1 .

In order to achieve our goal, i.e. maintaning a low cost while increasing the accuracy at each parareal step, we can now reuse information from previous parareal iterations for the starting guess. In [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF], two options are explored. The first is , if k ≥ 1.

In the first case, we take over the internal iterations at the point where they were stopped in the previous parareal iteration k -1. In addition to this, in the second case, the term u

N,n-1,J N,n,k k -u N,n-1,J N,n-1,k-1 k-1
tries to better take the dynamics of the process into account. Note that the use of solutions that have been produced in the previous parareal iterations is at the expense of additional memory requirements. It might also be at the cost of a certain increase in the complexity locally at certain times. However, [START_REF] Mula | Some contributions towards the parallel simulation of time dependent neutron transport and the integration of observed data in real time[END_REF] shows in a restricted setting that these starting guesses (in particular the second) have interesting potential to enhance the speed-up of the parareal algorithm. A general theory on this aspect will be presented in a forthcoming work.

Theorem 2 . 1 .

 21 u(T N ) ), and τ := C d T e -Cc∆T ε G . If G and δG satisfy Hypothesis (4), then,

eff

  AP/seq (η, [0, T ]) := speed-up AP/seq (η, [0, T ]) N . Assume that g N k and f N k are the numerical costs to realize G(T N , ∆T, y N k ) and [E(T N , ∆T, y N k ), ζ N k ]. Since the tolerances ζ N k decrease with k, we have

Figure 1 :

 1 Figure 1: Speed-up as a function of the number of processors N . Dashed lines: classical parareal. Continuous lines: Adaptive parareal.

Figure 2 :

 2 Figure 2: Convergence history of the errors. Top: classical parareal. Bottom: adaptive parareal. Left: Errors for t ∈ [0, 10]. Right: maximum error at each iteration k.

  n+1 , u N,n+1,j(t N,m , u N,m,j-1 k), 1 ≤ j, 0 ≤ n ≤ n -1, ≤ n ≤ n.To speed-up computations, one of the key elements is the choice of the starting guesses u N,n,0 k , 0 ≤ n ≤ n. Without entering into very specific details, the PFASST algorithm is a particular instantiation of the above scheme when J * = 1 and u N,n,0 k uses information produced at the previous parareal iteration k -1. Therefore, PFASST falls into the present framework in the sense that it produces uN,n,1 k = SDC(T N , ∆T, y N k ) with I N k = {T N , ∆T, y N k , P k-1 }and it is expected that uN,n,1 k = [E(T N , ∆T, y N k ); ζ N k

1 ,

 1 -1,J N,n-1,k k , if k = 0 ūN,n,0 k = ūN,n,J N,n,k-1 k-if k ≥ 1,and the second, less natural choice,

Table 1 :

 1 Impact of the cost of the coarse solver. Speed-up and efficiency with N = 50.

	Adaptive Parareal

Acknowledgment

This work was funded by the ANR project "Ciné-Para" (ANR-15-CE23-0019).


Nicolas F Bouché

Kevin Mcconway

Melatonin Levels and Low-Frequency Magnetic Fields in Humans and Rats: Reconciling Opposite Results with a Bayesian Regression

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since the epidemiological study of [START_REF] Wertheimer | Original Contributions[END_REF], concerns for an adverse health effect (in particular for childhood leukemia) due to electrical and magnetic fields (MFs) generated in the Extremely Low Frequency [ELF] regime (< 300 Hz, but mostly at 50-60 Hz) by power lines have been raised in the west and also from case-reports of electrical substation workers in the former Soviet Union [e.g. [START_REF] Zhadin | Review of russian literature on biological action of DC and low-frequency AC magnetic fields[END_REF]. This potential association between residential exposure to ELF Magnetic Fields (ELF-MF) and childhood leukemia has remained from the various pooled analysis of the numerous epidemiological studies [START_REF] Ahlbom | A pooled analysis of magnetic fields and childhood leukaemia[END_REF][START_REF] Savitz | Health Effects of Electric and Magnetic Fields: Are We Done Yet[END_REF][START_REF] Draper | Childhood cancer in relation to distance from high voltage power lines in England and Wales: a case-control study[END_REF][START_REF] Kheifets | Childhood leukemia and EMF: Review of the epidemiologic evidence[END_REF][START_REF] Schüz | Nighttime exposure to electromagnetic fields and childhood leukemia: An extended pooled analysis[END_REF][START_REF] Sermage-Faure | Childhood leukaemia close to high-voltage power lines -the Geocap study, 2002-2007[END_REF][START_REF] Kheifets | Epidemiologic study of residential proximity to transmission lines and childhood cancer in California: description of design, epidemiologic methods and study population[END_REF][START_REF] Schüz | Extremely low-frequency magnetic fields and risk of childhood leukemia: A risk assessment by the ARIMMORA consortium[END_REF] revealed that the relative risk for leukemia is 2× for MF of intensities ≥ 0.4 µT. This elevated risk for childhood leukemia has led to the World Health Organization to label ELF MF as possible carcinogen 'class 2B' based on International Agency for Research on Cancer (IARC) report on the subject [IARC, 2002], a conclusion recently re-affirmed by the IARC chair on non-ionizing radiation [START_REF] Schüz | Extremely low-frequency magnetic fields and risk of childhood leukemia: A risk assessment by the ARIMMORA consortium[END_REF].

Historically, [START_REF] Stevens | The Melatonin Hypothesis: Electric Power and Breast Cancer[END_REF] proposed the so-called melatonin hypothesis in the context of breast cancer involving ELF MFs discussed in the 90s [as reviewed in [START_REF] Brainard | The relationship between electromagnetic field and light exposures to melatonin and breast cancer risk[END_REF][START_REF] Kliukiene | Residential and occupational exposures to 50-hz magnetic fields and breast cancer in women: A population-based study[END_REF]. Under this hypothesis, the well-known melatonin (MLT) hormone produced by the pineal gland that controls the body's sleep/wake cycle [e.g. [START_REF] Reiter | Action Spectra, Dose-Response Relationships, and Temporal Aspects of Light's Effects on the Pineal Gland[END_REF][START_REF] Reiter | Melatonin: the chemical expression of darkness[END_REF], would be an intermediary agent where ELF MF would somehow impact MLT levels and this in turn would increase the risk of developing a disease or cancer. This hypothesis was put forward because (i) it was known by [START_REF] Stevens | The Melatonin Hypothesis: Electric Power and Breast Cancer[END_REF] that somehow the pineal gland responds to artificial EMF [since the 80s: [START_REF] Stemm | Effects of an earth-strength magnetic field on electrical activity of pineal cells[END_REF][START_REF] Wilson | Neuro-Endocrine Mediated Effects of Electromagnetic-Field Exposure -Possible Role of the Pineal-Gland[END_REF][START_REF] Wilson | Evidence for an Effect of Elf Electromagnetic-Fields on Human Pineal-Gland Function[END_REF][START_REF] Reiter | Alterations of the Circadian Melatonin Rhythm by the Electromagnetic Spectrum: A Study in[END_REF][START_REF] Reiter | Static and Extremely Low Frequency Electromagnetic Field Exposure: Reported Effects on the Circadian Production of Melatonin[END_REF][START_REF] Reiter | Melatonin suppression by static and extremely low frequency electromagnetic fields: relationship to the reported increased incidence of cancer[END_REF], and (ii) because MLT is an effective anti-oxident agent, free radical scavenger, and a potent oncostatic agent [e.g. [START_REF] Panzer | The validity of melatonin as an oncostatic agent[END_REF][START_REF] Allegra | The chemistry of melatonin's interaction with reactive species[END_REF][START_REF] Rodriguez | Regulation of antioxidant enzymes: a significant role for melatonin[END_REF][START_REF] Henshaw | Do magnetic fields cause increased risk of childhood leukemia via melatonin disruption?[END_REF][START_REF] Jung | Melatonin in Cancer Management: Progress and Promise[END_REF][START_REF] Reiter | Melatonin as an antioxidant: under promises but over delivers[END_REF]. Thus reduced MLT levels could lead to an increase risk of cancer [e.g. [START_REF] Guénel | Exposure to 50-Hz electric field and incidence of leukemia, brain tumors, and other cancers among French electric utility workers[END_REF][START_REF] Kliukiene | Residential and occupational exposures to 50-hz magnetic fields and breast cancer in women: A population-based study[END_REF][START_REF] Koeman | Occupational extremely low-frequency magnetic field exposure and selected cancer outcomes in a prospective Dutch cohort[END_REF] and other neurodegenerative illnesses [e.g. [START_REF] Feychting | Occupational magnetic field exposure and neurodegenerative disease[END_REF][START_REF] Huss | Source of Funding and Results of Studies of Health Effects of Mobile Phone Use : Systematic Review of Experimental Studies[END_REF][START_REF] Davanipour | Severe Cognitive Dysfunction and Occupational Extremely Low Frequency Magnetic Field Exposure among Elderly Mexican Americans[END_REF] by increasing the oxidative stress as described in [START_REF] Mevissen | Complex effects of long-term 50 Hz magnetic field exposure In Vivo on immune functions in female Sprague-Dawley rats depend on duration of exposure[END_REF] and reviewed in [START_REF] Consales | Electromagnetic Fields, Oxidative Stress, and Neurodegeneration[END_REF].

This hypothetical connection made by [START_REF] Stevens | The Melatonin Hypothesis: Electric Power and Breast Cancer[END_REF] between circadian rhythm disruption and certain illnesses has been revisited in the context of childhood leukemia by [START_REF] Henshaw | Do magnetic fields cause increased risk of childhood leukemia via melatonin disruption?[END_REF]. While this connection between MLT levels and ELF-MF lacked a clear mechanism, it seems to be related to the visual system since rats with severed optical nerves not longer respond to ELF-MF [START_REF] Olcese | Evidence for the involvement of the visual system in mediating magnetic field effects on pineal melatonin synthesis in the rat[END_REF]. The exact mechanism with magneto-receptors in the retina is now a plausible scenario in light of recent developments in the study of magneto-reception from behavioural [e.g. [START_REF] Kirschvink | Is Geomagnetic Sensitivity Real? Replication of the Walker-Bitterman Magnetic Conditioning Experiment in Honey Bees[END_REF][START_REF] Phillips | Behavioural evidence for use of a light-dependent magnetoreception mechanism by a vertebrate[END_REF][START_REF] Ritz | Resonance effects indicate a radical-pair mechanism for avian magnetic compass[END_REF][START_REF] Ritz | Magnetic Compass of Birds Is Based on a Molecule with Optimal Directional Sensitivity[END_REF][START_REF] Johnsen | The physics and neurobiology of magnetoreception[END_REF][START_REF] Gegear | Cryptochrome mediates light-dependent magnetosensitivity in Drosophila[END_REF][START_REF] Malkemper | Magnetoreception in the wood mouse (Apodemus sylvaticus): influence of weak frequency-modulated radio frequency fields[END_REF][START_REF] Yoshii | Cryptochrome mediates light-dependent magnetosensitivity of Drosophila's circadian clock[END_REF][START_REF] Bazalova | Cryptochrome 2 mediates directional magnetoreception in cockroaches[END_REF][START_REF] Wiltschko | Two Different Types of Light-Dependent Responses to Magnetic Fields in Birds[END_REF][START_REF] Winklhofer | Avian magnetic compass can be tuned to anomalously low magnetic intensities[END_REF][START_REF] Wiltschko | Sensing magnetic directions in birds: Radical pair processes involving cryptochrome[END_REF][START_REF] Wiltschko | Light-dependent magnetoreception in birds: the crucial step occurs in the dark[END_REF][START_REF] Sherrard | Low-intensity electromagnetic fields induce human cryptochrome to modulate intracellular reactive oxygen species[END_REF] and theoretical investigations [e.g. Ritz et al., 2010a;[START_REF] Hore | The Radical-Pair Mechanism of Magnetoreception[END_REF] where the cryptochrome CRY proteins discovered in the 90s [START_REF] Ahmad | HY4 gene of A. Thaliana encodes a protein with characteristics of a blue-light photoreceptor[END_REF][START_REF] Ahmad | Magnetic intensity affects cryptochrome-dependent responses in Arabidopsis thaliana[END_REF][START_REF] Ahmad | Seeing the world in red and blue: insight photoreceptors Margaret Ahmad[END_REF][START_REF] Chasmore | Cryptochromes: Blue Light Receptors for Plants and Animals[END_REF][START_REF] Chaves | The Cryptochromes: Blue Light Photoreceptors in Plants and Animals[END_REF][START_REF] Ahmad | Photocycle and signaling mechanisms of plant cryptochromes[END_REF] would provide the radical pair mechanism postulated by [START_REF] Schulten | A biomagnetic sensory mechanism based on magnetic field modulated coherent electron spin motion[END_REF] and be the (light-dependent) MF receptor [Ritz et al., 2010b;[START_REF] Liedvogel | Cryptochromes-a potential magnetoreceptor: what do we know and what do we want to know[END_REF][START_REF] Hore | The Radical-Pair Mechanism of Magnetoreception[END_REF][START_REF] Michael | Invited Review Animal Cryptochromes: Divergent Roles in Light Perception, Circadian Timekeeping and Beyond †[END_REF]. CRY proteins are widely expressed in cones, amacrine cells of the retina [e.g. [START_REF] Foley | Human cryptochrome exhibits light-dependent magnetosensitivity[END_REF][START_REF] Wong | differential roles for cryptochromes in the mammalian retinal clock[END_REF] and are thought to be the prime MF receptors involved in avian compass.

As discussed in [START_REF] Lagroye | ELF magnetic fields: Animal studies, mechanisms of action[END_REF], CRYs which are ubiquitous, and recently discovered (blue) light dependent magneto-photoreceptor, should be assessed as a plausible mechanism behind some of the bioeffects of ELF MFs. CRYs are also involved in the regulation of circadian biorhythms [e.g. [START_REF] Van Der Horst | Mammalian Cry1 and Cry2 are essential for maintenance of circadian rhythms[END_REF][START_REF] Yoshii | Cryptochrome mediates light-dependent magnetosensitivity of Drosophila's circadian clock[END_REF][START_REF] Ono | Cryptochromes are critical for the development of coherent circadian rhythms in the mouse suprachiasmatic nucleus[END_REF][START_REF] Wong | differential roles for cryptochromes in the mammalian retinal clock[END_REF], which led Vanderstraeten et al. [Vanderstraeten and Burda, 2012;Vanderstraeten et al., 2012[START_REF] Vanderstraeten | Could Magnetic Fields Affect the Circadian Clock Function of Cryptochromes? Testing the Basic Premise of the Cryptochrome Hypothesis (ELF Magnetic Fields)[END_REF][START_REF] Vanderstraeten | Magnetic Fields and Health: From Epidemiology To Cryptochrome Chemistry[END_REF] to revive the MLT hypothesis for childhood leukemia and to formulate the cryptochrome hypothesis in the context of the epidemiological results cited above [see also [START_REF] Lagroye | ELF magnetic fields: Animal studies, mechanisms of action[END_REF][START_REF] Juutilainen | Magnetocarcinogenesis : is there a mechanism for carcinogenic effects of weak magnetic fields ?[END_REF]. Under this hypothesis, weak MFs in the micro-tesla range disrupt the biorhythms, leading to disrupted MLT production rendering MLT as an effective marker to be used in relation to weak MFs. Moreover, it has been shown that pulsed MFs (PMF) can also stimulate a rapid accumulation of reactive oxygen species (ROS) -a metabolite implicated in stress response and cellular ageing-but only in insect cells expressing CRY [START_REF] Sherrard | Low-intensity electromagnetic fields induce human cryptochrome to modulate intracellular reactive oxygen species[END_REF] leading [START_REF] Landler | Cryptochrome : The magnetosensor with a sinister side ?[END_REF] to postulate that carcinogenesis associated with power lines, PMF-induced ROS generation, and animal magnetoreception share a common mechanism.

However, the epidemiological and laboratory studies on MLT levels and ELF MF are often contradictory [see reviews by [START_REF] Henshaw | Do magnetic fields cause increased risk of childhood leukemia via melatonin disruption?[END_REF][START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF][START_REF] Touitou | The effects of extremely low-frequency magnetic fields on melatonin and cortisol, two marker rhythms of the circadian system[END_REF][START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF][START_REF] Lewczuk | Influence of Electric , Magnetic , and Electromagnetic Fields on the Circadian System : Current Stage of Knowledge[END_REF]. In this paper, we use both a Bayesian analysis and a non-parametric approach on a compilation of 62 studies the evolution of MLT levels on humans and rats exposed under weak ELF MFs. Given that these studies are often inconsistent (in reporting variation or no changes in MLT levels), we are making sure to include both types of results.

Materials and Methods

Here, we present the compilation of 62 studies reporting MLT levels on humans ( § 2.1) and rats ( § 2.2) and our Bayesian methodology in ( § 2.3).

Melatonin Data on humans

Halgamuge [2013] compiled various studies on humans exposed to ELF published in the last 15-20 years where MLT levels -mostly urine (24hr)-were reported. These authors included both laboratory (short term) and epidemiological (long term) studies. From their collection of 33 studies, we noticed that some were duplicates, which were removed. We verified each entry listed in [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF] (their Table 4) regarding MF field strength and exposure duration, leading to some differences between their listing and this work.

Since our focus is on studying the putative effect of environmental/ambient (i.e. large scale) magnetic field on human MLT levels, mostly from power lines where the entire body is subject to the MF, we did not consider studies that involved very localized ELF such as those from electric blankets, video displays, nor cell phone usage. Furthermore, we did not consider those regarding geomagnetic storms, in-vitro studies, or involving static magnetic field.

Table 1 lists the studies from [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF] used in this analysis. Note that it includes the 14 studies listed in the review of [START_REF] Henshaw | Do magnetic fields cause increased risk of childhood leukemia via melatonin disruption?[END_REF], and we included a few studies not included in the original review of [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF] but in the review article of [START_REF] Touitou | The effects of extremely low-frequency magnetic fields on melatonin and cortisol, two marker rhythms of the circadian system[END_REF], such as [START_REF] Griefahn | Experiments on the effects of a continuous 16.7 Hz magnetic field on melatonin secretion, core body temperature, and heart rates in humans[END_REF]; [START_REF] Kurokawa | Acute Exposure to 50 Hz Magnetic Fields with Harmonics and Transient Components: Lack of Effects on Nighttime Hormonal Secretion in Men[END_REF]; [START_REF] Cocco | Urinary 6-sulfatoxymelatonin excretion in humans during domestic exposure to 50 hertz electromagnetic fields[END_REF]; [START_REF] Davis | Effects of 60-Hz Magnetic Field Exposure on Nocturnal 6-Sulfatoxymelatonin, Estrogens, Luteinizing Hormone, and Follicle-Stimulating Hormone in Healthy Reproductive-Age Women: Results of a Crossover Trial[END_REF] and [START_REF] Warman | Acute Exposure to Circularly Polarized 50-Hz Magnetic Fields of 200-300 µT Does Not Affect the Pattern of Melatonin Secretion in Young Men[END_REF]. The only study we rejected is that of [START_REF] Touitou | Magnetic fields and the melatonin hypothesis: a study of workers chronically exposed to 50-Hz magnetic fields[END_REF] which is based on a small sample (30) of individuals/electrical workers preselected not to have any sleep disturbances, i.e. is biased against finding any sleep/MLT perturbation from ELF MF. We note that last study is in contrast to the recent work of [START_REF] Liu | Occupational Electromagnetic Field Exposures Associated with Sleep Quality: A Cross-Sectional Study[END_REF] on 854 workers showing an increase of sleep disturbance in some utility workers [see also Monazzam et al., 2014, on this subject].

In Table 1, some studies claim that MLT levels are affected, but the change is invariably in the sense of a decrease of the MLT production or a phase shift. In contrast, other studies claim that the MLT level is not affected by MFs. The effect/no-effect outcome naturally leads to logistic modelling (described in § 2.3) appropriate for such binary situations [START_REF] Hosmer | Applied logistic regression[END_REF]. The logistic approach makes no implicit assumption and is simpler than invoking a model that assumes a linear relation between MLT levels, exposure duration, etc. Unfortunately, MLT studies are heterogeneous and there is no universal way to quantify the amount of decrease in MLT production across these studies. Thus, we assign the outcome of the studies listed in Table 1 a 1(0) depending on whether the original authors reported change (no change) in MLT levels, respectively. When the study reported 'some' change, we assign the outcome of the study a 0.5. This would correspond to, for instance, when changes were observed only for a sub-group of the study. [START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF] compiled various laboratory studies on the putative effect of ELF MF on rat MLT, whose list is reproduced in Table 2. We removed the entries that were not consistent with the original study, e.g. the entries with ID 13,14 and 15 from [START_REF] John | 60 Hz magnetic field exposure and urinary 6sulphatoxymelatonin levels in the rat[END_REF]. In addition, we added the study of [START_REF] Löscher | Effects of Weak Alternating Magnetic Fields on Nocturnal Melatonin Production and Mammary Carcinogenesis in Rats[END_REF] and [START_REF] Löscher | Exposure of Female Rats to a 100-µT 50 Hz Magnetic Field Does Not Induce Consistent Changes in Nocturnal Levels of Melatonin[END_REF]. As in § 2.1 we assign rat studies a 1(0) depending on whether the authors reported change (no change) in MLT levels, respectively. [START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF] investigated whether the MF exposure duration, MF polarization and other factors play a role. They concluded that the only factor that seemed to be the most significant is the duration of exposition to ELF MF, albeit with a P-value of 0.07 implying that this factor is not significant at more than >95% level, using a model linear with exposure duration and with field strength.

Melatonin Data on rats

A parametric bayesian analysis

As discussed in [START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF], logistic regression (LR) is a statistical technique commonly used to examine the possible relationship between a dichotomous-dependent variable (here the effect/non-effect on MLT excretion pattern) and independent variables (such as frequency, polarization, exposure duration, and MF). In general, the probability P to observe an effect (i.e. Y = 1) is given by the logistic function:

P (Y = 1) = 1 1 + exp(-t) = exp(t) 1 + exp(t) ≡ LR(t) (1)
where t is usually taken to be a linear combination of the dependent variables, X n , i.e.

t = α + β 1 X 1 + • • • + β n X n .
However, one should keep in mind that such a linear combination of dependent variables makes a critical assumption: namely that these variables are independent of one another. In other words, the probability to have an effect might depend on the field exposition duration and on the magnetic field strength, but the coefficient for each of these variables are assumed to be independent of one another.

In this work, we use a logistic function (L(t); Eq. 1) where t can be a non-linear function of the independent variables. L(t) gives the probability to observe an effect (p ≡ L(t)), and the observed realization is given by the Bernoulli (Bern) probability distribution since the observables are dichotomous, with values at 0 or 1, which can be written as (see Supplementary Material):

t = f (X i ; θ) (2) p = L(t) (3) O ∼ Bern(p) (4)
where O are the simulated observables.

The LR model is made robust to spurious data by including an (unknown) outlier fraction π, i.e. Eq. 3 becomes

p = π p out + (1 -π) L(t) (5)
where p out is the logistic probability for outliers and π is taken from a Uniform distribution from 0 to 0.5. We use uniform priors on π and p out .

Next, we will consider the following two parametric LR models. First, we use a model linear in exposure duration with log T as the single independent variables. Then, we will use a variant of the logistic model where the slope α is a function of the MF strength in a dichotomous fashion for reasons that will be clearer in § 3.4. The parametric models are Model A:

t ≡ α T (log 10 T -β T ) (6) Model B: t ≡ α (log 10 T -β) if B < B t γ (log 10 T -β) if B ≥ B t (7)
where α, γ are the linear slope, β the transition point of the logistic function, and B t is the threshold level for model B.

In order to find the best parameters θ for our model, we use a Monte Carlo Markov Chain (MCMC) algorithm. Because traditional MCMC algorithms are somewhat sensitive to the step size and the desired number of steps. In what follows, we use the No-U Turn Sampler (NUTS) of [START_REF] Hoffman | The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo[END_REF], a self-tuning variant of Hamiltonian Monte Carlo (HMC), except when the model is not continuous (as in model B) where we revert to the traditional Metropolis-Hasting sampling method. We typically use 2 MCMC chains per run and 15,000 iterations to 25,000 iterations per chain.

A non-parametric analysis

In order to investigate the inherent limitations of our parametric approach (as in any regression), we applied a non-parametric supervised classification algorithm to the data sets in order to determine whether there are robustly-defined regions in the parameter space that divide studies reporting a change in MLT levels with those that reported no change. We chose to apply the Support Vector Classification (SVC) algorithm [START_REF] Cortes | Support-Vector Networks[END_REF] implemented in the SVM module of the scikit-learn python package v0.19.1 [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Nonlinear regions were investigated, using a Gaussian RBF (Radial Basis Function) kernel which uses a Gaussian similarity measure between points in the parameter space. The use of the RBF kernel depends on two quantities, C, which describes the way in which smoothness of the boundaries of the classification regions in parameter space is traded off with misclassifications of the studies, and γ, which defines how much influence each individual study has. We used γ to be 1/number of features, and used a cross-validation technique to determine that C 3. Unlike the logistic regression methods, SVC does not directly give a probability of classification to each of the studies considered.

Results

Magnetic field strengths

Figure 1 shows the histogram of mean MF strengths for the studies compiled on humans in Table 1 (hatched) and on rats in Table 2 (solid). The strength of the static Earth magnetic field B ∼ 50µT is indicated with the vertical dotted line, but the local strength varies from ∼ 30 to 60µT, depending on the latitude.

This figure shows that human studies cover the range of MF of strength from 0.1 to 50µT, while rat studies are involving MF of higher strengths from 1 to 1000 µT. The MF distributions for human and rat studies appear to be significantly different as a KS-test indicate the two histograms are not drawn from the same parent population, with a P-value of 0.01. This difference is perhaps due to an implicit bias induced by researchers looking to bring out a signal in the lab, i.e. induced by a dose-response expectation [as in [START_REF] Warman | Acute Exposure to Circularly Polarized 50-Hz Magnetic Fields of 200-300 µT Does Not Affect the Pattern of Melatonin Secretion in Young Men[END_REF].

Results on human studies

Regarding model A (described in § 2.3), we use the following (uninformative) priors for the slope α and zero-point β:

α T ∼ N (0, 10)| up=5 lo=-5 β T ∼ U(-2, 3.5) (8)
where the N (µ, σ) is the normal distribution truncated on the interval [-5,5] and U is the Uniform distribution.

The best fit parameters of model A with their 95% credible intervals are β = 1.3 +1.6 -1.0 and α = 2.5 +2.3 -1.8

(Table 3).

Figure 2(left) shows the result from the LR model A applied on the 28 human studies reporting change or no change to human MLT levels ( § 2.1). The top panel shows the data in the plane log B-log T where the model predictive values are represented by the grey scale. The vertical line represents the best fit β parameter, i.e. where the probability to have an effect is modeled to be 0.5. The bottom panel shows the model prediction (red solid line) as a function of exposure duration T where vertical dotted-dashed lines indicating a day, a month, and a year. The shaded gray region represents the 95% posterior uncertainties, calculated using the [START_REF] Wilson | Probable inference, the law of succession, and statistical inference[END_REF] score confidence interval for binomial distributions, verified to be a continuous representation of the uncertainties found from the MCMC posteriors. This figure shows that ELF-MFs start to have an effect on MLT levels with a probability larger than 50% at around ∼22 days.

Figure 2(right) shows the same 28 human studies where we applied a non-parametric SVC algorithm and studies reporting change, partial change and no change on MLT levels are shown in red, yellow and blue respectively. This figure confirms that the studies reporting changes in MLT levels are predominantly in the region of parameter space with long exposure duration, supporting the results from the parametric LR shown on the left panel, and is not driven by a few rogue false data points.

Results on rat studies

We performed a similar analysis on the studies available on laboratory rats (described in § 2.2) and the results are listed in Table 4. One notable diff erence between studies involving humans or rats, is that the duration coefficient α T appears to be much weaker in the case of rat studies (α T 1.2) than in the case of humans (α T 2.5) and α T = 0 is much less significant for rats. However, we remind the reader that, as shown in Fig. 1, only a handful of human studies have MF strength above ∼50 µT, while about half of the studies on rats have MFs above this level.

Towards a unified frame work

Given that human and rat studies differ significantly in the field strengths, we show in Figure 3 the results for studies on laboratory rats when the MF strength is below (above) 45µT1 , shown in the bottom (top) panels, respectively. These two panels clearly show that the effect on MLT levels becomes random with respect to exposure duration T when the MFs are above ∼50µT. In both panels, the red solid line represents the best model (model A) obtained from the LR Bayesian analysis whose parameters are listed in Table 4.

Comparing Table 4 with the results of model A on humans Table 3, one sees that the statistics of change/no-change on MLT levels in rat and human studies are consistent with each other, only after restricting animal and human studies over the same range of MF strengths. The time-dependent factor is α T = 2.5 +2.3 -1.8 for human studies and α T 3.0 +1.9 -2.5 for rat studies. Furthermore the exposure duration where the MF exposure becomes significant (with a probability to affect MLT levels greater than 50%) is in both cases at β T ∼ 1.0, corresponding to ∼ 10 days.

Inspired by these results, we extended our LR model to include some (unknown) threshold MF, B t , i.e. model B introduced in § 2.3. Figure 4(left) shows the data in the log T -log B plane along with the model predictions respresented as the grey scale. Figure 4 (right) shows the non-parametric SVC analysis, and strongly supports the results from the Bayesian parametric LR. Figure 5 shows that posterior distribution on each of the parameters for model B, whose best parameters are listed in Table 4.

The best threshold value determined by the data is log B t = 1.4 +0.7 -0.5 (68% CL), i.e. the magnetic threshold is B t 10-65µT2 . We note that the transition field strength of ∼ 50µT corresponds to two different regimes, one where the ELF MF are a mere perturbation to the ambient static terrestrial MF, which has an amplitude of B 50µT, and the other where time varting ELF MF are the sole dominant contribution. We will return to this in the discussion ( § 3.6).

Model selection

In order to determine the model that performs the best, several information criteria have been proposed that penalize complex models with large number of parameters p [e.g. [START_REF] Raftery | Bayesian Model Selection in research[END_REF][START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF][START_REF] Sharma | Markov Chain Monte Carlo Methods for Bayesian Analysis in Astronomy[END_REF] such as: Schwartz's Bayesian Information Criterion (BIC= -2 ln L( θ) + p log n), Akaike's Information Criterion (AIC= -2 ln L( θ) + 2 p), the corrected (for small sample size) Akaike Information Criterion (AICc= -2 ln L( θ) + 2 p + 2p(p + 1)/Np -1) and the popular Deviance Information Criterion (DIC= -2 ln L( θ) + 2 p D ) introduced by [START_REF] Spiegelhalter | Bayesian Measures of Model Complexity anf Fit[END_REF]. As discussed in [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF] and [START_REF] Sharma | Markov Chain Monte Carlo Methods for Bayesian Analysis in Astronomy[END_REF], AIC and DIC favors the model with the best predictive performance while BIC favors the model with the smallest degrees of freedom.

However, these information criterion (BIC, AICs, DIC) are applicable only when the likelihood can be approximated by a continuous function (which is approximated by a Gaussian). When, this condition is not met, as in our case due to the discontinuity in model B, one can use the Widely Available Information Criteria [WAIC, from [START_REF] Watanabe | A Widely Applicable Bayesian Information Criterion[END_REF] or the Leave-One-Out (LOO) Cross-Validations [e.g. [START_REF] Vehtari | Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC[END_REF] technique which asymptotes to WAIC [START_REF] Watanabe | Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory[END_REF]. Table 5 lists the numerical values for these information criterion for the models explored in this paper and both WAIC and LOO-CV indicates that model B is slightly favored over model A.

Discussion

Our result of a threshold-dependent impact of man-made ELF MFs on MLT levels at intensities at or below B calls for a possible role of the geomagnetic field. Indeed, the amplitude of the static Earth MF |B | is not constant with time as there are fluctuations on a range of time-scales, from daily fluctuations to monthly, annual variations and up to time-scales of millions of years [see e.g. [START_REF] Courtillot | Time variations of the earth's magnetic field -From daily to secular[END_REF] due to complex interactions between the solar wind and the magnetosphere. The daily variations are of the order of 20 to a few hundreds of nT (i.e. 1000× smaller than the field strength) due to the impact of the solar wind pressure in the upper atmosphere [e.g. [START_REF] Hitchman | The quiet daily variation in the total magnetic field: Global curves[END_REF], and this led [START_REF] Liboff | Why are living things sensitive to weak magnetic fields?[END_REF] suggested that the biological genesis for interactions between living beings and weak ELF could originate from these tiny (∼50 nT) daily swing in the geomagnetic field because it is a remarkably constant effect exactly in phase with the solar diurnal change. Hence, as argued in [START_REF] Liboff | Why are living things sensitive to weak magnetic fields?[END_REF], the widespread sensitivity of biological systems to weak ELF magnetic fields is derived from the diurnal geomagnetic variations. However, while numerous studies show that MF can influence the circadian system, no study has experimentally established that the natural GMF variations can act as a reliable secondary zeitgeber.

In the context of our result of a threshold-dependent impact of man-made ELF MFs on MLT levels it is relevant to discuss the functional window discussed by [START_REF] Wiltschko | Sensing magnetic directions in birds: Radical pair processes involving cryptochrome[END_REF] in the case of the avian magnetic compass. The functional window at ∼ 50µT has been shown to be adaptable to variations in the static field. Indeed, [START_REF] Wiltschko | Sensing magnetic directions in birds: Radical pair processes involving cryptochrome[END_REF] and collaborators have shown that, after a few hours, migratory birds regain their magnetic sense at other intensities both low [e.g. Winklhofer et al., 2013, as low as 4µT] and high [Wiltschko et al., 2006, up to 92µT]. Note the coupling between such a weak field and biological organisms [e.g. [START_REF] Ritz | A model for photoreceptor-based magnetoreception in birds[END_REF][START_REF] Vanderstraeten | Theoretical Evaluation of Magnetoreception of Power-Frequency Fields[END_REF][START_REF] Vanderstraeten | Low-Light Dependence of the Magnetic Field Effect on Cryptochromes: Possible Relevance to Plant Ecology[END_REF][START_REF] Hore | The Radical-Pair Mechanism of Magnetoreception[END_REF][START_REF] Kattnig | The sensitivity of a radical pair compass magnetoreceptor can be significantly amplified by radical scavengers[END_REF] is far more complex than having an 'internal compass' in their beak and appears to involve chemical reactions on spin-correlated radical pairs, even though little is understood on the downstream signalling cascade mechanism(s) [as reviewed [START_REF] Nordmann | Magnetoreception -A sense without a receptor[END_REF].

Possible Limitations

Our study did not consider other possible parameters that may influence MLT excretion levels due to the lack of consistency in the parameters reported in MLT studies. In light of the main mechanism of interaction between MF and biological systems (discussed in § 3.6), such parameter might include (1) the MF polarization, (2) the amount of light and more importantly, whether or not the spectrum includes blue photons as magneto-reception appears to be blue-light dependent [e.g. [START_REF] Chasmore | Cryptochromes: Blue Light Receptors for Plants and Animals[END_REF][START_REF] Ritz | A model for photoreceptor-based magnetoreception in birds[END_REF][START_REF] Gegear | Cryptochrome mediates light-dependent magnetosensitivity in Drosophila[END_REF][START_REF] Chaves | The Cryptochromes: Blue Light Photoreceptors in Plants and Animals[END_REF][START_REF] Michael | Invited Review Animal Cryptochromes: Divergent Roles in Light Perception, Circadian Timekeeping and Beyond †[END_REF][START_REF] Vanderstraeten | Low-Light Dependence of the Magnetic Field Effect on Cryptochromes: Possible Relevance to Plant Ecology[END_REF], (3) the intensity of blue-light [as magneto-reception might be inversely proportional to the photon flux, e.g. [START_REF] Vanderstraeten | Low-Light Dependence of the Magnetic Field Effect on Cryptochromes: Possible Relevance to Plant Ecology[END_REF], (4) the time of exposure with respect to MLT rise [as suggested by [START_REF] Wood | Changes in human plasma melatonin profiles in response to 50 Hz magnetic field exposure[END_REF]Vanderstraeten and Burda, 2012], (5) the MF orientation with respect to the geomagnetic field since the radical pair (RP) mechanism involved in CRY might depend on the direction of the field line [START_REF] Wiltschko | Sensing magnetic directions in birds: Radical pair processes involving cryptochrome[END_REF][START_REF] Zhang | The Radical Pair Mechanism and the Avian Chemical Compass: Quantum Coherence and Entanglement[END_REF], ( 6) the (blue) light polarization [as discussed in [START_REF] Stoneham | A New Type of Radical-Pair-Based Model for Magnetoreception[END_REF][START_REF] Hore | The Radical-Pair Mechanism of Magnetoreception[END_REF]], (7) the possible adaptation time reported by [START_REF] Wiltschko | Sensing magnetic directions in birds: Radical pair processes involving cryptochrome[END_REF] for the avian magnetic compass, (8) the age [Vanderstraeten and Burda, 2012] and genetic factors as [START_REF] Fedrowitz | Significant Differences in the Effects of Magnetic Field Exposure on 7 , 12-Dimethylbenz ( a ) anthracene-Induced Mammary Carcinogenesis in Two Substrains of Sprague-Dawley Rats[END_REF] indicated that significant differences might occur from different substraints of rats.

Conclusions

From our analysis of 62 studies on the possible variations of MLT levels in humans and rats from [START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF]; [START_REF] Touitou | The effects of extremely low-frequency magnetic fields on melatonin and cortisol, two marker rhythms of the circadian system[END_REF]; [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF], we examined the possible relationship between a dichotomous dependent variable (corresponding to studies showing an effect or no effect on MLT excretion pattern) and independent variables such as exposure duration and magnetic field strength using a Bayesian approach and a simple logistic regression model. We find that :

• the MF exposure duration is the most significant parameter in causing changes in MLT levels both in human (Fig. 2) and rat (Fig. 4) studies, as others have reported [e.g. [START_REF] Savitz | Health Effects of Electric and Magnetic Fields: Are We Done Yet[END_REF][START_REF] Kurokawa | Acute Exposure to 50 Hz Magnetic Fields with Harmonics and Transient Components: Lack of Effects on Nighttime Hormonal Secretion in Men[END_REF][START_REF] Selmaoui | Sinusoidal 50-hz magnetic fields depress rat pineal nat activity and serum melatonin. Role of duration and intensity of exposure[END_REF][START_REF] Jahandideh | Comparing Performances of Logistic Regression and Neural Networks for Predicting Melatonin Excretion Patterns in the Rat Exposed to ELF Magnetic Fields[END_REF]Vanderstraeten et al., 2012];

• human and rat studies are entirely consistent with one another, but only after matching the MF strengths to similar ranges, i.e. B 50µT ;

• there seems to be no dose-dependence between any change in MLT levels with MF strengths ranging from 0.5 to 100µT as others have reported [e.g. [START_REF] Kato | Effects of exposure to a circularly polarized 50-Hz magnetic field on plasma and pineal melatonin levels in rats[END_REF][START_REF] Reiter | Static and Extremely Low Frequency Electromagnetic Field Exposure: Reported Effects on the Circadian Production of Melatonin[END_REF][START_REF] Pfluger | Effects of exposure to 16.7 Hz magnetic fields on urinary 6hydroxymelatonin sulfate excretion of Swiss railway workers[END_REF][START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF] 3 ;

• the impact of MF on MLT levels does, however, depend on the ELF MF strength, in the regime where ELF MFs are weaker than B t ∼ 30µT (Fig. 5). Such a window effect was already discussed in [START_REF] Löscher | Exposure of Female Rats to a 100-µT 50 Hz Magnetic Field Does Not Induce Consistent Changes in Nocturnal Levels of Melatonin[END_REF].

In light of these results, we suggest to perform additional research on rats with ELF MF with intensities in the range from 20nT to 20µT, while controlling the additional factors listed below ( § 3.7), because epidemiological studies have indicated that adverse effects on human health become noticeable at ∼0.4µT. But so far very few rat studies involved ELF MF with intensities below 5µT. This range 20nT to a few µT covers the regime experienced by humans in man-made and natural environments. Indeed, the natural variations of the geomagnetic field ranges from 20nT to a few hundred of nT [START_REF] Hitchman | The quiet daily variation in the total magnetic field: Global curves[END_REF].

Because MF strengths > 100µT are not found in nature, studies on rats with MF strengths > 100µT, or mT levels, might reveal a different (likely acute effect) than the duration-dependent effect discussed here, where perhaps one of the other factors listed below has become dominant. Table 1: Studies on the putative effect of ELF magnetic fields on MLT excretion in human subjects taken from [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF].

Notes: (1) survey identification number as in [START_REF] Halgamuge | Pineal Melatonin Levels Disruption on Human Due to Electromagnetic Fields and ICNIRP Limits[END_REF] Human studies Model A α 2.5 [0.7-4.9] (95%) β 1.3 [0.4-2.9] (95%) π 0.3 [0.0-0.5] (95%) p out 0.2 [0.0-0.9] (95%) Table 3: Results on the LR on MLT levels in humans from the logistic regression using our Bayesian analysis (see text). where the probability p for having an effect is 0.5. Right: Non-parametric Support Vector Classification (SVC) using a RBF kernel with penalty parameter C = 3.2 determined by cross-validation. In both panels, the x-values hasve been offseted by a small (random) amount to help distinguish overlapping data points.

Rat studies Model

A |B| <= 45µT |B| > 45µT α 1.2 [-0.3-4.7] (95%) 3.0 [0.5-4.9] (95%) -0.3 [-4.4-4.3] (95%) β 1.5 [-1.4-3.2] (95%) 1.2 [-0.0-2.2] (95%) 0.5 [-1.9-3.3] (95%) π 0.3 [0.0-0.5] (95%) p out 0.6 [0.0-1.0] (95%) Model B α 4.0 [1.4-7.1] (68%) β 1.1 [0.2-1.6] (68%) γ -0.1 [-2.9-2.5] (68%) log B t 1.4 [0.9-2.1] (68%) π 0.3 [0.2-0.5] (68%) p out 0.3 [0.1-0.7] (68%)

FigureFigure 2 :

 2 Figure Histogram of MF strength for studies involving human (hatched) and rats (solid). The vertical dotted line represents the geomagnetic field B at 50µT. A KS-test indicate the two histograms are not drawn from the same parent population, with a P-value of 0.01.

Figure 3 :Figure 4 :

 34 Figure 3: Bayesian Logistic Regression (model A) on rat studies with magnetic field strength B above (below) 45µT shown in the top (bottom) panel respectively

Table 2 :

 2 Studies on the putative effect of ELF magnetic fields on MLT excretion in rats.

	Notes: (1) survey identification number as in Jahandideh et al. [2010]; (2) ELF frequency (Hz); (3) EMF strength
	(µT); (4) polarization; (5) MF exposition in hours; (6) Impact on MLT excretion level; (7) Sample exposed/sham;
	(8) References.

Table 4 :

 4 Results from our Bayesian LR on MLT levels in rats (see text) with confidence intervals.

	model	t(θ)	Data-Set N WAIC LOO-CV
	(1)	(2)	(3)	(4)	(5)	(6)
	A A B	α(log T -β) α(log T -β) α, γ(log T -β) with log B ≶ B t	Humans 28 Rats 34 Rats 34	26.3 47.3 46.4	26.4 47.3 46.6

Table 5 :

 5 Information Criterion for the models considered in this paper. Notes: (1) Model number; (2) Functional form t in Eq. 2; (3) Date-set; (4) Number of studies; (5) WAIC: Widely Available Information Criteria from[START_REF] Watanabe | A Widely Applicable Bayesian Information Criterion[END_REF]; (6) LOO-CV: Leave-One-Out Cross-Validations from[START_REF] Vehtari | Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC[END_REF].

45µT is chosen to avoid the four studies which are at 50.0µT.

The threshold is somewhat loosely constrained owing to the spare sampling of MF strengths between 1 and 50 µT in published rat studies.

Note this conclusion is not in contradiction with the dose-response effect reported in epidemiological studies at the 0.1-0.5 µT levels.

Acknowledgments

This research did not receive any specific grant from funding agencies in the public, commercial or not-forprofit sectors. We are grateful to Prof. J. Vanderstraeten for very useful comments and perspectives on an earlier draft. The entire code used to generate the tables and figures in this paper is available on the CERN zenodo server at https://doi.org/10.5281/zenodo.3250993.

Software: This work made use of the following open-source software: Numpy [Van Der Walt et al., 2011], Scipy [Jones et al., 2001], Matplotlib [Hunter, 2007], Pymc3 [Salvatier et al., 2016] and Scikit-learn [Pedregosa et al., 2011].