Pierre-Emmanuel Leni

Rémy Laurent

Michel Salomon

Régine Gschwind

Libor Makovicka

Julien Henriet

Development of a 4D numerical chest phantom with customizable breathing

Keywords: Artificial Neural Network, Radiation Physics, Phantoms, Breath Simulation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The ICRP (International Commission on Radiological Protection) phantoms are continuously evolving and being refined through research. Indeed, various mathematical phantoms, such as MIRD [START_REF] Snyder | Medical internal radiation dose committee (MIRD) pamphlet no. 5[END_REF], have been used since 1969. Nevertheless, they were hermaphrodite and nothing existed for children. Therefore, new phantoms appeared, [START_REF] Kramer | The calculation of dose from external photon exposures using reference human phantoms and Monte Carlo methods, Part I: The male (Adam) and female (Eva) adult mathematical phantoms[END_REF][START_REF] Cristy | Specific Absorbed Fractions of Energy at Various Ages from Internal Photon Sources: 1, Methods[END_REF], refining the organ representation. Nowadays, the ICRP uses voxelized phantom references based on computed tomography images [START_REF] Menzel | ICRP Publication 110. Adult reference computational phantoms[END_REF], to resemble as closely as possible, the shapes, volumes and weights of various organs.

Email address: pierre_emmanuel.leni@univ-fcomte.fr (Pierre-Emmanuel LENI)

In parallel to ongoing research regarding phantoms, respiratory movement is also a growing research topic.

In 2005, Low et al. estimated respiratory movement using a 5-dimension modelization representing object coordinates, current volume and airflow [START_REF] Low | Novel breathing motion model for radiotherapy[END_REF][START_REF] Yang | 4D-CT motion estimation using deformable image registration and 5D respiratory motion modeling[END_REF][START_REF] Zhao | Characterization of free breathing patterns with 5D lung motion model[END_REF]. Other models also exist: in 2010, Eom et al. [START_REF] Eom | Predictive modeling of lung motion over the entire respiratory cycle using measured pressure-volume data, 4DCT images, and finite-element analysis[END_REF] developed a simulation based on the pressure-volume relationship linked to a finite element analysis; Vandemeulebrouke et al. [START_REF] Vandemeulebroucke | The POPI-model, a point-validated pixel-based breathing thorax model[END_REF][START_REF] Vandemeulebroucke | Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[END_REF] introduced the PoPi model, which relies heavily and a priori about the movement, since it is built using a 4DCT of the patient. The motion estimation is computed by comparing the images obtained at each respiratory phase with a simulated sequence (built using the spatial registration of images correlated with a cyclic model of breathing). It is also possible to simulate, almost in real time, respiratory movement customized to every patient, using neural networks [START_REF] Laurent | Respiratory lung motion using an artificial neural network[END_REF]; and to use external surrogate monitoring to drive a tumor movement model ( [START_REF] Martin | Building motion models of lung tumours from cone-beam CT for radiotherapy applications[END_REF]). Consequently, phantoms tend to take motion into account to obtain a simulation both accurate and as close as possible to physiologic reality. Segars integrated a module for cardiac deformation into his NCAT model [START_REF] Segars | Modeling respiratory mechanics in the MCAT and spline-based MCAT phantoms[END_REF][START_REF] Segars | Development and application of the new dynamic Nurbs-based Cardiac-Torso (NCAT) phantom[END_REF]. Furthermore, respiratory simulation was implemented, but remains quite limited: it mainly consists in lung deformations and translation/rotation of other organs, such as the liver and kidneys. Respiratory simulation also exists in the XCAT phantom, which relies mainly on lung deformation. This deformation is parametrized using two curves, which offers high flexibility but requires an expert to adapt it for a given patient [START_REF] Segars | 4D XCAT phantom for multimodality imaging research[END_REF].

Respiratory movement is useful for radiation treatment: a customized breathing phantom associated with patient monitoring during a treatment session allows for detailed monitoring of received doses and could lead to treatment adjustment. Nevertheless, from a radiation protection point of view, 4DCT provides movement information at the cost of increasing patient exposure to radiation (around 30 times the annual natural exposure). Moreover, due to technical limitations, the data may be inaccurate and contain artifacts (classification problem for example). In this context, using a customized simulation for breathing offers advantages for both radiation protection, by limiting the use of 4DCT that could be a significant step in avoiding small doses, whom effects are still not well documented; and for radiation treatment, by providing movement information for more accurate treatment planning as well as monitoring the dose actually delivered during a treatment session via the breathing cycle recorded from the patient.

Using Artificial Neural Networks (ANN), we proposed two methods: for the simulation of pulmonary movement [START_REF] Laurent | Respiratory lung motion using an artificial neural network[END_REF], and for the adaptation of the IRSN (French radiation protection and nuclear safety institute) phantoms [START_REF] Farah | Development of a new CBR-based platform for human contamination emergency situations[END_REF][START_REF] Henriet | EquiVox: an example of adaptation using an artificial neural network on a case-based reasoning platform[END_REF]. These two projects aim at developing adapted tools for a customized and accurate modelization of the human body, either in a radiation protection context (anthroporadiametry), or for delivering external beam treatment more accurately. In radiation therapy, ANN has recently been used to simulate the tumor movement according to external position trackers [START_REF] Seregni | Real-time tumor tracking with an artificial neural networks-based method: a feasibility study[END_REF].

In this paper, our contribution includes a new customized phantom model that simulates realistic respiratory motion of the lungs. The breathing cycle is simulated and validated using 4DCT data. The Dice indices are computed for every respiratory phase section 3. We also present some preliminary results for simulations of cardiac and oesophageal motion (section 4).

Method: ANN presentation

The neural network

The NEMOSIS (NEural NEtwork MOtion SImulation System) platform has been presented and validated in previous work [START_REF] Laurent | Respiratory lung motion using an artificial neural network[END_REF]: using an artificial neural network, it is possible to simulate pulmonary movement customized per patient. Nevertheless, customization was limited by a small amount of patient data, resulting in small learning sets. Moreover, only the internal anatomic structure of the lungs was considered. In this platform, a multi-layer perceptron with one hidden layer is trained using supervised learning. The synaptic weights and biases are optimized to minimize the mean squared error (MSE) using the gradient based Limited-memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) [START_REF] Hsieh | Machine learning methods in the environmental sciences: neural networks and kernels[END_REF] combined with the Wolfe linear search [START_REF] Nocedal | Numerical optimization[END_REF] to determine the optimal step size. Moreover, 1000 iterations are performed beyond the minimum MSE to prevent overfitting by checking that a degradation is observed on the interpolation of the generalization set. An incremental approach determines the optimal number of neurons in the hidden layer [START_REF] Sauget | Efficient domain decomposition for a neural network learning algorithm, used for the dose evaluation in external radiotherapy[END_REF]. The application of NEMOSIS to anthropomorphic phantoms requires several adjustments:

• Adding more patient data to the learning set via collaboration with Besançon CHRU,

• Adding the overall lung contours to the data (for 8 patients over 16),

• Modification of the ANN to improve the customization and compatibility with the IRSN phantoms : find relevant characteristics of the lung to be used as ANN entries.

The IRSN phantoms were generated by proportionally adapting an initial model for various heights [START_REF] Clairand | Improvement of internal dose calculations using mathematical models of different adult heights[END_REF][START_REF] Farah | Examples of Mech and NURBS phantoms to study the morphology effect over in vivo lung counting[END_REF]. Nevertheless, the analysis and the comparison between patient data and the phantoms demonstrated differences between real and theoretical morphologies. To represent the size of the lungs in the Superior-Inferior axis (SI axis), we take into account the size of the right lung (which is independant of the heart).

The measures on the Anterior-Posterior axis and Left-Right axis (AP and LR axes) are performed along the axial plane of the carina to allow reproducibility in all patients and phantoms. Using these two axes, we have approximated the lung bounding box by using a cylinder generated by an ellipse. The ellipse perimeter

P (P = Π 2(( dAP 2 ) 2 + ( dLR 2 ) 2 )
) is presented in the last column of the tables 1 and 2. This bounding box is used as an ANN entry.

Tables 1 and2 are sorted according to pulmonary volume. Table 1 shows that there is no correlation between patient height and lung volume. For example, a person measuring around 1.65 meters can have a pulmonary volume of 2.53 liters (163 cm) or 4.15 liters (165 cm). Similarly, the size of the lungs along the vertical axis (i.e. SI axis) is not a linear function of height. Nevertheless, pulmonary volume can be represented by the volume of the aforementioned cylinder. The ANN concept is extensively described in the literature [START_REF] Laurent | Simulation du mouvement pulmonaire personnalisé par réseau de neurones artificiels pour la radiothérapie externe[END_REF]. The entries of the ANN are raw data to avoid biasing the learning. Therefore, the existing entries (3D coordinates of a point (its location along the 3 axis X, Y and Z), pulmonary volume, and respiratory phase) have been complemented by three new ones, which represent an approximation of the real pulmonary volume by a cylinder, i.e. the dimensions along the SI, AP and LR axes. Unfortunately, we note that the dimensions along the AP axis for all the phantoms are inferior to the patient's. As already stated [START_REF] Laurent | Respiratory lung motion using an artificial neural network[END_REF], NEMOSIS is an interpolator inside the definition domain of the training data. Since the learning step is performed using patient data, the dimensions along the AP axis cannot be used for simulation on the phantoms. Therefore, the entries using the dimensions along AP and LR axes are replaced by the perimeter.

Moreover, the perimeter is correlated with chest size and thus can be "externally" measured. Only phantoms with a height greater than 173.18 cm are included in the definition domain; thus, further studies are limited to these phantoms. Ultimately, 7 variables constituted the entries of the NEMOSIS ANN:

• 3D point coordinates,

• Pulmonary volume,

• The cylinder perimeter,

• The dimension of the right lung (i.e. the cylinder height), and

• The respiratory phase.

In other words, given the external measurement of a person and the desired respiratory phase, the network can simulate the corresponding customized lung contours.

The training data

The learning set for pulmonary movement was constituted by 4DCT images of 15 patients acquired during free breathing and sorted into 5 phases by the 4DCT. Firstly, on every 4DCT, positions of anatomic points are manually defined by experts and tracked over various respiratory phases, for a total of 4680 points (10 phases per patient). These points are internal characteristic locations within the lungs that are recognizable over the phases (e.g. vessels, bronchus unusual junctions). Since the IRSN phantoms consist of the organ contours only, an algorithm (region-growing) has been implemented to automatically extract the lung contours of the patients on every 4DCT. Indeed, there is no visual reference to identify a point over phases, so we take advantage of parametric curves:

