INFLUENCE OF THE MODE OF REPRODUCTION ON DISPERSAL EVOLUTION DURING SPECIES INVASION

We consider a reaction-diffusion-reproduction equation, modeling a population which is spatially heterogeneous. The dispersion of each individuals is influenced by its phenotype. In the literature, the asymptotic propagation speed of an asexual population has already been rigorously determined. In this paper we focus on the difference between the asexual reproduction case, and the sexual reproduction case, involving a non-local term modeling the reproduction. This comparison leads to a different invasion speed according to the reproduction. After a formal analysis of both cases, leading to a heuristic of the asymptotic behaviour of the invasion fronts, we give some numerical evidence that the acceleration rate of the spatial spreading of a sexual population is slower than the acceleration rate of an asexual one. The main difficulty to get sharper results on a transient comes from the non-local sexual reproduction term. Résumé. Nous nous intéressons à un modèle de réaction-diffusion-reproduction d’une population, répartie de manière hétérogène. Le phénotype de chaque individu caractérise sa dispersion. Dnas la littérature, la vitesse asymptotique de propagation d’une population asexuée a déjà été rigoureusement déterminée. Dans cet article, nous comparerons la dispersion d’une population asexuée, avec celle d’une population sexuée, ce qui introduit un terme non-local modélisant la reproduction. Après une étude formelle des deux modes de reproduction, qui nous amènera vers une heuristique du comportement asymptotique du front d’invasion, nous effectuons des simulations numériques, montrant que la vitesse d’accélération d’une population sexuée est plus faible que celle d’une population asexuée. La principale difficulté pour obtenir des résultats, sur le régime transitoire, vient du terme non local modélisant la reproduction sexuée. 1 Unité de Mathḿatiques Pures et Applications (UMPA), Ecole Normale Supérieure de Lyon, France; e-mail: vincent.calvez@ens-lyon.fr 2 Institut de Mathématiques de Toulouse ; UMR5219, Université de Toulouse ; UPS IMT, F-31062 Toulouse Cedex 9 France; e-mail: joachim.crevat@math.univ-toulouse.fr 3 École Normale Supérieure de Lyon; e-mail: leonard.dekens@ens-lyon.fr 4 Institut Camille Jordan, Université Claude Bernard Lyon 1 in Lyon e-mail: fabreges@math.univ-lyon1.fr 5 Laboratoire de Mathématiques et de Physique Théorique de Tours; e-mail: frederic.kuczma@lmpt.univ-tours.fr 6 BioSP, INRA, 84914, Avignon, France; e-mail: florian.lavigne@inra.fr 7 Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France 8 Centre de Mathématiques Appliquées, École Polytechnique, Palaiseaue-mail: gael.raoul@cmap.polytechnique.fr © EDP Sciences, SMAI 2020 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Article published online by EDP Sciences and available at https://www.esaim-proc.org or https://doi.org/10.1051/proc/202067008 ESAIM: PROCEEDINGS AND SURVEYS 121


Introduction
Biological populations are often evolving in a heterogeneous environment, to which individuals can be more or less fitted, according to their phenotypes. Phenotypic distributions are shaped through processes like survival, reproduction and migration. Biological invasions are archetypal of a process where spatial heterogeneities are constitutive. Successful species manage to thrive due to some selective advantage in a new environment, and to expand in a new area. Different causes can promote the success of an invasion: wind or birds can disperse seeds (e.g. widespreading of lodgepole pine in western North America [12]) or climate change modifies the local environment and may influence species range (e.g. invasion of butterfly and bush cricket species in Britain [21]). Also, many species invasions result from an introduction, either voluntary or involuntary, by human beings (e.g. invasion of cane toads in Australia [1,19,20]), etc.
Recently, much attention was paid to evolutionary dynamics at the vanguard of species range expansion, see for instance [4,20,22,23]. Among other, dispersal evolution, that is the evolution of the dispersal abilities of individuals, was the subject of biological and mathematical studies in the past decade. In particular, the phenomenon of spatial sorting, i.e. the assembly of more dispersive individuals ahead of the range expansion by pure spatial effect (all other phenotypes being identical within the population), was the subject of several modelling and analytical studies [2,6,7,10,24]. In particular, it has been established that front acceleration results from sustained spatial sorting in the case of possibly unbounded dispersal rates (or equivalently in the transient regime before the physiological limits are reached) [3,5,9]. So far, analytical results were obtained based on the modeling assumption that reproduction is clonal (as for asexual reproduction) with possible mutations affecting the dispersal ability of offsprings. In particular, recent techniques borrowed from the approximation of geometric optics for front propagation were used successfully to compute quantitative features such as the asymptotic speed of propagation, or the rate of acceleration in the transient regime [9,24].
Besides, similar techniques were developed recently in a different context: the asymptotic description of equilibria in quantitative genetics models involving a sexual mode of reproduction in the regime of small variance. This methodology was recently completed for Fisher's infinitesimal model [14]. It is a simple model for the inheritance of quantitative traits in a population with a sexual mode of reproduction. It assigns to an offspring the average parental phenotypic trait up to a random normal deviation with constant variance [5].
Preliminary heuristics seem to show that this mode of reproduction significantly slows down the wave expansion as compared to an asexual reproduction. It is expected that, in the transient accelerating regime, the population spreads as t 5 4 under the infinitesimal model of inheritance, whereas it was previously shown that the spreading occurs at a rate of order t 3 2 with the asexual mode of reproduction [3,8].
In this paper, we present formal calculations in order to compare the rates of expansion for the two modes of reproduction under study. To validate these heuristics, we present also numerical results, which help us to catch a transient regime and its main features.

