Matthieu Alfaro

Quentin Griette

Denis Roze

Benoît Sarels

Part I: The case of stacked underdominant clines

Keywords: genetic incompatibilities, underdominance, quasi-linkage equilibrium, standing wave, traveling wave, perturbation analysis AMS Subject Classifications: 92D10 (Genetics and epigenetics), 35C07 (Traveling wave solutions), 35B20 (Perturbations in context of PDEs)

published or not. The documents may come

The spatio-temporal dynamics of interacting genetic incompatibilities.

Introduction

Genetic incompatibilities correspond to deleterious interactions between alleles (at the same locus or at different loci) within the same genome, and are the cause of the reduced fitness of hybrids between species [START_REF] Coyne | Speciation[END_REF], [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF]. Such incompatibilities may be revealed by crosses between divergent populations or species, which may be performed experimentally [START_REF] Fraïsse | The genetics of speciation: are complex incompatibilities easier to evolve?[END_REF], but may also occur naturally in hybrid zones resulting from secondary contacts between genetically divergent populations [START_REF] Barton | Analysis of hybrid zones[END_REF][START_REF] Barton | Adaptation, speciation and hybrid zones[END_REF]. Indeed, the offspring of such crosses carry a mix of alleles from the two parental populations, which may not function well together. Genetic incompatibilities may also be widespread within the same species, as suggested by recent data from Drosophila [START_REF] Corbett-Detig | Genetic incompatibilities are widespread within species[END_REF].

How several incompatibilities segregating within the same population interfere with each other has important consequences for the evolution of reproductive isolation, and the maintenance of distinct genetic entities after a secondary contact. Barton & de Cara [START_REF] Barton | The evolution of strong reproductive isolation[END_REF] showed that, in the case of a single population (sympatry), incompatibilities are expected to couple through the buildup of linkage disequilibrium among them, which may eventually lead to strong reproductive isolation between two distinct genetic backgrounds. In spatially structured populations with limited dispersal (parapatry), clines in allele frequencies may form between regions containing different sets of incompatible alleles [START_REF] Barton | The dynamics of hybrid zones[END_REF]. In this case, linkage disequilibria generated by dispersal generally tend to pull clines together [START_REF] Slatkin | Gene flow and selection in a two-locus system[END_REF], again leading to the coupling of genetic incompatibilities which then tend to reinforce each other. This process was explored by Barton [3] in the case of a continuous, linear habitat, and when genetic incompatibilities are generated by an arbitrary number of underdominant loci: at each locus, heterozygotes (say Aa) have a lower fitness than homozygotes, while the two homozygotes (AA and aa) have equal fitness. This form of symmetric selection (against heterozygotes) can maintain stable clines in allele frequencies [START_REF] Bazykin | A hypothetical mechanism of speciation[END_REF], [START_REF] Barton | The dynamics of hybrid zones[END_REF], separating regions where AA and aa individuals are abundant (Aa hybrids being generated between these regions). When two such clines overlap in space (one separating regions where AA and aa are abundant, and the other separating regions where BB and bb are abundant), they tend to attract each other until they coincide, as illustrated by the numerical simulations in Figure 1, and then become motionless.

In the asymmetric situation where one allele has a selective advantage over the other (i.e., one homozygous genotype, say AA, has a higher fitness than the other, aa), the cline will move in the direction of the less fit genotype [START_REF] Barton | The dynamics of hybrid zones[END_REF]. The interaction between several such asymmetric incompatibilities raises several questions that remain little explored to date, such as: when clines moving at different speeds come into contact, will they remain stacked (increasing the degree of reproductive isolation between the two sides of the clines) or not? If they do remain stacked, what will be the speed of the resulting front? Under which conditions may an asymmetric incompatibility escape from a hybrid zone in which several incompatibilities are segregating?

This article constitutes a first step in the exploration of the spatio-temporal dynamics of interacting asymmetric incompatibilities, focusing on a simple situation involving two coupled underdominant loci (with alleles A, a at the first locus, B, b at the second) with identical fitness effects. Notice that Barton [START_REF] Barton | Multilocus clines[END_REF] has considered the situation where heterozygotes present a cost in fitness and where the fitness of the homozygotes AB|AB and the one of ab|ab are the same, that is the symmetric case. In Section 3, we give a mathematical proof of the existence of such a cline in this symmetric situation, see Proposition 3.2. We also prove in Proposition 3.3 that such a stationary cline is stable, i.e. that small perturbations of the profile may shift its spatial position but essentially do not alter its shape. In other words, a perturbed stationary cline comes back to a possibly shifted stationary cline. We refer to [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF] for the he equation satisfied by such a cline.

When the fitness of the homozygote AB|AB becomes slightly larger than the one of ab|ab, it is not a priori obvious whether the stationary cline stays stationary or begins to move. Here we answer this question by showing that invasion does occur even if the difference in fitness between homozygotes has a lower order of magnitude (measured by 0 < ε 1) compared to the fitness cost of heterozygotes. By using a rather involved perturbation analysis, we show in Theorem 4.1 that a front traveling at a constant speed c ε > 0 emerges from the stationary cline u 0 solving [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF] when ε becomes positive. Such a traveling front is a solution to the reaction-diffusion equation involving nonlinear gradient terms [START_REF] Fraïsse | The genetics of speciation: are complex incompatibilities easier to evolve?[END_REF]. We give an explicit approximation of the speed c ε which is, from the modelling point of view, the main contribution of the present work. Among other implications, it reveals not only that recombination between the two loci tends to slow down the propagation of the front but also that the stacked clines always travel faster than one cline alone.

The organization of the paper is as follows. In Section 2 we derive the mathematical model, a PDE system involving nonlinear gradient terms. Through a phase plane analysis, we construct stationary solutions in Section 3. Then, in Section 4, we construct traveling fronts thanks to a perturbation argument. In Section 5, a trick enables to derive an explicit approximation of the speed c ε which sheds light on the original model. We conclude and present some perspectives in Section 6.

Derivation of the mathematical model 2.1 Biological assumptions

The population is homogeneous in a one-dimensional space, along which density is supposed uniform and large. We focus on the variations of the genetic composition of the population. We consider that the fitness of a (diploid) individual is affected by two loci: a first locus with two alleles A and a and a second locus with two alleles B and b. We assume that heterozygotes have the lowest fitness (underdominance), the fitness of the different genotypes at each locus being given in Table 1.

genotype fitness AA 1 + 2s A Aa 1 + s A -S A aa 1 genotype fitness BB 1 + 2s B Bb 1 + s B -S B bb 1
Table 1: Fitness of the different genotypes at each locus.

Here the constants s A , s B , S A , S B satisfy 0 ≤ s A < S A , 0 ≤ s B < S B . We then assume multiplicative effects among loci, so that the fitnesses W of two-locus genotypes are given in Table 2.

AA Aa aa BB (1 + 2s B )(1 + 2s A ) (1 + 2s B )(1 + s A -S A ) 1 + 2s B Bb (1 + s B -S B )(1 + 2s A ) (1 + s B -S B )(1 + s A -S A ) 1 + s B -S B bb 1 + 2s A 1 + s A -S A 1 Table 2: Fitness of diploid individuals.
Note that, because we will derive expressions to the first order in s A , s B , S A and S B , assuming additive effects among loci would lead to the same results. We assume that recombination occurs with probability r, so that AB|ab individuals may produce Ab and aB gametes. Throughout the paper, the population occupying the left-hand side of the linear habitat will consist mostly of AB|AB individuals, while the right-hand side will be mostly composed of ab|ab individuals. In the symmetric situation (s A = s B = 0), if the clines of A and B are shifted in space, linkage disequilibrium will develop between the two loci and will pull both fronts together until they are stacked [START_REF] Slatkin | Gene flow and selection in a two-locus system[END_REF], [START_REF] Barton | Multilocus clines[END_REF], as illustrated in Figure 1. This paper is concerned with the established regime where the fronts are stacked (this situation may also result from a secondary contact between two divergent populations, as considered by [START_REF] Barton | Multilocus clines[END_REF]). Note that in the general case (s A , s B = 0), shifted clines may not necessarily become stacked; however, we postpone the analysis of the conditions for stacking to future works.

