Segmenting documents into lines and words is a very critical step before the recognition task. It is even more difficult for ancient and calligraphic writings, as is often the case in Arabic manuscript documents. In this work, we propose a new attempt to segment documents into lines and words, using deep learning. For line segmentation, we use an RU-net which allows a pixel-wise classification, thus separating pixels of lines from the background pixels. For segmenting lines into words, not having a ground truth for the word segmentation (at the image level), we use the line transcription to find the words. A BLSTM-CTC is used to achieve this mapping directly between the transcription and line image, without any segmentation. A CNN precedes this sequence to extract the features and feeds the BLSTM with the essential of the line image. Tested on KHATT Arabic database, the system achieves good performance that is of the order of 96.7% correct lines and 80.1% correct words.

Introduction

The processing of handwritten documents for their recognition remains a big challenge. The cursive and ever-changing nature of handwritten lines makes segmentation into lines and then into words usually necessary steps. Among the difficulties encountered in this type of document, we find the irregular inclination of the lines, the occurrence of ascenders and descenders connecting the lines, and the change of style of writing according to the writer. The writing of words in several subwords introduces confusion in the separation of words between them, and finally the presence of diacritics adds an additional difficulty.

The works related to line segmentation can be divided in three categories: 1) top-down methods, 2) bottom-up methods and 3) machine learning methods. The two first methods generally require prior knowledge on the documents, such as line inter-spaces, document orientation, etc. Therefore, systems must combine many kinds of image processing methods to take into consideration all possible features.

In a particular way, top-down methods process the whole image and iteratively subdivide it into smaller blocks to isolate the desired part [START_REF] Papavassiliou | Handwritten document image segmentation into text lines and words[END_REF]. These methods include image processing techniques like projection profile, Hough transform or Gaussian filters.

Bottom-up methods, meanwhile, deal with noise problems and writing variation [START_REF] Belaïd | In: Segmentation of Ancient Arabic Documents[END_REF]. The connected components (CCs) based methods are the most popular methods used in this approach [START_REF] Boulid | Arabic handwritten text line extraction using connected component analysis from a multi agent perspective[END_REF][START_REF] Rajput | Text-line extraction from handwritten document images using histogram and connected component analysis[END_REF]. These methods generally use simple rules, analyzing the geometric relationships between neighboring blocks such as the distance or the overlap.

Conversely, machine learning methods treats the image as a whole without any prior knowledge. Some of the existing systems are end-to-end systems, with no further post or pre-processing [START_REF] Renton | Handwritten text line segmentation using fully convolutional network[END_REF][START_REF] Chen | Convolutional neural networks for page segmentation of historical document images[END_REF], while others use deep learning like one step among other processing [START_REF] Grüning | A two-stage method for text line detection in historical documents[END_REF][START_REF] Pastor-Pellicer | Complete system for text line extraction using convolutional neural networks and watershed transform[END_REF].

Often, line and word segmentation tasks are separated because of the difference between line and word features, so global processing for both tasks can be tough. Also, we noticed that deep learning is much less used in word segmentation than on line segmentation. The methods used in word segmentation are rather bottom-up methods, based on CC analysis [START_REF] Alkhateeb | Component-based segmentation of words from handwritten arabic text[END_REF][START_REF] Louloudis | Text line and word segmentation of handwritten documents[END_REF], structural feature extraction [START_REF] Aouadi | Word extraction and recognition in arabic handwritten text[END_REF] or even both [START_REF] Elzobi | Off-line handwritten arabic words segmentation based on structural features and connected components analysis[END_REF]. These approaches show interesting results on documents written in Latin or Germanic languages.

We propose a line segmentation system using an RU-net, and an end-to-end system for word segmentation, using a CNN (Convolutional Neural Network) followed by a BLSTM (bidirectional Long Short Term Memory) and finally a CTC function (Connectionist Temporal Classification) [START_REF] Graves | Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks[END_REF] that will automatically learn the alignment between text line images and the words in the transcription.

The paper is organized as follows: Section 0.1 gives an overview of some related works. Section 0.3 presents the proposed approaches for line and word segmentation. Section 0.4 reports the experimental results and section 0.5 concludes the work.

Related Work

In this section, line segmentation and word segmentation previous works will be presented.

Line Segmentation

As an example of bottom-up methods, Cohen & al. [START_REF] Cohen | Using scale-space anisotropic smoothing for text line extraction in historical documents[END_REF] starts their system by applying a multi-scale multi-scale anisotropic second derivative of Gaussian filter bank to enhance text regions, then apply a linear approximation to merge connected components from the same line using a K N N function. The authors proposed a generic system that may be used for any language, and obtained interesting results on Hebrew with 98% for text region detection.

As an example of top-down methods, the authors in [START_REF] Arvanitopoulos | Seam carving for text line extraction on color and grayscale historical manuscripts[END_REF] use a seam carving method. They first compute medial seams on the lines, using a projection profile matching approach, then compute separating seam using a modified version of seam carving procedure [START_REF] Avidan | Seam carving for content-aware image resizing[END_REF]. Their method showed good results on the Arabic dataset they used (99.9%), but coarse results on our dataset (see Figure 1).

Among the several works using deep learning for line segmentation, Renton & al. [START_REF] Renton | Handwritten text line segmentation using fully convolutional network[END_REF] use a modified version of the FCN (Fully Convolutional Network) [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF]. In [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF], the authors use FCN for semantic segmentation and introduce "skip" steps between first and last layers to avoid the coarse results due to many pooling layers, but the authors in [START_REF] Renton | Handwritten text line segmentation using fully convolutional network[END_REF] argue that this is not sufficient for the task of line segmentation which needs to be more precise. Thus, Renton & al. [START_REF] Renton | Handwritten text line segmentation using fully convolutional network[END_REF] propose a new architecture where the layer convolution + max pooling is replaced by dilated convolutions. Their network has been trained on x-height labeling, reaching up to 75% of F-measure on the cBad dataset [START_REF] Diem | ScriptNet: ICDAR 2017 Competition on Baseline Detection in Archival Documents (cBAD)[END_REF].

Grüning & al. [START_REF] Grüning | A two-stage method for text line detection in historical documents[END_REF] use a two stages method to extract lines from document images. The first stage is an ARU-net which is a U-net extended with an attention model (A) and a residual structure (R). The U-net is a variant of FCN [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] that introduces shortcuts between layers of the same spatial dimension, so we can have features from higher level and reduce the vanishing gradient problem. The authors add residual blocks [START_REF] He | Deep residual learning for image recognition[END_REF] to the U-net and an attention model that works simultaneously with the RU-net at multi-scale level. The output of the network is two maps, one related to the detected baselines and another to the beginning and end of the lines. In the second stage, they calculate a set of superpixels from the baseline map, then they estimate the state of these superpixels by computing their orientation and interline distance. Finally, using a Delaunay neighborhood system and calculating some functions like the projection profile, the data energy and the data cost, they find superpixels clusters that represent the lines. Their system reaches 95% of good detection of lines on the whole cBad dataset, but the many processing and computations used, especially on the second stage, makes it a heavy procedure.

Word Segmentation

Papavassiliou & al. [START_REF] Papavassiliou | Handwritten document image segmentation into text lines and words[END_REF] calculate an SVM-based gap metric for adjacent CCs within each text line, then apply a threshold to classify gaps as "within" and "between" words. They tested their approach on the datasets from the IC-DAR07 Handwriting Segmentation Contest [21] and reached an F-measure of 93%. Al Khateeb & al. [START_REF] Alkhateeb | Component-based segmentation of words from handwritten arabic text[END_REF] looked at the Arabic word segmentation problem with a component-based method where they analyze CCs with the help of baselines, and reach 85% of correct rate segmentation. In [START_REF] Al-Dmour | Word extraction from arabic handwritten documents based on statistical measures[END_REF], Al-Dmour & al. propose a method based on two spatial measures: CC length and gaps between them. Lengths are clustered to separate between the groups of letters, subwords and words. Gaps are clustered to figure out whether the gap occurs "between-words" or "within-a word". These measures are clustered using a Self-Organizing Map (SOM) algorithm [START_REF] Kohonen | The self-organizing map[END_REF]. The approach has been tested on the AHDB dataset [START_REF] Al-Ma'adeed | A data base for arabic handwritten text recognition research[END_REF] and achieved 86.3% of correct extraction rate.

Proposed Approach 0.3.1 Global Presentation

Inspired by [START_REF] Grüning | A two-stage method for text line detection in historical documents[END_REF], we choose to use an RU-net for our line segmentation, on an x-height labeling, followed by a simple post-processing that extracts baselines.

There are many representations in the literature for what is a text line in the line segmentation, like baselines or bounding boxes. The x-height represents the area corresponding to the core of the text without ascenders and descenders and seems to be a suitable text line representation for text recognition, while the baseline represents the main orientation of a text line and is mainly used for the evaluation. An example of x-height and baseline labeling is provided in the Figure 2. Figure 3 shows the structure of our system. The output of the RU-net is a prediction map with the same size of the input. The network is then followed by a post-processing step that extracts the baselines from the x-heights.

For word segmentation task, we chose a CNN followed by a BLSTM (forward and backward) and finally a CTC decoder. Figure 4 depicts the structure of the network.

Detailed Description of the Approach

Line Segmentation:

Our goal here is first to extract the x-heights of text lines, by providing a ground truth that separates 3 classes (background, paragraphs and lines' x-heights in each paragraph) as shown in Figure 8.b. From then on, our problem becomes similar to a semantic segmentation, which can be solved by an FCN [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF].

The FCN structure can be described as a symmetric encoder/decoder network where the encoder part convolves and downsamples the input with a number of convolution+pooling blocks, then the decoder part upsamples the result by the same factor of downsampling. So the output will be a map with the same size of the input where each pixel has been assigned to a probability of belonging to a class. We chose an U-net as FCN model, extended with residual connections. The U-net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] allows an easier combination of low level features and high level ones by introducing shortcuts between the same level blocks and the residual connections greatly reduce the vanishing gradient problem [START_REF] He | Deep residual learning for image recognition[END_REF].

We provide, as an input, paragraph images resized to 720 pixels for the largest side and keep the same ratio between height and width. The images are resized so we can feed the network smaller images, thus reduce the memory footprint. The encoder part of our network is composed of 3 convolution + maxpooling blocks and a 4 th convolution layer, and the decoder part of 3 convolution + unpooling blocks.