Our model
In this paper, we study the time evolution and the spatial and phenotypical propagation of population encowed with either an asexual or a sexual reproduction mode. More precisely, the population is distributed according to its location x ∈ R and its phenotypic trait θ ∈ (θ min , +∞), with θ min > 0, and is modeled with the following reaction-diffusion equation for all t > 0: where r 1 > 0 and K > 0 are fixed constants, f (t, x, θ) is the density of individuals presenting the trait θ at location x ∈ R at time t ≥ 0, and ρ(t, x) ∶= ∫ ∞ θ min f (t, x, θ) dθ is the population size at x ∈ R and time t > 0. The reaction term B[f ] will be detailed later. We also assume that initially, the density support is compact. Let us discuss the modeling interest of each term. Firstly, the term r 1 f 1 − K −1 stands for an adaptation force called selection. At point x ∈ R and at time t ≥ 0, any individual is in competition with the others for resources: when the density (t, x) at x is less than a threshold K, named the carrying capacity, there is enough resources for everyone so the population can grow; whereas, if (t, x) > K, then competition is involved, and the population is decreasing. The constants r 1 > 0 is therefore called competition rate.
Then, the term θ ∆ x f models the migration phenomenon. Individuals are assumed to diffuse through space at each time t, at a rate given by the phenotypic trait θ. For example, those with long legs or bigger wings can be faster to go into a new environment, and so help the invasion. In fact, the equation (1) without the last reaction term B[f ] is the generalized non-local Fisher-KPP equation, where the diffusion depends on the phenotypic traits, which can be seen as a given constant. In the case that is a convolution term with f , the previous equation has been studied in some different ways [13,[15][16][17].
Finally, the last reaction term B[f ] models the reproduction event. Let us turn to explain this new term. Asexual case. For asexual populations, reproduction is clonal: the offspring receives at birth the same trait as its sole parent, but mutations change the value of the trait. Assuming the variance of the mutation effects is very small, we can model the phenomenon by the following evolutionary equation, as in [2]: where µ > 0 is a constant depending on the mutation rate and on the variance of the mutation effects. Replacing f by θ min K −1 f (θ 2 min t µ, θ 3 min µ x, θ min θ) yields us to simplify the initial problem and to study: with θ ∈ (1, +∞), r > 0 and (t, In this case, the reproduction term is just Bf (t, x, θ) = ∆ θ f (t, x, θ).
Remark. The variation deduced by the mutation event is an approximation of an infinitesimal model [14,18], taking the limit of the equation when the variance of the mutation effect λ 2 > 0 tends to 0. In this infinitesimal model, this variation follows a centred normally distributed random variable. Without selection and migration, the mutations [18] can be modeled by the evolutionary equation: where r 2 > 0 is the mutation rate, J is the gaussian distribution of N (0, λ 2 ) and the "convolution term" x, θ ′ )dθ ′ . Linearising the previous equation, when λ is very small, we retrieve the same term as in (3), with µ = r 2 λ. Boundary condition for the asexual case. By definition of phenotypic trait, no individual can have a phenotype in (−∞, 1) or generate a descendent with such phenotype: there is no phenotypical exchange at θ = 1, which means that: Sexual case. During sexual reproduction, at the opposite of the asexual one, recombinations of the DNA happen, which changes the mathematical model and so the reproduction term in (2). We assume that individuals can breed only with those which are at the same location x ∈ R.
Let us turn to the sexual reproduction model, introduced by Bülmer in 1980's. At time t, an individual of trait θ 1 can find and copulate with a partner of trait θ 2 with probability density f (t, x, θ 2 ) (t, x). The offspring has originally the trait (θ 1 + θ 2 ) 2, which suffers directly from recombinations of variance λ 2 > 0: with r 3 > 0 the rate of reproduction/recombination. The term G θ θ − θ1+θ2 2 , symbolizing the stochasticity of the recombinations, is a normalized truncated centred gaussian, with variance λ 2 > 0, such that θ − θ1+θ2 2 > θ min . As a prospect, a further article could study the influence of the variance λ 2 of the recombination event; for the sake of simplicity, we take here λ 2 = 1 2. Replacing here f by θ min K −1 f (t r 3 , θ min r 3 x, θ min θ), and noting r ∶= r 1 K θ min , we can simplify the previous PDE into: with G θ is a normalized truncated centred gaussian, such that θ − θ1+θ2 2 > 1. By this simplification, the reproduction is:

Formal analysis
If we consider the asexual reproduction case without the competition, then the PDE that we obtained becomes linear and therefore calls for exponential solutions. Aligning thereby with others studies about diffusion of individuals and speed of the invasion (see for instance [8,9]), where the common method is to define the function u such that: where s is a time parametrization, and α and β are introduced to re-scale the spatial and phenotypic variables.
Finding (α, β) crystallizes the main issue tackled by the formal analysis, as it is the crucial key to get asymptotically a non trivial PDE satisfied by u (that is taking into account all the biological forces, namely selection, migration and reproduction), as it is shown in figure 1. Asexual case. Thanks to (3), we can check that: , y, η))dη. We choose the parametrization s = log(t) so that ts ′ (t) = 1. In this case, we obtain an viscous Hamilton-Jacobi EDP : Remark. In large time t, which is equivalent to taking the limit as s tends to +∞, the factor exp(− Cste s) in front of the laplacians terms will make the second order terms vanish. Nevertheless, the transformation (5) enables to capture the long time asymptotics. In fact, it transforms the nature of the problem, from a second order parabolic equation to a nonlinear viscous Hamilton-Jacobi, for which characteristic lines can be analytically computed. However, reviewing the rigorous derivation of the viscous Hamilton-Jacobi equation relies on the theory of viscosity solutions, which is beyond our scope, and we refer to [9] for details.
Formally, our purpose is to determine the asymptotic behaviour of (6) in the limit s → +∞, that is to derive a PDE in which all biological phenomena intervene, especially: ◇ the reproduction term: as the factor e (2−2β)s (∂ η u) 2 tends to +∞ or 0 depending on the sign of (2 − 2β), the reproduction term will either annihilate the contribution of the others or vanish. Thus the only relevant choice derives from 2 − 2β = 0. ◇ the diffusion term: the same considerations lead us to 2 − 2α + β = 0. We thereby get the following system for (α, β): which leads to: Sexual case. Although the PDE obtained without taking the competition in account is in this case non linear, it keeps the same homogeneity as the asexual one, for only the reproduction term differs between the two models. We thus proceed by analogy with the asexual case by looking for exponential solutions using the transformation (5). Taking again s(t) = log(t), equation (4) yields: with = e βs ∫ ∞ e −βs exp(−e s u(s, y, η))dη and: (8) As the sole difference between the asexual and sexual cases is the reproduction term, only the condition on the reproduction term will differ (and not the condition on the diffusion term that is coupling α and β). Since it dictated entirely the choice β in the previous case, we have to look precisely at the term B[f ] f in order to deduce formally a reasonable β. By the same heuristics as in the previous case, as we want reproduction term neither to vanish nor to undermine the contributions of the other terms as s goes to +∞, we need all the terms in the double integral to be treated equally with respect to s, which yields 1 − 2β = 0. Thereby, we get the following system: which leads to:

Schemes and numerical results
Here, we present some numerical schemes to approximate the solutions of the equations (3) and (4). During this CEMRACS project, we have simulated the solutions of the different PDEs with Fortran and Python. We have used the packages numpy, matplotlib and pythran.

Asexual case
In this section, we present two different numerical schemes to approximate the diffusion with asexual reproduction.
We consider x max ≥ 0 and θ max ≥ 1 so that we work with couples (x, θ) in the bounded domain [0, x max ] × [1, θ max ], discretized with the meshes (x i ) 1≤i≤Nx and (θ j ) 1≤j≤N θ , respectively of step ∆x > 0 and ∆θ > 0. As for the time discretization, let ∆t > 0 be a time step, and let us define for all n ∈ N, t n ∶= n ∆t.
In the following, we denote by A D x and A N x the matrix of the discrete Laplace operator in x of size N x respectively with Dirichlet and Neumann boundary conditions: We consider A D θ and A N θ the matrix of the discrete Laplace operator in θ of size N θ respectively with Dirichlet and Neumann boundary conditions with similar definitions. In the following, we choose to work with the fixed parameter r = 1.

Euler explicit scheme
Then, for all n ∈ N, we approach ( (t n , x i )) 1≤i≤Nx by a vector (̃ n i ) 1≤i≤Nx computed with an approximation of the integral in θ using the trapezoidal rule. We have chosen an initial truncated Gaussian distribution: First of all, we approximate the reaction-diffusion equation (3). We start by discretizing the diffusion terms in x and θ with the matrix of discrete Laplace operator with Dirichlet boundary conditions A D x and A D θ . Then, we approximate in time using an explicit Euler scheme that is for all n ∈ N, we compute the new value of F n+1 by: where D n ∶= diag ((̃ n i ) 1≤i≤Nx ) , and D θ ∶= diag ((θ j ) 1≤j≤N θ ) . In the following, the steps ∆t, ∆x and ∆θ have been chosen to satisfied the stability condition: In Figure 1 (a), we plot the density of population (t, ⋅) at different fixed times at regular intervals from the instant t = 10 to t = 150. As expected, the function seems to quickly converge towards an invasion front propagating towards right connecting 0 to 1. Moreover, we observe that the front of propagation in space seems to accelerate. Theoretically, according to [3,8], we should observe a propagation of the front of order t 3 2 in space and of order t in phenotype.
On Figure 1 (b), we verify the order of the spatial acceleration plotting the same curves with respect to the auto-similar variable y = xt −3 2 . We observe that the spatial densities converge towards an Heaviside distribution in y, centered at a constant y c which is approximately 1.25. Theoretically, according to [8], this constant is approximately y c ≈ 1.315.
Then, we focus on the shape of the transition front. We plot in Figure 1 (c) the spatial distribution with respect to a re-centered scale in X 1 2 (t), where (t, X 1 2 (t)) = 1 2 for all t ≥ 0. We can observe that the shape of the front flattens as time goes to infinity. In order to study this evolution, we display in Figure 1 (d) the same curves with respect to the re-scaled variable x − X 1 2 (t) t −1 2 . We can see that all curves are superposed. Therefore, the shape of the front seems to flatten at order t 1 2 .
Finally, in order to numerically study the accelerations in phenotype of the solution of the equation (3), we define for all t ≥ 0 the position of the front at time t with: and the main phenotype at the head of the frontθ(t) with: .
In Figure 1 (e), we show in log-log scale the evolution of the numerical and formal results ofθ with respect to time, and the identity function t ↦ t, in order to compare their slopes. We can observe that for high times, these two curves seem to be parallel. Using a linear regression without taking into account the values where t is too small in order not to consider an eventual transitory state, we find that our numerical approximation of the main phenotypeθ(t) is of order t 1.01 . Therefore, this explicit Euler scheme seems to be efficient to reproduce the expected acceleration phenomenon in the asexual case.