The PDEs describing the dynamics of two underdominant loci in a 1-dimensional continuous habitat can be obtained by combining the works [START_REF] Barton | The dynamics of hybrid zones[END_REF] and [START_REF] Barton | Multilocus clines[END_REF]. For the selfcontainedness of the present work, we present here a derivation of these equations, obtained by approximating a discrete-time model by a continuous-time model. In Section 2.2 we present the genetic model that drives the genetic dynamics. In Section 2.3 we introduce the spatial structure and the corresponding equations. Finally in Section 2.4, we precise our assumptions on the parameters and their respective magnitudes, as well as our precise objectives.

Recursions on gamete frequencies in a discrete in time setting

We start by considering a single population of diploid, hermaphroditic individuals with nonoverlapping generations. At the end of a generation (at time t), individuals release gametes and immediately die. The next generation, at time t + 1, is formed by the random fusion of gametes. Under these hypotheses, it is sufficient to follow the frequencies of gametes produced at each generation, which completely determine the next generation of individuals (by the law of large numbers). Denote by yA B , yA b , y a B , y a b the frequencies of the different types of gametes at generation t. The fusion at random of these four combinations gives birth to sixteen types of individuals ("ordered" in the sense that z i|j = z j|i for i = j)

z i|j , i, j ∈ { A B , A b , a B , a b },
with proportions p i|j . Notice that, for i = j, the fusion can be male-female or female-male so that we have p i|j = 2 × 1 2 y i y j , thus

p i|j = y i y j .
Each one of these individuals then produces gametes according to its fitness, providing the generation of gametes y A B , y A b , y a B , y a b at time t + 1. Here we assume that there is a probability of recombination 0 ≤ r ≤ 1 2 between the two loci. For each of the sixteen diploid genotypes, the process is as one of the three following examples:

• the individuals zA All these processes are weighted by the fitness of each type of individual, as in the above table. After a tedious but straightforward analysis, one obtains:

B | A B ,
y A B = 1 W (1 + 2s A )(1 + 2s B )y 2 A B + (1 + 2s A )(1 + s B -S B )yA B yA b + (1 + s A -S A )(1 + 2s B )yA B y a B + (1 -r)(1 + s A -S A )(1 + s B -S B )yA B y a b + r(1 + s A -S A )(1 + s B -S B )y a B yA b y A b = 1 W (1 + 2s A )y 2 A b + (1 + 2s A )(1 + s B -S B )yA b yA B + (1 + s A -S A )yA b y a b + (1 -r)(1 + s A -S A )(1 + s B -S B )yA b y a B + r(1 + s A -S A )(1 + s B -S B )yA B y a b y a B = 1 W (1 + 2s B )y 2 a B + (1 + s B -S B )y a B y a b + (1 + s A -S A )(1 + 2s B )y a B yA B (1 -r)(1 + s A -S A )(1 + s B -S B )y a B yA b + r(1 + s A -S A )(1 + s B -S B )yA B y a b y a b = 1 W y 2 a b + (1 + s A -S A )y a b yA b + (1 + s B -S B )y a b y a B + (1 -r)(1 + s A -S A )(1 + s B -S B )y a b yA B + r(1 + s A -S A )(1 + s B -S B )y a B yA b ,
where W is the average fitness:

W = i,j∈{ A B , A b , a B , a B } z i|j W i|j = (1 + 2s A )(1 + 2s B )y 2 A B + (1 + 2s A )y 2 A b + (1 + 2s B )y 2 a B + y 2 a b + 2(1 + 2s A )(1 + s B -S B )yA B yA b + 2(1 + s A -S A )(1 + 2s B )yA B y a B + 2(1 + s A -S A )(1 + s B -S B )yA B y a b + 2(1 + s A -S A )(1 + s B -S B )yA b y a B + 2(1 + s A -S A )yA b y a b + 2(1 + s B -S B )y a B y a b .
Notice that adding the four above equations, one can check

y A B + y A b + y a B + y a b = W W = 1.
For ease of notation, we now let

u := yA B , v := yA b , w := y a B , z := y a b , (1)
so that u + v + w + z = 1. (2)
As in [START_REF] Barton | Multilocus clines[END_REF], we shall rather work on the three components system satisfied by

p := u + v, q := u + w, D := uz -vw, (3)
where

• p measures the frequency of allele A,

• q measures the frequency of allele B,

• D stands for the linkage disequilibrium, measuring the association between alleles A and B within gametes (notice than, equivalently, D = u -pq).

Notice that

u = pq + D, v = p (1 -q) -D, w = (1 -p) q -D, z = (1 -p)(1 -q) + D. (4)
Next, we assume that s A , s B , S A , S B , r are small and of the same order of magnitude, that is

s A ← s A α, s B ← s B α, S A ← S A α, S B ← S B α, r ← rα, ( 5
)
for 0 < α 1. Taking into account (1), ( 2), (3), (4), [START_REF] Barton | Genetic analysis of hybrid zones[END_REF], one can perform straightforward (but tedious) computations and obtain to the first order in α:

               p = p + α (S A (2p -1) + s A ) p(1 -p) + (S B (2q -1) + s B ) D q = q + α (S B (2q -1) + s B ) q(1 -q) + (S A (2p -1) + s A ) D D = D -α r + (2p -1) (S A (2p -1) + s A ) + (2q -1) (S B (2q -1) + s B ) D. ( 6
)

Inserting a spatial structure and switching to continuous time

Finally we consider the associated problem with a spatial structure x ∈ R (corresponding to the position of individuals along space) and continuous time t ≥ 0. More precisely, we assume that gametes migrate according to a dispersal kernel centered on 0 and with variance σ 2 . In the diffusion limit, and from ( 6), the equations for the frequencies p = p(t, x), q = q(t, x) are

   p t = σ 2 2 p xx + (S A (2p -1) + s A ) p(1 -p) + (S B (2q -1) + s B ) D q t = σ 2 2 q xx + (S B (2q -1) + s B ) q(1 -q) + (S A (2p -1) + s A ) D,
where σ > 0. Notice that, since we assumed that the density of individuals is uniform and large, no advection term appears in the above system. As for the equation for the disequilibrium D = uz -vw, we have additional gradient terms (e.g. [START_REF] Barton | Genetic analysis of hybrid zones[END_REF], [START_REF] Barton | Multilocus clines[END_REF]) since

