Clara Benoit-Pilven

Camille Marchet

Emilie Chautard

Leandro Lima

Marie-Pierre Lambert

Gustavo Sacomoto

Amandine Rey

Audric Cologne

Sophie Terrone

Louis Dulaurier

Jean-Baptiste Claude

Cyril F Bourgeois

Didier Auboeuf
email: didier.auboeuf@inserm.fr

Vincent Lacroix
email: vincent.lacroix@univ-lyon1.fr

Complementarity of assembly-first and mapping-first approaches for alternative splicing annotation and differential analysis from RNAseq data

come

Complementarity of assembly-first and mapping-first approaches for alternative splicing annotation and differential analysis from RNAseq data

In the last 10 years, the prevalence of alternative splicing has been completely re-evaluated. Recent reports claim that more than 90% of multi-exon genes produce at least two splicing variants [START_REF] Pan | Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing[END_REF][START_REF] Wang | Alternative isoform regulation in human tissue transcriptomes[END_REF] . The depth at which transcriptomes can be sampled with next generation sequencing techniques opens the possibility not only to annotate splicing variants in various conditions, but also to detect which transcripts are differentially spliced across pathological and physiological conditions. This growing interest in splicing both as a fundamental process and because of its implication in pathologies [START_REF] Scotti | Rna mis-splicing in disease[END_REF][START_REF] Edery | Association of tals developmental disorder with defect in minor splicing component u4atac snrna[END_REF][START_REF] David | Alternative pre-mrna splicing regulation in cancer: pathways and programs unhinged[END_REF] has been accompanied by an increasing number of methods aiming at analyzing RNAseq datasets [START_REF] Trapnell | Differential gene and transcript expression analysis of rna-seq experiments with tophat and cufflinks[END_REF][START_REF] Wang | Mapsplice: accurate mapping of rna-seq reads for splice junction discovery[END_REF][START_REF] Robertson | De novo assembly and analysis of rna-seq data[END_REF] . The ultimate goal of these methods is to identify and quantify full-length transcripts from short sequencing reads. This task is particularly challenging and recent benchmarks show that all methods still make a lot of mistakes [START_REF] Steijger | Assessment of transcript reconstruction methods for rna-seq[END_REF] . The difficulty of reconstructing full-length transcripts (isoform-centric approaches) also prompted a number of authors to focus on identifying exons that are differentially included within transcripts (exon-centric approaches) [START_REF] Anders | Detecting differential usage of exons from RNA-seq data[END_REF][START_REF] Katz | Analysis and design of rna sequencing experiments for identifying isoform regulation[END_REF][START_REF] Shen | MATS: a Bayesian framework for flexible detection of differential alternative splicing from RNA-Seq data[END_REF][START_REF] Sacomoto | KISSPLICE: de-novo calling alternative splicing events from RNA-seq data[END_REF] .

Whether they are exon-or isoform-centric, methods to study splicing from RNAseq data can further be divided in two main categories [START_REF] Martin | Next-generation transcriptome assembly[END_REF] . The mapping-first approaches first map the reads to the reference genome and the mapped reads are then assembled into exons and eventually transcripts. In contrast, assembly-first approaches first assemble the reads based on their overlaps. The assembled sequences (corresponding to sets of exons) are then aligned to the reference genome.

Mapping-first approaches have been the most used so far, essentially because they were the first to be developed and because they initially required less computational resources. De novo assembly methods were also thought to be restricted to non-model species, where no (good) reference genome is available, and they seemed to be inadequate when an annotated reference genome is available.

Recent progress in de novo transcriptome assembly is clearly changing this view, and the argument of the heavier computational burden does not hold anymore.

The application of de novo assembly to human RNAseq datasets however still remains rare, although some studies have already shown its potential to detect novel biologically relevant splicing variants [START_REF] Dargahi | A pan-cancer analysis of alternative splicing events reveals novel tumor-associated splice variants of matriptase[END_REF][START_REF] Freyermuth | Splicing misregulation of scn5a contributes to cardiac-conduction delay and heart arrhythmia in myotonic dystrophy[END_REF] .

The generalization of de novo assembly approaches for studying splicing in human seems to be mostly impeded by the lack of a clear evaluation of its potential interest in comparison to more traditional mapping-based approaches.

This is the gap we aim at filling with the work presented here.

To achieve this goal, we performed a systematic evaluation of an assembly-first and a mapping-first approach on two RNAseq datasets.

As a first step, we compared pipelines that we developed in parallel, namely KisSplice and FaRLine, because we could easily control their parameters. Any difference between the predictions that is solely due to a parameter setting could be fixed easily, which enabled us to obtain a precise understanding of the irreducible differences between the two approaches.

In a second step, we confirmed the generality of our findings by benchmarking our methods against Cufflinks [START_REF] Trapnell | Differential gene and transcript expression analysis of rna-seq experiments with tophat and cufflinks[END_REF] , MISO [START_REF] Katz | Analysis and design of rna sequencing experiments for identifying isoform regulation[END_REF] and Trinity [START_REF] Grabherr | Trinity: reconstructing a full-length transcriptome without a genome from rna-seq data[END_REF] , which are widely used pipelines.

A significant part of our work has been to manually dissect a number of cases found by only one of the two methods. This enabled us to go beyond a simple qualitative description and provide the community with a precise understanding of which cases are overlooked by each type of method, and where new methods are needed.

All the software and step-by-step protocols presented in this work are freely available at http://kissplice.prabi. fr/pipeline_ks_farline. This should facilitate the reproducibility of our work, and applications to other datasets.

From a general point of view, the combination of approaches we propose should enable to improve splicing-related transcriptomic analyses in physiological and pathological situations.

Results

KisSplice and FaRLine. Figure 1 presents schematically the two pipelines that we developed and compared. A detailed description of each step is given in the Methods section. In the assembly-first approach, a De Bruijn graph is built from the reads. Alternative splicing events, which correspond to bubbles in this graph are enumerated and quantified by KisSplice. Each path is then mapped on the reference genome using STAR and the event is annotated by KisSplice2RefGenome, using the EnsEMBL r75 annotations as an evidence. Importantly, exons not present in the annotations can be identified by this approach. In the mapping-first approach, reads are aligned to the reference genome using TopHat2. Mapped reads are then analyzed by FaRLine, using the EnsEMBL r75 annotations as a guide.

We also tested STAR instead of TopHat2 for the mapping-first pipeline, and found that our main results were essentially unchanged (see Methods).

Quantification of splicing variation is performed similarly in the two pipelines. Only junction reads are considered. Exonic reads are not considered, for reasons exposed in Methods. For the inclusion isoform, there are two junctions to consider. We calculate the mean of the counts of these two junctions.

The differential analysis is performed by a common method for the two approaches: kissDE, which tests if the relative abundance of the inclusion isoform has changed significantly across conditions.

Overall, we developed and adapted jointly these two pipelines in order to minimize the discrepancies that could complicate the comparison.

The majority of frequent isoforms are identified by both approaches. Applying KisSplice and

FaRLine to the same RNAseq datasets generated by the ENCODE consortium (SK-N-SH cell lines treated or not with retinoic acid), we noticed that 68% of the alternatively skipped exons (ASE) identified by KisSplice were also identified by FaRLine and that 24% of ASEs identified by FaRLine were also identified by KisSplice (Fig. 2A). This observation highlights that the mapping-first approach predicts a much larger number of events. This difference in sensitivity is due to the fact that while mapping-first approaches require that each exon junction is covered by at least one read, assembly-first approaches require overlapping reads across the entire skipped exon. Therefore, it can be anticipated that low abundant isoforms, that are covered by few reads, will be reported by mapping, but not by the assembly-first approach. Supporting this prediction, we observed that for ASEs reported only by FaRLine, the number of reads supporting the minor isoform is much lower than in the other categories (Fig. 2 B). The same results were obtained using another RNAseq dataset representing MCF-7 cells expressing or not the DDX5 and DDX17 splicing factors (Supplementary Figure S1).

Having clarified that rare variants are better handled by the mapping-first approach, we decided to filter them out, in order to analyse other differences between the two approaches. Experimental validations by RT-PCR that we performed on rare variants stratified by read support enabled us to clarify that both an absolute and a relative cutoff on the number of reads are required to discriminate variants which can be validated from those which cannot. Indeed, out of the 48 tested cases, we were able to validate 41 (Supplementary Figure S9). The non validated cases indeed corresponded to cases supported by fewer reads. However, what really departed them from the validated cases was their lower relative abundance (Supplementary Figure S10, Supplementary Table 1). In the remaining of our work, we chose to use both criteria and we filtered variants supported by less than 5 reads, and less than 10% compared to the major isoform.

As expected, the proportion of candidates reported simultaneously by both methods increased significantly. Approximately 70% of predicted skipped exons were indeed found by both approaches after filtering lowly expressed isoforms. (Fig. 2C, Supplementary Figure S1C).

Furthermore, the estimation of their inclusion rates was consistent across the two approaches (R 2 > 0.9)).

Beyond the overall concordance of the two approaches in detecting common splicing events, a number of candidates remained reported by only one approach. Since many of them have a highly-expressed minor isoform (supported by more than 100 reads) (Fig. 2D, Supplementary S1D), the failure of one approach to detect them is likely not due to a lack of coverage.

For events only found by one approach, we patiently dissected the reasons why they could have been missed out by the other approach. This enabled us to define 4 main categories which cover 70% of the cases (Fig. 3A) The remaining 30% of cases did not fit into clearly defined biological categories. We however classified them using methodological criteria. The full list of categories is presented in Supplementary Table 2. For each of the 4 main categories, we selected cases to validate experimentally. All 34 RT-PCR validations were successful and are presented in Supplementary Figure S11 confirming that these events are not false positives. Some isoforms are systematically missed by one approach. The first category corresponds to cases that were missed out by the mapping-first approach and corresponds to alternative splicing events involving novel exons or novel combinations of existing exons. The two pipelines compared in this study: KisSplice and FaRLine. The first step of KisSplice is to assemble the reads and extract the splicing events. These events are then mapped back to the reference genome and classified by event type. The annotated and quantified events are then used for the differential analysis between the biological conditions. In contrast, the first step of FaRLine is to map the reads on the reference genome. From this mapping, annotated and quantified events are extracted. Finally, the differential analysis is done with the same method as in the KisSplice pipeline.

There are two reasons to explain why the mapping-first approach does not detect these events. First the mapper may fail to map the reads, or map them to an incorrect location, as junction discovery using short reads is a challenging task. Second, even in the case where the mapper succeeds, FaRLine may fail to report the event because it relies on annotations. Among these 1864 cases, we distinguished 3 sub-categories of errors due to the annotation. Either the exon is unannotated (30%), one of its flanking exon is unannotated (13%) or both exons are annotated but no transcript combining them was annotated (57%).

The assembly-first approach, KisSplice, does not consider annotations, and an interesting resulting advantage is that novel junctions have the same chance to be assembled as known junctions. Mapping assembled novel junctions to the genome is indeed less challenging than read mapping because the assembled sequences are longer.

