Lionel Roques

Florian Patout

Olivier Bonnefon

Guillaume Martin

Adaptation in General Temporally Changing Environments

We analyze a nonlocal PDE model describing the dynamics of adaptation of a phenotypically structured population, under the effects of mutation and selection, in a changing environment. Previous studies have analyzed the large-time behavior of such models, with particular forms of environmental changes, either linearly changing or periodically fluctuating. We use here a completely different mathematical approach, which allows us to consider very general forms of environmental variations and to give an analytic description of the full trajectories of adaptation, including the transient phase, before a stationary behavior is reached. The main idea behind our approach is to study a bivariate distribution of two 'fitness components' which contains enough information to describe the distribution of fitness at any time. This distribution solves a degenerate parabolic equation that is dealt with by defining a multidimensional cumulant generating function associated with the distribution, and solving the associated transport equation.

We apply our results to several examples, and check their accuracy, using stochastic individual-based simulations as a benchmark. These examples illustrate the importance of being able to describe the transient dynamics of adaptation to understand the development of drug resistance in pathogens.

Introduction and main assumptions

Understanding the impact of external factors on the dynamics of fitness distributions in asexuals is a fundamental issue in population genetics, with implications for the evolution of microbial pathogens such as viruses, bacteria and cancer cells. Drug resistance may occur when a pathogenic organism (e.g., a bacteria in presence of an antibiotic) manages to reach a positive growth rate (equivalently absolute fitness) due to genetic adaptation. Being able to describe the effect of various types of environmental changes on the trajectories of adaptation is therefore a crucial issue for the elaboration of drug resistance management strategies, to which mathematical models may help answer.

Recent models of asexual adaptation based on partial differential equations (PDEs) or integro-differential equations (IDEs) typically describe the dynamics of the distribution of a single phenotypic trait in a fixed environment. This trait can be fitness itself as in [START_REF] Alfaro | Explicit solutions for replicator-mutator equations: Extinction versus acceleration[END_REF][START_REF] Gil | Mathematical properties of a class of integrodifferential models from population genetics[END_REF][START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF][START_REF] Tsimring | RNA virus evolution via a fitness-space model[END_REF], or a given trait x ∈ R determining fitness, as in [START_REF] Alfaro | Replicator-mutator equations with quadratic fitness[END_REF][START_REF] Alfaro | Evolutionary branching via replicator-mutator equations[END_REF][START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF][START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF], leading to equations of the form: ∂ t q(t, x) = M[t, x, q(t, x)] + q(t, x) (m(x) -m(t)).

Here, M is a differential or an integral operator describing the effect of mutations on the distribution q(t, •) of the trait x. The last term q(t, x) (m(x) -m(t)) corresponds to the effects of selection, see e.g. [START_REF] Tsimring | RNA virus evolution via a fitness-space model[END_REF]: m(x) is a function which describes the relationship between the trait x and fitness, and m(t) is the mean fitness in the population at time t. The fitness that we consider in this work is a 'relative fitness'. It is connected to the Malthusian growth rate r via the formula: r(x) = r max + m(x) (m(x) ≤ 0 and r max > 0 is a constant corresponding to the growth rate of an optimum phenotype).

In Fisher's geometrical model (FGM), a multivariate phenotype at a set of n traits (a vector x ∈ R n ) determines fitness. The most widely used version assumes a quadratic form of the Malthusian fitness function m(x), which decreases away from a single optimum O 0 ∈ R n , [START_REF] Martin | The fitness effect of mutations across environments: Fisher's geometrical model with multiple optima[END_REF][START_REF] Tenaillon | The utility of Fisher's geometric model in evolutionary genetics[END_REF]:

m(x) = - x -O 0 2 2 , (1)
with • the Euclidian norm in R n . To describe the mutation effects on phenotypes, the standard 'isotropic Gaussian FGM' uses a normal distribution N (0, λ I n ) with λ > 0 the phenotypic mutational variance at each trait and I n the identity matrix [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Lande | The genetic covariance between characters maintained by pleiotropic mutations[END_REF]. Overall, assuming a constant mutation rate U per capita per unit time, the corresponding integro-differential equation describing the dynamics of the phenotype distribution q(t, x), under the combined effects of selection and mutation, is ∂ t q(t, x) = U (J q -q) + q(t, x) (m(x) -m(t)), t > 0, x ∈ R n , with m(t) the mean fitness in the population at time t, and the standard convolution product in R n and J the (Gaussian) probability density function associated with the normal distribution N (0, λ I n ).

In this work, we focus on the case of a changing environment: we assume that, due to an external factor (e.g., a drug dose, a temperature, etc), the phenotype to fitness relationship (1) is changed. We take this change into account through a moving optimum, i.e. we assume that

m(t, x) = - x -O(t) 2 2 , (2)
with

O(t) = O 0 + δ(t) u, (3)
with δ(t) ∈ C(R + ), δ(0) = 0 and u a unit vector in R n (without loss of generality, we assume in the sequel that O 0 = 0, and u = (1, 0, . . . , 0)). In such case, the equation describing the dynamics of the phenotype distribution becomes:

∂ t q(t, x) = U (J q -q) + q(t, x) (m(t, x) -m(t)), t > 0, x ∈ R n ,

with this time:

m(t) = R n m(t, x) q(t, x) dx. (5)
We approach the mutational effects U (J q -q) by a diffusion (Laplace) operator, leading to the main equation that is studied in this paper:

∂ t q(t, x) = µ 2 2 ∆q + q(t, x) (m(t, x) -m(t)), t > 0, x ∈ R n , (6)
with µ = √ U λ > 0 the mutation parameter; we refer to [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF] (Appendix) for further details on the derivation of this diffusion approximation. The regime where it applies corresponds to the 'Weak Selection Strong Mutation' (WSSM) regime, where a wide diversity of lineages accumulate mutations and co-segregate at all times.

The main goal of our work is to describe the dynamics of the mean fitness m(t) in the population for very general scenarios of environmental changes, i.e., with a general form for δ(t). The value of m(t) is fundamentally connected with the question of drug resistance, or evolutionary rescue in a broader context [START_REF] Gomulkiewicz | When does evolution by natural selection prevent extinction?[END_REF], as they occur when the mean growth rate (or equivalently mean absolute fitness) in the population r(t) := r max + m(t) becomes positive.

Several particular forms of environmental changes have already been considered. The case of an optimum shifting with a constant speed δ(t) = c t has inspired several developments. First, in models without adaptation, where a favorable region moves at a constant speed, as in the Fisher-KPP reaction-diffusion equations studied in [START_REF] Berestycki | Can a species keep pace with a shifting climate?[END_REF][START_REF] Berestycki | Forced waves of the fisher-KPP equation in a shifting environment[END_REF][START_REF] Berestycki | Reaction-diffusion equations for population dynamics with forced speed I -the case of the whole space[END_REF]. Then, [START_REF] Alfaro | The effect of climate shift on a species submitted to dispersion, evolution, growth, and nonlocal competition[END_REF] considered again an optimum shifting with a constant speed, in a model including both a 1D space variable and a 1D phenotypic trait. From a mathematical viewpoint, this corresponds to an equation of the form (6) in R 2 . The case of periodically fluctuating environments has also attracted much interest from mathematicians. In the 1D case, [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF] derived explicit Gaussian solutions of a PDE model describing the dynamics of a phenotype distribution with a periodically varying phenotype optimum. In [START_REF] Carrère | Influence of mutations in phenotypically-structured populations in time periodic environment[END_REF], comparable models were considered, with phenotypes in some bounded subset of R n ; the authors give conditions for the persistence of the population, based on the sign of the principal eigenvalue of a time-periodic parabolic operator, and study the large-time behavior of the solution. In [START_REF] Figueroa | Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments[END_REF], fitness functions m(t, x), periodic with respect to t and with x ∈ R n have also been considered (see also [START_REF] Figueroa | Selection and mutation in a shifting and fluctuating environment[END_REF]). The works [START_REF] Figueroa | Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments[END_REF][START_REF] Figueroa | Selection and mutation in a shifting and fluctuating environment[END_REF] are based on the method of constrained Hamilton-Jacobi equations, which has been developed to study the evolution of phenotypically structured populations, with integral or differential mutation operators (e.g., [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF][START_REF] O Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Gandon | A Hamilton-Jacobi method to describe the evolutionary equilibria in heterogeneous environments and with non-vanishing effects of mutations[END_REF][START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF][START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF]). This method assumes a small mutation parameter of order ε 1, and is based on a scaling t → t/ε. Thus, it typically describes asymptotic evolutionary dynamics, at large times and in a 'small mutation' regime. To the best of our knowledge, it cannot lead to explicit transient trajectories of adaptation. Temporally piecewise constant environments have also been studied with the same type of methods in the recent work [START_REF] Costa | Survival criterion for a population subject to selection and mutations; application to temporally piecewise constant environments[END_REF]. Note that the equations that were studied in [START_REF] Alfaro | The effect of climate shift on a species submitted to dispersion, evolution, growth, and nonlocal competition[END_REF][START_REF] Carrère | Influence of mutations in phenotypically-structured populations in time periodic environment[END_REF][START_REF] Figueroa | Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments[END_REF][START_REF] Figueroa | Selection and mutation in a shifting and fluctuating environment[END_REF][START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF] have the general form:

∂ t n(t, x) = µ 2 2 ∆n + n(t, x) (r(t, x) -ρ(t)), t > 0, x ∈ Ω ⊆ R n , (7)
with n(t, x) the total population density and ρ(t) its integral over Ω. The study of this equation is in fact equivalent to our problem (6): it is easily checked that q(t, x) = n(t, x)/ρ(t) satisfies (6) with r(t, x) = r max + m(t, x), see Appendix A. Compared to the above-mentioned works, we use here a completely different approach, which allows us: (i) to consider very general forms of environmental variations; (ii) to give an analytic description of the full trajectories of adaptation, including the transient phase, before a stationary behavior is reached. Our results are valid in any dimension n, and do not assume that the solution has a Gaussian form. The main ideas behind our approach is to study a bivariate distribution p(t, m 1 , m 2 ) of two 'fitness components' which contains enough information to describe the distribution of fitness at any time t. The distribution p solves a degenerate parabolic equation that is dealt with by defining a multidimensional cumulant generating function associated with the distribution, and solving the associated transport equation.

Our main results are presented in the next section. We begin in Section 2.1 with a preliminary standard existence and uniqueness results of the solution q(t, x) of the Cauchy problem associated with [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF]; then, in Section 2.2 we study the distribution p(t, •) and derive the equation solved by the cumulant generating function; in Section 2.3, we present our main results on the dynamics of the mean fitness in a general setting; in Section 2.4, we apply these results to particular forms of the function δ(t), and we compare our results with the existing literature. In Section 3, we compare our theoretical results with numerical simulations of a stochastic individual-based model. These sections are followed by a discussion. Proofs are presented in Section 5.

Main results

Existence and uniqueness of the solution of the Cauchy problem

The existence and uniqueness of the solution q(t, x) of ( 6) with initial condition q 0 does not follow from standard parabolic theory as the function m(t, x) is unbounded. However, they can easily be adapted from the results in [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF], in order to take into account the time-dependence of m(t, x). We recall here the main arguments that lead to these existence and uniqueness results.

We need the following assumptions on the initial distribution q 0 :

q 0 ∈ C 2+α (R n ), (8)
for some α ∈ (0, 1), that is, q 0 C 2+α (R n ) < +∞. Moreover, as q 0 is a (probability) distribution, we assume that: q 0 ≥ 0 and

R n q 0 (x)dx = 1. (9)
We also assume that q 0 has a fast decay rate as x → +∞, in the sense that there exists a non-increasing function g ∈ C(R + , R + ) (with R + = [0, +∞)) such that:

0 ≤ q 0 ≤ g( • ) in R n , and R n e b x g( x ) dx < +∞ for all b > 0. ( 10
)
We first recall a standard existence and uniqueness result for linear parabolic equations with unbounded coefficients in R n .

Theorem 2.1 ( [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF][START_REF] Chabrowski | Sur la construction de la solution fondamentale de l'équation parabolique aux coefficients non bornés[END_REF]). The problem

   ∂ t v(t, x) = µ 2 2 ∆v + m(t, x) v(t, x), t ≥ 0, x ∈ R n , v(0, x) = q 0 (x), x ∈ R n , (11)
admits a unique positive bounded solution v ∈ C 1,2 (R + × R n ).
Moreover, it follows from the same arguments as those in lemma 4.2 of [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF] that, with v defined in Theorem 2.1,

t → m v (t) := R n m(t, x) v(t, x) dx,
is real-valued and continuous in R + and, for every t ≥ 0, there holds:

1 + t 0 m v (s) ds = 1 + t 0 R n m(t, x)v(s, x) dx ds = R n v(t, x) dx > 0.
This allows us to define:

q(t, x) = v(t, x) 1 + t 0 m v (s) ds , (12)
for every (t, x) ∈ R + × R n . Arguing as in theorem 4.1 of [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF], it is straightforward to check that q is the unique solution of (6). More precisely,

Theorem 2.2 ([24]

). There exists a unique nonnegative solution q ∈ C 1,2 (R + × R n ) of (6) such that q ∈ L ∞ ((0, T ) × R n ) for all T > 0, and the function:

t → m(t) = R n m(t, x) q(t, x) dx,
is real-valued and continuous in R + . Moreover, we have:

∀ t ≥ 0, R n q(t, x) dx = 1.
Additionally, as m(t, x) ≤ 0 and m(t) is bounded, a standard comparison argument implies that: 0 ≤ q(t, x) ≤ B(t) K q 0 (x), for some positive bounded function B(t), and with K the heat kernel in dimension n:

K(x) = 1 (2 π t µ 2 ) n/2 e - x 2 2 t µ 2 .
With the assumption [START_REF] Carrère | Influence of mutations in phenotypically-structured populations in time periodic environment[END_REF], this implies that q is exponentially bounded at all times: for all t ≥ 0 for all b > 0,

R n e b x q(t, x) dx < +∞. ( 13
)

Fitness components and cumulant generating functions

In the previous work [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF], where the optimum O 0 remained constant, it was shown that the distribution of fitness, say p(t, m), satisfies a 1D degenerate parabolic PDE. Defining the cumulant generating function

C(t, z) = ln R p(t, s) e s z ds ,
associated with this distribution, an analytically tractable 1D transport equation for C(t, z) was obtained, leading to an explicit formula for m(t) = ∂ z C(t, 0) (which was consistent with the formula in [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF] in the isotropic case).

Here, due to the time-dependence of O(t), one cannot expect to obtain a single autonomous PDE for the fitness distribution p(t, m). Consider for instance two initial distributions of q(0, x) which are symmetric with respect to O 0 = 0: q 1 (0, x) = q 2 (0, -x); then the corresponding initial distribution of fitness is the same, as m(0, x) = m(0, -x): p 1 (0, m) = p 2 (0, m). However, it is natural to expect that p 1 (t, m) ≡ p 2 (t, m) for t > 0: if q 1 (0, •) is localised around the position of the optimum at later times, q 2 (0, •) is localised in the opposite direction.

Thus, instead of focusing on the fitness distribution, we define two time-independent 'components':

   m 1 (x) = u • x = x 1 , m 2 (x) = - x 2 2 , (14)
and we denote M(x) = (m 1 (x), m 2 (x)) for all x ∈ R n , see Fig. 1. We observe that, at any time t, the fitness associated with x is uniquely determined from its components by the following formula:

m(t, x) = - x -O(t) 2 2 = δ(t) m 1 (x) + m 2 (x) - δ(t) 2 2 . ( 15
)
We define p(t, m 1 , m 2 ) the (bivariate) distribution of the components (m 1 , m 2 ) at time t. More precisely, p is defined in the next theorem.

Theorem 2.3. There exists a unique nonnegative density function

p ∈ C 1 (R + , L 2 (R × R -)) that satisfies the following relationship R n q(t, x)ϕ(M(x))dx = R×R - p(t, m 1 , m 2 )ϕ(m 1 , m 2 )dm 1 dm 2 , (16)
for every test functions ϕ ∈ L 2 (R × R -) and all t ≥ 0. Using [START_REF] O Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF], we observe that the mean fitness in the population at time t is given by

m(t) = R n m(t, x) q(t, x) dx = R×R - p(t, m 1 , m 2 ) δ(t) m 1 + m 2 - δ(t) 2 2 dm 1 dm 2 . ( 17
)
Similarly, the fitness variance in the population,

V m (t) := R n m(t, x) 2 q(t, x) dx -m(t) 2 , (18)
is given by

V m (t) = δ(t) 2 V m 1 (t) + V m 2 (t) + 2δ(t)cov m (t), (19)
with, for

j = 1, 2 V m j (t) = R×R - p(t, m 1 , m 2 )m 2 j dm 1 dm 2 - R×R - p(t, m 1 , m 2 )m j dm 1 dm 2 2
, and:

cov m (t) = R×R - p(t, m 1 , m 2 )m 1 m 2 dm 1 dm 2 - R×R - p(t, m 1 , m 2 )m 1 dm 1 dm 2 R×R - p(t, m 1 , m 2 )m 1 dm 1 dm 2 .
Moreover, we can define the 'cumulant generating function' associated with p.

Theorem 2.4. The cumulant generating function of the components m 1 , m 2 :

C(t, z 1 , z 2 ) := ln R×R - p(t, m 1 , m 2 ) e m 1 z 1 +m 2 z 2 dm 1 dm 2 , ( 20
)
for all t ≥ 0, z 1 ∈ R and z 2 ∈ R + is well-defined and belongs to C 1 (R + × R × R + ). It satisfies the following equation, for t ≥ 0 and (z 1 , z 2 ) ∈ R × R + :    ∂ t C(t, z 1 , z 2 ) = a(t) • (∇C(t, z 1 , z 2 ) -∇C(t, 0, 0)) + k(z 1 , z 2 ) • ∇C(t, z 1 , z 2 ) + γ(z 1 , z 2 ), C(0, z 1 , z 2 ) = C 0 (z 1 , z 2 ), C(t, 0, 0) = 0, ( 21
)
where

a(t) = (δ(t), 1) ∈ R 2 and k(z 1 , z 2 ) = -µ 2 (z 1 z 2 , z 2 2 ), γ(z 1 , z 2 ) = µ 2 (z 2 1 /2 -n z 2 /2). (22)
The cumulant generating function leads to simple characterizations of the central moments of the distribution. For instance, the mean fitness and the variance can easily be computed from C. [START_REF] Figueroa | Selection and mutation in a shifting and fluctuating environment[END_REF], and ( 19), we get:

For j = 1, 2 ∂ j C(t, 0, 0) = R×R -m j p(t, m 1 , m 2 ) dm 1 dm 2 , ∂ jj C(t, 0, 0) = V m j (t), and ∂ 1,2 C(t, 0, 0) = cov m (t). Using
m(t) = δ(t) ∂ 1 C(t, 0, 0) + ∂ 2 C(t, 0, 0) - δ(t) 2 2 , (23)
V m (t) = δ(t) 2 ∂ 11 C(t, 0, 0) + ∂ 22 C(t, 0, 0) + 2δ(t)∂ 12 C(t, 0, 0). ( 24
)
2.3 General formulas for the mean fitness and the fitness variance

In order to solve the equation ( 21) satisfied by C, we first note that a simpler problem can be solved explicitly. Namely, we have the following proposition.

Proposition 2.5. Let β ∈ C 1 (R 3 + , R) with β(t, 0, 0) = 0 for all t ≥ 0, and Q 0 ∈ C(R 2 + , R). For t ≥ 0 and (z, z) ∈ R 2 + , the problem    ∂ t Q(t, z, z) = (∂ z Q + ∂ z Q)(t, z, z) -(∂ z Q + ∂ z Q)(t, 0, 0) + β(t, z, z), Q(0, z, z) = Q 0 (z, z), Q(t, 0, 0) = 0, (25)
admits a unique solution, which is given by the expression:

Q(t, z, z) = t 0 β(t -s, z + s, z + s) -β(t -s, s, s) ds + Q 0 (z + t, z + t) -Q 0 (t, t). (26)
Next, we look for a change of variables such that the function C in the rescaled variables solves a system of the form [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF]. This leads to our main theorem, which can be stated as follows.

Theorem 2.6. For each t ≥ 0, define

ϕ t : R 2 + → R × R + , by ϕ t (z, z) = (y 1 (t, z, z), y 2 (z)), with        y 1 (t, z, z) := z 0 δ(z + t -s) cosh(µ s) cosh(µ z) ds + (z -z) cosh(µ(z + t)) cosh(µ z) , y 2 (z) := tanh(µ z) µ .
Let β(t, z, z) := γ(ϕ t (z, z)) with γ defined in [START_REF] Gomulkiewicz | When does evolution by natural selection prevent extinction?[END_REF] and Q be defined by [START_REF] Lande | The genetic covariance between characters maintained by pleiotropic mutations[END_REF]. Then, for all t ≥ 0 and (z, z) ∈ R 2 + , the solution of (21) satisfies:

C(t, ϕ t (z, z)) = Q(t, z, z), (27)
This result leads to an explicit expression for m(t), as stated below.

Corollary 2.7. Let Q be defined as in Theorem 2.6 and R(t, z) := Q(t, z, z). The mean fitness in ( 6) is given by

m(t) = ∂ z R(t, 0) -δ(t) 2 /2, (28)
or, more explicitly,

m(t) = -µ n 2 tanh(µ t) - 1 2 (H δ (t) -δ(t)) 2 + R 0 (t), (29)
with

H δ (t) := µ t 0 δ(u) sinh(µ u) cosh(µ t) du and R 0 (t) = 1 cosh(µ t) (δ(t) -H δ (t)) ∂ 1 C 0 (ϕ 0 (t, t)) + (1 -tanh 2 (µ t)) ∂ 2 C 0 (ϕ 0 (t, t)). (30)
The term -µ (n/2) tanh(µ t) in ( 29) corresponds to the dynamics of m(t) with a steady optimum (δ ≡ 0), and is consistent with the results in [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF][START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF][START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF]. The second term -1 2 (H δ (t) -δ(t)) 2 is a sort of squared distance between the position of the optimum at time t, and a 'weighted history' of δ for u ∈ (0, t).

The dependence of the dynamics of m(t) with respect to the initial phenotype distribution clearly shows up in Corollary 2.7, through R 0 (t). Let C clonal be the solution of [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], with an initial condition C 0,clonal (z 1 , z 2 ) = 0 corresponding to a clonal population at the optimum at t = 0, i.e. with q 0 a Dirac mass at x = 0 (the cumulant generating function associated with a Dirac mass at x * is z 1 x * 1 -z 2 x * 2 /2). Even though this initial data does not satisfy the previous conditions ( 8)-( 9) on q 0 , the function C clonal can still be defined, at least locally, through the equality [START_REF] Liu | Monotonicity of the principal eigenvalue for a linear time-periodic parabolic operator[END_REF]. We denote by m clonal (t) the corresponding value of the mean fitness:

m clonal (t) = -µ n 2 tanh(µ t) - 1 2 (H δ (t) -δ(t)) 2 .
Then, the mean fitness in ( 6) is given by

m(t) = m clonal (t) + R 0 (t), (31)
with R 0 given by [START_REF] Martin | The fitness effect of mutations across environments: a survey in light of fitness landscape models[END_REF]. Another corollary of Theorem 2.6 gives a characterization of the variance of the fitness distribution at any time.

Corollary 2.8. The fitness variance V m (t) defined by ( 18) is given by

V m (t) = ∂ zz R(t, 0) + δ (t) cosh(µ t) ∂ z Q(t, 0, 0). ( 32
)
Other moments of the distribution of fitness could be computed as well, based on the result of Theorem 2.6. For instance, the third standardized moment of the distribution (skewness) is given by the formula (87) in Appendix B.

Explicit expressions for m(t) and V m (∞): some examples

In this section, we apply the general formula for m(t) derived in Corollary 2.7 to several particular forms of functions δ(t). Explicit but lengthy expressions for V m (t) can also be derived from Corollary 2.8. We only give here the formula for the asymptotic variance as t → +∞. These results are left without proof, as they are straightforward consequences of formulas ( 29) and [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF].

Optimum shifting with a constant speed. We make here the standard assumption (e.g., [START_REF] Alfaro | The effect of climate shift on a species submitted to dispersion, evolution, growth, and nonlocal competition[END_REF][START_REF] Figueroa | Selection and mutation in a shifting and fluctuating environment[END_REF]) of an optimum O(t) which moves at a constant speed. Proposition 2.9. Assume that δ(t) = c t for some c ∈ R. Then the mean fitness is given by:

m(t) = -µ n 2 tanh(µ t) - c 2 2 µ 2 tanh 2 (µ t) + R 0 (t), (33)
with

R 0 (t) = c µ tanh(µ t) cosh(µ t) ∂ 1 C 0 (ϕ 0 (t, t)) + (1 -tanh 2 (µ t)) ∂ 2 C 0 (ϕ 0 (t, t))
and

ϕ 0 (t, t) = c µ 2 1 - 1 cosh(µ t)
, tanh(µ t) µ .

Passing to the limit t → +∞ in [START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF], we observe that

m(∞) = -µ n 2 - c 2 2 µ 2 ,
is independent of the initial phenotype distribution. In the case c = 0 (steady optimum), m(∞) = -µ n/2. This quantity is the 'mutation load', i.e., the decrease in fitness (compared to the optimum 0), due to mutations. When c = 0, the additional negative term -c 2 /(2 µ 2 ) describes the 'lag load', which corresponds here to the decrease in fitness due to the shifting in the optimum. We observe that the mutation parameter has opposite effects on the mutation and lag loads: it tends to increase the mutation load and to decrease the lag load. This leads to an optimum value µ * = (2 c 2 /n) 1/3 which maximizes m(∞).

We recall that the growth rate of a population whose phenotype distribution satisfies (6) can be described by r(t) := r max + m(t), with r max > 0 a fixed constant corresponding to the growth rate of an optimum phenotype. Persistence of the population at large times is then equivalent to r(∞) := r max + m(∞) > 0 (see Appendix A). The results of Proposition 2.9 show that the critical speed c * for persistence (r(∞) > 0 if c < c * ) or extinction (r(∞) ≤ 0 if c ≥ c * ) of the population is given by:

c * = µ 2 r max -µ n. ( 34
)
Note that the condition 2 r max -µ n ≥ 0 is necessary for the survival of the population in a fixed environment.

A consequence of the results of [START_REF] Alfaro | The effect of climate shift on a species submitted to dispersion, evolution, growth, and nonlocal competition[END_REF] is that the critical speed satisfies

c * = 2 -λ ∞ µ 2 /2, with λ ∞ the principal eigenvalue of the operator ϕ → -µ 2 /2∂ xx ϕ -(r max -x 2 /2)ϕ in R, which is given here by λ ∞ = µ/2 -r max (principal eigenfunction: ϕ(x) = exp(-x 2 /(2 µ))). Our formula (34) is therefore fully consistent with the formula c * = 2 -λ ∞ µ 2 /2 in [1].
Using Corollary 2.8, we also obtain an explicit expression for the limit V m (∞) of V m (t) as t → +∞:

V m (∞) = µ 2 n 2 + c 2 µ .
Thus, the variance of the fitness distribution increases with the speed c: a higher speed leads to a flatter distribution. However, it is a nonmonotonic function of µ: contrarily to the case of a fixed environment, V m (∞) first decreases with µ, until a critical value which is reached at µ = (c 2 /n) 1/3 , and then increases with µ. Finally using the expression for the skewness (87) derived in Appendix B, we obtain that as t → +∞, the skewness converges to:

Skew m (∞) = - µ 3 n + 3 c 2 V m (∞) 3/2 .
This negative skewness implies that the distribution of fitness is asymmetrical, with a longer left tail. The skewness becomes even more negative when the speed c is increased, which therefore reinforces the asymmetry of the distribution.

Sub-and superlinear cases. We assume here that the position of the optimum O(t) is a sublinear, or superlinear function of t: δ(t) = c t α , with α > 0, α = 1. The general formulas [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF] and ( 32) can be applied to derive explicit expressions for m(t) and V m (t). As these expressions are rather complex, we only summarize some asymptotic properties below.

Proposition 2.10. Assume that δ(t) = c t α for some c ∈ R * and α > 0.

(i) If α < 1, then m(t) → -µ n/2 and V m (t) → µ 2 n/2, as t → +∞. (ii) If α > 1, then m(t) → -∞ and V m (t) → +∞, as t → +∞.
Thus, if α < 1 (sublinear case), the lag load is equal to 0: at large times, the population tends to be as well adapted as in the case of a steady optimum. On the other hand, if α > 1 (superlinear case) the lag load is infinite; this means that adaptation is not possible.

Periodically varying optimum. The case of an optimum O(t) with a periodic trajectory is particularly relevant in applications, when a population faces an external factor which is itself periodic (concentration of an antibiotic, temperature, . . . ). We consider here a particular case, to illustrate the global shape of m(t) in such situations: O(t) oscillates between the two points ±δ max u ∈ R n . Proposition 2.11. Assume that δ(t) = δ max sin(ω t) for some ω ∈ R * + and δ max > 0. The mean fitness is given by

m(t) = -µ n 2 tanh(µ t) - 1 2
δ max ω ω 2 + µ 2 2 (ω sin(ω t) + µ cos(ω t) tanh(µ t)) 2 + R 0 (t), ( 35
)
and the average value of m(t) over one period converges to

m ∞ := lim t→+∞ ω π t+π/ω t m(s) ds = -µ n 2 - δ 2 max ω 2 4ω 2 + 4µ 2 . ( 36
)
Thus, asymptotically in time, the mean fitness m(t) becomes periodic with period π/ω. Additionally, the formula (36) tells us that the average mean fitness at large times is a decreasing function of the frequency ω/π: higher frequencies tend to impede adaptation. In a rapidly oscillating environment, i.e. as ω → +∞, the average lag load converges to -δ 2 max /4, which means that the system behaves in average as if the phenotype distribution was at a distance δ max from the optimum. Conversely, in a slowly oscillating environment, i.e. as ω → 0, the average lag load is equivalent to -δ 2 max ω 2 /(4µ 2 ); in this case, the system behaves in average as in the case of a steadily moving optimum, with speed δ max ω/ √ 2 (see Proposition 2.9). Formula (36) also enables us to study the dependence of m ∞ with respect to the mutation parameter µ: it is convex until the inflexion point µ = ω/ √ 3 and then concave. At this inflexion point,

∂ µ m ∞ (µ = ω/ √ 3) = -n/2+3δ 2 max √ 3/(32 ω). Thus, if -n/2+3δ 2 max √ 3/(32 ω) ≤ 0, the average mean fitness m ∞ is a decaying function of µ; otherwise, if -n/2+3δ 2 max √
3/(32 ω) > 0, m ∞ reaches a minimum for some value of µ in (0, ω/ √ 3), and then becomes concave and reaches a maximum for some larger value of µ. This type of dependence, with the occurrence of an optimal mutation parameter (here, leading to a higher value of m ∞ ), has already been described for a fluctuating environment in [START_REF] Carrère | Influence of mutations in phenotypically-structured populations in time periodic environment[END_REF], based on numerical simulations (see their figure 2; in their case, the mean population over one period is represented).

In [START_REF] Figueroa | Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments[END_REF], the same example (with n = 1) δ(t) = δ max sin(ω t) was inspired by an experiment on the bacterial pathogen Serratia marcescens. The method in [START_REF] Figueroa | Long time evolutionary dynamics of phenotypically structured populations in time-periodic environments[END_REF] is based on large time small mutation limit, and therefore can only give an equivalent of m(t) at large times. Their work do not focus on the mean fitness, but on the mean trait x(t) and variance v 2 (t) (the general theory deals with n dimensions, but this particular example is in 1D):

x(t) := R x q(t, x) dx and v 2 (t) := R x 2 q(t, x) dx -x(t) 2 .
They show that:

x(t) ≈ ε δ max ω sin(ω t -π/2) and v 2 (t) ≈ ε √ 2.
In our framework, in the case n = 1, the mean fitness is given by:

m(t) = - 1 2 R (x -δ(t)) 2 q(t, x) dx = - 1 2 v 2 (t) + x(t) 2 + δ(t) x(t) - δ(t) 2 2 .
Thus, their results can be used to compute an approached formula for m(t) :

m(t) ≈ - µ 2 - 1 2 δ max ω 2 (ω sin(ω t) + µ cos(ω t)) 2 . ( 37
)
We note that this is fully consistent with the large time asymptotics given by our formula [START_REF] Tenaillon | The utility of Fisher's geometric model in evolutionary genetics[END_REF], in a small mutation regime (ω/(ω 2 + µ 2 ) is approached by 1/ω) and with n = 1.

Other forms of δ(t) could be considered as well, leading to more or less complex expressions for m(t). For instance, with δ(t) = δ max sin 2 (ω t), m(t) is given by [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF], with

H δ (t) = δ max 1 2 - 1 8ω 2 + 2µ 2 2 ω µ sin(2ωt) tanh(µt) + µ 2 cos(2ωt) + 4 ω 2 cosh(µ t) , (38)
and the average value of m(t) over one period converges to

m ∞ := lim t→+∞ ω π t+π/ω t m(s) ds = -µ n 2 - δ 2 max ω 2 16ω 2 + 4µ 2 . ( 39
)
Again, we observe that m ∞ is a decreasing function of the frequency ω/π. In light of formulas ( 36) and (39), a natural conjecture is that higher frequencies always tend to impede adaptation in periodically fluctuating environments. In the case of bounded domains, for equations of the form (7), theorem 1.2 in [START_REF] Carrère | Influence of mutations in phenotypically-structured populations in time periodic environment[END_REF] shows that, for r max large enough (so that persistence occurs), the average value of the total population ρ(t) over one period converges as t → +∞ towards the principal eigenvalue λ 1 (ω) of the time-periodic operator ϕ → ∂ t ϕ-µ 2 2 ∆ ϕ-r(t, x)ϕ, and ρ(t) converges to a periodic function. Using the relationship (86) (Appendix A) between ρ(t) and r(t), integrating over one period and passing to the limit t → +∞, we obtain that m ∞ = -λ 1 (ω) -r max . Then, theorem 1.1 in [START_REF] Liu | Monotonicity of the principal eigenvalue for a linear time-periodic parabolic operator[END_REF] implies that λ 1 (ω) is an increasing function of the frequency ω/π. Finally, this shows that in bounded domains, m ∞ is indeed a decreasing function of the frequency. In our case, proving this result would require further investigation of the general formula (28) of Corollary 2.7.

Shifting and periodically fluctuating optimum. Assume that δ(t) = δ 1 (t) + δ 2 (t), with δ 1 (0) = δ 2 (0) = 0. Formula (29) implies that:

m(t) = -µ n 2 tanh(µ t) - 1 2 (H δ 1 (t) -δ 1 (t) + H δ 2 (t) -δ 2 (t)) 2 + R 0 (t),
as H δ is linear with respect to δ. Thus, we get:

m(t) = -µ n 2 tanh(µ t) + L 1 (t) + L 2 (t) -(H δ 1 (t) -δ 1 (t)) (H δ 2 (t) -δ 2 (t)) + R 0 (t), with L i (t) := - 1 2 (H δ i (t) -δ i (t)) 2
the lag associated with δ i . Let us now consider the specific example of a periodically fluctuating and shifting optimum. Combing the results of Propositions 2.9 and 2.11, we obtain the following result.

Proposition 2.12. Assume that δ(t) = c t+δ max sin(ω t) for some c ∈ R, ω ∈ R * + and δ max > 0. The mean fitness is given by

m(t) = -µ n 2 tanh(µ t) - c 2 2 µ 2 tanh 2 (µ t) - 1 2
δ max ω ω 2 + µ 2 2 (ω sin(ω t) + µ cos(ω t) tanh(µ t)) 2 (40) - c µ tanh(µ t) δ max ω ω 2 + µ 2 (ω sin(ω t) + µ cos(ω t) tanh(µ t)) + R 0 (t).
The average value of m(t) over one period converges to

m ∞ := lim t→+∞ ω 2 π t+2 π/ω t m(s) ds = -µ n 2 - c 2 2 µ 2 - δ 2 max ω 2 4ω 2 + 4µ 2 .
(41)

The effect of the two simultaneous changes (constant speed shift and periodic oscillations) is therefore not additive: we observe the emergence of an additional term that changes sign. However, in average, over one period m ∞ is the sum of the mutation load and of the lag load induced by each movement independently.

With the growth rate r(t) := r max + m(t) and assuming that the population size satisfies ρ (t) = r(t) ρ(t) or ρ (t) = ρ(t)(r(t) -ρ(t)) (as in Appendix A), we can compute the critical shifting speed for persistence (by persistence, we mean that ρ(t) → 0 as t → +∞). Proposition 2.12 shows that r(t) converges to a periodic function, thus persistence occurs if and only if r ∞ := r max + m ∞ ≥ 0 (the inequality is strict if ρ (t) = ρ(t)(r(t) -ρ(t))). Applying formula (41), we obtain the following formula for the critical shifting speed:

c * = µ 2 r max -µ n - δ 2 max ω 2 2ω 2 + 2µ 2 , ( 42
)
provided that r max is large enough so that 2 r max -µ n -δ 2 max ω 2 2ω 2 +2µ 2 ≥ 0 (otherwise persistence never occurs). We observe that the critical speed decreases with the frequency and with the amplitude of the periodic fluctuations.

Numerical computations

In this section, we check the validity of our results, (i) compared to the initial integro-differential equation ( 4), and (ii) compared to stochastic individual-based simulations of a standard model of genetic adaptation. With these two comparisons, we test the accuracy of the diffusion approximation ((4) vs ( 6)) and the effect of neglecting the stochastic aspects of mutation and selection. See Supplementary Material (S1) for some details on the numerical computation of the solution of (4).

Description of a Wright-Fisher individual-based model (IBM) with moving optimum. We assume a constant population size N . Under the assumptions of the Fisher's geometrical model, each individual i = 1, . . . , N is characterized by a phenotype x i ∈ R n . Its relative Malthusian fitness at time t (exponential growth rate) is given by (2), i.e., m i = -x i -O(t) 2 /2 and its corresponding Darwinian fitness is exp(m i ) (geometric growth rate, a discrete time counterpart of the Malthusian fitness). We assume non-overlapping generations of duration δ t = 1. Each generation, selection and genetic drift are jointly simulated by the multinomial sampling of N individuals from the previous generation, each with weight given by their Darwinian fitnesses. Mutations are then simulated by randomly drawing, for each individual, a Poisson number of mutations, with rate U . We use a classic Gaussian FGM: each single mutation has a random phenotypic effect dx drawn into a multivariate Gaussian distribution: dx ∼ N (0, λI n ), where λ > 0 is the mutational variance at each trait, and I n is the identity matrix of size n × n. Multiple mutations in a single individual have additive effects on phenotype.

Parameter values. In all cases, we take n = 3 and λ = 0.005. Based on arguments in [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF][START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF], the WSSM (diffusion) approximation should apply for U U c := n 2 λ/4, see Supplementary Material (S2) for more details. We assume here that U = 10 U c (recall that µ = √ U λ). Smaller values of U are considered in Supplementary material (S2). In the individual-based model, we assume an initially clonal distribution of the phenotypes, at the optimum 0; to be consistent with this assumption, we take Q 0 ≡ 0.

Numerical results. Fig. 2a depicts the trajectory of mean fitness when the optimum moves at a constant speed. We observe a good agreement between the analytical result of Proposition 2.9, the numerical value of m(t) given by solving (4) and the mean value of m(t) averaged over 10 3 realisations of the IBM. Note that the value of the speed c was chosen here such that the lag load is equal to the mutation load.

The trajectories of mean fitness corresponding to periodically varying optimums are presented in Figs. 2b,c. Again, the theoretical formulas accurately describe the average dynamics of the IBM and of the integro-differential equation ( 4). In particular, they capture the transient dynamics of adaptation, before m(t) tends to become periodic. In Fig. 2c, the lowest value of m(t) is reached during this transient stage, which means that extinction (or evolutionary rescue) will mainly depend on the early adaptation of the population, and not on the ultimate periodic behavior. In these plots, the parameter values are chosen such that the lag load averaged over one period is approximately equal to the mutation load.

In Fig. 2d, we considered the case of a shifting and periodically fluctuating optimum, corresponding to the situation studied in Proposition 2.12. As expected, the trajectory of mean fitness is not just a combination of the trajectories of Figs. 2a,b: the extra term

- c µ tanh(µ t) δ max ω ω 2 + µ 2 (ω sin(ω t) + µ cos(ω t) tanh(µ t))
in ( 40) tends to lower some fluctuations and to increase others, leading to minimum values of m(t) much lower than expected by simply adding the trajectories in panels a,b. Lastly, we tested the accuracy of the general formula [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF] in the case of a stochastic moving optimum. We assumed here that δ(t) was an Ornstein-Uhlenbeck process:

dδ(t) = -νδ(t) dt + β dW t , (43)
with W t the Wiener process. Given a realization of this process, the formula (29) can still be used to compute the value of m(t) (though it requires a numerical evaluation of the integral in H δ ). The results are presented in Fig. 3. Again, the dynamics of the mean fitness simulated by the IBM are well-described by our theory. Note that all of the simulations were carried out based on a single realization δ(t) of the Ornstein-Uhlenbeck process. The comparison between Figs. 3 a) and b) illustrates the complex interplay between the environment and the mutation rate: the same environment leads to very different dynamics of adaptation depending on U (or equivalently µ).