Auto-similar variables
In this section, we focus on the approximation of the variable u(s, y, η) defined in (5), with reparametrized time s ∶= log(t), and auto-similar variables y ∶= x t −α and η ∶= θ t −β with α and β > 0. Thus, the equation (6) satisfied by u can be written as follows: where for all s > 0, and p ∈ R, the Hamiltonians in y and η are defined with: To be consistent with the last paragraph, we also choose a polynomial initial condition: We consider y max ≥ 0 and η max ≥ 1 so that we work with couples (y, η) in the bounded domain [0, y max ] × [1, η max ], discretized with the meshes (y i ) 1≤i≤Ny and (η j ) 1≤j≤Nη , respectively of step ∆y > 0 and ∆η > 0. In the following, we define A N y and A N η the matrix of the discrete Laplace operator with Neumann boundary conditions respectively in y of size N y and in η of size N η . We also define a time step ∆s > 0 and for all n ∈ N, s n = n ∆s. We want to approximate for all 1 ≤ i ≤ N y , 1 ≤ j ≤ N η , and n ∈ N the function u(s n , y i , η j ) with a matrix U n ij ∈ R Ny×Nη , defined for all n ∈ N with: where D η ∶= diag (η j ) 1≤j≤Nη , and for all 1 ≤ i ≤ N y , 1 ≤ j ≤ N η , wherẽ n i is the approximation of exp(s n ) , y i × s 5 4 n computed from U n . Furthermore, the matrix H 1 (s n , U n ) and H 2 (s n , U n ) ∈ R Ny×Nη are the respective approximations of the Hamiltonians in y and η. To choose good approximations H 1 and H 2 , we used the idea of the scheme of [11]. This method consists in rewriting the Hamiltonian H 1 (s, ⋅) (resp. H 2 (s, ⋅)) into the sum of an increasing function H + 1 (s, ⋅) and a decreasing function H − 1 (s, ⋅), which gives for all s > 0, y ∈ R, η > 0 and p ∈ R: Then, we approximate H 1 (s, ∂ y u) (resp. H 2 ) with: where ∂ − y u and ∂ + y u (resp. ∂ − η u and ∂ + η u) are the derivatives of u with respect to y (resp. η) respectively on the left and on the right. Thus, if we define the matrices D l y and D r y of size N y × N y and D l η and D r η of size N η × N η with: and: we approximate H 1 (s n , ∂ y u)(y i , η j ) and H 2 (s n , ∂ η u)(y i , η j ) with: At the end of an iteration, we compute F by the relation: In Figure 2 (a), we display the function −t log(f ) at fixed time t = 150 with respect to the re-scaled variable y = x t −3 2 and η = θ t −1 , where f is the solution of the reaction-diffusion (3) in the asexual case. We can compare this plot with the Figure 2 (b), where we display the function u at fixed re-parametrized time s = 5, with respect to y and η. First of all, the shape of these two functions seems to be roughly similar. Nevertheless, since f (t, x, θ) is close to 0 for large enough x and small enough θ, the function −t log(f ) drawn in Figure 2 (a) reaches high values for y = 1.6 and η = 0. On the contrary, the function u plotted in Figure 2 (b) does not explode for large y and small η. Therefore, the numerical scheme with auto-similar variables provides more information than the explicit Euler scheme explained in Section 3.1.1 on the behaviour of the solution at the head of the front.
Furthermore, in Figure 2 (c), we show the density population (s, y) computed from the solution u of the equation (6), with respect to the re-scaled variable y, at different fixed times from s = 0 to s = 7 at regular interval. We observe that the function seems to converge as expected towards an Heaviside function H(y c −y). This is coherent with the theoretical results explained previously in Section 2. Moreover, the numerical estimate for y c ≈ 1.32 is close to the theoretical value from [9] (approximately 1.315), and is consistent with theorem 1.2 from [8] which establishes y c ≤ 4 3 and an acceleration of order t 3 2 . Consequently, this numerical scheme enables us to approximate more accurately the order of acceleration of the front than the previous one explained in Section 3.1.1.
Besides, the position of the front y c depends on the value of the birth rate r. For small values of this parameter, we expect y c to be increasing with respect to r. Indeed, allowing individuals to mutate according to a stronger rate would lead to a larger range in the dispersal rates and build a faster front wave for the invasion. This phenomenon can be illustrated through a comparison with [3], in which the authors have analytically studied the same model with a different set of parameters. Up to a time scale, their model is equivalent to ours with r = 2 and we see that: y c (r = 2) = 2 3 2