• Every lung has 4 ridge lines (top and bottom most points along the LR and AP axis),

• A ridge line is a spline (a parametric curve). The control points of the splines are the extreme points of the lungs on every section,

• An identical parameter value over time defines the same point over different phases.

90 to 115 points per ridge line are computed, which leads to 720 to 920 contour points per patient. A total of 3248 points per phase are added to the dataset, which leads to 37160 points for 15 patients over 10 phases, including the anatomic points.

The Figure 1 illustrates the spline interpolation (in red) of the lung contours. For any number of points guiding the spline, a position along the curve can be defined by a parameter (a real value) varying between 0 and 1.In this example, three points (in yellow) are drawn, corresponding to the parameterization values 0.5, 0.7 and 0.8 (chosen for illustration purpose). We can observe that these locations can be represented at every phase, by using the same parameter value, even if the spline is not defined by the same number of points. The automatic incrementation algorithm of the number of neurons presented in [START_REF] Laurent | Respiratory lung motion using an artificial neural network[END_REF] is used. The dataset is split into two sets: a learning set and a generalization set. The generalization set is used to check the ANN interpolation performances during the learning step,to prevent over-fitting, and to stop the process when it starts diverging from the best solution. The results presented in the following section correspond to the best results : the learning set contains 12 patients, and the generalization set contains 3 patients. 20 neurons between the entries and the output layer have been determined as the best topology by the learning algorithm. Various constructions of the learning and generalization sets have been tested for the ideal data division.

Results: Lung movement and validation

Simulated movement

The number of phases to simulate can be adjusted. Two configurations are considered: 10 phases and 100 phases. To emphasize the movement gradients, we present the results obtained for 10 phases. 100 phases may allow the simulation of any lung movement (such as a sigh) or irregular breathing. Table 3 presents the simulation of respiratory movement using NEMOSIS on IRSN phantoms. The two phantoms represented are at the boundary of the definition domain (cylinder perimeter). The maximum exhalation point, i.e. the minimum pulmonary volume, is reached between the 50 to 60% respiratory phases. The variation between the various 4D TDM are due to the acquisition uncertainties and synchronization with the breath signal.

NEMOSIS computes the minimum at the phase 60% when 10 phases are simulated (as in Table 3), and 53% when 100 phases are generated.

HeightPhantom volume AN N 0% volume ∆V ol P hant-AN N AN N 60% volume ∆V ol 60%-0% respiratory phases: 0% and 60%, the minimum and maximum lung volumes, respectively. The ∆ measures represent the difference between the phantom and the RNA, and between the phases 60% and 0%.