Discussion

The approach developed in this paper led to an explicit characterization of the mean fitness and variance associated with the solution of ( 6) under very general assumptions on the movement of the optimum. Our results encompass in a single framework several examples that have been recently treated in the literature, and are consistent with these anterior results. As pointed out in [START_REF] Gomulkiewicz | Evolution and the duration of a doomed population[END_REF], most experimental studies are by nature restricted to finite time horizons. Contrarily to 'traveling wave', Hamilton-Jacobi, 'Gaussian solution' or spectral approaches, which deal with large time asymptotics, our framework leads to a description of the full dynamics of the mean fitness, since the initial state t = 0, which is therefore particularly relevant for the understanding of drug resistance. Besides, as illustrated in Section 3, the transient dynamics can be very different from the large time dynamics, even in a periodically fluctuating environment. As mentioned in the Introduction the most widely used version of the FGM assumes that the Malthusian fitness decreases quadratically away from the optimum [START_REF] Martin | The fitness effect of mutations across environments: Fisher's geometrical model with multiple optima[END_REF][START_REF] Tenaillon | The utility of Fisher's geometric model in evolutionary genetics[END_REF]. Several empirical tests showed that a deviation from the quadratic would actually lead to a lower fit of the data (see [START_REF] Fraïsse | The distribution of epistasis on simple fitness landscapes[END_REF] and fig. 5 in [START_REF] Martin | The fitness effect of mutations across environments: a survey in light of fitness landscape models[END_REF]). Thus, this assumption seems to be justified from a biological point of view. From a mathematical viewpoint, extending our approach to more general phenotype to fitness landscape model remains an interesting perspective. Being able to deal with anisotropic mutation effects, as in [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF] in the case of a fixed environment, could also lead to nontrivial results regarding the interaction between the direction of the environmental change and the type of anisotropy.