Sexual case
The most challenging problem for the simulation of the sexual case is the approximation of the reproduction term B[f ]. In this article, a straightforward approach is sufficient to display what we seek. Thus, for all n ∈ N, Similarly as in the asexual case, we discretize the time using an explicit Euler scheme, that is for all n ∈ N, we compute the matrix F n+1 using: with D n ∶= diag ((̃ n i ) 1≤i≤Nx ) , and D θ ∶= diag ((θ j ) 1≤j≤N θ ) . In Figure 3, we show some simulations of the solution of the equation (4). First of all, in Figure 3 (a) and (b), we display the contour lines of the function f respectively at fixed time t = 800 and t = 1, 600. We can observe the propagation both in space and phenotype of the function f as time passes. Moreover, the shape of the function f before the head of the front seems to be invariant with respect to time. Therefore, once the population is settled at a position x, it does not evolve any more.
Then, let us more precisely focus on the propagation in space. In Figure 3 (c), we plot the density function for different fixed times between t = 40 and t = 1600. As expected, if we do not take into account the values of x too close to the boundary, seems to quickly converge towards an invasion front. Moreover, this invasion front seems to accelerate. One can notice that for small x, the density function (t, x) seems to remain below 1 for all time. Indeed, unlike in the asexual case, this is possible since in the sexual case (4), the maximum principle is false.
To validate our formal analysis, we display in Figure 3 (d) the same functions with respect to the auto-similar variable y = xt −5 4 . As expected, as t increases, the function tends towards a stiff and stationary state. This shows that the acceleration of the invasion front showed in Figure 3 (c) seems to accelerate at the order y c t 5 4 . Nevertheless, similarly as in the asexual case, the approximation of the constant y c is not sharp.
Finally, as done in Section 3.1.1, let us study the deformation of the invasion front. In Figure 3 (e), we display the same functions as previously with respect to the re-scale and re-centered variable (x − X 1 2 (t)) t −1 4 . Since the graphs in Figure 3 seem to be merged, we can conjecture that the the typical width of the front in the sexual reproduction case is of order t 1 4 .
Finally, in order to study more precisely the acceleration of the invasion front in phenotype, we proceed as previously for the asexual case, computing for all time t ≥ 0 the front position X(t) and the mean phenotypē θ(t) at position X(t), respectively defined with (11) and (12). In Figure 3 (f), we display in log-log scale the numerical approximation of the mean phenotypeθ with respect to time, and the function t ↦ √ t, in order to compare the slope of the two curves for large times. A formal derivation and approximation of the multiplicative constant in the expression of the mean phenotype will be the objects of an upcoming work. We can observe that for t large enough, the two graphs seem to remain parallel. Proceeding with a linear regression, we find thatθ is approximately of order t 0.52 . In the same way that the explicit Euler scheme provides a good approximation of the asexual case (3), here the numerical results seem to fit with our formal analysis in Section 2.

