Aïda Ouangraoua

Frédéric Boyer
email: frederic.boyer@cea.fr

Andrew Mcpherson

Éric Tannier
email: eric.tannier@inria.fr.

Cedric Chauve

Prediction of Contiguous Regions in the Amniote Ancestral Genome

We investigate the problem of inferring contiguous ancestral regions (CARs) of the genome of the last common ancestor of all extant amniotes, based on the currently sequenced and assembled amniote genomes as ingroups and three teleost fish genomes as outgroups. We combine a methodological framework using conserved syntenies computed from whole genome alignments of amniote species together with double conserved syntenies (DCS) using gene families from amniote and fish genomes, to take into account the whole genome duplication that occurred in the teleost lineage. From these comparisons, ancestral genome segments are computed using techniques inspired by physical mapping. Due to the difficulty caused by the whole genome duplication and the large evolutionary distance to the closest assembled outgroup, very few methods have been published with a reconstruction of the amniote ancestral genome. This one is the first which is founded on a simple and formal methodological framework, whose good stability is shown and whose CARs cover large regions of the human and chicken genomes.

Introduction

The reconstruction of ancestral karyotypes and gene orders from homologies between extant species can help to understand the large-scale evolutionary mutations that differentiate the present genomes. It has been approached using cytogenetics methods and recently applied to mammalian genomes [START_REF] Wienberg | The evolution of eutherian chromosomes[END_REF]. Beyond this evolutionary distance, homologies are less visible and it is only with the recent availability of sequenced and assembled genomes that bioinformatics methods can predict the past of chromosomes. These methods address the problem at a much higher resolution, although with much less available genomes. The first results have been obtained on mammalian genomes [START_REF] Bourque | Reconstructing the genomic architecture of ancestral mammals: lessons from human, mouse and rat genomes[END_REF][START_REF] Murphy | Dynamics of mammalian chromosome evolution inferred from multispecies comparative maps[END_REF]16], and several reviews have been published [START_REF] Faraut | Addressing chromosome evolution in the whole-genome sequence era[END_REF][START_REF] Muffato | Paleogenomics, or the recovery of lost genomes from the mist of times[END_REF], analyzing the divergences with earlier cytogenetics results [START_REF] Froenicke | Are molecular cytogenetics and bioinformatics suggesting diverging models of ancestral mammalian genomes?[END_REF][START_REF] Bourque | The convergence of cytogenetics and rearrangement-based models for ancestral genome reconstruction[END_REF][START_REF] Rocchi | Ancestral genome reconstruction: An integrated, multi-disciplinary approach is needed[END_REF]. These methods can be divided into model-based methods, that compute complete evolutionary scenarios [START_REF] Bourque | Reconstructing the genomic architecture of ancestral mammals: lessons from human, mouse and rat genomes[END_REF][START_REF] Murphy | Dynamics of mammalian chromosome evolution inferred from multispecies comparative maps[END_REF] and model-free approaches that do not consider a precise rearrangement model, which are used by cytogeneticians and currently receive a lot of attention from computational biology (see [16,[START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF] and references there).

The application of such methods to more ancient genomes comes up against the difficulty to handle duplications and losses as evolutionary events. Yet the teleost fish genomes have undergone a whole genome duplication (WGD) [START_REF] Jaillon | Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF][START_REF] Muffato | Paleogenomics, or the recovery of lost genomes from the mist of times[END_REF] at an early stage of their evolution, and are currently the only available genomes that may serve as an outgroup to reconstruct amniote or tetrapod ancestral genomes. Two recent methods have been developed to reconstruct these ancestral genomes [START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF][START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF], and predict very divergent syntenic associations. Hence, while the reconstruction of the ancestral mammalian genome seems now to be close to a relative consensus, the reconstruction of the amniote ancestral genome looks as the next bottleneck on the way towards the ancestral proto-vertebrate genome.

In [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], a general model-free framework was introduced for the reconstruction of "Contiguous Ancestral Regions", or CARs (the terminology is borrowed from Ma et al. [16]) in an ancestral genome. It is inspired by genome physical mapping techniques, and roughly consists in two phases as follows. Given a set of "genomic markers", which are sets of orthologous positions in the ingroup genomes: (1) detect "ancestral syntenies", which are sets of genomic markers that are believed to have been contiguous in the ancestral genome, and (2) order the genomic markers into a set of "Contiguous Ancestral Regions" in which the ancestral syntenies are respected, discarding some of them if the whole set is not compatible with the formation of linear CARs. This second phase relies on combinatorial tools such as PQ-trees, that were introduced in computational biology for physical mapping of genomes. Indeed, our problem consists in the mapping of markers into ancestral chromosomes. This framework was applied in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF] for the reconstruction of contiguous ancestral regions of mammalian ancestors (ferungulates and boreoeutheria). It was shown to be very stable under different parameters for the computation of syntenies.

Our goal here is to apply this framework to compute CARs of the ancestral amniote genome. The method used to infer mammalian CARs needs to be extended to handle two main issues. First, the closest currently available sequenced and assembled outgroups are the teleost fishes, whose evolutionary distance to the ingroups (mammals and birds) is considerable. Hence it is impossible to define a high coverage set of genomic markers that appear once in each genome of this study. Moreover, the Whole Genome Duplication followed by massive gene losses and intensive rearrangements makes ancestral syntenies inaccessible by a classical comparison between amniote and fish genomes.

We handle these issues by using (1) genomic markers obtained from wholegenome alignments within amniote assembled genomes (chicken and mammals), and (2) gene families to compare amniote and fish genomes (teleost fishes) to construct ancestral syntenies. We rely on the Double Conserved Synteny (DCS) principle introduced by Kellis et al. [START_REF] Kellis | Proof and evolutionary analysis of ancient genome duplication in the yeast saccharomyces cerevisiae[END_REF] and Dietrich et al. [START_REF] Dietrich | The ashbya gossypii genome as a tool for mapping the ancient saccharomyces cerevisiae genome[END_REF], since then often used to detect syntenies in a WGD context [START_REF] Jaillon | Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF][START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF], and systematized by Van de Peer [START_REF] Van De Peer | Computational approaches to unveiling ancient genome duplications[END_REF]. Despite the principle is well known, the detection of such syntenies implies a non trivial methodological problem and formal descriptions of the expected signal, adapted to the highly rearranged fish genomes, are lacking. It is a contribution of this paper to propose a formal definition of DCSs.

We obtain a family of ancestral syntenies, and group them into Contiguous Ancestral Regions of the proto-amniote genome. This set can contain some "conflicting signal", which means that no linear ordering of the genomic markers can account for all the ancestral syntenies. While the ancestral syntenies computed from mammalian genomes presented very little conflict [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], the greater evolutionary distance mixes up the signal and we have to cope with more conflicting ancestral syntenies. That is why we produce here the sets of CARs with different sets of parameters used to compute the DCS and propose both a set of CARs obtained with stringent parameters and a set of consensus CARs obtained from several values of parameters. We compare our results with two other studies [START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF][START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF] that proposed a configuration of the amniote ancestral genome. These two present contradictory results, low coverage of the extant genomes by the reconstructed ancestral one, and no validation of the methods. We try to make significant progresses in these directions: our CARs are more numerous than in the previous studies but present a good coverage of the extant genomes. As in the previous studies, we find a proto-amniote genome that shows more similarity to the chicken genome than to mammalian genomes.

