F Hamel

F Lavigne

G Martin

L Roques

Dynamics of adaptation in an anisotropic phenotype-fitness landscape

We consider a diffusive 'replicator-mutator' equation to describe the adaptation of a large asexual population in a n-dimensional phenotypic space, under anisotropic mutation and selection effects. Though this equation has been studied previously when n = 1 or under isotropy assumptions, the n-dimensional anisotropic case remained unexplored. We prove here that the equation admits a unique solution, which is interpreted as the phenotype distribution, and we propose a new and general framework to the study of the quantitative behavior of this solution. When the initial distribution is Gaussian, the equation admits an explicit solution which remains normally distributed at all times, albeit with dynamic mean vector and variance-covariance matrix. In the general case, we derive a degenerate nonlocal parabolic equation satisfied by the distribution of the 'fitness components', and a nonlocal transport equation satisfied by the cumulant generating function of the joint distribution of these components. This last equation can be solved analytically and we then get a general formula for the trajectory of the mean fitness and all higher cumulants of the fitness distribution, over time. Such mean fitness trajectory is the typical outcome of empirical studies of adaptation by experimental evolution, and can thus be compared to empirical data. In sharp contrast with the known results based on isotropic models, our results show that the trajectory of mean fitness may exhibit (n-1) plateaus before it converges. It may thus appear 'non-saturating' for a transient but possibly long time, even though a phenotypic optimum exists. To illustrate the empirical relevance of these results, we show that the anisotropic model leads to a very good fit of Escherichia coli long-term evolution experiment, one of the most famous experimental dataset in experimental evolution. The two 'evolutionary epochs' that have been observed in this experiment have long puzzled the community: we propose that the pattern may simply stem form a climbing hill process, but in an anisotropic fitness landscape.

Introduction

Biological motivation

Understanding the adaptation of asexual organisms (such as viruses, bacteria or cancer cells) under the combined effects of selection and mutation is a fundamental issue in population genetics

In parallel, the development of experimental evolution (especially in microbes) has made it possible to compare observed dynamics, in the lab, with alternative models or infer various model parameters from the data (for a recent special issue on this subject see [START_REF] Rosenzweig | Experimental evolution: prospects and challenges[END_REF]). Still, the problem of modelling asexual evolutionary dynamics is inherently complex (discussed in e.g. [START_REF] Gerrish | Real time forecasting of near-future evolution[END_REF]): recurrent mutation changes the key parameters of the dynamical system constantly, competition between numerous, ever changing types must be described, and both mutational and demographic (birth/death) events are stochastic in nature.

To cope with this complexity, previous models of asexual adaptation dynamics have followed two main routes, both focusing on the dynamics of fitness, the expected reproductive output of a lineage. In a first class of models, only the expected mean fitness is followed, either directly by a deterministic ordinary differential equation (ODE) or by deriving the stationary regime of fitness change in a stochastic process of appearance and fixation of beneficial or deleterious mutations (that increase or decrease fitness, respectively). These so-called 'origin fixation models' [START_REF] Mccandlish | Modeling evolution using the probability of fixation: history and implications[END_REF] assume a limit of low diversity and large population (which approximately consists of a single genotype at all times): mean fitness follows a stochastic renewal-reward process which expectation dynamics at some stationary regime are studied. For a review of formal approaches see [START_REF] Sniegowski | Beneficial mutations and the dynamics of adaptation in asexual populations[END_REF]. Rigorous derivations of the limit behavior of the system (in regimes of low polymorphism), from a stochastic birth/death/mutation process, can be found in e.g. [START_REF] Casanova | An individual-based model for the Lenski experiment, and the deceleration of the relative fitness[END_REF] and [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF]. Some recent versions of these 'origin fixation models' have further allowed for simple forms of epistasis, namely the fact that the distribution of a mutation's effect on fitness depends on its ancestor's genetic background. In this case, a stationary regime typically does not exist but a limit ODE can still be derived for the mean fitness trajectory (e.g. [START_REF] Baake | Modelling and simulating Lenski's long-term evolution experiment[END_REF][START_REF] Kryazhimskiy | The dynamics of adaptation on correlated fitness landscapes[END_REF][START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF]).

A second class of models seeks to follow the dynamics of the full distribution of fitness within populations. These models relax the low polymorphism assumption, but in exchange for ignoring or simplifying the stochastic components of the dynamics. The resulting outcome is a partial differential equation (PDE) or an integro-differential equation (IDE). The PDE or IDE typically describes the dynamics of the distribution of a single 'trait' that may be fitness itself as in [START_REF] Tsimring | RNA virus evolution via a fitness-space model[END_REF], or a trait determining fitness (birth rate, phenotype), as in [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF]. The deterministic limit can be obtained from a microscopic stochastic model of mutation and demography (e.g. in [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF]), but the resulting equations are highly nonlinear and analytical progress requires further approximations (in a low polymorphism limit as above, using moment closure approximations). An alternative is to derive a PDE for a generating function of the fitness distribution, as in [START_REF] Gerrish | Real time forecasting of near-future evolution[END_REF][START_REF] Johnson | The approach to mutation-selection balance in an infinite asexual population, and the evolution of mutation rates[END_REF][START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF]. This PDE can also be obtained explicitly as a deterministic limit from the underlying stochastic microscopic model (e.g. in [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF]). In some cases, these PDEs for the generating functions of the fitness distribution can be solved explicitly, yielding a theory for the dynamics of the full fitness distribution under mutation and selection [START_REF] Gil | Mathematical properties of a class of integro-differential 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] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF]. In these approaches too, some forms of epistasis in the fitness effects of mutations have been considered [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].

One established find from empirical fitness trajectories is that epistasis must be accounted for to explain the data (e.g. [START_REF] Kryazhimskiy | The dynamics of adaptation on correlated fitness landscapes[END_REF]). More precisely, fitness trajectories tend to decelerate over time, implying a priori that beneficial mutations become weaker and/or rarer as the population adapts (ignoring deleterious ones, which is of course debatable). The question is then: which particular form of epistasis does explain/predict the data best? A common metaphor to explain this decelerating pattern in fitness trajectories is to invoke some form of 'fitness landscape' connecting genotypes or phenotypes with fitness, with one or several adaptive peak(s) where fitness is maximal (discussed in [START_REF] Elena | Climb every mountain?[END_REF]). In this view, deceleration in fitness trajectories stems from the hill climbing process of adaptation up a fitness peak. This view is appealing because of its intuitive/visual illustration, but also because of a particular form of landscape, Fisher's geometrical model (FGM; a single peak phenotype-fitness landscape), has shown promising potential when compared to empirical measures of fitness epistasis. For example, the FGM has shown to accurately predict observed distributions of epistasis among random and beneficial mutations in a virus and a bacterium [START_REF] Martin | Distributions of epistasis in microbes fit predictions from a fitness landscape model[END_REF], to accurately fit patterns of epistasis among beneficial mutations in a fungus [START_REF] Schoustra | Diminishing-returns epistasis among random beneficial mutations in a multicellular fungus[END_REF] or the pattern of re-adaptation (compensation) from different deleterious mutant backgrounds in a 'recovery' experiment with a bacterium [START_REF] Perfeito | Rates of fitness decline and rebound suggest pervasive epistasis[END_REF]. A review of the FGM and its empirical test can be found in [START_REF] Tenaillon | The utility of Fisher's geometric model in evolutionary genetics[END_REF].

On the other hand, this 'fitness peak' view is challenged by the observation [START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF] that, in the longest evolution experiment ever undertaken (Long Term Evolution Experiment LTEE, in the bacterium Escherishia coli ), fitness has not reached any maximum after more than 70 000 generations. It has been suggested [START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF] that this experiment actually shows a 'two epoch' dynamics with an initial (potentially saturating) fitness trajectory and a later and distinct non-saturating dynamics. A similar pattern could be invoked in another mid-term experiment with an RNA virus [START_REF] Novella | Exponential increases of RNA virus fitness during large population transmissions[END_REF]. Several other experiments did seem to yield a 'true' saturation (plateau) in fitness trajectories over several thousands of generations in E. coli (e.g. [START_REF] De Crécy-Lagard | Long term adaptation of a microbial population to a permanent metabolic constraint: overcoming thymineless death by experimental evolution of Escherichia coli[END_REF][START_REF] Fong | Metabolic gene-deletion strains of Escherichia coli evolve to computationally predicted growth phenotypes[END_REF][START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF]), but they may simply be on a too short timescale to identify subsequent fitness increases. In fact, the LTEE itself did seem to show a strongly saturating pattern after 10 000 generations [START_REF] Lenski | Dynamics of adaptation and diversification: a 10,000generation experiment with bacterial populations[END_REF].

Overall, these different insights on epistasis and adaptation trajectories are difficult to reconcile under a single theory: why would a single peak model show a good fit to mutation epistasis over single mutations, or short term fitness trajectories, yet be invalid over the longer term? Is there really a two 'epochs' dynamics in some long-term trajectories, when and why? A desirable model should reconcile both timescales, describe both epistasis among random single mutations and in long term fitness trajectories. Proposed models that do accurately fit non-saturating trajectories [START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF] have the drawback that they only focus on beneficial mutations and do not yield any prediction on the distribution of epistasis among random mutations. On the other hand, simple fitness landscapes like the FGM do predict some short term patterns and the full distribution of epistasis among random mutants, but cannot yield a never-ending fitness trajectory (a fitness plateau must be reached at the optimum).

Aim of the paper

In this article, we explore the possibility that an extended version of the FGM be able to fit multiple epochs fitness trajectories. The central extension proposed here is to consider anisotropy in the landscape, in that different underlying traits mutate and affect fitness differently. Using central limit and random matrix theory arguments [START_REF] Martin | Fisher's geometrical model emerges as a property of complex integrated phenotypic networks[END_REF], the FGM can be obtained as a limit of a much less constrained model where high-dimensional phenotypes integrate into a phenotypic network to a smaller set of traits that directly determine fitness, with an optimum. The resulting FGM, however, is not necessarily isotropic, it can also show a single dominant direction in phenotype space (that affects fitness much more) with all other directions remaining approximately equivalent. Our initial intuition is that, in these conditions, adaptation along the dominant direction will drive the early dynamics while a second 'epoch' will be visible when adaptation only proceeds along the remaining dimensions. Therefore, it seems natural to explore how such a model would fare when compared to the fitness trajectories of the LTEE. Yet, this requires deriving the fitness dynamics resulting from such a complex mutational model, in the presence of mutation and competition between asexual lineages. This is the aim of the present paper.

The existence of a phenotype optimum has been taken into account in recent PDE and IDE approaches (second class of models alluded to above). For instance, in [START_REF] Alfaro | Replicator-mutator equations with quadratic fitness[END_REF][START_REF] Alfaro | Evolutionary branching via replicator-mutator equations[END_REF], the fitness depends on a single phenotypic trait x ∈ R, through a function m(x) which admits an optimum. Extending these works to take into account the dependence of the fitness on several traits appears as a natural question; especially since we know that the number of traits affected by mutation and subject to selection critically affects the evolutionary dynamics [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF][START_REF] Waxman | Pleiotropy and the preservation of perfection[END_REF]. So far, and to the best of our knowledge, such mathematical models that take into account n-dimensional phenotypes together with the existence of a phenotype optimum always assume an isotropic dependence between the traits x ∈ R n and the fitness, and an isotropic effect of mutations on phenotypes (see [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF] for an IDE approach and [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF] for an approach based on the analysis of a PDE satisfied by a moment generating function of the fitness distribution). Our goal here is to extend these approaches in order to characterize the dynamics of adaptation in a n-dimensional phenotypic space, without making an isotropy assumption. For the sake of simplicity here, we ignore any stochastic component of the dynamics, but note that the same equations can be obtained explicitly as a deterministic limit of a stochastic model of mutation and growth [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF].

Model assumptions and definitions