a) δ(t) = c t (b) δ(t) = δmax sin(ω t) (c) δ(t) = δmax sin 2 (ω t) (d) δ(t) = c t + δmax sin(ω t)

Proofs

Proof of Theorem 2.3

For any couple of functions f and g in L 2 (R × R -), we denote by

•, • R×R -the usual inner product on L 2 (R × R -): f, g R×R -= R×R - f (s 1 , s 2 ) g(s 1 , s 2 ) ds 1 ds 2 . Similarly, •, • R n is the usual inner product on L 2 (R n ).
Let q(t, x) be the solution of (6) defined in Theorem 2.2 and M(x) : R n → R × R -the vector field defined by [START_REF] Coville | On generalized principal eigenvalues of nonlocal operators with a drift[END_REF]. For the sake of simplicity, for each t ≥ 0, we denote by q t (•) the function x → q(t, x). We consider T , the linear form defined over the Hilbert space L 2 (R × R -) by: T (ϕ) := q t , ϕ • M R n .

First, we ensure that T is a well-defined continuous linear form. Let ψ ∈ L 1 (R n ). Thanks to a change of variables, the following identity holds true:

R n ψ x 1 , - x 2 2 , x 3 , . . . , x n dx = 2 E 1,2   E n-2 ψ(y 1 , y 2 , . . . , y n ) 1 -2y 2 -y 2 1 -n i=3 y 2 i dy 3 . . . dy n   dy 1 dy 2 , (44)
where:

E 1,2 = (y 1 , y 2 ) ∈ R × R -/ -2y 2 -y 2 1 ≥ 0 , and
E n-2 = (y 3 , . . . , y n ) ∈ R n-2 / n i=3 y 2 i ≤ -2y 2 -y 2 1 for (y 1 , y 2 ) ∈ E 1,2 .
We will use this formula to show that

q t (•)ϕ•M(•) is a L 1 (R n ) function and that T is a continuous linear form. Let ϕ ∈ L 2 (R × R -).
The Cauchy-Schwartz inequality implies that, for any b > 0,:

|T (ϕ)| = R n q t (x)ϕ(M(x))dx ≤ R n ϕ x 1 , - x 2 2 2 e -b x dx 1 2 R n q t (x) 2 e b x dx 1 2 . ( 45
)
Using [START_REF] Costa | Survival criterion for a population subject to selection and mutations; application to temporally piecewise constant environments[END_REF] we know that q t is exponentially bounded, and so is q 2 t . Thus the last term of the above inequality is finite. We can now deal with the other part. Applying (44) with

ψ(x) = ϕ(x 1 , -x 2 /2) 2 exp(-b x ), we get R n ϕ x 1 , - x 2 2 2 e -b x dx = 2 E 1,2 ϕ(m 1 , m 2 ) 2 e -b √ -2m 2   E n-2 1 -2m 2 -m 2 1 -n i=3 y 2 i dy 3 . . . dy n   dm 1 dm 2 . (46)
Moreover, by a polar change of coordinates, we know that, for any y ≥ 0, and k ∈ N:

k i=1 y 2 i ≤y 1 y -k i=1 y 2 i dy 1 . . . dy k = √ y 0 1 y -r 2 Vol k-1 ( √ y)dr,
where Vol k-1 (s) is the volume of the ball of radius

s in dimension k-1, that is V k (s) = C k-1 s k-1 , with C k-1 a generic constant depending only upon the dimension k. Therefore, k i=1 y 2 i ≤y 1 y -k i=1 y 2 i dy 1 . . . dy k = π 2 Vol k-1 ( √ y), = C k y k-1 2 .
Plugging this computation into (46), we have found

R n ϕ x 1 , - x 2 2 2 e -b x dx = C n-2 E 1,2 ϕ(m 1 , m 2 ) 2 e -b √ -2m 2 (-2m 2 -m 2 1 )
n-3

2 dm 1 dm 2 , ≤ C n E 1,2 ϕ(m 1 , m 2 ) 2 dm 1 dm 2 ≤ C n R×R - ϕ(m 1 , m 2 ) 2 dm 1 dm 2 ,
for some generic constant C n that depends only on the dimension n. Using this estimate together with (45), we get:

|T (ϕ)| ≤ C q ϕ L 2 (R×R -) ,
where C q depends only on the bound on q established in ( 13) and on the dimension. Therefore, T is a well-defined continuous linear form on L 2 (R × R -).

From then, the existence and uniqueness of a function

p t ∈ L 2 (R × R -) satisfying T (ϕ) = p t , ϕ R×R - for all ϕ ∈ L 2 (R × R -)
is a straightforward application of the Riesz-Frechet representation theorem, see for instance theorem 6.19 in [START_REF] Rudin | Real and complex analysis[END_REF]. Defining p(t, •, •) := p t for each t ≥ 0, the regularity p

∈ C 1 (R + , L 2 (R × R -)
) is a straightforward consequence of the regularity of q.

Proof of Theorem 2.4

Fix z 1 , z 2 ∈ R × R + and consider an increasing sequence (ϕ k ) k∈N of nonnegative functions in

C ∞ c (R × R -) satisfying: ϕ k (m 1 , m 2 ) = exp(z 1 m 1 + z 2 m 2 ), if m 2 1 + m 2 2 < k, ϕ k (m 1 , m 2 ) = 0, if m 2 1 + m 2 2 > k + 1. (47)
Applying Theorem 2.3, one gets

R n q(t, x)ϕ k (M(x))dx = R×R - p(t, m 1 , m 2 )ϕ k (m 1 , m 2 )dm 1 dm 2 ,
with M(x) = (m 1 (x), m 2 (x)) : R n → R × R -defined by [START_REF] Coville | On generalized principal eigenvalues of nonlocal operators with a drift[END_REF]. Using the monotone convergence theorem, we can pass to the limit k → +∞ on each side of the above equality. This yields:

R n q(t, x) exp z 1 u • x -z 2 x 2 2 dx = R×R - p(t, m 1 , m 2 ) exp(z 1 m 1 + z 2 m 2 )dm 1 dm 2 .
According to [START_REF] Costa | Survival criterion for a population subject to selection and mutations; application to temporally piecewise constant environments[END_REF], q is integrable against exponential functions. Therefore, the left hand side is finite (and positive, since q ≡ 0), and so is the right hand side. Thus, the quantity

C(t, z 1 , z 2 ) := ln R×R - p(t, m 1 , m 2 ) e m 1 z 1 +m 2 z 2 dm 1 dm 2 ,
is well-defined. We now detail how to derive the equation ( 21) satisfied by C. First, differentiating ( 16) with respect to time, we find that for any compactly supported test function

ϕ ∈ C ∞ c (R × R -, R): ∂ t p, ϕ R×R -= ∂ t q, ϕ • M R n .
Plugging in (6), the equation solved by q, for any t > 0, one has

∂ t p(t, •, •), ϕ(•, •) R×R -= µ 2 2 ∆q(t, •), ϕ(M(•)) R n + q(t, •)(m(t, •) -m(t)), ϕ(M(•)) R n . ( 48
)
We deal separately with each term on the right hand side of (48). First, for the Laplace operator, we use a duality argument:

µ 2 2 ∆q, ϕ • M R n = µ 2 2 q, ∆(ϕ • M) R n , since ∆ is self-adjoint in L 2 (R n ).
We can then write:

µ 2 2 ∆(ϕ • M) = (Aϕ) • M, (49)
with A such that for all (m

1 , m 2 ) in R × R -, A(ϕ)(m 1 , m 2 ) = -µ 2 m 2 ∂ 22 ϕ(m 1 , m 2 ) + µ 2 2 ∂ 11 ϕ(m 1 , m 2 ) -µ 2 m 1 ∂ 21 ϕ(m 1 , m 2 ) - µ 2 n 2 ∂ 2 ϕ(m 1 , m 2 ). ( 50
)
With p defined by Theorem 2.3, we have:

µ 2 2 ∆q, ϕ • M R n = q, (Aϕ) • M R n , = p, A(ϕ) R×R -. (51)
Going back to (48), we notice that thanks to [START_REF] O Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF],

q(t, •)m(t, •), ϕ(M(•)) R n = q(t, •) δ(t) m 1 (•) + m 2 (•) - δ(t) 2 2 , ϕ(M(•)) R n
, which can be decomposed as

q(t, •)m(t, •), ϕ(M(•)) R n = δ(t) q(t, •), m 1 (•) ϕ(M(•)) R n + q(t, •), m 2 (•) ϕ(M(•)) R n - δ(t) 2 2 q(t, •), ϕ(M(•)) R n .
We apply Theorem (2.3) three times, to get

q(t, •)m(t, •), ϕ(M(•)) R n = p(t, m 1 , m 2 ), δ(t) m 1 + m 2 - δ(t) 2 2 ϕ(m 1 , m 2 ) R×R - . (52)
To deal with the last part of (48), we first notice that:

m(t) = q(t, •), m(t, •) R n = q(t, •), δ(t) m 1 (•) + m 2 (•) - δ(t) 2 2 R n = δ(t) p(t, m 1 , m 2 ), m 1 R×R -+ p(t, m 1 , m 2 ), m 2 R×R -- δ(t) 2 2 .
Again, the above equality involves Theorem (2.3), with ϕ = Id R×R -, which is made possible by considering increasing compactly supported approximations of identity, and a uniform limit as in (47). Thus,

