Dorothée Huchon

Frédéric Delsuc

François M Catzeflis

Emmanuel J P Douzery
email: douzery@isem.univ-montp2.fr

Armadillos exhibit less genetic

Keywords:

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

Armadillos exhibit less genetic polymorphism in North

America than in South America: nuclear and mitochondrial data confirm a founder effect in Dasypus novemcinctus (Xenarthra).

Dorothée Huchon, Frédéric Delsuc, François M. Catzeflis, Emmanuel J.P.

Douzery

To cite this version:

Dorothée Huchon, Frédéric Delsuc, François M. Catzeflis, Emmanuel J.P. Douzery. Armadillos exhibit less genetic polymorphism in North America than in South America: nuclear and mitochondrial data confirm a founder effect in Dasypus novemcinctus (Xenarthra).. Molecular Ecology, 1999Ecology, , 8 (10), pp.1743Ecology, -8. 10.1046Ecology, /j.1365Ecology, -294x.1999.00768..00768.x. halsde-00192969

INTRODUCTION

The nine-banded armadillo (Dasypus novemcinctus Linnaeus, 1758; Dasypodidae)

shows the largest distribution of all Xenarthra. This New World mammal ranges from Peru and Argentina to south-central United States [START_REF] Wetzel | The subgenera and species of long-nosed armadillos[END_REF]. Its expansion into the United States is unique among placentals in that it has been occurring since the mid-19th Century at a mean rate of 10 km a year (review in [START_REF] Taulman | Recent range expansion and distributional limits of the nine-banded armadillo (Dasypus novemcinctus) in the United States[END_REF]. Expanding from Texas, the territory of Dasypus novemcinctus now encompasses Kansas, Missouri, Illinois, Tennessee, South Carolina, and all the Gulf states. This fast migration may result from a low predation on adults, a lack of natural competitors, a weak homing ability (although it is rather sedentary with small home ranges), and human-induced translocations [START_REF] Mcbee | Dasypus novemcinctus[END_REF][START_REF] Taulman | Recent range expansion and distributional limits of the nine-banded armadillo (Dasypus novemcinctus) in the United States[END_REF][START_REF] Prodöhl | Genetic maternity and paternity in a local population of armadillos assessed by microsatellite DNA markers and field data[END_REF]. Dasypus novemcinctus is also characterized by an unusual reproductive mode, monozygotic polyembryony [START_REF] Mcbee | Dasypus novemcinctus[END_REF], wherein four genetically identical offspring are produced from a single egg [START_REF] Prodöhl | Molecular documentation of polyembryony and the micro-spatial dispersion of clonal sibships in the nine-banded armadillo, Dasypus novemcinctus[END_REF].

United States armadillos are almost completely monomorphic in enzymatic studies [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF][START_REF] Moncrief | Absence of genic variation in a natural population of nine-banded armadillos, Dasypus novemcinctus (Dasypodidae)[END_REF], and to display moderate microsatellite allelic variability [START_REF] Prodöhl | Molecular documentation of polyembryony and the micro-spatial dispersion of clonal sibships in the nine-banded armadillo, Dasypus novemcinctus[END_REF]. Preliminary results on one microsatellite locus also suggest that North American individuals may exhibit less variability than those from one Brazilian population [START_REF] Loughry | Polyembryony in armadillos[END_REF] Taulman & Robbins 1996), and/or by inbreeding possible reinforced by monozygotic polyembryony [START_REF] Moncrief | Absence of genic variation in a natural population of nine-banded armadillos, Dasypus novemcinctus (Dasypodidae)[END_REF]. Evaluation of these two hypotheses requires: (i) sampling of reference armadillos from South America, i.e. individuals from their primary habitat; (ii) estimation of the polymorphism of North and South American armadillos. The latter point involves screening by the biochemical systems of [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF] and [START_REF] Moncrief | Absence of genic variation in a natural population of nine-banded armadillos, Dasypus novemcinctus (Dasypodidae)[END_REF], and a confirmation using complementary (i.e.

nucleotide) markers, unlinked to nuclear protein loci.

Mitochondrial DNA (mtDNA) is a suitable tool for phylogeographic studies. The hypervariable 5'-(left)-end of its control region is of particular interest in revealing nucleotide variability (e.g., Wood & Phua 1996), and in evaluating the degree of polymorphism within populations [START_REF] Gonzalez | Conservation genetics of the endangered Pampas deer (Ozotoceros bezoarticus)[END_REF]. Moreover, the control region left domain of Dasypus novemcinctus contains tandemly repeated 80-nt blocks [START_REF] Arnason | Phylogenetic analyses of mitochondrial DNA suggest a sister group relationship between Xenarthra (Edentata) and Ferungulates[END_REF], which may exhibit a high within-population nucleotide variability (Wood & Phua 1996). According to the previous hypotheses, a low level of nuclear and mitochondrial polymorphism is expected over the whole geographic range of armadillos if polyembryony has an obvious impact on genetic diversity. To the contrary, if armadillos experienced persistent bottleneck(s) during migration to North America, higher polymorphism is predicted within their original range (South America).

MATERIALS AND METHODS

Protein data.

Horizontal starch gel electrophoresis was performed on blood extracts collected from 55 French Guiana individuals. Electrophoreses and protein staining were conducted according to [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF], with slight modifications. Only seven blood allozyme systems out of the 14 examined by [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF] Mitochondrial DNA control region data.

DNA was extracted from biopsies of 20 nine-banded armadillos collected in Mississippi, Texas, Louisiana, and French Guiana, and three great long-nosed armadillos (Dasypus kappleri Krauss, 1862) in French Guiana (Table 1). The 5'-(left)end (i.e. cytochrome b side) of the mitochondrial control region (581 base pairs) was amplified via polymerase chain reaction using primers L 0 [L15445] [START_REF] Douzery | The mitochondrial control region of Cervidae: evolutionary patterns and phylogenetic content[END_REF]) and E 3 [H15978] = 5'-ATGACCCTGAAGAAASAACCAG-3' (numbers defined by [START_REF] Arnason | Phylogenetic analyses of mitochondrial DNA suggest a sister group relationship between Xenarthra (Edentata) and Ferungulates[END_REF] for Dasypus novemcinctus mitochondrial genome). Amplification products were excised from agarose gels, and directly sequenced on both strands. L 0 and E 3 primed the sequencing reactions conducted with [α 33 P]ddNTP and the Thermo Sequenase radiolabeled terminator cycle sequencing kit (Amersham; Cleveland, OH, USA).

Sequences were manually aligned with the ED editor of the MUST package [START_REF] Philippe | MUST: a computer package of management utilities for sequences and trees[END_REF]. Haplotype and nucleotide diversities and their standard errors (SE)

were calculated according to [START_REF] Nei | Molecular evolutionary genetics[END_REF]. Phylogeographic relationships between ninebanded armadillo haplotypes were described using a minimum spanning network.

Phylogenetic relationships were reconstructed using maximum parsimony (MP) and maximum likelihood (ML), with respectively PAUP 3.1.1 [START_REF] Swofford | PAUP: Phylogenetic analysis using parsimony, version 3.1.1[END_REF]) and PUZZLE 4.0 [START_REF] Strimmer | Quartet puzzling: a quartet maximum-likelihood method for reconstructing tree topologies[END_REF]. Bootstrap and reliability estimates indicated robustness of the nodes, and the congeneric great long-nosed armadillo rooted nine-banded armadillo trees.

RESULTS

Protein polymorphism.

Monomorphism of 265 North American armadillos have been revealed by [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF] Relationships between control region sequences.

A minimum spanning network reveals a clear distinction between French Guiana and United States haplotypes (Fig. 1). Twenty substitutions separate the South American haplotype G from the North American haplotype Y, five of which are transversions.

French Guiana haplotypes cluster into two groups: A, C, D, H and I, vs. B, E, F, G and J (Fig. 1). The clear separation between North and South American populations is robustly confirmed by the ML reliability, MP bootstrap, and MP branch support (Fig. 2).

DISCUSSION

Founder effect and colonization of United States.

Nuclear and mitochondrial markers evidence contrasted polymorphism levels between South and North American nine-banded armadillos, with a severe loss of variability amongst individuals from Texas to Mississippi (including Louisiana). The latter observation confirms that of [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF]. Preliminary microsatellite studies also suggest a moderate to reduced variability in United States relative to one Brazilian population [START_REF] Loughry | Polyembryony in armadillos[END_REF]. Conversely, for a maximum distance range of 32 km, protein (see results) and mtDNA (comparative data not shown) polymorphisms in French Guiana fall within the range of usual placental values, suggesting that the polyembryony has no obvious consequence on the genetic diversity of Dasypus novemcinctus. Indeed, microsatellite surveys suggested that closely related adult armadillos do not often interact [START_REF] Prodöhl | Molecular documentation of polyembryony and the micro-spatial dispersion of clonal sibships in the nine-banded armadillo, Dasypus novemcinctus[END_REF][START_REF] Loughry | Polyembryony in armadillos[END_REF], limiting the opportunity of inbreeding. Our results therefore favour the founder effect hypothesis along colonization of North America. However, we cannot draw conclusions about the number and location of the bottlenecks without sampling key biogeographical intermediates -the Central American populations.

The rapid expansion of the nine-banded armadillo's range may have lead to establishment of populations from the front of colonization waves where the level of genetic variation decreases (e.g., [START_REF] Hewitt | The subdivision of species by hybrid zones[END_REF]. Moreover, recent human influence likely reinforced the founder effects with a hundred of releases leading to successful pioneer populations (see the example of Florida: [START_REF] Ramsey | Protein variation in populations of Dasypus novemcinctus and comparisons to D. hybridus, D. sabanicola and Chaetophractus villosus[END_REF][START_REF] Taulman | Recent range expansion and distributional limits of the nine-banded armadillo (Dasypus novemcinctus) in the United States[END_REF]. The genetic homogeneity of North American populations may thus reflect the combination of those phenomena.

Recovery of a mitochondrial polymorphism.

The mtDNA control region 5'-half of North American nine-banded armadillos is monomorphic for point mutations (except for one Louisiana individual), even for individuals separated by 1,000 km, but exhibits variable numbers of tandem TA repeats (VNTR: Table 1). Among placentals, microsatellites in the left domain of the control region are rare and always involve TA dinucleotides (in porpoises: [START_REF] Rosel | Phylogenetic relationships among true porpoises (Cetacea: Phocoenidae)[END_REF]pampas deer: Gonzalez et al. 1998). Measurement of VNTRs and levels of heteroplasmy will be developed to quantify the recovery of mitochondrial variability in United States armadillos.

Systematics of Dasypus novemcinctus.