In the following we first describe the method, following the framework of [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], and focusing on the novelty we introduce, which is the possibility of integrating duplicated syntenies in this framework. The definition and computation of the duplicated syntenies is discussed. Then we describe the CARs we obtain into details, comparing them to previous studies, showing some convergences and differences. Eventually we study the soundness of the proposed CARs by running the method under different sets of parameters, to understand its behavior, its stability and the confidence we may have on the proposed CARs.

Data and methods

Overview

We consider a dataset containing eleven amniote genomes (human, chimpanzee, orangutan, macaca, mouse, rat, dog, cow, horse, opossum, chicken) and three teleost fish genomes (tetraodon, stickleback, medaka) used as outgroups. We then proceed with the following steps.

1. We compute a set of "genomic markers" that are unique and universal in amniote genomes (i.e. each markers appears once and exactly once in each of the eleven genomes), using whole genome alignments available. 2. A first set of ancestral syntenies is generated by computing common intervals (as in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF]) of genomic markers between all pairs of amniote species whose evolutionary path goes through the amniote ancestor (here, chicken against every mammal). 3. A second set of ancestral syntenies is generated by computing "double conserved syntenies" (DCS) between each amniote and the three teleosts, using sets of gene families. This provides the coordinates on an amniote genome of a genomic segment which is likely to descend from a segment of the ancestral osteichthyes genome (the ancestor of teleost fishes and amniotes), and is then likely to have been present as an segment of the ancestral amniote genome. These coordinates provide a set of genomic markers (those which intersect the DCS) which is taken as the ancestral synteny. 4. We weight the ancestral syntenies from both sets according to their conservation pattern in the considered species tree. 5. We select from the ancestral syntenies a maximum weight subset such that alltogether they are compatible with the formation of linear CARs. This supposes that the least weight conflicting ancestral syntenies are more likely to be false positives. This phase relies on the "consecutive ones problem", widely used in physical mapping problems. The final result is a combinatorial structure (a PQ-tree) that allows to linearly represent the whole set of solutions to the consecutive ones problem. The children of the root of the PQ-tree are the amniote CARs.

Genomes and markers

We used data retrieved from the Compara database of Ensembl [8], comparing human (hg18), chimpanzee (CHIMP2.1), orangutan (PPYG2), macaca (Mmul 1) mouse (mm9), rat (RGSC 3.4), dog (CanFam2.0), cow (Btau 4.0), opossum (monDom5), platypus (Ornitorhynchus anatinus-5.0), chicken (WASHUC2), tetraodon (TETRAODON8.0), medaka (HdrR) and stickleback (BROAD S1).

For the phylogenetic relationships between the amniote species, we considered the species tree used by Compara to compute whole-genome alignments (see http://www.cecm.sfu.ca/~cchauve/SUPP/ISBRA09/ for data and results).

We construct a set of genomic markers by using the Pecan 12-amniotesvertebrates multiple alignments available in the release 50 of Compara with human genome as a reference. From the orthologous seeds on the 11 fully assembled genomes defined by the multiple alignments, we keep only the ones that have a minimum size (100b in the human reference genome). Two seeds are joined if they are distant from each other by less than 100Kb in all amniotes where they occur. A genomic marker is then an inclusionwise maximal set of linked seeds which spans more than 100Kb in all genomes and such that its seeds span at least 50% of its total span. In this way we expect to obtain a good set of orthologous markers, removing uncertain homologies as well as paralogies. We obtain 1101 non-overlapping genomic markers, spanning respectively 797Mb of the human genome (26% of its size) and 308Mb (29%) of the chicken genome.

Ancestral syntenies.

As described in Chauve et al. [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], ancestral synteny detection in the absence of large duplications (in mammalian and bird genomes) may be performed by the detection of groups of genomic markers contiguous in at least two genomes whose evolutionary path goes through the desired ancestor. But the WGD that is believed to have occurred in the lineage of the teleost genomes requires a more sophisticated treatment when a fish genome is involved in the comparison.

In the absence of large duplications: comparing two amniote genomes. First, when comparing the chicken genome to each mammalian genome, which are not separated by a whole genome duplication event, an ancestral synteny is the set of genomic markers intersecting (a) an inclusionwise maximal segment of the chicken genome that contains the same genomic markers as a segment of a mammalian genome (a maximal common interval) or (b) a segment of the chicken genome that contains only two markers that are also consecutive in a mammalian genome (an adjacency). We include adjacencies to balance the fact that there is no order information associated with common intervals.

The result of this process is a family of sets of markers, each of which covers amniote genomic segments that are believed to be ancestral.

In the presence of a WGD in the outgroups. For the comparisons between an amniote genome and a fish genome, the method described above can not be applied due to the numerous losses and intra-chromosomal rearrangements that shuffled the teleost fish chromosomes after the WGD. Indeed, massive losses of genes often follow a WGD, and two paralogous segments may present little similarity. It is possible to detect this paralogy indirectly by comparing the two segments with their common ortholog in a non duplicated genome. This principle is called the pivot method, and is now classical to detect chromosome segment homologies in a WGD context [START_REF] Kellis | Proof and evolutionary analysis of ancient genome duplication in the yeast saccharomyces cerevisiae[END_REF][START_REF] Dietrich | The ashbya gossypii genome as a tool for mapping the ancient saccharomyces cerevisiae genome[END_REF][START_REF] Jaillon | Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF][START_REF] Van De Peer | Computational approaches to unveiling ancient genome duplications[END_REF] (see Figure 1). Despite this principle is well known, no methodological discussion on the exact signal that should be detected using fish genomes has been published (due to the highly rearranged fish chromosomes, the approaches described in [START_REF] Van De Peer | Computational approaches to unveiling ancient genome duplications[END_REF] are not efficient in this case). So in the present method we propose a definition of a DCS. The proportion of conflicting signal in the whole set of DCS tends to show that there is still some space for improvement of the precision of this proposition.

We use a set of orthologous gene families constructed from the orthologies between genes of amniotes and fishes available in the Ensembl Compara database [8] 4 . When comparing an amniote genome to a fish genome, we use only genes that are annotated with coordinates on a chromosome on both species. A Fig. 1. A double synteny: an amniote chromosome (in the middle) is homologous to two fish segments (up and down), though paralogy between these two segments is not detectable through direct similarity (few genes are conserved in two copies). Clusters are much rearranged: this justifies the method of detection, which does not consider the order of the genes within a cluster.

double conserved synteny (DCS) is a segment S of the amniote genome that is orthologous to two paralogous segments S 1 and S 2 in two different chromosomes of the fish genome, i.e. satisfies the two following criteria.

1. S contains at least 20 genes, min prop percent of all genes having orthologs on S 1 or S 2 (we choose min prop = 95% and test the sensitivity to this parameter), 2. There is a minimum number of alternances (4 for this study) along the genes of S between those which have an ortholog on S 1 and those which have an ortholog on S 2 .

These conditions, inspired by the comparison between tetraodon and human genomes in [START_REF] Jaillon | Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF], were designed to retrieve genome segments whose gene content exhibits a clear signal for originating from a single genome segment pre-WGD. We considered here all comparisons between one of the three teleost genomes and each amniote. We obtained a list of amniote genome segments orthologous to two segments in two chromosomes of a fish, showing a signal for a "double synteny". The maximal set of genomic markers intersecting such a segment defines the corresponding ancestral synteny.

We also included ancestral syntenies from mammalian ancestral segments showing a DCS signal with fishes. These are sets of genomic markers intersecting maximal common intervals between two mammalian genomes whose evolutionary path goes through the boroeutherian ancestor, which are in addition included in a DCS in at least one of the two considered genomes. These segments are refined DCS, with a stronger conservation signal, as they are predicted to be ancestral by two methods, in the boreoeutherian ancestral genome and in the osteichtyes ancestral genome.