m(t) q(t, •), ϕ(M(•)) R n = δ(t) m1 (t) + m2 (t) - δ(t) 2 2 p(t, m 1 , m 2 ), ϕ(m 1 , m 2 )) R×R -, (53)
where mj is defined by mj (t) := p(t, m 1 , m 2 ), m j R×R -(j = 1, 2). Plugging (51), ( 52) and ( 53) into (48), we find that

∂ t p(t, m 1 , m 2 ), ϕ(m 1 , m 2 ) R×R -= p(t, m 1 , m 2 ), A(ϕ)(m 1 , m 2 ) R×R -+ p(t, m 1 , m 2 ), δ(t) (m 1 -m1 (t)) + m 2 -m2 (t) ϕ(m 1 , m 2 ) R×R - . (54)
To be able to compute an equation on the cumulant generating function C defined by [START_REF] Gil | Mathematical properties of a class of integrodifferential models from population genetics[END_REF], we first define the moment generating function, for all z 1 ∈ R and z 2 ∈ R + :

M (t, z 1 , z 2 ) := p, E z 1 ,z 2 R×R -with E z 1 ,z 2 (m 1 , m 2 ) := exp(z 1 m 1 + z 2 m 2 ).
Applying the relationship (54) to the sequence (ϕ k ) k∈N defined by (47) and passing to the limit k → +∞, we get:

∂ t p(t, •, •), E z 1 ,z 2 (•, •) R×R -= p(t, m 1 , m 2 ), A(E z 1 ,z 2 )(m 1 , m 2 ) + δ(t) (m 1 -m1 (t)) + m 2 -m2 (t) E z 1 ,z 2 (m 1 , m 2 ) R×R - . (55)
With the explicit expression of A given in (50), we can compute ∂ t M since

∂ t M (t, z 1 , z 2 ) = ∂ t p(t, •, •), E z 1 ,z 2 (•, •) R×R -.
We find,

∂ t M (t, z 1 , z 2 ) = p, -µ 2 m 2 z 2 2 E z 1 ,z 2 + µ 2 2 z 2 1 E z 1 ,z 2 -µ 2 m 1 z 2 z 1 E z 1 ,z 2 - µ 2 n 2 z 2 E z 1 ,z 2 R×R - + p, E z 1 ,z 2 m 2 -m2 + δ(t)(m 1 -m1 ) R×R - . (56)
Moreover, for j = 1, 2,

p, m j E z 1 ,z 2 R×R -= ∂ j M (t, z 1 , z 2 ) and mj (t) = p, m j R×R -= ∂ j M (t, 0, 0).
Therefore, (56) can be rewritten as

∂ t M (t, z 1 , z 2 ) =δ(t) ∂ 1 M (t, z 1 , z 2 ) -M (t, z 1 , z 2 )∂ 1 M (t, 0, 0) + ∂ 2 M (t, z 1 , z 2 ) -M (t, z 1 , z 2 )∂ 2 M (t, 0, 0) -µ 2 z 1 z 2 ∂ 1 M (t, z 1 , z 2 ) -µ 2 z 2 2 ∂ 2 M (t, z 1 , z 2 ) + µ 2 (z 2 1 /2 -n z 2 /2)M (t, z 1 , z 2 ).
Dividing this expression by M (t, z 1 , z 2 ), and since C(t, z 1 , z 2 ) = ln(M (t, z 1 , z 2 )), this shows that C satisfies the equation in Theorem 2.4.

Proof of Theorem 2.6

We begin with the proof of Proposition 2.5.

Proof. Let T > 0. For all t ∈ [0, T ] and z, z ∈ [t -T, +∞) 2 we define

W (t, z, z) = Q(t, z + T -t, z + T -t) -Q(t, T -t, T -t).
First, we observe that:

∂ t W (t, z, z) = (∂ t Q -∂ z Q -∂ z Q)(t, z + T -t, z + T -t) -(∂ t Q -∂ z Q -∂ z Q)(t, T -t, T -t). (57)
Using [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF], we then obtain:

∂ t W (t, z, z) = β(t, z + T -t, z + T -t) -β(t, T -t, T -t).
Integrating between 0 and t, and noting that W (0, z, z) = Q 0 (z + T, z + T ) -Q 0 (T, T ), we get: Computing this quantity at (z, z) = (t -T, t -T ) and using Q(t, 0, 0) = 0, we get: Conversely, it is straightforward to check that this expression solves [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF], and the proposition follows.