MtDNA sequences reveal marked differences between French Guiana and the United States populations: five transversions allow discrimination of South and North American haplotypes (Fig. 1). Transversions in the control region are rare events at the population level (Wood & Phua 1996), though they may occur at a few consecutive sites [START_REF] Douzery | The mitochondrial control region of Cervidae: evolutionary patterns and phylogenetic content[END_REF]. This emphasizes the deepness of the split between South and North American haplotypes (Figs. 1 and2), indicating a rupture of gene flow between French Guiana and United States armadillos, in conjunction with a possible fixing of a divergent haplotype (Y or Z) after severe drift. Seven subspecies of Dasypus novemcinctus have been recognized (e.g., [START_REF] Mcbee | Dasypus novemcinctus[END_REF], with morphological and chromosomal differences identified [START_REF] Wetzel | The subgenera and species of long-nosed armadillos[END_REF][START_REF] Jorge | The evolution and ecology of armadillos, sloths, and vermilinguas[END_REF].

The abundant mitochondrial divergence we observed may reflect a subspecies differentiation. Molecular and morphological analyses on individuals from more distant localities (e.g. Brazil and Central America) and from congeneric species should be undertaken to further investigate the phylogeography and systematics of Dasypus novemcinctus.

Figure 1 .

 1 Figure 1. Minimum spanning network based on the number of substitutions between 10

Figure 2 .

 2 Figure 2. Highest-likelihood phylogram (ln L = -1,029.4) depicting the phylogenetic

  ( a ): heteroplasmy detected on sequencing gels autoradiographies (n = 9 to 10 for DNO-T1672; n = 13 to 15 for individual DNO-T1674). ( b ): imperfect repeat of the motif due to a pyrimidine difference. ( c ): the exact origin of the sample from the study by Arnason et al. (1997) is unknown. ( d ): the outgroup, Dasypus kappleri (great long-nosed armadillo), displayed two haplotypes. Huchon, Delsuc, Catzeflis & Douzery -Figure 1.

  Both haplotype and nucleotide diversities are significantly lower (Student's t test: t = 4.90; d.f. = 19; P < 0.0001 and t = 5.48; d.f. = 19; P < 0.0001, respectively) in United States than in French Guiana.

	at 27 loci. We screened 8 of these loci for 55 French Guianan
	individuals, and found that HBB and EST were polymorphic (H e = 0.102; n = 8 loci; SE
	= 0.227). Heterozygosities were compared, locus by locus, to published data on
	numerous placentals. HBB and EST are each more polymorphic in South American
	armadillos than in 29 placental species, with respectively H e (HBB) = 0.167 versus 0.031
	(n = 34 populations; SE = 0.099), and H e (EST) = 0.645 versus 0.232 (n = 113; SE =
	0.246). For the six remaining loci, the monomorphism observed in French Guianan
	armadillos conforms with other placental surveys: distributions of H e (ALB) , H e (GPI) , H e (LDH
	1+2) , H e (MDH) , and H e (PGD) display each a null median (n = 67, 95, 226, 79 and 112
	populations, respectively).
	mtDNA control-region polymorphism.
	Comparison of 492 unambiguously aligned sites of the control region reveals two
	(Y and Z) and 10 (A to J) haplotypes in United States and French Guiana populations
	respectively (Table 1). Haplotype diversity is 0.22 (n = 9 individuals; SE = 0.17) in
	United States, 0.97 (n = 12; SE = 0.04) in French Guiana, and 0.86 (n = 21; SE = 0.07)

Control region sequences begin with a microsatellite area for all Dasypus novemcinctus (Table

1

). In United States, they exhibit a polymorphic dinucleotide TA motif, varying from seven to 14 tandem repeats, and with heteroplasmy for the two Mississippi individuals. In French Guiana, they display a different pattern with one to five polymorphic hexanucleotide repeats surrounding an imperfect motif.

overall. Nucleotide diversity is 0.046% (SE = 0.034) in United States, 0.990% (SE = 0.146) in French Guiana, and 3.045% (SE = 0.187) overall.

TABLETable 1 .

 1 Origin and properties of the armadillo mitochondrial control region (D-loop) DNA sequences. We give the tissue collection number (Collection of Tissues of the "Paléontologie,

	Paléobiologie and Phylogénie" Laboratory, Montpellier II University; curated by F. Catzeflis),
	DNA extraction number, geographical origin, accession number (EMBL / GenBank / DDBJ
	databases), haplotype assigned letter, and structure of the microsatellite area for each
	sequence. All French Guiana individuals but three (DNO-38, DNO-46, and DNO-49) were
	studied by both allozyme and mtDNA control region sequencing approaches.
	Mississippi, Texas, Louisiania, and French Guiana samples are respectively from
	Hancock county (30°24'N; 89°39'W), Kimble county (30°26'N; 99°58'W), Landry Parish
	(30°25'N; 92°03'W), and Petit-Saut (04°48'N; 52°56'W, along the Sinnamary river).
	Twelve haplotypes (A to J; Y and Z) were defined by nucleotide substitutions. The
	sequence of the microsatellite area is [TA] n [CATATA] p [CAYAYATATA] [CATATA] q ,
	with variable n, p and q repeat numbers. The [CATATA] repeats surround an imperfect motif
	which consensus sequence is CAYAYATATA (with CAYAYA lacking in haplotype J from
	individual DNO-166), and that is absent from the control regions of individuals DNO-38,
	DNO-155, and DNO-163. The variable number of tandem TA repeats (VNTR) started after
	nucleotide position 17 of the control region for North America individuals, and after position
	8 for South America individuals. No polymorphic repeats were found in control region of the
	three great long-nosed armadillos.
	Individual	DNA	Origin	Accession		Micro-		Haplo-
		sample		numbers		satellite		types
					n	p	q
	DNO-T1672	#5132	Mississippi	AJ010362	9 a	0	0	Z
	DNO-T1674	#5134	Mississippi	AJ010363 14 a	0	0	Z
	DNO-T1676	#5136	Texas	AJ010364	13	0	0	Z

ACKNOWLEDGMENTS

This study would not have been feasible without the valuable biological samples preserved, curated and managed by many individuals: Nancy D. Moncrief (Collection of Tissues; Louisiania State University) Richard Monk (Cryogenics Facility; Texas Tech University), Danny Devillier "guitariste de Louisiane", Jean-Christophe Vié and his team "Faune Sauvage" at Petit-Saut. Laboratory facilities were provided by François Renaud. We wish to thank Jacques Lagnel, Anne Lavergne, and Francis Santalla for their valuable bench advices, and Mark W. Chase and Robin Smith for English suggestions. We also thank Christine Chevillon, Jim Loughry, Robert Wayne and two anonymous referees for helpful comments on the manuscript. Research on French Guiana mammals was funded by EDF-CNEH (contract GP-0850) and Ministère de l'Environnement (contract SOFT-Guyane). This is contribution ISEM 99-CCC of the Institut des Sciences de l'Evolution de Montpellier (UMR 5554 -CNRS).

genus Dasypus L. In: Vertebrate ecology in the northern Neotropics (ed Eisenberg JF), pp. 43-63. Smithsonian Institution Press, Washington.

Wood NJ, Phua SH (1996) Variation in the control region sequence of the sheep mitochondrial genome. Animal Genetics, 27, 25-33.


Stef Graillat
email: graillat@univ-perp.fr

Philippe Langlois
email: langlois@univ-perp.fr

Nicolas Louvet
email: nicolas.louvet@univ-perp.fr

Choosing a Twice More Accurate Dot Product Implementation

The fused multiply and add (FMA) operation computes a floating point multiplication followed by an addition or a subtraction as a single floating point operation. Intel IA-64, IBM RS/6000 and PowerPC architectures implement this FMA operation. The aim of this talk is to study how the FMA improves the computation of dot product with classical and compensated algorithms. The latters double the accuracy of the former at the same working precision. Six algorithms are considered. We present associated theoretical error bounds. Numerical experiments illustrate the actual efficiency in terms of accuracy and running time. We show that the FMA does not improve in a significant way the accuracy of the result whereas it increases significantly the actual speed of the algorithms.

The fused multiply and add (FMA) operation computes a floating point multiplication followed by an addition or a subtraction as a single floating point operation. This means that only one final rounding (to the working precision) error is generated by a FMA whereas two occur in the classical implementation of x × y + z. Intel IA-64, IBM RS/6000 and PowerPC architectures implement this FMA operation. On the Itanium processor, the FMA operation enables a multiplication and an addition to be performed in the same number of cycles than one multiplication or one addition [START_REF] Markstein | IA-64 and elementary functions. Speed and precision[END_REF].

The FMA operation seems to be advantageous for both speed and accuracy. Indeed, it approximately halves the number of rounding errors in many numerical algorithms. This is the case for example within the computation of a dot product of two n-vectors where just n rounding errors occur instead of 2n -1 without FMA.

Moreover, it is well known that FMA yields an efficient computation of the rounding error generated by a floating point product. Such rounding error computation at the current working precision is a key task when implementing multi-precision libraries as double-double or quad-double ones [START_REF] Bailey | A Fortran-90 double-double library[END_REF] or even when designing compensated algorithms. Compensated algorithms implement inner computation of the rounding errors generated by the original algorithms and so provide more accurate results; [START_REF] Pichat | Correction d'une somme en arithmétique à virgule flottante[END_REF][START_REF] Ogita | Accurate sum and dot product[END_REF] are examples of compensated summation and dot product. The latter reference recently proved that these compensated implementations double the accuracy of the classical algorithm still running in the current working precision.

Here we study how the FMA can improve the computation of dot products in terms of accuracy and running time.

First, we consider the classical dot product computed at the working precision with or without FMA. We report the theoretical error analysis (worst case bounds) and some experimental results to show that the use of FMA only slightly improve the accuracy of the computed dot product, even if the number of rounding errors is halved.

Nevertheless, the accuracy provided by the classical dot product may not be sufficient when applied to ill conditioned dot products. Such cases appear for instance when computing residuals for ill conditioned linear systems. So we also consider accurate dot products whose computed result is as accurate as if computed in twice the working precision. Here we consider the classical dot product performed with double-double computation as it can be found in the XBLAS library [START_REF] Li | Design, implementation and testing of extended and mixed precision BLAS[END_REF] and the compensated dot product from [START_REF] Ogita | Accurate sum and dot product[END_REF] where the FMA is used to compute the rounding error generated by each product. We also present a new compensated dot product using a recent algorithm by Boldo and Muller [START_REF] Boldo | Some functions computable with a fused mac[END_REF] that computes the exact result of a FMA operation as the unevaluated sum of three floating point values. We present theoretical error bounds to prove that all these algorithms provide results as accurate as if computed in twice the working precision. Then we compare these implementations in terms of practical computing time to identify the best choice to double the computing precision. Our experimental results show that the compensated algorithms run both considerably faster than the one with double-double computation. Moreover, the most efficient approach to benefit from the availability of a FMA seems to be to compensate the rounding error generated by each (classical) product without using the FMA operation in the inner loop of the dot product.

Algorithm

Table 1 :

 1 Classical dot product with and without FMA: the FMA does not improve the actual accuracy. Measured computing times (Itanium 2, 1600 MHz, Intel C++ Compiler v9.0). Using the FMA to compensate the multiplication rounding error (Dot2FMA) is the best choice to double the accuracy.

								Brief description
					Dot		Dot product without FMA
					DotFMA		Dot product with FMA
					Dot2FMA		Compensated dot product with FMA
					DotThreeFMA Compensated DotFMA with FMA
					DotXBLASFMA XBLAS dot product with FMA
				Classical dot product with FMA (n=100, 720 samples)				Classical dot product without FMA (n=100, 720 samples)
		1								1
		0.01								0.01
								DotFMA							Dot
		1e-04	1/2 γ n cond							1e-04	1/2 γ n cond
	relative forward error	1e-12 1e-10 1e-08 1e-06							relative forward error	1e-12 1e-10 1e-08 1e-06
		1e-14								1e-14
		1e-16							eps	1e-16					eps
		1e-18								1e-18
			100000	1e+10	1e+15	1e+20	1e+25	1e+30	1e+35		100000	1e+10	1e+15	1e+20	1e+25	1e+30	1e+35
					condition number							condition number
	Figure 1: n		DotFMA Dot2FMA DotThreeFMA DotXBLASFMA
				Theoretical			1	10		19		22
				Measured	50		1.0	1.4		2.3		8.24
						100		1.0	1.29		2.37		8.98
						1000		1.0	1.24		2.63		10.46
						10000		1.0	1.25		2.63		10.5
						100000	1.0	1.07		1.76		6.27


Ying Tang
email: ying.tang@univ-lille.fr

Alessio Franci
email: afranci@ciencias.unam.mx

Romain Postoyan
email: romain.postoyan@univ-lorraine.fr

On-line detection of qualitative dynamical changes in nonlinear systems: the resting-oscillation case

Keywords: transition detection, qualitative methods, singular perturbation, singularity theory, neuroscience, nonlinear systems, Lyapunov methods

Motivated by neuroscience applications, we introduce the concept of qualitative detection, that is, the problem of determining on-line the current qualitative dynamical behavior (e.g., resting, oscillating, bursting, spiking etc.) of a nonlinear system. The approach is thought for systems characterized by i) large parameter variability and redundancy, ii) a small number of possible robust, qualitatively different dynamical behaviors and, iii) the presence of sharply different characteristic timescales. These properties are omnipresent in neurosciences and hamper quantitative modeling and fitting of experimental data. As a result, novel control theoretical strategies are needed to face neuroscience challenges like on-line epileptic seizure detection. The proposed approach aims at detecting the current dynamical behavior of the system and whether a qualitative change is likely to occur without quantitatively fitting any model nor asymptotically estimating any parameter. We talk of qualitative detection. We rely on the qualitative properties of the system dynamics, extracted via singularity and singular perturbation theories, to design low dimensional qualitative detectors. We introduce this concept on a general class of singularly perturbed systems and then solve the problem for an analytically tractable class of two-dimensional systems with a single unknown sigmoidal nonlinearity and two sharply separated timescales. Numerical results are provided to show the performance of the designed qualitative detector.