More importantly, the ability of KisSplice to identify novel splicing events comes from the fact that it introduces known annotations as late as possible in its pipeline (see Methods). Annotations are used as an evidence, not as a filter. AS events involving novel splice sites are clearly identified as such, and can be specifically tested and experimentally validated. More than 99% of the novel splice sites were canonical splice sites (GT-AG).

As an example, the HIRA gene contains a novel exon, whose inclusion is supported by at least 20 reads on each junction (Fig. 3B, Supplementary Figure S8A). This case was overseen by the mapping-first approach, FaRLine. The panel B of the Supplementary Figure S8 shows an example of an ASE not reported by FaRLine because the included exon was not present in the transcripts.

The second category of splicing events identified by only one approach corresponds to recent gene duplications. Untangling the relation between alternative splicing and gene duplication is a difficult topic, subject to debate [START_REF] Kopelman | Alternative splicing and gene duplication are inversely correlated evolutionary mechanisms[END_REF][START_REF] Roux | Age-dependent gain of alternative splice forms and biased duplication explain the relation between splicing and duplication[END_REF] . It is indeed difficult to assess the amount of alternative splicing that occurs within paralogous genes. With the mapping-first approach, the reads stemming from recent paralogs are classified as multi-mapping reads. FaRLine, like the vast majority of mapping-first pipelines, discards these reads for further analysis, as their precise location cannot be clearly established. This results in silently underestimating alternative splicing in recent paralog genes. Note that setting the mapper to keep multi-mapping reads in the analysis leads to overestimating alternative splicing, as all members of the family will be predicted as alternatively spliced. In opposition, de novo assembly can faithfully state that a family of recent paralogs collectively produce two isoforms that vary in their sequence. However, whether the two isoforms are produced from the same locus or from different loci remains undetermined. KisSplice detects these cases of putative AS in paralog genes. Figure 3C illustrates the case with The exon in intron 8 of the HIRA gene is an example of an exon not annotated in EnsEMBL r75. This event was identified by KisSplice but not by FaRLine. (C) RASA4 and RASA4B are 2 paralog genes. KisSplice detected 2 isoforms that could be produced by these 2 genes. FaRLine did not detect any event in either of these genes. The exon skipped is exon 18 in RASA4 (corresponding to exon 17 in RASA4B). The third band on the RT-PCR is the inclusion of another exon in the intron 18 of RASA4. (D) Exon 2 of the RAB5C gene is an example of exon skipping overlapping an Alu element identified only by FaRLine. The events in panel B to C were validated by RT-PCR. (E) The RPAIN gene contains a complex event with a lowly expressed isoform. This weakly expressed isoform was not identified by KisSplice, while the other isoforms were identified by both approaches.

genes RASA4 and RASA4B. Exon 18 in RASA4 (denoted as exon 17 in RASA4B) was detected to be skipped. The exclusion isoform is supported by 160 reads, while the inclusion isoform is supported by 400 reads. The mapping-first approach did not detect either of these isoforms at all. Another example from this category is presented in Supplementary Figure S2C.

The third category of splicing events identified by only one approach corresponds to cases that are missed out by the assembly-first approach. Out of the 1663 cases belonging to this category, a large fraction (21%) corresponds to cases where the skipped exon overlaps a repeat, notably Alu elements. Alu are transposable elements present in a very large number of copies in the human genome [START_REF] Batzer | Alu repeats and human genomic diversity[END_REF] . Most of these copies are located in introns and a number of them have been exonised [START_REF] Lev-Maor | The birth of an alternatively spliced exon: 3'splice-site selection in alu exons[END_REF][START_REF] Sorek | Minimal conditions for exonization of intronic sequences: 5′ splice site formation in alu exons[END_REF] . The reason why the mapping-first approach is able to identify these cases is because even though the reads partially map to repeated sequences, the boundaries of these exons are unique and annotated. Hence the mapper, if set properly, can map these reads to unique annotated exon junctions and is not confused by multiple mappings. Importantly, if the annotations are not provided to the mapper, it will be confused by multiple mappings and will not be able to map the read to the correct location (Supplementary Figure S7). Figure 3D and Supplementary Figure S2D represent two RT-PCR validated Alu-derived exons identified by the mapping-first approach. The assembly-based approach fails to detect most of these events. The reason is that, although they do form bubbles in the DBG generated by the reads, these bubbles are highly branching (supplementary figure http://kissplice.prabi.fr/sknsh/graph_RAB5C_distance_3.html23). Enumerating branching bubbles is computationally very challenging, and may take a prohibitive amount of time. In practice, we restrict our search to the enumeration of bubbles with at most 5 branches (Supplementary Figure S12A).

The fourth category of splicing events identified by only one approach corresponds to cases where more than two splicing isoforms locally coexist, and one of them is poorly expressed compared to the others. The RPAIN gene is a good illustration of such cases (Fig. 3E), as exons 5 and 6 of RPAIN may be skipped and the intron between exons 4 and 5 may be retained. While both methods successfully reported the skipping of exon 6, with exons 5 and 7 as flanking, FaRLine additionally reported the skipping of the same exon, but with exons 4 and 7 as flanking exons. The reason why KisSplice did not report this case is because the junction between exons 4 and 6 is relatively weakly supported. More specifically, this junction is supported by only 55 reads, which accounts for less than 2% of the total number of reads branching out from exon 4. Transcriptome assemblers, like KisSplice, usually interpret such relatively weakly supported junctions as sequencing errors or spurious junctions in highly-expressed genes, therefore disregarding them in the assembly phase (see Supplementary Methods). Supplementary Figure S2E shows another example of a complex event not correctly handled by KisSplice because there were locally more than 5 branches.

Comparison of the approaches after differential analysis.

Beyond the tasks of identifying exon skipping events, a natural question which arises when two conditions are compared is to assess if the exon inclusion rate significantly changed across conditions.

In order to test this, we took advantage of the availability of replicates for both the SK-N-SH cell line and the same cell line treated with retinoic acid. For each detected event, we tested with kissDE [START_REF] Lopez-Maestre | SNP calling from RNA-seq data without a reference genome: identification, quantification, differential analysis and impact on the protein sequence[END_REF] , whether we could detect a significant association between one isoform and one condition. Focusing on those condition-specific events, we again partitioned them in events reported by both methods, by FaRLine only and by KisSplice only. As shown in Fig. 4, the majority of condition-specific events were detected by both approaches. This is the case for instance of exon 22 of gene ADD3 which is clearly more included upon retinoic acid treatment (Fig. 4C), with a DeltaPSI of 27%. The estimation of the DeltaPSI is overall very similar across the two approaches (Fig. 4B) with a correlation of 0.94. The outliers essentially correspond to ASE with several alternative donor/acceptor sites. KisSplice considers these events as different exons while FaRLine considers them as an unique exon, and sums up all the incoming (resp. outgoing) junction counts. Hence, the read counts will differ. Supplementary Figure S8D gives an example.

When focusing on condition-specific events, the proportion of events predicted by only one method increased, for two main reasons. First, some ASE annotated by both approaches were predicted to be differentially included only by one method. This is again due to differences in the quantification of the inclusion rate, especially for ASE with multiple 5′ and 3′ splice sites. Second, some of the exons that were missed out by one method at the identification step happened to be condition specific. This is the case of an exon in NINL intron 5 (Fig. 4D), only identified by KisSplice because it was not annotated. This is also the case of SAR1B exon 3 (Fig. 4E), only identified by FaRLine because it overlaps with an Alu element. The analysis of the MCF-7 RNAseq dataset gave very similar results (Supplementary Figure S3).

The observation that many of the AS events that were annotated only by one method are differentially regulated across conditions confirms that these AS events should not be discarded from the analysis. Focusing only on AS events annotated by one approach may lead to miss splicing events which are central in the biological context.

Overlap with other methods. In a first step, we picked FaRLine and KisSplice as examples of a mapping-first and an assembly-first approach respectively. Clearly, there are other published methods in both categories. MISO is probably the most widely used to annotate AS events. We therefore ran it on the same datasets to check how its predictions overlapped with ours. As shown in Fig. 5A (SK-N-SH dataset), 77% of predictions made by MISO were common to both FaRLine and KisSplice, 18% were only common with FaRLine, 2% were only common to KisSplice and the remaining 3% were specific to MISO. The overlap between the different methods was very similar when the MCF-7 RNAseq dataset was used (Supplementary Figure S4A). Overall, almost all candidates predicted by MISO were also predicted by FaRLine. This large overlap with FaRLine was expected, because both are mapping-first approaches. This also shows that the differences between mapping-and assembly-first approaches reported above are not limited to one mapping-first approach.

Besides exon-centric approaches, which aim at finding the differentially spliced exons, there is also a number of published methods which are isoform-centric and have the more ambitious goal to reconstruct full-length transcripts at the expense of underestimating alternative splicing.

The most widely used mapping-first and isoform-centric approach is Cufflinks 6 that we compared to FaRLine using the same dataset. As shown in Fig. 5B (and Supplementary Figure S4B), we found that the vast majority of ASE were predicted by both approaches.

Finally, we compared KisSplice to one of the most widely used de-novo transcriptome assembler, Trinity [START_REF] Grabherr | Trinity: reconstructing a full-length transcriptome without a genome from rna-seq data[END_REF] . As shown in Fig. 5D (and Supplementary Figure S4D), most ASE found by Trinity were also found by KisSplice. However, KisSplice was significantly more sensitive. The goal of Trinity is to assemble the major isoforms for each gene, it therefore largely under-estimates alternative splicing, especially inclusion/exclusion of short sequences.

For completeness sake, we also provide an all-vs-all comparison (Supplementary Figure S5). An interactive version of this Figure is available at http://kissplice.prabi.fr/pipeline_ks_farline/. The list of events found by any used method can be retrieved from this interactive figure and analysed in IGV, to reproduce the sashimi plots of the paper. The general conclusions from these comparisons is that there is a clear distinction between mapping-first and assembly-first approaches, and between exon-centric and isoform-centric approaches, the latter being less sensitive.

Discussion

De novo assembly is usually applied to non-model species where no (good) reference genome is available. We show here that even when an annotated reference genome is available, using assembly offers a number of advantages. We named this approach "assembly-first" because it does use a reference genome, but as late as possible in the process, in order to minimize the a priori on which exons should be identified.

Using this strategy, we identified novel alternatively skipped exons, which were not identified by traditional read mapping approaches (Fig. 3 and Supplementary Figure S2). While it is believed that the human genome is fully annotated, it is important to underline that we have not yet established a final map of the parts of the genome that can be expressed. It can be anticipated that sequencing of single-cells from different parts of the body will lead to the discovery of a huge diversity of transcripts and that a substantial number of new exons will be discovered. An example is the case of unannotated skipped exons which overlap with repeat elements. We cannot exclude that this category is currently largely under-annotated.

We also showed that assembly-first approach has the ability to detect splicing variants within recently duplicated genes (Fig. 3 and Supplementary Figure S2). This is because mapping approaches discard reads which map to multiple genomic locations. Identification of such splicing variants produced from different genomic regions sharing sequence similarities (e.g. paralog genes, pseudogenes) is however very important, since splicing variants generated from paralogous genes but also from pseudogenes may have different biological functions [START_REF] Poursani | Differential expression of oct4 pseudogenes in pluripotent and tumor cell lines[END_REF] .