W (t,
-Q(t
In order to solve our main equation [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], we look for a function ϕ t (z, z) = (y 1 (t, z, z), y 2 (z)), with y 1 ∈ C 1 (R 3 + , R) and y 2 ∈ C 1 (R + , R + ), such that the function Q(t, z, z) = C(t, ϕ t (z, z)) (61) satisfies a problem of the form [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF]. In that respect, we first establish some conditions on the functions y 1 , y 2 .

Lemma 5.1. Assume that C is a solution of ( 21), and assume that y 1 ∈ C 1 (R 3 + , R) and y 2 ∈ C 1 (R + , R + ), satisfy, for (t, z, z) ∈ R 3 + :

   ∂ t y 1 -∂ z y 1 -∂ z y 1 = -δ(t) + µ 2 y 1 y 2 , ∂ z y 2 = 1 -µ 2 y 2 2
, y 1 (t, 0, 0) = 0, ∂ t y 1 (t, 0, 0) = 0, y 2 (0) = 0, (62) then the function Q(t, z, z) defined by (61) satisfies [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF], with β(t, z, z) = γ(y 1 (t, z, z), y 2 (z)) and Q 0 (z, z) = C 0 (y 1 (0, z, z), y 2 (z)).

Proof of Lemma 5.1. Let Q be defined by (61). Then, one can note that:    ∂ t Q(t, z, z) = ∂ t C(t, ϕ t (z, z)) + ∂ t y 1 (t, z, z)∂ 1 C(t, ϕ t (z, z)), ∂ z Q(t, z, z) = ∂ z y 1 (t, z, z)∂ 1 C(t, ϕ t (z, z)) + ∂ z y 2 (z)∂ 2 C(t, ϕ t (z, z)), ∂ z Q(t, z, z) = ∂ z y 1 (t, z, z)∂ 1 C(t, ϕ t (z, z)).

(63) Thus, Q satisfies [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF], with β(t, z, z) = γ(ϕ t (z, z)) if and only if ∂ t C(t, ϕ t (z, z)) = (∂ z y 1 + ∂ z y 1 -∂ t y 1 )∂ 1 C(t, ϕ t (z, z)) + (∂ z y 2 ) ∂ 2 C(t, ϕ t (z, z))

-(∂ z y 1 + ∂ z y 1 )(t, 0, 0) ∂ 1 C(t, ϕ t (0, 0)) -∂ z y 2 (0) ∂ 2 C(t, ϕ t (0, 0)) +γ(ϕ t (z, z)).

(64)

Using (62), we just have to check that the coefficients in front of the differential terms in (64) correspond to those in [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], computed at (t, ϕ t (z, z)) to conclude the proof. Note that, at z = 0, ∂ t y 1 -∂ z y 1 -∂ z y 1 = -δ(t), and ∂ t y 1 (t, 0, 0) = 0 thus implies that (∂ z y 1 + ∂ z y 1 )(t, 0, 0) = δ(t).

Our goal is now to find some functions y 1 , y 2 , satisfying the conditions of Lemma 5.1. In order to solve the system (62), we first note that y 2 (z) = tanh(µ z)/µ (65)

satisfies the second equation in the system and the condition y 2 (0) = 0. Then, fix T > 0 and define, for t ∈ [0, T ) and z, z ∈ [t -T, +∞)

which can be solved explicitly. Namely, for any function B in C 1 (R 2 ), a solution is given by: h(t, z, z) = F (t, z + T -t) B(z + T, z) -t 0 δ(s)F (-s, z + T ) ds , with F (t, z) = cosh(µ(z + t)) cosh(µ z) .

This leads to the following expression for y 1 (t, z, z) = h(t, z + t -T, z -z): y 1 (t, z, z) = F (t, z) B(z + t, z -z) -t 0 δ(s)F (-s, z + t) ds .

(67)

The function B must be such that y 1 (t, 0, 0) = 0. We chose B(z + t, z -z) = z -z + z+t 0 δ(s)F (-s, z + t) ds, and finally get:

y 1 (t, z, z) = z 0 δ(z + t -s) cosh(µ s) cosh(µ z) ds + (z -z) F (t, z). (68)
Finally, it is immediate to check that y 2 and y 1 respectively defined by (65) and (68) satisfy the conditions of Lemma 5.1. Now, let C ∈ C 1 (R + × R × R+, R) be defined by [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], ϕ t (z, z) := (y 1 (t, z, z), y 2 (z)) : R + → R × R + and set Q(t, z, z) = C(t, ϕ t (z, z)) for all t ≥ 0 and z, z ∈ R 2 + . Then, Lemma 5.1 implies that Q(t, z, z) satisfies [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF] with Q 0 (z, z) = C 0 (ϕ 0 (z, z)) and β(t, z, z) = γ(ϕ t (z, z)). Proposition 2.5 implies that Q(t, z, z) is given by the expression [START_REF] Lande | The genetic covariance between characters maintained by pleiotropic mutations[END_REF]. This proves the result of Theorem 2.6.

Figure 1 :

 1 Figure 1: Schematic illustration of the fitness components. Given O(t), the fitness of a phenotype x, m(t, x) = -x -O(t) 2 /2, is uniquely defined by the couple (m 1 (x), m 2 (x)).

(

Figure 2 :

 2 Figure 2: Trajectories of mean fitness. Blue curve: theoretical value of m(t) (given by Proposition 2.9 in panel a, by Proposition 2.11 in panel b, by formula (29) with H δ given by (38) in panel c and Proposition 2.12 in panel d). Black dashed curve: numerical value obtained by solving the integro-differential equation (4); red circles: mean value of the mean fitness, averaged over 10 3 replicate individual-based simulations (N = 10 4 individuals in panel a; N = 10 3 individuals in panels b,c,d); pink shading: interval between the 0.025 and 0.975 quantiles of the distribution of m(t) obtained from the individual-based model; thin horizontal red line: mutation load -µ n/2; thin horizontal black line (panels b,c): asymptotic average value of the mean fitness, m ∞ given by (36) and (39); green lines: δ(t). Parameter values: panel a: c = n µ 3 ; panel b: δ max = √ 31 λ and ω = µ π; panel c: δ max = 10 √ λ and ω = µ π; panel d: same values as in panels a,b.

  (a) U = Uc (b) U = 10 Uc

Figure 3 :

 3 Figure 3: Trajectories of mean fitness with a stochastically varying optimum. Blue curve: theoretical value of m(t) given by (29); red circles: mean value of the mean fitness, averaged over 10 3 replicate individual-based simulations with N = 10 3 individuals; pink shading: interval between the 0.025 and 0.975 quantiles of the distribution of m(t) obtained from the individual-based model; thin horizontal red line: the mutation load -µ n/2. We assumed here that the position of the optimum was given by O(t) = δ(t) u, with u = (1, 0, 0) a unit vector in R 3 . Green line: δ(t), a realization of the Ornstein-Uhlenbeck process (43), with ν = 0.01 and β = 0.1.

z, z) = t 0 β

 0 (s, z + T -s, z + T -s) -β(s, T -s, T -s) ds + W (0, z, z),which leads to:Q(t, z + T -t, z + T -t) -Q(t, T -t, T -t) = t 0 β(s, z + T -s, z + T -s) -β(s, T -s, T -s) ds + W (0, z, z). (58)

  , T -t, T -t) = -s, t -s) -β(s, T -s, T -s) ds + W (0, t -T, t -T ). (59) Combining (58) and (59), we obtain:

	t
	β(s, t
	0

Q(t, z + T -t, z + T -t) = t 0 β(s, z + T -s, z + T -s) -β(s, t -s, t -s) ds + W (0, z, z) -W (0, t -T, t -T ), (60) which implies that Q(t, z, z) = t 0 β(s, z + t -s, z + t -s) -β(s, t -s, t -s) ds + Q 0 (z + t, z + t) -Q 0 (t, t).

Acknowledgements

This work was supported by the French Agence Nationale de la Recherche (ANR-18-CE45-0019 "RESISTE"). The authors thank the reviewers for valuable comments and suggestions.

Proof of Corollary 2.7

Set R(t, z) = Q(t, z, z). We have ∂ z R(t, 0) = ∂ z Q(t, 0, 0) + ∂ z Q(t, 0, 0). Using (63), we observe that ∂ z R(t, 0) = (∂ z y 1 + ∂ z y 1 )(t, 0, 0)∂ 1 C(t, 0, 0) + ∂ z y 2 (0)∂ 2 C(t, 0, 0).

From (62), we have (∂ z y 1 +∂ z y 1 )(t, 0, 0) = δ(t) and ∂ z y 2 (0) = 1. Finally, ∂ z R(t, 0) = δ(t) ∂ 1 C(t, 0, 0)+ ∂ 2 C(t, 0, 0), and using [START_REF] Gomulkiewicz | Evolution and the duration of a doomed population[END_REF] we obtain the general formula [START_REF] Lorenzi | Dissecting the dynamics of epigenetic changes in phenotype-structured populations exposed to fluctuating environments[END_REF] for m(t) :

To derive a more explicit expression, we begin by observing that:

and

Using (62), and since y 2 (u) = tanh(µ u)/µ, we get:

and differentiating ỹ1 with respect to t and integrating by parts, we get:

Combining (72), ( 73) and (74), we get:

Next, integrating by parts and using standard trigonometric formulas, we note that:

with

Similarly,

Integrating (75) between 0 and t, and using ( 76) and (77), we get:

Using (71), we get:

and coming back to (70), this shows that:

This proves the formula [START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF]. Lastly, we note that

and from formula (74),

Finally, this yields:

with ϕ 0 (t, t) = (y 1 (0, t, t), y 2 (t)). This concludes the proof of Corollary 2.7.

Proof of Corollary 2.8

In order to simplify the computations, we introduce in this section the function ŷ1 (t, z) := y 1 (t, z, z).

Notice that:

and, finally with R defined by R(t, z)

Then, from a straightforward computation, we get:

Thanks to (84),

Now coming back to the formula for the variance we established thanks to the CGF in [START_REF] Hamel | Dynamics of adaptation in an anisotropic phenotype-fitness landscape[END_REF],

Next, from Theorem 2.6, differentiating with respect to z, we get:

Since ∂ z y 1 (t, 0, 0) = -cosh(µt), we obtain the expression of the variance given in Corollary 2.8.

Appendix A: relationship between the equations ( 6) and (7)

Consider q(t, x) the solution of ( 6) and set r(t, x) = r max + m(t, x), with r max > 0 and m(t, x) defined by [START_REF] Alfaro | Explicit solutions for replicator-mutator equations: Extinction versus acceleration[END_REF]. Define the 'total population' at time t as the solution of

with r(t) := R n r(t, x) q(t, x) dx, the mean growth rate in the population at time t. Then, setting n(t, x) := ρ(t) q(t, x) (the population density) we observe that

Thus, if m(t) has a limit m(∞) as t → +∞ (e.g., in the case of a linearly or sublinearly moving optimum, see Propositions 2.9 and 2.10 (i)) the population size ρ(t) converges to r max + m(∞).

In particular, large-time persistence is equivalent to r max + m(∞) > 0.

Conversely, consider n(t, x) a positive solution of (7) with ρ(t) = Ω n(t, x) dx > 0 and r(t, x) ≤ r max , with r max a positive constant. Define q(t, x) = n(t, x)/ρ(t). Integrating the equation [START_REF] Berestycki | Can a species keep pace with a shifting climate?[END_REF] over Ω ⊆ R n , and provided that Ω ∆n = 0 (which means that the mutations do not change the total mass) we note that ρ(t) satisfies (86). Thus,

and finally, q solves an equation of the form [START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF], with m(t, x) = r(t, x) -r max .

Appendix B: Skewness

As the fitness satisfies [START_REF] O Diekmann | The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach[END_REF], m(t, x) = δ(t) m 1 (x) + m 2 (x) -δ(t) 2 /2, the third central moment M 3 (t) of the variable m(t, x) is equal to the third central moment of δ(t) m 1 (x) + m 2 (x). Let us define S(t, z) := C(t, δ(t) z, z), t ≥ 0, z ≥ 0.

The function S(t, z) satisfies:

and therefore corresponds to the cumulant generating function of δ(t) m 1 (x) + m 2 (x). Its third central moment is therefore given by: M 3 (t) = ∂ zzz S(t, 0). The skewness of the distribution of fitness is therefore given by the formula:

.

By differentiating, one finds that

with ŷ1 defined by (83). Therefore, the skewness is equal to: Next, straightforward computations yield ∂ z z Q(t, 0, 0) = ∂ 11 C(t, 0, 0) cosh(µt) 2 , ∂ z z Q(t, 0, 0) = -(δ(t) + cosh(µ t))∂ 11 C(t, 0, 0) cosh(µt) -∂ 12 C(t, 0, 0) cosh(µt) -µ sinh(µt)∂ 1 C(t, 0, 0), ∂ zzz y 2 (0) = -2µ 2 , ∂ zzz ŷ1 (t, 0) = -2µ 2 δ(t) + δ (t).

Finally, one finds that Skew m (t) = 1 V m (t) 3/2 ∂ zzz R(t, 0) + 2µ 2 ∂ z R(t, 0) + δ (t) ∂ z Q(t, 0, 0) cosh(µt) +3δ (t) ∂ z z Q(t, 0, 0) + ∂ z z Q(t, 0, 0) -µ tanh(µt)∂ z Q(t, 0, 0) cosh(µt) . (87)