This network is followed by a simple post-processing step on the output of the network (a map of 3 classes): A dilation is performed on the x-heights, so the over-segmented lines can be merged, then an erosion, to recover the original shape of x-heights. After that, polynomial regression is applied on the connected components (x-heights here), so the baselines can be computed. Although this post processing is well adapted to our case, it would not be for lines written on several columns, for example.

Word Segmentation:

The network is fed by line images extracted from the previous network. All images have a normalized size of 48×1600. To deal with the images, we chose to use a CNN, that will extract the most important features of the input. Every convolutional block is followed by a batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] that greatly reduces the vanishing gradient problem and makes the use of dropout unnecessary [START_REF] Li | Understanding the disharmony between dropout and batch normalization by variance shift[END_REF]. The output of the CNN is then passed sequentially to a bidirectional LSTM (100 neurons for each LSTM) which is followed by a CTC function [START_REF] Graves | Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks[END_REF]. Let's consider that we have two classes (word: 1, space: 0). The ground truth provided is then the Unicode transcription of the text lines where each word is labeled by 1 and each space by 0, as shown in Figure 5, and the CTC decoder's output will be the sequence word-space found. No study has proven which one is better [START_REF] Chung | Empirical evaluation of gated recurrent neural networks on sequence modeling[END_REF], and it all depends on one's case, so the 2 networks have been tested for comparison. Figure 6 shows the difference between LSTM's and GRU's accuracy through epochs, and LSTMs showed the better results.

On the other hand, the use of a BLSTM instead of a simple one may provide additional context for the network, so the learning will be faster and better.

CTC function:

As there is no obvious alignment between the network's output and the ground truth, the CTC function is used to perform it automatically. Graves & al. [START_REF] Graves | Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks[END_REF] initially introduced the CTC to address the problem of alignment in speech recognition. Given the output of an RNN (a sequence of probabilities), the CTC loss function computes the probability of an alignment per each time-step using a dynamic programming algorithm. A beam search decoder is then used to extract the top paths, i.e. the paths with the greatest probability values.

Results & Discussion

Line Segmentation:

From the KHATT database, we labeled 50 paragraph images following the pattern shown in Figure 8.b, 35 for the training set and 15 for the testing set. For all our experiments, a learning rate of 10 -4 and the Adam Optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] are used. The network has been trained for 100 epochs. Table 1 summarizes the parameter settings of the line segmentation architecture.

Image pre-processings: 50 paragraph images: 35 for training and 15 for testing. Images are normalized to a size of 720px for the largest side. No further processing.

RU-net:

See Figure 3 for the detailed number of filters, kernels size and pool size. Strides : 2. Dropout : 0.5 (after every convolutional layer).

Training settings:

Initial weights: 1.0 for the 3 classes. Initial learning rate: 10 -4 . Optimizer : Adam. Initial number of epochs: 100. Table 1: Parameter settings of the line segmentation architecture. The 4 nearest graph modeling further connected these prices to form a structure.

The metric used for the evaluation of the baseline detection is presented in [START_REF] Grüning | Read-bad: A new dataset and evaluation scheme for baseline detection in archival documents[END_REF] and computes the F-measure on the extracted baselines. tested on KHATT dataset. The system proposed in [START_REF] Grüning | A two-stage method for text line detection in historical documents[END_REF] has not been tested, since the code for the second stage is not provided. The results in [START_REF] Renton | Handwritten text line segmentation using fully convolutional network[END_REF] are evaluated on the x-heights, while those in [START_REF] Grüning | A two-stage method for text line detection in historical documents[END_REF] and RU-net are evaluated on the baselines. Figure 8 shows the results obtained by the used system. The results obtained by the used system are encouraging, especially since it is quite simple and requires less parameters and processing steps. Furthermore, we trained it on a small set, a bigger training set may provide better results, but having to label the images by hand could be the brake of this system.

Word Segmentation:

For our experiments, we used 5000 line images, 4800 for the training set and 200 for the testing set. A learning rate of 10 -4 with an Adam optimizer are also used. Table 3 summarizes the parameter settings of the word segmentation architecture. Figure 9 shows an example of the results obtained by our network and Table 4 summarizes the F-measure values of our approach and other works. Since the prediction is done on a downscaled image, the output needs to be upscaled to fit the original image. The F-measure was calculated manually by comparing the words on the transcription and the resulting segmentation. The most common segmentation errors are misplaced segmentations, like shown in Figure 10 where the segmentation is between sub-words.

Compared to other works on the Arabic word segmentation, ours gives less good measure but doesn't use any shape specification of the language, which is probably the reason of these results. Furthermore, the databases are not the same and KHATT may provide more difficult writings.

Conclusion

In this paper, line and word segmentation approaches based on deep learning models has been proposed. The line segmentation system uses an RU-net to extract x-heights from text images, then a post-processing step extracts baselines. The word segmentation approach uses a CNN with a BLSTM, then a CTC to find the alignment between the line transcription and the line image. The results of the two approaches are promising, but still need to be improved, especially the word segmentation approach. A future work would be to add some post-processing on the connected components for example, to correct the wrong segmentation.

Figure 1 :

 1 Figure 1: Results of [15]: a) on their dataset(Al-Majid), b) on our dataset.

Figure 2 :

 2 Figure 2: X-height (green) and baseline (red) of a text line.

Figure 3 :

 3 Figure 3: Structure of the RU-net. Blue boxes denote multi-channel feature maps. The number of channels is provided on the top of each box while sizes are provided on the bottom-left of the boxes.

Figure 4 :

 4 Figure 4: Proposed word segmentation network.

Figure 5 :

 5 Figure 5: Ground truth provided to the CTC (1: word; 0: space).

Figure 6 :

 6 Figure 6: Obtained results using LSTM and GRU.

Figure 7 :

 7 Figure 7: Example of a text line and it's transcription in the KHATT database.

Figure 8 :

 8 Figure 8: Results obtained for the line segmentation: a) the original image, b) the ground truth composed of three classes (background: red; paragraph: blue; x-height: green), c) the output of the RU-net, d) the final result after post-processing.

  Image pre-processing: 5000 line images: 4800 for training and 200 for testing. Images are normalized to a size of 48 × 1600. No further processing. CNN+BLSTM+CTC See Figure 4 for detailed description of the CNN used, strides: 1, no dropout. The BLSTM has 100 neurons for each direction. See Figure 5 for an example of the ground truth provided to the CTC function. Beam size: 100. Training settings: Weights initializer: Xavier. Initial learning rate: 10 -4 . Optimizer: Adam. Number of epochs: 100.

Figure 9 :

 9 Figure 9: Results of the proposed approach. The lower image is the output of the BLSTM for the class space (0), the 200 probabilities has been upscaled ×8 to recover the original shape. The upper image is the segmented line after the projection of the probabilities.

Figure 10 :

 10 Figure 10: Example of wrong segmentations. The red box denotes one word and an over-segmentation is made on the sub-words.

Table 2

 2 below gives the results of the machine learning based methods for line extraction,

	Architecture Dataset F-measure on line extraction
	[5]	KHATT cBad	81% 75%
	[7]	cBad	95%
	RU-Net(ours)	KHATT	96.71%

Table 2 :

 2 Results on line extraction.

Table 3 :

 3 Parameter settings of the word segmentation architecture.

	Architecture	Dataset	F-measure
	[9]	IFN/ENIT	85%
	[22]	AHDB	86.3%
	CNN+BLSTM+CTC	KHATT	80.1%

Table 4 :

 4 Results of our approach.


Keywords: . hal-01912694 Navier-Stokes equation, isokinetic curves, -time, bounded polynomial iterations, Weierstrass-Mandelbrot functions

Isokinetic solutions of the Navier-Stokes equation.

solution, at such a point that this existence constitutes one of the millennium's problems posed by Fefferman. But all the methods used to solve the question are too much dependant of the linearization of the equation. An oversee of the question can be found in the Lemarie-Rieusset's book. Fefferman [6], writer of the problem of the equation of Navier for the Clay's foundation, seems very embarrassed to formulate it: in fact, he proposes four problems A, B, C, D. Carefully, he seems accept periodical solutions to represent the physically reasonable movement. As we will see, we are not sure it is true. On other hand, he exposes very clearly many ambiguities of the question. But, his paper is singularly imprisoned by many physical considerations, with many parameters and derivatives in relation with the non linearity of the equations. The Navier's equations concern the unknown velocity vector and the pressure of an incompressible fluid in movement. The variables are the position x and the time t: x=(x, ); and . At time , the velocity vector is function of (x, ):

-The millennium presentation

We have as many equations as unknown functions (n+1): , with the null divergence:

and the known function: at the initial time: .

-Vector notation

-Let For , we note:

And the vectors:

Then:

So, the NSE is: ∂ -Let And: ∂

-Main hypothesis

In order to simplify the problem, we suppose that the fluid is located in a bounded open domain such as its frontier is empty of fluid: so, all quantities concerning the fluid are null on . This frontier is smooth and don't move with the time. That means : at the time 0 and during all the time , as is open, all the fluid is far from the frontier : for this finite time , especially at time ; the pressure is null on and verifies the Dirichlet's conditions. We can have problems with the non autonomous externally forces So, we suppose free of : . is at least of class. If reach , the problem change completely.

-Isokinetic curve

Let at time 0, the velocity is known for the position . From this starting point at time 0, we suppose that the NSE gives a velocity at time .

Definition

An isokinetic curve is such as: . We search in this paper all the isokinetic curves except near the frontier . This approach gives some new reflexions about the conditions and the understanding of the existence of the solutions of the Navier-Stokes equation (NSE). Note that the problem of the isokinetic curves is different from the wellknown stationary solutions where ∂ . (see Leray).

To find the isokinetic curves, we decompose this equation in two problems when the external forces don't depend on the time:

-The first problem is the existence of the initial pressure defined by a Poisson's equation when the velocity is known; this question seems well-known under Dirichlet's conditions, especially if the frontier is smooth and bounded. -The second is the resolution of an ordinary differential equation defining the position of the velocity in function of the time. This equation is completely defined by initial conditions. The Dirichlet's conditions imply that the domain is bounded and relatively compact. The technics of iteration with the Borh's theorem give conditions to obtain an interesting solution.

-Computation of the pressure