D t = σ 2 2 u xx + • • • z + u σ 2 2 z xx + • • • - σ 2 2 v xx + • • • w -v σ 2 2 w xx + • • • = σ 2 2 (D xx + 2(-u x z x + v x w x )) + • • • = σ 2 2 (D xx + 2(p x q x ) + • • •
where we have used the identity

p x q x = (u + v) x (u + w) x = u x (u x + v x + w x ) + v x w x = -u x z x + v x w x .
Hence, from (6), the equation for D = D(t, x) is

D t = σ 2 2 D xx + σ 2 p x q x -[r + (2p -1) (S A (2p -1) + s A ) + (2q -1) (S B (2q -1) + s B )] D.

Conclusion and goals

Hence the system for the allele frequencies p = p(t, x), q = q(t, x) and the linkage disequilibrium D = D(t, x) is written

           p t = σ 2 2 p xx + (S A (2p -1) + s A ) p(1 -p) + (S B (2q -1) + s B )D q t = σ 2 2 q xx + (S B (2q -1) + s B ) q(1 -q) + (S A (2p -1) + s A )D D t = σ 2 2 D xx + σ 2 p x q x -[r + (2p -1) (S A (2p -1) + s A ) + (2q -1) (S B (2q -1) + s B )] D, (7)
where σ > 0, r > 0, s A > 0, s B > 0, s A > 0 and S B > 0 are given parameters. Observe that, starting from D ≡ 0 (no disequilibrium) the dynamics of p and q are decoupled but the gradient terms p x and q x in the D-equation cause disequilibrium and thus interaction [START_REF] Barton | Multilocus clines[END_REF], see Figure 1. 2 = 1. Left column: p, q and D; the clines are initially uncoupled; next, in a transitory regime, they are driven closer to each other, and eventually become stacked. Right column: the original unknowns, that is the frequencies of gametes u, v, w, z. Remark: the partial differential system in (u, v, w, z) is a reaction-diffusion system for which a standard Strang splitting method was used; numerical simulations were done in python 3.9.2 with the NumPy package version 1.20.1.

Assuming that recombination r is sufficiently large relative to the strength of selection against heterozygotes (S A , S B , determining the gradients in allele frequencies, e.g. [START_REF] Barton | The dynamics of hybrid zones[END_REF]), one expects that D approximately follows

D t ≈ σ 2 2 D xx + σ 2 p x q x -rD.
In the sequel, we use a quasi-linkage equilibrium approximation [START_REF] Barton | Multilocus clines[END_REF], meaning that the dynamics on D is much faster than the one of p and q. As a result,

σ 2 2 D xx + σ 2 p x q x -rD ≈ 0,
whose solution is given by (one may use the Fourier transform to see it)

D(t, x) ≈ σ 2 r ρ σ * (p x q x (t, •))(x), ρ α (x) := 1 2 2r σ 2 e -2r σ 2 |x|
.

For σ sufficiently small, the kernel ρ α "approaches" the Dirac delta function, and thus

D ≈ σ 2 r p x q x . ( 8
)
As a result, using (8) and writing (p, q)(t, x) = (p, q) t,

√ 2
σ x , we reach a simplified version of system [START_REF] Barton | Adaptation, speciation and hybrid zones[END_REF], namely

   pt = pxx + S A f (p) + s A g(p) + 2 r (S B (2q -1) + s B )p x qx , qt = qxx + S B f (q) + s B g(q) + 2 r (S A (2p -1) + s A )p x qx , where f (u) := u(2u -1)(1 -u), g(u) := u(1 -u).
For ease of notation in the mathematical analysis, we now drop the tildes but keep in mind that, when returning to the original model, the traveling waves speeds we will find have to be multiplied by the factor σ √ 2 . Last, we assume that

S A = S B = S, s A = s B = s =: ε, ( 9
)
and thus focus on the system

   p t = p xx + Sf (p) + εg(p) + 2 r (S(2q -1) + ε)p x q x , q t = q xx + Sf (q) + εg(q) + 2 r (S(2p -1) + ε)p x q x . ( 10
)
Notice that f is a balanced bistable nonlinearity, which is slightly unbalanced by the term εg.

In the sequel, our goal is to inquire on the situation where the A cline, measured by p, and the B cline, measured by q, remain stacked together. To do so we look after u = p = q solving the nonlinear equation

u t = u xx + Sf (u) + εg(u) + 2 r (S(2u -1) + ε)u 2 x . ( 11
)
We suspect the existence of a stationary solution connecting 1 to 0 for ε = 0 and that of a front connecting 1 to 0 and traveling at a speed c ε ∼ c 1 ε for some c 1 > 0 and 0 < ε 1. These facts are proved in Section 3 and 4, while c 1 is explicitly identified in Section 5.

Standing together (ε = 0)

In this section, we construct a stationary solution connecting 1 to 0 in [START_REF] Fraïsse | The genetics of speciation: are complex incompatibilities easier to evolve?[END_REF] when ε = 0, and then prove its stability. [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]. In red, the nullcline x = 0, in green the nullcline y = 0, in brown dashed the linear unstable manifold at (1, 0), in blue (an approximation of) the heteroclinic orbit from (1, 0) to (0, 0). Left: the parameters are S = 0.6, r = 0.25 so that ( 14) holds. Right: the parameters are S = 0.85, r = 0.15 so that ( 14) does not hold.

Construction of the standing wave

We are here looking after a u 0 : R → R solving

     u 0 + Sf (u 0 ) + 2 r S(2u 0 -1)(u 0 ) 2 = 0 on R, u 0 (-∞) = 1, u 0 (+∞) = 0. ( 12
)
Lemma 3.1 (A priori estimates). Any standing wave solution of (12) has to satisfy 0 < u 0 < 1 and u 0 (±∞) = 0.

Proof. If u 0 ≤ 1 is not true then, from the boundary conditions, u 0 has to reach a maximum value strictly larger than 1 at some point but, testing the equation at this point, this cannot hold. Hence u 0 ≤ 1 and, from the strong maximum principle, u 0 < 1. Similarly u 0 > 0.

From the equation and the boundary condition, u 0 > 0 in some (A, +∞), so that u 0 is increasing on (A, +∞). As a result u 0 has a limit in +∞, which has to be zero since u 0 is bounded. Similarly u 0 (-∞) = 0.

Using a phase plane analysis (x, y) = (u 0 , u 0 ), the equation in ( 12) is recast

   x = y y = -Sf (x) -2 r S(2x -1)y 2 . ( 13
)
The phase plane analysis is depicted in Figure 2. The equilibria (0, 0) and (1, 0) are saddle points, the eigenvalues of the Jacobian matrix at these points being ± √ S, whereas the equilibrium ( 12 , 0) is a center, the eigenvalues of the Jacobian matrix at this point being ±i S 2 . At equilibrium (1, 0) the linear unstable manifold is the line y = √ S(x -1). To prove the existence of a heteroclinic orbit from (1, 0) to (0, 0), we consider the orbit leaving (1, 0) along the unstable manifold. As long as it has not reached x = 1 2 this trajectory satisfies x < 0 and y < 0 (south west trajectory). In order to prove that the trajectory does cross the vertical line x = 1 2 , we need to construct a barrier, from below, preventing the situation x → l ≥ 1 2 , y → -∞. We choose the line y = α(x -1) with α > 0 to be selected large enough. Choosing α > √ S insures that the trajectory is above the barrier in a neighborhood of (1, 0). We thus need to show that

|y | |x | < α on the points (x, y) such that y = α(x -1), 1 2 ≤ x < 1.
After some straightforward computations, this is recast

ϕ(x) := (2x -1) 1 - 2α 2 r x + 2α 2 r < α 2 S , for all 1 2 ≤ x < 1.
Assuming 1 -2α 2 r < 0, and evaluating the maximum of ϕ on [ 1 2 , 1], we reach

2α 2 r + 1 2 8 2α 2 r -1 < α 2 S ,
which can be obtained with α sufficiently large provided

S < 4r. ( 14
)
Notice that, from the modelling point of view, assumption ( 14) is consistent with the asymptotics "S small" performed in Section 2 (quasi-linkage equilibrium approximation).

On the other hand, even if [START_REF] Slatkin | Gene flow and selection in a two-locus system[END_REF] does not hold, the (right) phase plane analysis of Figure 2 suggests that the heteroclinic orbit joining (1, 0) to (0, 0) still exists, but the above argument does not apply. As a result, under assumption [START_REF] Slatkin | Gene flow and selection in a two-locus system[END_REF], the orbit touches the line x = 1 2 at some point ( 12 , -β) for some β > 0. Since the problem is symmetric with respect to x = 1 2 , we conclude that the orbit then converges to the equilibrium (0, 0) along the stable manifold, the linear stable manifold being given by y = -√ Sx. This trajectory provides a positive and decreasing solution u 0 to [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF].

In other words, we have (nearly) proved the following.

Proposition 3.2 (Stationary solution for ε = 0). Let us assume [START_REF] Slatkin | Gene flow and selection in a two-locus system[END_REF]. Then there is a unique u 0 : R → R solving [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF] and satisfying the normalization condition u 0 (0) = 1 2 . Moreover, u 0 is positive, decreasing, symmetric in the sense that

u 0 (-x) = 1 -u 0 (x) for all x ∈ R,
and has the asymptotics

1 -u 0 (x) ∼ Ce √ Sx as x → -∞, u 0 (x) ∼ Ce - √ Sx as x → +∞, ( 15
)
for some C > 0.

Proof. From the above phase plane analysis, we are already equipped with a positive, decreasing and symmetric u 0 solving [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF]. The asymptotics ( 15) is rather classical but, for the convenience of the reader, we sketch a short and direct proof. We work as x → +∞. We know from the phase plane analysis that u 0 (x) ∼ -√ Su 0 (x) so that

u 0 (x) = e - √ Sx+o(x) . ( 16
)
Now, from the nonlinear ODE, we have, for some K > 0,

-Ku 2 0 (x) ≤ u 0 (x) -Su 0 (x) ≤ Ku 2 0 (x).
Multiplying this by u 0 (x) < 0 and integrating from x to +∞, we have,

- K 3 u 3 0 (x) ≤ - 1 2 (u 0 ) 2 (x) + S 2 u 2 0 (x) ≤ K 3 u 3 0 (x),
so that, for some M > 0,

-M u 2 0 (x) ≤ u 0 (x) + √ Su 0 (x) ≤ M u 2 0 (x). ( 17
)
From this and ( 16) we deduce that e

√ Sx (u 0 (x) + √ Su 0 (x)) = d dx e √ Sx u 0 (x) must be integrable in +∞. As a result there is C ≥ 0 such that e √ Sx u 0 (x) → C as x → +∞. Now the left inequality in (17) implies - √ S ≤ u 0 u 0 + M √ S u 2 0 = u 0 u 0 - M √ S u 0 1 + M √ S u 0 .
Integrating this from 0 to x provides

u 0 (0) 1+ M √ S
u 0 (0) as a positive lower bound for e √ Sx u 0 (x) so that C > 0 and we are done with [START_REF] Volpert | Fredholm theory of elliptic problems in unbounded domains[END_REF]. It remains to prove uniqueness. We use a sliding method argument. Let v 0 be "another" solution such that v 0 (0) = 1 2 . For K ≥ 0, define the shifted function v K (x) := v 0 (x -K). Since v 0 must also have some asymptotics of the form [START_REF] Volpert | Fredholm theory of elliptic problems in unbounded domains[END_REF], say with some constant C > 0 instead of C, we see that u 0 ≤ v K on R for K > 0 sufficiently large. As a result the real number

K 0 := inf {K ∈ R : u 0 (x) ≤ v K (x), ∀x ∈ R}
is well defined and nonnegative. Assume by contradiction that K 0 > 0. Then there is a point x 0 ∈ R where u 0 (x 0 ) = v K 0 (x 0 ) and u 0 (x 0 ) = v K 0 (x 0 ) so that, from Cauchy-Lipschitz theorem, u 0 ≡ v K 0 on R, which is excluded by the normalization conditions. As a result K 0 = 0 and thus u 0 ≤ v 0 . Similarly v 0 ≤ u 0 and we are done.

Stability of the standing wave

We prove here that the standing wave constructed in Proposition 3.2 is linearly stable in the L ∞ norm. More precisely the following holds.

   v t (t, x) = v xx (t, x) + Sf (v(t, x)) + 2 r S(2v(t, x) -1)(v x (t, x)) 2 , t > 0, x ∈ R, v(0, x) = u 0 (x) + εh(x), x ∈ R.
Then there is λ 0 > 0 such that, for any 0 < λ < λ 0 , the following holds: for sufficiently small ε, there is a continuous function γ(ε) satisfying

γ(0) = R h(x)u 0 (x)e 4S r (u 2 0 (x)-u 0 (x)) dx,
and a constant K > 0 such that, for all t > 0,

v(t, •) -u 0 (• + εγ(ε)) C 1 b (R) ≤ Ke -λt .
Proof. We aim at applying a result of Sattinger, namely [13, Theorem 4.1]. To do so, we need to show that the linear operator (obtained by linearizing [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF] around the solution u 0 )

Lh := h + 4S r (2u 0 -1)u 0 h + S f (u 0 ) + 4 r (u 0 ) 2 h,
satisfies the assumptions (i) and (ii) of [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 3.4]. Since equation ( 12) is a scalar quasilinear second-order differential equation set on R and with a smooth nonlinearity, the assumption (ii) of [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 3.4] can be readily checked thanks to [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 5.4]. As for the assumption (i) of [13, Lemma 3.4], we point out that [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Corollary 5.7] does not apply to our situation, and we thus need to determine the spectrum of L.

The liner operator L admits u 0 as principal eigenvector with eigenvalue 0. We remark that L can be written as

Lh = e -2S r (u 2 0 -u 0 ) M he 2S r (u 2 0 -u 0 ) ,
where

M k := k + 2S 2 r (2u 0 -1)f (u 0 ) + Sf (u 0 ) k =: k + Since the weight function e 4S r (u 2 0 -u 0 )
is bounded and uniformly positive, the operators L and M can be considered as acting on the same space C 0 b (R). In particular, λI -L admits a bounded inverse if and only if λI -M does (where I is the identity mapping on C 0 b (R)), and we have (λI

-L) -1 = e -2S r (u 2 0 -u 0 ) (λI -M ) -1 e 2S r (u 2 0 -u 0 ) .
Below, by following ideas of [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF], we analyze, for g ∈ C 0 b (R), the set of solutions to the resolvent equation (λI

-M )k = -k + (λ -c(x))k = g(x), ( 18
)
and then determine the spectrum of M .

System of fundamental solutions to the homogeneous equation:

we first look for a system of fundamental solutions to

-k + (λ -c(x))k = 0, (19)
whose behaviour near ±∞ can be determined (see [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 5.1] for related arguments) for λ ∈ C such that λ + S ∈ R -. Near +∞, this is performed by substituting ϕ 1 (x) = z 1 (x)e -γ + x in (19), where γ + ∈ C solves γ 2 + = λ + S and Re γ + > 0. We obtain

-z 1 + 2γ + z 1 -(S + c(x))z 1 = 0, (20)
which is recast -(z 1 e -2γ + x ) -(S + c(x))z 1 e -2γ + x = 0, so that, assuming z 1 (+∞) = 0,

z 1 (x) = +∞ x e -2γ + (y-x) (S + c(y))z 1 (y)dy, (21)
and thus, assuming z 1 (+∞) = 1,

z 1 (x) = 1 + +∞ x e 2γ + (x-y) -1 2γ + (S + c(y))z 1 (y)dy. ( 22
)
Hence z 1 is written as the solution of a fixed-point problem (22) set on C 0 b (R + ). Notice that the asymptotic behaviour (15) of u 0 implies y → S+c(y) ∈ L 1 (R + ). As a result, for a given x 0 > 0, the right-hand side operator appearing in ( 22 b ([x 0 , +∞)) for x 0 sufficiently large, and z 1 can be extended to (-∞, x 0 ) by solving the adequate Cauchy problem associated with (20). We have therefore constructed a solution ϕ 1 (x) = z 1 (x)e -γ + x to (19) with z 1 ∈ C 0 b (R + ), z 1 (+∞) = 1. By the same procedure, but integrating on [x 0 , x] instead of [x, +∞) in (21), we can a solution ϕ 2 (x) = z 2 (x)e γ + x to (19) with z 2 ∈ C 0 b (R + ) provided by the fixedpoint problem

z 2 (x) = 1 + x x 0 1 -e -2γ + (x-y) 2γ + (S + c(y))z 2 (y)dy.
By the continuous dependence of the fixed-point with respect to the parameter x 0 [17, Proposition 1.2], and by selecting x 0 sufficiently large, z 2 (x) can be made arbitrarily close to 1. Indeed z 2 (x + x 0 ) is the unique fixed point of the operator

T x 0 z(x) := 1 + x 0 1 -e -2γ + (x-y) 2γ + (S + c(x 0 + y))z(y)dy,
and T x 0 converges uniformly to the constant operator T +∞ z ≡ 1 as x 0 → +∞:

T x 0 z -1 C 0 b ([0,+∞)) ≤ 1 2|γ + | +∞ x 0 |S + c(y)|dy z C 0 b ([0,+∞)) ----→ x 0 →∞ 0.
Therefore we have found a system of fundamental solutions (ϕ 1 , ϕ 2 ) to (19) whose behaviour near +∞ is known. We can proceed similarly near -∞ and find another system of fundamental solutions (ψ 1 , ψ 2 ) whose behaviour near -∞ is known.

Summarizing, for each λ ∈ C \ (-∞, -S], we have

ϕ 1 (x) ≈ +∞ e -γ + x , ϕ 2 (x) ≈ +∞ e γ + x , ψ 1 (x) ≈ -∞ e γ + x , ψ 2 (x) ≈ -∞ e -γ + x , (23) ϕ 1 (x) ≈ +∞ e -γ + x , ϕ 2 (x) ≈ +∞ e γ + x , ψ 1 (x) ≈ -∞ e γ + x , ψ 2 (x) ≈ -∞ e -γ + x , ( 24
)
where

A(x) ≈ +∞ B(x) means 0 < lim inf x→+∞ |A(x)| B(x) ≤ lim sup x→+∞ |A(x)| B(x) < +∞.
Notice that, if λ is not an eigenvalue of M , we further know that ϕ 1 is unbounded as x → -∞ (or else it would be an eigenvector), and ψ 1 is unbounded as x → +∞. Notice also that the constants involved in the above estimates are locally uniform in λ.

Solving equation

(18) if λ ∈ C \ (-∞, -S]
is not an eigenvalue of M : from the behaviours near -∞, the functions ϕ 1 and ψ 1 are linearly independent. Therefore, up to redefining ϕ 2 = ψ 1 , we may consider that (ϕ 1 , ϕ 2 ) is a system of fundamental solutions satisfying

ϕ 1 (x) ≈ +∞ e -γ + x , ϕ 2 (x) ≈ +∞ e γ + x , ϕ 1 (x) ≈ -∞ e -γ + x , ϕ 2 (x) ≈ -∞ e γ + x , ϕ 1 (x) ≈ +∞ e -γ + x , ϕ 2 (x) ≈ +∞ e γ + x , ϕ 1 (x) ≈ -∞ e -γ + x , ϕ 2 (x) ≈ -∞ e γ + x .
We use the method of variation of constants to solve (18) and straightforwardly reach

k(x) = C 1 - 1 W x -∞ ϕ 2 (y)g(y)dy ϕ 1 (x) + C 2 - 1 W +∞ x ϕ 1 (y)g(y)dy ϕ 2 (x),
where C 1 and C 2 are arbitrary constants and W is the constant Wronskian W = W (x) = ϕ 1 (x)ϕ 2 (x) -ϕ 1 (x)ϕ 2 (x). Therefore, there is a unique bounded solution k(x), which corresponds to

C 1 = C 2 = 0.
Hence, for each g ∈ C 0 b (R) there exists a unique k ∈ C 2 b (R) such that (λI -M )k = g. By the open mapping theorem, the operator λI -M has a bounded inverse (λI -M ) -1 :

C 0 b (R) → C 2 b (R) → C 0 b (R). In particular, if λ ∈ C \ (-∞, -S]
is not an eigenvalue of M , then λ is in the resolvent set of M. 23), the eigenvector must be proportional to both ϕ 1 and ψ 1 , hence ϕ 1 and ψ 1 are not linearly independent. Hence the Wronskian ϕ 1 ψ 1 -ϕ 1 ψ 1 must vanish. Since the Wronskian is analytic in λ (see [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 5.2]) and not identically zero, the eigenvalues of M in C \ (-∞, -S] are isolated.

The eigenvalues in

C \ (-∞, -S] of M : if λ ∈ C \ (-∞, -S] is an eigenvalue of M from (
Let λ ∈ C \ (-∞, -S] be an eigenvalue of M . Then the associated eigenvector ϕ is a solution to (18) and the former analysis applies. In particular, ϕ and ϕ converge exponentially fast to 0 near ±∞ (at rate ∓γ + , Re γ + > 0) and therefore ϕ ∈ H 1 (R). Since M is symmetric on H 1 (R), we have in fact λ ∈ R. Reproducing the argument of [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Theorem 5.5], we see that there are no positive eigenvalues of M .

We conclude from the above analysis that the eigenvalues of M in C \ (-∞, -S] form a sequence (λ n ) n∈N (with λ 0 = 0) of isolated values in (-S, 0]. As a result the spectrum of M satisfies

σ(L, C 0 b (R)) = σ(M, C 0 b (R)) ⊂ (-∞, -S] ∪ {λ n , n ≥ 0}.
This shows that the assumption (i) of [START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]Lemma 3.4] holds in our case and concludes the proof of Proposition 3.3.

Traveling together (0 < ε 1)

In this section, we construct a traveling front connecting 1 to 0 in [START_REF] Fraïsse | The genetics of speciation: are complex incompatibilities easier to evolve?[END_REF], when 0 < ε 1, through a perturbation argument from the case ε = 0 studied above.

We are here looking after a nonnegative profile u : R → R and a speed c ∈ R solving

   u + cu + Sf (u) + εg(u) + 2 r (S(2u -1) + ε)(u ) 2 = 0 on R, u(-∞) = 1, u(+∞) = 0. (25)
Observe that, from the strong maximum principle we have u > 0. Also, as in the proof of Lemma 3.1, we have u < 1. Hence, we a priori know 0 < u < 1.

We use a perturbation technique and look for u in the form

u = u 0 + h,
where u 0 is provided by Proposition 3.2 and with, typically, h(±∞) = h (±∞) = 0. Plugging this ansatz into the equation, we see that we need F(ε, c, h) = 0, where

F : R × R × E → E is defined by F(ε, c, h) := h + cu 0 + ch + S(f (u 0 + h) -f (u 0 )) + εg(u 0 + h) + 2 r (S(2u 0 + 2h -1) + ε) (u 0 + h ) 2 - 2 r S(2u 0 -1)(u 0 ) 2 . ( 26
)
As for the function spaces, we choose the weighted Hölder spaces

E := C 2,α µ (R), E := C 0,α µ (R), 0 < α < 1, (27)
where, for k ∈ N,

C k,α µ (R) := f ∈ C k (R) : f C k,α µ (R) < +∞ , f C k,α µ (R) := x → e µ √ 1+x 2 f (x) C k,α (R)
, for well-chosen µ ≥ 0. Here, C k,α (R) denotes the Hölder space consisting of functions of the class C k , which are continuous and bounded on the real axis R together with their derivatives of order k, and such that the derivatives of order k satisfy the Hölder condition with the exponent 0 < α < 1 . The norm in this space is the usual Hölder norm.

Our main result in this section then reads as follows.

Theorem 4.1 (Traveling waves for 0 < ε

1). Let 0 ≤ µ < √ S be given. Let F : R × R × C 2,α µ (R) → C 0,α µ (R) be defined as in (26). Then there is ε 0 > 0 such that, for any 0 ≤ ε ≤ ε 0 , there exists (c ε , h ε ) ∈ R × E such that F(ε, c ε , h ε ) = 0. Moreover the map ε → (c ε , h ε ) is continuous, the speed c ε satisfies c ε = - R g(u 0 ) + 2 r (u 0 ) 2 u 0 e 4S r (u 2 0 -u 0 ) R (u 0 ) 2 e 4S r (u 2 0 -u 0 ) ε + o(ε), as ε → 0, ( 28
)
whereas the perturbation profile h ε satisfies R h ε u 0 = 0, for all 0 ≤ ε ≤ ε 0 . ( 29
)
In what follows we aim at applying the Implicit Function Theorem A.1 to the operator F defined in (26), see [START_REF] Apreutesei | Travelling waves for integro-differential equations in population dynamics[END_REF] for a related argument. We straightforwardly compute the derivatives with respect to c and h at the origin (0, 0, 0):

∂ c F(0, 0, 0)(c) = cu 0 ,
and

Lh := ∂ h F(0, 0, 0)(h) = h + 4S r u 0 (2u 0 -1)h + S f (u 0 ) + 4 r (u 0 ) 2 h. ( 30
)
We need to show that ∂ c,h F(0, 0, 0) given by (c, h) → Lh + cu 0 is bijective from and to a well-chosen pair of function spaces. Our strategy is as follows.

In subsection 4.1, thanks to some results of [START_REF] Volpert | Topological degree for elliptic operators in unbounded cylinders[END_REF], [START_REF] Volpert | Fredholm theory of elliptic problems in unbounded domains[END_REF] (recalled in Appendix), we show that L is a Fredholm operator and compute its index (which depends on the choice of µ). Next, in subsection 4.2, we determine the kernel of L. In particular u 0 is the only bounded solution. We also determine the kernel of L * thanks to an algebraic symmetric formulation in a well-chosen weighted L 2 space, from which we deduce the surjectivity of ∂ c,h F(0, 0, 0). Then we conclude the proof of Theorem 4.1 in subsection 4.3.

Fredholm property

Lemma 4.2 (Fredholm property). The operator L :

C 2,α µ (R) → C α µ (R), defined in (30), is Fredholm if µ = √ S and we have ind L = 0 if 0 ≤ µ < √ S, -2 if µ > √ S.
Proof. In view of Remark A.4 it suffices to study the limiting operators (L µ ) ± associated with L µ defined as in (41), namely

(L µ ) ± h = h ∓ 2µh + (µ 2 -S)h,
thanks to Theorem A.3. First since -ξ 2 ∓ 2µiξ + µ 2 -S = 0, corresponding to (39), has no real solution, L is Fredholm. Next, the associated characteristic equation, corresponding to (40), writes

X 2 ± 2µX + (µ 2 -S) = 0,
and has the following roots:

X + 1,2 = -µ ± √ S, X - 1,2 = +µ ± √ S.
If 0 ≤ µ < √ S we deduce that κ + = 1 and κ -= 1 (in the notations of Theorem A.3), hence ind L = 0; if √ S < µ we have κ + = 0 and κ -= 2, hence ind L = -2. This completes the proof of Lemma 4.2.

Kernels of L, L

* and surjectivity of ∂ c,h F(0, 0, 0)

Lemma 4.3 (The kernel of L). Two linearly independent solutions to the linear homogeneous ordinary differential equation

Lh := h + 4S r u 0 (2u 0 -1)h + S f (u 0 ) + 4 r (u 0 ) 2 h = 0 (31)
are given by u 0 and v 0 : x → u 0 (x)

x 0 1 (u 0 ) 2 (z) e -4S r (u 2 0 (z)-u 0 (z)) dz.
Among the two, u 0 is the only bounded solution. As a result, for 0 ≤ µ < √ S, the kernel of the operator L acting on the space

C 2,α µ (R) into C 0,α µ (R) is given by ker L = span u 0 .
Proof. We investigate the solutions h to (31). This is a second-order linear homogeneous ordinary differential equation, and we already know a solution u 0 (as seen by differentiating ( 12)). In this case a second solution v 0 can be sought in the form v 0 (x) = z(x)u 0 (x). Indeed plugging this ansatz into (31) yields the following first order linear ordinary differential equation for z :

z + 2 u 0 u 0 + 4S r (2u 0 -1)u 0 z = 0, or, equivalently, z + ln((u 0 ) 2 ) + 4S r (u 2 0 -u 0 ) z = 0.
As a result, we can select the solution

z (x) = 1 (u 0 ) 2 (x) e -4S r (u 2 0 (x)-u 0 (x)) ,
which we integrate to reach z(x), and thus

v 0 (x) = u 0 (x) x 0 1 (u 0 ) 2 (z) e -4S r (u 2 0 (z)-u 0 (z)) dz. (32)
Now, from the analysis in Section 3, we know that, for some C > 0,

u 0 (z) ∼ Ce - √ Sz , as z → +∞. ( 33
)
Since u 0 (+∞) = 0, the integrand in (32

) is equivalent to 1 C 2 e 2 √
Sz as z → +∞, and thus

v 0 (x) ∼ 1 C2 √ S e √ Sx , as x → +∞. ( 34
)
Thus v 0 is unbounded and, in particular, v 0 ∈ C 2,α µ (R). Since solutions to (31) are the linear combinations of On the other hand, if µ > √ S then

u 0 ∈ C 2,α µ (R) when 0 ≤ µ < √ S, v 0 / ∈ C 2,α µ (R), and since L : C 2,α µ (R) → C α µ (R), we conclude that ker L = span u 0 when 0 ≤ µ < √ S.
ker L * = span u 0 e 4S r (u 2 0 -u 0 ) , v 0 e 4S r (u 2 0 -u 0 ) ,
where v 0 (x) := u 0 (x)

x 0 1 (u 0 ) 2 (z) e -4S r (u 2 0 (z)-u 0 (z)) dz is as in Lemma 4.3.
Proof. Our starting point is to notice that the coefficient of the first-order term in the definition of L, that is u 0 (2u 0 -1), is the derivative of u 2 0 -u 0 so that

Lh = h + 4S r (u 2 0 -u 0 ) h + S f (u 0 ) + 4 r (u 0 ) 2 h,
from which we deduce the formulation

Lh = h e 4S r (u 2 0 -u 0 ) e -4S r (u 2 0 -u 0 ) + S f (u 0 ) + 4 r (u 0 ) 2 h,
which is symmetric in the adequate weighted L 2 space:

R k(Lh)e 4S r (u 2 0 -u 0 ) = - R k h e 4S r (u 2 0 -u 0 ) + R S f (u 0 ) + 4 r (u 0 ) 2 hke 4S r (u 2 0 -u 0 ) = R (Lk)he 4s r (u 2 0 -u 0 ) .
Remark 4.6. We present here an alternate way to prove that u 0 / ∈ rg L remains true when µ ≥ √ S. To do so, let us solve the second-order linear ordinary differential equation

w + 4S r u 0 (2u 0 -1)w + S f (u 0 ) + 4 r (u 0 ) 2 w = u 0 . ( 35
)
Recall that the solutions of the associated homogeneous equation are spanned by u 0 and v 0 provided by Lemma 4.3. To find a particular solution to (35), we use the method of variation of constants. We see that ϕ(x) := λ 1 (x)u 0 (x) + λ 2 (x)v 0 (x) solves (35) as soon as

u 0 λ 1 + v 0 λ 2 = 0 u 0 λ 1 + v 0 λ 2 = u 0 , which yields λ 2 u 0 v 0 -u 0 v 0 u 0 = u 0 , λ 1 = - v 0 u 0 λ 2 .
Since u 0 v 0 -u 0 v 0 is nothing else than the Wronskian, it is equal to θ -1 e -4S r (u 2 0 -u 0 ) for some θ = 0, and thus

   λ 2 (x) = θ(u 0 ) 2 (x)e 4S r (u 2 0 (x)-u 0 (x)) ∼ θC 2 e -2 √ Sx λ 1 (x) = -θv 0 (x)u 0 (x)e 4S r (u 2 0 (x)-u 0 (x)) ∼ -θ 2 √ S ,
where the equivalents are taken as x → +∞ and where we have used (33) and (34). Hence, we can select

     λ 2 (x) = -+∞ x θ(u 0 ) 2 (z)e 4S r (u 2 0 (z)-u 0 (z)) dz ∼ θC 2 2 √ S e -2 √ Sx λ 1 (x) = -x 0 θv 0 (z)u 0 (z)e 4S r (u 2 0 (z)-u 0 (z)) dz ∼ -θ 2 √ S x.
Hence the solutions to (35) are

w(x) = (C 1 + λ 1 (x))u 0 (x) + (C 2 + λ 2 (x))v 0 (x) for any C 1 ∈ R, C 2 ∈ R.
If C 2 = 0 then, from all the above asymptotic, w is unbounded.

If C 2 = 0 then, from all the above asymptotics,

w(x) ∼ - θC 2 √ S xe - √ Sx , as x → +∞.
This above asymptotics shows that w / ∈ C 2,α µ (R) when µ ≥ √ S, and thus u 0 / ∈ rg L.

Construction of traveling waves

We are now in the position to complete the proof of Theorem 4.1, that is the construction of traveling waves for (25) when 0 < ε 1.

Proof of Theorem 4.1. Assume 0 ≤ µ < √ S. Let us recall that F : R × R × C 2,α µ (R) → C 0,α µ (R) is given by (26). It is Fréchet differentiable (even of the class C 1 ) with respect to each of its variables, and we have

∂ ε F(0, 0, 0) = g(u 0 ) + 2 r (u 0 ) 2 , ∂ c F(0, 0, 0) = u 0 L = ∂ h F(0, 0, 0) : h → Lh = h + 4S r u 0 (2u 0 -1)h + S f (u 0 ) + 4 r (u 0 ) 2 h.
We have shown, in Lemma 4.2, that L is a Fredholm operator with indice 0 and, in Lemma 4.3, that the kernel of L is span u 0 in the considered weighted Hölder space.

Our concern is the derivative ∂ c,h F(0, 0, 0) : (c, h) → Lh + cu 0 . It has been shown in Lemma 4.5 that it is surjective. It is not difficult to show that ker ∂ c,h F(0, 0, 0) = {0} × span u 0 , and that the restriction of ∂ c,h F(0, 0, 0) to R × N , where

N := f ∈ C 2,α µ (R) : R f u 0 = 0
is a topological complement of ker L, is injective and still surjective. Therefore we can apply the Implicit Function Theorem A.1 to the restriction of F to R × R × N . We deduce the existence of a branch (c ε , h ε ), 0 ≤ ε 1, of solutions with ε → (c ε , h ε ) continuous and h ε satisfying (29).

It remains to prove (28). Since F is C 1 in all its variables we deduce from F(ε, c ε , h ε ) = 0 and the chain rule that

∂ ε F(ε, c ε , h ε ) + dc ε dε ∂ c F(ε, c ε , h ε ) + ∂ h F(ε, c ε , h ε ) dh ε dε = 0,
which we evaluate at ε = 0 to get and integrating over R, we reach

g(u 0 ) + 2 r (u 0 ) 2 + dc ε dε ε=0 u 0 + L dh ε dε ε=0 = 0.
dc ε dε ε=0 = - R g(u 0 ) + 2 r (u 0 ) 2 u 0 e 4S r (u 2 0 -u 0 ) R (u 0 ) 2 e 4S r (u 2 0 -u 0 ) > 0,
which yields (28) and concludes the proof of Theorem 4.1.

The speed of the traveling stacked clines

In this section, we obtain an explicit form for c 1 = c 1 (r, S) appearing in the asymptotic formula for the speed c ε = c 1 ε + o(ε) as given in (28). This in turn provides valuable insights on the model for coupled underdominant clines.

For convenience let us temporarily denote u = u 0 the standing wave solution constructed in Proposition 3.2. From [START_REF] Gavrilets | Fitness landscapes and the origin of species[END_REF], we see that v := u 2 solves the linear first order

ODE v + 4S r (u 2 -u) v = -2Su (2u -1)u(1 -u),
which is solved as

u 2 (x) = e -4S r (u 2 -u)(x) C -2S x 0 (u 2 -u) (t)e 4S r (u 2 -u)(t) (u -u 2 )(t)dt ,
for some constant C. We integrate by parts and, up to changing the value of the constant C, reach

u 2 (x) = e -4S r (u 2 -u)(x) C - r 2 e 4S r (u 2 -u)(x) (u -u 2 )(x) + x 0 e 4S r (u 2 -u)(t) (u 2 -u) (t)dt = e -4S r (u 2 -u)(x) C - r 2 e 4S r (u 2 -u)(x) (u -u 2 )(x) + r 4S e 4S r (u 2 -u)(x)
.

Letting for instance x → +∞ enforces C = r 2 8S and, returning to the notation u 0 , we finally obtain

u 0 2 (x) = r 2 8S e 4S r (u 0 -u 2 0 )(x) - r 2 (u 0 -u 2 0 )(x) - r 2 8S . ( 36
)
The fact that u 0 can be expressed in terms of u 0 , already observed in [START_REF] Barton | Multilocus clines[END_REF], enables to obtain an explicit form c 1 = c 1 (r, S) appearing in c ε = c 1 ε + o(ε) as given in (28). Indeed, using (36) and recalling that u 0 < 0, we obtain

c 1 = - R r 4S u 0 (x) 1 -e -4S r (u 0 -u 2 0 )(x) dx - R u 0 (x) r 2 8S e 4S r (u 0 -u 2 0 )(x) - r 2 (u 0 -u 2 0 )(x) - r 2 8S 1 2 e -4S r (u 0 -u 2 0 )(x) dx .
Performing the change of variable u = u 0 (x) this is recast

c 1 = 1 0 r 4S 1 -e -4S r (u-u 2 ) du 1 0 r 2 8S e 4S r (u-u 2 ) - r 2 (u -u 2 ) - r 2 8S
1 2

e -4S r (u-u 2 ) du .

Expanding with respect to S r 1, we reach, after a straightforward computation,

c 1 = 1 √ S 1 + 4 15 S r + 2 45
S 2 r 2 + • • • . ( 37
)
In the sequel we denote

c * 1 := 1 √ S 1 + 4 15 S r , (38)
the first order term of expansion (37).

To verify the accuracy of our previsions, we ran simulations of the full system (7), the one established before simplification thanks to the quasi-linkage equilibrium approximation. We numerically estimate the instantaneous speed by following the movement of the center of the fronts. The comparison with the theoretical speed εc * 1 = sc * 1 (let us recall that s appearing in the original model is nothing else than ε, see ( 9)) is shown in Figure 3. We observe that formula (38) gives a very good approximation of the instantaneous speed when S r is not too small (right of the figure). This validates a posteriori the quasi-linkage equilibrium approximation. When r = 0.5 (free recombination), the linkage disequilibrium stays small and the coupled clines move at a speed which is close to the one each cline would have if travelling alone, that is s/ √ S (or sσ/ √ 2S in the original spatial scale, as obtained by Barton [2] in a single-locus model). Indeed, without interaction, we are left with u t = u xx +Sf (u)+sg(u), which is nothing else than the bistable equation (0 < s S < 1)

u t = u xx + Su(1 -u) 2u -1 + s S ,
whose traveling wave, explicitly computed as

1 2 -1 2 tanh √ S 2 (x -s √ S t) , has speed s/ √ S.
At the other extreme, when r = 0 (no recombination) the system becomes equivalent to a single locus where one allele has a fitness advantage 2s and with a cost for heterozygotes 2S, leading to a bistable wave speed of 2s/ √ 2S. When r ∈ (0, 0.5), the speed of the coupled clines decreases monotonously as recombination r increases.

Our concluding remark is as follows: whatever the values of the parameters, interacting and eventually stacked clines travel faster than one cline alone.

Conclusion and perspectives

In this paper we have investigated the solutions of equation [START_REF] Fraïsse | The genetics of speciation: are complex incompatibilities easier to evolve?[END_REF], describing the dynamics of two coupled, asymmetric genetic incompatibilities (underdominant loci) with identical fitness effects, in a quasi-linkage equilibrium regime. The two main results are as follows: first, we have shown that when ε = 0, there is a unique standing wave u 0 under a normalization condition; then, in Section 4, we have shown that when ε > 0 is small enough, there exists a traveling wave u ε defined as a perturbation of u 0 , and we obtained a simple approximation for its speed.

Those results were obtained under a series of assumptions that we recall here for discussion:

s A , s B < S (H1) s A , s B , S r (H2) S A = S B , s A = s B (H3) p = q.
(H4) Assumption (H1) is the frame of this work which was devoted to the heterozygote inferior case. It is therefore not a hypothesis we want to discuss per se.

Assumption (H2) expresses that we are in the case of small selective advantages. When it does not hold, D may not be small, in which case the quasi-linkage equilibrium approximation (that allowed us to reduce the number of variables) is no longer valid. It can easily be seen that -1 4 ≤ D ≤ 1 4 always holds, and that, as shown by the D equation in [START_REF] Barton | Adaptation, speciation and hybrid zones[END_REF], positive D is generated whenever p and q travel in the same direction (that is p x q x > 0), while negative D is generated otherwise. These facts help to understand the kind of contribution D makes to the coupling between p and q in [START_REF] Barton | Adaptation, speciation and hybrid zones[END_REF].

Assumption (H3) is basically a hypothesis of exchangeability between loci. Although this allowed us to simplify the algebra, different loci should have different fitness effects, and it would thus be of interest to relax this hypothesis.

Last but not least, assumption (H4) conveys the strong argument that the A cline and the B cline were stacked in the past and will remain stacked forever in the future. This is indeed a good starting point from a mathematical perspective. Nevertheless, in the context of population genetics, more interesting questions arise when (H4) does not hold. In such a situation, the coupling in [START_REF] Barton | Adaptation, speciation and hybrid zones[END_REF] can give rise to non-standard behaviours, such as adaptation of the speed. The questions that arise are such as: can a traveling front be pinned by a standing front? Will a front traveling at a large speed crossing a slower traveling front adapt its speed so as to remain stacked with the slower one? A preliminary numerical exploration has shown that there can be a vast zoology of situations. We hope to present them in a future work.

A Some useful results and tools

We recall the Implicit Function Theorem, see [START_REF] Zeidler | Nonlinear functional analysis and its applications[END_REF]Theorem 4.B] for instance.