We assume that a phenotype is a set of n biological traits, which is represented by a vector x ∈ R n . To describe the evolution of the phenotypic distribution we begin with an intuitive description of the concept of fitness. In an asexual population of K types of phenotypes x 1 , . . . , x K , with r max > 0 the growth rate of the optimal phenotype, the Malthusian fitness, relative to the optimal genotype is denoted m j = m (x j ) ≤ 0 for a phenotype x j ∈ R n . It is defined by the following relation: N ′ j (t) = (r max + m j ) N j (t), with N j (t) the abundance of the phenotype x j at time t. When we sum these equations over all indexes j = 1, . . . , K, we obtain an ordinary differential equation for the total population size

N (t) = ∑ K j=1 N j (t) at time t: N ′ (t) = (r max + m(t)) N (t), where the quantity m(t) = ∑ K j=1 m j N j (t) N (t)
is the mean relative fitness in the population at time t. Now if we turn to the distribution of the phenotype frequencies q(t, x j ) = N j (t) N (t), we get the partial differential equation:

∂ t q (t, x j ) = (m(x j ) -m(t)) q (t, x j ). (1)
This equation can be generalized to a continuous distribution of phenotype frequencies (see e.g. [START_REF] Tsimring | RNA virus evolution via a fitness-space model[END_REF]), as we assume in the sequel, with:

m(t) = R n m(x) q(t, x) dx. (2)
If m(x) < m(t), then (1) implies that the frequency of the phenotype x decreases, whereas if m (x) > m(t), then the frequency increases.

To define the fitness function x ↦ m(x), we use the standard Fisher's Geometric Model or FGM [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], which assumes that there is a unique optimal phenotype O ∈ R n -without loss of generality, we assume that O = (0, . . . , 0) -and that the fitness m(x) decreases quadratically with respect to each trait:

m(x) = n i=1 α i m i (x), where m i (x) = - (x i -O i ) 2 2 = - x 2 i 2 ,
and the coefficients α i are positive. 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.

In the present paper, we assume that the mutation effects can be anisotropic i.e. the mutation effects on phenotypes are described with a distribution N (0, Λ), where Λ = diag (λ 1 , . . . , λ n ) is any given positive diagonal matrix.

Assuming a small mutational variance max λ i << 1, the mutational effects can be approximated by an elliptic operator ∑ n i=1 (µ 2 i 2)∂ ii , where µ i = √ U λ i > 0 and ∂ ii denotes the second order partial derivative with respect to the i-th coordinate of x (or later m as in (15) below). We refer to Appendix A for further details.

Finally the corresponding PDE describing the dynamics of the phenotype distribution q(t, x), under the combined effects of selection and mutation, is:

∂ t q(t, x) = n i=1 µ 2 i 2 ∂ ii q(t, x) + (m(x) -m(t)) q(t, x), t > 0, x ∈ R n , (3)
with m(t) defined by [START_REF] Alfaro | Bistable travelling waves for nonlocal reaction diffusion equations[END_REF]. Recent studies [START_REF] Alfaro | Replicator-mutator equations with quadratic fitness[END_REF] have already treated the case n = 1 (see also [START_REF] Alfaro | Evolutionary branching via replicator-mutator equations[END_REF] for more general fitness functions). In this paper, we consider the general case n ≥ 1. Without loss of generality, we may assume in the sequel that α 1 = . . . = α n = 1, up to a scaling of the phenotype space (q(t, x) = (∏ n i=1 α

1 2 i ) q(t, x 1 √ α 1 , . . . , x n √ α n ))
. This leads to:

m (x) = - x 2 2 ,
with ⋅ the standard Euclidian norm in R n . Similarly, we could remove the anisotropy in the mutation effects, up to an other scaling of the phenotype space (defined this time by q(t, x) = (∏ n i=1 µ -1 i )q(t, µ 1 x 1 , . . . , µ n x n )), but in this case the coefficients α i would be replaced by α i √ µ i and hence cannot be taken all equal to 1. Note that, for the sake of simplicity, we assumed a diagonal coraviance matrix Λ and that m(x) was a linear combination of the fitness components m i (x) = -x 2 i 2. More general covariance matrices Λ and quadratic forms m(x) could have been considered as well using the transformations of the phenotype space presented in [START_REF] Martin | A general multivariate extension of Fisher's geometrical model and the distribution of mutation fitness effects across species[END_REF].

Compared to standard reaction-diffusion equations, the mathematical analysis of (3) is more involved, due to the nonlocal term m(t) (see e.g. [START_REF] Alfaro | Bistable travelling waves for nonlocal reaction diffusion equations[END_REF][START_REF] Berestycki | The non-local Fisher-KPP equation: travelling waves and steady states[END_REF][START_REF] Fang | Monotone wavefronts of the nonlocal Fisher-KPP equation[END_REF][START_REF] Faye | Modulated traveling fronts for a nonlocal Fisher-KPP equation: a dynamical systems approach[END_REF][START_REF] Genieys | Pattern and waves for a model in population dynamics with nonlocal consumption of resources[END_REF][START_REF] Gourley | Travelling front solutions of a nonlocal Fisher equation[END_REF][START_REF] Hamel | On the nonlocal Fisher-KPP equation: steady states, spreading speed and global bounds[END_REF] for other reaction-diffusion equations with nonlocal terms). In the one-dimensional case, an equation of the form (3) with a general nonlocal reaction term of the form m(x) -∫ R m(z)q(t, z)dz has been studied in [START_REF] Alfaro | Evolutionary branching via replicator-mutator equations[END_REF]. Under the assumption that m(x) tends to -∞ as x tends to ±∞, the authors have established a formula for q involving the eigenelements of the operator H = -µ 2 2 d 2 dx 2 -m(x). The formula can be made more explicit for our choice of m, m(x) = -x 2 2 [START_REF] Alfaro | Replicator-mutator equations with quadratic fitness[END_REF]. However, the method used in [START_REF] Alfaro | Replicator-mutator equations with quadratic fitness[END_REF], which consists in reducing the equation (3) to the heat equation thanks to changes of variables based on Avron-Herbst formula and generalized Lens transform, cannot be directly applied in our n-dimensional anisotropic framework. We use here a completely different approach.

In order to understand the dynamics of adaptation, an important quantity is of course the fitness distribution p(t, m), such that p(t, m) dm is the pushforward measure of the measure q(t, x) dx by the map x ↦ -x 2 2, and the mean fitness m(t). In the isotropic case [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], the authors directly focused on the equation satisfied by p(t, m) and not by q(t, x). Whereas (3) is quadratic into the variables x 1 , . . . , x n , the equation satisfied by p(t, m) is linear with respect to m, which makes possible the derivation of some PDEs satisfied by generating functions of p(t, m). Here, due to the anisotropy, the dynamics of p(t, m) is not summarized by a single one-dimensional PDE. Instead, we define the fitness components m = (m 1 , . . . , m n ) and the joint distribution of the fitness components p(t, m) such that p(t, m) dm is the pushforward measure of the measure q(t, x) dx by the map x ↦ (-x 2 1 2, ⋯, -x 2 n 2). As fitness is the sum of the fitness components, the distributions p and p are linked by:

p(t, m) = R n -∩{m≤∑ n-1 i=1 m i } p t, m 1 , . . . , m n-1 , m - n-1 i=1 m i dm 1 . . . dm n-1 , (4)
where

R -= (-∞, 0].
Thus the mean fitness m(t) can be easily connected with these distributions:

m(t) = R n m(x)q(t, x)dx = n i=1 R n - m i p(t, m)dm = R- mp(t, m)dm,
see also [START_REF] Fong | Metabolic gene-deletion strains of Escherichia coli evolve to computationally predicted growth phenotypes[END_REF] below. This paper is organized as follows. The results are presented in Section 2, and their proofs are presented in Section 4. More precisely, Section 2.1 is dedicated to the analysis of the time-dependent problem (3). We begin with the existence and uniqueness of the solution of the Cauchy problem. We then give an explicit formula for q(t, x) in the particular case of a Gaussian initial distribution of phenotypes q 0 (x). Then, we derive a nonlocal degenerate parabolic equation satisfied by p(t, m), and the equation satisfied by its cumulant generating function (CGF). Solving the equation satisfied by the CGF, we derive an explicit formula for m(t). Then, in Section 2.2, we study the long time behavior and the stationary states of (3). We shall see that the distribution of the fitness components p(t, m) converges towards a distribution p ∞ (m) as t → +∞, and we give an explicit formula for p ∞ (m). Lastly, in Section 2.3 we show that including anisotropy in the models may help to understand experimental trajectories of fitnesses, such as those obtained in the famous experiment of Richard Lenski [START_REF] Lenski | Experimental evolution in Escherichia coli. I. Adaptation and divergence during 2,000 generations[END_REF][START_REF] Lenski | Dynamics of adaptation and diversification: a 10,000generation experiment with bacterial populations[END_REF][START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF].

Main results

The time-dependent problem

Solution of the Cauchy problem associated with equation (3) for q(t, x)

We first show that the Cauchy problem admits a unique solution. We need the following assumption on the initial distribution q 0 :

q 0 ∈ C 2+α (R n ), (5)
for some α ∈ (0, 1), that is, q 0 C 2+α (R n ) < +∞. Additionally, we assume that:

R n q 0 (x)dx = 1, (6)
and there exists a function

g ∶ R + → R + (with R + = [0, +∞)) such that g is non-increasing, 0 ≤ q 0 ≤ g( ⋅ ) in R n , x ↦ m(x)g( x ) is bounded in R n and R n m(x) g( x ) dx < +∞. (7)
These assumptions are made throughout the paper, and are therefore not repeated in the statements of the results below.

We can now state an existence and uniqueness result for the distribution of phenotypes.

Theorem 2.1. There exists a unique nonnegative solution

q ∈ C 1,2 (R + × R n ) of: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t q(t, x) = n i=1 µ 2 i 2 ∂ ii q(t, x) + (m(x) -m(t)) q(t, x), t ≥ 0, x ∈ R n , q(0, x) = q 0 (x), x ∈ R n , (8)
such that q ∈ L ∞ ((0, T ) × R n ) for all T > 0, and the function:

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

∀ t ≥ 0, R n q(t, x) dx = 1.
The next result gives an explicit solution of (8) in the particular case where the phenotypes are initially Gaussian-distributed.

Corollary 2.2. Assume that the initial distribution of phenotype frequencies is Gaussian, that is,

∀ x ∈ R n , q 0 (x) = (2π) -n 2 n i=1 (s 0 i ) -1 2 exp - n i=1 (x i -q 0 i ) 2 2s 0 i , (9)
for some parameters q 0 i ∈ R and s 0 i > 0. Then the solution q(t, x) of the Cauchy problem ( 8) is Gaussian at all time:

∀ t ≥ 0, ∀ x ∈ R n , q(t, x) = (2π) -n 2 n i=1 (s i (t)) -1 2 exp - n i=1 (x i -q i (t)) 2 2s i (t) , (10)
with:

q i (t) = µ i q 0 i µ i cosh(µ i t) + s 0 i sinh(µ i t)
, and

s i (t) = µ i µ i sinh(µ i t) + s 0 i cosh(µ i t) µ i cosh(µ i t) + s 0 i sinh(µ i t) . (11)
Moreover, we have:

m(t) = - n i=1 q 2 i (t) + s i (t) 2 .
We also note, in Corollary 2.2, that the distribution q(t, x) converges to a Gaussian distribution with mean q ∞ = 0 and variances s ∞,i = µ i , as t → +∞.

The determination of an explicit formula for q(t, x) becomes more involved when the initial distribution q 0 is not a Gaussian. In this case, we study the equation satisfied by the distribution of the fitness components p(t, m) (see Introduction) and, as a by-product, we derive an explicit formula for m(t).

A degenerate parabolic PDE satisfied by p(t, m)