(cm) (l) (l) (l) (l) (l
Figure 2 presents the breath simulation obtained for the 185 cm phantom using NEMOSIS. The phantom is rendered at 5 respiratory phases (10, 30, 50, 70 and 90%). Since the number of points within the contour is constant, a color code is used to demonstrate the amplitude of movement of a point from its current position to its position at the next phase. The initial lung contour (phase 0%) is represented on every phase to visualize the global movement amplitude. In this case, the maximum magnitude is 12.54 mm. The mean and standard deviations are computed for all the point variations at every phase. The main direction of the movement gradient is along the SI axis, which corresponds to reality. Indeed, despite the numerous muscles involved in respiration (intercostals, scalenes, diaphragm), the main one is the diaphragm. Its contraction deforms the lungs, including their contour as well as their internal anatomy.

Validation

Validation is conducted by comparing ANN simulations with the 4DCT images of 3 patients, phase by phase. The patients had lung dimensions and contour coordinates included in the learning set minimum and maximum values. These data were not used during the training step.

Table 4 shows the pulmonary volumes computed using the simulated contours on the test patients. In this case, since the number of points is not the same for every phase, it is not possible to compute the error per point. Instead, the volumes are compared using the Dice index, which evaluates the similarity between two sets (see equation 1, V 4D and V N EM representing the measured volume on 4DCT and simulated using NEMOSIS, respectively). Unfortunately, some 4DCT were not available for some phases so the corresponding comparisons were not computed. The Figure 3 illustrates lung contours for patient A extracted from 4DCT data at the 70% respiratory phase, and the corresponding contours simulated by NEMOSIS. We note that the contours are not distorted.

D(V 4D , V N EM ) = (2card(V 4D ∩ V N EM )) (card(V 4D ) + card(V N EM )) (1)
On this example, the Dice index is equal to 0.97.

4. Discussion: limits of current approach and extension to other organs

Simulated lung movement

We propose an ANN-based approach to simulate the organ contour at any respiratory phase. Our solution has been tested on both phantoms, and it can be noted that NEMOSIS simulates the largest variations during the breathing cycle for the tallest phantom. This behavior is related to the phantoms: they are generated from a model by applying a proportionality coefficient along the 3 dimensions, without any morphologic difference. To validate the approach, we compare the ANN results with 4DCT data using Dice index. Table 4 shows that for patient A, the maximum exhalation point is obtained at 53% using NEMOSIS, but measured on the real patient at 60%: while the real patient is still inhaling, the simulation has already started exhaling. Thus, for the phases 60% to 90%, NEMOSIS is "in advance" and provides higher volumes than reality. Despite this, the correlation between NEMOSIS and reality is very promising, providing a Dice index higher than 0.94. By construction, the ANN can only perform the interpolation for the input data included inside the learning set. For example, if the patient heights for training are between 160cm and 165cm, to be able to use ANN for a 170cm person, the learning sets has to be extended.

Cardiac and OEsophageal movement

Only three 4DCT are currently available with heart and oesophageal information. Moreover, these data are incomplete, with missing phases and/or incoherent data. Indeed, these organ contours are manually defined. Nevertheless, since the data are acquired using breath synchronization (and thus 10 phases), the main movement measured on the heart and oesophagus is likely due to the breathing and lung deformation.

Therefore, for a preliminary study, we simulate the movement of both organs using the ANN used for lungs.

Although it is not an ANN dedicated to these organs, the training using the lungs allows for presenting interesting results, but requires further improvements. It certainly illustrates the feasibility and flexibility of our approach. To compare the movement simulated by NEMOSIS with the 4DCT data, only the sections (i.e. CT images along the SI axis) having mutual contour definitions are retained. Table 5 presents the oesophagus and the heart dice index, between the simulated volumes using NEMOSIS and the 4DCT data.

that the contour definitions (on 4DCT) are consistent over all phases. The Dice index of patient 3 is lower than 0.75, due to poor contour definition on the 4DCT.

D(V 4D , V N

Figure 1 :

 1 Figure 1: Spline interpolation of the lung contours.

Figure 2 :

 2 Figure 2: Simulation of the lung movement using NEMOSIS on the 185 cm phantom.

Figure 3 :

 3 Figure 3: Comparison between 4DCT data (in blue) and NEMOSIS (in red) at the phase 70%.

Figure 4

 4 Figure4illustrates these results for the patient 1 in the worst case. All the sections available either through ANN or 4D are represented. We can notice that the shape is preserved, and the anatomic characteristics (such as the curvature at the bottom) are properly located.

Figure 4 :

 4 Figure 4: Comparison of patient 1 oesophageal shape at the 50% phase. In black: expert's contour definition, in blue: NEMOSIS contours.

Table 1 :

 1 Measures (in mm) of the lungs for various patients, obtained from computed tomography images, besides the perimeters, computed using the measures along the AP and LR axes.

	Height Age/Gender	Pulmonary volume Right lung (liter) dimension dimension dimension LR Axis AP Axis	Representative ellipse perimeter
	1620		88/F	2.53		185	224.90	149.30	599.67
	1630		-/F	2.60		190	204.10	148.44	560.63
	1700	83/M	3.11		185	248.05	166.99	664.26
	1680	74/M	3.71		227.5	230.90	148.10	609.37
	1540		73/F	3.85		240	208.99	157.27	580.97
	1660	77/M	4.04		200	244.14	185.55	681.20
	1800	64/M	4.05		220	244.14	161.13	649.81
	1650	78/M	4.15		247.5	226.85	166.88	625.60
	1740		-/M	4.42		220	221.20	195.68	656.06
	Height	Pulmonary volume Right lung (liter) dimension dimension dimension LR Axis AP Axis	Representative ellipse perimeter
	1650			2.64	224.82	208.07	118.14	531.52
	1680.3		2.81	229.46	212.36	120.58	542.49
	1731.8		3.10	236.92	219.26	124.50	560.12
	1758			3.24	240.54	222.62	126.40	568.69
	1783.1		3.38	243.91	225.74	128.17	576.66
	1807.1		3.51	247.08	228.67	129.84	584.15
	1850			3.75	252.50	233.68	132.68	596.94

Table 2 :

 2

Measures (in mm) of the lungs for the IRSN female phantoms, obtained through mesh measurements, besides the perimeters.

Table 3 :

 3 NEMOSIS applied on 2 phantoms at the boundary of the definition domain. The measured are performed at two

Table 4 :

 4 Simulation of the lung movements on the test patients.

	Phase (%)

Table 5 :

 5 Simulation of the oesophagus and heart movement using NEMOSIS.

Conclusions and perspectives

We have presented an approach to reproducing breath movement on 3D phantoms. Our contribution consists of the training of ANN, the NEMOSIS application, to simulate the position and volume of the lungs, the oesophagus and the heart at every phase of the respiratory cycle. The ANN are trained using organ contours extracted from 4DCT acquisitions of real patients. The interpolation of the lung contours is validated by comparing the ANN results with 4DCT, on 3 patients, not used during network training. The Dice index is used to measure the superposition of interpolated and 4DCT lung contours, and the obtained results (Dice indices greater than 0.94) show the accuracy of the simulation. The aim of this approach is to improve the patient radiation protection by limiting the 4DCT acquisitions, which could be substituted or sustained by respiratory simulation from one 3DCT acquisition.

Nevertheless, the results presented for the heart and the oesophagus cannot yet be validated using the same approach, and adding new data for cardiac and oesophageal motion is a priority. Moreover, when enough data is available, one ANN will be trained for every organ (3 in total). Therefore, it is important to ensure that no superpositions nor collisions exist between the interpolated organ contours. Several solutions may be considered, such as working directly with voxel data. Finally, other movements should be taken into account, such as the heart beat and swallowing, which is complex: 4D TDM are synchronized with breathing but these other movements cannot be acquired using this tool.

Acknowledgements

The authors would like to thank the Besançon university hospital (CHRU) for providing patient data. These works were supported by a funding obtained during a call for projects in the INSERM's PhysiCancer plan, which we would like to sincerely thank. Finally, the authors would like to thank the reviewers for their constructive remarks, and Dr. G. Hruby (Staff specialist, Sidney) for carefully reading this article.


Serge Moulin
email: serge.moulin@univ-fcomte.fr

Nicolas Seux

Stéphane Chrétien

Christophe Guyeux

Emmanuelle Lerat

Simulation-based estimation of branching models for LTR retrotransposons

Motivation: LTR retrotransposons are mobile elements that are able, like retroviruses, to copy and move inside eukaryotic genomes. In the present work, we propose a branching model for studying the propagation of LTR retrotransposons in these genomes. This model allows us to take into account both the positions and the degradation level of LTR retrotransposons copies. In our model, the duplication rate is also allowed to vary with the degradation level. Results: Various functions have been implemented in order to simulate their spread and visualization tools are proposed. Based on these simulation tools, we have developed a first method to evaluate the parameters of this propagation model. We applied this method to the study of the spread of the transposable elements ROO, GYPSY, and DM412 on a chromosome of Drosophila melanogaster.

Introduction

A transposable element (TE) is a DNA sequence able to move from one location to another inside a genome. These sequences, discovered during the 50's by Barbara McClintock (?) exist in almost all living organisms and are the source of a huge number of mutations. They are considered as a major cause of genetic disease in human [START_REF] Belancio | Mammalian non-ltr retrotransposons: for better or worse, in sickness and in health[END_REF] or in Drosophila where they are responsible for more than 80% of the spontaneous mutations [START_REF] Green | Mobile dna elements and spontaneous gene mutation[END_REF]. DNA sequences derived from these TEs can represent a large part of a genome. For example, they represent about 45% of the human genome [START_REF] Lander | Initial sequencing and analysis of the human genome[END_REF] and over 70% of the corn genome [START_REF] Sanmiguel | Evidence that a recent increase in maize genome size was caused by the massive amplification of intergene retrotransposons[END_REF]. Fortunately, most of these sequences correspond to fragments or "dead" elements that have lost their ability to move in the genome due to several lethal mutations or are controlled, especially via epigenetic mechanisms.

TEs have two possible ways to move in a genome, according to their type (?) [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF]. The first class of mobile elements are cut from their original place to move to another one, and are called "DNA transposons" or "Class II transposable elements". The other class of mobile elements, called "retrotransposons" or "Class I transposable elements", use an RNA intermediate to duplicate themselves, the new copy being inserted into another location of the genome. Two orders are identified among the retrotransposons according to the presence or absence of Long Terminal Repeat (LTR) sequences at their extremities. The LTR retrotransposons are similar in structure to retroviruses such as HIV. In both classes, TEs can be classified as either "autonomous", if they encode the enzymes that will allow them to move, or "non autonomous" if they use the enzymes produced by other elements. In an assembled genome, the various sequences corresponding to TE insertions can be found using different bioinformatic approaches (see [START_REF] Lerat | Comparative analysis of transposable elements in the melanogaster subgroup sequenced genomes[END_REF] for a review), which allow us to determine the exact number and positions of each TE insertion. In this article, we focused on the important problem of inferring the history of the spreading of LTR retrotransposons. For this purpose, we modeled the evolution using a branching process where each element (i.e., a copy of a given TE) can randomly evolve via duplication or mutation.

Instances of branching processes have already been proposed in the literature, as putative models for the propagation of TEs. However, most of these studies focus on the evolution of the host population, and not on the propagation of the TEs in the host. The "subject" of these branching models (i.e., the entity able to change or duplicate) is generally the host, while in our case it is the TE itself. For instance, Michael E. Moody [START_REF] Moody | A branching process model for the evolution of transposable elements[END_REF] has used a branching model, where the studied variable was the number of individuals owning i copies of a given TE. Sawyer et al. [START_REF] Sawyer | Distribution and abundance of insertion sequences among natural isolates of escherichia coli[END_REF] produced almost the same model, in order to study the distribution and abundance of insertion sequences. Kaplan et al. proposed a model where TEs can be either of wild type (i.e., non mutated) or of mutant one, which is a little closer to our proposal. When a host gives birth to its child, wild copies can mutate or be deleted, whereas mutant ones can only be removed. New copies can be additionally created. This number of new created copies is supposed to decrease with the proportion of mutants. More recently, interesting models have been proposed that take into account the location of TEs. For instance, Drakos and Wahl [START_REF] Drakos | Extinction probabilities and stationary distributions of mobile genetic elements in prokaryotes: The birth-deathdiversification model[END_REF] suggested a model of mobile promoter evolution, where the probabilities for promoters to duplicate inside or outside their region is potentially not the same.

In the present work, the objective is to propose a new approach for the propagation of LTR retrotransposons that combines a location-dependent model with the fact that LTR retrotransposons can face degradation (i.e., mutations, recombination, etc.), which may decrease their duplication rate, that is, their potentiality to copy and insert elsewhere in the genome. Then, we have developed a first method to evaluate the model parameters: average distance traveled by the TEs before insertion, location of the original copy, average time between two degradations (mutations, recombination, etc.), average impact of a degradation, and the impact of degradations undergone by copies on their duplication speed. This method requires to define a distance between the results of the simulations and the observed chromosome, which is based on the Hungarian method [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. This method has been applied to the spreading of the LTR retrotransposons ROO, DM412, and GYPSY on the chromosome 3L of Drosophila melanogaster. The parameters associated to each TE are computed and a branching tree is proposed in each case. Our results show that, according to our model and method, the roots of ROO, DM412, and GYSPSY on the chromosome 3L could correspond to the annotated copies FBti0059644, FBti0061034, and FBti0062705 respectively.

System and methods

The branching model

2.1.1 The branching tree An example of branching tree is shown in Figure 1.

This branching tree represents the spread of the LTR retrotransposon "ROO" between times 0 and 5.601. At time 0 there is only one copy, called the "root" in this article. At time 1.488, the root duplicates itself to give birth to its first "child".

Finally, the process ends at time T obs = 5.601 with 32 copies. The branching tree represents only the duplications, but copies are also subject to degradations as explained in Section 2.1.2. The state of the tree (i.e., number of TE copies, copy positions, degradations undergone by the copies) at time T obs = 5.601 is named "final state" of this tree. The working principle of our estimation method is to simulate trees, in order to determine in which conditions the final states of simulated trees match well with the observed chromosome. To compare the final state of a simulated tree with the observed chromosome, a distance was build, as detailed in Section 2.2.2.

The general model

The branching model is constructed as follows.

1. The spread starts with only one copy, called the root, at time zero in a location X 0 to be determined. 2. Each copy can either be duplicated or undergo degradations (i.e., mutations, recombinations, etc.) at any time. 3. The number of copies increases due to duplications. When a new copy is created, it receives an index equal to the number of existing copies at the time of its birth, including itself. In the remainder of this article, let T i be the birth date of the i th copy and let τ i,k be the time when the i th copy faces its k th degradation. 4. The time interval τ i,k+1 -τ i,k between two degradations is supposed to be independent and identically distributed (i.i.d.) with an exponential distribution E( 1 µ ), where µ (i.e., average time between two degradations) must be determined. 5. Each copy is also associated to its Needleman-Wunsch [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] similarity to the original state of the root. This similarity is a value between 0 and 1. Let us denote S i (t) the similarity between the original state of the root and the state of the i th copy at time t. This similarity decreases as a function of time, due to degradation effects. In addition, we defined the state of deterioration by D i (t) = 1 -S i (t). 6. At each degradation, the similarity to the root is divided by 1+E( 1 β ), where β has to be determined. In other terms, S i (τ i,k+1 ) =

S i (τ i,k ) 1+E( 1 β )
.

7. At time t, for the i th copy, conditionally on D i (t), we assume that the time before the next duplication follows a distribution E(

1 1+p×D i (t) )
, where p > 0 is a parameter to be determined. In other words, the time before the next duplication is longer when the copy is far from the original state of the root (in terms of Needleman-Wunsch distance). Note that, in this article, "duplication rate" is just the inverse of the "average time before the next duplication". In other words, duplication rate = 1 + p × D i (t). 8. Moreover, each copy is also associated to its position in the chromosome. This position is denoted by X i for the i th copy. This position is constant with time. We assume that each child j of a copy i satisfies X j = X i + χ i,j , where χ i,j follows a distribution U {-1, 1} × E( 1 L ), in which U represents the uniform law (i.e., the probability to choose -1 or 1 is the same) and L is a parameter to be determined. 9. We also take into account the host structure (i.e., position of host genes) to insert or not the child in the chromosome. on the real chromosome. During simulations, when the child moves into a gene, it can be inserted with a probability equal to "density", otherwise the child position is recomputed (if "density" > 1, the child can always insert itself). Furthermore, the child position is also relaunched when it goes outside the chromosome (cf. Part 3.1.3).

Our goal is thus to estimate the parameters of this model, i.e., X 0 , µ, β, p, L and T obs . Note that the duplication speed of the nondegraded root is set to 1 and it does not need to be determined. Indeed, this duplication speed is redundant with µ and p. In addition, note that the parameter L is not really informative about the mean distance traveled by the child before its insertion, due to putative relaunching processes. This is why we also provide the mean traveled distance in simulations, that is, the mean jump, which is denoted by J. See Part 3.1 for an example of its computation.

The estimation method

As explained in Section 2.1.1, the working principle of our estimation method is to simulate trees in order to determine in which conditions the final states of simulated trees match well with the observed chromosome.

Trees are simulated according to the model defined in Section 2.1.2. The stopping criterion of these simulated trees depends on the number of copies in the observed chromosome. Actually, the simulation was constrained to stop at the birth date of the n + 1 th copy, where n is the number of copies in the observed chromosome, see Section 3.1.2 for further details.

In order to improve computation speed, parameters estimation has been split in three parts, i) we estimate µ, β, and L (cf. Section 2.2.1); ii) we estimate X 0 and L (cf. Section 2.2.2), which requires to define a distance between simulated trees and the observed chromosome; iii) we estimate J and T obs (cf. Section 2.2.3).

Estimation of µ, β, and p

The objective here is to estimate the parameters µ, β, and p. We can note that these parameters only affect the distribution of deterioration states. They have no direct influence on copy positions. Thus, the goal at this step is to minimize the differences between the distribution of states of deterioration in the simulated trees, and the states of deterioration distribution in the observed chromosome. More precisely, we want to minimize

D1(T r, C) = n i=1 (D i,T r -D i,C ) 2 where D 1,T r ...D n,T r
is the sorted distribution of the states of deterioration for the simulated tree, D 1,C ...D n,C is the sorted distribution of the states of deterioration for the observed chromosome, and n is the number of TEs in the observed chromosome.

For this purpose, a 3-dimensional grid has been constructed, where each point of this grid represents a triplet (µ, β, p), while X 0 and L are set to predefined values (they do not matter at this stage). A score