We obtain 2745 ancestral syntenies containing more than one genomic marker. These syntenies are then weighted according to the pattern of conservation they present in the phylogeny, using the formula described in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF] that accounts for a species tree and the branch lengths.

Assembling ancestral syntenies

The output of the phase described above is a set L of n genomic markers, and a family S = {S 1 , . . . , S m } of m subsets of L, where each subset is a set of genomic markers that are believed to be contiguous in the ancestral genome of interest. Following [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], we use the approach traditionally applied to physical mapping problems [START_REF] Alizadeh | Physical mapping of chromosomes using unique probes[END_REF]. It is based on the consecutive ones property (C1P) and PQ-trees. We encode S by an m × n 0/1 matrix M where row i represents S i as follows: M[i, j] = 1 if marker j belongs to S i and 0 otherwise. Ordering markers into CARs consists in finding a permutation of the columns of the matrix M, such that all 1's entries in each row are consecutive (also called a C1P ordering for M). Finding such an order of the columns of M is not always possible, in particular if there are false positives in S, that is groups of markers that were not contiguous in the ancestral genome. Moreover, if there exists a C1P ordering of the columns of M, there are often several possible orderings that make all 1's consecutive on each row, that represent different ancestral genome architectures. An equivalent representation of T (M) that highlights all ancestral genome architectures that correspond to C1P orderings for M: each row corresponds to a chromosomal segment represented by a child of the root, two glued blocks have to be adjacent in any ancestral genome architecture and sets blocks that float in the same box have to be consecutive in any genome architecture but their order is not constrained. Here for example we see three ancestral chromosomal segments and the second one contains markers 5 to 8, with only constraint that markers 6 and 7 are adjacent; hence, 5 6 7 8 is a possible order for this last segment, but not 5 6 8 7. All 13824 possible C1P orderings (possible ancestral orderings) are visible on this representation, that we use to present the amniote CARs in Figure 3.

If M is not C1P, then we know that some sets of markers in S are false positives and were not contiguous in the ancestral genome. Following [16,[START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], we clear ambiguities by computing a maximal subset of S that is C1P, using a branchand-bound algorithm described in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF] that finds an exact solution. Then, given this C1P subset of S, all C1P orderings can be represented in a compact way, using the PQ-tree of the resulting matrix M ′ , denoted T (M ′ ), that contains three kinds of nodes: leaves (labeled by L), P-nodes and Q-nodes. Computing T (M ′ ) can be done efficiently [START_REF] Mcconnell | A certifying algorithm for the consecutive-ones property[END_REF]. See Figure 2 for an illustration. T (M ′ ) encodes in a compact way all possible C1P orderings of the columns of M and then all genome architectures we can deduce from S: the root of T (M) is a P-node, the children of the root represent CARs, where Q-nodes describe fixed orderings, up to a reversal, while P-nodes, but the root, describe subsets of markers that have to be contiguous but where there is no information to fix a relative order (see Figure 2 for an illustration). Two markers that are consecutive children of a Q-node are said to define an adjacency.

Results and discussion

An ancestral amniote genome architecture

We applied the described method to propose a genomic architecture of the ancestral species of all amniotes, using the data presented in the previous section. Of the 2745 ancestral syntenies, 372 had to be removed during the optimization phase, in order to have some solutions to the C1P problem. This resulted in an ancestral amniote genome architecture composed of 79 CARs, 63 of them containing more than one genomic marker. In these 63 CARs, 983 of 1101 genomic markers are included in adjacencies, which indicates that there is little ambiguity due to P-nodes in the PQ-tree, although more than in the boreoeutherian CARS of [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF]. These CARs define 282 segments that are strictly colinear with segments of the chicken genome, and cover 75% of the chicken genome. Similarly, these CARs define 225 segments that are colinear with the human genome and cover 67% of this genome. Although these numbers are smaller than for reconstruction of mammalian ancestors described in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF], they are much larger than the amniote CARs inferred in [START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF]. The 51 CARs which span more than 1 Mb of the chicken genome are illustrated on Fig. 3, with their correspondence with the chicken chromosomes. 1. Recovered syntenic associations between chicken chromosomes, for three different methods. Numbers refer to pieces of chicken chromosomes. CARs which contain markers from only one chicken chromosome are not mentioned here.

Chicken syntenic associations. A "syntenic association" is the presence in a single CAR of genomic markers from two different chromosomes of an extant amniote. Here we may observe several chicken syntenic associations and compare them with other published methods [START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF][START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF]. These are, up to our knowledge, the only two methods that lead to the proposition of an architecture for the amniote genome. These propositions are very divergent: not only the number of chromosomes varies between 18 [START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF] and 26 [START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF], but the observed syntenic associations between chicken chromosomes are not always compatible. Of 13 syntenic associations found by Kohn et al. [START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF] and 6 found by Nakatani et al. [START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF], only two are common (Z-17 and 2-9). The reason is probably the absence of a formal framework, that we tend to fill here. We give a summary of the differences in Table 1, together with the syntenic associations we find in this study. We find one of the common syntenic associations (17-Z), plus two associations from Kohn et al. [START_REF] Kohn | Reconstruction of a 450-my-old ancestral vertebrate protokaryotype[END_REF] and none additional from Nakatani et al. [START_REF] Nakatani | Reconstruction of the vertebrate ancestral genome reveals dynamic genome reorganization in early vertebrates[END_REF].

Stability of the method and sensitivity to parameters

The advantage of using a general framework for ancestral genome reconstruction is the possibility, to a certain extent, of assessing the quality and robustness of the results. We have a simple support to every pair of adjacent markers in the CARs: the existence of an ancestral synteny that contains these markers. So every adjacency in CARs may be examined using the data, independently of the methodology.

The choices we have to make are the parameters used in the computations of ancestral syntenies. When two amniotes are compared, the ancestral synteny construction requires no parameter and its stability has been assessed in [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF].

When an amniote genome and a fish genome are compared, we rely on a novel method that is less tried and tested, that is the computation of DCS. The principle itself is well known and employed, but we are not aware of any formal study on a reliable implementation of this principle. The two parameters that are used to construct the DCS are inclusive: strengthening both parameters will improve the specificity. Any DCS which is found for a set of parameters will also be found by less stringent parameters. The question is which parameters are stringent enough to assure a good specificity. We think an amniote segment with at least 20 genes, covering 95% of the genes annotated in this segments is a sufficient proof for an double orthology signal, and it is confirmed by the comparison with the map of [START_REF] Jaillon | Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF] made from the same principle with visual expertise. The optimization step is also a source of possible instability of the method, especially as, from our experiments with the stringent criterion min prop = 95, at least 10% of the possible ancestral syntenies are false positive or result from convergent evolution, and have to be discarded during the optimization phase.

To assess these two possible sources of instability, we inferred DCS with the following values of the parameter min prop : 80%, 83%, 86%, 89%, 92% and 95%. To measure the stability, we considered the set of adjacencies (pairs of markers that are consecutive children of a Q-node) defined by each set of CARs. An adjacency is said to be conserved between two sets of CARs if it is found in both sets of CARs. An adjacency of a given set of CARs is said to be weakly conserved in another set of CARs if it is absent in the latter one but the two markers that define it belong to the same CAR. We show in Table 2 the characteristics of the set of CARs we computed.