Our objective here is to derive an equation for p(t, m) that only involves linear dependencies with respect to the coefficients m i , and holds for general initial phenotype distributions q 0 , which may not be Gaussian.

First, we express the distribution of the fitness components p(t, m) in terms of the distribution of phenotypes q(t, x) given in Theorem 2.1.

Proposition 2.3. For all t ≥ 0 and m = (m 1 , . . . , m n ) ∈ (R * -) n , there holds:

p(t, m) = 2 -n 2 m 1 ⋯m n ε=(ε 1 ,...,εn)∈{±1} n q(t, x ε (m)), (12)
with x ε (m) = (ε 1 √ -2m 1 , . . . , ε n √ -2m n ) ∈ R n
. Furthermore, we have:

∀ t ≥ 0, R n - p(t, m) dm = 1 and m(t) = n i=1 R n - m i p(t, m) dm, (13)
where all above integrals are convergent.

It also turns out that the expression (12) becomes simpler when q satisfies some symmetry properties. In the sequel, for a given function

f ∈ C(R n ), we define its #-symmetrization f # ∈ C(R n ) by: ∀ x = (x 1 , . . . , x n ) ∈ R n , f # (x) = 2 -n ε=(ε 1 ,...,εn)∈{±1} n f (ε 1 x 1 , . . . , ε n x n ).
From the symmetries inherent to (3), it is easy to check that, if q is the solution of (3) with initial condition q 0 satisfying the conditions of Theorem 2.1, then q # is the solution of (3) with initial condition q # 0 , and q(t, ⋅) and q # (t, ⋅) have the same mean fitness m(t) at every time t ≥ 0. Furthermore, using the expression [START_REF] Faye | Modulated traveling fronts for a nonlocal Fisher-KPP equation: a dynamical systems approach[END_REF], we observe that p(t, m) can be described in terms of the #-symmetrization of q(t, ⋅):

p(t, m) = 2 n 2 m 1 ⋯m n q # (t, x 1 1 (m)), (14)
for every t ≥ 0 and m ∈ (R * -) n , with:

x 1 1 (m) = ( √ -2m 1 , . . . , √ -2m n ) ∈ R n + .
This function p(t, m) satisfies a nonlocal degenerate parabolic equation, as the following result shows:

Theorem 2.4. The distribution function of the fitness components p is a classical C 1,2 (R + × (R * -) n ) solution of: ∂ t p(t, m) = n i=1 µ 2 i m i ∂ ii p(t, m) - 3 2 n i=1 µ 2 i ∂ i p(t, m) + n i=1 m i -m(t) p(t, m), (15)
for t ≥ 0 and m ∈ (R * -) n , with initial condition:

p 0 (m) = 2 n 2 m 1 ⋯m n q # 0 (x 1 1 (m)). (16)
In the isotropic case (µ i = µ > 0 for all i), it is also possible to derive a scalar equation for the distribution of fitness p(t, m), defined by (4).

Theorem 2.5. (Isotropic case) If µ i = µ for all 1 ≤ i ≤ n, then the fitness distribution p(t, m) is a classical C 1,2 (R + × R * -) solution of: ∂ t p(t, m) = -µ 2 m ∂ mm p(t, m) + µ 2 n 2 -2 ∂ m p(t, m) + (m -m(t)) p(t, m), (17)
for t ≥ 0 and m < 0, with initial condition,

p 0 (m) = R n-1 -∩{m≤∑ n-1 i=1 m i } p 0 m 1 , . . . , m n-1 , m - n-1 i=1 m i dm 1 . . . dm n-1 . (18)
As expected, the equations ( 15) and ( 17) only involve linear combinations of the coefficients m i . This allows us to derive simpler equations satisfied by the generating functions of p(t, m) and p(t, m).

Generating functions

We define the moment generating functions (MGFs) M p and M p of p and p and their logarithms -the cumulant generating functions (CGFs) -C p and C p by:

M p (t, z) = R n - e z ⋅ m p(t, m) dm, M p (t, z) = 0 -∞ e zm p(t, m) dm, (19)
and:

C p (t, z) = log M p (t, z), C p (t, z) = log M p (t, z), (20)
for t ≥ 0, z ∈ R n + and z ≥ 0. The integrals are well defined because the components m i are all negative. Furthermore, it follows from ( 13) and the nonnegativity of p that, for each t ≥ 0, the functions M p (t, ⋅) and

C p (t, ⋅) are of class C ∞ ((R * + ) n ) ∩ C 1 (R n + ), while the functions M p (t, ⋅) and C p (t, ⋅) are of class C ∞ (R * + ) ∩ C 1 (R + ).
The following result gives the equation satisfied by C p .

Theorem 2.6. The cumulant generating function C p of p is of class C 1,1 (R + × R n + ) 1 and it solves:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂ t C p (t, z) = A(z) ⋅ ∇C p (t, z) -b(z) -m(t), t ≥ 0, z ∈ R n + , C p (0, z) = C p 0 (z), z ∈ R n + , (21)
1. This means that the partial derivatives of C p with respect to the variables t and z exist and are continuous in

R + × R n + .
where:

A(z) = (1 -µ 2 1 z 2 1 , . . . , 1 -µ 2 n z 2 n ) ∈ R n , b(z) = n i=1 µ 2 i 2 z i ∈ R, m(t) = n i=1 ∂ i C p (t, O), (22)
and ∇C p (t, z) denotes the gradient of C p with respect to the variable z.

As a corollary of this theorem, we derive an equation satisfied by C p in the isotropic case.

Corollary 2.7. (Isotropic Case) If µ i = µ for all 1 ≤ i ≤ n, then the cumulant generating function C p of the fitness distribution p is a C 1,1 (R + × R + ) solution of: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂ t C p (t, z) = (1 -µ 2 z 2 )∂ z C p (t, z) - n 2 µ 2 z -m(t), t ≥ 0, z ∈ R + , C p (0, z) = C p 0 (z), z ∈ R + , (23)
where p 0 is the initial fitness distribution given in [START_REF] Gil | Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF], and

m(t) = ∂ z C p (t, 0).
Note that the equation ( 23) is directly obtained as a by-product of ( 21), without using the equation ( 17) satisfied by p(t, m). Each of these two equations ( 21) and ( 23) has a unique solution which can be computed explicitly, leading to a formula for m(t). We refer to Section 2.3 for an application of this result.

In the Introduction, we mentioned that the coefficients µ 2 i can be interpreted as the product between the mutation rate U and the variance λ i at the i-th trait. In the isotropic case, µ 2 i = µ 2 = U λ. Thus we have retrieved the equation mentioned in [32, Eq. (E5) in Appendix E].

The last two results of this section provide some explicit expressions of C p (t, z) and C p (t, z) when z and z are close enough to O and 0.

Proposition 2.8. The cumulant generating function C p of p is given by, for all t ≥ 0 and z

= (z 1 , . . . , z n ) ∈ [0, 1 µ 1 ) × ⋯ × [0, 1 µ n ), C p (t, z) = 1 2 n i=1 log cosh(µ i t) cosh(atanh(µ i z i )) cosh(µ i t + atanh(µ i z i )) + C p (0, ψ(t, z)) -C p (0, ψ(t, O)), (24)
with:

ψ(t, z) = (ψ 1 (t, z), . . . , ψ n (t, z)) and ψ j (t, z) = 1 µ j tanh (µ j t + atanh(µ j z j )) . (25)
Moreover, for all t ≥ 0, we have:

m(t) = n i=1 (1 -tanh 2 (µ i t)) ∂ i C p (0, ψ(t, O)) - µ i 2 tanh(µ i t) . ( 26
)
Corollary 2.9. (Isotropic case) If µ i = µ for all 1 ≤ i ≤ n, then the cumulant generating function C p of p is given by:

C p (t, z) = n 2 log cosh(µ t) cosh(atanh(µ z)) cosh(µ t + atanh(µ z)) + C p (0, ϕ(t, z)) -C p (0, ϕ(t, 0)), (27)
for t ≥ 0 and 0 ≤ z < 1 µ, with ϕ(t, z) = (1 µ) tanh µ t + atanh(µ z) . Moreover, we have:

m(t) = (1 -tanh 2 (µt)) ∂ z C p 0, tanh(µt) µ - nµ 2 tanh(µt). ( 28
)

Long time behavior and stationary states

We are here interested in the long time behavior of the solutions of ( 3) and [START_REF] Genieys | Pattern and waves for a model in population dynamics with nonlocal consumption of resources[END_REF]. We begin with a result on the convergence of the solution of (15) at t → +∞. Theorem 2.10. Let p and m(t) be as in the previous section. Then: (i) p(t, ⋅) weakly converges in (R * -) n to p ∞ as t → +∞, where:

p ∞ (m) = 1 π n 2 √ µ 1 ⋯µ n m 1 ⋯m n exp n i=1 m i µ i for all m ∈ (R * -) n , (29)
in the sense that

(R * -) n p(t, m) φ(m) dm → (R * -) n p ∞ (m) φ(m) dm as t → +∞ for every test function φ ∈ C ∞ c ((R * -) n ); (ii) m(t) → m ∞ = - n i=1 µ i 2 as t → +∞ and m ∞ = n i=1 R n - m i p ∞ (m)dm; (iii) the function p ∞ is a classical C 2 ((R * -) n ) solution of: 0 = n i=1 µ 2 i m i ∂ ii p ∞ (m) - 3 2 n i=1 µ 2 i ∂ i p ∞ (m) + n i=1 m i -m ∞ p ∞ (m). (30)
In the isotropic case, we retrieve the result of [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF] in the WSSM case (Weak Selection and Strong Mutation), which says that the fitnesses are asymptotically distributed according to the symetrized Gamma distribution -Γ(n 2, µ), with µ = √ U λ:

Corollary 2.11. (Isotropic case) If µ i = µ for all 1 ≤ i ≤ n, then p(t, ⋅) weakly converges in R * -to p ∞ as t → +∞, where:

p ∞ (m) = m n 2 -1 Γ(n 2) µ n 2 exp m µ for all m < 0,
and Γ(x) = +∞ 0 t x-1 e -t dt is the standard Gamma function.

Thanks to the previous two results, we get the asymptotic behavior of the phenotype distribution in the symmetric case.

Corollary 2.12. If q 0 is #-symmetric in the sense that q 0 = q # 0 , then q(t, ⋅) = q # (t, ⋅) weakly converges to q ∞ as t → +∞ on R n where, for all x ∈ R n ,

q ∞ (x) = 1 (2π) n 2 √ µ 1 ⋯µ n exp - n i=1 x 2 i 2µ i , (31)
and

m(t) → m ∞ = R n m(x) q ∞ (x) dx as t → +∞.
We note that q ∞ is a classical positive stationary state of (3), i.e., it satisfies the following equation:

n i=1 µ 2 i 2 ∂ ii q ∞ (x) + m(x) q ∞ (x) = m ∞ q ∞ (x) for all x ∈ R n .
Thus, as already observed in [START_REF] Alfaro | Evolutionary branching via replicator-mutator equations[END_REF] in the 1D case, the distribution q(t, x) and the mean fitness m(t) converge to the principal eigenfunction (resp. eigenvalue) of the operator

H = ∑ n i=1 µ 2 i 2 ∂ ii + m(x)
. By analogy with their results, we expect that the convergence of q(t, ⋅) also occurs in L p (R n ), for all 1 ≤ p ≤ +∞, and that the convergence result remains true when the initial condition q 0 is not #-symmetric.

Effect of anisotropy: numerical computations and connection with Escherichia coli long-term evolution experiment

The objective of this section is to illustrate the importance of taking anisotropy into account when modelling adaptation trajectories. Isotropic models [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF] lead to regularly saturating trajectories of m(t) with a plateau, i.e. a single 'epoch'. Here, we show that, in the presence of anisotropy, the trajectory of m(t) can exhibit several plateaus before reaching a stable level close to m ∞ . Thus, the dynamics of adaptation can show several evolutionary 'epochs', as those observed in the E. coli long-term evolution experiment [START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF], corresponding to different time-scales at which adaptation occurs.