Theorem A.1 (Implicit Function Theorem). Let X, Y and Z be three Banach spaces. Suppose that:

(i) The mapping F : U ⊂ X × Y → Z is defined on an open neighbourhood U of (x 0 , y 0 ) ∈ X × Y and F(x 0 , y 0 ) = 0.

(ii) The partial Fréchet derivative of F with respect to y exists on U and F y (x 0 , y 0 ) : Y → Z is bijective.

(iii) F and F y are continuous at (x 0 , y 0 ).

Then, the following properties hold:

(a) Existence and uniqueness. There exist r 0 > 0 and r > 0 such that, for every x ∈ X satisfying x -x 0 ≤ r 0 , there exists a unique y(x) ∈ Y such that y -y 0 ≤ r and F(x, y(x)) = 0.

(b) Continuity. If F is continuous in a neighbourhood of (x 0 , y 0 ), then the mapping x → y(x) is continuous in a neighbourhood of x 0 .

(c) Higher regularity. If F is of the class C m , 1 ≤ m ≤ ∞, on a neighbourhood of (x 0 , y 0 ), then x → y(x) is also of the class C m in a neighbourhood of x 0 .

In Section 4 we apply Theorem A.1 to the operator F defined in (26), with X = R, x = ε, x 0 = 0, Y = R × C 2,α µ (R), y = (c, h), y 0 = (0, 0), and Z = C 0,α µ (R).