Introduction

In neurosciences, to detect on-line the current activity type of a neuron (spiking, bursting etc.) or of a population of neurons (resting, ictal, interictal, slow/fast oscillations etc.) is of fundamental importance. For epilepsy for instance, such algorithms could be used to detect or even predict seizures. It seems that this problem is very challenging from a control-theoretic viewpoint for several reasons. First, disparate combinations of biophysical parameters are known to lead to the same activity pattern at the cellular level [START_REF] Goldman | Global structure, robustness, and modulation of neuronal models[END_REF], and the same degenerated parametrization property propagates at the A preliminary version of this work was presented at the 20th IFAC World Congress, Toulouse, France, 2017 [START_REF] Tang | Qualitative parameter estimation for a class of relaxation oscillators[END_REF]. This work was supported by the ANR under the grant SEPI-COT (ANR 12 JS03 004 01), by the French "Région Grand Est" through a fellowship grant 2016-2017, and by DGAPA-UNAM under the grant PAPIIT RA105518. neuronal circuit level [START_REF] Marder | Variability, compensation, and modulation in neurons and circuits[END_REF]. Second, biophysical parameters change smoothly over time but by doing so they induce sharp, almost discontinuous, transitions between qualitatively different activity types (spiking or bursting, healthy or epileptic, etc.) at the crossing of critical parameter sets. Third, the few available models are often subject to large uncertainties. As a consequence, quantitative modeling and fitting of experimental neuronal data generally constitute an ill-posed problem: a new viewpoint on the problem is needed, which is the purpose of this work.

A first major point is to extract the ruling parameters governing neuronal dynamics and their critical values. In spite of the nonlinear, degenerate structure of the biophysical parameter space, neuronal systems typically exhibit only a few distinct qualitative behaviors. For instance, the activity can be classified as resting, spiking, bursting at the single cell level, and as resting, low amplitude/high frequency, high amplitude/low frequency os-cillations at the network level. Predicting the structure of the critical parameter set, where switches in behavior occur, remains an open problem in general. Recently, based on singularity and singular perturbation theories, it was shown that the quantitative parameter space of biophysical neuronal models can be mapped to a small number of lumped parameters. The lumped parameters define low-order polynomial models, capturing the qualitative neuronal behavior at the single cell level [START_REF] Drion | Dynamic input conductances shape neuronal spiking[END_REF][START_REF] Franci | Modeling the modulation of neuronal bursting: a singularity theory appraoch[END_REF]. The resulting geometric framework allows a local, qualitative sensitivity analysis of neuronal dynamics [START_REF] Franci | Modeling the modulation of neuronal bursting: a singularity theory appraoch[END_REF][START_REF] Drion | Neuronal behaviors: a control perspective[END_REF], which predicts the effect of biophysical parameter variations on the qualitative behavior and the qualitative shape of the critical parameter set. Motivated by the above considerations and inspired by the qualitative sensitivity analysis in [START_REF] Franci | Modeling the modulation of neuronal bursting: a singularity theory appraoch[END_REF][START_REF] Drion | Neuronal behaviors: a control perspective[END_REF], our first objective is to extract analytically the parameters that rule the activity type of a highly uncertain, redundantly parametrized system and then to detect on-line the activity type and whether a dynamical change is likely to occur.

We formulate these ideas on a general class of nonlinear singularly perturbed control systems that embrace, for instance, the Hodgking-Huxley model for neuronal spiking [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] and the Wilson-Cowan model for neuronal population activity [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF]. The singularly perturbed nature of the model is justified by neuronal biology, where sharply separated timescales govern the neuronal dynamics [START_REF] Franci | Modeling the modulation of neuronal bursting: a singularity theory appraoch[END_REF][START_REF] Drion | Dynamic input conductances shape neuronal spiking[END_REF][START_REF] Drion | Neuronal behaviors: a control perspective[END_REF]. Because the proposed problem is very general and challenging, we then focus on solving it for a two dimensional class of nonlinear systems with a single unknown sigmoidal nonlinearity and two sharply separated timescales. This class of systems captures the qualitative features of the center manifold reduction of higher-dimensional biophysical models, as we explain. In order to cope with the peculiarities of neuronal dynamics, we assume that the exact functional form of the sigmoid nonlinearity is unknown as well as the exact timescale separation. We firstly show that, independently of the particular expression of the sigmoid nonlinearity, the system exhibits either global exponential stability (resting) or relaxation oscillations depending on a single ruling parameter. Guided by singularity theory [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF][START_REF] Franci | Realization of nonlinear behaviors from organizing centers[END_REF], we subsequently design an on-line qualitative detector that infers the activity type of the system and determines whether a qualitative dynamical transition is close to occur. By qualitative, we mean that the detector neither needs to quantitatively fit any model nor to asymptotically estimate any parameter. Because activity transitions are governed by bifurcations of the underlying vector field, the detection problem put forward here is closely related to the problem of steering a system toward an a priori unknown bifurcation point [START_REF] Moreau | Balancing at the border of instability[END_REF][START_REF] Moreau | Feedback tuning of bifurcations[END_REF].

The main contributions of the present work are summarized as follows.

•Inspired by neurosciences, we formulate the problem of detecting on-line qualitative changes in nonlinear sys-tems for which an exact model is not available.

•We provide a solution for a class of two-dimensional nonlinear systems, to detect whether the system is close to a transition between resting and oscillatory activity. The stability and the robustness of the detector are ensured via Lyapunov techniques and geometric singular perturbation arguments.

•We provide numerical evidences that the same detector performs well when tested on a classical, high dimensional, biophysical neuronal model of resting and spiking oscillations due to Hodgkin and Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF].

The paper is organized as follows. In Section 2, we first formulate the idea on a general nonlinear systems. We then introduce and analyze a class of singularly perturbed nonlinear systems with a single unknown sigmoidal nonlinearity. In Section 3, we design and analyze the qualitative detector for the latter. Numerical simulations are presented in Section 4. Section 5 concludes the paper. The proofs are given in the appendix.

Notation. Let R 0 := [0, ∞), R >0 := (0, ∞), and N := {1, 2, . . . , }. The usual Euclidean norm is denoted by | • |. For x, y, z ∈ R, the vector [x, y, z] is denoted by (x, y, z). For a function h : R 0 → R n , the associated infinity norm is denoted by h [0,∞) = sup s∈[0,∞) |h(s)|,
when it is well-defined. We use sgn(x) to denote the sign function from R to {-1, 0, 1} with sgn(0) = 0. For any function f : R → R, we denote the range of f as

f (R) = {z : z = f (x) for some x ∈ R}. Let A, B be two non-empty subsets in R n , their Hausdorff dis- tance is noted by d H (A, B) = max{sup a∈A inf b∈B |a - b|, sup b∈B inf a∈A |a-b|}. The point to set distance from p ∈ R n to A ⊆ R n is denoted by |p| A = inf a∈A |p -a|.
2 Problem statement

General formulation

We first introduce the main ideas on the following general class of nonlinear singularly perturbed systems

   ẋf = f (x f , h s (x s ), α, u), ẋs = εg(h f (x f ), x s ), y = (h f (x f ), h s (x s )), (1a) (1b) (1c) where x f ∈ R n , x s ∈ R m are the state variables, α ∈ R p is an unknown parameter vector, u ∈ R is a constant input, 0 < ε
1 is an unknown parameter, h s : R m → R, h f : R n → R and all functions are smooth.

System (1) is said to be redundantly parametrized if it depends on numerous parameters, i.e. p 1, but the possible qualitatively distinct dynamical behaviors are only a few. Parametric redundancy naturally arises under two conditions: timescale separation and the presence of an organizing center, as explained below. System (1) evolves according to two time scales, since ε is small. Under suitable stability and monotonicity conditions of the fast dynamics (see e.g. [START_REF] Gedeon | Oscillations in multistable monotone systems with slowly varying feedback[END_REF]), the trajectories of the fast variable x f are forced toward its instantaneous quasi-steady state x * f (t), which is defined as follows. For all t ∈ R ≥0 , α ∈ R p , u ∈ R, let x s (t) be the slow state at time t. Then x * f (t) is a solution of the quasi-steady state equation f (x * f (t), h s (x s (t)), α, u) = 0, and the fast output y f = h f (x f ) will change accordingly. Now, for u fixed, define the map