For the sake of simplicity of the computations, and although the existence and uniqueness results of Section 2.1 were only obtained with continuous initial distributions of phenotypes, we assume a Dirac initial distribution of the phenotypes. Namely, we assume that q 0 = δ x 0 with x 0 = (x 0,1 , . . . , x 0,n ) ∈ R n . This corresponds to an initially clonal population. In this case, the expression (26) in Proposition 2.8 simplifies to:

m(t) = n i=1 x 2 0,i 2 tanh 2 (µ i t) -1 - µ i 2 tanh (µ i t) . ( 32
)
Trajectory of adaptation in the presence of anisotropy: an illustrative example

We take n = 3 and µ 1 > µ 2 > µ 3 . The corresponding trajectory of m(t) is depicted in Figure 1. After a brief initial decay which was already observed in the isotropic case [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF], m(t) rapidly increases and reaches a first plateau (of value close to m ∞,1 ∶= -x 2 0,2 2 -x 2 0,3 2 -µ 1 2). Then, m(t) rapidly increases again to reach a second plateau (of value close to m ∞,2 ∶= -x 2 0,3 2 -µ 1 2 -µ 2 2). Finally, m(t) increases again and stabilises around m ∞ = -µ 1 2 -µ 2 2 -µ 3 2. Interestingly, although the ultimate value m ∞ does not depend on the initial phenotype, the intermediate plateaus depend on x 0 . Their values approximately correspond to the fitnesses associated with the successive projections of x 0 on the hyperplanes {x 1 = 0} and {x 1 = x 2 = 0} minus the mutation load (we recall that the optimal phenotype was fixed at x = (0, . . . , 0), see Introduction).

More generally speaking, for n ≥ 2 and µ 1 > . . . > µ n , the trajectory exhibits (n -1) The function m(t) is given by formula [START_REF] Martin | The non-stationary dynamics of fitness distributions: Asexual model with epistasis and standing variation[END_REF], and the approximations of the values of the intermediate plateaus are

m(t) m ∞ m ∞,2 m ∞,1
m ∞,1 = -x 2 0,2 2 -x 2 0,3 2 -µ 1 2 and m ∞,2 = -x 2 0,3 2 -µ 1 2 -µ 2 2. The mutational parameters are µ 1 = √ 2, µ 2 = √ 2 ⋅ 10 -2 and µ 3 = √ 2 ⋅ 10 -4 .
The other parameter values are x 0 = (3 2, 1, 1). plateaus (before the final one of value m ∞ ), of respective values:

m ≈ m ∞,k ∶= - n i=k+1 x 2 0,i 2 - k i=1 µ i 2 ,
for k = 1, . . . , (n -1). Note that these values may not be ordered by increasing order, depending on the parameter values, possibly leading to nonmonotone trajectories of m(t).

The plateaus are more visible when the µ i 's have different orders of magnitude. More precisely, we show in Section 4.5 that, given T > 0, m(t) remains around each plateau of value m ∞,k at least during a time interval of duration T , for a good choice of the parameters µ i .

Long term evolution experiment with Escherichia coli

The long term evolution experiment (LTEE) has been carried by Lenski and his collaborators since 1988 [START_REF] Lenski | Experimental evolution in Escherichia coli. I. Adaptation and divergence during 2,000 generations[END_REF]. Twelve populations of E. coli have been founded from a single common ancestor, and are still evolving after more than 70,000 generations. The fitness evolved rapidly during the first 2,000 generations, and then remained nearly static between generations 5, 000 and 10, 000 [START_REF] Lenski | Dynamics of adaptation and diversification: a 10,000generation experiment with bacterial populations[END_REF], which would at least phenomenologically advocate for the existence of a phenotype optimum. However, more recent data (after generation 10,000) indicate that the mean fitness seems to increase without bounds [START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF]. Our goal here is not to propose a new explanation of the LTEE data, but simply to check whether the anisotropic model (3) leads to a better fit than an isotropic model.