positive ancestral syntenies, as the ratio of discarded syntenies drops from 29% to 14%. However, the results obtained are very stable in terms of adjacencies, as most adjacencies of a given set of CARs are consistent with the adjacencies in the previous set of CARs. The increase in the number of CARs is expected as less ancestral syntenies are obtained with more stringent parameters to compute DCS. This seems to indicate that, on this particular dataset, the optimization phase behaves very consistently with various sets of DCS and seems to conserve a subset of DCS that lead to very similar sets of CARs. Note also that the proportion of genomic markers that do not belong to adjacencies (i.e. are children of a P-node) is low, which indicates that the proposed ancestral genomic architectures contain very few segments where the order of the markers is not fixed.

Additional results regarding the support of each adjacencies by the different types of ancestral syntenies -common intervals between ingroups, DCS and mammalian syntenies included in DCS -are available on the companion website. They show that with these data, that contain only one species on one branch from the ancestor (birds and reptiles), DCS are fundamental to detect and support a significant number of adjacencies: 112 adjacencies are not supported by any ancestral synteny chicken-mammalian.

Conclusion

We extended a general framework for reconstructing ancestral genome architectures [START_REF] Chauve | A methodological framework for the reconstruction of contiguous regions of ancestral genomes and its application to mammalian genome[END_REF] in order to handle WGD events. We apply our method to reconstruct the architecture of the ancestral amniote genome. While we put a lot of attention to the specificity of the method, not to infer doubtful ancestral syntenies, sensitivity is not sufficient to provide the exact set of chromosomes of the amniote ancestor. The small number of genomes used, as well as the fact that the chicken is the only available genome among birds and reptiles and has several very small chromosomes, is also a reason for which there are certainly more CARs than ancestral chromosomes. The definitive amniote ancestral genome is still an open problem, but with this general, simple and formal method, some of its characteristics are accessible and valuable for further studies.

From a methodological point of view, several avenues can be explored. One of the issues is the fact that the C1P framework requires that each genomic marker appears exactly once in the ancestral genome architecture. This forced us to define our genomic markers from whole-genome alignments. But already at the level of the amniotes, these markers span around 30% of the chicken genome. Another approach could have been to use genes as genomic markers. However, in order to apply the C1P framework, this requires to compute the gene content of the ancestral amniote genome using the gene trees/species tree reconciliation approach. This problem is still a hard problem, that is very sensitive for example to errors in computing gene trees [START_REF] Hahn | Bias in phylogenetic tree reconciliation methods: implications for verte-brate genome evolution[END_REF]. The other main issue is the computation of the DCS. The results we obtain with several values of the parameter that defines DCS clearly show that many of the DCS we compute are probably not ancestral syntenies. They are very likely to intersect or contain genome segments that really originate from an ancestral amniote genome segment, but due to the lack of flexibility of the C1P framework, they are considered as false positive. Indeed, an amniote segment that only overlaps, even on a large part, a complete segment derived from the amniote common ancestor, will induce conflict with respect to the C1P framework. In the same time, without DCS, a significant number of adjacencies are not supported, which makes DCS instrumental in our method. This motivates the problem, open up to now, of the design and study of formal methods for the reliable and precise detection of DCS.

Fig. 2 .

 2 Fig. 2. (a) A matrix M with the consecutive ones property. (b) A PQ-tree T (M). (c)An equivalent representation of T (M) that highlights all ancestral genome architectures that correspond to C1P orderings for M: each row corresponds to a chromosomal segment represented by a child of the root, two glued blocks have to be adjacent in any ancestral genome architecture and sets blocks that float in the same box have to be consecutive in any genome architecture but their order is not constrained. Here for example we see three ancestral chromosomal segments and the second one contains markers 5 to 8, with only constraint that markers 6 and 7 are adjacent; hence, 5 6 7 8 is a possible order for this last segment, but not 5 6 8 7. All 13824 possible C1P orderings (possible ancestral orderings) are visible on this representation, that we use to present the amniote CARs in Figure3.

Fig. 3 .

 3 Fig. 3. The PQ-tree of amniote CARs, with their correspondence in the chicken genome. All ancestral architectures are represented, in the format described in Figure 2, with few chromosomic segments in which the order of the markers is not fixed.

We use release

51, which is based on the same Genome assemblies than release 50, with improved gene annotation, but does not contain Pecan alignments.

Acknowledgments

Aïda Ouangraoua is funded by the ANR BRASERO (ANR-06-BLANC-0045) and SFU. Eric Tannier is funded by the ANR (ANR-08-GENO-003-01 and NT05-3 45205) and by the CNRS. Cedric Chauve is funded by NSERC and SFU.


Abdelkader Akhenak

Mohammed Chadli
email: mohammed.chadli@u-picardie.fr

José Ragot

Didier Maquin

* Irsee

Design of a fuzzy observer for Takagi-Sugeno fuzzy

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

Design of observers for Takagi-Sugeno fuzzy models for Fault Detection and Isolation

A. Akhenak * M. Chadli * * J. Ragot * * * D. Maquin * * *

INTRODUCTION

A system for automatic fault detection and isolation (FDI) in a technical process can bring several benefits to the operator of the plant. Indeed, it provides him with information on the state of health of the supervised system at each moment. This allows the operator to follow the evolution of possible system degradations, and to predict when maintenance will be needed (see for example Patton et al. [1998], Kinnaert [1999]). A typical system for FDI is made of three parts: a residual generator, a residual evaluation system, a decision-making part.

Many standard observer-based techniques exist in the literature providing different solutions to both the theoretical and practical aspects of FDI problem for linear and nonlinear system (Frank [1996], Frank et al. [2000]). Many of these procedures are based on the design of an unknown input observer robust with respect to the disturbances (Staroswiecki and Varga [2001], Kinnaert [2003], Yan and Edwards [2007]). If the disturbances and modelling errors are not properly taken into account in the estimation process, it is then likely that any attempt in monitoring the system's health based on the observer leads to numerous false alarms.

Over the past decades, many researchers have paid attention to the problem of state estimation of dynamic linear systems subjected to both known and unknown inputs (Darouach et al. [1994], Wang et al. [1995], Sename [1997]). Many approaches have developed full and/or reduced order unknown input observers to estimate the state of linear time-invariant dynamical system driven by both known and unknown inputs (Xiong and Saif [2003], Edwards [2004]).

However, the real physical systems are often nonlinear. As it is delicate to synthesize an observer for an unspecified nonlinear system, it is preferable to represent this system with the Takagi-Sugeno (T-S) fuzzy model (Takagi and Sugeno [1985]). The idea of the T-S model approach is to apprehend the global behaviour of a system by a set of local models (linear or affine), each of them characterizing the behaviour of the system in a particular zone of operation. The local models are then aggregated by means of an interpolation mechanism. This approach includes the multiple model approach (Murray-Smith [1997], Chadli et al. [2003]) and Polytopic Linear Differential Inclusions (Boyd et al. [1994]) .

In this paper, for state and unknown input estimation, the suggested technique consists in associating to each local model a local unknown input high gain observer. The considered observer is then a convex interpolation of these local observers. This interpolation is obtained throughout the same activation functions as the T-S fuzzy model. Our contribution here lies in the design of the unknown input fuzzy observer by eliminating the unknown inputs from the dynamics of the state estimation error. The synthesis conditions of fuzzy observer are expressed in LMI terms. The designed observer is then applied for sensors and actuators FDI.