Conversely, we showed that some ASE were detected only by the mapping-first approach. As shown in Fig. 2 (and Supplementary Figure S1), we observed that the mapping-first approach has a better ability to detect lowly-expressed splicing variants. Although such lowly-expressed splicing variants are often considered as "noise" or biologically non relevant, caution must be taken with such assumptions for several reasons. First, mRNA expression level is not necessarily correlated with protein expression level. Second, as observed from single-cell transcriptome analyses, some mRNAs can be expressed in few cells, within a cell population (e.g. they are expressed at a specific cell cycle step) and may therefore appear to be expressed at a low level in total RNAs extracted from a mixed cell population [START_REF] Bacher | Design and computational analysis of single-cell rna-sequencing experiments[END_REF] . Therefore, computational analysis should not systematically discard lowly-expressed splicing variants and filtering these events should depend on the biological questions to be addressed.

We also observed that the mapping-first approach better detects exons corresponding to annotated-repeat elements (Fig. 3 and Supplementary Figure S2). While it has been assumed for a long time that repeat elements are "junk", increasing evidences support important biological functions for such elements. For example, repeat elements like Alu can evolve as exons and the presence of Alu exons in transcripts has been shown to play important regulatory functions [START_REF] Sorek | Minimal conditions for exonization of intronic sequences: 5′ splice site formation in alu exons[END_REF][START_REF] Shen | Widespread establishment and regulatory impact of alu exons in human genes[END_REF] .

When two methods give non-overlapping predictions, the temptation could be to focus on exons found by both approaches and to discard the others. We argue that this is not the best option, because approach-specific cases can be validated experimentally, and also because many of them correspond to regulated events, i.e. the inclusion isoform is significantly up or down regulated depending on the experimental condition.

In conclusion, combining mapping-and assembly-first approaches allows to detect a larger diversity of splicing variants. This is very important towards the in depth characterization of cellular transcriptome although other approaches are further required to analyze their biological functions.

From a computational perspective, a number of challenges are still ahead. The co-development of two approaches enabled us to narrow down the list of difficult instances not properly dealt with by at least one approach, but we cannot exclude that some categories are still missed out by both approaches. The categories of challenging cases that we defined in Fig. 3: lowly-expressed variants, exonised Alu, complex splicing variants, paralogs have been overlooked up to now. Possibly because they are much harder to detect, they have been assumed to play a minor role in transcriptomes, but more recent studies however argues the opposite.

For exonised ALUs, paralog genes and genes with complex splicing patterns, the possibility to sequence longer reads with third generation techniques [START_REF] Tilgner | Defining a personal, allele-specific, and single-molecule long-read transcriptome[END_REF][START_REF] Bolisetty | Determining exon connectivity in complex mRNAs by nanopore sequencing[END_REF] should prove very helpful. The number of reads obtained with these techniques is however currently much lower than with Illumina, thereby preventing their widespread use for differential splicing, for which the sequencing depth, and not so much the length of the reads, is the critical parameter which conditions the statistical power of the tests. In the coming years, methods combining second and third generation sequencing should enable to obtain significant advances in RNA splicing.

Material and Methods

FaRLine and KisSplice. Figure 1 shows the two pipelines that we are comparing. While STAR and TopHat are third-party softwares, we developed the other methods ourselves. KisSplice was first introduced in Sacomoto et al. [START_REF] Sacomoto | KISSPLICE: de-novo calling alternative splicing events from RNA-seq data[END_REF] . The novelty here is that its usage is now possible in the case where a reference genome is available, which required specific methodological developments implemented in the newly released KisSplice2RefGenome software. kissDE was first introduced in Lopez-Maestre et al. [START_REF] Lopez-Maestre | SNP calling from RNA-seq data without a reference genome: identification, quantification, differential analysis and impact on the protein sequence[END_REF] in the context of SNPs for non-model species. We present here its extension for alternative splicing. FaRLine is a new mapping-first pipeline, that we introduce in this paper. It is the RNAseq pipeline associated to the FasterDB database [START_REF] Mallinjoud | Endothelial, epithelial, and fibroblast cells exhibit specific splicing programs independently of their tissue of origin[END_REF] and was already successfully applied to the analysis of the effect of metformin treatment on myotonic dystrophy type I (DM1) with a validation rate of 95% [START_REF] Laustriat | In Vitro and In Vivo Modulation of Alternative Splicing by the Biguanide Metformin[END_REF] . Specifically, 20 cases of ASE regulated by the metformin treatment were tested, and 19 were validated. In this paper, we provide additional validations of FaRLine with similar validation rates (36 out of 38), Supplementary Figure S19.

For the sake of self-containment, we explain all methods here.

KisSplice. KisSplice is a local transcriptome assembler. As most short reads transcriptome assemblers [START_REF] Robertson | De novo assembly and analysis of rna-seq data[END_REF][START_REF] Grabherr | Trinity: reconstructing a full-length transcriptome without a genome from rna-seq data[END_REF][START_REF] Schulz | Oases: robust de novo rna-seq assembly across the dynamic range of expression levels[END_REF] , it relies on a De Bruijn graph (DBG). Its originality lies in the fact that it does not try to assemble full-length transcripts. Instead, it assembles the parts of the transcripts where there is a variation in the exon content. By aiming at a simpler goal, it can afford to be more exhaustive and identify more splicing events. The key concept on which KisSplice is built is that variations in the nucleotide content of the transcripts will correspond to specific patterns in the DBG called bubbles (Supplementary Figure S13). KisSplice's main algorithmic step therefore consists in enumerating all the bubbles in the graph built from the reads. Examples of bubbles in the DBG and explanation of the parameters used to filter out sequencing errors and repeat-induced bubbles are given in Supplementary Methods.

Annotating the events with KisSplice2RefGenome. KisSplice outputs bubbles in the form of a pair of fasta sequences. Clearly, such information is insufficient to analyse alternative splicing for model species. KisSplice2RefGenome enables to provide for each bubble: the gene name, the AS event type, the genomic coordinates and the list of splice sites used (novel or annotated). Bubbles found by KisSplice are mapped to the reference genome using STAR, with its default settings, which means that in the case of multi-mappings, STAR reports all equally best matches. The mapping results are then analysed by KisSplice2RefGenome. Bubbles are classified in sub-types depending on the number of blocks obtained when mapping each path of the bubble to the genome (Supplementary Figure S14). For exon skipping, the longer path of the bubble corresponds to 3 blocks, while the lower path corresponds to 2 blocks. The splice sites are located and compared to the annotations. Events with novel splice sites are reported explicitly as such in the output of the program.

In the case where the bubble corresponds to a genomic insertion or deletion, it exhibits a specific pattern in terms of block numbers (one block for one path and two blocks for the other) and is reported separately.

The criterion of the number of blocks is discriminative in most cases. However, there is a possible confusion between intron retentions and genomic deletions, since in both cases, the longer path will map into one block and the lower path in two blocks. In this case, we also use the distance between the blocks, and introduce a user-defined threshold, which we set to 50nt, below which the bubble is classified as a genomic deletion, and above which it is classified as an intron retention.

In the special case where the exon flanking the AS event is very short (less than k nt), the number of blocks is increased for both paths, but the difference of number of blocks remains unchanged.

In the special case where there is a genomic polymorphism located less than k nt apart from the AS event, KisSplice will report several bubbles (possibly all combinations of genomic and transcriptomic variants). This redundancy is removed in KisSplice2RefGenome where the primary focus in on splicing.

In the case where the bubble maps to two locations on the genome, a distinction is made between the case of exact repeats where both paths map to both locations and inexact repeats where each path maps to a distinct location (Supplementary Figure S12B). The cases of exact repeats correspond to recent gene duplications.

FaRLine. FasterDB EnsEMBL r75 annotation. FasterDB RNAseq Pipeline, FaRLine, uses the FasterDB-based EnsEMBL r75 annotation database. FasterDB is a database containing all annotated human splicing variants [START_REF] Mallinjoud | Endothelial, epithelial, and fibroblast cells exhibit specific splicing programs independently of their tissue of origin[END_REF] .

Each transcripts present in the FasterDB, is composed of a succession of exons, that we call transcript exons (represented in blue in Supplementary Figure S15). The genomic exons (represented in red in Supplementary Figure S15) are defined by projecting the transcript exons. First, the transcript exons are grouped by position. Then each group of exons defines a projected exon with the following rules:

• The start is the leftmost start of the non-first-exon of the group.

• The end is the rightmost end of the non-last-exon of the group that ends before the start of the next group of exons.

When the most frequent event annotated in the transcripts is an intron retention, the projected genomic exon is defined as a combination of the two exons flanking the retained intron. In Supplementary Figure S15, the exons 5 and 6 and the intron 5 are considered as one unique exon. As events included within one exon are not tested, this results in some events being missed.

Mapping. The first step of FaRLine is to map the reads to a reference genome. This step is done using Tophat-2.0.11 [START_REF] Trapnell | Differential gene and transcript expression analysis of rna-seq experiments with tophat and cufflinks[END_REF] . tophat-min-intron-length 30-max-intron-length 1200000\-p 8 [-solexa1.3-quals for Sknsh_rep1 and Sknsh_rep2]\-transcriptome-index A transcriptome index has been built by TopHat using EnsEMBL r75 annotations in gtf format. When a transcriptome index is used, the mapping steps are modified: instead of aligning first to the genome, which is the default behavior, TopHat uses Bowtie to align the reads to the transcript sequences first, then align the remaining unmapped reads to the genome. Minimal and maximal intron lengths have been modified (default 70 and 500000) to maximize the number of junctions detected, according to the statistics provided by FasterDB EnsEMBL r75 annotations.

The resulting alignment files have been filtered using samtools 0.1.19 [START_REF] Li | The sequence alignment/map format and samtools[END_REF] . Samtools view -F 260 -f 1 -q 10 -b With this step, only the primary alignments are kept. The minimum read alignment quality was set up so that multi-mapping reads were removed from the alignment file.

Annotation and quantification of alternative splicing events. For each gene, all the reads with at least one base overlapping the gene from the start to the end coordinates are retrieved. CIGAR strings are then used to find the alignments blocks. Junction reads are identified by the presence of at least one'N' letter in the CIGAR. Junction reads were filtered if:

• More than 10% of soft-clipping was detected in the alignment (it should not be the case with TopHat).

• An indel was close to the junction site, as it would make the junction position uncertain.

Junction read alignments are then processed block by block sequentially from left to right. Alignment blocks under 4 bp on read extremities are removed from the reads as we considered it is not sufficient to identify correctly the mapping localization. Then each block is compared to FasterDB annotations to check if the block boundaries correspond to known exons annotated in FasterDB, or to a putative new acceptor or donor site. First and last alignment blocks for each read must overlap one and only one exon for a read to be considered. For the inner blocks, if alignment blocks map to a succession of exons and introns, it is considered as an intron retention. For the acceptors and donors, we also added a supplementary filter. If a new donor is identified within a junction, we check if the junction also has an acceptor identified of the same length +/-1bp on the other side of the junction, showing most probably a problem of mapping. Once all the blocks are identified, the block annotations are used to annotate putative alternative splicing events: alternative skipped exon, multiple exon skipping, acceptor, or donor sites.