First, we recall: 1 -Lemma 1 If is a known function of ( , the condition implies that the pressure follows the Poisson's equation: , With the Dirichlet's conditions. Then, the pressure is the unique solution in the space H( of the infinitely derivable functions with compact support in .

■ The divergence operator commutes with all derivations. If we apply this operator to the first n equations defining ∂ , we obtain: ;

With: =∂

=∂

And:

So:

The pressure is only function of and of . If is known, the unicity of depends on the unicity of the solution of the Poisson's equation:

But, the solution of the Poisson's equation is given modulo a function satisfying the Laplace's equation. Limit conditions such as the Dirichlet's conditions are necessary to prove the existence and the uniqueness of this solution. This condition doesn't appear clearly in the text of Fefferman, but many authors, as Tao, know it. ■ The Dirichlet's condition implies that is bounded.

-Lemma 2

If, at any arbitrary fixed time 0, the velocity is a ( known function of the position : and if frontier is smooth and don't move with the time, then, the pressure is the unique solution in the space H( , which is the adherence of the infinitely derivable functions with compact support in . If the velocity is a smooth known ( function of the position and smooth on , then, the pressure is ( . ■ We obtain the same Poisson's equation with the same limit conditions: So, at time 0, the existence and the uniqueness of the solution is valid. ■ 3 -Isokinetic curves in a small neighbourhood of 0 First, we present some properties of the isokinetic curves if they exist.

Lemma 3 A differentiable isokinetic curve can exist only under the condition:

Then, we have:

■ If is a differentiable isokinetic curve, we have = .
That means:

. Then an unique isokinetic curve can exist only if is invertible. So, we must have . ■

-Hypothesis H0 at a fixed time =0:

is a known smooth function of .

-The jacobian matrix is invertible. If , that means:

, the domain of is (Dirichlet) bounded under the previous condition for . We note H if H is true at each time between and .

-Proposition 1

If the isokinetic solution of the NSE for exists under H , for all , , and this solution is: where is the known initial velocity, then must be the solution of the ordinary differential equation with the known function of :  . We suppose that is locally lipchitzian with respect to . Then, the Cauchy's problem:

)= for each
with the initial condition , has unique solution in a neighbourhood of . If is of class , is and maximal.

-Remark

If , the corresponding surface is a singularity for the NSE and we have risk of turbulences all along this surface. So, before to do anything, we have to note the points such as , then, we study the behaviour of the process far from these points as possible. So, we define the time:

). But, we note that remains bounded. We will clarify this condition in the next proposition.

-Proposition 3

Let be the known smooth function of the position at time . In a neighbourhood of 0 :

, under H , the unique maximal solution of the NSE, is for each initial condition :

+ is solution of the ONSE with known for an initial position : ); -Then: =

Where: With : = More, it exists a local flow in a neighbourhood of . ■ -We take 0 and write under H for all , If Is solution of the NSE, must be solution of:

Where: And is the unique solution = If we note the initial position at time 0, the position at time will be . ■

-General solution far from the turbulences and the frontiers