The paper is organized as follows. In Section 2, the general structure of the considered T-S fuzzy model is presented. In Section 3, the proposed structure of unknown input fuzzy observer is described and the main result is presented. The derived conditions ensuring the global asymptotic convergence of the estimation error are given as a set of LMI terms. The proposed T-S observer is used for detection and isolation of faults which can affect a nonlinear models in T-S representation. The validity of the proposed methodology is illustrated by estimating the yaw rate and faults of automatic steering vehicle.

Notation: Throughout the paper, the following useful notation is used: X T denotes the transpose of the matrix X, X>0m e a n st h a tX is a symmetric positive definite matrix, I N = {1, 2, ..., N } and . represents the Euclidean norm for vectors and the spectral norm for matrices.

TAKAGI-SUGENO FUZZY MODEL REPRESENTATION

Many physical systems are very complex in practice so that rigorous mathematical model can be very difficult to obtain, if not impossible. However, many of these systems can be expressed in some form of mathematical model locally or as an aggregation of a set of mathematical models. Here, using the T-S fuzzy dynamic model (Takagi and Sugeno [1985]), we consider a complex nonlinear system with unknown inputs. Then, the following T-S fuzzy model is adopted :

     ẋ = N i=1 µ i (ξ) A i x + B i u + R i ū + D i y = Cx (1)
with

N i=1 µ i (ξ)=1, 0 ≤ µ i (ξ) ≤ 1 ∀i ∈ I N
where x ∈ R n is the state vector, u ∈ R m the input vector, ū ∈ R q , q<n , contains the unknown inputs and y ∈ R p the measured outputs. Matrices A i ∈ R n×n and B i ∈ R n×m denote respectively the state matrix, the input matrix associated with the ith local model. Matrices R i ∈ R n×q are the distribution matrices of unknown inputs. D i ∈ R n is introduced to take into account the operating point of the system. At last, ξ is the so-called decision vector which may depend on some subset of the known inputs and/or measured variables to define the operating regimes.

Matrices A i , B i , R i and C can be obtained by using the direct linearization of an a priori nonlinear model around operating points, or alternatively by using an identification procedure (Gasso et al. [2001], Johansen and Babuska [2003], Angelov and Filev [2004], Murray-Smith [1997]).

In the following it is assumed that the vector ξ depends on measurable variables.

STRUCTURE OF THE UNKNOWN INPUT FUZZY OBSERVER

In this paper, we consider the state and unknown input estimation of a T-S fuzzy model perturbed by unknown inputs. The proposed unknown input fuzzy observer is based on a nonlinear combination of local unknown input observer. The proposed structure of the T-S fuzzy observer has the following form:

ẋ = N i=1 µ i (ξ) A i x + B i u + D i + G i (y -C x)+R i ûi (2a) ûi = γW i y -C x (2b) û = M i=1 µ i (ξ) ûi (2c)
T h ea i mo ft h ed e s i g ni st o determine gain matrices G i ∈ R n×p , W i ∈ R q×p and the positive scalar γ, that guarantee the asymptotic convergence of x towards x. Let us note that ûi can be considered as variables which compensate the errors due to the unknown inputs.

Assumption: In this paper, we consider that the inputs ū(t) are bounded: ū <ρ (3) where ρ is a positive scalar.

The set of residuals is defined as follows:

r = y -ŷ (4a) = Ce (4b
) where e represents the state estimation error: e = x -x (5) The dynamic of the state estimation error is given as follows:

ė = N i=1 µ i (ξ) (A i -G i C)e + R i ū -R i ûi (6)
Replacing ûi by its expression (2b) we get

ė = N i=1 µ i (ξ) (A i -G i C)e + R i ū -γR i W i Ce (7)
The following result proposes a solution to design the gain parameters G i ∈ R n×p , W i ∈ R q×p and the positive scalar γ of the proposed observer (2).

Theorem 1: the state estimation of the unknown input fuzzy observer (2) can asymptotically estimate with any desired degree of accuracy ε>0, the state of the T-S fuzzy model (1), if there exist a symmetric positive definite matrices P and Q and the gain matrices G i and W i which satisfies the following constraints ∀ i ∈ I N :

(A i -G i C) T P + P (A i -G i C) < -Q (8a) W i C = R T i P (8b)
The fuzzy observer (2) is then completely defined by choosing:

γ ≥ 1 2 λ min (P -1 Q) λ min (P ) ε 2 -1 ρ 2 (9)
and the input estimation is given by

û = γ M i=1 µ i (ξ) W i y -C x ( 10
)
where ρ is defined in (3).

Proof: Consider the Lyapunov function V (e)=e T Pe (11) Its derivative with respect to time is

V = N i=1 µ i (ξ) e T (A i -G i C) T P + P (A i -G i C) e+ 2e T PR i ū -2γe T PR i W i Ce ≤ N i=1 µ i (ξ) e T (A i -G i C) T P + P (A i -G i C) e+ 2ρ (PR i ) T e -2γe T PR i R T i Pe
For any positive scalar β, we have the following inequality:

2ρ (PR i ) T e(t) ≤β -1 ρ 2 (PR i ) T e(t) 2 + β (12) Then V ≤ N i=1 µ i (ξ) e T (A i -G i C) T P + P (A i -G i C) e+ β -[2γ -β -1 ρ 2 ] (PR i ) T e 2 Choosing γ ≥ 1 2 β -1 ρ 2 (13) V ≤ N i=1 µ i (ξ) e T (A i -G i C) T P + P (A i -G i C) e + β (14)
Using the inequality (8a), the derivative of the Lyapunov function becomes as follows V (e) ≤-e T Qe + β (15) with α = 1 2 λ min (P -1 Q) > 0, we can easily deduce V (e) ≤-2 αV(e)+β (16) Tacking account the expression of V (e) (11) and the fact that α>0, we get when t →∞the following inequality e ≤ 1 λ min (P )

β 2α (17)
To guarantee for any desired ǫ>0that e <ǫ, it suffices to choose

γ ≥ 1 2 λ min (P -1 Q) λ min (P ) ε 2 -1 ρ 2 (18)
which end the proof.

Remark 1: The inequalities (8a) are nonlinear in P and G i . To linearize these inequalities, the following change of variables is used

K i = PG i (19)
We obtain a linear matrix inequalities in P and K i that can be easily solved by the means of LMI tools:

A T i P + PA i -C T K T i -K i C<-Q (20a) W i C = R T i P (
20b) Finally, the matrix gains G i are computed as follows

G i = P -1 K i (21) Remark 2:
In case of presence of noise measurement on the system output, choosing a large gain γ leads to the amplification of noise and thus a bad estimation state and unknown inputs. [2008]). Here, we have chosen to consider the coupling model of longitudinal and lateral motions of a vehicle. This model, already used in Zhang and Xu [2002], is strongly nonlinear and is given by the following equations:

u = vr -fg + (fk 1 -k 2 ) M u 2 + c f v + ar Mu δ + T M (22a) v = -ur - c f + cr Mu v + bcr -ac f Mu r + c f δ + Tδ M (22b) ṙ = bcr -ac f Izu v - b 2 cr + a 2 c f Izu r + aT δ + ac f δ Iz (22c)
where, u, v and r are the longitudinal velocity, the lateral velocity and the yaw rate, respectively, δ is the steering angle, T is the traction and/or braking force. The nonlinear vehicle dynamics can be written as follows:

ẋ = F (x, w) (23a) y = Cx (23b)
where C = 100 001 , F is a nonlinear function of the state vector x =[uvr ], w gathers the two inputs δ and T and y(t) gathers the two inputs y 1 = u and y 2 = r.A si ti s delicate to synthesize an observer for a nonlinear system, we preferred to represent this system with a T-S fuzzy model. Then, we propose to linearize the nonlinear model ( 23) around some operating points [x (i) w (i) ]. Next, we integrate the set of the linear models in a T-S fuzzy model. The proposed T-S model is described as follows (Akhenak et al. [2007]):