Conclusion
Dispersal evolution at the vanguard of species range expansion can bring intricated dynamics such as spatial sorting [20] and front wave acceleration [3,8]. As mentioned in the Introduction, some heuristics predict that the sexual reproduction slows down the wave expansion in comparison to the asexual reproduction. In this paper, we have developed numerical schemes to bring elements of confirmation to this phenomenon.
As seen in section 4, the numerical simulations supports the formal analysis in confirming the predicted asymptotic order of the propagation speed in space and in the phenotypic range, for both the asexual and sexual type of reproduction (respectively (t 3 2 , t) and (t 5 4 , t 1 2 )). Therefore, it does appear that the mode of reproduction does bear a strong influence on the speed of the invasion, as well as on the genetic mixing: a species using a sexual reproduction will be slower to expand on both levels than a species using an asexual reproduction.
However satisfactory bringing a numerical confirmation to this question is, the formal analysis performed for the sexual case in Section 3 could not be considered as complete, as it does arise some new remarks. The knowledge of the values of the rescaling constants (α, β) = (5 4, 1 2) allows us to take a new look at the equation involving the re-scaled variables, obtained by the transformation (7): with = e s 2 ∫ ∞ e −s 2 exp(−e s u(s, y, η))dη and: The presence of factors like e s , only in terms where η is variable at a given point of space y, implies that there is a hierarchy between the different equilibria in the variables η and y on how fast they are reached. More precisely, it suggests that we can consider that the equilibrium towards the variable η is always reached at each spatial point y. Moreover, its mathematical description should clearly state its dependence towards the space variable y, which surely calls for a more thorough analysis.
Numerically, we have observed in the asexual case that a stable code computing u in terms of the re-scaled variables (s, y, η) was far more efficient to capture the complexity of the asymptotic state. Thus, the isolation of such subtle phenomena of hierarchy between phenotypic and space equilibria requires a similar code in the sexual case. Both of these are prospects of an upcoming work. to the re-centered variable x − X 1 2 (t), and (d) to the re-scaled variable (x − X 1 2 (t))t −1 2 , where for all t ≥ 0, (t, X 1 2 (t)) = 1 2. (e) Plot of the main phenotype θ(t) (see (12)) at the front position with respect to time (blue curve) and of the function t → t (red curve), in log − log scale.  , and (e) to the re-scaled and re-centered variable (x − X 1 2 (t))t −1 4 , where for all t ≥ 0, (t, X 1 2 (t)) = 1 2. (f) Plot of the main phenotype θ(t) (see (12)) at the front position with respect to time in log − log scale.