As the solution of the ONSE defining is bounded, we have proved in our book the following result: -The solution of the ONSE is reduced to the study of the zero of ) ( and the eigen values at these points. More, if we suppose that can be uniformly approximated by polynomials, at least locally near fixed points ( , especially if they are hyperbolic, the solution is a power series of Weirestrass -Mandelbrot functions. The equality of Perceval is satisfied.

■ Hint: If the space of the ONSE is bounded, the image is relatively compact and we can use the Bohr -Cauchy's representation with almost periodic functions. (see my little book). Then, under the hypothesis of the hyperbolicity of the fixed points, the solution is developed in power series of Weierstrass -Mandelbrot's functions. But, here the fixed points and the eigen values at these fixed points are determined by the function Anyway, we have seen that the solutions are local near each fixed points and we can meet many other perturbations when the process cross from a basin of a fixed point to an another as in the Lorenz iteration. ■ But, we can have many other turbulences when .

-Construction of near by the singularities: the turbulences

We stay far from the frontiers Exercise 1 In , if is not invertible, we have with .

-1 -1

That means belongs in . Then, in , the solution exists. So we work in with .

-Computation when

If , the corresponding surface is a singularity for the NSE and we have turbulences all along this surface. So, before to do anything, we have to note the points such as: . In this case, for all :

,we have at least a vector such as and . (We don't examine here the other situations). So, , and then, in a neighbourhood of 0:

-Construction in a neighbourhood of

Now, study the behaviour of the process near these points. We need this result:

Exercise 2 Let and a -lipchitzian function of , with independent of , then, for all , the Cauchy,s problem of with >0has a unique solution. ■ With the change of variable , the proof is achieved. We note that is a point of singularity. More, at this point , is not lipchitzian. ■ Notation . We denote , and we can write with a coordinate of and all the other coordinates of . is smooth and we apply the implicit theorem . To simplify the problem, we choose such as will be invertible in a very large domain.

We write the ONSE with the comatrix of :

That means:

So, except at time where , the ONSE is completely modified.

-Computation in a neighbourhood of

It is easy to check the coherence of the solution , at every time as we have done for all . But we have to consider all the turbulences where at each time .

-Conclusion

The paradigm of the problem has completely moved: the equations of the movement depend entirely on the given initial situation . But the good solution is local near by the initial when there is no turbulence. Even in this case, the situation is complicate because the solution is developed in power series of Weierstrass -Mandelbrot's functions and not with periodic functions. The mapping of the space with isokinetic curves seem useful to study the NSE, but we have to study an infinity of equations with possible turbulences instead of the NSE. We note the importance of the isobaric curves, well-known in meteorology, more easy to operate than the velocity because it is a one dimensional parameter. The analysis of the isokinetic curves are more complicated when the problem is not autonomous: that is if . But, the method can be used to integrate in the equation more parameters like the temperature…

New generalised equation

We see easily that we can take more complicated derivable function . The method will be the same. In any functions (smooth functions, derivatives, even unique solution of an equation as the pressure of .

Definition

  initial condition: The pressure verifies the Poisson equation at : If does not depend of , the isokinetic curves are isobaric curves. ■ At time the pressure must satisfy the Poisson's equation for fixed and the known function at : = If has a good smooth frontier (as under Dirichlet's conditions), the solution is unique. And is a known function of . Now, we study solution of the NSE for and NSE is reduced to an ordinary differential equation (ONSE) near with the initial position . Suppose that be the solution of ONSE for the initial position at time . The velocity at time for the known initial position is: this solution exists. As the velocity is constant along an isokinetic curve: = ; then the pressure obeys the same Poisson's equation at time and at time 0: The unicity of the solution implies that the pressure remains constant. ■ function defined on OI, where O is an open set of and I an open set of .  . We recall well-known Cauchy-Lipchitz results for the ODE: 1 -Proposition 2 Let a continuous function defined on OI, where O is an open set of and I an open set of .

  Let be a smooth function of, and of (any smooth functions, derivatives, even unique solution of a functional or integro-differential equation of) isokinetic curves as curves such as: Theorem Then, the isokinetic curves solution of:Is, for all : , where is the solution of the ordinary differential equation with function of)■ Hints: As is function of every smooth functions, derivatives and unique solutions of a functional or integro-differential equation of can be explained as functions of and we can write . The only important things to check are: -the jacobian , -the unicity of the solution of -the veracity of H at each time t, and especially the point, -the border of the domain of . And we have to solve the ordinary differential equation: ■


Hélène Leman
email: helene.leman@polytechnique.edu

H Leman Cimat

A stochastic model for reproductive isolation under asymmetrical mating preferences

Keywords: mating preference, asymmetrical preference, birth-death stochastic model, dynamical system, long-time behaviour, perturbation method Mathematics Subject Classification (2000) 92D40, 37N25, 60J27

More and more evidence shows that mating preference is a mechanism that may lead to a reproductive isolation event. In this paper, a haploid population living on two patches linked by migration is considered. Individuals are ecologically and demographically neutral on the space and differ only on a trait, a or A, affecting both mating success and migration rate. The special feature of this paper is to assume that the strengths of the mating preference and the migration depend on the trait carried. Indeed, patterns of mating preferences are generally asymmetrical between the subspecies of a population. I prove that mating preference interacting with frequency-dependent migration behavior can lead to a reproductive isolation. Then, I describe the time before reproductive isolation occurs. To reach this result, I use an original method to study the limiting dynamical system, analyzing first the system without migration and adding migration with a perturbation method. Finally, I study how the time before reproductive isolation is influenced by the parameters of migration and of mating preferences, highlighting that large migration rates tend to favor types with weak mating preferences.

Introduction

Understanding the mechanisms of speciation and reproductive isolation is a major issue in evolutionary biology. There is now strong evidence that sexual preferences and speciation are tied [START_REF] Lande | Models of speciation by sexual selection on polygenic traits[END_REF][START_REF] Boughman | Divergent sexual selection enhances reproductive isolation in sticklebacks[END_REF]. Initially, the role of 'magic' or 'multiple effect' traits, which associate both adaptation to an ecological niche and a mate preference, has been studied deeply. It has been shown that such traits may lead to speciation, using direct experimental evidence [START_REF] Merrill | Disruptive ecological selection on a mating cue[END_REF] or theoretical works [START_REF] Lande | Ecological speciation by sexual selection[END_REF][START_REF] Van Doorn | Sympatric speciation and extinction driven by environment dependent sexual selection[END_REF]. Then, studies focused on the particular role of mating preference during a speciation event [START_REF] Gavrilets | Models of speciation: Where are we now?[END_REF], highlighting that (i) it may impede reproductive isolation [START_REF] Servedio | Limits to the evolution of assortative mating by female choice under restricted gene flow[END_REF][START_REF] Servedio | Limits to the evolution of assortative mating by female choice under restricted gene flow[END_REF]Bürger, 2014, 2015), or, (ii) it may promote reproductive isolation. This promoting role may be secondary or primary. For example, the initial divergence in traits may be the result of natural selection in order to decrease hybridization and then be subjected to mating preference [START_REF] Panhuis | Sexual selection and speciation[END_REF], producing speciation by reinforcement [START_REF] Gregorius | Extinction times in multitype markov branching processes[END_REF]. Other studies illustrate the direct and promoting role of assortative mating, using numerical simulations [START_REF] Kondrashov | On the origin of species by means of assortative mating[END_REF][START_REF] 'gonigle | Sexual selection enables long-term coexistence despite ecological equivalence[END_REF][START_REF] Smadi | Looking for the right mate in diploid species: How does genetic dominance affect the spatial differentiation of a sexual trait?[END_REF], or theoretical results [START_REF] Rudnicki | Model of phenotypic evolution in hermaphroditic populations[END_REF][START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF].

The studies mentioned above focus on a symmetrical sexual preference, assuming that all individuals express the same sexual preference. Numerous observations and studies though do not support this assumption and describe examples of species that express different patterns of preference (See [START_REF] Panhuis | Sexual selection and speciation[END_REF] for examples). [START_REF] Smadja | Asymmetrical reproductive character displacement in the house mouse[END_REF] describe such an example between two subspecies of the house mouse. The subspecies Mus musculus musculus is characterized by a stronger assortative preference than the subspecies Mus musculus domesticus [START_REF] Smadja | Strong premating divergence in a unimodal hybrid zone between two subspecies of the house mouse[END_REF]. A mechanism for subspecies recognition mediated by urinary signals occurs between these two taxa and seems to maintain reproductive isolation. Another example comes from Drosophila melanogaster populations where strong sexual isolation with an asymmetrical pattern of sexual preference was observed [START_REF] Wu | Sexual isolation in drosophila melanogaster: a possible case of incipient speciation[END_REF][START_REF] Hollocher | Incipient speciation by sexual isolation in drosophila melanogaster: variation in mating preference and correlation between sexes[END_REF]. The Zimbabwe female lines of Drosophila melanogaster have a nearly exclusive preference for males from the same locality over the males from other regions or continents; the reciprocal mating is also reduced but to a lesser degree.

In this paper, I was thus interested in asymmetrical patterns of preference. I focused on the cases where mating preference promotes sexual selection, and I was interested in two main problematic: (i) studying the influence of an asymmetric mating preference pattern on speciation mechanisms, and (ii) understanding the effects of migration on mating preference advantages. To do so, I aimed to generalize the model of [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] to account for asymmetrical sexual preferences. A haploid population divided in two demes but connected by migration is considered. Following the seminal papers of [START_REF] Bolker | Using moment equations to understand stochastically driven spatial pattern formation in ecological systems[END_REF]; [START_REF] Dieckmann | Relaxation projections and the method of moments. The Geometry of Ecological Interactions: Symplifying Spatial Complexity[END_REF]; [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF], I used a stochastic individual-based model with competition and varying population size. Individuals are assumed not to express any local adaptation. Their parameters do not depend on their location. Individuals, however, are characterized by a mating trait, encoded by a bi-allelic locus, and which has two consequences: (i) individuals of the same type have a higher probability to mate and give an offspring, and (ii) the migration rate of an individual increases with the proportion of individuals caring the other trait in its deme. Finally, the two alleles may not have identical effects, in the sense that strengths of mating preferences and of migration depend on the allele carried by the individual.

Using convergence to the large population limit, I first connected the microscopic model to a macroscopic and deterministic model. Then, studying both models, I established the main result of the paper, which ensures that the mechanism of mating preference combined with a negative typedependent migration is sufficient to entail reproductive isolation and which gives the time needed before reproductive isolation. Here, unlike [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF], the approximation of this time depends on both mating preference parameters and both migration parameters, related to both alleles. I finally conducted an extensive study on the influence of migration and preference parameters on this time showing that large migration rates can favor types with a weak mating preference. The proof of the main result is based on a fine analysis of the deterministic limiting model. In particular, global results on the dynamical system are established such that the dynamics of almost all trajectories can be predicted. To do so, I developed an original method based on a perturbation theory of the migration parameters, which strongly differs from the method used by [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. The asset of this method is that it can easily be adapted to other dynamical systems.

The paper is organized as follows. In Section 2, the stochastic model is introduced and motivated from a biological perspective. Section 3 presents the results of the paper. In particular, the main results on the deterministic limiting model and on the stochastic process are stated in Section 3.1. Section 3.2 presents the main result in the case without migration between both patches. In Section 3.3, the influence of migration on the time before reproductive isolation is analyzed. Section 4 establishes the proof of the key result using perturbation theory. Finally, Section 5 gives a discussion about the paper results. Proofs of the case without migration will be found in Appendix A. Proofs of the probabilistic parts of the main result will be found in Appendix B.

Model

The population is divided into two patches. The individuals are haploid and characterized by a diallelic locus (a or A) and a position (1 or 2 depending on the patch in which they are). The set E := {(α, i), α ∈ {a, A}, i ∈ {1, 2}} is used to characterize the individuals. The population dynamics follows a multi-type birth and death process with competition in continuous time. In other words, the dynamics follows a Markov jump process in space N E , where the rates are described below. At any time t, the population is represented by the following vector of dimension 4 in N E :

N K (t) = (N K A,1 (t), N K a,1 (t), N K A,2 (t), N K a,2 (t)) ∈ N E
where N K α,i (t) denotes the number of individuals with genotype α in the deme i at time t. K is an integer parameter associated to the concept of carrying capacity and accounting for the quantity of available resources or space (see also [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] for more details). Consequently, it is a scaling parameter for the size of the community. It is assumed to give the order of magnitude of the initial population, in the sense that the initial number of individuals divided by K converges (in probability) when K goes to infinity. The competition for resources is also scaled with 1/K, as presented below.

In what follows, if α denotes one of the alleles, notation ᾱ denotes the other allele and if i denotes one of the demes, ī denotes the other one.

The birth, death and migration rates of each individual are now described. At a rate B > 0, a given individual with trait α ∈ {a, A} encounters uniformly at random another individual of its deme. Then it mates with the latter and transmits its trait with probability bβ α /B ≤ 1 if the other individual carries also the trait α, and with probability b/B ≤ 1 if the other individual carries the trait ᾱ. That is to say, after encountering, two individuals that carry the same trait α have a probability β α -times larger to mate and give birth to a viable offspring than two mating individuals with different traits. Hence, as the current state of the population is denoted by

N K ∈ N E , the total birth rate of α-individuals in patch i is bN K α,i β α N K α,i + N K ᾱ,i N K α,i + N K ᾱ,i .
(2.1)

Note that two parameters, β a and β A , are used to model the sexual preference depending on the trait carried by the individual. The limiting case where β A = β a was studied by [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. Here, I was interested in the case where β a = β A although the result of the limiting case can be rediscovered with our calculation. As presented in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF], Formula (2.1) models an assortative mating by phenotypic matching or recognition alleles (Blaustein, 1983;[START_REF] Jones | Mate choice and sexual selection: What have we learned since darwin?[END_REF]. Note that, in the present model, preference modifies the rate of mating and not only the distribution of genotypes, unlike what is usually assumed in classical generational models (O 'Donald, 1960;[START_REF] Lande | Models of speciation by sexual selection on polygenic traits[END_REF][START_REF] Kirkpatrick | Sexual selection and the evolution of female choice[END_REF][START_REF] Gavrilets | Fitness landscapes and the origin of species (MPB-41)[END_REF][START_REF] Bürger | Intraspecific competitive divergence and convergence under assortative mating[END_REF][START_REF] Servedio | Limits to the evolution of assortative mating by female choice under restricted gene flow[END_REF]. The present model can be compared with these classical ones by computing the probabilities that an individual of trait α in the deme i gives birth after encountering an individual of the same trait (resp. of the opposite trait) conditionally on the fact that this individual gives birth at time t, and we find

β α N K α,i β α N K α,i + N K ᾱ,i resp. N K ᾱ,i β α N K α,i + N K ᾱ,i
.

Note that these terms correspond to the ones presented in the supplementary material of [START_REF] Servedio | Limits to the evolution of assortative mating by female choice under restricted gene flow[END_REF], or in [START_REF] Gavrilets | On the evolution of premating isolation after a founder event[END_REF]. A extended discussion between these two types of models (discrete time versus continuous time models) can be found in Section 2 of [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF].

The death rate of a given individual is composed of a natural death rate and a competition death rate. Individuals compete for resources or space against all individuals of its own deme. The competitive death rate of each individual is thus proportional to the total population size of its deme. Finally, the total death rate of α-individuals in patch i is

d + c K N K α,i + N K ᾱ,i N K α,i , (2.2)
where d models the natural death and c models the competition for resources. As presented previously, K is the scaling parameter that scales the amount for resources. Hence, the larger K is, the smaller the strength of competition between two individuals, c/K, is. Finally, individuals can migrate from one patch to the other one. Following [START_REF] Payne | Sexual selection, space, and speciation[END_REF][START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF][START_REF] Smadi | Looking for the right mate in diploid species: How does genetic dominance affect the spatial differentiation of a sexual trait?[END_REF], I use a frequency-dependent migration rate in such a way that individuals are more prone to move if they do not find a suitable mate. This hypothesis is relevant for all organisms with active mate searching [START_REF] Uy | Complex mate searching in the satin bowerbird ptilonorhynchus violaceus[END_REF][START_REF] Jugovic | Movement, demography and behaviour of a highly mobile species: A case study of the black-veined white, aporia crataegi (lepidoptera: Pieridae)[END_REF]. The migration term of an individual is proportional to the proportion of individuals carrying the other allele in its deme, and to a parameter m α which depends on the trait of the individual. Hence, the alleles code for the strength of the mating preference and simultaneously, the speed of migration. The total migration rate of α-individuals from patch 1 to patch 2 finally is

m α N K ᾱ,1 N K α,1 + N K ᾱ,1 N K α,1 . (2.3)
Note that the migration rate does not depend on the other deme composition.

In what follows, the following statements on the parameters are assumed:

β A > 1, β a > 1, b > d > 0, c > 0, m A ≥ 0, m a ≥ 0.
3 Results

Time needed before reproductive isolation

In this section, I present the main result of the paper that gives the time needed for the process N K to reach reproductive isolation. This time depends on K, the carrying capacity of the process.

To this aim, let us first give the average behavior of the process using convergence to the large limit population. Precisely, Lemma 1 below ensures that the sequence of re-scaled processes

(Z K (t)) t≥0 = N K (t) K t≥0 converges when K goes to infinity to                          d dt z A,1 (t) = z A,1 b β A z A,1 + z a,1 z A,1 + z a,1 -d -c(z A,1 + z a,1 ) -m A z a,1 z A,1 + z a,1 + m A z A,2 z a,2 z A,2 + z a,2 d dt z a,1 (t) = z a,1 b β a z a,1 + z A,1 z A,1 + z a,1 -d -c(z A,1 + z a,1 ) -m a z A,1 z A,1 + z a,1 + m a z A,2 z a,2 z A,2 + z a,2 d dt z A,2 (t) = z A,2 b β A z A,2 + z a,2 z A,2 + z a,2 -d -c(z A,2 + z a,2 ) -m A z a,2 z A,2 + z a,2 + m A z A,1 z a,1 z A,1 + z a,1 d dt z a,2 (t) = z a,2 b β a z a,2 + z A,2 z A,2 + z a,2 -d -c(z A,2 + z a,2 ) -m a z A,2 z A,2 + z a,2 + m a z A,1 z a,1 z A,1 + z a,1 . (3.1)
Lemma 1 Assume that the sequence (Z K (0)) K≥0 converges in probability to the deterministic vector

z 0 ∈ R E . Then, for any T ≥ 0, lim K→∞ sup s≤T Z K (s) -z (z0) (s) = 0 in probability, (3.2)
where . denotes the L ∞ -Norm on R E and (z (z0) (t)) t≥0 denotes the solution of (3.1) with initial condition z 0 ∈ R E This result can be deduced from a direct application of Theorem 2.1 p. 456 by [START_REF] Ethier | Markov processes: Characterization and convergence[END_REF].

A direct computation implies that the following four points are stable equilibria of the system (3.1):

equilibria with fixation of an allele (where only an allele is maintained in both patches)

(ζ A , 0, ζ A , 0), (0, ζ a , 0, ζ a ), (3.3)
equilibria with maintenance of each allele in a different patch

(ζ A , 0, 0, ζ a ), (0, ζ a , ζ A , 0), (3.4) with ζ α := bβα-d c
, α ∈ {A, a}. These four equilibria describe states of reproductive isolation: Once reaching one of these equilibria, migration rates equal zero and individuals do not migrate anymore. More specifically, observe that equilibria (3.4) are of particular interest to our problematic. Indeed, once reaching one of these equilibria, even if a small basal migration (i.e. constant migration) is added, the mating preferences and the frequency-dependent migration terms will prevent the populations of both demes to mix again, leading to migration-selection balance [START_REF] Karlin | Polymorphisms for genetic and ecological systems with weak coupling[END_REF] but where selection is due to sexual selection and not to natural selection. Precisely, if an A-individual travels because of basal migration from patch 1 to patch 2, which is filled with a-individuals, its probability to reproduce will be significantly reduced in patch 2 and its migration rate to come back will be so high that it is quite unlikely that its offspring establish in patch 2. This reasoning, however, fails with equilibria (3.3).

Our aim is then to understand the long-time behavior of trajectories of the dynamical system and more specifically to detail the set of initial states that lead to one of these equilibria, this corresponds exactly to the basin of attraction of an equilibrium. With this aim in mind, let us define the weighted sums

Σ i := (β A -1)z A,i + (β a -1)z a,i , for i = 1, 2, Σ := Σ 1 + Σ 2 = (β A -1)(z A,1 + z A,2 ) + (β a -1)(z a,1 + z a,2 ),
and the compact set

S := z ∈ R E , Σ i ≥ (β min -1)(b -d) 4c , for i = 1, 2, and Σ ≤ 4(β max -1)(bβ max -d) c , (3.5)
where β min = min(β a , β A ) and β max = max(β a , β A ). Next Lemma ensures that we can focus on trajectories starting from S. Indeed any trajectory reaches it in finite time.

Lemma 2 Assume that

m A (β A -1) + m a (β a -1) ≤ 2(b -d)(β min -1). (3.6)
S is a positively invariant set for the dynamical system (3.1), in the sense that any trajectories starting from this set will stay in it when t grows to +∞. Moreover, any trajectory solution of (3.1) hits S after a finite time.

The aim is thus to study trajectories in the compact set S.

Theorem 1 1. Assume that m A = 0 if and only if m a = 0. There exists m 0 > 0 such that for all m A ≤ m 0 and m a ≤ m 0 , there exist four open subsets (D α,α m A ,ma ) α,α ∈A of S that are the basins of attraction in S of the four equilibria (3.3) and (3.4) of the system (3.1), and such that the closure of ∪ α,α ∈A D α,α m A ,ma is equal to S.

2. In the case m A = m a = 0, the basins of attraction are exactly

D α,α 0,0 = z ∈ R E , (β α -1)z α,1 > (β ᾱ -1)z ᾱ,1 and (β α -1)z α ,2 > (β ᾱ -1)z ᾱ ,2 ∩ S.
Theorem 1 ensures that any trajectory starting from S, except from a set with empty interior which corresponds to the boundary of the basins of attraction, reaches one of the steady states (3.3) or (3.4). In particular, coexistence of both alleles in a single deme does never occur. In other words, assortative mating combined with negative type-dependent migration entails reproductive isolation.

The assumption on the migration rate is essential to obtain this result. Different results are deduced in models with frequency-independent migration [START_REF] Servedio | The counterintuitive role of sexual selection in species maintenance and speciation[END_REF][START_REF] Smadi | Looking for the right mate in diploid species: How does genetic dominance affect the spatial differentiation of a sexual trait?[END_REF].

In particular, reproductive isolation may be prevented. Smadi et al ( 2018) study a similar model as the one used here, but individuals are diploid. A mechanism of mating preference interacting with frequency-dependent migration is studied. In Section 3.4 of this paper, the frequency-dependent migration term is replaced by a constant migration term. Then, polymorphic equilibria with both alleles in demes can be observed, but only if the migration rate is sufficiently large. This highlights that, although using other kind of migration prevents reproductive isolation, the mechanism that would prevent reproductive isolation is migration and not assortative mating, in their case as in the one presented here. Theorem 1 is, furthermore, a key result to deduce the next theorem, which gives the time before reproductive isolation. It can be compared to Theorem 2 of [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] which gives same results in the symmetrical case (β A = β a and m A = m a ). In the latter, the equilibrium reached is given by the alleles that make up the majority initially in each patch. In our case, the dynamics is more involved. Without migration, the equilibrium reached depends on the initial number of individuals of each type and of the mating preference strengths. Then, when m A and m a are small, the basin of attraction D α,α m A ,ma is a continuous deformation of D α,α 0,0 . I drew such an example in Section 3.3. Note that no basin of attraction is empty, since the four equilibria are stable equilibria.

The asymmetrical sexual preferences make the long-time behavior more involved than in the symmetrical case and proofs here use completely different mathematical techniques. I used perturbation theory to deduce Theorem 1 : the system is first studied in the particular case where m A = m a = 0, then one makes m A and m a grow up to deduce the result for positive migration rates. Unfortunately, I was not able to give an explicit formulation for the sets D α,α m A ,ma unlike in the symmetrical case.

Let us now state the main result. It describes the random time T K B A,a,ε that is the first time when the population process N K reaches the set

B A,a,ε := [(ζ A -ε)K, (ζ A + ε)K] × {0} × {0} × [(ζ a -ε)K, (ζ a + ε)K],
with ε > 0 and when K is large. In other words, it is the random time before (1) all a-individuals in patch 1 and all A-individuals in patch 2 get extinct, and (2) the population size in patch 1 is approximately Kζ A and the one in patch 2 is approximately Kζ a . In the light of the previous discussion about equilibrium (ζ A , 0, 0, ζ a ), it thus corresponds to the time before reproductive isolation occurs.

Theorem 2 Assume that Assumptions of Theorem 1 holds and that m A ≤ m 0 and m a ≤ m 0 . Let ε 0 > 0 and assume also that Z K (0) = N K (0)/K converges in probability to a deterministic vector z 0 ∈ D A,a m A ,ma such that (z 0 a,1 , z 0 A,2 ) = (0, 0). Then there exist C 0 > 0, M > 0, and V > 0 depending only on (M, ε 0 ) such that, for any ε ≤ ε 0 ,

lim K→∞ P T K B A,a,ε log K - 1 ω(A, a) ≤ C 0 ε, N K T K B A,a,ε + t ∈ B A,a,M ε ; ∀t ≤ e V K = 1, (3.7)
where for all α, α ∈ A,

ω(α, α ) = 1 2 b(β α -1 + β α -1) + m α + m α - b(β α -β α ) + (m α -m α ) 2 + 4m α m α .
(3.8) Similar results hold for the three other equilibria of (3.3) and (3.4).

Theorem 2 gives the first-order approximation of the time before reproductive isolation. The latter is proportional to log(K), which is short compared to K, the order of magnitude of the population size. Comparatively, the time scale needed for random genetic drift to cause the end of gene flow between two populations is of order K in many models (as for example in Wright-Fisher model, see chapter 7 in Hartl et al (1997)). Hence, Theorem 2 implies that reproductive isolation due to mating preference is much shorter. Note also that Theorem 2 not only gives the time before reproductive isolation but also it ensures that once the equilibrium is reached, the population sizes of both patches stay around Kζ α during at least a long time of order e KV . Secondly, as ω(α, α) = b(β α -1), the time before reaching one of equilibria (3.3) does not depend on migration parameters unlike the time before reaching one of equilibria (3.4). I studied more specifically the influence of migration parameters on this time in Section 3.3.

The assumption on initial condition ((z 0 a,1 , z 0 A,2 ) = (0, 0)) is only needed to obtain the lower bound on the time T K B A,a,ε given in (3.7). Otherwise, this time would be faster. Finally, note that, assuming β A = β a and m A = m a , Theorem 3 of Coron et al ( 2018) is rediscovered.

Study of the system without migration

The proofs of Theorems 1 and 2 require a full understanding of the dynamics without migration. Hence before proceeding with the proofs, I present a complete study of the dynamical system when m A = m a = 0. Since both patches evolve independently in this case, only the dynamics of patch 1 is studied and, for the sake of simplicity, the dependency on patches in notation is dropped. From (3.1), we find that

       d dt z A (t) = z A b β A z A + z a z A + z a -d -c(z A + z a ) d dt z a (t) = z a b β a z a + z A z A + z a -d -c(z A + z a ) (3.9)
The equilibria of the system will be written with the following quantities

χ α := (β ᾱ -1)χ, where χ := b(β a -1)(β A -1) + (b -d)(β A -1 + β a -1) c(β A -1 + β a -1) 2 ,
and where ᾱ is the complement of α ∈ A. A direct computation implies that there exist exactly four fixed points of the dynamical system (3.9):

(0, 0), (ζ A , 0), (0, ζ a ), and (χ A , χ a ).

These equilibria represent respectively the extinction of the population, the loss of allele a or allele A, or the long-time coexistence of both alleles.

Let us now describe their stability and the long time behavior of any solution of (3.9).

Lemma 3 -(ζ A , 0) and (0, ζ a ) are stable equilibria, (0, 0) is unstable and (χ A , χ a ) is a saddle point.

-The set

D A 0 := (z A , z a ) ∈ R 2 , (β A -1)z A > (β a -1)z a (3.10)
is a positively invariant set under the dynamical system (3.9). Moreover, any solution starting from D A 0 converges to (ζ A , 0) when t converges to +∞.

-The set

D a 0 := {(z A , z a ) ∈ R 2 , (β A -1)z A < (β a -1)z a }
, is a positively invariant set under the dynamical system (3.9). Any solution starting from D a 0 converges to (0, ζ a ) when t converges to +∞.

-Finally, {(z A , z a ) ∈ R, (β A -1)z A = (β a -1)z a } is also a positively invariant set and any solution starting from this set converges to (χ A , χ a ) when t grows to +∞.

In other words, the system is bi-stable: All trajectories converge to (ζ A , 0) or (0, ζ a ), except the trajectories starting from a line (see Fig. 1). A direct consequence of this Lemma is that the basin of attraction D α,α 0,0 are exactly the ones described by Theorem 1.

Influence of parameters on the time before reproductive isolation

In this section, the model under study is the initial one with two demes. I used functional studies and simulations to explore the influence of migration rates and mating preference parameters on the process. The simulations below were computed with the following demographic parameters: Assume that the process starts from a state z 0 ∈ D A,a m A ,ma . Then, according to Theorems 1 and 2, the trajectory will reach a neighborhood B A,a,ε of equilibrium (ζ A , 0, 0, ζ a ) after a time of magnitude log(K)ω(A, a) -1 . Direct functional studies ensure that the constant of interest, ω(A, a) -1 , is a decreasing function with respect to β A and to β a whatever the other parameters are (see Fig. 2,left). Hence, the stronger the sexual preference is, the faster the reproductive isolation is.

β A = 2, β a = 1
Then, I focus on how the constant depends on m A and m a . It may be natural to consider that m A and m a can be rewritten using three positive parameters γ A , γ a and m as follows:

m A := γ A m and m a := γ a m.
In this way, both migration parameters change simultaneously with m. Once again, a direct functional study ensures that ω(A, a) -1 is a non-increasing function with respect to m (see Fig. 2, right). Hence, increasing both migration rates at the same time accelerates the reproductive isolation, in the same way as when mating preference parameters increase. Moreover, the migration parameters used in the model are frequency-dependent terms such that individuals are more prone to migrate when they do not find suitable mates in their deme. With this in mind, the first conclusion is that a large migration rate seems to strengthen the homogamy.

The result is then improved by studying how constant ω(A, a) -1 changes with respect to m A and m a separately. A direct computation shows that ω(A, a) -1 is a decreasing function with respect to m A if β A > β a and it is an increasing function with respect to m A if β A < β a . In other words, if A-individuals have a stronger sexual preference than a-individuals (β A > β a ), the bigger their migration rate is when they are in contact with too much a-individuals, the shorter the time before reaching the equilibrium is. Once again, it highlights that the effects of migration and sexual preference are similar. However, assuming again that A-individuals have a stronger sexual preference than a-individuals (β A > β a ), the bigger the a-individuals migration rate is, the longer the time before reproductive isolation is. This is more surprising. In particular, it highlights that a large migration rate does not only reflect a strong sexual preference but implies more involved behavior. This will be corroborated in what follows.

Basins of attraction : I then explored how basins of attraction are modified when migration parameters increase. To simplify the study, I assumed here that m := m A = m a .

Figure 3 presents the trajectories of some solutions of dynamical system (3.1) in both phase planes which represent both patches. The trajectories are drawn for the initial condition z A,1 (0) = 4, z a,1 (0) = 10, z A,2 (0) = 8.5 and z a,2 (0) = 15, and for three different values of m: 0, 1 and 5. It is important to notice that the equilibrium reached depends not only on the initial condition but also on the value of m, unlike the symmetrical case. Indeed, on the example of Figure 3, when m is small, the trajectory converges to (0, ζ a , ζ A , 0). When m is larger, only a-individuals survive, the trajectory converges to (0, ζ a , 0, ζ a ). Hence, a large migration rate m can favor allele a, which codes for the weakest of both mating preferences (β a < β A ), to invade both patches.

Then, an example of basins of attraction D α,α m,m is given in the case of a large migration parameter (m = 5). Figure 4 presents the projections of the four sets on six different planes. More specifically, each graph (a-f) represents the equilibrium reached with respect to the initial condition in patch 1 for a couple of initial conditions in patch 2, which is plotted on graph (g). In order to compare results for m = 5 and m = 0, I plotted the line solution of (β A -1)z A,1 -(β a -1)z a,1 = 0 on all planes. Indeed, according to Lemma 3, without migration any trajectory with initial conditions in patch 1 above (resp. below) this line converges to a patch filled with a-individuals (resp. A-individuals). Generally, on the planes characterized by the values of (z A,2 (0), z a,2 (0)) given in captions. On each plane, the four sets from white to dark grey corresponds to initial conditions with convergence to (ζ A , 0, ζ A , 0), (ζ A , 0, 0, ζa), (0, ζa, ζ A , 0) and (0, ζa, 0, ζa) respectively. The black line is the solution of (β A -1)z A,1 -(βa -1)z a,1 = 0. (g): The black diamond points correspond to the initial conditions in patch 2 taken to obtain plots (a) to (f).

(a) z A,2 (0) = 4, z a,2 (0) = 5 (b) z A,2 (0) = 4, z a,2 (0) = 10 (c) z A,2 (0) = 4, z a,2 (0) = 15 (d) z A,2 (0) = 8, z a,2 (0) = 5 (e) z A,2 (0) = 8, z a,2 (0) = 10 (f) z A,2 (0) = 8,
observe that when the number of a-individuals is large in patch 1, these individuals are favored by a large migration rate. Thus, the conclusion here is that a large migration parameter m favors the allele coding for the weakest mating preference by mixing the populations of both patches.

Minimal number of individuals for invasion : In the following simulations, it is assumed that, initially, each patch is filled with a density of ζ a a-individuals and A-individuals are introduced in patch 1. To corroborate previous observations, I computed the minimal number of A-individuals that is needed to be introduced such that they can survive, i.e. such that the dynamical system converges to a stable equilibrium with a positive number of A-individuals. I computed this minimal number, denoted by N min (β A , m), for a range of values of β A (β A ∈ (1, 2]) and m (m ∈ [0, 2]), other parameters are defined by (3.11). On the left part of Figure 5, the number N min (β A , m) is drawn using a logarithmic color scale. Note that the minimal number of A-individuals required for survival decreases when β A increases. For example, when β A is large (β A = 2), observe that the minimal number of A-individuals needed for survival, is only half (resp. two-thirds) of ζ a = 20 when m = 0 (resp. m = 2). Moreover, if β A and m are sufficiently large (β A ≥ 2.9 and m ≥ 1.9 (data not shown)), the A-population replaces the resident a-population in both patches as soon as the initial number of A-individuals is equal to N min (β A , m). This suggests that individuals with a higher mating preference have a selective advantage.

Secondly, to better understand how m affects N min (β A , m), I computed the scaling difference A-individuals, especially if D(β A , m) is large. The opposite conclusion holds for negative value of D(β A , m). Here, when β A is smaller than β a = 1.5, D(β A , m) is positive and increases with migration m whereas, when it is smaller than β a , it is decreasing with m. Hence, migration seems here again to favor the allele with the weakest mating preference.

D(β A , m) := N min (β A , 0) -N min (β A , m) N min (β A ,

Proofs

This part is devoted to the proof of Theorem 1. The main idea is to start from the results without migration, then use a perturbation method to make m A and m a grow up and deduce results for some positive migration parameters. However, this perturbation technique will only apply on a bounded set of R E excluding 0. Thus, let us first prove Lemma 2, which allows us to restrict the study of the dynamical system (3.1) to the compact set S. To help with proofreading, we recall here the definitions of the weighted sums :

Σ i := (β A -1)z A,i + (β a -1)z a,i , for i = 1, 2, Σ := Σ 1 + Σ 2 = (β A -1)(z A,1 + z A,2 ) + (β a -1)(z a,1 + z a,2 ).
Proof (Proof of Lemma 2) The proof is based on the equations satisfied by Σ 1 , Σ 2 and Σ. From (3.1), we find

d dt Σ 1 = Σ 1 b Σ 1 z A,1 + z a,1 -2b(β A -1)(β a -1) z a,1 z A,1 (z A,1 + z a,1 )Σ 1 + b -d -c(z A,1 + z a,1 ) -m A (β A -1) + m a (β a -1) z A,1 z a,1 z A,1 + z a,1 - z A,2 z a,2 z A,2 + z a,2 . (4.1) Since Σ 2 1 -2(β A -1)(β a -1)z a,1 z A,1 ≥ 0 and Σ 1 ≥ (β min -1)(z a,1 + z A,1 ), we obtain d dt Σ 1 ≥ Σ 1 b -d - c (β min -1) Σ 1 -m A (β A -1) + m a (β a -1) z A,1 z a,1 (z A,1 + z a,1 )Σ 1 . (4.2)
We then find an upper bound of z A,1 za,1

(z A,1 +za,1)Σ1 : Σ 1 (z A,1 + z a,1 ) = (β A -1)z 2 A,1 + (β a -1)z 2 a,1 + (β A + β a -2)z A,1 z a,1 ≥ (β min -1)[z 2 A,1 + z 2 a,1 + 2z A,1 z a,1 ] ≥ 4(β min -1)z A,1 z a,1 .
In addition with (3.6) and (4.2), we deduce

d dt Σ 1 ≥ Σ 1 b -d 2 - c (β min -1) Σ 1 .
Hence, as soon as

Σ 1 < (β min -1)(b -d)/2c, its derivative is strictly positive. In other words, if Σ 1 (0) ≤ (β min -1)(b -d)/4c
, there exists t 1 > 0 such that for all t ≥ t 1 , Σ 1 (t) is higher than this threshold. Moreover, if Σ 1 (t 2 ) is higher than this threshold, for all t ≥ t 2 , Σ 1 (t) remains higher than it. The same conclusion holds for Σ 2 .

Let us now deal with Σ. From equations (4.1) satisfied by Σ 1 and Σ 2 , we find

d dt Σ = i=1,2 Σ i b β A z A,i + β a z a,i z A,i + z a,i -d -c(z A,i + z a,i ) -2b(β A -1)(β a -1) z A,i z a,i z A,i + z a,i ≤ i=1,2 Σ i bβ max -d - c β max -1 Σ i ≤ Σ bβ max -d - c 2(β max -1) Σ .
Using a reasoning similar to the previous one, we conclude that there exists a time after which Σ(t) remains lower than 4(β max -1)(bβ max -d)/c. Finally any trajectory hits S after a finite time and S is a positively invariant set. That ends the proof of Lemma 2.

Lemma 2 implies that the study of the dynamical system (3.1) can be restricted to the study of trajectories belonging to S. Note that when m A = m a = 0, Subsection 3.2 ensures that the dynamical system (3.1) has exactly 9 equilibria which belong to S:

(ζ A , 0, 0, ζ a ), (ζ A , 0, ζ A , 0), (0, ζ a , ζ A , 0), (0, ζ a , 0, ζ a ), (4.3) (χ A , χ a , ζ A , 0), (χ A , χ a , 0, ζ a ), (0, ζ a , χ A , χ a ), (ζ A , 0, χ A , χ a ). (4.4) (χ A , χ a , χ A , χ a ). (4.5)
Equilibria (4.3) are stable fixed point whereas equilibria (4.4) (resp. (4.5)) are unstable with a local stable manifold of dimension 3 (resp. 2), i.e. there exists a set of dimension 3 (resp. 2) such that any trajectory starting from this set converges to equilibria (4.4) or (4.5).

In order to simplify the notation of the proofs, let us write migration rates m A and m a using three parameters γ A ∈ [0, 1], γ a ∈ [0, 1] and m ≥ 0 as m A := mγ A and m a := mγ a .

We consider that γ A and γ a are fixed parameters and we will make m grow up in the following proof. We can rewrite the dynamical system (3.1) considering m as a parameter

d dt z(t) = F (z(t), m). (4.6)
The solution of (4.6) with initial condition z 0 is written t → ϕ m,z 0 (t). Our goal is to understand the dynamics of the flow ϕ m,z 0 associated to the vector field F (z, m) using ϕ 0,z 0 (without migration) which is entirely described in Subsection A. Theorem 1 can be rewritten as follows using the notion of flow.

Theorem 3 (Theorem 1') There exists m 0 > 0 such that for all m ≤ m 0 , we can find four open subsets (D α,α m ) α,α ∈A of S with the following properties:

-The closure of ∪ α,α ∈A D α,α m is equal to S.
-For all z 0 ∈ D A,a m , the flow ϕ m,z 0 (t) converges to (ζ A , 0, 0, ζ a ) when t tends to +∞. Similar results hold for the three other equilibria (4.3).

Proof The first step is to construct a neighborhood around each equilibrium (4.3)-(4.5) which includes a unique equilibrium of the dynamical system (3.1) with m > 0. Let us first focus our study on the equilibrium (ζ A , 0, 0, ζ a ). Subsection (A) implies that, when m = 0, the equilibrium (ζ A , 0, 0, ζ a ) is an attractive stable equilibrium. The first derivative

D z F evaluated at (z, m) = ((ζ A , 0, 0, ζ a ), 0) is     -(bβ A -d) -b(β A -1) -(bβ A -d) 0 0 0 -b(β A -1) 0 0 0 0 -b(β a -1) 0 0 0 -b(β a -1) -(bβ a -d) -(bβ a -d)     .
(4.7)

Since matrix (4.7) is invertible and F is smooth on S × R + , the Implicit Function Theorem implies that there exists m 1 and a neighborhood V 1 of (ζ A , 0, 0, ζ a ) in S such that there is a unique point y 1 (m) ∈ V 1 satisfying F (y 1 (m), m) = 0 for all m < m 1 . And m → y 1 (m) is regular and converges to (ζ A , 0, 0, ζ a ) when m converges to 0. A simple computation ensures that F (y 1 (0), m) = F ((ζ A , 0, 0, ζ a ), m) = 0, for any m > 0. Since y 1 (m) is unique, we deduce that y 1 (m) = y 1 (0). Moreover, from Theorem 6.1 and Section 6.3 of [START_REF] Ruelle | Elements of differentiable dynamics and bifurcation theory[END_REF] (see also Appendice B of [START_REF] Collet | A rigorous model study of the adaptive dynamics of mendelian diploids[END_REF][START_REF] Hoppensteadt | Singular perturbations on the infinite interval[END_REF]), we conclude that if m 1 and V 1 are small enough, any solution ϕ m,z 0 with z 0 ∈ V 1 and m < m 1 converges uniformly to ϕ 0,z 0 when m converges to 0. In other words, y 1 (0) attracts all the orbits ϕ m,z 0 starting from V 1 . Similarly, we find (m i ) i=2,3,4 and (V i ) i=2,3,4 neighborhoods around the three other equilibria of (4.3), denoted by (y i (0)) i=2,3,4 , such that, for i ∈ {2, 3, 4}, for all m < m i , y i (0) attracts all solutions ϕ m,z 0 with z 0 ∈ V i and m < m i . Theorem 6.1 and Section 6.3 of [START_REF] Ruelle | Elements of differentiable dynamics and bifurcation theory[END_REF] ensure also the stability of the local stable and unstable manifolds of a hyperbolic non-attractive fixed points. Thus, we find m 5 , .., m 9 and V 5 , .., V 9 , neighborhoods around equilibria (4.4) and (4.5) that satisfy the following properties. For all i ∈ {5, .., 9}, for all m < m i , there exists a unique fixed point y i (m) ∈ V i invariant by F (., m) which repulses all orbits solution associated to F (., m), except the orbits that start from a surface of dimension 3 or 2, depending on whether we are focused on equilibria (4.4) or (4.5) respectively. These surfaces are the stable manifolds of (y i (m)) i=5,..,9 in (V i ) i=5,..,9 respectively. Without loss of generality, we assume that these nine neighborhoods are disjoint sets.

The second step is to deal with trajectories outside these nine neighborhoods. Let ε > 0 and for i = 1, .., 9, we define

V ε i = B(y i (0), R i ) = {z ∈ S, z -y i (0) ≤ R i }, where R i = max{r > 0, B(y i (0), r + ε) ⊂ V i },
which is a neighborhood of y i (0) slightly smaller than V i .

Recall that the five neighborhoods (V ε i ) i=5,..,9 attracts some solutions ϕ 0,z 0 . Thus, we set

W =   i=5,..,9 z 0 ∈V ε i (ϕ 0,z 0 ) -1 ([0, +∞))   S, (4.8)
which is a neighborhood of the union of all stable manifolds of unstable equilibria (4.5) and (4.4) assuming m = 0. We denote the complement of W in S by W c . Let us first deal with the trajectories starting from W c . According to Appendix A, all trajectories ϕ 0,z 0 starting from W c converge to a stable equilibrium, i.e. they reach any neighborhood of set {y i (0), i = 1, ., 4} in finite time. Since W c is compact, there exists a finite time t 1 > 0 such that ϕ 0,z 0 (t 1 ) ∈ ∪ 4 i=1 V ε i , for all z 0 ∈ W c . Moreover, from Theorem 1.4.7 by [START_REF] Berger | Géométrie différentielle: variétés, courbes et surfaces. Presses Universitaires de France Blaustein AR[END_REF], the flow ϕ is uniformly continuous with respect to m, to z 0 and to t. We can thus find m 10 < min i=1,..,9 m i such that for every m ≤ m 10 , z 0 ∈ W c ϕ 0,z 0 (t 1 ) -ϕ m,z 0 (t 1 ) ≤ ε.

Then, by definition of (V i ) i=1,..,4 and (V ε i ) i=1,..,4 , we deduce that for all m ≤ m 10 , all z 0 ∈ W c and all t ≥ t 1 ,

ϕ m,z 0 (t) ∈ 4 i=1 V i .
Then we deal with the trajectories starting from W. According to the definition of W (4.8), all trajectories ϕ 0,z 0 starting from W reach one of the five neighborhoods (V ε i ) i=5,..,9 in finite time. Thus, by reasoning as above, we can find m 11 ≤ m 10 and t 2 > 0 such that for all m ≤ m 11 and all z 0 ∈ W, there exists t ≤ t 2 , with

ϕ m,z 0 (t) ∈ 9 i=5 V i .
Let us fix m ≤ m 11 , z 0 ∈ W and assume that ϕ m,z 0 (t 3 ) ∈ V i . We have then three possibilities:

(i) If ϕ m,z 0 (t) ∈ V i for all t ≥ t 3 , then z 0 belongs to the stable manifold of y i (m) in S. Since we have a global diffeomorphism on S, we can find the stable manifold of y i (m) by iterating the Implicit Function Theorem and deduce that this stable manifold is a set of dimension 3 or 2, depending on which equilibrium is considered. (ii) Otherwise, there exists t 4 ≥ t 3 such that ϕ m,z 0 (t 4 ) ∈ V i . If ϕ m,z 0 (t 4 ) ∈ W c , the flow will converge to one of the four equilibria (4.3) according to previous reasoning.

(iii) The last possibility is ϕ m,z 0 (t 4 ) ∈ W \∪ 9 i=5 V i . Thus, the flow (ϕ m,z 0 (t)) t≥t4 will reach again one of the neighborhoods (V j ) j=5,..,9 . It would have a problem if the trajectory went from a neighborhood to an other without living W as t → +∞. Thus, let us show that this is not possible. Indeed, the flow goes out of V i by following the unstable manifold of y i (m) which is close to the unstable manifold of y i (0) (according to the continuity of the unstable manifolds with respect to m, cf Theorem 6.1 by [START_REF] Ruelle | Elements of differentiable dynamics and bifurcation theory[END_REF]). Since ϕ m,z 0 leaves V i by staying in W, the intersection of the unstable manifold of y i (0) and W is not empty. From the definition of W (4.8) and Appendix A, it is possible if and only if y i (0) = y 9 (0) = (χ A , χ a , χ A , χ a ) and if ϕ m,z 0 leaves V 9 through the neighborhood of the stable manifold of one of the equilibria (4.4). Thus, the flow (ϕ m,z 0 (t)) t≥t4 will reach one of the neighborhood (V j ) j∈{5,6,7,8} . Then, only the two previous possibilities (i) or (ii) are possible.

Finally, we have shown that any solution ϕ m,z 0 of (4.6) starting from S and with m ≤ m 11 converges to one of the equilibria (4.3), except if it starts from a set with empty interior which is the union of the global stable manifolds of the equilibria (y i (m)) i=5,..,9 .

Finally, m 0 := m 11 , D A,a m = ∪ z 0 ∈V1 ϕ m,z 0 -1 ([0, +∞)), and D A,A m , D a,A m , D a,a m are defined in a similar way using sets V 2 , V 3 and V 4 respectively. We have shown that for all m ≤ m 0 , the four non empty interior sets (D α,α m ) α,α =A,a satisfy Theorem 3.

Discussion

The first conclusion of the study is that a population with a large mating preference has selective advantages: (1) the larger the mating preference strength is, the shorter is the time before reaching an equilibrium where this allele is maintained, and (2) a population with a strong mating preference can invade a resident population with a weak preference even if its initial number of individuals is small. Same kind of conclusion is drawn by [START_REF] Smadja | Strong premating divergence in a unimodal hybrid zone between two subspecies of the house mouse[END_REF]. In the latter, the authors predict that the asymmetrical mating preference observed between two species of mouse could lead to the replacement of the subspecies with the weakest mating preference (M. m. domesticus) by the other subspecies (M. m. musculus), if no other mechanism was involved. This conclusion is a substantial added value compared to [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] where only the symmetric case (β A = β a , m A = m a ) is considered.

Accounting for asymmetrical preference gave the possibility to better understand advantages of a strong mating preference. Migration has a more involved impact on the system dynamics than mating preference, although the frequency-dependent term I used for migration seemed only to mimic mating preferences. More precisely, there exists a trade-off between two phenomena [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]: (1) large migration rates can help individuals to escape disadvantageous patches [START_REF] Clobert | Dispersal[END_REF] but (2) large migration rates entail also risks of moving to unfamiliar patches (i.e. filled with not-preferred individuals) and thus may increase the time before reproductive isolation. This is understandable since the migration terms only focus on the departure patch. More surprisingly, large migration rates seem to favor alleles with reduced mating preferences. This tendency was not noticed by [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] and could be linked to the effects of migration on habitat specialization [START_REF] Brown | Evolution in heterogeneous environments: effects of migration on habitat specialization[END_REF][START_REF] Cuevas | Evolution of rna virus in spatially structured heterogeneous environments[END_REF][START_REF] Elena | The evolution of viruses in multi-host fitness landscapes[END_REF]. In these articles, the authors highlight that migration may prevent the local specialization of subpopulations and favor generalist species. Hence, in both cases, large migration rates tend to avoid specialized behaviors in terms of ecological niche adaptation or mating partner adaptation.

any α-individual gives birth to a α-individual at rate b, any α-individual gives birth to a ᾱ-individual at rate mᾱ, any α-individual dies at rate bβᾱ + mα.

The goal is thus to estimate the extinction time of such a sub-critical two types branching process. Let M (t) be the mean matrix of the multitype process, that is, where ca, c A are two positive constants and r is the largest eigenvalue of the matrix G. With a simple computation, we find that r = -ω(A, a). From (B.2), we deduce that the extinction time is of order ω(A, a) -1 log K when K tends to +∞ by arguing as in step 2 of Proposition 4.1's proof of [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. This concludes the proof of Lemma 4.

M (t) =     E E
Finally, this gives all elements to induce Theorem 2.

Fig. 1

 1 Fig. 1 Example of phase portrait of a single patch dynamics. The black line is the set of initial conditions for which trajectories converge to the unstable equilibrium (χ A , χa). Above (resp. below) this line, trajectories converge to (0, ζa) (resp. (ζ A , 0)). The red and the blue curves are examples of trajectories. The black and the gray points represent respectively the stable and the unstable equilibria.

  .5, b = 2, d = 1 and c = 0.1, (3.11) unless stated otherwise. For these parameters, ζ A = 30 and ζ a = 20.

Fig. 2

 2 Fig.2Graphs of the constants in front of the times before reproductive isolation, ω(A, a) -1 (blue line), ω(A, A) -1 (red dashed line), ω(a, a) -1 (red dashed-dotted line), with respect to β A (left) and to m (right). The demographic parameters are defined by (3.11), γ A = 1, γa = βa -1 = 0.5 and m = 2 on the left and β A = 2 on the right.

Fig. 3

 3 Fig.3Plots of the trajectories in the phase planes which represent the patch 1 (left) and the patch 2 (right) for t ∈ [0, 10] and for three values of m: m = 0 (red), m = 1 (blue), m = 5 (green). The initial condition is(4, 10, 8.5, 15), represented by the black dots. The dark line is the solution of (β A -1)z A -(βa -1)za = 0

  Fig. 4 (a-f): Projections of sets D α,α 5

Fig. 5

 5 Fig.5Left: Minimal number of initial A-individuals in patch 1 that is needed for a long time survival when starting from two patches filled with ζa a-individuals; a logarithmic color scale is used. Right: Scaling differences between the minimal number of A-individuals needed without migration (i.e. (βa -1)ζa/(β A -1)) and the one computed on the left plot. Parameters are defined by (3.11) where βa = 1.5.

  Na(t) (Na(0), N A (0)) = (1, 0) E E N A (t) (Na(0), N A (0)) = (1, 0) E E Na(t) (Na(0), N A (0)) = (0, 1)E E N A (t) (Na(0), N A (0let G be the infinitesimal generator of the semigroup {M (t), t ≥ 0}. From the book of Athreya and Ney (1972) p.202, we deduce a formula of G which isG = -b(β A -1) -ma m A ma -b(βa -1) -m A .Applying Theorem 3.1 of Heinzmann et al (2009), we find thatP (Na(t), N A (t)) = (0, 0) (Na(0), N A (0)) = (z 0 a,1 K, z 0 A,2 K) = (1 -cae rt ) z 0 a,1 K (1 -c A e rt ) z 0 A,2 K , (B.2)

Acknowledgements I would like to thank Pierre Collet for his help on the theory of dynamical systems. I also thank two anonymous reviewers for their constructive comments, which helped me to improve the manuscript. This work was partially funded by the Chair "Modélisation Mathématique et Biodiversité" of VEOLIA-Ecole Polytechnique-MNHN-F.X.

A Dynamical system without migration

In this appendix, we will prove the results of Section 3.2, which is related to the case without migration. To this aim, we use the two following weighted quantities Ω(t) := (β A -1)z A (t) -(βa -1)za(t), Σ(t) := (β A -1)z A (t) + (βa -1)za(t).

From (3.9), we find that

Proof (Proof of Lemma 3) We start by studying the stability of equilibrium (0, 0). Assume that Σ(0) > 0. From (A.2), we derive

.

Hence, as long as Σ ∈]0, (b -d)(β min -1)/c[, Σ(t) is increasing. Thus (0, 0) is an unstable equilibrium.

The stability of the three other equilibria, (ζ A , 0), (0, ζa) and (χ A , χa), can be deduce by a direct computation of Jacobian matrices at these points, which we do not detail.

Finally, let us study the long time behavior of any solution. Equation (A.1) implies that the sign of Ω(t) is equal at all time and, that D A 0 is a positively invariant set under dynamical system (3.9). Moreover, there exists only a stable equilibrium that belongs to the set D A 0 , which is (ζ A , 0). We consider the function W :

From (A.1) and (A.2), we deduce that

Moreover for any (z A , za) ∈ D A 0 , W (z A , za) = 0 if and only if za = 0. W (z A , za) converges to +∞ when (β A -1)z A -(βa -1)za converges to 0 and dW dt is non-positive on D A 0 and is equal to zero if and only if za = 0. It ensures that W is a Lyapunov function for (3.9) on the set D A 0 which cancels only on D A 0 ∩ {za = 0}. Furthermore, a simple computation gives that the largest invariant set in 1960) is thus sufficient to conclude that any solution of (3.9) with initial condition in D A 0 converges to (ζ A , 0) when t tends to +∞. Similarly, we prove that any solution with initial condition in D a 0 converges to (0, ζa). Finally, assume that Ω(0) = 0. Then, Ω(t) = 0 for all t ≥ 0 according to (A.1) and, in addition with (3.9), we derive for all α ∈ A,

We deduce the last point of Lemma 3 easily.

B Extinction time

This subsection is devoted to the proof of Theorem 2 following ideas similar to the ones of the proof of Theorem 3 and Proposition 4.1 in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. Hence, we do not give all details, but explain only parts that are different. Assume that m A ≤ m 0 , ma ≤ m 0 and that Z K (0) converges in probability to a deterministic vector z 0 belonging to D A,a m A ,ma , Lemma 3 and Theorem 1 ensure that (Z K (t), t ≥ 0) reaches a neighborhood of the equilibrium (ζ A , 0, 0, ζa) after a finite time independent from K. Indeed, the process dynamics is close to the one of the limiting deterministic system (3.1).

To prove Theorem 2, it remains to estimate the time before all a-individuals in patch 1 and all A-individuals in patch 2 disappear. We denote it by

and we assume that the process is initially close to equilibrium (ζ A , 0, 0, ζa). The estimation is deduced from the following Lemma.

Lemma 4 There exist two positive constants ε 0 and C 0 such that for any ε ≤ ε 0 , if there exists

Proof Following the first step of Proposition 4.1's proof given by [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF], we prove that as long as the population processes Z K a,1 (t) and Z K A,2 (t) have small values, the processes Z K A,1 (t) and Z K a,2 (t) stay close to ζ A and ζa respectively.

Then, by bounding death rates, birth rates and migration rates of (Z K a,1 (t), t ≥ 0) and (Z K A,2 (t), t ≥ 0), we are able to compare the dynamics of these two processes with the ones of

where (Na(t), N A (t)) ∈ N {a,A} is a two-types branching process with types a and A and for which