G : R n+1 × R × R p → R n+1 as G(x, λ, α) := f (x f , λ, α, u) y f -h f (x f )
, where x = (x f , y f ) and λ is the bifurcation parameter.

For system (1) to be redundantly parameterized, we assume that the map G is organized by a k-codimensional singularity, with k p, in the sense of [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]Section III.1]. This means that the qualitative shape of the graph of the, possibly multivalued, mapping λ → x f defined by solving the equation G(x, λ, α) = 0 is given by one of the persistent bifurcation diagrams [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]Sections III.5,[START_REF] Drion | Neuronal behaviors: a control perspective[END_REF] in the universal unfolding of the organizing singularity. This shape is exactly what determines the quasi-steady state evolution of the fast variable. Under suitable stability and monotonicity properties of the slow dynamics, (for instance that they provide slow adaptation -negative feedback -on the fast dynamics), it follows that the possible qualitative dynamic output behaviors of system (1) are fully determined by the persistent bifurcation diagrams of the organizing singularity. Hence, since k p, the possible different persistent bifurcation diagrams, and therefore the possible output behaviors are much less than the number of model parameters: parametric redundancy occurs.

When system (1) exhibits parametric redundancy, we do not need to know the full vector of parameters α ∈ R p to detect what is the current dynamic behaviour of the system and whether a qualitative change is prone to appear. The projection of α onto the k-dimensional unfolding space of the organizing suffices to determine the model behavior. The sensitivity analysis proposed in [START_REF] Drion | Dynamic input conductances shape neuronal spiking[END_REF] is fully based on this idea. In this paper, we explore this viewpoint to formulate the following qualitative detection problem: Let f , g, h f , h s be unknown. Let the organizing singularity be known. Let the slow and fast output be measurable. Let the dimension p of the parameter space be unknown. Can we detect on-line in which region of the organizing singularity unfolding space the parameter vector lies and whether this vector is approaching a transition variety, where a qualitative dynamical transition occurs? This problem is new, general and very challenging. That is the reason why we solve it for a particular yet relevant class of systems of the form of (1) in the following.

A tractable class of systems

In the following, we concentrate on systems of the form

ẋf = -x f + S(βx f + u -x s ), ẋs = ε(x f -x s ), (2a) (2b)
where x f , x s ∈ R are the state variables, β ∈ R is an unknown parameter, u ∈ R is the input, which is assumed to be known and constant, and 0 < ε 1 is small unknown parameter. The mapping S : R → R is an unknown sigmoid function, which is assumed to satisfy the following properties. Due to the small parameter ε, system (2) evolves according to two time scales. We follow the standard approach of singular perturbation theory [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF][START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] to decompose system (2) into two subsystems, which represent the fast dynamics and the slow dynamics, respectively, called the layer and the reduced subsystems. By setting ε = 0 in (2b), we obtain the layer dynamics

ẋf = -x f + S(βx f + u -x s ), ẋs = 0. (3a) (3b)
Because of (3b), the variable x s can be treated as a parameter in (3a).

To account for the slow variations of x s , we rescale the time as τ = εt, hence d dt = ε d dτ . Then, system (2) becomes

εx f = -x f + S(βx f + u -x s ), x s = x f -x s , (4a) (4b)
where stands for d dτ . Setting ε = 0 in (4a), we obtain the reduced dynamics

     0 = -x f + S(βx f + u -x s ) =: F f (x f , u -x s , β), x s = x f -x s =: F s (x f , x s ). (5a) (5b)
The reduced dynamics evolves in the slow time τ and is an algebro-differential equation. In particular it defines a one-dimensional vector field over the critical manifold

S 0 u,β := {(x s , x f ) ∈ R 2 : F f (x f , u -x s , β) = 0}. (6)
The critical manifold S 0 u,β depends on u and β. However, since we assume u and β constant, for simplicity we omit the index u and β in the rest of the paper. The critical manifold will therefore simply be denoted by S 0 . The idea behind singular perturbations is that dynamics (2) are well approximated by layer dynamics far from S 0 and by trajectories of the reduced dynamics close to S 0 . The shape of S 0 therefore plays a key role in determining the behavior of system (2).

Let

β c := 1 S (0) , (7)
which is well-defined and strictly positive according to item c) of Assumption 1. The parameter β c is crucial in shaping the critical manifold S 0 , as illustrated in Figure 1. For β < β c the shape of S 0 is given by graph of a monotone decreasing function. Whereas for β > β c the shape of S 0 is S-shaped. We emphasize that β c is unknown since so is the sigmoid function S. The qualitative properties of the critical manifold fully determine the model dynamical behaviors, as formalized in the next proposition, whose proof can be found in [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF].

Proposition 1 Consider system (2), the following holds.

1) For any β -β c < 0 and any constant input u, there exists ε > 0 such that, for all ε ∈ (0, ε], system (2) has a globally exponentially stable fixed point. Moreover, all trajectories converge in an O(ε)-time to an O(ε)neighborhood of the critical manifold S 0 defined in (6).

2.a) For all 0 < β -β c < 1, there exists ū > 0, such that, for any constant input u in (-ū, ū), there exists ε > 0 such that, for all ε ∈ (0, ε], system (2) possesses an almost globally asymptotically stable1 and locally exponentially stable periodic orbit P ε . 2.b) Let β, ū, u, ε, ε be as in item 2.a), and let T ε be the period of P ε and p : R → P ε , with p(t) = p(t+T ε ) for all t ∈ R, be a periodic solution of (2). Let P 0 be the singular limit of P ε . The Hausdorff distance between P ε and P 0 verifies d H (P ε , P 0 ) = O(ε 2/3 ). Moreover, for all t 0 ∈ R, there exists

δT ε ⊂ [t 0 , t 0 +T ε ) and |δT ε | = O(ε)T ε , such that |p(t) -S 0 | < O(ε) for all t ∈ [t 0 , t 0 + T ε ) \ δT ε . 2
Proposition 1 shows that system (2) either has a globally exponentially stable fixed point (item 1)), which we call the resting activity, or almost all its solutions converge to a limit cycle (item 2)), in which case we talk of oscillating activity. The nature of the activity is solely determined by the sign of β -β c , independently of the particular expression of S.

Remark 1 Proposition 1 states stability results for β < β c + 1. When β > β c + 1, system (1) exhibits three equilibria, two of which are asymptotically stable and one is unstable. Then trajectories converge to the two stable equilibria. 2

Remark 2 Higher-dimensional control systems of the form (1) may be reduced to the simpler lowerdimensional system (2) via center manifold reduction [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF]. We sketch the ideas as follows. The transition from globally asymptotically stable fixed point to almost globally asymptotically, locally exponentially stable oscillations in model ( 2) is via a Hopf bifurcation, as it can easily be verified by linearization at the origin (for u = 0). Assume that the same bifurcation exists in the full model [START_REF] Tang | Qualitative parameter estimation for a class of relaxation oscillators[END_REF]. Then, close to this bifurcation, by the center manifold theorem, the dynamics of model ( 1) can be reduced to a two-dimensional invariant manifold: the center manifold of the Hopf bifurcation. The same assumptions as for (2) can then be formulated for the full model ( 1) by restricting them to the Hopf center manifold. Classical two-dimensional reduction of the Hodgkin-Huxley model [START_REF] Fitzhugh | Impluses and physisological states in theoretical models of nerve membrane[END_REF][START_REF] Rinzel | Excitation dynamics: Insights from simplified membrane models[END_REF] satisfy these conditions with the sodium maximal conductance ḡNa as unfolding parameter and the applied current I app as constant input. To illustrate this, numerical results for Hodgkin-Huxley model are provided in Section 4.2 2

In the next section we reformulate the qualitative detection problem in the specific context of model ( 2) and we provide a solution.

3 Qualitative detector

Objective and design

Our objective is to detect on-line the current type of activity of system (2), i.e. whether the solution converges to a fixed point or to a limit cycle, and how close the system is to a change of activity. In view of Proposition 1, activity switch in system ( 2) is ruled by β only. Hence, the most natural idea is to estimate β. For this purpose, we might approximate the sigmoid function S by a Taylor expansion up to some order and quantitatively estimate the expansion coefficients together with ε and β. However, due to the generic impossibility of quantitatively modeling and fitting experimental neuronal data, this approach would not overcome, in a real experimental setting, the limitation of quantitative/asymptotic estimation. Instead, we aim at qualitatively detection the different dynamical activities of the system: whether the system exhibits oscillation or has a globally exponentially stable fixed point and whether it is near to a change of activity.

Proposition 1 proves that (almost) all trajectories of system (2) converge to an O(ε)-neighborhood of the critical manifold S 0 in both the fixed point and the limit cycle cases, at least for most of the time in the latter case. Based on singularity theory [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF][START_REF] Franci | Realization of nonlinear behaviors from organizing centers[END_REF] and in view of Assumption 1, the shape of the critical manifold S 0 is, modulo an equivalence transformation, the same as that of the set

{(x f , x s ) : -x 3 f + (β -β c )x f + u -x s = 0}, ( 8
)
whose definition is independent of S. We exploit this information to design and analyze our detector in the following.

We first make the following assumption.

Assumption 2 Both x f and x s variables are measured in system (2).

Remark 3 When only the fast dynamics is available by measurement, we may reconstruct the slow variable using a state estimator of the form ẋs = ε(x f -xs ), for some small ε > 0. The analysis of this scenario is not addressed in this paper and is left for future work. 2

We propose the nonlinear qualitative detector

β = -kx f (-x 3 f + βx f + u -x s ) =: f ( β, x f , u -x s ), (9)
where k > 0 is a design parameter. The definition of f does not require the knowledge of S or ε, it only relies on the information provided by [START_REF] Wilson | Excitatory and inhibitory interactions in localized populations of model neurons[END_REF].

In the following, we define β * the steady state of β, and we analyze its relationship with respect to β -β c . We use the terminology steady-state with some slight abuse as we will see that β * is not always a constant. We then investigate the incremental stability properties of ( 9), which we finally exploit to prove the stability of ( 9) and the convergence of β to β * .

Steady state properties of the qualitative detector

Steady state of β with x f = 0 satisfies -x 3 f + βx f + u -x s = 0 in view of [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]. We thus implicitly define the steady state function β * (x f , x s , u) by

-x 3 f + β * (x f , x s , u)x f + u -x s = 0, (10)
for u, β ∈ R and x f , x s ∈ S 0 , i.e. when x f , x s are related via the critical manifold equation

-x f + S(βx f + u -x s ) = 0. ( 11
)
Recalling that the sigmoid function S is invertible on S(R) according to item c) of Assumption 1, we can explicitly invert (11) as follows