ẋ = N i=1 µ i (y 1 )(A i x + B i w + D i ) (24a) A i = ∂F ∂x x=x (i) w=w (i) B i = ∂F ∂w x=x (i) w=w (i) (24b)
D i = F (x (i) ,w (i) ) -A i x (i) -B i w (i) (24c) with N i=1 µ i (u)=1 and µ i (u) ≥ 0 ∀i ∈ I M A i = A 11i A 12i A 13i A 21i A 22i A 23i A 31i A 32i A 33i A 11i = 2(fk 1 -k 2 ) u i M - c f (v i + ar i ) Mu 2 i δ i A 12i = r i + c f δ i Mu i A 13i = v i + ac f δ i Mu i A 21i = -r i + (c f + c r ) Mu 2 i v i - (bc r -ac f ) Mu 2 i r i A 22i = - (c f + c r ) Mu i A 23i = -u i + (bc r -ac f ) Mu i A 31i = - (bc r -ac f ) I z u 2 i v i + b 2 c r + a 2 c f I z u 2 i r i A 32i = (bc r -ac f ) I z u i A 33i = - b 2 c r + a 2 c f I z u i B i =        c f v i + ar i Mu i 1 M c f + T i M δ i M aT i + ac f I z aδ i I z        D i = F (x i ,δ i ,T i ) -A i x i -B i δ i T i
where x i =( u i v i r i ) T ,w i t hi = {1, 2, 3}, x i , δ i and T i are the operating points. Three local models are used to approximate the nonlinear model ( 22). Using a quadratic criterion of the error between the state variables of the nonlinear model ( 22) and state variables of the T-S fuzzy model ( 24), we obtain the following operating points:

Three local models were chosen for this application. This number gives a good compromise between the quality of the obtained model and its complexity. The activation functions given in figure 2 depend only on the longitudinal velocity u(t). The numerical values of the different matrices A i , B i , D i are: The membership functions µ 1 (u), µ 2 (u)a n dµ 3 (u)a r e given in figure 2.

A 1 =    0
The system ( 22) is simulated using the steering angle δ and attraction force T depicted in figure 1. Figures 2 and3show that the state variables of the nonlinear model ( 22) and its approximate by T-S fuzzy model ( 24) are superimposed. In order to check the good accuracy of the T-S fuzzy model ( 24), its behaviour and that of the nonlinear model ( 22) have been simulated. Figure 1 shows the time evolution of the two inputs δ, the steering angle and T the tractionbraking force.

Figures 2 and3 show the superposition of the state vector of the nonlinear model ( 22) and their approximation by the T-S fuzzy model ( 24). Clearly the T-S fuzzy model is able to represent the nonlinear behaviour of the considered system.

The following section exploits the obtained T-S model and the observer designed above to propose a diagnosis method for the steering vehicle model ( 22).

FAULT DETECTION AND ISOLATION

The objective of this part is to generate residuals that reflect the faults acting on the system (24). An ideal residual signal should remain zero in the fault-free case and non-zero when fault occurs. Once a fault has been detected, it must be estimated. The fault estimation will specify the type of fault, its duration, its amplitude and eventually its probable evolution. In the literature, there are several fault detection techniques. They are generally based on the change detection of the average and the variance. In this FDI study, we will not deal with the detection thresholds of residuals. We will confine ourselves only to the detection and localization of sensor (subsection A) and actuator faults (subsection B) taking into account the uncertainties modeling.

Sensor fault detection and isolation

In order to identify the sensor fault, we consider that the actuators are faultless ( w = 0) while the output vector y is corrupted by the sensor fault ∆y. Then the system (24) becomes:

     ẋ = M i=1 µ i (y 1 ) A i x + B i w + D i y = Cx +∆y + ν (25)
Firstly we can easily checked that the following necessary conditions are satisfied: ∀i ∈{ 1, 2, 3}, j ∈{ 1, 2} rank(A i ,C(j, :)) = 3

Three T-S observers are designed, the first is based on the longitudinal velocity y 1 = u, the second is based on the yaw rate y 2 = r and the last is based on the two outputs u and r. The sensor fault detection and localization is based on the analysis of the residuals r y ik = y i -ŷi,k , k ∈{ 1, 2, 3}, i ∈{ 1, 2} generated by three observers (figure 4) which depend on two inputs δ and T applied to the system (22). The longitudinal velocity observer 1 and the yaw rate observer 2 use respectively only one output u and r.T h e global observer 3 uses two outputs u and r.

Figures 5 shows the additive signal that represents sensor failure, the fault has been added to sensor 2 output y 2 = r between 5 and 10s.

Fault detection using unknown input fuzzy observer

The problem is to detect the occurrence of two fault signals m 1 and m 2 such as the wind force, which disturb the motion of the vehicle. To this end, we suppose the existence of a fault signal on the system (nonlinear model ( 22)). For that, we use the unknown input fuzzy observer (2) developed previously in this paper. In this case, we consider the nonlinear model with defect describes as follows: ẋ = F (x, δ, T )+Rm ( 26)

with m = m 1 m 2 ,R = 10 0 -1 10
The considered unknown input fuzzy observer that estimates the state and the defect (regarded as an unknown input) of the nonlinear model ( 22) is described by:

ẋ = 2 i=1 µ i (u) A i x + B i w + D i + G i (y -C x)+R mi (27a) mi = γW i y -C x (27b) m = 2 i=1 µ i (u)m i (27c)
For simulation the following initial conditions are considered: The synthesis conditions lead to the resolution of an LMI problem. Moreover, the estimation of unknown inputs of the system is considered. The proposed observer is then used for state estimation and for detection and isolation of faults which can affect nonlinear models. The effectiveness of the proposed methodology is illustrated by estimating the yaw rate and faults of automatic steering vehicle. The considered structure of unknown input T-S observer can be also benefic for fault detection and isolation of a fault actuator.