Next, we quote some results on Fredholm operators. Let us recall that the operator L has the Fredholm property with index 0 if ker L has a finite dimension, rg L is closed and has finite codimension and ind L := dim ker L -codim rg L = 0.

In particular, since its range is closed, such an operator is normally solvable: Finally, let us define the limiting operators

L ± u := a ± u + b ± u + c ± u,
and assume that for any λ ≥ 0, the equation

L ± u -λu = 0
has no nontrivial solution in C 2,α (R).

Then L is Fredholm with index 0.

Let us also recall a Fredholm property result for second-order ordinary differential equations, see the monograph of Volpert [15, Chapter 9, Theorem 2.4 p. 366].

Theorem A.3 (Fredholm property for second-order ODEs). With the notations of Theorem A.2, the operator L is Fredholm provided the two equations

-a ± ξ 2 + b ± iξ + c ± = 0 ( 39
)
has no real solution ξ ∈ R. In this case the index of L is given by the formula

ind L = κ + -κ -,
where κ ± is the number of complex solutions to the characteristic equation

a ± X 2 -b ± X + c ± = 0 (40)
which have a positive real part.

Figure 1 :

 1 Figure 1: Numerical solutions with parameters s A = s B = 0 (symmetric case), S A = S B = r = 0.1 and σ 22 = 1. Left column: p, q and D; the clines are initially uncoupled; next, in a transitory regime, they are driven closer to each other, and eventually become stacked. Right column: the original unknowns, that is the frequencies of gametes u, v, w, z. Remark: the partial differential system in (u, v, w, z) is a reaction-diffusion system for which a standard Strang splitting method was used; numerical simulations were done in python 3.9.2 with the NumPy package version 1.20.1.