x s = -S -1 (x f ) + βx f + u. ( 12
)
Replacing this expression for x s in [START_REF] Franci | Realization of nonlinear behaviors from organizing centers[END_REF], we obtain the explicit expression of β * (x f , x s , u)

β * (x f , x s , u) = x 3 f -S -1 (x f ) + βx f x f . ( 13
)
Because β * depends only on x f and β in (13), we write β * (x f , β) in the sequel.

The next theorem provides important properties of β * (x f , β). In particular, it ensures that β * (x f , β) is well-defined when x f = 0. Furthermore, it formally characterises in which sense β * provides qualitative information about β -β c .

Theorem 1 The following holds. 1) β * is smooth on S(R) × R.

2) For all β < β c and x f ∈ S(R),

∂ β * (x f ,β) ∂β > 0.
3) For all β ∈ (0, β c + 1) and x f in a neighborhood of the origin,

∂ β * (x f ,β) ∂β > 0. Moreover, let x + f old > 0 and x - f old < 0 be the two distinct solutions to -1 + βS (S -1 (x f )) = 0. If, for x f old ∈ {x + f old , x - f old }, 2x 2 f old + S -1 (x f old ) x f old -(S -1 ) (x f old ) > 0, ( 14
)
then ∂ β * (x f ,β) ∂β > 0 for x f ∈ (-∞, x - f old ) ∪ (x + f old , +∞). 4) β * (x f , β) = β -β c + O(x 2 f ) for any x f ∈ S(R) and β ∈ R. 2
Theorem 1 reveals that β * provides information about the qualitative trend of β. According to item 2) of Theorem 1, an increase in β * indicates that β is increasing in the resting activity (i.e. β < β c ). If a change of activity is prone to appear (i.e. β is close to β c ), then necessarily x f will lie in a neighborhood of the origin (for u small), which relates to the point of infinite slope, corresponding to the hysteresis singularity of the critical manifold. In this case, β * will change its sign from negative to positive by virtue of item 4) of Theorem 1, and then it will keep increasing in β in the oscillation activity (i.e. β c < β < β c + 1) in view of item 3) of Theorem 1.

The same holds when β * goes from positive to negative values. Hence, properties of β * ensure that we can detect online when β is about to cross, either from above or from below, the critical value β c , at least for small u. Condition ( 14) can be verified off line for different sigmoid functions and might be instrumental in tracking variations of β for β -β c strictly positive and large.

Incremental stability properties of the detector

If the value of β converges to β * as time tends to infinity, then [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] can inform us about the activity of system (2) as explained above. To prove this, we first show that system (9) exhibits incremental stability properties [START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF] if the input signal x f to the detector is bounded and a persistency of excitation (PE) condition holds. We then prove that under mild conditions, the solutions to system (2) do satisfy the aforementioned PE requirement.

We use the following notion of PE, which is similar to the one in [START_REF] Anderson | Exponential stability of linear equations arising in adaptive identification[END_REF].

Definition 1

The piecewise continuous signal y : R 0 → R is (T, µ)-PE, if there exist T, µ > 0 such that t+T t y 2 (τ ) dτ µ for all t 0. 2

Proposition 2 For any ∆, M, T, µ > 0, there ex-

ist 1 (∆, M, T, µ), 2 (∆, M, T, µ), 3 (∆, M, T, µ) > 0 such that, if (i) x f,i is (T, µ)-PE, i ∈ {1, 2}, and (ii) x f,i , x s,i , u i are piecewise continuous with max( x f,i ∞ , x s,i ∞ , u i ∞ ) M , i ∈ {1, 2}, then for any | βi (0)| ∆, i ∈ {1, 2}
, the corresponding solution βi to system (9) satisfies for all t 0

| β1 (t) -β2 (t)| 1 (∆, M, T, µ)e -2(∆,M,T ,µ)t | β1 (0) -β2 (0)| + 3 (∆, M, T, µ) x f 1 -x f 2 [0,t) + x s1 -x s2 [0,t) + u 1 -u 2 [0,t) . ( 15
)
2 Proposition 2 states that system (9) satisfies a semiglobal input-to-state incremental stability property, when its input signals verify conditions (i) and (ii).

In the next lemma, we analyze under what conditions the solution to system (2) verifies these two conditions. Lemma 1 Consider system (2). 1) For any β -β c < 0, any ∆ > δ > 0 and any constant input u satisfying δ < |u| < ∆, there exist ε, T (∆, δ), µ(∆, δ), C(∆, δ) > 0, such that for all ε ∈ (0, ε] and initial conditions x f (0) and x s (0)

with |x f (0)|, |x s (0)| ∆, the corresponding solution to system (2) is such that (i) its x f -component is (T (∆, δ), µ(∆, δ))-PE and (ii) max( x f [0,∞) , x s [0,∞) , u [0,∞) ) C(∆, δ).
2) For any 0 < β -β c < 1, any ∆ > δ > 0 there exists ū ∈ (0, ∆) such that for any u ∈ (-ū, ū), there exist ε, T (∆, δ), µ(∆, δ), C(∆, δ) > 0, such that for all ε ∈ (0, ε] and initial conditions x f (0) and x s (0) with

|x f (0)|, |x s (0)| ∆ and |x f (0)-p * f | δ, |x s (0)-p * s | δ, where (p * f , p * s )
is the unique unstable fixed point of ( 2), the corresponding solution to system (2) satisfies conditions (i) and (ii) as in the above item 1).

2

Lemma 1 means that the conditions of Proposition 2 are verified semiglobally by solutions to system (2) whenever β < β c + 1, under mild conditions on the input u.

Main result

To finalize the analysis, we consider the overall system, which is a cascade of system (2) with ( 9)

     ẋf = -x f + S(βx f + u -x s ), ẋs = ε(x f -x s ), β = -kx f (-x 3 f + βx f + u -x s ). (16a) (16b) (16c)
We state the stability property of system [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF] and the relationship between β and β * in the following two theorems respectively.

Theorem 2 Consider system [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF]. 1) For any β < β c , ∆ > δ > 0, and any constant input u satisfying δ < |u| < ∆, there exists ε > 0 such that, for any ε ∈ (0, ε], system (15) has a globally asymptotically stable and locally exponentially stable fixed point.

2) For any 0 < β -β c < 1, ∆ > 0, there exist ε > 0 and ū ∈ (0, ∆) such that for any u ∈ (-ū, ū), and for any ε ∈ (0, ε], system (15) possesses an almost globally asymptotically stable and locally exponentially stable periodic orbit Q ε with period T ε . 2

Theorem 2 shows that the solution to the overall system (16) converges either to a stable fixed point or, almost all of these, converge to a limit cycle depending on the value of β. In the next theorem, we specify the relationship between β and β * .

Theorem 3 Consider system [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF]. 1) For any β < β c , ∆ > δ > 0, and any constant input u satisfying δ < |u| < ∆, there exists ε > 0 such that, for any ε ∈ (0, ε] and any initial conditions

|x f (0)|, |x s (0)|, | β(0)| ∆, β(t) converges to β * (x * f , β
) as time goes to infinity, where x * f is x f -component of the globally asymptotically stable and locally exponentially stable fixed point of ( 16), as stated in item 1) of Theorem 2.

2) For any 0 < β -β c < 1, ∆ > δ > 0, there exist ε > 0 and ū ∈ (0, ∆) such that for any u ∈ (-ū, ū), any ε ∈ (0, ε] and any initial conditions

|x f (0)|, |x s (0)|, | β(0)| ∆ and |x f (0) -p * f | δ, |x s (0) -p * s | δ, where (p * f , p * s ) is the unique unstable fixed point of system (16a)-(16b), β(t) converges in an O(ε)-neighborhood of β * (x lc f (t + θ), β) for all t ∈ [t 0 , t 0 + T ε ) \ δT ε with |δT ε | = O(ε)T ε , where x lc f (t + θ) is x f -component of the T ε -
periodic almost globally asymptotically stable and locally exponentially stable periodic orbit of ( 16), as stated in item 2) of Theorem 2 with phase θ ∈ [0, 2π).

2

When β < β c , item 1) of Theorem 3 shows that the value of β is β * when time goes to infinity. For 0 < β -β c < 1, item 2) of Theorem 3 indicates that β is in an O(ε)-neighborhood of β * for all the time, after a sufficiently long time, except during jumps of length

|δT ε | = O(ε)T ε .
We see that β * is not a constant in this case but a state-dependent signal. This is not an issue as all what matters is the sign of β * , which is informative according to Theorem 1, and not its actual value. Various filtering techniques, like simple averaging over a sliding window, can be used to extract an approximated constant value of β * if needed in this case, see Section 4.1 for instance.

Remark 4 In view of [START_REF] Sontag | New characterizations of input to state stability[END_REF], the stability property of system [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF] in Theorem 2 implies its local input-to-state stability with respect to any small exogenous disturbance. In other words, the detector is robust to small measurement noises as well as small perturbations in its dynamics. 2

4 Numerical illustrations

System (2)

We consider several sigmoid functions, namely

S 1 : x → tanh(x), S 2 : x → e -e -x -e -1 , S 3 : x → 1 1+e -x -1 2 .
The corresponding critical values are β c1 = 1, β c2 = e and β c3 = 4. We choose ε = 0.001, u = -0.01 and k = 5 in (9). Figure 2 illustrates the fact that, in spite of the different sigmoid nonlinearities, when β -β c < 0, the state x f converges to a constant value as time increases, as well as x s , which is consistent with item 1) of Proposition 1. When 0 < β -β c < 1, both x f and x s converge to an oscillatory behavior, which is in agreement with item 2) of Proposition 1.

To check that the proposed algorithm is able to detect a change of activity of the system, we use a varying signal for β and we consider S = S 1 and the same data as above. In particular, β(t) = 0.5 on t ∈ [0, 300], it increases with slope 0.005 on t ∈ [300, 425], and β(t) = 1.2 on t ∈ [425, 1000], which leads to a change of sign of β-β c at t = 400, since β = β c1 = 1 in this case. We observe in Figure 3a that, when β -β c < 0, β converges to a constant value as time increases, corresponding to the resting activity. When 0 < β -β c < 1, β tends to a periodic function, which is strictly positive. Hence, it indicates the oscillation activity. The spikes seen in Figure 3a are due to the "jump" of the solution to system (2) from one stable branch of the critical manifold to the other. This phenomenon is captured by item 2) of Theorem 3, as β is guaranteed to be close to β * for all time except, periodically, over interval of length |δT ε | = O(ε)T ε . These spikes can simply be removed by using a low pass filter on β as illustrated in Figure 3b. Moreover, the value of β evolves from negative to positive when the change of sign of β -β c occurs. In other words, the qualitative detector is able to detect the current activity and whether a change is occuring. We have tested different values of k in [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]. The simulations indicate that the speed of convergence of β increases with k. However, this also leads to bigger spikes at "jumps", which may provide wrong information during a very short interval, which may again be moderated by a low pass filter.