Once all the junction reads are processed, the alternative splicing events identified are pooled and the reads participating to each event are quantified, as well as the known exon-exon junction. If an exon-exon junction is annotated with multiple known acceptors and/or donors, all the possible junction reads are quantified and summed up. To fasten the quantification step, a junction coordinate file with the corresponding read numbers is produced from the read alignment using the same filters than described above and will be used for all the quantification tools: junction, exon skipping, acceptor and donor.

A challenge in defining the alternative skipped exon events is to identify the flanking exons. In the first version of FaRLine, these flankings exons were defined as the closest annotated genomic exons. This rule led to miss a lot of ASE events. Therefore, to define the flanking exons, we now use the information contained in the transcripts and in the reads. We consider each junction which skips an exon and is covered by at least one read. If this junction is annotated in the transcripts, we extract all annotated events containing this junction. Else, we annotate the event with the longest covered inclusion isoform. It allows FaRLine to be more robust to the incompleteness of the annotation compared to other methods, like MISO (Supplementary Figure S6). Panel C of Supplementary Figure S8 gives an example of an ASE reported by FaRLine but not by MISO because the exclusion isoform is not annotated in the transcripts.

Comparison with STAR. We also mapped the reads with STAR, ran FaRLine on these alignments and compared the predicted skipped exons with KisSplice. The main results are similar to what we found with TopHat. Indeed, without any filter, 69% of ASE annotated by KisSplice are also found by FaRLine and 24% of FaRLine's event by KisSplice (compared to 68% and 24% respectively for the mapping with TopHat). When we filter out the events with an unfrequent variant, we show that approximately 70% of predicted ASE are found by both approaches.

Quantification and differential analysis. Both pipelines perform ASE detection and quantification. The quantification step was done similarly in the two pipelines where only the junction reads were taken into account. To evaluate if using exonic reads in the quantification could increase the accuracy of our methods, we ran KisSplice on the MCF-7 dataset with the option -exonic reads set to on. In doing so, only the inclusion rate of the AS events changes. When comparing usage of only junction reads to usage of both junction and exonic reads, we observed that the p-values calculated strongly correlate as shown in Supplementary Figure S16. We found that some AS events became significant upon the addition of exonic reads but the opposite also happened. Inspection of these events revealed that many are borderline cases, where the p-value is close, but slightly above 5%. A manual inspection of the AS events with a very different p-value upon addition of exonic reads revealed that they correspond to exons overlapping alternative first or last exons (see STARD4, Supplementary Figure S17A) or novel exons located in poorly spliced introns (see PANK2 and PRRC2B, Supplementary Figure S17 B andC). Overall, we concluded that exonic reads can bring some statistical power in cases where the skipped exon does not overlap with any other event. In case of more complex events, exonic reads tend to "pollute" the pairwise comparison.

The last step of the pipelines is the differential analysis of the expression levels of the variants. This task is performed using the kissDE [START_REF] Lopez-Maestre | SNP calling from RNA-seq data without a reference genome: identification, quantification, differential analysis and impact on the protein sequence[END_REF] R package, which takes as input a table of read counts as in Supplementary Figure S18, and outputs a p-value and a DeltaPSI (Percent Spliced In).

Our statistical analysis adopted the framework of count regression with Negative Binomial distribution. We considered a 2-way design with interaction, with isoforms and experimental conditions as main effects. Following the Generalized Linear Model framework, the expected intensity of the signal was denoted by λ ijk and was decomposed as:

λ μ α β αβ = + + + log ( ) (1)
ijk i j i j
where μ is the local mean expression of the gene, α i the contribution of splicing variant i on the expression, β j the contribution of condition j to the total expression, and (αβ) ij the interaction term. The target hypothesis was αβ =

{ }

H : ( ) 0 ij 0 i.e. no interaction between the variant and the condition. If this interaction term is not null, a differential usage of a variant across conditions occurred. The test was performed using a Likelihood Ratio Test with one degree of freedom. To account for multiple testing, p-values were adjusted with a 5% false discovery rate (FDR) following a Benjamini-Hochberg procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] .

In addition to adjusted p-values, we report a measure of the magnitude of the effect. The measure we provide is based on the Percent Spliced In (PSI): If counts for a variant are below a threshold, then the PSI is not calculated. This prevents from over-interpreting large magnitudes derived from low counts. When several replicates are available for a condition, then a PSI is computed for each replicate, and we calculate their mean.

Finally, we output the DeltaPSI: unless one of the mean PSI of a condition could not be estimated. The higher the DeltaPSI, the stronger the effect.

In practice, we consider only DeltaPSI larger than 0.1, a threshold below which it is difficult to perform any experimental validation.

SK-N-SH dataset. We downloaded a total of 959 M reads from http://genome.crg.es/encode_RNA_dashboard/hg19/35. They correspond to long polyA+ RNAs generated by the Gingeras lab, and are also accessible with the following accession numbers (ENCSR000CPN -SRA: SRR315315, SRR315316 and ENCSR000CTT -SRA: SRR534309, SRR534310). For cell lines treated by retinoic acid, the reads were 76nt long, while they were 100nt long for the non treated cells. Hence we trimmed all reads to 76nt.

MCF-7 dataset. MCF-7 were transfected (two biological replicates) with siRNA targeting both DDX5 and DDX17 RNA helicases, and total RNA were extracted as described previously [START_REF] Dardenne | RNA Helicases DDX5 and DDX17 Dynamically Orchestrate Transcription, miRNA, and Splicing Programs in Cell Differentiation[END_REF] . cDNA synthesis was made using the TruSeq Stranded Total RNA protocol after Ribo-Zero Gold-mediated elimination of ribosomal RNA (Beckman Coulter Genomics). High throughput sequencing (2 × 125 bp) was carried out on an Illumina HiSeq 2500 platform (Beckman Coulter Genomics), generating between 45 and 50 millions of paired-end pairs of reads. Raw datasets are available on GEO under the accession number GSE94372. Reads were trimmed according to standard quality control filters using prinseq [START_REF] Schmieder | Quality control and preprocessing of metagenomic datasets[END_REF] and adapter were removed using cutadapt [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF] . The resulting reads had length between 25 and 125nt. Because MISO is unable to deal with reads of unequal length, we selected only reads with length larger than 100nt (87% of the reads) and trimmed longer reads to 100nt.

Figure 1 .

 1 Figure1. The two pipelines compared in this study: KisSplice and FaRLine. The first step of KisSplice is to assemble the reads and extract the splicing events. These events are then mapped back to the reference genome and classified by event type. The annotated and quantified events are then used for the differential analysis between the biological conditions. In contrast, the first step of FaRLine is to map the reads on the reference genome. From this mapping, annotated and quantified events are extracted. Finally, the differential analysis is done with the same method as in the KisSplice pipeline.

Figure 2 .

 2 Figure 2. Comparison of the ASE identified by the assembly-first and mapping-first pipelines. (A) Venn diagram of ASEs identified by the two pipelines. FaRLine detected many more events than KisSplice. 68% of ASE found by KisSplice were also found by FaRLine and 24% of ASE detected by FaRLine were also found by KisSplice. (B) Boxplot of the expression of the minor isoform in the 3 categories defined in the Venn diagram of panel A: ASE identified only by FaRLine, ASE identified by both pipelines and ASE identified only by KisSplice. The number of reads supporting the minor isoform of the ASE identified by FaRLine is overall much lower. Many isoforms are supported by less than 5 reads. (C) Venn diagram of ASEs identified by the two pipelines after filtering out the poorly expressed isoforms (less than 5 reads, or less than 10% of the number of reads supporting both isoforms). The common events represent a larger proportion than before filtering: 77% of the ASE identified by FaRLine and 69% of the ASE identified by KisSplice. (D) Boxplot of the expression of the minor isoform in the 3 categories defined in the Venn diagram of panel C: ASE identified only by FaRLine, ASE identified by both pipelines and ASE identified only by KisSplice. The distribution of the number of reads supporting the minor isoform is similar for the 3 categories with highly expressed variants in each category.

Figure 3 .

 3 Figure 3. (A)Main categories explaining why some exons are detected by only one method. (B) The exon in intron 8 of the HIRA gene is an example of an exon not annotated in EnsEMBL r75. This event was identified by KisSplice but not by FaRLine. (C) RASA4 and RASA4B are 2 paralog genes. KisSplice detected 2 isoforms that could be produced by these 2 genes. FaRLine did not detect any event in either of these genes. The exon skipped is exon 18 in RASA4 (corresponding to exon 17 in RASA4B). The third band on the RT-PCR is the inclusion of another exon in the intron 18 of RASA4. (D) Exon 2 of the RAB5C gene is an example of exon skipping overlapping an Alu element identified only by FaRLine. The events in panel B to C were validated by RT-PCR. (E) The RPAIN gene contains a complex event with a lowly expressed isoform. This weakly expressed isoform was not identified by KisSplice, while the other isoforms were identified by both approaches.

Figure 4 .

 4 Figure 4. (A) Condition-specific variants identified by FaRLine, KisSplice or both methods. Within dashed lines are events identified by both approaches but detected as condition-specific by only one approach. (B) DeltaPSI as estimated by KisSplice and FaRLine, for events identified by both methods. The red dots represent complex events for which KisSplice found at least 2 'bubbles' . (C) Exon 22 of the ADD3 gene is an example of regulated ASE identified by both approaches. (D) A new exon in intron 5 of NINL gene is identified only by KisSplice. The inclusion of this exon is differentially regulated between the 2 experimental conditions. (E) Because exon 3 of the SAR1B gene is an exonised Alu element, only FaRLine identified this event. Moreover this exon is significantly more included in the treated cells (SK-N-SH RA) compared to the control cells.

Figure 5 .

 5 Figure 5. (A) 77% of ASE identifed by MISO are also annotated by FaRLine and KisSplice. 18% of MISO's ASE are also annotated by FaRLine while only 2% of MISO's ASE are also annotated by KisSplice. Finally, only 3% of these ASEs are only annotated by MISO. (B) Most of the events annotated by Cufflinks are identified by FaRLine. (C) GTF2I exon 13 is an example of an ASE annotated by FaRLine but not by Cufflinks. Indeed, Cufflinks only identified the isoform corresponding to the exon inclusion. (D) Most of the events annotated by Trinity are also annotated by KisSplice. But half of the ASE annotated by KisSplice are not annotated by the global assembler Trinity. (E) KisSplice annotates an ASE in the RFWD2 gene, while Trinity only identified the isoform corresponding to the exon inclusion. The events in panels C and E have been validated by RT-PCR.

Scientific REPoRTS | (2018) 8:4307 | DOI:10.1038/s41598-018-21770-7

Acknowledgements