The interpretation of the fitness data from the LTEE is quite subtle (see the comments in Richard Lenski's Experimental Evolution website http://myxo.css.msu.edu/ ecoli/srvsrf.html). For the sake of simplicity, as our objective is to check if the trajectories given by ( 32) can be qualitatively consistent with the fitness data from the LTEE, we make the following simplifying assumptions: (1) a time unit t in our model corresponds to 1 cycle of the experiment (1 day), corresponding to ≈ 6.64 generations [START_REF] Lenski | Experimental evolution in Escherichia coli. I. Adaptation and divergence during 2,000 generations[END_REF]; [START_REF] Alfaro | Bistable travelling waves for nonlocal reaction diffusion equations[END_REF] the link between m(t) and the mean fitness values given by the LTEE data, which are of the form w(gen) (gen is the generation number and w the Darwinian mean fitness, which is related to the Malthusian fitness through an exponential) is made by assuming that m(t+2)-m(2) = ln [w(6.64t)]. Thus, we assume (arbitrarily) that 2 cycles were necessary to build the founding colony from the single ancestor (hence the term t + 2). Additionally, in the data, the fitness is a relative fitness against the ancestor, which implies that w(0) = 1; this is why the quantity m(2) was subtracted to m(t + 2). As mentioned above, the data are available for 12 populations. Here, we only use the data from one population (the population Ara-1, see [START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF]), for which measurements were obtained at 100-generation intervals during the first 2,000 generations, and then at 500-generation intervals.

We carried out a fit (using Matlab R ○ Curve Fitting Toolbox R ○ , with a nonlinear least squares method) with the function

f (t) = m(t + 2) -m(2).
For the sake of simplicity, we assumed a two-dimensional phenotypic space (n = 2). The only parameters to estimate are µ 1 , µ 2 and x 0 . We compared the result of this fit with a fit of the isotropic model (µ 1 = µ 2 ). The results are depicted in Figure 2. A graphical analysis shows that the anisotropic model gives a far better fit. This is confirmed by the adjusted R 2 : 0.89 for the anisotropic model versus 0.57 for the isotropic model (R 2 =1 indicates that the fitted model explains all of the variability in the data). In the anisotropic model, the fitted mutational parameters have two different orders of magnitude: µ 1 = 1.3 ⋅ 10 -2 and µ 2 = 3.0 ⋅ 10 -4 . This leads to a first plateau until about 1,000 cycles (6, 640 generations) followed by a second increase of the fitness. As expected, the isotropic model cannot explain this type of pattern.

Discussion

We considered a natural n-dimensional extension of the standard diffusive 'replicatormutator' equation, to describe the dynamics of a phenotype distribution under anisotropic mutation and selection effects, and in the presence of a phenotype optimum. We proved that the corresponding Cauchy problem was well-posed (existence and uniqueness of the solution) and we proposed a new and general framework to the study of the quantitative behavior of the solution q(t, x). This framework enabled us to derive a formula for the mean fitness in the population at all time (equation ( 26)), an important quantity to describe the dynamics of adaptation. The case of an initially Gaussian distribution q 0 of the phenotypes is simpler, as the phenotype distribution q(t, x) remains Gaussian (though anisotropic) at all time. However, when q 0 is not Gaussian, this result obviously breaks down. In this case, the method that we proposed consists in two steps: (1) to derive a degenerate nonlocal parabolic equation satisfied by the distribution p(t, m) of the fitness components, i.e., the joint distribution of the m i = -x 2 i 2. This equation is degenerate at 0, but has the advantage of being linear with respect to m; (2) to derive a nonlocal transport equation satisfied by the cumulant generating function of p(t, m). This last equation can be solved analytically, and its solution leads to explicit formulae for all the moments of p.

Conversely, the methods that are developed in this paper could be applied to solve more general degenerate parabolic PDEs. The idea would be firstly to transform the degenerate equation into a non-degenerate equation, through a change of function of the type ( 14) (rewriting q # in terms of p) to get an existence result, and secondly to consider PDEs satisfied by moment generating functions in order to obtain uniqueness and more quantitative properties of the solution of the degenerate PDE.

A natural idea to solve (3) could be to consider directly the cumulant generating function associated with q:

C q (t, z) = ln R n e z⋅x q(t, x) dx . (33)
One would then have to see for which z this quantity makes sense, since now one integrates with respect to x ∈ R n . Due to the nonlinear term m(x) in (3), the equation satisfied by C q would then be a nonlocal second-order viscous Hamilton-Jacobi type equation:

∂ t C q (t, z) = - 1 2 ∆C q (t, z) - 1 2 ∇C q (t, z) 2 + n i=1 µ 2 i 2 z 2 i -m(t),
whereas equation [START_REF] Hamel | On the nonlocal Fisher-KPP equation: steady states, spreading speed and global bounds[END_REF] for C p was a nonlocal first-order transport equation. Finally, the results of Section 2.3 illustrate the importance of taking anisotropy into account from an applied perspective, as it can open up further explanations of experimental data. In sharp contrast with the known results based on isotropic models, our results show that the trajectory of adaptation may exhibit (n -1) plateaus before it reaches the final one. In particular, using our analytic formulae for the dynamics of m(t), we obtained a very good fit of one of the most famous experimental dataset in experimental evolution, for which several evolutionary 'epochs' had already been observed [START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF]. This suggests that the FGM, in its anisotropic form, can generate a diversity of behaviours that may reconcile the various outcomes observed in long vs. short term experimental evolution. However, whether this versatility may lead to overparameterization remains an open statistical question.

Proofs

This section is devoted to the proofs of the results announced in Section 2. Section 4.1 is concerned with problem (8) satisfied by the distribution of phenotype frequencies q(t, x), while Section 4.2 deals with the fitness frequencies p(t, m) and the proofs of Proposition 2.3 and Theorem 2.4. In Section 4.3, we carry out the proofs of Theorem 2.6 and Proposition 2.8 on the cumulant generating functions, and their corollaries (Theorem 2.5 and Corollaries 2.7 and 2.9) in the isotropic case. Lastly, Sections 4.4 and 4.5 are concerned with the stationary states and the existence of plateaus for the mean fitness.

Proofs of Theorem 2.1 and Corollary 2.2 on the Cauchy problem (8)

Before considering the nonlocal problem (8), we begin with the local problem:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ t v(t, x) = n i=1 µ 2 i 2 ∂ ii v(t, x) + m(x) v(t, x), t ≥ 0, x ∈ R n , v(0, x) = q 0 (x), x ∈ R n , (34)
where q 0 satisfies ( 5)- [START_REF] Casanova | An individual-based model for the Lenski experiment, and the deceleration of the relative fitness[END_REF]. As the fitness function:

m(x) = - x 2 2 ,
is unbounded, standard parabolic theory with bounded coefficients does not apply. However, some properties of (34) can be for instance be found in [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF][START_REF] Lunardi | Schauder theorems for linear elliptic and parabolic problems with unbounded coefficients in R n[END_REF][START_REF] Protter | Maximum Principles in Differential Equations[END_REF], which in particular lead to the following result.

Theorem 4.1. The problem (34) admits a unique bounded solution

v ∈ C 1,2 (R + × R n ).
Additionally, we have:

∀T > 0, ∃ S > 0, ∀ t ∈ [0, T ], v(t, ⋅) C 2+α (R n ) ≤ S q 0 C 2+α (R n ) , (35)
and:

0 < v(t, x) < (K t * q 0 )(x) for all (t, x) ∈ R * + × R n , (36)
with:

K t (x) = 1 (2πt) n 2 µ 1 ⋯µ n exp - n i=1 x 2 i 2µ 2 i t . ( 37
)
Proof. Let us fix a time T > 0. Theorem 2 in [START_REF] Lunardi | Schauder theorems for linear elliptic and parabolic problems with unbounded coefficients in R n[END_REF] implies that (34) admits a unique bounded solution v ∈ C 1,2 ([0, T ] × R n ) and this solution satisfies [START_REF] Perfeito | Rates of fitness decline and rebound suggest pervasive epistasis[END_REF]. Theorem III in [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF] further implies that this solution is nonnegative. As T was chosen arbitrarily, these existence, uniqueness and nonnegativity results extend to t ∈ (0, +∞), with local boundedness in t.

Let us set h(t, x) ∶= (K t * q 0 )(x) for t > 0, and h(0, x) = q 0 (x). The function h satisfies:

∂ t h(t, x) = n i=1 µ 2 i 2 ∂ ii h(t, x),
for all t > 0 and x ∈ R n . Let ψ(t, x) ∶= v(t, x) -h(t, x). We see that, for all t > 0 and

x ∈ R n , ∂ t ψ(t, x) - n i=1 µ 2 i 2 ∂ ii ψ(t, x) = m(x)v(t, x) ≤ 0, ( 38
)
and ψ(0, x) = 0. By the Phragmèn-Linderlöf principle [36, Theorem 10, Chapter 3], we get that

ψ ≤ 0 in R + × R n , i.e., v(t, x) ≤ h(t, x) = (K t * q 0 )(x). We therefore infer that 0 ≤ v(t, x) ≤ (K t * q 0 )(x) in (0, +∞) × R n . Since q 0 is bounded, this also implies that v is bounded in R + × R n .
By the standard strong parabolic maximum principle, we conclude that the first inequality is strict, i.e., 0 < v(t, x), for all (t, x) ∈ (0, +∞) × R n , since v(0, ⋅) = q 0 is continuous, nonnegative and not identically equal to 0. Furthermore, since the inequality in ( 38) is then strict for all (t, x) ∈ (0, +∞) × (R n ∖ {O}), we get that ψ(t, x) < 0, i.e., v(t, x) < (K t * q 0 )(x), for all (t, x) ∈ (0, +∞) × R n .

In order to connect [START_REF] Novella | Exponential increases of RNA virus fitness during large population transmissions[END_REF] and ( 8), we need the following lemma.

Lemma 4.2. The function:

t ↦ m v (t) ∶= R n m(x) v(t, x) dx, (39)
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(x)v(s, x) dx ds = R n v(t, x) dx > 0. ( 40
)
Proof. First of all, denote:

µ ∶= min(µ 1 , . . . , µ n ) > 0 and µ ∶= max(µ 1 , . . . , µ n ) > 0. ( 41
)
It follows from the assumptions on q 0 that m v (0) is a nonpositive real number. Consider now any t > 0 and let us check that m v (t) defined in ( 39) is a nonpositive real number. The function x ↦ m(x) v(t, x) is nonpositive and continuous in R n . Furthermore, it follows from Theorem 4.1 and the assumptions on q 0 that:

R n m(x) v(t, x) dx ≤ 1 2(2πt) n 2 µ n R n R n exp - x -y 2 2µ 2 t g( y ) x 2 dy dx = µ n 2π n 2 µ n R n R n exp(-z 2 ) g( x -µ √ 2tz ) x 2 dz dx ≤ µ n 2π n 2 µ n R n z ≤ x (2µ √ 2t) exp(-z 2 ) g( x 2) dz + z > x (2µ √ 2t) exp(-z 2 ) g(0) dz x 2 dx ≤ µ n 2π n 2 µ n π n 2 R n g( x 2) x 2 dx + g(0) R n h t ( x ) x 2 dx , where h t (r) ∶= z ≥r (2µ √ 2t)
exp(-z 2 ) dz = O(exp(-r)), as r → +∞.

Since the integral ∫ R n x 2 g( x ) dx converges by assumption, one concludes that the integral ∫ R n m(x) v(t, x) dx converges as well, hence m v (t) is a nonpositive real number. Furthermore, since the quantities h t (r) are non-decreasing with respect to t > 0, the same arguments together with Lebesgue's dominated convergence theorem imply straightforwardly that the function m v is continuous in R + .

The convergence of the integral defining m v (t) for every t ≥ 0, together with the nonnegativity and continuity of v and x ↦ -m(x) = x 2 2, immediately implies that the integral ∫ R n v(t, x) dx is a nonnegative real number for each t ≥ 0. Moreover, since v(t, x) > 0 for all (t, x) ∈ R * + × R n , we infer that:

R n v(t, x) dx > 0 for all t > 0.
As in the previous paragraph, the function t ↦ ∫ R n v(t, x) dx is also continuous in R + .

Fix now an arbitrary t > 0. For R > 0, denote B R = {x ∈ R n , x < R}, ν the outward unit normal on ∂B R and dσ(x) the surface measure on ∂B R . For every ε ∈ (0, t) and every R > 0, the integration of ( 34) over [ε, t] × B R yields:

B R v(t, x) dx - B R v(ε, x) dx = 1 2 t ε ∂B R n i=1 µ 2 i ν i ∂ i v(s, x) dσ(x) ds + t ε B R m(x) v(s, x) dx ds.
Since the first term of the right-hand side converges to 0 as R → +∞ from standard parabolic estimates (see [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF]), one gets that:

R n v(t, x) dx - R n v(ε, x) dx = t ε m v (s) ds,
by passing to the limit R → +∞. The limit ε → 0 + then yields:

R n v(t, x) dx -1 = R n v(t, x) dx - R n q 0 (x) dx = t 0 m v (s) ds, (42)
which gives the desired result (40) for t > 0. Formula (40) for t = 0 is trivial since v(0, ⋅) = q 0 has unit mass. The proof of Lemma 4.2 is thereby complete.

Proof of Theorem 2.1. Let v ≥ 0 be the unique classical solution of [START_REF] Novella | Exponential increases of RNA virus fitness during large population transmissions[END_REF] given in Theorem 4.1. From Lemma 4.2, the function m v defined by ( 39) is continuous in R + . Let us then set:

q(t, x) = v(t, x) 1 + ∫ t 0 m v (s) ds , (43)
for every (t, x) ∈ R + × R n . We recall that the denominator in the right-hand side of ( 43) is positive from Lemma 4.2. From Theorem 4.1 and Lemma 4.2, the function q is nonnegative and of class C 1,2 (R + × R n ). Furthermore, q(0, ⋅) = v(0, ⋅) = q 0 in R n and it follows from Lemma 4.2 that:

∀ t ≥ 0, R n q(t, x) dx = 1 and m(t) ∶= R n m(x) q(t, x) dx = m v (t) 1 + ∫ t 0 m v (s) ds ,
hence the function m is real-valued, nonpositive and continuous in R + . Lastly, it is immediate to see that q obeys (8). Furthermore, since v is bounded in R + × R n , the function q is bounded in (0, T ) × R n for every T > 0.

To show the uniqueness, assume now that we have two nonnegative classical solutions q 1 and q 2 of (8

) in C 1,2 (R + × R n ) ∩ L ∞ ((0, T ) × R n ) (for every T > 0)
, with the same initial condition q 0 satisfying ( 5)-( 7), and such that the functions:

m 1 (t) ∶= R n m(x) q 1 (t, x) dx and m 2 (t) ∶= R n m(x) q 1 (t, x) dx
are real-valued and continuous in R + . Define:

v i (t, x) = q i (t, x) exp t 0 m i (s)ds , for i = 1, 2, and (t, x) ∈ R + × R n .
The two functions v 1 and v 2 satisfy (34) and are bounded in (0, T ) × R n for every T > 0. From Theorem 4.1, we get v 1 ≡ v 2 in R + × R n . Furthermore, for all i = 1, 2 and t ∈ R + , there holds:

m v i (t) = R n m(x) v i (t, x) dx = m i (t) exp t 0 m i (s)ds = d dt exp t 0 m i (s)ds .
Hence, for all (t, x) ∈ R + × R n , we get:

q 1 (t, x) = v 1 (t, x) 1 + ∫ t 0 m v 1 (s) ds = v 2 (t, x) 1 + ∫ t 0 m v 2 (s) ds = q 2 (t, x).
The proof of Theorem 2.1 is thereby complete.

Proof of Corollary 2.2. It is a straightforward calculation to check that the function q defined by ( 10)-( 11) is a classical bounded solution of ( 8) with initial condition given by ( 9). The conclusion then follows from the uniqueness part of Theorem 2.1.

A degenerate parabolic PDE satisfied by p(t, m)

Proof of Proposition 2.3. We recall that the phenotypes are represented by n traits, and so by a vector in R n and that we have a constant optimum O, which is fixed to (0, . . . , 0) up to translation. We define, for each ε = (ε 1 , . . . , ε n ) ∈ {±1} n , the subset:

Q ε = {x ∈ ε 1 R * + × ⋯ × ε n R * + } ⊂ R n .
For any time t ≥ 0, we get from the law of total probability that:

R n - p(t, m) dm = R n q(t, x) dx = 1,
and that, for any

m ∈ R n -, p(t, m) = ε∈{±1} n p(t, m x(m) ∈ Q ε ) P(x(m) ∈ Q ε ), (44)
with p(t, m x(m) ∈ Q ε ) the conditional density of the fitness vector m, given that the associated phenotype x(m) is in Q ε , and

P(x(m) ∈ Q ε ) the probability that x(m) ∈ Q ε , i.e., P(x(m) ∈ Q ε ) = Qε q(t, y)dy.
In the above formula [START_REF] Waxman | Pleiotropy and the preservation of perfection[END_REF], we also use the fact that ∫ H q(t, x) dx = 0 with:

H = ⋃ 1≤i≤n {x ∈ R n , x i = 0}, since q(t, ⋅) is continuous in R n . As the fitness function x ∈ Q ε ↦ m(x) = (m 1 (x), . . . , m n (x)) with m i (x) = -x 2 i 2 is one-to-one from Q ε to (R * -) n , with inverse m ↦ x ε (m) = (ε 1 √ -2m 1 , . . . , ε n √ -2m n ), we infer that, for every t ≥ 0 and m ∈ (R * -) n , p(t, m x(m) ∈ Q ε ) = 1 det J ε ⋅ q(t, x ε (m)) P(x(m) ∈ Q ε ) = 2 -n 2 m 1 ⋯m n q(t, x ε (m)) P(x(m) ∈ Q ε ) , with J ε = diag(-ε 1 √ -2m 1 , . . . , -ε n √ -2m n ).
Finally, we get:

p(t, m) = 2 -n 2 m 1 ⋯m n ε∈{±1} n q(t, x ε (m)),
and:

m(t) = n i=1 R n - x 2 i 2 q(t, x)dx = n i=1 R n ∖H - x 2 i 2 q(t, x)dx = n i=1 ε∈{±1} n {ε 1 x 1 >0,...,εnxn>0} - x 2 i 2 q(t, x)dx = n i=1 ε∈{±1} n (R * -) n m i q(t, x ε (m)) 2 -n 2 m 1 ⋯m n dm = n i=1 (R * -) n m i p(t, m)dm.
Notice that all integrals in the last sum converge since all integrands are nonpositive and the sum of these integrals is a real number. Observe lastly that

(R * -) n m i p(t, m)dm = R n - m i p(t, m)dm, for every 1 ≤ i ≤ n, since p(t, ⋅) is an L 1 (R n -) function. The proof of Proposition 2.3 is thereby complete.
Proof of Theorem 2.4. Formula [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF] implies that the function p is of class C 1,2 (R + ×(R * -) n ) with initial condition p 0 given by [START_REF] Gerrish | Real time forecasting of near-future evolution[END_REF]. Furthermore, it is straightforward to check that, for all t ≥ 0 and m ∈ (R * -) n :

2 -n 2 ∂ i p(t, m) = - q # (t, x 1 1 (m)) 2m i m 1 ⋯m n - 1 2 m i ∂ i q # (t, x 1 1 (m)) m 1 ⋯m n , 2 -n 2 ∂ ii p(t, m) = 3 4 q # (t, x 1 1 (m)) m 2 i m 1 ⋯m n + 3 2
∂ i q # (t, x 1 1 (m)) m i 2 m i m 1 ⋯m n - 1 2m i ∂ ii q # (t, x 1 1 (m)) m 1 ⋯m n , with x 1 1 (m) = ( √ -2m 1 , . . . , √ -2m n ) ∈ (R * + ) n
. Hence, we have:

n i=1 µ 2 i m i ∂ ii p(t, m) + 3 2 ∂ i p(t, m) = - n i=1 2 n 2 m 1 ⋯m n µ 2 i 2 ∂ ii q # (t, x 1 1 (m)) = -2 n 2 ∂ t q # (t, x 1 1 (m)) -(m(x 1 1 (m)) -m(t)) q # (t, x 1 1 (m)) m 1 ⋯m n = -∂ t p(t, m) + n i=1 m i -m(t) p(t, m).
Theorem 2.4 is thereby proven.

Proof of Theorem 2.5. Since by definition p(t, m) dm is the pushforward measure, at each time t ≥ 0, of the measure q(t, x) dx by the map x ↦ -x 2 2, it follows from the layer-cake formula that:

p(t, m) = (2 m ) n 2-1 S n-1 q t, 2 m σ dσ = (2 m ) n 2-1 Q t, 2 m , for all (t, m) ∈ R + × R * +
, where S n-1 denotes the unit Euclidean sphere of R n and:

Q(t, r) = S n-1 q(t, r σ) dσ, for (t, r) ∈ R + × R + . Since q is of class C 1,2 (R + × R n ),
it is easy to see that the function:

(t, x) ↦ q(t, x) ∶= S n-1 q(t, x σ) dσ, is of class C 1,2 (R + × R n ) too, hence the function Q is of class C 1,2 (R + × R + ) and p is of class C 1,2 (R + × R * -). Furthermore, q(0, ⋅) is of class C 2+α (R n ) since q(0, ⋅) = q 0 is of class C 2+α (R n ).
Since q solves (8), which is invariant by rotation in the present isotropic case (µ i = µ for every 1 ≤ i ≤ n) and since m(t) = ∫ R n m(x) q(t, x) dx = ∫ R n m(x) q(t, x) dx, for every t ≥ 0, it follows that q solves (8) as well, with initial condition q(0, ⋅). As a consequence, Q satisfies:

∂ t Q(t, r) = µ 2 2 ∂ rr Q(t, r) + n -1 r ∂ r Q(t, r) + - r 2 2 -m(t) Q(t, r), for all (t, r) ∈ R + × R * + . But since: Q(t, r) = r 2-n p t, - r 2 2 ,
for all (t, r) ∈ R + × R * + , it is then straightforward to check that p solves [START_REF] Gil | Mathematical properties of a class of integro-differential models from population genetics[END_REF] in R + × R * -. Lastly, the formula ( 18) is an immediate consequence of the definitions of p and p, and the proof of Theorem 2.5 is thereby complete.

Generating functions

Proof of Theorem 2.6. Given q 0 satisfying ( 5)-( 7), we call v the unique bounded nonnegative C 1,2 (R + × R n ) solution of (34) defined in Theorem 4.1 with initial condition q # 0 . Notice that the function q # 0 satisfies the same assumptions ( 5)-( 7) as q 0 . By uniqueness and symmetry of [START_REF] Novella | Exponential increases of RNA virus fitness during large population transmissions[END_REF] with respect to the change of variable x i into -x i , for any 1 ≤ i ≤ n, it follows that v(t, x) = v # (t, x) and, as in [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF],

q # (t, x) = v # (t, x) 1 + ∫ t 0 m v (s) ds = v(t, x) 1 + ∫ t 0 m v (s) ds , (45)
for all (t, x) ∈ R + × R n + . As already noticed in Section 2.1, from [START_REF] Fong | Metabolic gene-deletion strains of Escherichia coli evolve to computationally predicted growth phenotypes[END_REF] and the nonnegativity of p, the functions (t, z) ↦ M p (t, z) and (t, z) ↦ C p (t, z) given in ( 19)-( 20) are well defined in R + × R n + . Let us start by proving the continuity of these functions M p and C p in R + × R n + . Owing to the relations ( 14) and ( 45), we see that:

M p (t, z) = 2 n 2 1 + ∫ t 0 m v (s)ds R n - e z⋅m m 1 ⋯m n v(t, x 1 1 (m)) dm, (46)
for all (t, z

) ∈ R + × R n + , with x 1 1 (m) = ( √ -2m 1 , . . . , √ -2m n ).
Hence, we have:

M p (t, z) = 2 n 1 + ∫ t 0 m v (s)ds R n + exp - n i=1 z i x 2 i 2 v(t, x) dx = 1 1 + ∫ t 0 m v (s)ds R n exp - n i=1 z i x 2 i 2 v(t, x) dx. (47)
Notice that the function t ↦ ∫ t 0 m v (s)ds is continuous in R + . Furthermore, the function

(t, z) ↦ exp(-∑ n i=1 z i x 2 i 2) v(t, x) is also continuous in R + × R n + , for every x ∈ R n .
Lastly, as in the proof of Lemma 4.2, it follows from ( 7), ( 36)-( 37) and ( 41) that, for any z ∈ R n + , T > 0, t ∈ (0, T ] and x ∈ R n , there holds:

0 ≤ exp - n i=1 z i x 2 i 2 v(t, x) ≤ v(t, x) ≤ 1 (2πt) n 2 µ n R n exp - x -y 2 2µ 2 t q # 0 (y) dy ≤ µ n π n 2 µ n R n e -y ′ 2 g( x -µ √ 2ty ′ ) dy ′ ≤ µ n π n 2 µ n y ′ ≤ x (2µ √ 2t)
e -y ′ 2 g( x 2) dy ′ +g(0)

y ′ > x (2µ √ 2t) e -y ′ 2 dy ′ ≤ µ n π n 2 µ n π n 2 g( x 2) + g(0) y ′ > x (2µ √ 2T )
e -y ′ 2 dy ′ .

(48)

Call h(x) the quantity given in the right-hand side of the last inequality. Since by [START_REF] Casanova | An individual-based model for the Lenski experiment, and the deceleration of the relative fitness[END_REF] the function

g( ⋅ ) is in L ∞ (R n ) ∩ L 1 (R n ), the function h belongs to L 1 (R n ),
∂ i M p (t, z) = R n - m i e z⋅m p(t, m) dm, (49)
hence, as in ( 46)-(47),

∂ i M p (t, z) = 2 n 2 1 + ∫ t 0 m v (s)ds R n - m i e z⋅m m 1 ⋯m n v(t, x 1 1 (m)) dm = - 1 1 + ∫ t 0 m v (s)ds R n x 2 i 2 exp - n i=1 z i x 2 i 2 v(t, x) dx (50)
from ( 14) and ( 45). On the one hand, the function t ↦ ∫ t 0 m v (s)ds is continuous in R + and so is the function

(t, z) ↦ x 2 i exp(-∑ n i=1 z i x 2 i 2) v(t, x) in R + × R n + , for every x ∈ R n .
On the other hand, as in the previous paragraph, it follows from (36)-( 37) and (41) that, for any 1 ≤ i ≤ n, z ∈ R n + , T > 0, t ∈ (0, T ] and x ∈ R n , there holds:

0 ≤ x 2 i exp - n i=1 z i x 2 i 2 v(t, x) ≤ x 2 v(t, x) ≤ µ n π n 2 µ n π n 2 x 2 g( x 2) + g(0) x 2 y ′ > x (2µ √ 2T )
e -y ′ 2 dy ′ .

Call h(x) the quantity given in the right-hand side of the last inequality. Since by [START_REF] Casanova | An individual-based model for the Lenski experiment, and the deceleration of the relative fitness[END_REF] the function x ↦ x 2 g( x ) is in L 1 (R n ), the function h belongs to L 1 (R n ), and is independent of z ∈ R n + and t ∈ (0, T ]. One then infers from (50) and Lebesgue's dominated convergence theorem that the function

∂ i M p is continuous in R + × R n
+ , and so is the function

∂ i C p = ∂ i M p M p .
In this paragraph, we are interested in the differentiation of M p with respect to t. By [START_REF] Novella | Exponential increases of RNA virus fitness during large population transmissions[END_REF], we already know that:

∂ t v(t, x) ≤ n i=1 µ 2 i 2 ∂ ii v(t, x) + m(x) v(t, x) , (51)
for all (t, x) ∈ R + × R n . Fix T > 0 and let S > 0 be the constant given in [START_REF] Perfeito | Rates of fitness decline and rebound suggest pervasive epistasis[END_REF]. Thus, for all (t, x) ∈ [0, T ] × R n , there holds:

∂ ii v(t, x) ≤ v(t, ⋅) C 2+α (R n ) ≤ S q # 0 C 2+α (R n ) ≤ S q 0 C 2+α (R n ) . (52)
Let us now focus on the boundedness of the second term of the right-hand side of (51), that is, the boundedness of the function (t, x) ↦ m(x)v(t, x) in [0, T ] × R n . Since this function is continuous in R + × R n , let us show its boundedness in (0, T ] × R n . Thanks to [START_REF] Casanova | An individual-based model for the Lenski experiment, and the deceleration of the relative fitness[END_REF] and ( 36)-( 37), we get, as in (48), that:

m(x)v(t, x) ≤ R n µ n x 2 2π n 2 µ n q # 0 x -µ √ 2ty e -y 2 dy,
for all (t, x) ∈ (0, T ] × R n . Thus, we have:

m(x)v(t, x) ≤ µ n π n 2 µ n R n x -µ √ 2ty 2 + 2tµ 2 y 2 g x -µ √ 2ty e -y 2 dy ≤ C µ n π n 2 µ n + 2µ n+2 T g(0) π n 2 µ n R n y 2 e -y 2 dy,
where the constant C is such that x ′ 2 q # 0 (x ′ ) ≤ x ′ 2 g( x ′ ) ≤ C, for all x ′ ∈ R n . Therefore, the function (t, x) ↦ m(x)v(t, x) is bounded in [0, T ] × R n for any T > 0, and so is ∂ t v by (51)-( 52). Together with [START_REF] Perfeito | Rates of fitness decline and rebound suggest pervasive epistasis[END_REF], [START_REF] Wiser | Long-term dynamics of adaptation in asexual populations[END_REF] and the continuity of m v in R + , it follows that the function ∂ t q # is bounded in [0, T ] × R n , for every T > 0. Finally, [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF] implies that for all (t, z)

∈ [0, T ] × R n + and m ∈ (R * -) n , e z⋅m ∂ t p(t, m) ≤ 2 n 2 ∂ t q # L ∞ ([0,T ]×R n ) e z⋅m m 1 ⋯m n .
Since the integrals:

R n - e z⋅m m 1 ⋯m n dm = 2 n 2 R n + exp - n i=1 z i x 2 i 2 dx,
converge for all z ∈ (R * + ) n , it then easily follows from the previous estimates and from Lebesgue's dominated convergence theorem that the function M p is differentiable with respect to t in R + × (R * + ) n , with:

∂ t M p (t, z) = R n - e z⋅m ∂ t p(t, m) dm, (53)
and that the function

∂ t M p is itself continuous in R + × (R * + ) n . So is the function ∂ t C p = ∂ t M p M p . The continuity of the functions ∂ t M p and ∂ t C p in the closure R + × R n + of R + × (R *
+ ) n will be obtained as a consequence of the equations satisfied by these two functions, which shall be established below.

Let us then turn to find an equation satisfied by M p , in order to derive the equation (21) satisfied by C p . Fix (t, z) ∈ R + × (R * + ) n . Thanks to (15), we have:

e z⋅m ∂ t p(t, m) = - n i=1 µ 2 i m i e z⋅m ∂ ii p(t, m) - 3 2 n i=1 µ 2 i e z⋅m ∂ i p(t, m) + n i=1 m i -m(t) e z⋅m p(t, m), (54)
for all m ∈ (R * -) n . We are now going to integrate (54) over R n -. To do so, let us first focus on the first two terms of the right-hand side of (54). Fix an index i ∈ 1, n and consider the cubes:

B ε = m = (m 1 , . . . , m n ) ∈ R n -, -ε -1 < m j < -ε, for each j ∈ 1, n = (-ε -1 , -ε) n , with 0 < ε < 1. Denote B = (-ε -1 , -ε) n-1 and: m = (m 1 , . . . , m i-1 , m i+1 , . . . , m n ), m ρ = (m 1 , . . . , m i-1 , -ρ, m i+1 , . . . , m n ),
for ρ ∈ R. By using Fubini's theorem and integrating by parts with respect to the variable m i , one infers that:

Bε m i e z⋅m ∂ ii p(t, m) + 3 2 e z⋅m ∂ i p(t, m) dm = Bε z i 2 + m i z 2 i e z⋅m p(t, m) dm + Bε z i ε e z⋅m ε p(t, m ε ) -z i ε -1 e z⋅m 1 ε p(t, m 1 ε ) d m + Bε e z⋅m ε p(t, m ε ) 2 -ε∂ i p(t, m ε ) -e z⋅m 1 ε p(t, m 1 ε ) 2 -ε -1 ∂ i p(t, m 1 ) d m. (55)
Let us pass to the limit as ε → 0 in the three integrals of the right-hand side of (55). Firstly, since p(t, ⋅) is nonnegative and the functions m ↦ p(t, m) and m ↦ m i p(t, m) are in L 1 (R n -), it follows from Lebesgue's dominated convergence theorem together with ( 19) and (49) that:

Bε z i 2 +m i z 2 i e z⋅m p(t, m) dm → R n - z i 2 +m i z 2 i e z⋅m p(t, m) dm = z i 2 M p (t, z)+z 2 i ∂ i M p (t, z),
as ε → 0. Secondly, by denoting:

Ĉε = √ 2ε, √ 2ε -1 n-1 , ẑ = (z 1 , . . . , z i-1 , z i+1 , . . . , z n ),
and:

x = (x 1 , . . . , x i-1 , x i+1 , . . . , x n ), x ρ = (x 1 , . . . , x i-1 , 2ρ, x i+1 , . . . , x n ),
for ρ ≥ 0, it follows from (14) that:

Bε z i ε e z⋅m ε p(t, m ε ) d m = z i √ ε e -z i ε Bε 2 n 2 e ẑ⋅ m q # (t, x 1 1 (m ε )) m 1 ⋯m i-1 m i+1 ⋯m n d m = z i √ ε e -z i ε 2 n-1 2 Ĉε exp - j≠i z j x 2 j 2 q # (t, x ε ) dx.
Since the continuous function q # (t, ⋅) is bounded in R n by ( 35) and ( 45), and since z ∈ (R * + ) n , one then gets that:

Bε z i ε e z⋅m ε p(t, m ε ) d m → 0, as ε → 0.
Similarly, we prove that:

Bε z i ε -1 e z⋅m 1 ε p(t, m 1 ε ) d m = z i √ ε -1 e -z i ε 2 n-1 2 Ĉε exp - j≠i z j x 2 j 2 q # (t, x 1 ε ) dx → 0,
as ε → 0. Thirdly, from the computations done in the proof of Theorem 2.4, we already know that, for all m ∈ (R * -) n :

1 2 p(t, m) + m i ∂ i p(t, m) = 2 n 2 ∂ i q # (t, x 1 1 (m)) 2 m 1 ⋯m i-1 m i+1 ⋯m n .
Hence, we have:

Bε e z⋅m ε p(t, m ε ) 2 -ε∂ i p(t, m ε ) d m = e -z i ε Bε 2 n 2-1 2 e ẑ⋅ m ∂ i q # (t, x 1 1 (m ε )) m 1 ⋯m i-1 m i+1 ⋯m n d m = e -z i ε 2 n-1 Ĉε exp - j≠i z j x 2 j 2 ∂ i q # (t, x ε ) dx.
Since the function ∂ i q # (t, ⋅) is continuous and bounded in R n by ( 35) and ( 45), since ∂ i q # (t, x 0 ) = 0, by #-symmetry of q # (t, ⋅), and since z ∈ (R * + ) n , one then infers from Lebesgue's dominated convergence theorem that:

Bε e z⋅m ε p(t, m ε ) 2 -ε∂ i p(t, m ε ) d m → 0, as ε → 0.
Furthermore, the integral:

Bε e z⋅m 1 ε p(t, m 1 ε ) 2 - √ ε -1 ∂ i p(t, m 1 ε ) d m = e -z i ε 2 n-1 Ĉε exp ⎛ ⎝ - j≠i z j x 2 j 2 ⎞ ⎠ ∂ i q # (t, x 1 ε ) dx,
converges to 0 as ε → 0. Coming back to (55) and passing to the limit as ε → 0, it follows from the previous estimates that:

Bε m i e z⋅m ∂ ii p(t, m) + 3 2 e z⋅m ∂ i p(t, m) dm → ε→0 z i 2 M p (t, z) + z 2 i ∂ i M p (t, z). (56)
Let us finally remember (54) and that the functions:

m ↦ e z⋅m ∂ t p(t, m), m ↦ m i e z⋅m p(t, m) and m ↦ e z⋅m p(t, m), are in L 1 (R n -) with integrals given by (53), ( 49) and [START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF], respectively. Together with (55) and (56), one concludes that:

∂ t M p (t, z) = n i=1 (1 -µ 2 i z 2 i )∂ i M p (t, z) - 1 2 n i=1 µ 2 i z i M p (t, z) -m(t)M p (t, z), (57)
for every (t, z) ∈ R + × (R * + ) n . Since the right-hand side of the above equation is continuous

in R + × R n + , one infers that the function ∂ t M p is extendable by continuity in R + × R n + Hence, m(t) → -∑ n i=1 µ i 2 =∶ m ∞ , as t → +∞. It is also straightforward to check that: n i=1 R n - m i p ∞ (m) dm = m ∞ ,
and that p ∞ is a classical C 2 ((R * -) n ) solution of (30) (this property is also a consequence of the fact that p satisfies (15) and the distributions p(t, ⋅) weakly converge to p ∞ as t → +∞). The proof of Theorem 2.10 is thereby complete.

Proof of Corollary 2.11. By the same arguments as in the proof of Theorem 2.10, thanks to [START_REF] Lunardi | Schauder theorems for linear elliptic and parabolic problems with unbounded coefficients in R n[END_REF], we see that, for all z ∈ [0, 1 µ), C p (t, z) → -(1 2) log(1 + µ z) as t → +∞. This limiting function corresponds to the cumulant generating function of a random variable distributed according to -Γ(n 2, µ). Since a distribution is uniquely determined by its cumulant generating function, this implies that p ∞ is the probability density function of this random variable, i.e., for all m < 0,

p ∞ (m) = 1 Γ(n 2) µ n 2 m n 2 -1 exp m µ ,
where Γ(x) = ∫ +∞ 0 t x-1 e -t dt is the standard Gamma function.

Proof of Corollary 2.12. We assume here that q 0 is #-symmetric, that is, q 0 = q # 0 . As already emphasized in Section 2.1, the uniqueness for problem [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF] implies that q(t, ⋅) is also #-symmetric at each time t ≥ 0. Proposition 2.3 (or formula ( 14)) yields: q(t, x) = 2 -n x 1 ⋯x n p t, - q ∞ (x) φ(x) dx, as t → +∞, with x 1 1 (m) = ( √ -2m 1 , . . . , √ -2m n ) and:

q ∞ (x) = 2 -n x 1 ⋯x n p ∞ - x 2 1 2 , . . . , - x 2 n 2 .
The above formula corresponds to [START_REF] Martin | The fitness effect of mutations across environments: Fisher's geometrical model with multiple optima[END_REF]. Furthermore, since ∫ R n q ∞ (x) dx = 1 and since ∫ R n q(t, x) dx = 1, for every t ≥ 0, it then easily follows from the previous estimates that R n q(t, x) φ(x) dx → t→+∞ R n q ∞ (x) φ(x) dx for every φ ∈ C ∞ c (R n ). In other words, the distributions q(t, ⋅) weakly converge in R n to q ∞ as t → +∞. Lastly, the formula lim t→+∞ m(t) = m ∞ = ∫ R n m(x) q ∞ (x) dx is a consequence of the previous arguments and Theorem 2.10. The proof of Corollary 2.12 is thereby complete.

Plateaus

We show in this section that, given an initial phenotype x 0 = (x 0,1 , . . . , x 0,n ) ∈ R n , a value µ 1 > 0, and a duration T > 0, we can choose some positive real numbers µ 2 , . . . , µ n (here, n ≥ 2) such that the mean fitness:

t ↦ m(t) ∶= n i=1 x 2 0,i 2 tanh 2 (µ i t) -1 - µ i 2 tanh(µ i t) ,
is close to each of the plateaus:

m ∞,k = - n i=k+1 x 2 0,i 2 - k i=1 µ i 2
, for k = 1, . . . , n -1, at least during a time interval of duration T . We recall that the above formula for m(t) corresponds to the limit of formula [START_REF] Lenski | Dynamics of adaptation and diversification: a 10,000generation experiment with bacterial populations[END_REF], when the initial conditions q 0 approach the Dirac distribution δ x 0 . More precisely, we are given x 0 = (x 0,1 , . . . , x 0,n ) ∈ R n , T > 0, ε > 0 and µ 1 > 0, and we shall prove the existence of (µ 2 , . . . , µ n ) ∈ (R * + ) n-1 such that, for each k ∈ 1, n -1 , the set:

t ≥ 0, m(t) -m ∞,k ≤ ε , (62)
contains an interval of length at least equal to T . In that respect, we firstly define some functions:

s i ∶ (µ, t) ∈ R * + × R + ↦ x 2 0,i 2
S k (t) × (1 -tanh(µ k τ k )) ≤ ε n + k n , (63)
◇ a real number µ k+1 ∈ (0, µ k ) such that:

∀ i ∈ 1, k , ∀ t ∈ [τ i , τ i + T ], s k+1 (µ k+1 , t) tanh(µ k+1 (τ i + T )) ≤ ε 2n .
Note that the last property implies that:

∀ k ∈ 1, n -1 , ∀ i ∈ k + 1, n , ∀ t ∈ [τ k , τ k + T ], s i (µ i , t) tanh(µ i (τ k + T )) ≤ ε 2n . (64)
Fix now an index k ∈ 1, n -1 and a time t ∈ [τ k , τ k + T ] (⊂ R * + ). There holds: Thus, with this choice of (µ 1 , . . . , µ n ) ∈ (R * + ) n , each set defined in (62) contains an interval of length at least equal to T .

m(t) -m ∞,k = k i=1 x 2

A A formal derivation of the diffusive approximation of the mutation effects

The goal of this appendix is to give a formal justification of the diffusion term in (3). The case n = 1 is classical and can be found e.g. in [START_REF] Roques | Modèles de réaction-diffusion pour l'écologie spatiale[END_REF][START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF]. The anisotropic case n ≥ 2 is less standard, but it will easily follow from the same arguments.

Namely, we assume that the mutation effects on phenotypes follow a normal distribution N (0, Λ), with Λ = diag (λ 1 , . . . , λ n ) and λ i > 0 for each i ∈ 1, n , and that these mutations occur with a rate U > 0. In other words, the dynamics of the phenotype distribution under the mutation effects only (i.e., without selection) can be described by an integro-differential equation: ∂ t q(t, x) = U (J ⋆ q -q)(t, x), t ≥ 0, x ∈ R n , where ⋆ is the standard convolution product in R n defined by: (J ⋆ q)(t, x) = R n q(t, x 1 -y 1 , . . . , x n -y n ) J(y 1 , . . . , y n ) dy 1 . . . dy n , and J the (Gaussian) probability density function associated with the normal distribution N (0, Λ).

Formally, by writing a Taylor expansion of q(t, xy) at x:

q(t, x -y) = ∞ k 1 ,...,kn=0 (-1) k 1 +⋯+kn y k 1 1 ⋯y kn n k 1 ! ⋯ k n ! ∂ k 1 +⋯+kn q ∂x k 1 1 ⋯∂x kn n (t, x),
and by defining the central moments of the normal distribution:

ω k 1 ,...,kn = R n y k 1 1 ⋯y kn n J(y 1 , . . . , y n ) dy 1 . . . dy n , we then get that:

(J ⋆ q)(t, x) = ∞ k 1 ,...,kn=0
(-1) k 1 +⋯+kn ω k 1 ,...,kn k 1 ! ⋯ k n ! ∂ k 1 +⋯+kn q ∂x k 1 1 ⋯∂x kn n (t, x).

Since ω k 1 ,...,kn = 0 if at least one of the k i 's is odd, and since

ω k 1 ,...,kn = n i=1 k i ! 2 k i 2 (k i 2)! × n i=1 λ k i 2 i
, otherwise, one infers in particular that the second-order moments with even indexes are such that ω 0,...,0,k i =2,0,...,0 = λ i . Assuming that max 1≤i≤n λ i << 1, we may formally neglect the moments of order k 1 + ⋯ + k n ≥ 4, leading to:

(J ⋆ q)(t, x) ≈ q(t, x) + n i=1 λ i 2 ∂ ii (t, x).
Finally, setting µ 2 i = U λ i , we obtain:

U (J ⋆ q -q)(t, x) ≈ n i=1 µ 2 i 2
∂ ii q(t, x).

Figure 1 -

 1 Figure 1 -Trajectory of adaptation in the presence of anisotropy, with n = 3. The function m(t) is given by formula (32), and the approximations of the values of the intermediate plateaus are m∞,1 = -x 2 0,2 2 -x 2 0,3 2 -µ 1 2 and m ∞,2 = -x 2 0,3 2 -µ 1 2 -µ 2 2. The mutational parameters are µ 1 = √ 2, µ 2 = √ 2 ⋅ 10 -2 and µ 3 = √ 2 ⋅ 10 -4 .The other parameter values are x 0 = (3 2, 1, 1).

Figure 2 -

 2 Figure 2 -Trajectory of adaptation, anisotropic and isotropic model vs LTEE data. The function f (t) = m(t + 2) -m(2), with m(t) given by (32) was fitted to the LTEE data. The functions f (t) with the parameter values corresponding to the best fit are depicted in blue. The red crosses correspond to the LTEE data (population Ara-1, ln [w(6.64t)], with w the original data). The values leading to the best fit are (a, anisotropic model) x 0 = (0.73, 0.76), µ 1 = 1.3 ⋅ 10 -2 and µ 2 = 3.0 ⋅ 10 -4 ; (b, isotropic model) x 0 = (0.30, 0.09) and µ 1 = µ 2 = 5.3 ⋅ 10 -3 .

2 ,

 2 for all (t, x) ∈ R + × (R * ) n and, for each functionφ ∈ C ∞ c ((R * ) n ), R n q(t, x) φ(x) dx = 2 -n R n x 1 ⋯x n p t, -, m) φ # (x 1 1 (m)) dm → R n p ∞ (m) φ # (x 1 1 (m)) dm = R n

tanh(µ t) - µ 2 for i = 1 ,◇

 1 . . . , n -1.Secondly, by iteration for k = 1, . . . , n -1, we can then define:◇ a function S k ∶ t ∈ R + ↦ a time τ k > τ k-1 + T , (with τ 0 = -T ) such that: max t∈[τ k ,τ k +T ]

  and is independent of z ∈ R n + and t ∈ (0, T ]. One then infers from (47) and Lebesgue's dominated convergence theorem that the function M p is continuous in R + × R n + . As C p = log M p , the cumulant generating function C p is also continuous in R + × R n + . Let us then check that M p and C p are C 0,1 (R + × R n + ), meaning that the functions ∂ i M p = ∂Mp ∂z i and ∂ i C p = ∂Cp ∂z i exist and are continuous in R + × R n + , for every 1 ≤ i ≤ n. As a matter of fact, since p(t, ⋅) is a probability density function in R n -for any t ≥ 0 and since the integral ∫ R n -m p(t, m) dm converges by formula (13) in Proposition 2.3, it easily follows from Lebesgue's dominated convergence theorem that ∂ i M p (t, z) exists for all (t, z) ∈ R + × R n + and 1 ≤ i ≤ n, with:

  As t ∈ [τ k , τ k + T ], we have 0 < 1 -tanh(µ i t) ≤ 1 -tanh(µ k τ k ) for every i ∈ 1, k (remember that 0 < µ k < µ k-1 < ⋯ < µ 1 ), whereas tanh(µ i t) ≤ tanh(µ i (τ k + T )) for every i ∈ k + 1, n . It then follows from (63)-(65) that, for every k ∈ 1, n -1 and t ∈ [τ k , τ K + T ], m(t)-m ∞,k ≤ S k (t) 1-tanh(µ k τ k ) +

				n i=k+1 s i (µ i , t) tanh(µ i (τ k +T )) ≤ ε	n + k 2n	+ i=k+1 n	ε 2n	= ε.
			0,i 2	tanh 2 (µ i t) -1 -	µ i 2	tanh(µ i t) -1
	≤	k i=1	x 2 0,i 2	+ + tanh(µ i t) + 1 -n i=k+1 x 2 0,i 2 n i=k+1 0,i 2 x 2	tanh 2 (µ i t) -tanh(µ i t) -µ i µ i 2 2 2 1 -tanh(µ i t) tanh(µ i t) tanh(µ i t). µ i	(65)

Acknowledgments

This work was supported by the French Agence Nationale de la Recherche (ANR-13-ADAP-0016 "Silentadapt" to G.M., ANR-13-ADAP-0006 "MeCC" to L.R., ANR-14-CE25-0013 "NONLOCAL" to F.H. and L.R. and ANR-18-CE45-0019 "RESISTE" to F.H., G.M. and L.R.). This work was fostered by stimulating discussions with Thomas Lenormand.

and (57) holds in R + × R n + . Owing to the definition C p = log M p , one concludes that ∂ t C p is continuous in R + × R n + (finally, C p is of class C 1,1 (R + × R n + )) and:

for all (t, z) ∈ R + × R n + , where A and b are as in [START_REF] Johnson | The approach to mutation-selection balance in an infinite asexual population, and the evolution of mutation rates[END_REF]. Therefore, [START_REF] Hamel | On the nonlocal Fisher-KPP equation: steady states, spreading speed and global bounds[END_REF] holds in R + × R n + and the proof of Theorem 2.6 is thereby complete.

Before going into the proof of the remaining results, let us first observe that, in ( 21)- [START_REF] Johnson | The approach to mutation-selection balance in an infinite asexual population, and the evolution of mutation rates[END_REF], m(t) = 1 ⋅ ∇C p (t, O), with:

It turns out that, if A(z) = 1, then equation ( 21) can be solved explicitely by the method of characteristics, as the following lemma shows (this lemma is used later in the proof of Proposition 2.8 in the general case A(z) given in [START_REF] Johnson | The approach to mutation-selection balance in an infinite asexual population, and the evolution of mutation rates[END_REF]).

Lemma 4.3. The Cauchy problem:

) solution, which is given by the expression:

Proof. First of all, it is immediate to check that the function Q given by (59

and obeys:

. It follows from the equation satisfied by Q that, for all s ∈ [0, t], there holds:

Proof of Proposition 2.8. In order to derive a general formula for the C 1,1 (R + × R n + ) solution C p of ( 21), we make a substitution of the spatial variable and use the previous special case described in Lemma 4.3. To do so, we set, for t ≥ 0 and z ∈ R n + ,

where y(z) = (y 1 (z), . . . , y n (z)) and

where A is given in [START_REF] Johnson | The approach to mutation-selection balance in an infinite asexual population, and the evolution of mutation rates[END_REF]. As 13) and [START_REF] Good | The impact of macroscopic epistasis on long-term evolutionary dynamics[END_REF], it follows from (21) that:

Therefore, Lemma 4.3 implies that:

for all (t, z) ∈ R + × R n + . Consider now any t ∈ R + and z = (z 1 , . . . , z n ) ∈ [0, 1 µ 1 ) × ⋯ × [0, 1 µ n ). Set:

and observe that y(z ′ ) = z. Hence, we have:

which leads straightforwardly to the formulae ( 24)- [START_REF] Lenski | Experimental evolution in Escherichia coli. I. Adaptation and divergence during 2,000 generations[END_REF]. Furthermore, for every t ∈ R + , the formula m(t) = 1 ⋅ ∇Q(t, O) together with (60) easily yields [START_REF] Lenski | Dynamics of adaptation and diversification: a 10,000generation experiment with bacterial populations[END_REF]. The proof of Proposition 2.8 is thereby complete.