The same tests have been done for S 2 and S 3 . We emphasize that, even though the nonlinearities are different, the detector remains the same as defined in [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]. Figure 4 shows the relationship between β and β -β c for different input u and perturbation parameter ε. When β -β c < 0, the value of β provided in Figure 4 is the constant value to which it converges as seen in simulations. When 0 < β -β c < 1, the value of β reported in the figure corresponds to the average value of β. We observe that for β -β c < 0, β is negative and it increases as β -β c increases. Moreover, because the fixed point of system (2) is near zero for small u, β evolves almost linearly with respect to β -β c . When 0 < β -β c < 1, β is positive. In addition, when β -β c = 0, β is around the origin. Hence β is able to detect the current activity type and whether a change is likely to occur. Figure 4 illustrates the efficiency of the approach.

To further evaluate the robustness of the scheme, we have added small additive measurement noises given by d out (t) = 0.005 sin(50t) for x f , x s and d in = 0.008 for input u. It is found from Figure 5 that the detector still provides good results in this case, which is in agreement with Remark 4.

Hodgkin-Huxley model

We consider Hodgkin-Huxley model

[7]                    C V = -g K n 4 (V -V K ) -g N a m 3 h(V -V N a ) -g l (V -V l ) + I app , τ m (V ) ṁ = -m + m ∞ (V ), τ n (V ) ṅ = -n + n ∞ (V ), τ h (V ) ḣ = -h + h ∞ (V ), (17a) (17b) (17c) (17d)
where the variable V is the membrane potential. The sodium fast activation is m and its slow inactivation is h. The potassium slow activation is n. The input is the applied current I app . The parameters V K , V N a are the equilibrium potentials for the potassium and sodium ions and V l is the potential at which the leakage current is zero. The constant C is the membrane capacity. The sodium, potassium and ionic conductances are denoted by g N a , g K , and g l . The gating variable time constants are τ m , τ n , τ h and their steady state characteristics are m ∞ , n ∞ , h ∞ , which are all monotone sigmoidal functions. We consider the same parameter values as in [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] unless otherwise specified.

As explained in Remark 2, Hodgkin-Huxley model can be reduced to two dimensions by assuming that m ≡ m ∞ (V ) and that h ≡ l(n), for a suitable linear function l, which is fine close to the Hopf bifurcation point via center manifold reduction [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF], see e.g., [START_REF] Fitzhugh | Impluses and physisological states in theoretical models of nerve membrane[END_REF][START_REF] Rinzel | Excitation dynamics: Insights from simplified membrane models[END_REF]. In its twodimensional reduced form, Hodgkin-Huxley model falls into the class of systems (1), with V the fast variable and n the slow variable. The fast nullcline can also be shown to be either monotone or cubic depending on g N a and, since n ∞ is monotone, the slow variable nullcline can locally be approximated as being linear. We therefore expect the detector designed for the class of system (2) to work locally for the Hodgkin-Huxley model with u = I app and β = g N a .

In the following, we apply the proposed algorithm to the full Hodgkin-Huxley model in [START_REF] Angeli | An almost global notion of input-tostate stability[END_REF] to detect changes due to the ruling parameter g N a . We assume that the state variables V, n are measurable. The input I app is a constant and is assumed to be known. From ( 9), the detector for the Hodgkin-Huxley model is

ġNa = -kV (-V 3 + ĝNa V + I app -m), (18)
where k > 0.

In the simulations, the state variables V, m are centered at zero by using a high-pass filter. We choose k = 4, I app = 20. The values of g N a are selected from the set {23, 25, 27, 29, 31}, which gives g N a -g N ac {-4, -2, 0, 2, 4} corresponding to an approximate critical value g N ac 27. Figure 6 shows the asymptotic value of ĝNa as a function of g N a -g N ac , just like we did to generate Figures 4-5 in Section 4.1. The blue line corresponds to the nominal case I app = 20.

Red and purple lines correspond to the perturbed cases I app = 20 ± 2. When g N a -g N ac > 0, ĝNa is selected as the average value over a period of the asymptotic periodic function to which it converges. We observe that in the nominal case, ĝNa is negative for g N a -g N ac < 0 and it increases as g N a -g N ac increases. The value of ĝNa is positive when g N a -g N ac > 0. Moreover, ĝNa is around zero when g N a = g N ac . Hence ĝNa provides information about the actual activity of the model and whether the model is prone to change its activity.

The perturbed cases also show that the detector performance is robust to small input uncertainties. These preliminary numerical results highlight the potentiality of the proposed approach beyond the academic example thoroughly analyzed in the present paper.

Conclusion

We have introduced the concept of qualitative detection, as the problem of informing on-line the qualitative dynamical behavior of a multiple-time scale nonlinear system, independently of large uncertainties on the system nonlinearities and without using any quantitative fitting of measured data. This is achieved by first extracting the system ruling parameter(s) and by subsequently designing a qualitative detector, which determines the system activity type and whether a qualitative change in the activity is close to occur. We have presented this idea on a general class of nonlinear singularly perturbed systems. As a first application, we have focused on a class of two-dimensional nonlinear systems with two time-scales and a single nonlinearity either exhibiting resting or relaxation oscillation behaviors. We have illustrated the extension of the proposed detector to real-word application through numerical simulations. Future extensions will include generalizing the detector design problem to systems exhibiting more than two possible qualitatively distinct behaviors. According to [START_REF] Gedeon | Oscillations in multistable monotone systems with slowly varying feedback[END_REF]

, β * (x f , β) = x 2 f - S -1 (x f ) x f
+ β, which is well-defined and differentiable on S(R) \ {0} × R. The only critical term for the smoothness of β * with respect to x f is

S -1 (x f )
x f when x f = 0. In view of items a) and c) of Assumption 1, the function S -1 is differentiable. Then the Taylor expansion of S -1 (x f ) at the origin is, for any

x f ∈ S(R) S -1 (x f ) =S -1 (0) + (S -1 ) (0)x f + (S -1 ) (0) 2 x 2 f + O(x 3 f ). ( 19
)
In view of item b) of Assumption 1 and (7), since S(S -1 (0)) = 0, it holds that S -1 (0) = 0, (S -1 ) (0) =

S (0) = β c > 0. Due to item d) of Assumption 1, (S -1 ) (0) = -S (0) (S (0)) 2 = 0. We deduce from (19), S -1 (x f ) = βx f + O(x 3 f ). Then the term S -1 (x f ) x f be- comes S -1 (x f ) x f = β c + O(x 2 f ), 1
which is smooth with respect to x f . This concludes the proof of the first item.

2) Let β ∈ R, we first consider that β < β c and x f = x s ∈ S 0 . By implicitly differentiating [START_REF] Moreau | Feedback tuning of bifurcations[END_REF] with respect to β, for x f ∈ S(R),

∂x f ∂β = x f -(β -1) + (S -1 ) (x f ) . (21)
Using ( 13) and ( 21), we obtain

∂ β * (x f , β) ∂β = ∂ β * (x f , β) ∂β + ∂ β * (x f , β) ∂x f ∂x f ∂β = 1 + 2x f - (S -1 ) (x f ) x f + S -1 (x f ) x 2 f x f -(β -1) + (S -1 ) (x f ) = 1 -β + 2x 2 f + S -1 (x f ) x f -(β -1) + (S -1 ) (x f ) . (22)
By item c) of Assumption 1, (S -1 ) (

x f ) = 1 S (S -1 (x f )) 1 S (0) = β c . Since β -β c < 0, we have -(β - 1) + (S -1 ) (x f ) β c -β + 1 > 0. Then the term 1 -(β-1)+(S -1 ) (x f )
in the right-hand side of ( 22) is strictly positive. According to the mean value theorem, it holds

S -1 (x f )-S -1 (0) x f -0 = S -1 (x f ) x f = (S -1 ) (c) for some c such that c ∈ (x f , 0) if x f < 0 and c ∈ (0, x f ) if x f > 0. Using again item c) of Assumption 1, it holds (S -1 ) (c) = 1 S (c) 1 S (0) = β c
. Hence, we deduce from [START_REF] Anderson | Exponential stability of linear equations arising in adaptive identification[END_REF] with the fact that β c -β > 0,

∂ β * (x f ,β) ∂β (1+βc-β+2x 2 f ) -(β-1)+(S -1 ) (x f ) > 0.
3) Let β ∈ (0, β c + 1), for any fixed x s ∈ R on S 0 and u ∈ R, by implicitly differentiating [START_REF] Moreau | Feedback tuning of bifurcations[END_REF] with respect to β, we have

∂x f ∂β = x f
-β+(S -1 ) (x f ) , which is well-defined as long as x f = {x - f old , x + f old }. We deduce from [START_REF] Gedeon | Oscillations in multistable monotone systems with slowly varying feedback[END_REF] and

∂x f ∂β = x f -β+(S -1 ) (x f ) , ∂ β * (x f , β) ∂β = ∂ β * (x f , β) ∂β + ∂ β * (x f , β) ∂x f ∂x f ∂β = 1 + 2x 2 f -(S -1 ) (x f ) + S -1 (x f ) x f -β + (S -1 ) (x f ) . (23)
We use again

S -1 (x f ) x f = (S -1 ) (c) for some c such that c ∈ (x f , 0) if x f < 0 and c ∈ (0, x f ) if x f > 0, and write (S -1 ) (c) = (S -1 ) (x f ) + (S -1 ) (c) -(S -1 ) (x f ).
Substituting this expression into [START_REF] Sontag | New characterizations of input to state stability[END_REF], we obtain

∂ β * (x f , β) ∂β = 1 + (S -1 ) (c) -(S -1 ) (x f ) + 2x 2 f -β + (S -1 ) (x f ) , (24) since
(S -1 ) (c) = β c + O(c 2 ), (S -1 ) (x f ) = β c + O(x 2 f ) and |c| < |x f |, it follows from (24) that ∂ β * (x f ,β) ∂β = 1 + O(x 2 f
). This proves the first part of item 3).

We now prove the second part of item 3). Let

x f old ∈ {x + f old , x - f old }. By definition of x + f old , x - f old in Lemma 1, it holds -1 + βS (S -1 (x f old )) = 0, which is equivalent to β = (S -1 ) (x f old ), (25) as
(S -1 ) (x f old ) = 1 S (S -1 (x f old )) .
We deduce from (23) by using (25),

∂ β * (x f , β) ∂β = 1 + 2x 2 f -(S -1 ) (x f ) + S -1 (x f ) x f -β + (S -1 ) (x f ) = -β + 2x 2 f + S -1 (x f ) x f -β + (S -1 ) (x f ) = -(S -1 ) (x f old ) + 2x 2 f + S -1 (x f ) x f -β + (S -1 ) (x f ) . (26)
By virtue of items c) and d) of Assumption 1, the following holds, for