This work was performed on the computing facilities of the computing center LBBE/PRABI and the PSMN (Pole Scientifique de Modelisation Numerique) computing center of ENS de Lyon. This work was funded by the ANR-12-BS02-0008 (Colib'read) by the ABS4NGS ANR project (ANR-11-BINF-0001-06), Action n3.6 Plan Cancer 2009-2013, Fondation ARC (Programme Labellisé Fondation ARC 2014, PGA120140200853) and INCa (2014-154). Doctoral fellowships from ARC 1 -Région Rhône-Alpes (C.B.P), Science Without Borders -CNPq -Brazil (L.L. -grant process number 203362/2014-4), ARS Rhône-Alpes (A.R.) and post-doctoral fellowships from Fondation ARC (M.P.L).

Computational requirements, software availability and reproducibility of the results. FaRLine took 45 hours and 10 Go of RAM. The time-limiting step was TopHat2, which took 41 hours, even parallelised on 8 cores. When STAR was tested instead of TopHat2, it took 4 hours, but 30 Go of RAM. KisSplice took 30 hours and 10 Go of RAM. The RAM-limiting step was STAR which took 30Go of RAM. All the steps of the pipelines can be reproduced using the following tutorial: http://kissplice.prabi.fr/pipeline_ks_farline.

Experimental Validation. SK-N-SH cells were purchased from the American Type Culture Collection (ATCC) and cultured using EMEM medium (ATCC) complemented with 10% FBS (Thermo Fisher Scientific). Cells were differentiated for 48 h using 6 μM of all-trans retinoic acid (Sigma-Aldrich). After harvesting, total RNA were extracted using Tripure isolation reagent (Sigma-Aldrich), treated with DNase I (DNAfree, Ambion) for 30 min at 37 °C and reverse-transcribed (RT) using M-MLV reverse transcriptase and random primers (Invitrogen). Before PCR, all RT reaction mixtures were diluted at 2.5 ng μL of initial RNA. PCR reactions were performed using GoTaq polymerase (Promega).

MCF7 cells were cultured as described in [START_REF] Dardenne | RNA Helicases DDX5 and DDX17 Dynamically Orchestrate Transcription, miRNA, and Splicing Programs in Cell Differentiation[END_REF] . RT-PCRs were performed using the same protocol as for SK-N-SH cells.

Author Contributions

Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-21770-7.

Competing Interests:

The authors declare no competing interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


a continuous-time system with uncertain but bounded (by the maximum sampling interval) time-varying delay in the control input. We verify that the LMI sufficient conditions for stability of [START_REF] Fridman | An improved stabilization method for linear timedelay systems[END_REF] are valid also in the case of piecewisecontinuous delay and derive LMIs for the feedback gain. The conditions which we obtain are robust with respect to different samplings with the only requirement that the maximum sampling interval is not greater than h. As a by-product we show that for h → 0 the conditions coincide with the necessary and sufficient conditions for the continuous-time stabilization. Such convergence in H 2 framework and related results were proved in [START_REF] Yu | Asymptotic analysis of digital control systems[END_REF], [START_REF] Chen | h2 optimal sampled-data control[END_REF], [START_REF] Fridman | Use of models with aftereffect in the problem of design of optimal digital control[END_REF], [START_REF] Osborn | An exact treatement of the achievable closed-loop h2-performance of sampled-data controllers: From continuous-time to open-loop[END_REF], [START_REF] Trentelman | Sampled-data and discrete-time h 2 optimal control[END_REF] and [START_REF] Oishi | A bound of conservativeness in sampled-data robust stabilization and its dependence on sampling periods[END_REF].

For the first time the new approach allows to develop different robust control methods for the case of sampled-data control. The LMIs are affine in the system matrices and thus for the systems with polytopic type uncertainty the quadratic stabilization conditions, where the common Lyapunov functional for different vertices of the polytope is used, readily follow. We derive a parameter dependent solution, where different Lypunov functionals are used for different vertices by modifying results of (Fridman & Shaked, 2003). We also consider the regional stabilization by sampled-data saturated state-feedback, where we give an estimate on the domain of attraction. For continuous-time stabilization of state-delayed systems by saturated-feedback see e.g. [START_REF] Dambrine | Feedback control of time-delay systems with bounded control and state[END_REF], [START_REF] Tarbouriech | Synthesis of controllers for continuoustime delay systems with saturating controls via lmi's[END_REF], [START_REF] Cao | Stability analysis of linear time-delay systems subject to input saturation[END_REF] and [START_REF] Fridman | Regional stabilization and control of timedelay systems with saturating actuators[END_REF].

Notation: Throughout the paper the superscript 'T ' stands for matrix transposition, R n denotes the n dimensional Euclidean space with vector norm | • |, R n×m is the set of all n × m real matrices, and the notation P > 0, for P ∈ R n×n means that P is symmetric and positive definite. Given ū = [ū 1 , ..., ūm ] T , 0 < ūi , i = 1, ..., m, for any u = [u 1 , ..., u m ] T we denote by sat(u, ū) the vector with coordinates sign(u i )min(|u i |, ūi ). By stability of the system we understand the asymptotic stability of it.

where x(t) ∈ R n is the state vector, u(t) ∈ R m is the control input.

We are looking for a piecewise-constant control law of the form

u(t) = u d (t k ), t k ≤ t < t k+1 , where u d is a discrete-time control signal and 0 = t 0 < t 1 < • • • < t k < • • • are the sampling instants. Our objective is to find a state-feedback controller given by u(t) = Kx(t k ), t k ≤ t < t k+1 , (3)
which stabilizes the system. We represent a piecewise-constant control law as a continuous-time control with a time-varying piecewise-continuous (continuous from the right) delay τ (t) = t -t k as given in [START_REF] Astrom | Adaptative Control[END_REF]. We will thus look for a state-feedback controller of the form:

u(t) = Kx(t -τ (t)). (4)
Substituting ( 4) into (2), we obtain the following closed-loop system:

ẋ(t) = Ax(t) + BKx(t -τ (t)), τ (t) = t -t k , t k ≤ t < t k+1 . (5)
We assume that A1 t k+1 -t k ≤ h ∀k ≥ 0. From A1 it follows that τ (t) ≤ h since τ (t) ≤ t k+1 -t k . We will further consider (5) as the system with uncertain and bounded delay.

Stability of the Closed-Loop System

Similarly to [START_REF] Fridman | An improved stabilization method for linear timedelay systems[END_REF], where the continuous delay was considered, we obtain for the case of piecewise-continuous delay the following result: Lemma 1. Given a gain matrix K, the system (5) is stable for all the samplings satisfying A1, if there exist n×n matrices 0 < P 1 , P 2 , P 3 , Z 1 , Z 2 , Z 3 and R > 0 that satisfy the following LMIs:

Ψ 1 < 0, and R [0 K T B T ]P * Z ≥ 0, (6)
where

P = P 1 0 P 2 P 3 , Z = Z 1 Z 2 * Z 3 , Ψ 1 = Ψ 0 + hZ + 0 0 0 hR , Ψ 0 = P T 0 I A + BK -I + 0 I A + BK -I T P.
Proof is based on the following descriptor representation of (5) [START_REF] Fridman | New lyapunov-krasovskii functionals for stability of linear retarded and neutral type systems[END_REF]:

ẋ(t) = y(t), 0 = -y(t)+(A+BK)x(t) -BK t t-τ (t) y(s)ds, (7)
which is valid in the case of piecewise-continuous delay τ (t) for t ≥ 0. Given a matrix K and initial condition x(t) = φ(t)(t ∈ [-h, 0]), where φ is a piecewise continuous function, x(t) satisfies (5) for t ≥ 0 iff it satisfies [START_REF] Fridman | Use of models with aftereffect in the problem of design of optimal digital control[END_REF]. Note that the descriptor system [START_REF] Fridman | Use of models with aftereffect in the problem of design of optimal digital control[END_REF] has no impulsive solutions since in (7) y(t) is multiplied by the nonsingular matrix I [START_REF] Fridman | Stability of linear descriptor systems with delay: A lyapunov-based approach[END_REF]. We apply the Lyapunov-Krasovskii functional of the form:

V (t) = V 1 + V 2 , (8)
where

x(t) = col{x(t), y(t)}, E = I n 0 0 0 , P = P 1 0 P 2 P 3 , P 1 = P T 1 > 0, (9a-d) and V 1 = xT (t)EP x(t), V 2 = 0 -h t t+θ y T (s)Ry(s)dsdθ, (9e-f)
which satisfies the following inequalities

a|x(t)| 2 ≤ V (t) ≤ b sup s∈[-h,0] |x(t + s)| 2 , a > 0, b > 0. ( 10
)
Differentiating V (t) along the trajectories of ( 7) for t ≥ h we find (see [START_REF] Fridman | An improved stabilization method for linear timedelay systems[END_REF]) that

V (t) < x (t) T Ψ 1 x (t) < -c|x(t)| 2 , c > 0, (11)
Provided that (6a,b) hold. Integrating [START_REF] Fridman | An improved stabilization method for linear timedelay systems[END_REF] we have

V (t) -V (h) ≤ -c t -h |x(s)|ds (12)
and, hence, [START_REF] Fridman | Regional stabilization and control of timedelay systems with saturating actuators[END_REF] Hale & Lunel, 1993, p168) and thus ẋ, defined by the right-hand side of ( 5), satisfy