Figure 2 :

 2 Figure2: Phase plane analysis for[START_REF] Sattinger | On the stability of waves of nonlinear parabolic systems[END_REF]. In red, the nullcline x = 0, in green the nullcline y = 0, in brown dashed the linear unstable manifold at (1, 0), in blue (an approximation of) the heteroclinic orbit from (1, 0) to (0, 0). Left: the parameters are S = 0.6, r = 0.25 so that (14) holds. Right: the parameters are S = 0.85, r = 0.15 so that (14) does not hold.

Proposition 3 . 3 (

 33 Stability of standing waves). Let u 0 be the standing wave constructed in Proposition 3.2. Let h ∈ C 1 b (R) be given. Let v solve the parabolic Cauchy problem

  ) is globally Lipschitz continuous on C 0 b ([x 0 , +∞)) with Lipschitz constant 1 2|γ + | +∞ x 0 |S + c(y)|dy. Hence, equation (22) has a unique solution z 1 on C 0

Lemma 4 . 4 (r (u 2 0

 442 The kernel of L * ). If 0 ≤ µ < √ S then the kernel of the adjoint operator L * is ker L * = span u 0 e 4S -u 0 ) .

  Since rg L = (ker L * ) ⊥ = span u 0 e

Figure 3 :

 3 Figure 3: Comparison of the theoretical speed sc * 1 and of the numerically estimated speed of the stacked fronts.