x f ∈ (-∞, x - f old ) ∪ (x + f old , +∞) -1 + βS (S -1 (x f )) < 0. (27)
Dividing both sides of the above inequality by S (S -1 (x f )), which is strictly positive, and using (S -1 ) (x f ) = 1 S (S -1 (x f )) , we obtain

-β + (S -1 ) (x f ) > 0. ( 28
)
This implies that the denominator of ( 26) is strictly positive. We denote g(x f ) =

(S -1 )(x f ) x f
. The derivative of g(x f ) is given by

g (x f ) = x f (S -1 ) (x f ) -S -1 (x f ) x 2 f . (29)
We show below that the growth of g(x f ) increases with

|x f |. According to item c) of Assumption 1, it holds (S -1 ) (x f ) = 1 S (S -1 (x f )) > 0.
Hence, the function S -1 is strictly increasing. Moreover, item b) of Assumption 1 implies S -1 (0) = 0. We first consider the case where x f > 0. It holds that S -1 (x f ) > 0 for all x f > 0. Using (S -1 )(x f )

x f = (S -1 ) (c) as previously done in this proof, for some c such that c ∈ (0, x f ) and in view of items c) and d) of Assumption 1, S (S -1 (x f )) < S (S -1 (c)) for x f > c > 0, and thus (S -1 ) (

x f ) = 1 S (S -1 (x f )) > 1 S (S -1 (c)) = (S -1 ) (c). It holds that (S -1 )(x f ) x f = (S -1 ) (c) < (S -1 ) (x f ). Then, we deduce x f (S -1 ) (x f ) -(S -1 )(x f ) > 0. (30)
In view of ( 29) and (30), we have g (x f ) > 0 for all x f > 0, which means that

S -1 (x f ) x f increases with x f when x f > 0. Similar arguments show that S -1 (x f )
x f increases with |x f | when x f < 0. Therefore, the term

S -1 (x f ) x f
, as well as 2x 2 f , increases when 26) is strictly positive, we obtain the desired result.

x f ∈ (-∞, x - f old ) ∪ (x + f old , +∞). Under condition (14), -(S -1 ) (x f old ) + 2x 2 f + S -1 (x f ) x f is strictly posi- tive for x f ∈ (-∞, x - f old ) ∪ (x + f old , +∞). Recall that 1 -β+(S -1 ) (x f ) in (
4) Combining ( 13) and ( 20), we obtain for any

x f ∈ S(R) and β ∈ R, β * (x f , β) = β -β c + O(x 2 f
). This concludes the proof of the last item.

Proof of Proposition 2. Two steps are used to prove this proposition. We first show a boundedness property for system [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF]. Then we prove that (15) holds. Let ∆, M, T, µ > 0 and x f , x s , u be such that conditions (i) and (ii) of Proposition 2 hold. We rewrite [START_REF] Golubitsky | Singularities and Groups in Bifurcation Theory[END_REF] as follows

β = -kx 2 f β + kx 4 f + kx f (x s -u) . (31)
System (31) can be interpreted as a linear time-varying system with input kx 4 f + kx f (x s -u). We first consider the following nominal system β = -kx 2 f β.

(32)

Solutions to system (32) are given for t 0 by β(t) = Φ(t, 0) β(0), where Φ(t, s) := exp t s -kx 2 f (τ ) dτ and |Φ(t, 0)| 1 for all t 0. Hence the origin is stable for system (32). Let | β(0)| ∆ and consider the following Lyapunov function candidate for system (31)

V (t, β) := p(t) β2 , (33)
where p(t):= ∞ t Φ 2 (s, t)ds. We show below that V (t, β) is positive definite and radially unbounded with respect to β, uniformly in t. For 0 t s, we write

p(t) = ∞ t Φ 2 (s, t) ds = ∞ t exp s t -2kx 2 f (τ ) dτ ds. (34)
Since x f is (T, µ)-PE by assumption, for t s < t + T ,

µ t+T t x 2 f (τ )dτ = s t x 2 f (τ )dτ + t+T s x 2 f (τ )dτ . Hence, s t x 2 f (τ )dτ µ - t+T s x 2 f (τ )dτ. ( 35
)
Under assumption max( x f ∞ , x s ∞ , u ∞ ) M , we deduce from (34) and ( 35)

p(t) = ∞ t exp s t -2kx 2 f (τ ) dτ ds ∞ t exp -2kµ + 2k t+T s x 2 f (τ ) dτ ds exp(-2kµ) ∞ t exp 2kM 2 (t + T -s) ds exp(2kM 2 T ) 2kM 2 . ( 36
)
We also have from ( 34)

p(t) = ∞ t exp s t -2kx 2 f (τ ) dτ ds ∞ t exp -2kM 2 (s -t) ds 1 2kM 2 . (37)
Hence, V (t, β) is lower and upper bounded as

1 2kM 2 | β| 2 V (t, β) exp(2kM 2 T ) 2kM 2 | β| 2 . ( 38
)
The time-derivative of V (t, β) along the solution to (31) is, for any t 0,

∂V ∂t + ∂V ∂ β f ( β, x f , u -x s ) = ṗ(t) β2 + 2p(t) β -kx 2 f β + kx 4 f + kx f (x s -u) = ṗ(t) -2kx 2 f p(t) β2 + 2 βp(t) kx 4 f + kx f (x s -u) . (39)
Since Φ(t, t) = 1 and dΦ(s,t) dt = Φ(s, t)kx 2 f , we deduce after some calculations that

ṗ(t) -2kx 2 f p(t) = -1. (40)
Using (36) and (40)

, ∂V ∂t + ∂V ∂ β f ( β, x f , u -x s ) -| β| 2 + 2| β| exp(2kM 2 T ) C(M ) 2kM 2
. According to [24, Theorem 4.18], we conclude the boundedness property of ( 9) from (38) and

∂V ∂t + ∂V ∂ β f ( β, x f , u -x s ) -| β| 2 + 2| β| exp(2kM 2 T ) C(M ) 2kM 2
. That is, there exists η(∆, M, T ) > 0 such that solution β to system (31) holds | β(t)| η(∆, M, T ), for all t 0. To prove [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], let x f , x s , u be such that conditions (i) and (ii) of Proposition 2 holds, we define for t 0,

W (t, β1 , β2 ) := p 1 (t)( β1 -β2 ) 2 , (41)
where p 1 (t) := ∞ t Φ 2 1 (s, t) ds, and Φ 1 (t, 0) := exp

t 0 -kx 2 f 1 (τ )dτ . For | β1 (0)|, | β2 (0)| ∆, similar
computations in (36)-(37) shows, for any t 0

1 2kM 2 | β1 -β2 | 2 W (t, β1 , β2 ) exp(2kM 2 T ) 2kM 2 | β1 -β2 | 2 . ( 42
)
Following the similar computations for (40), we obtain ṗ1 (t)-2kx 2 f 1 p 1 (t) = -1. The following holds along to the solutions to (31) for any t 0, where the time arguments are omitted,

∂W ∂t + ∂W ∂ β1 f ( β1 , x f 1 , u 1 -x s1 ) + ∂W ∂ β2 f ( β2 , x f 2 , u 2 -x s2 ) = ṗ1 (t)( β1 -β2 ) 2 + 2p 1 (t)( β1 -β2 ) -kx 2 f 1 β1 + kx 4 f 1 +kx f 1 (x s1 -u 1 )) + kx 2 f 2 β2 -kx 4 f 2 + kx f 2 (x s2 -u 2 ) = ṗ1 (t) -2kx 2 f 1 p 1 (t) ( β1 -β2 ) 2 + 2kp 1 (t)( β1 -β2 ) × -β2 (x 2 f 1 -x 2 f 2 ) + (x 4 f 1 -x 4 f 2 ) + x s2 (x f 1 -x f 2 ) +x f 1 (x s1 -x s2 ) -x f 1 (u 1 -u 2 ) -u 2 (x f 1 -x f 2 ) . ( 43
)
Under condition (ii), ṗ1 (t) -2kx 2 f 1 p 1 (t) = -1, and the fact that x → x 4 and x → x 2 are locally Lipschitz, we derive that there exists ρ(∆, M, T ) > 0 such that

∂W ∂t + ∂W ∂ β1 f ( β1 , x f 1 , u 1 -x s1 )+ ∂W ∂ β2 f ( β2 , x f 2 , u 2 -x s2 ) -| β1 -β2 | 2 + ρ(∆, M, T )| β1 -β2 |(|x f 1 -x f 2 | + |x s1 - x s2 | + |u 1 -u 2 |
), we then obtain the desired result by following the similar analysis as in the first step of this proof.

Proof of Lemma 1. We first prove item 1). Let ∆ > δ > 0. For any β -β c < 0, any constant input δ < |u| < ∆, according to item 1 of Proposition 1, there exists ε such that for any ε ∈ (0, ε] and any

|x f (0)|, |x s (0)| ∆, it holds |(x f (t)-x * f ), (x s (t)-x * s )| c(∆, δ)e -σ(∆,δ)t |(x f (0) -x * f ), (x s (0) -x * s )
|, for t 0, where the fixed point (x * f , x * s ) is different from (0, 0) (because u = 0) and c(∆, δ), σ(∆, δ) > 0. Hence, condition (ii) of item 1) of Lemma 1 holds. We also deduce that there exists T

1 (∆, δ) > 0 such that |x f (t) -x * f | |x * f | 2
for any t T 1 (∆, δ). It holds that for t T 1 (∆, δ),

|x f (t)| = |x f (t) -x * f + x * f | |x * f | -|x f (t) -x * f | |x * f | - |x * f | 2 = |x * f | 2 . Let T (∆, δ) > T 1 (∆, δ), we deduce, for any t 0 t+T (∆,δ) t x 2 f (τ )dτ = t+T1(∆,δ) t x 2 f (τ )dτ + t+T (∆,δ) t+T1(∆,δ) x 2 f (τ )dτ t+T (∆,δ) t+T1(∆,δ) |x f (τ )| 2 dτ |x * f | 2 2 (T (∆, δ) -T 1 (∆, δ)).
Thus, condition (i) of item 1) of Lemma 1 holds with

µ(∆, δ) = |x * f | 2 2 × (T (∆, δ) -T 1 (∆, δ)).
We next prove item 2). Let ∆ > δ > 0, for any 0 < β -β c < 1 and for any input u ∈ (-ū, ū), where ū ∈ (0, ∆], according to item 2) of Proposition 1, there exists ε such that for all ε ∈ (0, ε], and

|x f (0)|, |x s (0)| ∆ and |x f (0) -p * f |, |x s (0) -p * s |
δ, all trajectories of system (2) converge to an asymptotic periodic orbit P ε with period T ε and it is locally exponentially stable. Thus solutions x f , x s to system (2) are bounded. Therefore, condition (ii) of item 2) of Lemma 1 holds. Moreover, let

|x lc f | be the absolute value of x f -component of P ε , which is T ε -periodic. The average value of |x lc f | over a period is denoted by xf = 1 T ε t+T ε t |x lc f (τ )
|dτ for all t 0 and it is strictly positive. Due to the globally attractivity of P ε , for any η > 0 there exists T