yields |x(t)| 2 ≤ V (t)/a ≤ V (h)/a < b/a sup s∈[-h,0] |x(h + s)| 2 . Since sup s∈[-h,0] |x(h + s)| ≤ c 1 sup s∈[-h,0] |φ(s)|,c 1 > 0 (cf.
sup s∈[-h,0] | ẋ(h + s)| ≤ c 2 sup s∈[-h,0] |φ(s)|, c 2 > 0, we obtain that |x(t)| 2 ≤ c 3 sup s∈[-h,0] |φ(s)| 2 , c 3 > 0. (13)
Hence ( 5) is stable (i.e. x(t) is bounded and small for small φ). To prove asymptotic stability we note that x(t) is uniformly continuous on [0, ∞) (since ẋ(t) defined by the right-hand side of (5 is uniformly bounded). Moreover, [START_REF] Fridman | Parameter dependent stability and stabilization of uncertain time-delay systems[END_REF] yields that |x(t)| 2 is integrable on [0, ∞). Then, by Barbalat's lemma, x(t) → 0 for t → ∞ Consider now the continuous state-feedback

u(t) = Kx(t) (14)
and the closed-loop system (2), ( 14)

ẋ(t) = (A+BK)x(t). ( 15
)
It is clear that the stability of the latter system is equivalent to the stability of its equivalent descriptor form

ẋ(t) = y(t), 0 = -y(t)+(A+BK)x(t), (16)
which coincides with [START_REF] Fridman | Use of models with aftereffect in the problem of design of optimal digital control[END_REF] for h = 0. It is well-known ( [START_REF] Takaba | A generalized lyapunov theorem for descriptor systems[END_REF]) that the stability of the latter system is equivalent to the condition Ψ 0 < 0. If there exists P of the form (9c,d) which satisfies Ψ 0 < 0, then for small enough h > 0 LMIs of Lemma 2.1 are feasible (take e.g. Z = I 2n and R = [0 K T B T ]P T P [0 K T B T ] T ). We, therefore, obtain the following result: Corollary 1. If the continuous-time state-feedback ( 14) stabilizes the linear system (2), then the sampled-data state-feedback (3) with the same gain K stabilizes (2) for all small enough h.

In the case where the matrices of the system are not exactly known, we denote Ω = A B and assume that Ω ∈ Co{Ω j , j = 1, ...N }, namely,

Ω = N j=1 f j Ω j for some 0 ≤ f j ≤ 1, N j=1 f j = 1, ( 17
)
where the N vertices of the polytope are described by

Ω j = A (j) B (j) .
In order to guarantee the stability of (2) over the entire polytope one can use the result of Lemma 2.1 by applying the same matrices P 2 and P 3 for all the points in the polytope and solving (6a,b) for the N vertices only. A quadratic stability type criterion is then obtained:

Corollary 2. Given a gain matrix K, the system ( 5) is stable, over the entire polytope Ω, if there exist n×n matrices 0 < P

(j) 1 , P 2 , P 3 , Z (j) 1 , Z (j) 2 , Z (j) 3
and R (j) > 0 that satisfy the following LMIs:

Ψ (j)
1 < 0, and

R (j) [0 K T B (j)T ]P (j) * Z (j) ≥ 0, j = 1, ..., N, (18)
where

P (j) = P (j) 1 0 P 2 P 3 , Z (j) = Z (j) 1 Z (j) 2 * Z (j) 3 , Ψ (j) 1 = Ψ (j) 0 + hZ (j) + 0 0 0 hR (j) , Ψ (j)
0 = P (j)T 0 I A (j) + B (j) K -I + 0 I A (j) + B (j) K -I T P (j) .

Quadratic Stabilization

LMIs of Lemma 2.1 are bilinear in P and K. In order to obtain LMIs we use P -1 . It is obvious from the requirement of 0 < P 1 , and the fact that in ( 6)

-(P 3 + P T
3 ) must be negative definite, that P is nonsingular. Define

P -1 = Q = Q 1 0 Q 2 Q 3 and ∆ = diag{Q, I} (19)
Applying Schur formula to the term hR in (6a), we multiply (6a,b) by ∆ T and ∆, on the left and on the right, respectively. Denoting R = R -1 and Z = Q T ZQ we obtain, similarly to [START_REF] Fridman | An improved stabilization method for linear timedelay systems[END_REF], the following Theorem 1. The control law of (3) stabilizes ( 2) for all the samplings with the maximum sampling interval not greater than h and for all the system parameters that reside in the uncertainty polytope Ω, if there exist:

Q 1 > 0, Q (j) 2 , Q (j) 3 , R, Z(j) 1 , Z(j) 2 , Z(j)
3 ∈ R n×n , Ȳ ∈ R q×n that satisfy the following LMIs:

   Q (j) 2 +Q (j)T 2 + h Z(j) 1 Ξ(j) hQ (j)T 2 * -Q (j) 3 -Q (j)T 3 + h Z(j) 3 hQ (j)T 3 * * -h R    < 0,
and nonlinear matrix inequalities

   Q 1 R-1 Q 1 0 Ȳ T B (j)T * Z(j) 1 Z(j) 2 * * Z(j) 3   ≥ 0, (20)
where

Ξ(j) = Q (j) 3 -Q (j)T 2 +Q 1 A (j)T + h Z(j) 2 + Ȳ T B (j)T , j = 1, 2, ..., N. (21)
The state-feedback gain is then given by

K = Ȳ Q -1 1 . (22)
For solving [START_REF] Niculescu | Delay effects on stability: A robust control approach[END_REF] there exist two methods. The first uses the assumption

R = Q 1 , > 0, (23)
and thus leads to 2N LMIs with tuning parameter :

   Q (j) 2 +Q (j)T 2 + h Z(j) 1 Ξ(j) hQ (j)T 2 * -Q (j) 3 -Q (j)T 3 + h Z(j) 3 hQ (j)T 3 * * -hQ 1    < 0, (24)
   Q 1 0 Ȳ T B (j)T * Z(j) 1 Z(j) 2 * * Z(j) 3   ≥ 0, (25)
where Ξ(j) and j are given by [START_REF] Oishi | A bound of conservativeness in sampled-data robust stabilization and its dependence on sampling periods[END_REF].

Similarly to Corollary 2.2 we can show that if the system (2) is quadratically stabilizable by a continuous-time state-feedback [START_REF] Hale | Introduction of functional differential equations[END_REF], then for all small enough h the latter LMIs are feasible and the sampled-data state-feedback with the same gain stabilizes the system.

The second method for solving the matrix inequalities of Lemma 2.4 is based on the iterative algorithm developed recently by Gao and Wang (2003). This method is preferable in the cases of comparatively large h, since it leads to less conservative results. However it may take more computer time due to iterative process. In the sequel we shall adopt the first method for solving the matrix inequalities of Lemma 2.4.

Example 1. We consider (2) with the following matrices:

A = 1 0.5 0 -1 , B = 1 -1 .
It is verified by using Theorem 2.4 that the system is stabilizable by a sampleddata state-feedback with the maximum sampling interval h ≤ 0.69. Thus, for h = 0.69 the resulting K = [-1.048 0.2511] (with = 0.34). Simulation results (for uniform samplings with the sampling interval less than 0.7) show that the closed-loop system is stable.

Parameter Dependent Stabilization

The requirement for the quadratic stabilization imposes a serious constraint on the solution, where the same matrices Q 1 should satisfy the matrix inequalities in all the vertices of the polytope. To alleviate this difficulty a parameter dependent solution with different matrices

Q (j) 1
was derived in [START_REF] Fridman | Parameter dependent stability and stabilization of uncertain time-delay systems[END_REF] for the case of state delay. We modify the results of Corollary 4 of [START_REF] Fridman | Parameter dependent stability and stabilization of uncertain time-delay systems[END_REF] for the case of input delay by assuming that R = G 1 and obtain Theorem 2. Consider the system (2) and assume that its parameters lie in the polytope Ω. The system is stabilized, over the entire polytope Ω, by the controller of (3), for all the samplings with the maximum sampling interval not greater than h, if for some tuning positive scalar parameters and α there exist 2n×2n matrices: Q j , G j and H j , of the form

Q j = Q (j) 1 0 Q (j) 2 Q (j) 3 , G j = G 1 0 G (j) 2 G (j) 3 , H j = αG 1 0 H (j) 2 H (j) 3
.

and Zj , j = 1, ..., N , a m × n matrix Ȳ and n × n matrix R that satisfy the following LMIs.

     M j Q (j)T 2 Q (j)T 3 Q T j -G T j + Ā(j) H j + 0 αB (j) Ȳ 0 * -h -1 G 1 0 * * -H T j -H j      < 0 (26a) G 1 [0 Ȳ T B (j)T ] * Zj > 0, j = 1, ..., N. (26b)
where

M j = G T j Ā(j)T + Ā(j) G j + Ȳ T 0 0B (j)T + 0 B (j) Ȳ 0 +h Zj , Ā(j) = 0 I A (j) -I .
The state-feedback gain that stabilizes the system over Ω is then given by K = Ȳ G -1 1 .

Example 2. We consider (2) with the following matrices taken from [START_REF] Fridman | Regional stabilization and control of timedelay systems with saturating actuators[END_REF], where h = 0:

A = 0 1 -1 + g 1 -.5 , B = -1 + g 2 1
and where |g 1 | ≤ 0.53 and |g 2 | ≤ 1.7. It was shown in [START_REF] Fridman | Parameter dependent stability and stabilization of uncertain time-delay systems[END_REF] that the system is not quadratically stabilizable by continuous state-feedback (i.e. for h = 0). It is verified by using Theorem 2.5 that the system is stabilizable by a sampleddata state-feedback with the maximum sampling interval h ≤ 0.299. Thus, for h = 0.299 the resulting K = [0.0821 -0.1487] (with = 3.56 and α = 1.1). Simulation results (see e.g. Fig. 1 for the case of g 1 = g 1 (t) = 0.53 sin t, g 2 = g 2 (t) = 1.7 cos t, t k+1 -t k = 0.299 ∀k ≥ 0 and the initial condition x(0) = [5 -5] T ) show that the resulting closed-loop solutions converge to origin.

3 Regional Stabilization by Sampled-Data Controller with Saturation

Problem Formulation

Consider the system (2) with the sampled-data control law (3) which is subject to the following amplitude constraints

|u i (t)| ≤ ūi , 0 < ūi , i = 1, ..., m (27)
Denote by x(t, x(0)) the state trajectory of (2) with the initial condition x(0) ∈ R n . Then the domain of attraction of the origin of the closed-loop system (2), (3) is the set Lemma 3. Given β > 0, assume that there exists H in R m×n such that |h i x| ≤ ūi for all x(t) ∈ X β . Then for x(t) ∈ X β the system (30) admits the following representation.

ẋ(t) = Ax(t) + 2 m j=1 λ j (t)A j x(t -τ (t)) (32)
where

A j = B(D j K + D - j H) j = 1, ..., 2 m , 2 m j=1 λ j (t) = 1, 0 ≤ λ j (t), ∀ 0 < t, (33)
We denote

Ω α = 2 m j=1 λ j Ω j for all 0 ≤ λ j ≤ 1, 2 m j=1 λ j = 1 ( 34
)
where the vertices of the polytope are described by Ω j = A j , j = 1, ..., 2 m . The problem becomes one of finding X β and a corresponding H such that |h i x| ≤ ūi , i = 1, ...2 m for all x ∈ X β and that the state of the system

ẋ(t) = Ax(t) + A j x(t-τ (t)), τ (t) = t -t k , t k ≤ t < t k+1 , (35)
remains in X β .

Regional Stabilization

Applying the descriptor model transformation and the Lyapunov-Krasovskii functional of (8) by using the first method for solving the stabilization matrix inequalities (with tuning parameter ) , we obtain the following result:

Theorem 3. Consider the system (2) with the sampled-data control law (3) which is subject to the constraints (27). The system is stable with X β inside the domain of attraction for all the samplings with the maximum sampling interval not greater than h, if there exist 0

< Q 1 , Q (j) 2 , Q (j) 3 , Z (j) 1 , Z (j)
2 , Z

∈ R n×n , Y, G ∈ R m×n and β > 0 that satisfy the following set of inequalities:

   Q (j) 2 +Q T (j) 2 +hZ (j) 1 Σ j hQ (j) 2 * -Q (j) 3 -Q T (j) 3 +hZ (j) 3 hQ (j) 3 * * -hQ 1    < 0, j = 1, ..., 2 m (36a)    Q 1 0 (Y T D j + G T D - j )B T * Z (j) 1 Z (j) 2 * * Z (j) 3    ≥ 0 (36b) β g i * ū2 i Q 1 ≥ 0, i = 1, ..., m, (37)
where

Σ j = Q (j) 3 -Q T (j) 2 + Q 1 A T +(Y T D j +G T D - j )B T +hZ (j) 2 . ( 38
)
The feedback gain matrix which stabilizes the system is given by K = Y Q -1 1 .

Proof: For V given by ( 8) conditions are sought to ensure that V < 0 for any x(t) ∈ X β . As in [START_REF] Fridman | Regional stabilization and control of timedelay systems with saturating actuators[END_REF], the inequalities (37) guarantee that |h i x| ≤ ūi , ∀x ∈ X β , i = 1, ..., m, where g i ∆ = h i Q 1 , i = 1, ..., m and Q 1 ∆ = P -1 1 , and the polytopic system representation of (35) is thus valid. Moreover, (36a,b) guarantee that V < 0.

From V < 0 it follows that V (t) < V (0) and therefore for the initial conditions of the form (31)

x T (t)P 1 x(t) ≤ V (t) < V (0) = x T (0)P 1 x(0) ≤ β -1 .

(

) 39
Then for all initial values x(0) ∈ X β , the trajectories of x(t) remain within X β , and the polytopic system representation (35) is valid. Hence x(t) is a trajectory of the linear system (35) and V < 0 along the trajectories of the latter system which implies that lim t→∞ x(t) = 0.

Example 3. We consider (2) with the following matrices (taken from [START_REF] Cao | Stability analysis of linear time-delay systems subject to input saturation[END_REF], where h = 0):

A = 1.1 -0.6 .5 -1 , B 1 = 1 1
and where ū = 5. Applying Theorem 3.3 a stabilizing gain was obtained for all samplings with the maximum sampling interval h ≤ 0.75. In order to 'enlarge' the volume of the ellipse we minimized the value of β (to improve the result we also added the inequality Q 1 > αI and chose such α > 0 that enlarged the resulting ellipse). The ellipse volume increases when h decreases (see Figure 2). , α = 1) and we show (see Figure 3) that a trajectory starting on

Fig. 1 .

 1 Fig. 1. State trajectories of the closed-loop system for h = 0.299, g1 = 0.53 sin t, g2 = 1.7 cos t and x(0) = [5 -5] T

  For, say, h = 0.75 we obtain K = [-1.6964 0.5231] (with = 0.325, β = 0.1261, P 1 = 0.9132 -0.2816 -0.2816 0.0868

Fig. 2 .

 2 Fig. 2. Ellipsoidal bounds on the domain of attraction: line corresponds to h = 0.1; slash line to h = 0.75; point line to h =

A = {x(0) ∈ R n : lim t→∞ x(t, x(0)) = 0}.

We seek conditions for the existence of a gain matrix K which leads to a stable closed-loop. Having met these conditions, a simple procedure for finding the gain K should be presented. Moreover, we obtain an estimate X β ⊂ A on the domain of attraction, where

and where β > 0 is a scalar and P 1 > 0 is an n × n matrix. We represent the state-feedback in the delayed form

Reducing the original problem to the problem with input delay, we solve it by modifying derivations of [START_REF] Fridman | Regional stabilization and control of timedelay systems with saturating actuators[END_REF], where the case of state delay was considered.

A Linear System Representation with Polytopic Type Uncertainty

Applying the control law of (29) the closed-loop system obtained is

Though the closed-loop system has a delay, we keep in mind that in the case of sampled-data control the initial condition is defined in the point t = 0 and not on the segment [-h, 0]. That is why for the estimation of the domain of attraction we can restrict ourself to the following initial functions φ(s), s ∈ [-h, 0]:

Denoting the i-th row by k i , we define the polyhedron

If the control and the disturbance are such that x ∈ L(K, ū) then the system (30) admits the linear representation. Following [START_REF] Cao | Stability analysis of linear time-delay systems subject to input saturation[END_REF], we denote the set of all diagonal matrices in R m×m with diagonal elements that are either 1 or 0 by Υ , then there are 2 m elements D i in Υ , and for every i = 1, ...,

The following is obtained from Lemma 3.1. the periphery of the ellipse (for the case of the uniform sampling with the sampling period t k+1 -t k = 0.75) never leaves this ellipse and converges to the origin, while a trajectory starting not far from the ellipse remains outside the ellipse.

Conclusions

A new method for robust sampled-data stabilization of linear continuous-time systems is introduced. This method is based on the continuous-time model with time-varying input delay. Under assumption that the maximum sampling interval is not greater then h > 0, the h-dependent sufficient LMIs conditions for stabilization of systems with polytopic type uncertainty and for regional stabilization of systems with sampled-data saturated state-feedback are derived via descriptor system approach to time-delay systems.

The new approach solves the problems for comparatively small h and leads to sufficient conditions only, however these conditions are simple. The method may be applied to a wide spectrum of robust sampled-data control problems (e.g. to guaranteed cost or to H ∞ control).


François Veysset
email: francois.veysset@onera.fr

Lot Belkoura
email: lotfi.belkoura@univ-lille1.fr

Patricia Coton
email: patricia.coton@onera.fr

Jean-Pierre Richard
email: jean-pierre.richard@ec-lille.fr

DELAY SYSTEM IDENTIFICATION APPLIED TO THE LONGITUDINAL FLIGHT OF AN AIRCRAFT THROUGH A VERTICAL GUST

Keywords: Dynamics of ight, Delay systems, Identication

This paper deals with modelling and identication of aircraft dynamic entering a vertical gust. The identication approach initiated in (Fliess 2003) falls under a prospect for identication from tests carried out in the Flight Analysis Laboratory of the DCSD of ONERA in Lille. The plane is considered into various elements which consist in the fuselage, the wing and the tail. The model incorporates delays linked to the aircraft passage through the atmospheric turbulence.

INTRODUCTION

Flight through atmospheric turbulence is one of the signicant subjects of research in aeronautics. This topic concerns at the same time manufacturers of planes and operators of civil transport. The development of a representation of the unsteady aerodynamic phenomena as well as the establishment of command laws for gust alleviation are signicant objectives to ensure safety in critical phases (takeo and landing) but also to ensure the comfort of the passengers.

Our study relates to the mathematical representation of the behaviour of the plane in atmospheric turbulence, and to the identication of the dynamics induced by the gust. A modelling based on ight mechanics will make it possible to simulate its ight. The fast identication technique proposed in [START_REF] Fliess | An algebraic framework for linear identication[END_REF] is used here as a rst step of the approach for identication based on tests carried out in the Flight Analysis Laboratory of the DCSD of ONERA in Lille. During those experiments, a model of civil aircraft equipped with an embedded instrumentation will be catapulted and will cross, during its free ight, a turbulence generated by a vertical blower.

Various approaches of the modelling of the ight of a plane in atmospheric turbulence and of the representation of the unsteady eects were considered [START_REF] Van Staveren | Analyses of aircraft responses to atmospheric turbulence[END_REF], [START_REF] Klein | Modeling of aircraft unsteady aerodynamic characteristics[END_REF]). Among them, a description including terms of delays was introduced [START_REF] Coton | Contribution au DTP Turbulence -Modélisation du vol en turbulence[END_REF]), [START_REF] Coton | DTP Turbulence -Identication du modèle aérodynamique instationnaire[END_REF], [START_REF] Jauberthie | Méthodologies de planication d'expériences pour systèmes dynamiques[END_REF], in order to get a better accuracy of the representation of the eects of the penetration of the aircraft through the gust. This article falls under the continuity of the work initiated by [START_REF] Coton | Contribution au DTP Turbulence -Modélisation du vol en turbulence[END_REF], and a more detailed attention is related to the modelling and identication of the contribution of the gust to the aerodynamic coecients. It is known that delays create a specic diculty for the identication (see [START_REF] Richard | Time delay systems : An overview of some recent advances and open problems[END_REF][START_REF] Belkoura | Identiability and identication of linear systems with delays[END_REF][START_REF] References Belkoura | Identiability of systems described by convolution equations[END_REF]). Such distributed and delayed actions of the turbulence on the lifting surfaces is particularly considered here.