S 1 = N 1 i=1 D1(T r i , C
) has been associated to each of these triplets. In this formula, T r i is the i th tree simulated with the parameter set, C represents the observed chromosome, and N 1 is a parameter chosen by the user of the optimization method (for instance, N 1 = 10,000 in the case study of Section 4). The best of these points is selected, and a smaller grid is constructed around it. This iterative process is continued until the precision chosen by the user of the optimization method has been obtained.

Distance between trees, estimation of X 0 and L

The first step, in the estimation of X 0 and L, is to define a distance between the final state of a simulated tree and the observed chromosome. For this purpose, we first need to design a distance between a copy of the simulated tree and a copy of the observed chromosome. Let us name R i the i th copy of the simulated tree, X i its position, and D i (T obs ) its state of deterioration at the end of the process. Similarly, R j is the j th copy of the observed chromosome, X j its position, and D j (T obs ) its state of deterioration. The distance between two copies has been designed as follows:

D2(R i , R j ) = (X i -X j ) 2 w 1 + (D i (T obs )-D j (T obs )) 2 w 2 where w 1 = i=1...n,j=1...n ( X i -X j ) 2 and w 2 = i=1...n,j=1...n ( D i (T obs ) -D j (T obs )) 2
. This weighting by w 1 and w 2 allows us to give the same weight to positions and states of deterioration.