(u 0 v 0 r 0 ) = ( 1 500 )a n d( x1 (0) x2 (0) x3 (0)) = (16 1 1). The gain matrices G 1 , G 2 , G 3 , W 1 ,

  4. APPLICATION TO AUTOMATIC STEERING OF VEHICLE4.1 Representation of the vehicle model by a T-S fuzzy model Different models related to automatic steering of vehicle have been studied in the literature (see for example Zhang and Xu [2002], El Hajjaji and Bentalba [2003], Moriwaki [2005], Chadli et al. [2008], Oudghiri et al.

Fig. 3 .

 3 Fig. 1. δ: steering angle T : attraction force

Fig. 4 .

 4 Fig. 4. Block diagram of the banc observer-based FDI

Fig. 5 .Fig. 7 .

 57 Fig. 5. Sensor failure ∆y 2

  Figures 9 and 10 represent the comparison between the fault signals m 1 and m 2 affecting the nonlinear model (22) and their estimates by the unknown input fuzzy observer (27). These figures clearly show the occurrence of the fault signals m 1 and m 2 . Figures 11 and 12 present the comparison between the output of the nonlinear model with the faults m 1 and m 2 (26) and their estimates with the designed unknown input observer. The output variables and their estimates are superimposed except in the vicinity of the origin (choice of the initial conditions).

Fig. 12 .

 12 Fig. 9. Fault m 1 and its estimate m1 by (27)

  Table 1 lists the parameters of the above vehicle model.

		Parameters of the vehicle system
	M	Mass of the full vehicle	1480 kg
	Iz	Moment of inertia	2350 kg.m 2
	g	Acceleration of gravity force	9.81 m/s 2
	f	Rotating friction coefficient	0.02
	a	Distance from front axle to CG 1	1.05 m
	b	Distance from rear axle to CG	1.63 m
	c f	Cornering stiffness of front tyres	135000 N/rad
	cr	Cornering stiffness of rear tyres	95000 N/rad
	k 1	Lift parameter from aerodynamics	0.005 Ns 2 /m 2
	k 2	Drag parameter from aerodynamics	0.41 Ns 2 /m 2

  W 2 and W 3 are obtained by solving the constraints (20). We get

		9.22 -3.88		10.78 -4.10
	G 1 =	0.45 -1.02	G 2 =	6.64 0.55
		22.51 -11.92		27.38 -16.19
			8.49 -4.83
		G 3 =	4.40 1.36	γ =78.12
			20.27 -17.07

A. Akhenak, M. Chadli, D. Maquin, and J. Ragot. Design of a sliding mode fuzzy observer for uncertain takagisugeno: application to automatic steering of vehicle.

International Journal of Vehicle Autonomous Systems, 5(3/4):288 -305, 2007. P.P. Angelov and D.P. Filev. An approach to online identification of Takagi-Sugeno fuzzy models. IEEE Transactions on Systems, Man and Cybernetics, Part B, 34(1):484 -498, 2004. S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in System and Control Theory. Philadelphia: SIAM, 1994. M. Chadli, D. Maquin, and J. Ragot. Multiple observer for discrete-time multiple model. 5th IFAC Symposium on Fault Detection Supervision and Safety for Technical Processes, Safeprocess, Washington, USA, June 9-11, 2003. M. Chadli, A. Elhajjaji, and M. Oudghiri. Robust output fuzzy control for vehicle lateral dynamic stability improvement. International Journal of Modelling, Identification and Control, 3(3):247 -257, 2008. M. Darouach, M. Zasadzinski, and S. J. SXu. Full-order observers for linear systems with unknown inputs. IEEE Transactions on Automatic Control, 39(3):606 -609, 1994. C. Edwards. A comparison of sliding mode and unknown input observers for fault reconstruction. IEEE Conference on Decision and Control, 5:5279 -5284, 2004. A. El Hajjaji and S. Bentalba. Fuzzy path tracking control for automatic steering of vehicles. Robotics and Autonomous Systems, 43(3):203 -213, 2003. P. M. Frank. Analytical and qualitative model-based fault diagnosis -a survey and some new results. European Journal of Control, 2:6 -28, 1996. P. M. Frank, S. X. Ding, and T. Marcu. Model-based fault diagnosis in technical processes. Trans. of the Inst. of Measurement and Control, 22(1):57 -101, 2000. K. Gasso, G. Mourot, and J. Ragot. Structure identification in multiple model representation: elimination and merging of local models. IEEE Conference on Decision and Control, 3:2992 -2997, 2001. T. A. Johansen and R. Babuska. Multiobjective identification of Takagi-Sugeno fuzzy models. IEEE Transactions on Fuzzy Systems, 11(6):847 -860, 2003. M. Kinnaert. Robust fault detection based on observers for bilinear systems. Automatica, 35:1829 -1842, 1999. M. Kinnaert. Fault diagnosis based on analytical models for linear and nonlinear systems. Proceedings of IFAC Safeprocess, Washington, USA, pages 37 -50, 2003. K. Moriwaki. Autonomous steering control for electric vehicles using nonlinear state feedback H1 ∞ control. Nonlinear Analysis, 63(5-7):e2257-e2268, 2005. R. Murray-Smith. Multiple model approaches to modelling and control. Taylor and Francis, 1997. M. Oudghiri, M. Chadli, and A. El Hajjaji. Robust observer-based fault tolerant control for vehicle lateral dynamics. International Journal of Vehicle Design, 48(3/4), 2008. R. J. Patton, J. Chen, and C. J. Lopez-Toribio. Fuzzy observer for nonlinear dynamic systems fault diagnosis. IEEE Conference on Decision and Control, 1:84 -89, 1998. O. Sename. Unknown input robust observer for time delay system. IEEE Conference on Decision and Control,2 : 1629 -1630, 1997. M. Staroswiecki and G.C. Varga. Analytical redundancy relations for fault detection and isolation in algebraic dynamic systems. Automatica, 37:687 -699, 2001. J. Takagi and M. Sugeno. Fuzzy identification of systems and its application to modelling and control. IEEE Transactions on Systems Man and Cybernetics, 15:116 -132, 1985. S. H. Wang, E. J. Davison, and P. Dorato. Observing the states of systems with unmeasurable disturbances. IEEE Transactions on Automatic Control, 20:716 -717, 1995. Y. Xiong and M. Saif. Unknown disturbance inputs estimation based on a state functional observer design. Automatica, 39:1389 -1398, 2003. X.G. Yan and C. Edwards. Nonlinear robust fault reconstruction and estimation using a sliding mode observer. Automatica, 43(9):3 -40, 2007. J.R. Zhang and S.J. Xu. Nonlinear observer design for automatic steering of vehicles. IEEE 5th International Conference on Intelligent Transportation Systems, pages 179 -182, 2002.


Boleslaw Mielniczuk
email: boleslaw.mielniczuk@univ-montp2.fr

Tomasz Hueckel
email: hueckel@duke.edu

Moulay Said

El Youssoufi

Micro-scale Testing of Capillary Bridge Evolution due to Evaporation

Capillary bridge evolution between two fixed glass spheres during its natural convective evaporation is examined experimentally. For comparison extension tests were also carried out. The calibrated balance recording and digital image processing allow monitoring of a number of key variables of the process: the resultant capillary force, the water mass loss, radii of the bridge curvature. On that basis evaporating surface area, suction and surface tension force, interparticle force, axial stress vs (relative) volumetric mass loss are calculated. Testing shows a gradual decrease of suction within bridges down to zero and into a positive pressure range before a two step failure including a formation of a water thread according to a traditional Rayleigh instability pattern followed by a simultaneous rupture at two points of the lowest (negative) total (Gauss) curvatures of the bridge surface.

Introduction

The liquid loss from a porous (granular) material (i.e. soil) during evaporation plays a key role in the onset of drying cracks. Studies of desiccation shrinkage undertaken i.e. by Hu et al. [1], Kowalski and Mielniczuk [2], Peron at al. [3] revealed that the shrinkage in drying soil develops when soil is still saturated, while cracks appear shortly after the air entry, and formation and rupture of liquid bridges. On the other hand, the apparent cohesion of unsaturated soils is attributed to the generation of suction and capillary forces. Suction hardening and softening (wetting collapse) of granular materials are related to understanding of the evolution and rupture of liquid bridges.

To identify the mechanisms and critical variables that control evolution of the water body (liquid bonds) during evaporation we focus on the evolution and rupture of elementary liquid bridges between two, three and more particles. The important step of these studies is the analysis of the changes of the liquid bridge geometry and interparticle (capillary) forces during desiccation of the smallest clusters of grains at the local scale. In this paper, capillary bridges are addressed between two grains only. Such systems are not only the most elementary ones, but also, it has been seen that more articulated systems break down eventually to a number of two-grain arrangements.

Most of the past studies were conducted in static conditions or during extension of the bridge with a fixed liquid volume [4][5][6][7][8]. Instead, in the present experiments, the driving force of the investigated reponse is evaporation in a constant relative humidity environment, and at a constant predetermined separation of grains.

Materials and Methods

Liquid bridge is produced by placing a 4 µl of distilled water between two glass spheres (radius R=4mm, grade 25), with the use of microsyringe. Next, the liquid bridges are exposed to evaporation at constant ambient conditions, at constant sphere separations D from 0.1 to 2.0 mm, or are stretched at a constant speed of 16 µm/s, at several constant water volumes, from 0.2 to 10 µl. The lower sphere is placed on a precision balance (with precision of 10 -4 g). The upper sphere, attached to the micrometer, is suspended from the support, without any contact with the balance (Fig. 1a). The attractive force between spheres is calculated as F C =(m dry -m exp )•g, where m dry is the mass of dry lower sphere, m exp is a mass supported by the balance including the lower sphere and the liquid bridge during drying, g is the gravity. The evaporating mass of water is obtained in the same conditions, but also with the upper sphere supported by the balance. The changes in the bridge geometry were recorded by a photo camera with macro lens. The moment of the bridge rupture was recorded with a high-speed digital camera, with 27000 frames per second. Images were processed with Fiji and Mathlab software, to determine gorge radius r g and external radius r ext (see Fig. 1b). It was hence possible to calculate via Laplace law suction and surface tension force components of the capillary force, and the average axial stress at the bridge neck.

Results

Variations of bridge volume (evaporation) or bridge length (extension) are the main factors affecting the radii of the meniscus surface. During evaporation the profile of drying bridge evolves with the receding lateral meniscus surface of revolution, following the classical sequence of Plateau, from the surface of the negative unduloid, through catenoid to negative nodoid (Fig. 2) [9]. Measured value of r g decreases smoothly as a convex parabola, and external radius r ext decreases proportionally with the decreasing liquid volume. The r g for the bridge in extension decreases as a concave parabola with the increasing separation, with almost constant r ext (see [10] for details). The duration of the process till the bridge rupture varies from 9 to 67 min for evaporation and up to 150 seconds for extension.

Most notable is a relatively slow change of the bridge shape from a catenoid into a nodoid, which then developing into a very thin water wire. The wire ruptures simultaneously near both ends within 0.3-0.9 msec and then transitions into a small droplet of water, which next free-falls in the gravity field. Eventually the mass of water is separated into two drops attached to each of the spheres (Fig. 2).

Fig. 2 Evolution of evaporating capillary bridge between two spheres at separation of 1.3mm with the initial total water volume 4 µl, at decreasing water volume (in µl), detail of the rupture process (fast camera, time in milliseconds)

The bridge evolution due to evaporation at D= 1.3 mm is shown in Fig. 2. The bridge ruptures at the water volume V= 2.28 µl., with the water wire 90 µm in diameter, and 130 µm in length (see [11] for details).

Based on the measured evolution of the meniscus radii it is possible to estimate capillary forces, as a sum of a vertical resultant of the surface tension acting at the perimeter of the neck and the Laplace pressure ∆p (difference between ambient and internal bridge pressure) expressed as ∆p=γ(r g -1 -r ext -1 ) acting at the crosssection of the neck, F C =2πr g γ-πr g 2 ∆p [4,6,8], where γ= 0.0496 N/mis the experimentally determined surface tension of the actual water. At our calculations hydrostatic pressure of water bridge acting at the neck level is also considered. The surface tension force is positive and decays gradually following the neck perimeter length with volume loss or increasing separation. Its initial value is the highest at small separations and large volumes, as shown on Fig. 3a.

Laplace pressure arises as a consequence of the total bridge curvature, which can be either positive or negative, depending on gorge (neck) and external radii (Fig. 2). It is initially negative (suction) at small separations, but it may be positive at large separations (Fig. 3b). Such cases are described by Kralchevski and Nagayama [12]. Most notable is the evolution of suction during both evaporation and extension showing a consistent decrease until zero, and becoming eventually a positive pressure. Hence, pressure at bridge rupture is positive in all cases (except for the largest volume extension. The resultant suction force however remains attractive (due to surface area factor), for small separations, while reaching negligible positive values (repulsion) near rupture (see Fig. 3c). At large distances the resultant force is consistently repulsive.

The total capillary force (Fig. 4a) decreases consistently in both processes. For evaporation at large separations the force decreases almost linearly with an early sudden jump to zero at the bridge rupture. At separations D smaller than 0.7 mm there is an initial increase of F C up to a maximum, followed by a decay to zero. At separation D=2 mm the maximum (initial) force value F C = 0.14•10 -3 N (0.14 mN), with force jump of 0.09 mN at 23% of water loss. At D=0.1 mm the maximum force F C = 0.85 mN at 33% of water loss, and at rupture the force jump is 0.02 mN, at 96% of water loss.

Also during extension the total capillary force F C decreases to zero, with a slight initial increase at higher volumes.

Based on on the capillary force it was possible to visualize stress evolution for both tests, as presented in Fig. 4b. The values of tensile stress in water in both processes are of the same order of magnitude, with rupture occurring between 200 and 600 Pa, depending on separation D and total liquid volume V.

Conclusions

The behavior of an evaporating liquid bridge resembles tensile failure of a ductile steel rod, while no external tensile force is applied to the system. The process is controlled by a decreasing volume of the liquid bridge (at constant separation). The moment of rupture is dependent on separation of spheres. For larger separations the bridge rupture alters a smooth desiccation process early, whereas for very small separations, the rupture blends with the vanishing of all the liquid mass.

Notable findings concern the evolution of liquid suction, which monotonically decreases during evaporation (and extension) to zero, and becomes a positive pressure prior to rupture. The total capillary force between the grains remains attractive due to the surface tension effect. However, it also has a decreasing trenda and approaches very small values near rupture.

The intergranular bridge rupture is a two-stage process. The first stage includes a formation of a water thread in the central part of the bridge according to a traditional Plateau-Rayleigh instability pattern, followed by a rupture at simultaneously two points of the lowest (negative) total (Gauss) curvatures of the bridge surface. The later subject is pursued in Mielniczuk et al. [13]. The rupture of the bridge is critically important as it marks the end of suction controlled strengthening of unsaturated granular media. As seen it may occur very early in the process of evaporation in a loose medium with large grain separation. It also determines how much of water is left out after the process of bridge formation and rupture. Notably, the stage where two-grain bridges dry in isolation is the last stage of earlier processes when suction growth plays a crucial role in the removal of water from the porous medium (Hu et al.,[1]).

Fig. 1

 1 Fig. 1 a) Apparatus configuration to measure capillary force, b) Capillary bridge parameters

Fig. 3

 3 Fig. 3 a) Surface tension force for evaporation and extension tests, b) Laplace pressure changes for both experiments, c) Resultant suction force for evaporation

Fig. 4

 4 Fig. 4 a) Capillary force evolution during evaporation and extension tests, b) Axial stress vs decreasing volume for evaporation, vs increasing separation for extension tests The tensile stress at rupture at D= 1.3 mm (V= 2.28 µl) caused by evaporation is 220 Pa. For an ad hoc extension test performed at this particular liquid volume, the bridge ruptures at D= 1.76 mm, with tensile stress of 243 Pa. Notably, the axial stress at rupture is different for every test.

Acknowledgements. The authors acknowledge the partial support of US NSF grant #0700294 of CMMI Division, Geomechanics and Geomaterials and Program (T. Hueckel) and support of CNRS and MIST laboratory (B. Mielniczuk, M.S. El Youssoufi).