The paper is organized as follows. Section 2 introduces the approach via a delay model. Contribution of the various elements of the plane to the aerodynamic coecients is presented in section 2.1, followed by the description of the dynamic induced by the gust in 2.2. An identication method for the gust transfer function based on simulated data is presented in 3.

MODELLING

The general model describing the behaviour of a plane in a longitudinal ight results from the fundamental principle of dynamics and is governed by:

m V = -mg sin (θ -α) - 1 2 ρSV 2 C x , mV α = mg cos (θ -α) + mV q - 1 2 ρSV 2 C z , B q = 1 2 ρSlV 2 C m , (1)
θ = q,
where V , α, q et θ represent respectively the speed, the kinematic angle of attack, the pitch rate and the pitch attitude of the plane. The constant parameters m, B, g, ρ, l and S correspond respectively to the mass of the plane, its inertia, the constant of gravity, the density of the air, a length and a surface of reference. Lastly, C x , C z and C m respectively represent the aerodynamic coecients of drag, lift and pitching moment. On the assumption of a ight with weak variation of angle of attack, the drag coecient generally results from the lift by a relation of the form:

C x = C x0 + C x1 C z + C x2 C 2 z , (2)
with known coecients C xi , and the main point of our study will be concentrated on the analysis of the lift and pitching coecients C z and C m . The contribution of the gust to these global coecients is considered separately through the relations:

C z = C plane z + C gust z , (3)
C m = C plane m + C gust m .
(4)

2.1 Contribution of the elements of the plane to the aerodynamic coecients Within the framework of the longitudinal ight, the plane is split up into three elements which consist in the fuselage, the wing and the horizontal tail. Expressions of the aerodynamic lift coecients (gust not considered) are obtained from integration, along each element, of local angles of attack, respectively weighted by the section variation of the fuselage (denoted dS dx ) in the case of the fuselage, and by local coecients depending on the surface of reference S, the local chord c (y) and the lift gradient C zα (y) of the surface considered (these coecients are noted C wing zα (y)

for the wing and C ht zα (y) for the horizontal tail). A similar method is followed for the pitching moment coecients where preceding weightings are multiplied by x l with x the X-coordinate of the element at the point of integration and l a length of reference. These various operations reveal a behaviour closely connected to the state in the form (see e.g [START_REF] Etkin | Dynamics of ight[END_REF]))

C plane z = C z0 + C z1 α + C z2 q/V, (5)
C plane m = C m0 + C m1 α + C m2 q/V. (6)

Contribution of the gust

In order to obtain an accurate description of the the aerodynamic coecients induced by the gust, terms of delays are introduced as shown in Figure 1. A distributed delay τ (y) is introduced in order to take the sweepback wing into account, while the delay τ 1 corresponds to the delayed eect of the gust on the horizontal tail. The fuselage eect is assumed negligible and we consider separately the wing and horizontal tail contributions to the aerodynamic coecients. We therefore write, with obvious notations,

C gust z = ζ wing z + ζ ht z , (7)
C gust m = ζ wing m + ζ ht m . (8)
2.2.1. Contribution on the wing The method adopted in this section consists in considering the distributed action of the gust on the wing induced by its sweep. In a matter of space saving, only lift is analyzed and relations relating to pitching moment are deduced as shown in section 2.1, by replacing the gradient C wing zα (y) by x l C wing zα (y)

with x the X-coordinate of the point considered and l the length of reference.