From these distances between two copies, we can now create a matrix of distances W , verifying W ij = D2(R i , R j ). Then, the distance between final states of two trees has been defined as the best possible adjustment between copies, using the so-called Kuhn-Munkres algorithm, also named the Hungarian method [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF][START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. The Hungarian method is an algorithm that allows us to minimize the sum of n elements of a n × n matrix, under the condition that there is only one element by row and only one element by column. In our case, the Hungarian method allows us to assign exactly one copy of the simulated tree to each copy of the real chromosome while minimizing the sum of distances between paired copies. Let us name D3 this distance created by this way.

Once this distance between trees has been defined, we use it to estimate X 0 and L with the same type of process as in the previous step. In other words, a 2-dimensional grid has been constructed, where each point of this grid represents a couple of parameters (X 0 , L) while µ, β, and p are fixed to the values found in the previous stage. The score

S 2 = N 1 i=1 D3(T r i , C)
is then computed for each of these points, and a smaller grid is recursively built around the best point.

Estimation of J and T obs

Conversely to these µ, β, p, X 0 , and L, which are inputted in our simulation algorithm, J (mean jump) and T obs are outputs. It is thus easier to estimate them. In this step, we run N 2 simulations where µ, β, p, X 0 , and L are set to the values found in the two previous steps. The estimations of J and T obs are then the mean results of the output J and T obs of these N 2 simulations (for instance N 2 = 20,000 in the case study of Section 4).

Algorithm

Our proposal has been implemented using Python1 . A short application programming interface is detailed thereafter.

TreeBuild

This main function is used to build branching trees following the model defined in Section 2.1.2. Its halt condition is the targeted number of copies. Its prototype meets the following canvas:

(S, T, T obs ) = T reeBuild(X, µ, β, p, L, n, genes, density),

where n is the desired number of copies, while X 0 , µ, β, p, and L are the model parameters as defined in Section 2.1.2. Moreover genes is the positions of each gene in the observed chromosome, and density is the value defined in Part 2.1.2. Concerning the outputs, S is a n × 1 vector representing states of deterioration while T obs is the propagation time. Finally T is a n × 3 matrix containing, for each copy: its position, its birth date, and the row of its mother, like in Table 1.

In this example, the mother of the copy located in 0.1832 is the root. The mother of the copy located in 0.6644 is the copy located in 0.1832. Other details regarding this main function are provided thereafter. The mean jump J is not a direct output of TreeBuild, but it can be easily computed with T . In this example, J = |0.19-0.5|+|0.18-0.5|+|0.66-0.18|+|0.48-0.19|+|0.14-0.19| 5 .

Multiple clocks management

The working principle of TreeBuild can be summarized as follows: it determines the next event (deterioration or duplication) and executes it until the stopping criterion is satisfied. To determine the next event means to know its nature (deterioration or duplication), its time, and in which of the available copies it happens. Let j be the of available copies at time t 1 . The easiest way to determine the next event is to simulate 2 × j exponential laws, one for each possible deterioration or duplication. The minimum of these 2 × j simulations can thus provide the time, the nature, and the copy related to the next event.

Actually, TreeBuild does not really simulate 2 × j exponential laws, as two properties of this law have been used to shorten computations. Indeed,

∀(p 1 , ..., p 2j ) ∈ R 2j , ∀(Y 1 , ..., Y 2j ) ∼ (E(p 1 ), ..., E(p 2j )), we have: 1. min (Y 1 , ..., Y 2j ) ∼ E   2j i=1 p i   , 2. ∀i ∈ 1...2j, P (Y i = min (Y 1 , ..., Y 2j )) = p i 2j k=1 p k .
Hence, due to the first property, the time of the next event can be simulated by a single exponential law. The second property, for its part, allows us to determine the nature and the copy affected by the next event using a single uniform law.

Stopping criterion

As stated before, the stopping criterion of TreeBuild is related to n (the number of copies of the observed chromosome). But when a chromosome is observed, there is no way to detect that a new duplication has just occurred. Thus, the program cannot stop exactly at the birth of the n th copy. Actually, TreeBuild must run until the T n+1 iteration (the birth date of the n + 1 th copy), and then the propagation time T obs can be determined by:

T obs = Tn+T n+1 2
Furthermore, each value taken by S and T between Tn and T n+1 is kept in memory. Thereby, the values of T and S returned by TreeBuild are values of T and S at time T obs .

The management of copy locations

Copy positions in the chromosome are in the interval [0,1]. The distance traveled by a TE before insertion is assumed to follow an exponential law, but this latter can send the new copy outside the interval [0,1]. The solution chosen in this case is to launch again the computation of the new copy position.

In addition, the copy position is also relaunched with a probability "density" if its position falls into a host gene, as explained in Part 2.1.2. In other words: while (X child / ∈ [0, 1] or (X child ∈ gene and U 1 < density)) :

X child = X mother + U 2 × Y where U ∼ U {-1, 1}, U 2 ∼ U {-1, 1}, and Y ∼ E( 1 L ).

Critical situations

When TreeBuild is launched for each point of the grid of parameters, some critical situations can happen, which may induce a significant slowdown of the program. In particular, when µ is small and β is large, the probability for an event to be a duplication rather than a deterioration becomes very low. Thus, TreeBuild executes an inordinate number of deteriorations before reaching the desired number of copies. To solve this issue, we have decided that when the similarity to the root becomes lower than 0.03, then this copy cannot be degraded anymore.

Estimation method

In the available package, the estimation of the branching model parameters is realized by the Optim function. Its prototype is as follows:

(Best, Score) = Optim(Grid, Case, n1, N 1, N 2, genes).

Here, Grid is a 5 × 4 matrix of settings defined exactly as in Section 2.

Case, for its part, is a 2 × n matrix containing locations and state of deterioration for each copy of the observed chromosome. N 1 and N 2 are settings defined in Sections 2.2.1 and 2.2.3, while n1 indicates how the grid is shrunk at each step after obtaining the best point (cf. the following section). Finally, genes are the positions of each gene in the observed chromosome. The output Best is the parameter set (X 0 , µ, β, p, L, M J, T obs ) returned by the Optim function, while Score is the sum of N 2 differences between simulations and the observed chromosome (this Score is useful if we relaunch Optim several times).

Interval reduction

As explained in Sections 2.2.1 and 2.2.2, the estimation method works with a grid where each point represents a parameter set. When the best point of the grid is found, a new grid is constructed around this point. Note that the new grid is not necessarily included in the previous one, in order to provide a larger degree of freedom of the parameters (in particular, when the latter are close to zero). For instance, in the case of parameter L, the minimum of the new interval is min L min 2 , L best -L delta 2×n1 , where L min and Lmax are the minimum and maximum of the previous interval, L delta = Lmax -L min , L best is the L coordinate of the best parameter set, and n1 is the reduction parameter selected by the user. Thus, the minimum value of the test interval is divided by two at each time the best point of the grid is close enough to zero. The maximum value of the new interval is simply L best -L delta 2×n1 . These formulas, written for L, are also valid for β, µ, and p.

Location in the chromosome

Unlike the other parameters for which we scan a continuous interval, in the case of X 0 , we only consider the positions of the TE copies on the observed chromosome. we scan the interval of integers 1, .., n. However, we process in the same way as with the other parameters (excepted that we use rounded values), and we never get out of the original interval 1, .., n in this case.

For instance, if we choose to test this interval of integers in four points (as we do in our case studies): at the first step, X 0 is tested in the first, round(1 + n-1

3 )-th, round(1 + 2×(n-1) 3

)-th, and nth position of TE copies in the observed chromosome. In the next step, the new test interval thus becomes max(1, round(X best -X delta 2×n1 )...min(n, round(X best -X delta 2×n1 )) where X delta is n -1 here.

Module and package dependencies

The Hungarian method has been applied using the "munkres" module, implemented in 2008 by Brian M. Clapper (Brian M. Clapper, munkres 1.0.7 for Python, https://pypi.python.org/pypi/munkres/).

Results and Discussion

The data

This proposal has been applied to the spread of the LTR retrotransposons ROO, DM412, and GYPSY on the euchromatin part of the chromosome 3L of the Drosophila melanogaster genome. This sequence corresponds to the left arm of the chromosome 3, which is the largest autosomal chromosome of D. melanogaster. This is also the most prolific chromosome for each of the LTR retrotransposons we considered, this is why it has been chosen for this case study. ROO corresponds to the LTR retrotransposon in Drosophila melanogaster with the largest number of copies [START_REF] Kaminker | The transposable elements of the drosophila melanogaster euchromatin: a genomics perspective[END_REF][START_REF] Lerat | Sequence divergence within transposable element families in the drosophila melanogaster genome[END_REF][START_REF] De La Chaux | Evolutionary dynamics of the ltr retrotransposons roo and rooa inferred from twelve complete drosophila genomes[END_REF]. DM412 is supposed to have been recently acquired by D. melanogaster through horizontal transfer from a close relative species [START_REF] Bartolomé | Widespread evidence for horizontal transfer of transposable elements across drosophila genomes[END_REF][START_REF] Lerat | Comparative analysis of transposable elements in the melanogaster subgroup sequenced genomes[END_REF][START_REF] Modolo | A new genome-wide method to track horizontally transferred sequences: application to drosophila[END_REF]. Finally, GYPSY is an older and likely well regulated LTR retrotransposon [START_REF] Lerat | Comparative analysis of transposable elements in the melanogaster subgroup sequenced genomes[END_REF].

Chromosome 3L contains 32 copies of ROO (with a mean nucleotide identity of 68.82%), 16 copies of DM412 (mean nucleotide identity of 60.24%), and six copies of GYPSY (mean nucleotide identity of 13.6%).

Three databases have been used during the experiments. The first one contains positions and nucleotide sequences for each TE copy annotated in D. melanogaster (flybase website 2 version number 5.51 of the D. melanogaster genome [START_REF] Adams | The genome sequence of drosophila melanogaster[END_REF][START_REF] Smith | The release 5.1 annotation of drosophila melanogaster heterochromatin[END_REF]). The second database has been downloaded from the RepBase website 3 and contains the consensus sequences for each TE corresponding to reference elements. The Needleman-Wunsch distance between each TE copy (from 2 http://flybase.org/ 3 http://www.girinst.org/repbase/ the first database) and its reference (from the second database) has been calculated, in order to obtain the deterioration states.

Finally, the third database comes from flybase too. This is the position of all the annotated genes in the euchromatin part of chromosome 3L, in version number 5.51 of the D. melanogaster genome.

In this case study, the estimation method described in Section 2.2 has actually been applied not only once but 40 times in each situation, in order to check the consistency of the obtained parameter sets. The best parameter set of each case study, considering the output "score" (cf. Part 3.2), is presented in Section 4.3. The whole obtained parameter sets are presented in supplementary data with their descriptive statistics. Some indications about consistency of these results are provided in Section 4.4.

Settings

Let us first recall that X 0 , which represents the root position in the chromosome, is inside the interval [0 , 1]. In other words, copy positions have been divided by the chromosome size. For the euchromatic part of chromosome 3L, this size has been set at 24,543,557 base pairs (bp) in the version 5.51. [START_REF] Smith | The release 5.1 annotation of drosophila melanogaster heterochromatin[END_REF].

In Table 2, each row represents the beginning and the end of the test interval, the number by which the test interval has been divided, and the final desired accuracy regarding the parameter. In particular, in the third column, the value associated to each is 3. This latter means that these parameters have been tested at the beginning, in the first third, in the second one, and at the end of the test interval.

Finally, at each iteration, the grid used for µ, β, and p estimations contains 4 3 = 64 points while, at each iteration, the grid used to estimate X 0 and L contains 4 2 = 16 points.

The other parameters are:

• n1 = 1.5: at each step, after the best point has been found, the grid's dimensions have been divided by 1.5.

• N 1 = 10,000: each point has been tested 10,000 times during estimation of X 0 , µ, β, p, and L.

• N 2 = 20,000: J and T obs estimations are the average values of 20,000 simulations. The output "score" is computed based of these simulations.

Furthermore, we can notice that in the case of ROO, 18 of the 32 copies are located inside genes while genes hold 72.7% of the studied part (euchromatin). Thus, density = 18 72,7 14 27,2 = 0.482. In the cases of DM412 and GYPSY, the densities are respectively equal to 0.625 ( 10 16 ) and 0.187 ( 2 6 ).

Results

The obtained parameters are summarized in Table 3.

If we consider for instance the spread of ROO, the obtained parameters can be interpreted as follows:

• X 0 = 29.
The root is on the 29 th position in chromosome 3L. This is the copy FBti0059644, which is located between the 21, 954, 331 th and the 21, 954, 698 th nucleotide. • µ = 2.396. The average time between two degradations is 2.396, where 1 is the average time before duplication of the root. Degradations are thus less frequent than duplications. Please note that this estimation of µ is without time unit: it is related to the duplication speed of the root. It allows us to estimate duplication speed when the deterioration speed is known, and vice versa. • β = 0.351. Each degradation causes a division by 1 + E( 1 0.351 ) of the similarity. Thus the similarity is divided by 1.351 on average at each degradation.

• p = 0.051. p allows us to determine how many degradations led to a decrease in the duplication speed. For example, in this case, if the identity between a copy and the reference is 0.75 (i.e., state of deterioration = 0.25), then the duplication speed of this copy is reduced by 1.275% (indeed 0.25 × 0.051 = 0.01275). This looks like a low effect. • L = 1.331. The distance traveled by the TE before insertion follows a distribution E( 1 1.331 ), with the relaunching process: (1) when the position of the child is out of the chromosome (cf. Section 3.1.3) and (2) possibly when the position of the child is inside a gene. As explained in Part 2.1.2, the parameter J is better to represent the average distance before insertion. However, we can notice that L = 1.331 is a pretty large number (larger than 1). This value tends to suggest that the child inserts itself more or less anywhere in the case of ROO (compared to DM412 or GYPSY).

• J = 0.300. The mean distance traveled by the TE before insertion is 0.300. • T obs = 4.070. The time of ROO propagation is 4.070 larger than the non-deteriorated root. In addition, 4.070/2.396 = 1.699, thus the root has faced 1.699 degradations on average. From a global perspective, we can note that in each case, the root corresponds to the border elements. They are in position 29 th over 32, 14 th over 16, and 5 th position over 6 respectively for ROO, DM412, and GYPSY. The latter correspond to positions 0.894, 0.965, and 0.969 when the chromosome is considered as a [0,1] interval. This can be due to a larger density of TEs in this area. In addition, we can notice that p is really close to 0 in each case. This should imply no (linear) effect of the degradation on the duplication rate. This is an unexpected result. Indeed, the fact that the degradation undergone by the copies reduces their ability to duplicate sounds natural. This is why the parameter p was added in the model. Finally, we can notice that in the case of ROO, the results suggest a few big degradations, while in the case of GYPSY, they suggest a lot of little ones.

The fact that some of the obtained parameters are outside the test interval chosen at the beginning of the program (for instance, L = 1.331 in the case of ROO or µ = 0.050 in the case of GYSPY) is a desired effect, to let a larger freedom to the parameters. In particular, the aim was to let parameters to be as close as possible to zero if required (cf. Section 3).

In each of these three cases, one billion trees have been simulated with the obtained parameter set. The best of these trees is shown in Figure 1 for ROO and in supplementary data for DM412 and GYPSY.

Focusing on the roots

According to these models and methods, the root of ROO could correspond to the FBti0059644 copy. This is an incomplete copy (368 bp, compared to the reference which has 9,112 bp) that corresponds to a solo-LTR, a remnant from a LTR-LTR recombination. This copy is thus no longer active, but it is quite recent since its divergence to the reference is rather low (95.71% of identity on the aligned part of the sequences).

The root of DM412 could correspond to the FBti0061034 copy. This is a very degraded copy that has 88 bp in length, corresponding to an internal portion of the reference element, whose length is 7440 bp. The copy is old since it is very divergent compared to the reference (80.90% of identity).

Finally, the root of GYPSY could correspond to the FBti0062705 copy. This copy is an incomplete and very degraded element of 1,282 bp length (7,471 bp for the reference) with a very high divergence to the reference (70.59% identity). This copy corresponds to a piece of the inner part of the gypsy element, and this is a very old and degraded copy that is not currently active.

Consistency of results

As explained previously, the optimization method has been actually applied 40 times for each TE. The descriptive statistics for these three cases are summarized in the supplementary data. In addition, for each parameter that has been estimated by scanning an interval (i.e., X 0 , µ, β, p, and L), quotients Standard deviation of the results Test interval have been computed in each case, in order to assess the consistency of the results. These quotients are reproduced in Table 3. Standard deviations of the output J and T obs are also set in the same table.

Several versions of the code have been implemented in this project in order to increase the consistency of these results. In a previous version, X 0 , µ, β, p, and L were estimated all together (i.e., steps 2.2.1 and 2.2.2 were merged). This approach implied to work in a 4 5 = 1024 points grid, thus, it provided a lower number of trials by point, which reduced consistency. In our very first version, the duplication rate was a parameter to be estimated while T obs = 1 was the stopping criterion of TreeBuild. The difference between n and the number of copies at the end of the simulation (which was not necessarily n in this setup) was penalized. However, in this setup, the duplication rate was the only parameter to be consistent.

Finally, in the setup presented here, the consistency looks acceptable in most of the cases (i.e., excepted for L in the case of ROO or for µ in the case of GYPSY). Nevertheless, this requires a large number of simulations implying that the program runs several days.

Conclusion and future perspectives

In this article, a model has been proposed for the propagation of LTR retrotransposons in a genome. Various functions have been implemented to simulate this spread as well as graphic representations. Finally, a first method for estimating the parameters of this propagation model has been proposed and applied to the spread of TEs corresponding to the ROO, GYPSY, and DM412 elements in a chromosome of Drosophila melanogaster. However, this work can be improved in various directions, some of them being listed below.

The first point is that the model of propagation should be applied to the full genome instead of a single chromosome of Drosophila melanogaster. Indeed, a copy inserted in a given chromosome can produce a child that will not necessarily inserts itself in the same chromosome. An idea to extend the model to the full genome could be to let the position of the child copy following a δ(a)U + δ(1 -a)E( 1L ) law. In other words, the child copy inserts itself anywhere in the genome (uniform distribution), with a probability a, or it inserts itself on the same chromosome than the mother copy (with the distance to the mother following an exponential law) with probability 1 -a, where a is a new parameter to determine. Nevertheless, this approach raises various questions that must be answered in a further study. Firstly, should E( 1L ) represent a number of nucleotides or a proportion of the considered chromosome? How to redefine the distance built in Part 2.2.2? etc.

Secondly, as explained in Part 4.4, the estimation method needs a large compilation time to give acceptably consistent results. This compilation time increases quickly when the sample size (number of TEs) increases or if we want to estimate more parameters. Therefore, we will test other ways (likelihood maximization, neighbor joining, or Bayesian estimation) to improve our proposal. These types of methods would also allow us to produce some confidence intervals for the estimated parameters. Our long-term objective is to create a useful tool for estimating consistently both parameters and the branching process itself. In other word, our goal is to produce a tool close to phylogenetic tree estimation but adapted to TE constraints. Another possibility of improvement could be to consider the possibility of several roots. For instance, a method of unsupervised classification like Gaussian Mixture model could be applied in order to detect the number of clusters.

In this project, we used RepBase consensus sequences, based on a lot of TE sequences as root sequence. Another possibility could also be to produce an ancestral reconstruction based only on the sequences of the case study. In this way, it would not be necessary to search the data in two different databases. We can also note that in this work, we consider all modifications (i.e., mutations, LTR recombination, and so on) as one single and global deterioration effect. It could be interesting to try to distinguish each effect. Finally the effect of an epigenetic regulation that can affect TE behaviour even if they do not face sequence degradation could be taken into account.

Fig. 1 .

 1 Fig. 1. ROO spread

Table 1 .

 1 Example of the output T

	[[ 0.5	0.	0. ]
	[ 0.19031606 1.83699228 0. ]
	[ 0.18321005 11.25706728 0. ]
	[ 0.66442132 17.61532334 2. ]
	[ 0.48479738 25.45993783 1. ]
	[ 0.13876928 28.11662473 1. ]]

Table 2 .

 2 Setting tableparameter starting point end point interval division desired accuracy

Table 3 .

 3 Results and consistency

	Best parameter sets
	parameter	ROO DM412 GYPSY
	X0	29	14	5
	µ	2.396 1.135	0.050
	β	0.351 0.235	0.093
	p	0.051 0.001	0.007
	L	1.331 0.336	0.016
	J	0.300 0.216	0.013
	T obs	4.070 3.353	2.379
	Consistency indicators
	parameter	ROO DM412 GYPSY
	X0	0.104 0.096	0.392
	µ	0.020 0.017	0.002
	β	0.040 0.064	0.044
	p	0.001 0.000	0.000
	L	0.262 0.036	0.033
	J	0.010 0.013	0.024
	T obs	0.042 0.014	0.031

Available at https://github.com/SergeMOULIN/retrotransposonsspread

Funding

Computations have been performed on the supercomputer facilities of the Mésocentre de calcul de Franche-Comté. This work has been supported by the Transposable Elements project of the Franche-Comté region. The work of S.C. was also sponsored by NMS/IRD project 117480 at NPL.

Our proposal has been implemented using Python software. Source code is freely available on the web at https://github.com/SergeMOULIN/retrotransposons-spread.


Jean-Luc Baril
email: barjl@u-bourgogne.fr

Richard Genestier
email: richard.genestier@femto-st.fr

Alain Giorgetti
email: alain.giorgetti@femto-st.fr

Armen Petrossian
email: armen.petrossian@u-bourgogne.fr

Rooted planar maps modulo some patterns

Keywords: Rooted planar map, shuffle of Dyck words, equivalence relation, pattern

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and notations

A map is defined topologically as a 2-cell imbedding of a connected graph, loops and multiple edges allowed, in a 2-dimensional surface. A rooted planar map is a map on a sphere with a distinguished edge, called the root, assigned with an orientation. Figure 1(a) shows a rooted planar map with seven edges. We refer to [START_REF] Arquès | Une relation fonctionnelle nouvelle sur les cartes planaires pointées[END_REF][START_REF] Tutte | A census of planar maps[END_REF][START_REF] Tutte | On the enumeration of planar maps[END_REF] for the enumeration of rooted planar maps with respect to the number of edges. Planar maps have been widely studied for their combinatorial structure and their links with other domains such as theoretical physics where they appear in some models of 2-dimensional quantum gravity for instance. From a combinatorial point of view, it is proved in [START_REF] Lehman | Counting rooted maps by genus II[END_REF] that rooted planar maps are in one-to-one correspondence with shuffles of Dyck words that avoid a specific pattern. Then, it becomes natural to extend to maps statistical studies on Dyck words (see [START_REF] Deutsch | Dyck path enumeration[END_REF][START_REF] Mansour | Statistics on Dyck words[END_REF][START_REF] Merlini | Some statistics on Dyck words[END_REF][START_REF] Mansour | Counting peaks at height k in a Dyck path[END_REF][START_REF] Mansour | Counting Dyck paths according to the maximum distance between peaks and valley[END_REF][START_REF] Sapounakis | Counting strings in Dyck paths[END_REF] for instance). It is one of the main purpose of this paper.

A Dyck word on the alphabet {a, ā} is a word generated by the contextfree grammar S → ε | aSāS where ε is the empty word. Let D a be the set of all Dyck words on the alphabet {a, ā}, and let D b be the set of all Dyck words on the alphabet {b, b}. It is well known that Dyck words of semilength n are counted by Catalan numbers (A000108 in the on-line encyclopedia of integer sequences [START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF]), and that any non-empty Dyck word w ∈ D a has a unique decomposition of the form w = aαāβ where α and β are two Dyck words in D a (see [START_REF] Deutsch | Dyck path enumeration[END_REF]). Also, a word w on the alphabet {a, ā} belongs to D a if and only if the following two properties hold: (i) the number of letters a is equal to the number of letters ā in w, and (ii) in any prefix of w, the number of letters a is greater than or equal to the number of letters ā.

We say that an occurrence of the letter a matches an occurrence of ā located to its right in w ∈ D a , whenever the subword of w consisting of all letters between these two occurrences also belongs to D a . In this case, the pair (a, ā) is called a matching in w. For instance, if w = aāaaāaaāāā, then the second occurrence of the letter a matches the last occurrence of ā since aāaaāā is a Dyck word.

A shuffle of two Dyck words v ∈ D a and w ∈ D b is a word s on the alphabet {a, ā, b, b} where s is constructed by interspersing the letters of v and w. Let S be the set of all shuffles of two Dyck words of D a and D b . Shuffles of semilength n, n ≥ 0, in S are enumerated by the sequence A005568 in [START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF] A shuffle s of two Dyck words v ∈ D a and w ∈ D b will be called crossing whenever there exists a matching (a, ā) in v and a matching (b, b) in w such that s can be written s = αbβaγ bδāη where α, β, γ, δ and η belong to {a, ā, b, b} * . Then the occurrence ba bā will be called a pair of crossing matchings. Let N CS ⊂ S be the subset of non-crossing shuffles in S,

(b) (a)
aaabbāā bā bab bā Figure 1: A rooted planar map with 7 edges and its non-crossing shuffle representation.

i.e., shuffles with no pair of crossing matchings. For instance, aāb baāb baā is in N CS, while aābb ba bā is not in N CS. Notice that non-crossing shuffles are called canonical parenthesis-bracket systems in [START_REF] Walsh | Combinatorial enumeration of non-planar maps[END_REF]. The shuffles of semilength n ≥ 0 in N CS are enumerated by the sequence A000168 in [START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF] whose first values for n ≥ 0 are 1, 2, 9, 54, 378, 2916, 24057, 208494. They are in one-to-one correspondence with the rooted planar maps with n edges [START_REF] Bernardi | Bijective counting of tree-rooted maps and shuffles of parenthesis systems[END_REF][START_REF] Cori | Enumeration des graphes planaires à l'aide des séries formelles en variables non commutatives[END_REF][START_REF] Lehman | Counting rooted maps by genus II[END_REF][START_REF] Walsh | Combinatorial enumeration of non-planar maps[END_REF]. See Figure 1 for an example of rooted planar map with its representation as a non-crossing shuffle in N CS. This one-to-one correspondence is obtained by the following process. Starting with the root edge, we follow or cross all edges of the map by making its tour in counter-clockwise direction. Each edge in the map must be reached twice. If the final vertex of the encountered edge has not yet been considered, then we follow this edge and we write the letter a; otherwise, if the edge is reached for the first time, then we write the letter b and we cross it; in the other cases the edge is reached for the second time, and we write ā (resp. b) and we follow the edge (resp. we cross the edge) whenever the first assignment of the edge was the letter a (resp. b). Table 1 gives the correspondence between some patterns of length at most two in a shuffle with their meaning in terms of map. For instance, a pattern b b in a shuffle of N CS corresponds to an empty loop on a vertex of the corresponding rooted planar map. Also, the pattern a b cannot occur in the shuffle representation of a rooted planar map because it necessarily creates a pair of crossing matchings.

In a recent paper [START_REF] Baril | Equivalence of Dyck paths modulo some statistics[END_REF] (for equivalence classes of permutations see [START_REF] Baril | Equivalence classes of permutations modulo excedances[END_REF]), the authors introduced an equivalence relation on the set of Dyck words where two Dyck words are equivalent whenever the positions of the occurrences of a given pattern are the same in both words. They provided generating functions for the numbers of equivalence classes whenever the patterns are of length two: aa, aā, āā and āa. See also [START_REF] Baril | Equivalence classes of Motzkin paths modulo a pattern of length at most two[END_REF] for a study of this equivalence two shuffles with the same semilength are π-equivalent whenever they coincide on the positions of occurrences of the pattern π.

For instance, the shuffle s = aaba bābā bā is bā-equivalent to s = aāab bāab bā since the occurrences of bā appear in positions 5 and 9 in s and s .

In Section 2, we show that the problem of the enumeration of the πequivalence classes in N CS is the same as in S whenever π is a pattern of length at most two that does not belong to the set { ba, bā, a b, āb}. In Sections 3-6, we present enumerative results by providing generating functions for the number of π-equivalence classes when π / ∈ { ba, bā, a b, āb}. Using the one-toone correspondence between rooted planar maps and non-crossing shuffles of Dyck words, this induces enumerative results for equivalence classes of rooted planar maps relatively to the positions of some patterns. See Table 2 for an overview of these results. Notice that the pattern a b does not appear in any shuffle of a rooted planar map, and that we did not succeed to obtain the number of π-equivalence classes in S and N CS for π ∈ {āb, ba, bā}. So, we leave these cases as open problems.

Pattern Sequence Sloane a n , 1 ≤ n ≤ 9 {a}, {ā}, {b}, { b} 1 √ 1-4x
Central binomial coeff.

A000984

2, 6, 20, 70, 252, 924, 3432, 12870, 48620

{aā}, {b b} 1-x 1-3x+x 2

Shift of A001519

2, 5, 13, 34, 89, 233, 610, 1597, 4181

{āa}, { bb} 1-2x 1-3x+x 2 A001519
1, 2, 5, 13, 34, 89, 233, 610, 1597

{aa}, {āā}, {bb}, { bb } 1+x- √ 1-2x-3x 2 2x 2 +3x-1+ √ 1-2x-3x 2 Bisection of A191385 1, 2,

Some preliminary results

In this section, we prove that for some specific patterns π the numbers of π-equivalence classes in S and N CS are the same.

Lemma 1 Let w = αāβγ be a Dyck word in D a with βγ in {a, ā} * . Then the word w = αβāγ is also in D a .

Proof. The Dyck word w is obtained from w by a shift on the right of a letter ā. So, the number of a and the number of ā remains unchanged in w . Moreover, in any prefix of w the number of a is greater than or equal to the number of ā and this property remains satisfied in w . Using the characterization of a Dyck word given in the introduction, w is in D a . 2

Lemma 2 Let w = αaβāγ be a Dyck word in D a . If β = β 1 β 2 . . . β k is in D a then for any i, 1 ≤ i ≤ k, the word w = αaβ 1 . . . β i-1 āβ i . . . β k γ is also in D a .
Proof. Since w is in D a , (i) the number of a and the number of ā are equal in w, and in any prefix (ii) the number of a is greater than or equal to the number of ā. Since β is in D a , the number of a in any prefix of αaβ 1 . . . β i-1 is strictly greater than the number of ā. This means that αaβ 1 . . . β i-1 ā satisfies Condition (ii). Then, w satisfies (ii) which implies that w belongs to D a .

2

Of course, the last two lemmas also hold if we replace the letter a with b. So, we use them indifferently for a and b. Theorem 1 Let s be a shuffle in S and π be a pattern of length at most two not in {a b, āb, ba, bā}. Then, there exists a shuffle s in N CS with the same semilength as s so that s and s are π-equivalent.

Proof. Let s be a shuffle in S that does not lie in N CS, i.e. s contains a pair of crossing matchings. Then, s can be decomposed s = αbβaγ bδāη so that α, β, γ, δ and η are in {a, ā, b, b} * and the crossing matching ba bā is chosen the leftmost possible in s.

We distinguish three cases: (i) π ∈ {a, aa, ab, ba}; (ii) π = aā; and (iii) π = āa. The other cases are easily obtained using classical symmetries on shuffles ( a ↔ b and mirror). Case (i): Let us consider the word s = αbβaγāδ bη obtained from s by an exchange of ā and b. With Lemma 1 and Lemma 2, s is also a shuffle of two Dyck words. The two shuffles s and s belong to the same π-equivalence class since this operation cannot create or delete a pattern in {a, aa, ab, ba}. Moreover the leftmost pair of crossing matchings in s is shifted on the right compared to s. We repeat the process for as long as required, and we obtain a shuffle in N CS with the same semilength as s, which is π-equivalent to s. Case (ii): Let us assume that s = s 1 s 2 . . . s 2k for some k ≥ 1 with s i ∈ {a, ā, b, b}, 1 ≤ i ≤ 2k. We obtain s from s by the following process: we preserve the positions of all occurrences of the pattern aā; the other letters of s are chosen so that the restriction of s to them is the Dyck word of the form b b for some ≥ 0. It is straightforward to see that s ∈ N CS and that the two shuffles s and s belong to the same aā-equivalence class. Case (iii): Let us assume that s = s 1 s 2 . . . s 2k for some k ≥ 1 with s i ∈ {a, ā, b, b}, 1 ≤ i ≤ 2k. We obtain s from s by the following process: we preserve the positions of all occurrences of the pattern āa and we set s 1 = a and s 2k = ā; the other letters of s are chosen so that the restriction to them is the Dyck word of the form a ā for some ≥ 0. It is straightforward to see that s ∈ D a ⊂ N CS and that the two shuffles s and s belong to the same āa-equivalence class.

2

Corollary 1 Let π be a pattern of length at most two not in { ba, bā, a b, āb}.

For each word length the number of π-equivalence classes in S is also the number of π-equivalence classes in N CS.

This corollary allows us to limit our study of the general case to shuffles in S and patterns π not in { ba, bā, a b, āb}. All results in the following sections will hold for both sets S and N CS.

3 Equivalence modulo π ∈ {a, b, ā, b}

The results for π ∈ {b, ā, b} are deduced from the ones for π = a by using the classical symmetries (mirror and a ↔ b) on shuffles. So, we set π = a in this section.

Let A be the set of shuffles in S defined by the grammar A → D a | D a bD a bA where D a is the set of Dyck words on the alphabet {a, ā}.

Lemma 3

The set A and the set of a-equivalence classes of S are in lengthpreserving one-to-one correspondence.

Proof. For k ≥ 1 let s = s 1 s 2 . . . s 2k be a shuffle in S\A with s i ∈ {a, ā, b, b} for 1 ≤ i ≤ 2k. Let us prove that there exists a shuffle s ∈ A (with the same semilength as s) such that s and s belong to the same a-equivalence class. We define the word s by performing the following process on s for i from 1 to 2k:

-If s i = a then we set s i = a; -Otherwise, we distinguish three cases: Since the process preserves the occurrences of the letter a in s and do not introduce other letters a in s , s and s belong to the same a-equivalence class. Moreover, s necessarily lies in A. Indeed, the process sets s i = ā whenever w a (s 1 . . . s i-1 ) is not a Dyck word of D a ; in all other cases, it sets s i = b or s i = b alternatively. If s does not contain any occurrence of b, then s ∈ D a ⊂ A; otherwise, the previous construction ensures that s has a prefix of the form αbβ b with α and β in D a . This means that s satisfies the grammar A → D a | D a bD a bA. Thus, we have s ∈ A. Now, it suffices to prove that two distinct shuffles s and s with the same semilength in A are not a-equivalent. For a contradiction, let us assume that s and s belong to the same class.

(i) If w a (s 1 . . . s i-1 ) / ∈ D a then we set s i = ā; (ii) If w a (s 1 . . . s i-1 ) ∈ D
We distinguish two cases:

-If s ∈ D a or s ∈ D a with the same positions of the occurrences of a then we have s = s .

-Otherwise, s can be written s = αbβ bγ with α and β in D a and γ in A. Also, s can be written s = α bβ bγ with α and β in D a and γ in A.

Since α and α belong to D a and have the same positions of occurrences of the letter a, we necessarily have α = α . The same argument for β implies that β = β . We complete the proof by induction on the length for γ and γ in A and we conclude s = s .

2

Theorem 2 The generating function for the set of a-equivalence classes in S (or in N CS) with respect to the semilength is given by

1 √ 1 -4x = n≥0 2n n x n .
Proof. By Lemma 3, it suffices to provide the generating function A(x) of the set A, with respect to the semilength. Every nonempty shuffle s ∈ A is obtained by the grammar A → D a | D a bD a bA where D a is the set of Dyck words on the alphabet {a, ā}. This induces the functional equation:

A(x) = D(x) + xD(x) 2 A(x)
where

D(x) = 1- √ 1-4x 2x
is the Catalan generating function for the set D a , which implies the required result. 2

4 Equivalence modulo π ∈ {aā, b b, āa, bb}

The results for π ∈ {b b, bb} are deduced from the ones for π = aā and π = āa by using the classical symmetries. Then, we study in this section the two cases π = aā and π = āa. Let B be the subset of S defined by the grammar B → aāB | bB bB | ε where B → aāB | ε. It is worth noticing that the set B consists of words of the form (aā) k for all k ≥ 0. Moreover, it is straightforward to see that for any s ∈ B we have w b (s) = (b b) i for some i ≥ 0 and w a (s) = (aā) j for some j ≥ 0. Finally, if X = aā then the set B consists of all shuffles of words X i , i ≥ 0, with Dyck words of the form (b b) j , j ≥ 0, where the occurrences of X = aā are never split. For instance, baāaā baāb b belongs to B, while baāa bāaāb b / ∈ B.

Lemma 4

The set B and the set of aā-equivalence classes of S are in lengthpreserving one-to-one correspondence.

Proof. For k ≥ 1 let s = s 1 s 2 . . . s 2k be a shuffle in S\B with s i ∈ {a, ā, b, b} for 1 ≤ i ≤ 2k. Let us prove that there exists a shuffle s ∈ B (with the same semilength as s) such that s and s belong to the same class. We define the word s by performing the following process on s: -all occurrences of the pattern aā in s are preserved in s ; -for i from 1 to 2k such that s i does not belong to a pattern aā in s, we distinguish two cases:

(i) if w b (s 1 . . . s i-1 ) = (b b) j for some j ≥ 0, then we set

s i = b, (ii) if w b (s 1 . . . s i-1 ) = (b b) j b for some j ≥ 0, then we set s i = b. For all i (1 ≤ i ≤ 2k) we have either w b (s 1 . . . s i ) = (b b) j or w b (s 1 . . . s i ) = (b b) j b
for some j ≥ 0. Moreover, we have w a (s 1 . . . s i ) = (aā) for some ≥ 0 so that a and ā appear necessarily adjacent in s in an occurrence aā. This means that s is a shuffle that lies in B. Since the process preserves the occurrences of the pattern aā, s and s belong to the same aā-equivalence class.

Less formally, s is obtained from s by fixing the occurrences of aā and by replacing all other letters by a Dyck word of the form (b b) j for some j ≥ 0. Therefore, the definition of B ensures us that s is the unique element in B that lies into the same aā-equivalence class as s.

2

Theorem 3 The generating function for the set of aā-equivalence classes in S with respect to the semilength is given by

1 -x 1 -3x + x 2 = n≥0 F 2n+1 x n
where F m is the m-th Fibonacci number (defined by F m = F m-1 +F m-2 with F 0 = 0 and F 1 = 1).

Proof. By Lemma 4, it suffices to provide the generating function B(x) of the set B, with respect to the semilength. Every shuffle s ∈ B is obtained by the grammar B → aāB | bB bB | ε where B → aāB | ε. This induces the functional equation:

B(x) = 1 + xB(x) + xB(x)B (x)
where B (x) = 1 + xB (x) is the generating function for the set B .

A simple calculation provides the result. 2

Now we consider the equivalence relation for the pattern π = āa. Let B be the set B = aBā | ε where B is defined at the beginning of this section.

Lemma 5

The set B and the set of āa-equivalence classes of S are in length-preserving one-to-one correspondence. Proof. Let s = s 1 s 2 . . . s 2k (k ≥ 1) be a nonempty shuffle in S\B with s i ∈ {a, ā, b, b} for 1 ≤ i ≤ 2k. Let us prove that there exists a shuffle s ∈ B (with the same semilength as s) such that s and s belong to the same āa-equivalence class.

We define the shuffle s as follows:

-we set s 1 = a and s 2k = ā; -all occurrences of the pattern āa in s are preserved in s ; -for i from 2 to 2k -1 such that s i does not belong to a pattern āa, we distinguish two cases:

(i) if w b (s 1 . . . s i-1 ) = (b b) j for some j ≥ 0, then we set

s i = b, (ii) if w b (s 1 . . . s i-1 ) = (b b) j b for some j ≥ 0, then we set s i = b.
For all i, 1 ≤ i ≤ 2k we have either w b (s 1 . . . s i ) = (b b) j or w b (s 1 . . . s i ) = (b b) j b for some j ≥ 0. Moreover, we have w a (s 1 . . . s i ) = (aā) for some ≥ 0 so that a and ā appear necessarily adjacent in s in an occurrence āa, except for s 1 = a and s 2k = ā. This means that s belongs to B . Since the process preserves the occurrences of the pattern āa, s and s belong to the same āa-equivalence class.

This process ensures us that s is the unique element in B that lies into the same āa-equivalence class as s.

2

Theorem 4 The generating function for the set of āa-equivalence classes in S with respect to the semilength is given by

1 -2x 1 -3x + x 2 = n≥0 F 2n-1 x n
where F m is the m-th Fibonacci number (with F -1 = 1).

Proof. By Lemma 5, it suffices to provide the generating function for

B = aBā | ε which is B (x) = xB(x) + 1 = 1-2x 1-3x+x 2 . 2
5 Equivalence modulo π ∈ {aa, āā, bb, bb }

The results for π ∈ {āā, bb, bb } are deduced from the ones for π = aa by using the classical symmetries.

Let C be the subset of S defined by the grammar C → C | C bC bC where C → aaāC āC | a(C \ε)āC | ε. In fact, the set C consists of Dyck words in D a such that every occurrence of the letter a is contiguous with another occurrence of a.

Lemma 6

The set C and the set of aa-equivalence classes of S are in lengthpreserving one-to-one correspondence. Proof. Let s = s 1 s 2 . . . s 2k (k ≥ 1) be a shuffle in S\C with s i ∈ {a, ā, b, b} for 1 ≤ i ≤ 2k. Let us prove that there exists a shuffle s ∈ C (with the same semilength as s) such that s and s belong to the same class. We define the word s by performing the following process on s:

-all occurrences of the pattern aa in s are preserved in s ; -for i from 1 to 2k such that s i does not belong to a pattern aa, we distinguish three cases:

(i) If w a (s 1 s 2 . . . s i-1 ) ∈ D a and w b (s 1 s 2 . . . s i-1 ) = (b b) j b for some j ≥ 0, then we set s i = b;

(ii) If w a (s 1 s 2 . . . s i-1 ) ∈ D a and w b (s 1 s 2 . . . s i-1 ) = (b b) j for some j ≥ 0, then we set s i = b;

(iii) If w a (s 1 s 2 . . . s i-1 ) / ∈ D a then s i = ā. For all i, 1 ≤ i ≤ 2k we have either w b (s 1 s 2 . . . s i ) = (b b) j or w b (s 1 s 2 . . . s i ) = (b b) j b for some j ≥ 0; w a (s 1 s 2 . . . s i ) is a prefix of a Dyck word in D a , w a (s 1 s 2 . . . s 2k ) ∈ D a and w b (s 1 s 2 . . . s 2k ) = (b b) j for some j ≥ 0. So, the word s is a shuffle in S. Moreover, any occurrence of the letter a in s is always contiguous with another occurrence of a. Let i 1 ≥ 1 (resp. i 2 > i 1 ) be the position of the leftmost b (resp. b) in s ; then the prefix w a (s 1 s 2 . . . s i 1 -1 ) (resp. w a (s 1 s 2 . . . s i 2 -1 )) lies into D a which implies that the shuffle s is of the form αbβ bγ with α and β in C and γ in C; thus we have s ∈ C.

The process preserves the positions of the occurrences of aa. So, s and s belong to the same aa-equivalence class.

The proof for the unicity of s in C is obtained mutatis mutandis from the proof of Lemma 3.

2

Theorem 5 The generating function for the set of aa-equivalence classes in S with respect to the semilength is given by

1 + x - √ 1 -2x -3x 2 2x 2 + 3x -1 + √ 1 -2x -3x 2 .
This sequence corresponds to the values of even ranks in A191385 ( Let C (x) be the generating function for the set C . From the above grammar, we deduce the functional equation:

C (x) = x 2 C (x) 2 + x(C (x) -1)C (x) + 1 and C(x) = C (x) + xC (x) 2 C(x) which provides C (x) = x+1- √ 1-2x-3x 2 2x(x+1) and C(x) = 1+x- √ 1-2x-3x 2 2x 2 +3x-1+ √ 1-2x-3x 2 . 2 6
Equivalence modulo π ∈ {ab, ba, āb , bā}

The results for π ∈ {ba, āb , bā} are deduced from the ones for π = ab by using the classical symmetries.

Let E be the set of shuffles in S defined by the grammar E → E | E bE bE where E → abE bE āE | ε. In fact, the set E consists of shuffles s in S such that every occurrence of a letter x ∈ {a, b} appears in a pattern ab, and there is no matchings (a, ā) and (b, b) such that abā b appears in s.

Lemma 7

The set E and the set of ab-equivalence classes of S are in lengthpreserving one-to-one correspondence. Proof. Let s = s 1 s 2 . . . s 2k (k ≥ 1) be a shuffle in S\E with s i ∈ {a, ā, b, b} for 1 ≤ i ≤ 2k. Let us prove that there exists a shuffle s ∈ E (with the same semilength as s) such that s and s belong to the same class. We define the shuffle s by performing the following process on s:

-all occurrences of the pattern ab in s are preserved in s ; -for i from 1 to 2k such that s i does not belong to a pattern ab, we distinguish two cases:

(i) If w b (s 1 s 2 . . . s i-1 ) ∈ D b and w a (s 1 s 2 . . . s i-1 ) ∈ D a , then we set s i = b; (ii) If w b (s 1 s 2 . . . s i-1 ) / ∈ D b or w a (s 1 s 2 . . . s i-1
) / ∈ D a then there is a letter x ∈ {a, b} which is not matched in the prefix s 1 s 2 . . . s i-1 (we choose x the rightmost possible). If x = a, we set s i = ā, otherwise we set s i = b. Now, let us prove that this process produces a shuffle in E. We distinguish two cases.

-Let us assume that any letter x ∈ {a, b} belongs to an occurrence ab in s . If s i / ∈ {a, b} then s i = ā (resp. s i = b) whenever the rightmost nonmatched x ∈ {a, b} in s 1 s 2 . . . s i-1 is x = a (resp. x = b). So, s 1 s 2 . . . s i-1 s i can be written either abα bβā or abα b according to the value of s i (ā or b), such that w a (α) ∈ D a , w a (β) ∈ D a , w b (α) ∈ D b and w b (β) ∈ D b . Moreover, α and β also satisfy the property that any x ∈ {a, b} belongs to an occurrence ab. Then, s can be generated by the grammar E → abE bE āE | ε.

-Now, let us assume that there exists an occurrence b that does not lie in a pattern ab (the process ensures that any occurrence of the letter a lies into a pattern ab in s ). We choose the leftmost b with this property. Then, s can be written s = αbβ bγ where α ∈ E , β ∈ E and γ ∈ {a, b, ā, b} * . By induction, we have γ ∈ E and s can be generated by the grammar E → E bE bE.

Considering the two cases, s can be generated by the grammar E → E | E bE bE where E → abE bE āE | ε.

It is clear that two shuffles s and s in E lying into the same abequivalence class are necessarily identical. Using a recursive argument, this induces the unicity of a shuffle in E for a given ab-equivalence class. Let E (x) be the generating function for the set E . From the above grammar, we deduce the functional equation

E (x) = 1 + x 2 E (x) 3 and E(x) = E (x) + xE (x) 2 E(x).
The generating function E (x) is known (see A001764 in [START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF]) and given by The positions of the occurrences of the pattern π involve our equivalence relation. So, it would be interesting to study the weaker equivalence relation where two maps are equivalent when they have the same number of occurrences of a given pattern in their shuffle representations.

More generally, such a study could be made on maps using other representations and patterns.

  . The first values for n ≥ 0 are 1, 2, 10, 70, 588, 5544, 56628. For instance, s = aabāā bab bā is a shuffle of the two Dyck words aaāāaā ∈ D a and b bb b ∈ D b . For any shuffle s, we denote by w a (s) (resp. w b (s)) the Dyck word in D a (resp. D b ) obtained from s by deleting the letters b and b (resp. a and ā). In the following, we will extend the definition of w a (s) and w b (s) for any word in {a, ā, b, b} * . In particular, if s is a prefix of shuffles, w a (s) (resp. w b (s)) becomes a prefix of a Dyck word in D a (resp. D b ). For instance, if s = aabāā bab bā, then w a (s) = aaāāaā and w b (s) = b bb b; and if s = aabā, it is a prefix of a shuffle, and w a (s) = aaā and w b (s) = b are prefixes of Dyck words. Then, a word s is a shuffle in S if and only if w a (s) ∈ D a and w b (s) ∈ D b .

  a and w b (s 1 . . . s i-1 ) / ∈ D b then we set s i = b; (iii) If w a (s 1 . . . s i-1 ) ∈ D a and w b (s 1 . . . s i-1 ) ∈ D b then we set s i = b. It is worth noticing that w a (s 1 . . . s i ) (resp. w b (s 1 . . . s i )) is a prefix of a Dyck word in D a (resp. D b ) for all i, 1 ≤ i ≤ 2k. Moreover, at the end of the process we necessarily have w a (s 1 . . . s 2k ) ∈ D a and w b (s 1 . . . s 2k ) ∈ D b which means that s belongs to S.

  [START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF]), and the first values for n ≥ 1 are 1, 2, 5, 12, 31, 81, 216, 583, 1590. Proof. By Lemma 6, it suffices to provide the generating function C(x) of the set C, with respect to the semilength. Every nonempty shuffle s ∈ C is obtained by the grammar C → C | C bC bC where C → aaāC āC | a(C \ε)āC | ε.

  Our study focused on the enumeration of the π-equivalence classes of rooted planar maps in shuffle representation for a pattern π of length at most two not in {a b, āb, bā, ba}. Since the pattern a b does not occur in any shuffle of a rooted planar map, only three patterns are left as open enumeration problems.

Table 1

 1 So, we define the equivalence relation on S and N CS for any pattern π of length at most two:

	: Correspondence between patterns in a shuffle and patterns in a
	rooted planar map.
	relation on Motzkin words.
	The motivation of this paper is to present a similar study for shuffles of
	Dyck words and for rooted planar maps considered with their representation
	by shuffles in N CS.

Table 2 :

 2 Enumeration of π-equivalence classes for shuffles and rooted planar maps.

2

  Theorem 6 The generating function for the set of ab-equivalence classes in S with respect to the semilength is given by This is the sequence A138164 in[START_REF] Sloane | The On-line Encyclopedia of Integer Sequences[END_REF], and the first values for n ≥ 1 are 1, 2, 4, 9, 20, 47, 109, 262, 622. Proof. By Lemma 7, it suffices to provide the generating function E(x) of the set E, with respect to the semilength. Every nonempty shuffle s ∈ E is obtained by the grammar E → E | E bE bE where E → abE bE āE | ε.

	-2 4 sin 1/3 arcsin 3/2 √ 3 sin 1/3 arcsin 3/2 √ 3x	√	3x 2 -3 x	.