1 (∆, δ, η) > 0 such that |x f (t) -xf | η for t T 1 (∆, δ, η). Let T (∆, δ) = T 1 (∆, δ, η) + T ε , the following holds t+T (∆,δ) t x 2 f (τ )dτ= t+T1(∆,δ,η) t x 2 f (τ )dτ + t+T (∆,δ) t+T1(∆,δ,η) x 2 f (τ )dτ t+T (∆,δ) t+T1(∆,δ) x2 f dτ + t+T (∆,δ) t+T1(∆,δ) (x 2 f (τ ) -x2 f )dτ x2 f T ε - t+T (∆,δ) t+T1(∆,δ) |x 2 f (τ ) -x2 f |dτ. (44) Since x → x 2 is locally Lipschitz, it holds |x 2 f (t) -x2 f | ψ(∆, δ)|x f (t) -xf |. We deduce from (44) t+T (∆,δ) t x 2 f (τ )dτ x2 f T ε - t+T (∆,δ) t+T1(∆,δ) |x 2 f (τ ) -x2 f |dτ x2 f T ε -ψ(∆, δ)ηT ε . (45)
Let choose T 1 (∆, δ, η) sufficiently large such that η = xf 2ψ(∆,δ) . Then, the above inequality follows

t+T (∆,δ) t x 2 f (τ )dτ x2 f T ε -ψ(∆, δ)ηT ε = x2 f 2 T ε .
Thus, condition (i) of item 2) of Lemma 1 holds with

µ(∆, δ) = x2 f 2 T ε .
Proof of Theorem 2. In the following, we denote the (x f , x s )-component of the solutions to (16) as x. We first prove item 1). Let ∆ > δ > 0, β < β c and constant input δ < |u| < ∆, according to item 1) of Proposition 1, there exists ε > 0 such that for any ε ∈ (0, ε] and any |x(0)| ∆, the corresponding solution x to subsystem (16a)-(16b) satisfies for t 0,

|x(t) -x * | c 1 (∆, δ)e -c2(∆,δ)t |x(0) -x * |, (46)
where x * := (x * f , x * s ) is the fixed point of subsystem (16a)-(16b), and c 1 (∆, δ), c 2 (∆, δ) > 0.

We next consider for t 0,

W (t, β) := p(t)( β -β) 2 , ( 47
)
where β is the equilibrium point of (16c) associated with constant input (x * , u), and p(t):= ∞ t Φ 2 (s, t) ds as in the proof of Proposition 2. For any | β(0)| ∆, by the similar computations in (36)-(37) for any t 0

M (∆, δ)| β -β| 2 W (t, β) M (∆, δ)| β -β| 2 , (48)
where M (∆, δ), M (∆, δ) > 0. Following similar lines as in (40), we obtain ṗ(t) -2kx 2 f p(t) = -1, and we derive that for any t 0, where the time arguments are omitted,

∂W ∂t + ∂W ∂ β f ( β, x f , u -x s ) -| β -β| 2 + ν(∆, δ)| β -β| |x f -x * f | + |x s -x * s | , (49)
where ν(∆, δ) > 0. We deduce from the above inequality

|β(t) -β| 2 c 3 (∆, δ)e -c4(∆,δ)t | β(0) -β| 2 +c 5 (∆, δ)e -c4(∆,δ)t |x(0) -x * | 2 , (50)
with c 3 (∆, δ), c 4 (∆, δ), c 5 (∆, δ) > 0. This property together with (46) imply that the semiglobal exponential stability of system [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF], which, in turn, implies the global asymptotic stability, according to [START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF]Proposition 3.4]. Thus, the desired result follows.

We next prove item 2). Let ∆ > 0. For any 0 < ββ c < 1, according to item 2) of Proposition 1, there exist ε > 0, ū ∈ (0, ∆) such that for any ε ∈ (0, ε], for any u ∈ (-ū, ū), there exists m > 0 such that for any |x(0)| P ε m except the unique unstable fixed point of (16a)-(16b), the corresponding solution x to subsystem (15a)-(15b) satisfies for t 0,

|x(t)| P ε c 6 (m)e -c7(m)t |x(0)| P ε , (51)
where P ε is the periodic solution of subsystem (16a)-(16b) and c 6 (m), c 7 (m) > 0. By virtue of Lemma 2 and Proposition 2, subsystem (16c) is semiglobally incrementally ISS. Then according to [START_REF] Angeli | A Lyapunov approach to incremental stability properties[END_REF]Proposition 4.4], there exists initial condition βlc (0) such that βlc (t) is a periodic solution to (16c) with periodic input (P ε , u). We choose W (t, β) := p(t)( ββlc ) 2 . Then using similar arguments as above we obtain the local exponential stability of system (16) with respect to a limit cycle. Proof of Theorem 3. We start by proving item 1). By Theorem 2, if β < β c then (x f (t), x s (t)) → (x * f , x * s ) and β(t) → β as time goes to infinity. Note that, since (x * f , x * s ) is a steady-state of (16a)-(16b), it must satisfies (x * f , x * s ) ∈ S 0 . Moreover, u = 0 implies x * f = 0. Because -kx * f (-x * 3 f + βx * f + u -x * s ) = 0 and (x * f , x * s ) ∈ S 0 , it follows that β =

x * 3 f -S -1 (x * f )+βx * f x * f = β * (x * f , β).
We next prove item 2). At the singular limit = 0, the one dimensional estimate critical manifold R 0 of system ( 16) is defined by the following two equations 0 = -x f + S(βx f + u -x s ), 0 = -x 3 f + βx f + u -x s , which can explicitly be solved as R 0 = (x f , x s , β) :

(53)

x s = φ(x f , u, β) := -S -1 (x f ) + βx f + u,

β = x 3 f -u + φ(x f , u, β) x f = β * (x f , β) .
Note that the (x f , x s )-projection of R 0 is S 0 . The layer dynamics reads

     ẋf = -x f + S(βx f + u -x s ), ẋs = 0, β = -kx f (-x 3 f + βx f + u -x s ). (54a) (54b) (54c)
To study the normal hyperbolicity of R 0 , we compute the Jacobian matrix of (54) on R 0 as for x s = φ(x f , u, β), β = β * (x f , β). The two nonzero eigenvalues are λ 1 = -kx 2 f , λ 2 = -1 + βS (S -1 (x f )). For 0 < β -β c < 1, it holds that β > β c = 1 S (0) , which is equivalent to 1 β < S (0) and therefore, invoking item d) of Assumption 1, the equation -1 + βS (S -1 (x f )) = 0 has exactly one positive and negative roots, x + f old and x - f old , respectively. Moreover,

J=     -1 + βS (βx f + u -x s ) -S (βx f + u -x s ) 0 0 0 0 4kx 3 f -2k βx f -k(u -x s ) kx f -kx 2 f    
λ 2 < 0, if x f < x - f old or x f > x + f old , > 0, if x - f old < x f < x + f old .
Since λ 1 < 0 for x f = 0, it follows that two branches of R 0 given by

R 0 + := R 0 ∩ {x f > x + f old }, R 0 -:= R 0 ∩ {x f < x - f old },
are locally exponentially attractive branches of the estimate critical manifold R 0 . Let β+ f old = β * (x + f old , β), βf old = β * (x - f old , β),

x + s,f old = φ(x + f old , u, β), x - s,f old = φ(x - f old , u, β). The end points of R 0 + and R 0 -are given by F + := (x + f old , β+ f old , x + s,f old ), F -:= (x - f old , βf old , x - s,f old ) respectively, are fold points of the estimate critical manifold R 0 . See Figure 7.

The slow flow on R 0 is defined by the reduced dynamics

     0 = -x f + S(βx f + u -x s ), x s = x f -x s , 0 = -kx f (-x 3 f + βx f + u -x s ). (55a) (55b) (55c)
The flow on R 0 + , as well as on R 0 -, with respect to the reduced dynamics of (55) is given as follows. It holds that x s < 0(x f < x s < 0) on R 0 -and, similarly, x s > 0 on R 0 + . Noticing that the x s -projection of R 0 + (resp. R 0 -) is the semiline (-∞, x + s,f old ) (resp. (x - s,f old , ∞)), it follows that all trajectories on R 0 + eventually reach the fold point F + . Conversely, all trajectories on R 0 -eventually reach F -(see Fig. 7).

At the folds, we connect the slow flow with the (twodimensional) fast flow. We claim that the fast flow brings the trajectory on the opposite branch of the critical manifold. Indeed, by monotonicity of (54a), x f converges to equilibrium. By the cascade structure of (54a)-(54c) and incremental ISS property of (54c), β also converges to equilibrium. Because the only equilibria of (54a)-(54c) are on the estimate critical manifold, the result follows.

We have therefore constructed a candidate singular periodic orbit Q 0 . The persistence of this singular orbit for ε > 0 follows as in the proof of Proposition 1 and is omitted. Let Q ε be the resulting periodic orbit of [START_REF] Tang | Parameter estimation without fitting: a qualitative approach[END_REF]. By Theorem 2, almost all the trajectories of system ( 16) asymptotically converge to Q ε . Following again the same arguments as the proof of Proposition 1, it holds that trajectories along Q ε spend only an O(ε)-fraction of the limit cycle period outside an O(ε) neighborhood of the estimate critical manifold. Because β(t) = β * (x lc f (t + θ), β) on the estimate critical manifold, the result follows.

Assumption 1

 1 The following properties holds. a) S is smooth on R. b) S(0) = 0. c) S (a) > 0 for all a ∈ R (monotonicity), and argmax a∈R S (a) = 0 (sector-valued). d) sgn (S (a)) = -sgn(a) for all a ∈ R. 2 Standard sigmoid functions such as a → c 1 tanh(c 2 a), a → c1 1+e -c 2 a -c1 2 , a → e -e -c 1 a -c 2 with c 1 , c 2 ∈ R >0 , verify Assumption 1. Moreover, when S = n i=1 S i where S 1 , • • • , S i with (i ∈ N) are sigmoid functions satisfying Assumption 1, then so does S.

Fig. 1 .

 1 Fig. 1. Qualitative shape of the critical manifold S 0 .

2 S 1 : 2 ( 1 Fig. 2 .Fig. 3 .Fig. 4 . 6 S 1 : 2 Fig. 5 .Fig. 6 .

 212123461256 Fig. 2. States x f , xs of system (2) for different nonlinearities.

Due to [ 21 ,

 21 Proposition 4.5], for any initial condition x(0) except the unique unstable fixed point x * of subsystem (16a)-(16b), the corresponding solution to system (16c) holds lim t→∞ | β(t)| βlc = 0 as lim t→∞ |x(t)| P ε = 0 . In addition, local exponential stability implies local stability, then the desired result follows.

  Fig. 7.

We adopt the notion of almost global asymptotical stability defined in[START_REF] Angeli | An almost global notion of input-tostate stability[END_REF]