The modelling of the unsteady eects induced on a airfoil by its penetration through a vertical gust was introduced by Küssner [START_REF] Kussner | Stresses produced in airplane wings by gusts[END_REF] and Theodorsen [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of utter[END_REF]. Within this framework, the contribution of the gust to the lift on the wing is expressed through the relation:

ζ wing z = b/2 -b/2 C wing zα (y) k(t, y) * w(t, y) V (t) dy (9)
where b is the span of the wing, * denotes the convolution product, w(t, y) is the local gust value, and k(t, y) is the Küssner function. The latter may be regarded as the impulse response of a transfer function K(s, y) which depends on the chord c(y) of the wing (see gure 1), the average velocity V m and a dominant mode a according to the relation:

2K(s, y) = 1 + aV m /c(y) s + aV m /c(y) (10)
The distributed delay τ (y) = 2la

Vmb y = τ y (see gure 1) describing the sweepback wing eects is particularly considered. The approximation of the speed V (t) by its average value V m results in a input w(t, y) = w(t -τ y) and allows us to formulate the gust contribution to the lift coecient of the wing as:

V m ζ wing z (s) = (Φ w (s) + Φ w z (s)) w(s), (11)
where

Φ w (s) = b/2 0 C wing zα (y) e -τ.ys dy, (12)
and

Φ w z (s) = b/2 0 C wing zα (y)
e -τ.ys 1 + c(y) aVm s dy.

(

) 13
The main diculty lies in the fact that Φ w z (s)

does not admit an explicit primitive. Even in the restrictive case where the sweepback wing is neglected and the local coecient C wing zα (y) is assumed constant, integration of (13) results in a transfer function of the form log( 1+λ1 s 1+λ2 s ) not easily workable in identication or control perspectives. Two successive approximations are then carried out.

(a) The gust is supposed to be relatively smooth in comparison to the dimensions of the wing so the distributed delay is approached by e -τ.ys ≃ 1 -τ y s.

(b) The Laplace transform of the impulse response Φ w z (t) is also approximated by a Pade's development of order 1/1 with a corrected static gain.

Let us note that although approximation (a) is getting worst as one moves away from the fuselage, it should be attenuated by the decrease of the local chord c (y) and the lift gradient C zα (y), like generally conrmed by wind tunnel tests. The correction listed in (b) is due to the development in 1/s (and not in s) for which we wish to preserve the static gain Φ w z (0) = Γ w . The obtained approximations are therefore given by:

Φ w (s) ≃ Γ w -sΓ 1 , (14)
Φ w z (s) ≃ Γ w - α 0 s s + α 2 /α 1 , (15)
where results from experiments carried out in wind tunnel. This result consolidates the approximation carried out, although a more accurate approximation is possible by complexifying the model.

α i = b/2 0 aV m c (

Contribution on the horizontal tail Al-

though the method is similar to that used for the wing, reduced dimensions of the tail allow a simpler description of the dynamics induced by the gust. In particular, the chord c (y) is estimated at its average value and the sweepback tail is neglected. This yields the delayed transfer function:

V m ζ ht z (s) = Φ ht z (s) w(s), (17)
Φ ht z (s) = Γ ht s + 2aV m /c m s + aV m /c m e -τ1s . (18)
2.2.3. First order approach In order to evaluate the performances of our identication technique, we consider, in a rst step, a rst order approximation of relation ( 11) for the contribution on the wing. This expression can be reformulate, considering K(s) = 1 ks+1 and using approximation (a), as : to a crenel gust considering the dominant mode a (Eq.11) and the rst order (Eq.19). It clearly conrms the approximation carried out.

V m ζ wing z (s) = -s2Γ 1 /k + 2Γ w /k s + 1/k w(

IDENTIFICATION

During the tests in the Flight Analysis Laboratory of the DCSD of ONERA, the state and state derivative in equation ( 1) are available using on the one hand, embedded gyrometers, accelerometers and Kalman-Rauch ltering, and on the other hand, an optical trajectory system based on two sets of ten video cameras spaced along the facility. Moerover, by means of wind tunnel experiments, the aerodynamic coecients associated to the elements of the plane are also estimated. Therefore, the gust aerodynamic coecients are available.

To avoid repetition, only the transfer function of lift coecient C gust z is considered for the identication, and by virtue of the developments of the previous section, it admits an expression of the form:

y = k w0 + k w1 s 1 + τ w s + k t0 + k t1 s 1 + τ t s e -τ s w (20)
where for ease of notations we have denoted:

y = C gust z = C z -C z0 -C z1 α -C z2 q/V, k w0 = 2Γ w /V m , k w1 = -2Γ 1 /V m , k t0 = 2Γ ht /V m , k t1 = Γ ht c m /(aV 2 m ), τ w = k, τ t = c m /(aV m ).
The main diculty in this problem lies in the presence of the unknown delay τ . We shall focus in this section on the identication of the parameters τ , τ w , and τ t from step responses. As we shall see, the remaining coecients need not be known and could be identied in a second step. The approach used here is based on the work initiated in [START_REF] Fliess | An algebraic framework for linear identication[END_REF] and extended in our paper to parameters and delay identication.

Mathematical framework

Functions are considered through the distributions they dene and are therefore indenitely dierentiable. If y is a continuous function except at a point a with a nite jump σ a , its derivative dy/dt writes dy/dt = ẏ + σ a δ a (21) where ẏ is the distribution dened from the usual derivative of y. Derivation, integration and translation can be formed from the convolution products

ẏ = δ (1) * y, y = H * y, y(t-τ ) = δ τ * y (22)
where δ (1) is the derivative of the Dirac distribution, and H denotes the Heaviside function. With a slight abuse of notations, we shall write H k y the iterated integration of y and more generally T k the iterated convolution product of order k. A distribution is said to be of order r if it acts continuously on C r -functions but not on C r-1functions. Measures and functions are of order 0.

The multiplication of two distributions (say α and T ) always make sense when one of the two terms is a smooth function. Particularly, when T is a Dirac derivative of order n, on has:

α δ (n) = n k=0 (-1) (n-k) C k n α (n-k) δ (k), (23)
The next Theorem is the key result from which most of the parameters (including the delays) can be identied from step input responses.

Theorem 1. [START_REF] Schwartz | Théorie des distributions[END_REF] If a distribution T has a compact support K and is of order m (necessarily nite), α T = 0 whenever α and its derivatives of order ≤ m vanish on K.

The following examples illustrate this statement in case α is a polynomial and T a singular distribution. Note that, in forming the product α T , the delay τ involved in the argument T (t -τ ) now appears also as a coecient.

t δ = 0, (1 -e -γt ) δ = 0, (1 -e -γ(t-τ ) ) δ τ = 0. (24)
We shall make use of another property involving both multiplication with e -γt and the convolution product, in case one of the two distributions (S or T ) has a compact support.

e -γt (S * T ) = e -γt S * e -γt T.

(

) 25
With S = δ (p) and T = y this equation allows us to transform terms of the form e -γt y (p) into linear combinations of derivatives of products e -γt y.

Denoting z = e -γt y, one has for example, e -γt y (2) = γ 2 z + 2γ z (1) + z (2) .

(26)

Note that integrating twice this expression by considering H 2 e -γt y (2) results in nothing but the integration by parts with available data z.

Application to the gust transfer identication

With a step gust w = w 0 H, a rst order derivation of the dierential equation induced by (20) results in a relation of the form

a 2 y (3) + a 1 y (2) + y (1) = ϕ 0 + β 0 δ + β 1 δ (1) + β 2 δ τ + β 3 δ (1) τ ( 27
)
where a 2 = τ t τ w , a 1 = τ t + τ w , ϕ 0 (of order 2 and support {0}) contains the initial condition terms (discontinuities of y at t = 0) that naturally appear in the distributional framework, and β i are combinations of the unknown parameters.

By virtue of Theorem 1, the right hand side of equation ( 27) can be cancelled by means of a multiplication with a function α such that α (k) (0) = α (k) (τ ) = 0 for k = 0, 1, 2. The choice of the function

α(t) = (1 -e -γt ) 3 (1 -e -γ(t-τ ) ) 3
(28) results in

(1-e) 3 (1-λ e) 3 (a 2 y (3) +a 1 y (2) +y (1) ) = 0, (29) where for ease of notations we denoted e = e -γt and λ = e γτ . Note that the latter equation no longer requires the knowledge of the gains k w0 , k w1 , k t0 and k t1 . As an equality of singular distributions, this relation doesn't make sense for any t (otherwise we would have τ = t). However, k ≥ 1 successive integrations (or a convolution with H k ) result in functions equality from which the delay τ and the coecients τ w , τ t become available. More precisely, equation ( 29) combined with integrations leads to the following formulation:

(A 0 + λ A 1 + λ 2 A 2 + λ 3 A 3 )   a 2 a 1 1   = 0, (30)
where the matrices entries described below are realized by means of integration by part formula described previously.

A 0 (i, j) = +1H i+3 (1 -e) 3 y (4-j)
A 1 (i, j) = -3H i+3 e(1 -e) 3 y (4-j)

A 2 (i, j) = +3H i+3 e 2 (1 -e) 3 y (4-j)

A 3 (i, j) = -1H i+3 e 3 (1 -e) 3 y (4-j)

The identication problem is therefore transformed into a generalized eigenvalue problem (30) for which the delay is deduced from one of the eigenvalues (i.e. τ = log(λ)/γ), while the coecients a 2 and a 1 are obtained from the corresponding normalized eigenvector.

Figure 4 shows in simulation the step response of transfer to be identied with the numerical values τ w = 0.6, k w0 = 2, k w1 = 0.5, τ t = 0.4, k t0 = 0.7, k t1 = 0.1, τ = 0.5 and a time scale of 1/10. Figure 5 shows the generalized eigenvalues of (30) (and more precisely log(λ)/γ) solved using the polyeig Matlab function and the parameter (in (28)) γ = 0.2. The result clearly shows one constant eigenvalue corresponding to the unknown delay τ . This value appears after a transitory phase for which the delay is not identiable. Figure 6 shows the second and third component of the corresponding normalized eigenvector. These clearly converge to the desired values a 2 = τ w τ t and a 1 = τ w + τ t . Moreover, and although these coecients are not identiable for t < τ , they seem to converge in a rst step to the pair (a 2 , a 1 ) for which τ w = 0. Note also that a singularity occurred for t ≃ 1s. This paper illustrated new identication and modelling approaches and a rst step towards the identication of the aerodynamics induced by a gust in the longitudinal ight of an aircraft. Robustness issues, higher order of the unsteady aerodynamics, fuselage eects as well as additional delays taking into account the deection term in the expression of the angle of attack of the horizontal tail are under investigation.

Fig. 1 .

 1 Fig. 1. Description of the delays in the modelling.

  with the constant terms Γ w = b/2 0 C wing zα (y) dy and Γ 1 = b/2 0 C wing zα (y) τ y dy and

Figure 2 Fig. 2 .

 22 Figure 2 represents in simulation the theoretical (Eq.13) and approached (Eq.15) responses of the transfer Φ w z (s) to a crenel gust. The data correspond to the model and the experiments under consideration at the Flight Analysis Laboratory of the DCSD of ONERA: b = 2.23 m, V m = 23.25 m/s, τ = 0.023s.

Fig. 3 .

 3 Fig. 3. Responses of V m ζ wing z

Fig. 4 .

 4 Fig. 4. Step response of the gust transfer

Fig. 5 .

 5 Fig. 5. Eigenvalues of problem (30)

  Fig. 6. Eigenvector of the constant eigenvalue