Proof of Corollary 2.7. Let p be the fitness distribution, that is, p(t, m) dm is the pushforward measure of q(t, x) dx by the map x ↦ -x 2 2. Let:

e zm p(t, m) dm, be the moment generating function of p. As the fitness m ∈ R -is the sum of the fitness components (m 1 , . . . , m n ) ∈ R n -, we have:

Thanks to the equations ( 21)-( 22) satisfied by C p , it follows that

for all (t, z) ∈ R + × R + . This is the desired result and the proof of Corollary 2.7 is thereby complete.

Proof of Corollary 2.9. We have seen in the proof of Corollary 2.7 that C p (t, z) = C p (t, z1) for all (t, z) ∈ R + × R + . Thus, formulae ( 24)-( 26) straightforwardly yield ( 27)-( 28) for t ≥ 0 and z ∈ [0, 1 µ).

Stationary states

Proof of Theorem 2.10. Formulae ( 24)-( 25) in Proposition 2.8 imply that, for every z ∈ [0, 1 µ 1 ) × ⋯ × [0, 1 µ n ), one has:

as t → +∞. It then follows from the generalization of the Curtiss theorem [START_REF] Yakymiv | A generalization of the Curtiss theorem for moment generating functions[END_REF] that, if the limit as t → +∞ of the cumulant generating functions C p (t, ⋅) is the cumulant generating function of p ∞ given by ( 29) in some subset of R n -with non-empty interior, then the distributions p(t, ⋅) weakly converge to p ∞ as t → +∞. So let us compute the CGF of p ∞ . For all z ∈ R n + , Fubini's theorem yields that:

Hence, the CGF of p ∞ is equal to the function z = (z 1 , . . . , z n ) ↦ -(1 2) ∑ n i=1 log (1 + µ i z i ), that is, the limit in (61). As a consequence, the distributions p(t, ⋅) weakly converge to p ∞ in R n -as t → +∞. On the other hand, thanks to Proposition 2.8, we also know that:


In this work, linear stability analysis and Direct numerical simulations (DNS) are used to investigate the coupling between hydrodynamic instabilities, membrane transfer and osmotic pressure in Taylor-Couette configuration. The emphasis is on characterizing the effect of the osmotic pressure related to the concentration boundary layer forming near the membrane on the structure and dynamics of Taylor vortices. We consider a Taylor-Couette cell with two semi-permeable membranes totally rejecting the solute transported by a newtonian fluid filling the gap. For fixed operating conditions, linear stability analysis shows that the osmotic pressure tends to alter centrifugal instabilities as a result of an original self-sustained mechanism coupling the advection of the concentration boundary layer by the vortices, molecular diffusion and osmotic pressure driving a transmembrane flow fostering the vortices. This mechanism can induce a substantial reduction of the critical rotation rate above which vortices are observed. Furthermore, stability analysis shows that critical conditions are also impacted by the radius ratio. These analytical results are compared to recent DNS. *Financial support by Thomas Jefferson Fund.