∃u = 0 ,

 0 Lu = f ⇔ ∀φ ∈ (rg L) ⊥ , φ(f ) = 0, and remark that (rg L) ⊥ = ker L * . We recall below a theorem from Volpert, Volpert and Collet [16, Theorem 2.1 and Remark p787].Theorem A.2 (Fredholm property on the line). For 0 < α < 1, consider the operatorL : C 2,α (R) → C α (R) defined by Lu := a(x)u + b(x)u + c(x)u,where the coefficients a(x), b(x), c(x) are smooth, and a(x) ≥ a 0 for some a 0 > 0. Assume further that the coefficients a(x), b(x), and c(x) have finite limits as x → ±∞ and denote a ± := lim x→±∞ a(x), b ± := lim x→±∞ b(x), c ± := lim x→±∞ c(x).

whose proportion is y 2

	A B , whose proportion is yA , release gametes A B in proportion 1. B yA , release gametes A B in proportion 1 2 b b in proportion 1 • the individuals zA B | A b and gametes A 2 . • the individuals zA B | a b , whose proportion is yA y a b , release gametes A B and a b both in B proportion 1-r 2 (no recombination), and gametes A b and a B both in proportion r 2 (recombi-
	nation).

In particular, for any k ∈ C 2,α µ (R), we have

Then we deduce from Lemma 4.2 and Lemma 4.3 that dim ker L * = -ind L + dim ker L = 0 + 1 = 1, and therefore we do have ker

This time, the asymptotics for v 0 being given in (34), terms Lemma 4.5 (Surjectivity of ∂ c,h F(0, 0, 0)). Let 0 ≤ µ < √ S be given. Then, the application

Proof. We check that u 0 is not in the range of L.

Since rg L has codimension 1 by Lemma 4.2 and 4.3, we have C 0,α µ (R) = rg L ⊕ span u 0 . This shows that ∂ c,h F(0, 0, 0) is surjective. Remark A.4 (Fredholm property in weighted Hölder spaces). We cannot directly apply Theorem A.2 and Theorem A.3 to our situation since we consider the operator L acting from C 2,α µ (R) into C α µ (R), and not from C 2,α (R) into C 0,α (R). To circumvent this, we consider the operator L µ : C 2,α (R) → C α (R) defined by:

µ (R) is continuously invertible, the map L = T -1 µ L µ T µ shares the same Fredholm property and index as L µ . As a result, if L µ satisfies the assumptions of Theorem A.2, or Theorem A.3, then L is a Fredholm operator with the same index as that of L µ .