P Loiseau

D Caurant

N Baffier

K Dardenne

J Rothe

M Denecke

S Mangold

Structural characterization of Nd-doped calcium aluminosilicate glasses designed for the preparation of zirconolite (CaZrTi 2 O 7 ) -based glass-ceramic

This work concerns glasses belonging to SiO 2 -Al 2 O 3 -CaO-ZrO 2 -TiO 2 -Nd 2 O 3 system that lead to zirconolite crystallization in their bulk after nucleation + crystal growth thermal treatments and that could find application as nuclear waste form. The understanding of crystallization processes in glasses implied to investigate their structure. The environment around Ti, Zr (nucleating agents) and Nd was characterized for various Nd 2 O 3 loadings. Electron spin resonance study of the small amount of Ti 3+ occurring in glasses enabled to identify two types of sites for titanium. EXAFS showed that Zr occupied a quite well defined 6-7-fold site coordinated by oxygen at 2.20 Å, with second neighbors that could correspond to Ca/Ti and Zr around 3.5 Å. This short range order presents some similarities with zirconolite, which could predispose glass to zirconolite nucleation. Nd environment was probed by optical spectroscopies, ESR and EXAFS. Results showed that the environment around Nd was very constrained by the glassy network and significantly differed from the one in zirconolite. Nd occupied a highly distorted 8-9-fold coordinated site in glass with oxygen atoms at 2.53 Å. No second neighbors were clearly identified by EXAFS around Nd, but the study of Nd optical fluorescence decays suggested a strong interaction between Nd ions.

INTRODUCTION

Zirconolite (CaZrTi 2 O 7 ) is a naturally occurring phase which exhibits an excellent chemical durability and a good actinide containment capacity. This is a waste form of choice for the specific immobilisation of minor actinides or plutonium. A simpler alternative to the ceramic route is to prepare zirconolite crystals by controlled crystallization (nucleation+growth) of a parent glass. In addition, such glass-ceramic offers a greater chemical flexibility against waste composition fluctuations due to the presence of a residual glass surrounding zirconolite crystals. The basic parent glass composition which was studied is (mol. %): 48.8 SiO 2 , 8.5 Al 2 O 3 , 25.3 CaO, 11.3 TiO 2 , 5.0 ZrO 2 , 1.1 Na 2 O 1 . Various loadings of Nd 2 O 3 (trivalent actinide surrogate) were added up to 2.2 mol. % (10 wt. %). Such glasses led to the only crystallization of zirconolite by devitrification (nucleation+crystal growth) 2 . To get insight into crystallization processes, the environment around three elements was investigated: Ti and Zr (main components of zirconolite playing a key role in nucleation processes), and Nd. EXPERIMENTAL 40g batch glasses were prepared following a standard method comprising two stages of melting and casting at 1550°C with intermediate grinding. A similar glass transformation temperature around 760 °C was found for all the glasses, so that the addition of Nd 2 O 3 up to 2.2 mol. % has no significant effect on the mean bond strength of the glassy network. Ti 3+ Electron Spin Resonance (ESR) measurements were performed at 120 K using a Bruker ESP 300e spectrometer operating at X-band. Room temperature EXAFS (Extended X-ray Absorption Fine Structure) measurements were performed at Zr-K edge (17998 eV) and Nd-L 3 edge (6208 eV) at ANKA synchrotron (Germany, Karlsruhe) in transmission mode through quenched glass samples containing 1.3 mol. % of Nd 2 O 3 (6 wt. %). EXAFS signals were simulated using the UWXAFS package 3 . Nd 3+ optical properties were characterised at 15 K by absorption from 400 to 900 nm and by emission from 860 to 960 nm ( 4 F 3/2  4 I 9/2 channel). At 15 K, only the lowest energy Kramers doublet of the 4 I 9/2 or of the 4 F 3/2 multiplet is populated for absorption or emission measurements respectively.

RESULTS AND DISCUSSION

ESR study of native Ti 3+ ions in parent glass

During glass meting, a Ti 4+  Ti 3+ equilibrium was established and then was frozen by quenching at room temperature. Therefore, a small quantity of Ti 3+ formed. Their amount was quantified by ESR at 120 K (Fig. 1) with a DPPH standard: the proportion of Ti at the trivalent oxidation state was below 0.7 ppm (case of the undoped glass) and decreased with increasing Nd 2 O 3 concentration (no Ti 3+ ion detected for the 10 wt. % Nd 2 O 3 glass). This evolution could result from an increase of the glass basicity 4 with Nd concentration, if we admit that Nd acts as a network modifier. Ti 3+ ESR signal was simulated with the assumption of an axial symmetry for simplicity reasons. Without any distribution of g factor, it was necessary to consider two contributions (with gaussian line shape) to simulate the signal and the shoulder occurring around 344 mT (Fig. 1 and Table I). This gives evidence of two types of sites for Ti 3+ ions in the glass, referred as site 1 and site 2. Site 2 contributes only to 8.5 % of the overall signal. The g factors give pieces of information about the symmetry of Ti 3+ environment in the glass [START_REF] Wertz | Electron Spin Resonance: Elementary theory and practical applications[END_REF] . As Ti 3+ ESR signal can be detected at room temperature, several simple geometries can be excluded: tetrahedron, octahedron and tetragonally elongated octahedron. Moreover, a four-fold coordination is quite unlikely for a cation as big as Ti 3+ (r VI (Ti 3+ )=0.67 Å> r VI (Ti 4+ )=0.605 Å 6 ). However, ESR results can be fully explained by considering a compressed octahedral symmetry, highly distorted to have the order g // < g  (more distorted for site 2). This geometry could correspond either to a C 4v or a D 4h symmetry (5 to 6 coordination number). It is interesting to notice that the C 4v symmetry is the one of the square pyramidal O=TiO 4 frequently proposed for Ti 4+ environment in glasses 7 .

EXAFS characterisation of Zr 4+ environment in parent glass

The short range order around Zr in the 6 wt. % Nd 2 O 3 parent glass was determined from the EXAFS spectrum (Fig. 2). Zr occupies a well defined site (debye-waller factor  2 = 0.007  0.001 Å -2 ) of 6.5  0.7 coordination number with a mean Zr-O distance d(Zr-O) = 2.15  0.02 Å. This environment is relatively close to the one in a depolymerised glass such as R7T7 glass, with a coordination number around 6 and a mean Zr-O distance ranging from 2.07 to 2.10 Å 8 . The slightly higher coordination number and Zr-O distance in our glass could be due to the higher polymerisation of the glassy network. For 100 moles of undoped glass, assuming that Al, Ti and Zr adopt respectively the geometries (AlO 4 ) -, (O=TiO 4 ) 2-and (ZrO 6 ) 2-, 24.7 moles of (CaO+Na 2 O) would be mobilised for their charge compensation, so that only 1.7 moles of (CaO+Na 2 O) would be free to depolymerize the glassy network. Then, oxygens are less polarisable and a higher coordination number is required for Zr 4+ charge shielding. In zirconolite, Zr coordination number is 7 and d(Zr-O) = 2.20(1) Å. This short range order similarity with the parent glass could predispose to zirconolite nucleation. The EXAFS analysis of the next nearest neighbours around Zr is complicated (Fig. 2): at least two shells (likely three) with overlapping shapes are necessary for fitting, and they partly interfere destructively. Nevertheless, the Fourier peak around 3.2 Å seems to be associated to Ca/Ti on the one hand and Zr on the other hand, both occurring around 3.5 Å from Zr.

Spectroscopic investigation of Nd 3+ environment in parent glass

The fluorescence spectra corresponding to the 4 F 3/2  4 I 9/2 Nd transition were recorded at 15 K by exciting the Nd-doped glasses around 806.5 nm. These spectra allow to estimate the total splitting of the 4 I 9/2 ground multiplet which is about 500 cm -1 , against 775 cm -1 for Nddoped zirconolite samples. The 4 I 9/2 splitting is a measure of the mean crystal field experienced by Nd 3+ ion in a host 9 . The greater crystal field of Nd in zirconolite as compared to parent glass indicates not only that the environments are different between the two hosts, but also that the mean Nd-O distance is shorter in zirconolite.

Nd optical absorption spectra were recorded for the Nd-doped glasses at 15 K. They are all very similar, even though a slight broadening of absorption lines occurs with increasing Nd 2 O 3 concentration from 0.5 to 10 wt. %. Three absorption bands from 4 I 9/2 state are particularly interesting: towards 2 P 1/2 ( 23200 cm -1 ), 4 F 3/2 ( 11400 cm -1 ) and 2 G 7/2 -4 G 5/2 ( 17200 cm -1 ). The first two bands ( 2 P 1/2 and 4 F 3/2 ) notably give information about the variety of Nd environment because of the low degeneracy of the excited state, whereas the latter band ( 2 G 7/2 -4 G 5/2 ) is described as hypersensitive and is very sensitive to Nd environment. For instance, the FWHM of the 4 I 9/2  2 P 1/2 band is approximately 150 cm -1 in Nd-doped parent glasses, much more than the one measured for zirconolite (30 cm -1 ), which reflects a higher distribution of environment in the glasses. Moreover, by comparison with literature 10 , the shape of the 4 I 9/2  2 G 7/2 -4 G 5/2 is structure-less and is characteristic either of highly polymerised glassy network or of glasses containing high field strength network modifiers. All these observations agree with an irregular Nd environment partly imposed by the glassy network topology, contrarily to Zr.

The environment of Nd was probed by EXAFS at L 3 edge (Fig. 3). Nd occupies a highly distorted site: the EXAFS simulation implies a cumulant analysis of the data (non-gaussian asymmetric Nd-O pair distribution 11 ) with a large  2 = 0.029  0.001 Å -2 . Nd-O coordination number is 8.6  0.8 and the mean Nd-O distance is about 2.53  0.02 Å in the glass. In zirconolite, Nd coordination by O is 8 and, by EXAFS, Nd-O bond distance was found to be 2.440.01 Å with  2 =0.0100.001 Å -2 . So, the environment around neodymium in the glass shows striking differences with the one in zirconolite and is much more distorted. Therefore, the incorporation of neodymium into zirconolite crystals during devitrification will imply strong reconstructive phenomena which are assumed to occur at high temperature, when the viscosity of the supercooled liquid is lowered. In these conditions, neodymium can hinder zirconolite nucleation with increasing diffusionnal difficulties.

CONCLUSION

ESR study of Ti 3+ ions formed during glass melting revealed that their environment could correspond to a square pyramidal O=TiO 4 or to a six-fold coordinated octahedron axially compressed. The coordination numbers involved in such geometries are close to the ones encountered in zirconolite (6 and 5) which is aimed to crystallize by devitrification. Zr-K edge EXAFS spectra showed that Zr 4+ ions occupied well defined 6-7 coordinated sites, close to the ones of zirconolite, able to explain the nucleating effect of Zr in the glass. Contrarily to Zr 4+ ion which can be described as a cross-linking agent that reinforces the glassy network, Nd 3+ ion must rather be considered as a network modifier or a charge compensator which exhibits a lower ability to impose its own environment. All the optical spectroscopies and EXAFS agree to describe Nd environment as being constrained by the glassy network with coordination number and Nd-O distance significantly greater than the ones of zirconolite: Nd 3+ ion cannot act as a nucleating agent for zirconolite crystallization.

Fig. 1 .

 1 Fig. 1. Experimental and simulated Ti 3+ ESR signal for the undoped glass

Fig. 2 .

 2 Fig. 2. Fit (solid lines) to the EXAFS (open circles) of the 6 wt. % Nd 2 O 3 parent glass (Zr-K edge, room temperature).

Fig. 3 .

 3 Fig. 3. Fit (solid lines) to the EXAFS (open circles) of the 6 wt. % Nd 2 O 3 parent glass (Nd-L 3 edge, room temperature).

TABLE I .

 I Results of the ESR simulation.

		g 	FWHM	g 	FWHM
	Site 1	1.946 6.5 mT	1.88	12 mT
	Site 2	1.973	4 mT	1.89	4 mT


Miguel Navascués

Frantz Depaulis

Brent C Emerson
email: b.emerson@uea.ac.uk

Combining contemporary and ancient DNA in population genetic and phylogeographic studies

Keywords: phylogeography, serial coalescence, heterochronous data, ancient DNA approximate Bayesian computation

The analysis of ancient DNA in a population genetic or phylogeographic framework is an emerging field, as traditional analytical tools were largely developed for the purpose of analyzing data sampled from a single time point. The analysis of heterochronous sequence data from closed panmictic populations has received attention with Markov chain Monte Carlo (MCMC) approaches, but attributing genetic differences between temporal samples to mutational events between time points requires the consideration of other factors that may also result in genetic differentiation. Geographic effects are an obvious factor for species exhibiting geographic structuring of genetic variation, and departures such as this from a closed panmictic model require researchers to either exploit software developed for the analysis of isochronous data, take advantage of simulation approaches using algorithms developed for heterochronous data, or explore approximate Bayesian computation. Here we review statistical approaches employed and available software for the joint analysis of ancient and modern DNA, and where appropriate we suggest how these may be further developed.

Introduction

Non-contemporaneous, or ancient DNA, is providing biologists with new and exciting opportunities to investigate evolutionary pattern and process over a range of temporal scales, from decades (e.g. [START_REF] Harper | Analysis of museum specimens suggests extreme genetic drift in the adonis blue butterfly (Polyommatus bellargus)[END_REF][START_REF] Martinez-Cruz | Population fragmentation leads to spatial and temporal genetic structure in the endangered Spanish imperial eagle[END_REF] to hundreds of thousands of years (e.g. [START_REF] Willerslev | Ancient biomolecules from deep ice cores reveal a forested Southern Greenland[END_REF]. For the purposes of this paper we define ancient DNA (hereafter aDNA) as DNA recovered from non-ideal biological material -that is to say material that was not preserved or maintained in a manner typically associated with downstream DNA analysis, for which the host organism is no longer alive. Such material includes subfossil remains (typically bones and teeth), archaeological remains, coproliths, mummies, naturally (i.e. not laboratory) frozen remains, ice cores, sediments, museum and herbarium tissues. DNA extracted from such material is typically low quality, with constraints upon both the amount and integrity of the DNA that can be obtained. Nevertheless, there is an increasing volume of studies successfully obtaining samples of aDNA, and frequently these data are being analyzed in a temporal and geographic context (see [START_REF] Ramakrishnan | Using phylochronology to reveal cryptic population histories: review and synthesis of 29 ancient DNA studies[END_REF] for a review). Thus aDNA is becoming increasingly accessible for both population genetic and phylogeographic analysis, offering the possibility of using temporal samples of DNA to characterize population history. However, the analysis of such data is an emerging field, as traditional population genetic and phylogeographic tools were largely developed for the purpose of analyzing data sampled from a single time point, or at the most a sampling interval that spanned no more than a few generations. In this paper we review some of the analytical approaches employed for the joint analysis of ancient and modern DNA, and offer some suggestions for future directions. While we focus on studies utilizing aDNA for temporal sampling, it is important to point out that temporal sampling is not restricted to aDNA. Organisms with short generation times such as viruses offer similar opportunities (e.g. [START_REF] Norja | Rapid sequence change and geographical spread of human parvovirus B19: Comparison of B19 virus evolution in acute and persistent infections[END_REF][START_REF] Shackleton | JC evolution and its association with human populations[END_REF]), but without constraints on DNA quality. Perhaps unsurprisingly, the most widely employed marker for phylogeographic and population genetic analyses with temporal sampling over more distant time scales (from hundreds, to thousands of years) is the mitochondrial genome, in particular the fast evolving control region (e.g. [START_REF] Barnett | Phylogeography of lions (Panthera leo ssp.) reveals three distinct taxa and a late Pleistocene reduction in genetic diversity[END_REF][START_REF] Lambert | Rates of evolution in ancient DNA from Adelie penguins[END_REF][START_REF] Shapiro | Rise and fall of the beringian steppe bison[END_REF][START_REF] Valdiosera | Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears[END_REF][START_REF] Valdiosera | Surprising migration and population size dynamics in ancient Iberian brown bears (Ursus arctos)[END_REF]. Given recent evidence that mtDNA is not a reliable marker for demography [START_REF] Bazin | Population size does not influence mitochondrial genetic diversity in animals[END_REF], nuclear data would be a welcome addition, and next generation sequencing technologies offer some promise in this direction. Over more recent time scales of decades, microsatellite markers gain greater prominence, where their high variability provide more potential to reveal demographic events over narrower time scales (e.g. [START_REF] Harper | Analysis of museum specimens suggests extreme genetic drift in the adonis blue butterfly (Polyommatus bellargus)[END_REF][START_REF] Martinez-Cruz | Population fragmentation leads to spatial and temporal genetic structure in the endangered Spanish imperial eagle[END_REF]). However, the utility of microsatellites is not restricted to decadal analyses, and recent efforts have seen the combination of nuclear microsatellites with mtDNA to investigate demographic changes spanning several thousands of years [START_REF] Keyser-Tracqui | Population origins in Mongolia: Genetic structure analysis of ancient and modern DNA[END_REF]. Single nucleotide polymorphism (SNP) data offers similar potential and has been used to assess genetic change over a 4,000 year period for cattle [START_REF] Svensson | Tracing genetic change over time using nuclear SNPs in ancient and modern cattle[END_REF]). As in classical population genetics, the choice of genetic marker for studies incorporating aDNA should be appropriate for the question(s) to be addressed, and in common to all genetic markers are the fundamental population genetic processes that may shape differentiation between heterochronous samples. Population samples of DNA from different time points may vary from subtle changes in allele frequency through to allele loss or gain, due to the processes of mutation, selection, genetic drift and migration. Recent interest has focused on sampling DNA sequences from "measurably evolving" populations [START_REF] Drummond | Measurably evolving populations[END_REF] for which there are sufficiently long or numerous sampled sequences and a fast mutation rate relative to the available range of sequence sampling times. Sequences from such populations have the potential to enable the analysis of temporal changes in the size, structure and substitution rates of populations. Problems have been pointed out about the reliability of aDNA for the estimation of population history, as postmortem DNA damage may act to inflate diversity estimates [START_REF] Axelsson | The effect of ancient DNA damage on inferences of demographic histories[END_REF]. Distinguishing what is a genuine mutation and what is the consequence of postmortem degradation is thus important for the reliable estimation of population history. A conservative approach is to remove all sites that might represent false, or post mortem induced degradation. Applying this to a data set of steppe bison (Bison bison/Bison priscus) previously analyzed by [START_REF] Shapiro | Rise and fall of the beringian steppe bison[END_REF] results in the removal of signal for population size change observed in the original analyses [START_REF] Axelsson | The effect of ancient DNA damage on inferences of demographic histories[END_REF]). However, it has been argued that this approach is too conservative, and a subsequent reanalysis implementing a model where sequence variation is the result of a joint process of mutation and postmortem DNA damage is consistent with previous conclusions of population size change over time [START_REF] Rambaut | Accommodating the effect of ancient DNA damage on inferences of demographic histories[END_REF]. Clearly there is a need to accommodate mutational artifacts when analyzing aDNA in a population genetic and phylogeographic framework, and recent methods developed to predict errors due to postmortem degradation are a welcome addition [START_REF] Mateiu | Bayesian inference of errors in ancient DNA caused by postmortem degradation[END_REF]. Phylogeographic analysis of temporally sampled DNA sequences provides for the direct quantification of population turnover within species, with particular reference to climatically mediated regional extinction and recolonization [START_REF] Benton | How did life become so diverse? The dynamics of diversification according to the fossil record and molecular phylogenetics[END_REF][START_REF] Stewart | Refugia revisited: individualistic responses of species in space and time[END_REF]. This is perhaps the most powerful contribution of aDNA to the field of phylogeography, with recent studies suggesting regional population extinction and subsequent recolonization [START_REF] Barnes | Dynamics of Pleistocene population extinctions in Beringian brown bears[END_REF][START_REF] Barnett | Phylogeography of lions (Panthera leo ssp.) reveals three distinct taxa and a late Pleistocene reduction in genetic diversity[END_REF][START_REF] Hofreiter | Sudden replacement of cave bear mitochondrial DNA in the late Pleistocene[END_REF][START_REF] Leonard | Megafaunal extinctions and the disappearance of a specialized wolf ecomorph[END_REF]. Temporally sampled DNA sequences also provide the potential for the estimation of substitution rates and divergence times without paleontological calibrations, and there has been much recent interest in the idea that substitution rates within species may be much higher than previously thought [START_REF] Ho | Elevated substitution rates estimated from ancient DNA sequences[END_REF][START_REF] Ho | Molecular clocks: when times are a-changin[END_REF][START_REF] Ho | Time dependency of molecular rate estimates and systematic overestimation of recent divergence times[END_REF][START_REF] Ho | The effect of inappropriate calibration: Three case studies in molecular ecology[END_REF][START_REF] Penny | Evolutionary biology -relativity for molecular clocks[END_REF]. Bayesian Markov chain Monte Carlo (MCMC) analyses of ancient data sets has consistently generated substitution rate estimates exceeding those from the literature. While this acceleration of molecular rates may partially be explained by evolutionary processes [START_REF] Ho | Time dependency of molecular rate estimates and systematic overestimation of recent divergence times[END_REF]) it has been shown that such rate estimates can be biased due to the low information content of aDNA data [START_REF] Debruyne | Time dependency of molecular rates in ancient DNA data sets, a sampling artifact[END_REF] or demographic model misspecification [START_REF] Navascués | Elevated substitution rate estimates from ancient DNA: model violation and bias of Bayesian methods[END_REF]. These caveats should be taken into account when estimating molecular rates and in their use to date historical events. The concerns raised by [START_REF] Axelsson | The effect of ancient DNA damage on inferences of demographic histories[END_REF] and the results of [START_REF] Navascués | Elevated substitution rate estimates from ancient DNA: model violation and bias of Bayesian methods[END_REF] highlight the need for caution when analyzing aDNA in a population genetic or phylogeographic context. Both studies reveal that variation among DNA sequences that originates from processes not accounted for within an underlying evolutionary model may lead to incorrect conclusions regarding evolutionary history. [START_REF] Axelsson | The effect of ancient DNA damage on inferences of demographic histories[END_REF] demonstrate that when nucleotide differences between temporal samples arise from postmortem DNA damage they can misleadingly contribute to demographic inference. [START_REF] Navascués | Elevated substitution rate estimates from ancient DNA: model violation and bias of Bayesian methods[END_REF] demonstrate that when nucleotide differences between temporal samples arise not from mutation, but from other population genetic processes such as genetic drift, immigration or selection, they can bias divergence and substitution rate estimation. The take home message is model violation should be given due consideration as a potential explanation for the result of a model-based analysis. This is not something unique to population genetic/phylogeographic analyses incorporating aDNA (e.g. [START_REF] Becquet | Learning about modes of speciation by computational approaches[END_REF][START_REF] Strasburg | How robust are "Isolation with Migration" analyses to violations of the IM model? A simulation study[END_REF]. New DNA sequencing technologies promise to deliver greater amounts of aDNA sequence information for a greater number of taxa (e.g. [START_REF] Allentoft | Identification of microsatellites from an extinct moa species using high-throughput (454) sequence data[END_REF][START_REF] Briggs | Targeted Retrieval and Analysis of Five Neandertal mtDNA Genomes[END_REF]). While the majority of studies to date have relied on mtDNA sequence variation, nuclear microsatellites and nuclear coding sequence variation is also within reach (e.g. [START_REF] Keyser-Tracqui | Population origins in Mongolia: Genetic structure analysis of ancient and modern DNA[END_REF][START_REF] Krause | The derived FOXP2 variant of modern humans was shared with neandertals[END_REF][START_REF] Lalueza-Fox | Bitter taste perception in Neanderthals through the analysis of the TAS2R38 gene[END_REF][START_REF] Lalueza-Fox | Genetic characterization of the ABO blood group in Neandertals[END_REF]. It is thus both timely and appropriate to evaluate how best to combine ancient DNA with contemporary DNA in a population genetic and phylogeographic context. What follows is a review of the statistical tools available, and suggestions for their future development.

Testing for heterochrony

For the analysis of heterochronous data, one may, as a first step, test if the data show measurable evolution. If they do not, data may be pooled to proceed with further analyses using available tools for isochronous data. PATH-O-GENE (Table 1), is a simple tool that may be used with an input tree, typically a phylogeny reconstructed without assuming a molecular clock, to compute the correlation between tip to root distances and sampling times. This provides associated estimates of the mutation rate and a dated root (the slope and abscissa intercept of the regression) and an output tree assuming a molecular clock, but taking into account sampling times. It should however be noted that rigorous testing of such a correlation would require an additional randomization procedure, since the data do not obey parametric assumptions.

To take into account the potentially confounding effects of factors such as geography (genetic differentiation resulting from spatial rather than temporal separation) it is possible to partition genetic variation through the AMOVA [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data[END_REF] implemented in ARLEQUIN (Table 1). Although originally intended to assess geographic population structure, AMOVA can also be applied to heterochronous data partitioned into sampling times, or time bins when pooling is appropriate. This analysis tests the sampling time effect and estimates the proportion of genetic variation explained by sampling times. Formally, this is equivalent to testing F st between sampling times. Indeed, the F st , commonly used to describe differentiation among populations, is defined as a fixation index and can be seen as a measure of temporal divergence of populations by drift and not only as a (static) distance between equilibrium populations. The AMOVA implemented in Arlequin has been applied to several aDNA data sets [START_REF] Bramanti | Genetic discontinuity between local hunter-gatherers and central Europe's first farmers[END_REF][START_REF] Dalen | Ancient DNA reveals lack of postglacial habitat tracking in the arctic fox[END_REF][START_REF] Malmström | Ancient DNA reveals lack of continuity between neolithic hunter-gatherers and contemporary Scandinavians[END_REF], although these analyses did not allow for the distinction of population effects from temporal ones. In a different approach, [START_REF] Valdiosera | Surprising migration and population size dynamics in ancient Iberian brown bears (Ursus arctos)[END_REF] applied the analyses separately on different geological layers to assess geographical structure through time, though the confounding effect of sampling time may not be entirely overcome since some layers may be more heterochronous than others (e.g. Pleistocene vs modern). It should be noted that the AMOVA framework was originally designed for contemporaneous (modern) data and presents several limitations for the analysis of heterochronous data within the ARLEQUIN implementation. Ideally factors should be nested, and including both sampling time and geographic location in an analysis would be more appropriate if a set of sampling times (or time bins) is found only in a specific location (or conversely if samples from one population all belong to one chronological layer). An additional constraint is that the available implementation of AMOVA does not quantitatively take into account sampling time differences. [START_REF] Excoffier | Analysis of population subdivision[END_REF] notes that a population delimitation issue arises prior to the application of AMOVA and refers to associated attempts through aggregation techniques. This issue of population delimitation also applies to time bins. The qualitative feature of the analyses makes it strictly applicable only for pairwise analyses between two sampling time bins. It would be appropriate for some studies of modern vs ancient data if the latter do not span a substantial time range, but frequently this is not the case (e.g. [START_REF] Barnes | Dynamics of Pleistocene population extinctions in Beringian brown bears[END_REF][START_REF] Lambert | Rates of evolution in ancient DNA from Adelie penguins[END_REF][START_REF] Shapiro | Rise and fall of the beringian steppe bison[END_REF]. More quantitatively, it is possible to correlate the pairwise genetic distance with a pairwise time difference through the Mantel test or to both time difference and geographical distance using a partial Mantel test using IBD (Table 1). Finally, there is a generalization of AMOVA (GAMOVA) where sampling time can potentially be treated as a quantitative explanatory cofactor together with other factors [START_REF] Nievergelt | Generalized analysis of molecular variance[END_REF], although the procedure has not been applied for this purpose yet. Beyond these statistical tools to test for heterochrony and assess the proportion of genetic variation explained by sampling time differences, parametric (model based) statistics offer a powerful complement. [START_REF] Liu | Test of genetical isochronism for longitudinal samples of DNA sequences[END_REF] proposed two classes of standardized summary statistics aimed at detecting measurable evolution between two subsets at different times: the first one, D c , is related to a [START_REF] Nei | Molecular evolutionary genetics[END_REF] net distance D a (see Fig. 1). The other class of statistic, T c , quantifies the excess of polymorphism exclusive to one of the two subsets when compared to the isochronous case (i.e. the null hypothesis being tested). The tests proceed through either (i) parametric simulation if it is desirable to take into account the stochasticity of the evolutionary process, or (ii) permutation between the subsets, which may partly lessen confounding historical effects and would be appropriate to test if samples could be pooled for subsequent analyses. Both approaches should in general be complementary, the former being generally more robust. Such statistics are largely inspired from coalescent theory [START_REF] Kingman | The coalescent[END_REF] which is discussed in the next section, in particular its extension to heterochronous data.

Making inferences under the standard neutral model

The serial coalescent Coalescent theory models the genetic history of a sample by probabilistic genealogies where the nodes represent most recent common ancestors and where the lengths of the branches -or lineages -are proportional to time (see [START_REF] Wakeley | Coalescent Theory: An Introduction[END_REF] for a review). The serial (heterochronous) coalescent is a simple extension of the isochronous standard case [START_REF] Rodrigo | Coalescent Approaches to HIV Population Genetics[END_REF]. Heterochronous data can be described as a list of subsets, each defined by a sampling time and the associated subset of sequences. Events occurring on such a genealogy (e.g. coalescence, mutation, recombination, migration) are limited to lineages that co-occur at the same time point. Thus, at a given time, the number of extant lineages governs the rate of those events. Analytical derivations on the serial coalescent may help to predict the heterochrony effect or to interpret data, to suggest new statistics aimed at investigating heterochrony (see previous section) and to provide necessary corrections for most available statistics so as to enable comparison of data sets with different time sampling schemes. [START_REF] Fu | Estimating mutation rate and generation time from longitudinal samples of DNA sequences[END_REF] has proposed a diversity-based estimator of mutation rates based on serial samples; alternatively this method can be used to estimate sample age when an independent estimate of the mutation rate is available. This approach was subsequently refined by [START_REF] Liu | Summary statistics of neutral mutations in longitudinal DNA samples[END_REF]. For simplicity [START_REF] Liu | Summary statistics of neutral mutations in longitudinal DNA samples[END_REF] focused on a two sampling time case for both the standard neutral model of Wright-Fisher (WF;[START_REF] Wright | Evolution in Mendelian populations[END_REF]), and a model involving any deterministic population size change starting after the first sampling time. From this approach mean, variance and covariance is derived for the number of mutations affecting i individuals in a temporally older sampled subset and j individuals in a more recent sampling. Although the derivation relates to mutation counts, or segregating sites in the infinitely many site model (IMSM; [START_REF] Watterson | On the number of segregation sites[END_REF], more general cases with other mutational models can be estimated with corrected distances. From these counts commonly used summary statistics are derived (see Fig. 1 for examples) such as the total number of mutations in the sample, and hence the related [START_REF] Watterson | On the number of segregation sites[END_REF] θ W estimator of polymorphism as well as [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF] π diversity estimator. Testing the standard model Neutrality test statistics (see Nielsen 2005 for a review) can potentially be corrected for heterochrony within a general time sampling scheme with an arbitrary number of sampling times. Assuming WF and IMSM, [START_REF] Forsberg | Tree measures and the number of segregating sites in time-structured population samples[END_REF] derived the full probability distribution for the time to the most recent common ancestor of the total data set, the total length of the genealogy and the length of lineages exclusively ancestral to ancient samples (see lineages marked with dotted lines in Fig. 1). Such time quantities (time to MRCA, length of genealogy, length of lineages exclusively ancestral to ancient samples) can be related to the number of segregating sites, the number of mutations affecting only ancient samples (note that this may include some possible fixed differences between ancient and modern samples) and the mutations shared between ancient and modern samples. Thus, it is possible to perform the numerical computation of the probabilities to obtain a given summary statistic value in order to perform a neutrality test. Although the analytical derivation of means and variances of summary statistics is typically feasible for simple models, deriving their full distribution (as in the work described in the previous paragraph) is generally more challenging. Simulation approaches provide a sensible alternative to test evolutionary models or for parameter inference. Coalescent simulations are commonly used to empirically investigate the distribution of polymorphisms under flexible evolutionary scenarios, to test those scenarios through comparison with data, or make inferences for associated parameters in a more refined approach such as MCMC and ABC (see latter section). Simulations also allow for flexibility of model specification and provide efficient simulation schemes because one need only consider explicitly the history of a sample, and only at times where events modify that history, such as common ancestry, mutation, recombination or migration. Simulation within the serial coalescent algorithm is a straightforward minor modification of the contemporaneous case [START_REF] Achaz | A robust measure of HIV-1 population turnover within chronically infected individuals[END_REF][START_REF] Anderson | Serial SimCoal: a population genetics model for data from multiple populations and points in time[END_REF][START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF]) and provides a powerful tool to assess heterochrony effects, test evolutionary scenarios, design appropriate sampling schemes, and investigate the statistical properties of available methods. In the context of assessing the effects of heterochrony on neutrality test statistics, coalescent simulation has recently been used to reveal that heterochrony can introduce substantial bias to parameter estimation from summary statistics [START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF]. Heterochrony increases coalescence times by increasing the difference between sampling times, since the ancestral lineage of the youngest individual must first reach the sampling point of the ancient one before the two lineages can coalesce. This leads to lengthening of a genealogy, and the overestimation of polymorphism with classical estimators such as diversity π (Fig. 1), indicating that direct comparison of polymorphism between data sets with sampled from different time points is inappropriate. Using a large aDNA data set for cave bears, heterochrony has been shown to strongly influence the conclusions of several neutrality tests [START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF]. A straightforward correction factor can however be easily implemented for the diversity estimator and between population distances [START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF]. For moderate heterochrony (sampling time differences below N e generations), a likely scenario for most available aDNA data, terminal branches (i.e. directly leading to sampled tips) tend to be proportionally more affected. This leads to star like trees, an excess of mutations with low frequency in the sample, as revealed by negative Tajima's D, and a deficit of associations between mutations (linkage disequilibrium statistics, being sensitive to mutation frequency and the shape of the tree, provide information even in the absence of recombination). The result of this is the spurious mimicry of typical signatures of demographic processes such as population expansion. Both [START_REF] Achaz | A robust measure of HIV-1 population turnover within chronically infected individuals[END_REF] and [START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF] have shown through simulation that heterochrony can lead to substantial genetic differentiation between sampling points or, more generally, between data sets with different time sampling schemes. A useful consequence of this is that such differentiation may be used to estimate the rate of drift and thus, the effective population size [START_REF] Achaz | A robust measure of HIV-1 population turnover within chronically infected individuals[END_REF].

Making inferences under complex models Likelihood and Bayesian inference using MCMC genealogy sampling

The likelihood of a given parameter value is a function proportional to the probability of the data given the parameter value, L(θ|D)=αP(D|θ) (we will ignore the constant of proportionality α from this point). This function is used in statistical modelling to make inferences about the parameter values of a model. One approach to make inferences is to quantify the parameter value that maximizes the likelihood function (maximum likelihood estimate) and use the likelihood function profile around that value to determine the confidence intervals for the estimate. A Bayesian alternative approach is to combine the likelihood with some prior probability (information or belief prior to the experiment) to obtain the posterior probability distribution, and use this distribution to obtain a point estimate and credible intervals of the model parameter (see [START_REF] Beaumont | The Bayesian revolution in genetics[END_REF] for further details). It is possible to calculate the probability of a genealogy given a demographic model, P(G|θ), by using the coalescent as a model [START_REF] Felsenstein | Estimating effective population size from samples of sequences: inefficiency of pairwise and segregating sites as compared to phylogenetic estimates[END_REF]. Additionally, calculating the probability of the data (genetic state of a sample) given a genealogy, P(D|G) has been established from the field of phylogenetics [START_REF] Felsenstein | Evolutionary trees from DNA sequences: a maximum likelihood approach[END_REF]. The genealogy of the sample is unknown; thus, in order to derive the likelihood, it is necessary to integrate over all possible genealogies, L(θ

| D) = P(D | G)P(G |θ) G ∑ .
This will be the likelihood of full data. In practice, only a selection of random genealogies will be used in this integration, because the number of possible genealogies is infinite (a large number of topologies and an infinite combination of branch lengths). However, most of the possible genealogies will be very unlikely and will contribute little to the estimation of the likelihood. Therefore, it is desirable to sample the genealogies in an efficient way, i.e. favouring plausible genealogies. Currently, the most popular framework to sample genealogies efficiently is the Markov chain Monte Carlo (MCMC) approach, in an implementation originally proposed by [START_REF] Kuhner | Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling[END_REF]. This method consists of a random walk over the space of the model parameter values and over possible genealogies, and it is designed in a way that it should visit each point of that space proportionally to its likelihood (see [START_REF] Wakeley | Coalescent Theory: An Introduction[END_REF] for further details and a review of various implementations). This approach has been implemented in several software packages targeting different demographic scenarios (see Kuhner 2009 for a review); among them only BEAST (Table 1) allows for the sampling of genealogies of genes collected at different time points [START_REF] Drummond | Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data[END_REF] or a model of DNA damage for aDNA [START_REF] Rambaut | Accommodating the effect of ancient DNA damage on inferences of demographic histories[END_REF]). The initial interest in developing BEAST was to take advantage of heterochronous data to estimate mutation rates [START_REF] Drummond | Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data[END_REF]. In the case of isochronous data, substitution rates are frequently estimated from a phylogeny using fossil or geological data to calibrate node dates. For heterochronous data, it is possible to estimate the number of mutations occurring in the time interval between samples (and, thus, the mutation rate) without using any external data to calibrate the genealogy (see [START_REF] Drummond | Measurably evolving populations[END_REF] for further details). This potentially allows one to obtain species-specific or population-specific molecular rate estimates that cannot usually be obtained from phylogenies because of the difficulty of assigning fossils to a particular lineage below the genus level. The MCMC method to estimate the molecular rate implemented in BEAST was first used in the analysis of Adélie penguin (Pygoscelis adeliae) modern and ancient DNA samples [START_REF] Lambert | Rates of evolution in ancient DNA from Adelie penguins[END_REF] and has subsequently been applied to a number of other aDNA data sets (e.g. [START_REF] Ho | Elevated substitution rates estimated from ancient DNA sequences[END_REF]), usually providing rate estimates higher than those obtained from phylogenies (this has lead to the debate mentioned in the introduction). Analysis of complex historical change in population size can be addressed in BEAST using the Bayesian skyline plot framework [START_REF] Drummond | Bayesian coalescent inference of past population dynamics from molecular sequences[END_REF]. This method is based on the classical skyline plot analysis [START_REF] Pybus | An integrated framework for the inference of viral population history from reconstructed genealogies[END_REF], in which population size is estimated from consecutive time intervals separated by coalescent events from a given genealogy (in practice, a maximum likelihood estimate of the genealogy). An estimate of population size is obtained for each time interval from its length, based on the expected coalescence time for the number of lineages present at in that time interval. The graphical representation of this piecewise demographic model yields a plot that evokes the skyline of a city (Fig. 2b). The Bayesian skyline plot is an extension of this approach with the following modification: (i) time intervals of the piecewise demographic model do not need to be defined by consecutive coalescent events (generalised skyline plot by [START_REF] Strimmer | Exploring the demographic history of DNA sequences using the generalized skyline plot[END_REF]); (ii) a prior is used to make population sizes of consecutive time intervals correlated; and (iii) the uncertainty of the genealogy is taken into account with the MCMC algorithm for sampling genealogies [START_REF] Drummond | Bayesian coalescent inference of past population dynamics from molecular sequences[END_REF]. Although the model parameters are the population sizes at each interval, and the times defining the intervals, posterior probabilities for them are ignored in the output of BEAST. Instead, a plot of the estimated population size and the associated credibility interval as a function of time is presented (Fig. 2c), which is more informative about the past demographic history than the model parameter values. The analysis of heterochronous steppe bison mtDNA data remains an enlightening example of this analysis and its development [START_REF] Drummond | Bayesian coalescent inference of past population dynamics from molecular sequences[END_REF][START_REF] Rambaut | Accommodating the effect of ancient DNA damage on inferences of demographic histories[END_REF][START_REF] Shapiro | Rise and fall of the beringian steppe bison[END_REF]). The coalescent model implemented within BEAST was (until recently) that of a single population, and, thus, did not allow one to address the study of spatially structured populations. Potential problems of imposing spatially structured heterochronous data onto such a model for the estimation of mutation rate have been pointed out [START_REF] Navascués | Elevated substitution rate estimates from ancient DNA: model violation and bias of Bayesian methods[END_REF], and although it has not been evaluated, similar issues may exist for demographic inference. However, two recent developments have opened the door to phylogeographical analyses with BEAST. First, [START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF] have introduced a Bayesian modelling of character evolution for the inference of ancestral states. Using geographical locations as character states, [START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF] advocate the inference of posterior probabilities for the location (state) of ancestral nodes and migration events (changes of state). This is a significant improvement over classical phylogeographical studies (e.g. [START_REF] Hofreiter | Evidence for reproductive isolation between cave bear populations[END_REF][START_REF] Leonard | Population genetics of Ice Age brown bears[END_REF] for two reasons: (i) the uncertainty on the genealogy reconstruction is taken into account (classical phylogeography is frequently based on a single 'phylogeny') and (ii) there is a rigorous statistical interpretation of the results (classical phylogeography is frequenctly reduced to a visual description of the geographical distribution of lineages). Nevertheless, it must be noted that the use of this model of character evolution does not change the model of the coalescent, which remains constrained to that of a single panmictic population and it is not clear whether this would be robust to some population structure scenarios. The second development is the implementation of the multispecies coalescent within BEAST by [START_REF] Heled | Bayesian inference of species trees from multilocus data[END_REF]. In this new method (named *BEAST) a set of species is considered to diverge (by successive bifurcations, following a birth-death model) from a common ancestral species. Divergence times and effective population sizes for each contemporary and ancestral species (i.e. the 'species tree') determine the coalescent probabilities used in the sampling of genealogies. The MCMC scheme works by integrating over gene genealogies and over species tree topologies to obtain estimates of the posterior probability distribution of the species tree (gene genealogies have only the role of a 'nuisance parameter' and are ignored in the output). The purpose of the method is to reconstruct the phylogeny from multiple loci. In particular, it addresses the problem of incongruent gene genealogy topologies among loci and with the species tree topology due to incomplete lineage sorting by using a coalescent model. Although the method is focused on phylogeny reconstruction, it can be applied to a set of populations to determine the 'population tree' and use it to make phylogeographical inferences, i.e. inferring and dating the vicariance or colonization events responsible of the divergence of those populations. Estimating the likelihood distribution from full data with the MCMC algorithm is in principle the best way to extract most of the information contained within the data. BEAST is currently the only implementation of this method for heterochronous data but unfortunately it only includes single population demographic models, divergence models without migration (*BEAST) and coalescent models without intragenic recombination. Regarding the lack of intragenic recombination in the coalescent model, this has not been considered a problem in the past, as most aDNA studies have targeted mitochondrial DNA (see [START_REF] Ramakrishnan | Using phylochronology to reveal cryptic population histories: review and synthesis of 29 ancient DNA studies[END_REF] for a review). However, new sequencing technologies are a promising tool for the characterization of nuclear genetic diversity in aDNA [START_REF] Millar | New developments in ancient genomics[END_REF]) and future studies would likely benefit from the implementation of recombination. These features may be available in future versions or new programs of MCMC coalescent sampling, but efficiently incorporating recombination into MCMC presents many challenges -in particular likelihood surfaces become very rugged and difficult to explore. As an alternative researchers may use alternative methods for the estimation of the likelihood based on summary statistics rather than on full data. Approximating the likelihood using summary statistics Estimating the likelihood of full data, L(θ|D), is computationally intensive, the methods are difficult to implement in a computer program, and these difficulties increase with model complexity and the size of a given data set. Alternatively, it is possible to estimate the likelihood of a subset of the information within the data contained in summary statistics, S: L(θ|S). This approach takes advantage of the ease of simulating pseudo-samples using the coalescent model and computing summary statistics for those pseudo-samples. The general idea consist in simulating pseudosamples under a range of parameter values and rejecting those simulations yielding summary statistics very different (using a predetermined threshold) to the summary statistics of the target sample. The proportion of accepted simulations is used as an approximation of the likelihood. This approximation of likelihoods is most frequently used in the Bayesian framework in what has been termed Approximate Bayesian Computation (ABC). Specifically, the model parameter values are taken from some prior distributions for performing the simulations. Simulations are rejected or accepted in the same way as described above. The distribution of parameter values from the accepted simulations is used as an approximation of the posterior distribution (e.g. see [START_REF] Pritchard | Population growth of human Y chromosomes: a study of Y chromosome microsatellites[END_REF]. The estimates obtained from this method are highly dependent on the arbitrarily chosen threshold value for the rejection step: values that are too large will increase the error by accepting pseudo-samples far from real sample; values that are too small will require a prohibitively large number of simulations in order to accept enough simulations for the approximation of the posterior distribution. In order to address this problem, [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF] proposed an additional regression step after the rejection. In this step a linear regression of the parameter values as a function of the summary statistics is calculated from the accepted simulations. The parameter values are then adjusted using the regression and the posterior distribution is estimated from the adjusted values. This adjustment is termed the regression algorithm and has been responsible for the recent success of ABC in population genetics because it allows accurate inferences with an affordable number of simulations (Fig. 3). ABC methods have been little used in aDNA phylogeographic and population genetics. One of the reasons for this could be the lack of user-friendly software to perform all the necessary steps of the analysis (but see DIY ABC; Table 1). However, there are at least three coalescent simulators (SERIAL SIMCOAL, COMPASS and NETRECODON; Table 1) that perform coalescent simulations of heterochronous data. With a minimal knowledge of programming, they can be used for the simulation step of the ABC and their output can be analysed in R (with the functions available from Mark Beaumont or SERIAL SIMCOAL websites, see Table 1) for the rejection and regression steps. This was the procedure followed by [START_REF] Chan | Bayesian estimation of the timing and severity of a population bottleneck from ancient DNA[END_REF] to characterize a bottleneck for the rodent Ctenomys sociabilis from aDNA and by Ghirotto et al. (in press) to study the genetic continuity of bronze-age and modern human populations of Sardinia; both works using SERIAL SIMCOAL and R to perform the analyses. A simple example for using COMPASS and R for ABC analysis can be found as supplementary material (it will reproduce the analyses presented in Fig. 3). With the demographic and mutational model flexibility of the currently available coalescent simulators, there is much potential in an ABC approach for the analysis of heterochronous DNA sequence data. This is particularly relevant for addressing models outside the scope of BEAST, such as structured populations with migration or genes with intragenic recombination, for which no other alternative is currently available. However, a note of caution is necessary. With so much flexibility one might be tempted to fit a complex model to data with low information content. As in any model-based inference, special care should to be taken when choosing the proposed model and the posterior and prior distributions should be compared. It is also advisable to evaluate only a few different models and use a procedure of model selection (such as that proposed with an R function by [START_REF] Beaumont | Joint determination of topology, divergence time, and immigration in population trees[END_REF]) to estimate the posterior probability for each model.

Model testing and model selection

Estimating the parameters of a single demographic model will not suffice in most study cases for which several demographic models are plausible and should be considered. Thus, it is necessary to be able to test the models or measure the goodness of fit of the models to the data. Coalescent simulation may be used to test particular evolutionary scenarios within a range of parameter values of interest for heterochronous data. Several aDNA studies have made used of SERIAL SIMCOAL to test different demographic hypothesis: genealogic continuity between ancient and modern populations (e.g. Belle et al 2006[START_REF] Bramanti | Genetic discontinuity between local hunter-gatherers and central Europe's first farmers[END_REF]), demographic changes (e.g. Ramakrishan et al. 2005;[START_REF] Valdiosera | Surprising migration and population size dynamics in ancient Iberian brown bears (Ursus arctos)[END_REF] or population structure (Fabre et al 2009(Fabre et al (doi: 10.1371/journal.pone.0005151) /journal.pone.0005151) ). The general procedure consist in choosing a set of models, fixing a single value for each parameter of the model (or few combinations of values) and simulating a large number of pseudosamples for each model. The testing of each model (for a fixed combination of parameter values) is performed by estimating p-values from the distributions of the summary statistic (see Guimaraes et al. 2009 (doi: 10.1093/molbev/msp126) for a test statistic combining the p-values estimated from different summary statistics). An important drawback of this approach is that it requires a large number of simulations to estimate the p-value for each combination of parameter values considered, which severely limits the range of parameter values that can be explored. In order to select the model with the best fit to the data, Belle et al (2009, doi: 10.1038/hdy.2008.103) propose to rank the models by the number of summary statistic for which the estimated p-value is lower than the threshold for significance. Although this ranking gives an idea of the plausibility of the different models it is a poor measure of goodness of fit because some summary statistics are expected to be correlated and the 'measure' is not quantitative to the departure from the model. Alternatively, Ramakrishnan and Hadley (2009) propose to calculate the Akaike information criterion using a crude approximation of the likelihood based on the estimated pvalues. Despite their past popularity, these methods suffer from serious limitations. It is our view that this kind of approach will be superseded by the ABC methods described in the previous section, as they allow to explore each model for a continuous range of parameter values (instead of a single or few fixed values) and offer more statistically rigorous procedures for model selection (see Ghirotto et al. in press) as an example for model selection procedures in ABC).

Concluding remarks

Ancient DNA is being incorporated in population genetic and phylogeographic analyses with increasing frequency. As such, careful consideration is required by researchers regarding available statistical approaches, and the appropriateness of these for the data to be analysed. Approaches are available to assess the influence of evolutionary change on sequences across sampling intervals, and to control for other variables, such as geographic location, that may contribute to genetic differentiation among samples of sequences. While these approaches currently present some limitations, it would seem that these could probably be overcome with further development of existing tools. The analyses available to researchers are perhaps less limited than they might first appear, but to increase flexibility researchers are required to explore options beyond stand alone analytical packages. For the analysis of single populations conforming to panmixia, a full-data likelihood approach as implemented in BEAST offers a versatile tool to make demographic inferences, and recent developments open the door for phylogeographic analysis with multiple population demographic inference [START_REF] Heled | Bayesian inference of species trees from multilocus data[END_REF][START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF]. For aDNA studies where sampled populations deviate from models available within BEAST (e.g. structured populations, gene flow, recombination), summary statistics and their approximation in a Bayesian framework provide for an alternative approach. However for all model-based analyses researchers are advised to take a cautious approach to the interpretation of their results by comparing posterior and prior parameter values, considering model violation, and the potential consequences of lack of information within heterochronously sampled data. We echo previous calls for the need for rigor in aDNA analysis [START_REF] Cooper | Ancient DNA: do it right or not at all[END_REF], but extend this call to the downstream analysis of this data for historical inference.

Figure 1. The coalescent and summary statistics. Three coalescent trees are presented. Mutations are indicated in bold on the resulting sequence alignments (right). There are S=4 segregating sites in all cases. Lineages exclusively ancestral to ancient DNA samples are marked with a dotted line. (A): isochronous case, average difference between sequences (total diversity) π=1.8. (B): moderate heterochrony, mutation 1 is shared between ancient and modern sequences, mutation 2 is exclusive to ancient sequences while mutation 3 and 4 are exclusive to modern sequences; π=2.0; Nei's net distance: D a =0. (C): large heterochrony, mutations 1, 2 and 3 lead to fixed differences between modern and ancient sequences while mutation 4 is exclusive to modern sequences; π=2.2, D a =3. However, estimates corrected for heterochrony [START_REF] Depaulis | Using classical population genetic tools with heterochronous data: time matters! PLoS ONE[END_REF]) are: π h =1.8 (A), π h =1.54, D ah =-1 (B) and π h =1.16, D ah =2.33 (C). Figure 2. Skyline plot. The classical skyline plot analysis uses the reconstructed genealogy of the sample (a) to obtain estimates of the population size for each time interval defined by consecutive coalescent events, resulting in a plot shape similar to the skyline of a city (b). The Bayesian skyline plot (c) implemented in BEAST (Table 1) takes into account the uncertainty within the genealogy of the sample; the final (smooth) plot represents the posterior probability density for the population size (typically the median and 95% highest posterior density interval) calculated from the skyline plots (in grey) for the genealogies sampled with the MCMC. Figure 3. Likelihood approximation through summary statistics. A fictitious data set of ten sequences (five mDNA and five aDNA sampled 0.2×N generations ago) containing 20 segregating sites is analysed. The objective of the analyses is to estimate the parameter θ=2Nµ of a constant population size model and an IMSM of mutation. Prior and posterior distributions (estimated with rejection and regression algorithms) from an approximate Bayesian computation analysis are presented (50% of simulations rejected from a total of 30 000 simulations; a particularly high rejection threshold has been chosen for a better illustration of the improvement by the regression algorithm). Simulations were performed with COMPASS (Table 1) under the infinite site model and ABC was performed with the R functions of Mark Beaumont (Table 1). See supplementary files for a script in R to reproduce the analyses represented in this figure.

Table 1 . Software discussed in the main text for the analysis of ancient DNA.

 1

	Program	Purpose	Method	Models	Web	Reference
	ARLEQUIN	Test heterochrony &	Nested AMOVA,
		population structure	Mantel test


V B Kazantsev

V I Nekorkin

S Binczak

S Jacquir

J M Bilbault

Spiking dynamics of interacting electronic neurons

Spiking sequences emerging from dynamical interaction in a pair of electronic neurons is investigated theoretically and experimentally. The model comprises two unidirectionally coupled FitzHugh-Nagumo units with modified excitability (MFHN). The The first (master) unit exhibits a periodic spike sequence with a certain frequency. The second (slave) unit is in its excitable mode and responds on the input signal with a complex (chaotic) spike trains. We analyze the dynamic mechanisms underlying different response behavior depending on interaction strength. The spiking phase maps describing the response dynamics are investigated. Complex phase locking and chaotic sequences are analyzed. We show how the response spike trains can be effectively controlled by the interaction parameter and discuss the possibility of neuronal information encoding.

I. INTRODUCTION

One of the key problems in the understanding of brain functions is how the neurons represent and bind sensory information converging to the brain from different sensory channels [START_REF]Principles of Neural Science[END_REF][START_REF] Winfree | The Geometry of Biological Time[END_REF][START_REF] Llinas | I of the Vortex: From neurons to self[END_REF]. Neurons exhibit electrical activity and communicate using spike (electrical pulse) sequences so that all the information is represented (encoded) in the characteristics of the sequences. The neurons are subserved by the synaptic interactions (chemical and electrical) providing the activity modifications from one cell to another. The resulting space-time patterns of electrical activity is believed to be a neuronal tool to represent and process the information [START_REF] Llinas | I of the Vortex: From neurons to self[END_REF]. A number of recent neurophysiological studies of the activity in different brain areas have shown that neuronal information "code" deals with the control of spike timing. Take, for instance, spiking synchrony and reset properties in the olivo-cerebellar dynamics [START_REF] Llinas | I of the Vortex: From neurons to self[END_REF][START_REF] Sugihara | [END_REF][5], δ-oscillation phase response in the thalamo-cortical system [6], spike-timing in cerebellar stellate cells [7], θ-oscillation phase precessing in hippocampus [8,9]. It appears that not only spiking rate but also spiking phase is important. Many recent works have been devoted to the study of "spike-communication-code" [10][11][12][13]. It has been found that integro-resonant properties of neurons as nonlinear dynamical systems underly the inter-neuron communications. Getting external stimuli the neuron may respond differently depending on both characteristics of the stimulus and internal dynamical state. ¿From neuron modeling point of view the basic concept of inter-neuron communication deals with the neuron response on incoming information messages. A number of recent works have contributed a lot on this issue [START_REF] Glass | From Clocks to Chaos: The Rhythms of Life[END_REF][15][16]. Among others the synchronous and chaotic response have been found in the different models. By varying the input frequency one can achieve various spike frequency locking modes and chaotic areas. It has been shown that the response can be characterized by the phase maps describing possible inter-spike timings.

In this paper we report on spike dynamics in the neuronal system comprising two units with a unidirectional interaction. The unit is described by the FitzHugh-Nagumo model [17] with modified excitability (MFHN) [13]. For some particular parameters its dynamics is similar to a Morris-Lecar system [18,[START_REF] Rinzel | Methods in neuronal modelling[END_REF]. In contrast with classical FitzHugh-Nagumo model typically used for cardiac cells modeling [15,16] the MFHN dynamics involves a number of separatrix loop bifurcations what will be essential for the response dynamics.

We assume that the input spike sequence is fixed with a certain frequency and come from the master MFHN unit in its oscillatory mode. The second (slave) unit is in the excitable mode and responds only when the external stimulation is strong enough. The unidirectional coupling can be seen as a simplest tool to model a synaptic interaction. It is activated when an input pulse comes and provides a dynamic delay of the response. We focus on the analysis of response sequences depending coupling parameter. We provide theoretical results on the dynamic (bifurcation) mechanisms of different response and experimental results obtained using MFHN electronic circuits [START_REF] Binczak | [END_REF].

The paper is organized as follows. In Sect. II we describe the model and its implementation with analog electronic circuit. In Sect. III we define 1D point map describing the dynamics of spiking sequences. Section IV is devoted to the map analysis and possible attractors of the map. We analyze basic bifurcation scenarios leading to complex (chaotic) spike dynamics. We also present the experimental study of the phase maps obtained with two MFHN electronic circuits in master-slave configuration. In Sect. V we develop the spike number analysis that could be a tool for "spike-communication-code". Section VI contains brief discussion of the results.

II. MODEL

The dynamics of two udidirectionally coupled MFHN electronic neurons is described by the following equations

u m = f (u m ) -v m , v m = (g(u m ) -v m -I m ), us = f (u s ) -v s + du m , vs = (g(u s ) -v s -I s ). (1)
The u m,s -variable describes the evolution of the neuron membrane potential, v m,s mimics the dynamics of ionic currents of the master and slave units, respectively. The function f has a cubic shape, f (u) = u -u 3 /3, the function g is taken piece-wise linear, g(u) = αu, if u < 0, and g(u) = βu, if u ≥ 0. The parameters α and β control the shape and the location of the v-nullcline, hence the excitation threshold. We fix the values α = 0.5, β = 2.

Parameter defines the relative time scales of u and v variable. Parameter I m,s controls the neuron depolarization level. The coefficient d describes the interaction, which is excitatory for d > 0 and inhibitory for d < 0.

O 2 O 3 O 1 u v O 2 O 3 O 1 u v (a) (b)
The dynamics of MFHN units (d = 0) is illustrated on the phase plane shown in Fig. 1.

The master unit (Fig. 1 2). In this circuit, the voltage U corresponds to variable u by setting u = R 0 U where R 0 is equal to R 2 R 4 /R 3 . The currents I 1 and I 2 are related to the variable v in Eq. ( 1) by v = γR 0 (I 1 + I 2 ) where γ is a fitting parameter so that the current voltage characteristic of part (A) is expressed by

I N L = 1 R 0 U -γ 2 U 3 R 5 R 6 (B) + - R R 3 4 R 2 R 1 U C D D D 1 2 3 4 D 5 L R 7 2 1 L 1 E 2 E D J (A) FIG. 2: MFHN electronic circuits. Electrical parameters : R 1 = 500 Ω, R 2 = 2 kΩ, R 3 = 8 kΩ, R 4 = 4 kΩ, R 5 = 1.2 kΩ, leading to R 0 = 1010 kΩ and γ = 1.138 V -1 . L 1 = 10.2 mH, L 2 = 3.5 mH, R 6 = 2021 Ω, R 7 = 690 Ω and E 2 = -0.

III. DESCRIPTION OF RESPONSE DYNAMICS A. Spiking phase sequences

Let us consider the case of the excitatory interaction when the excitable unit gets a periodic spike train with d > 0. Then, the dynamics of the slave unit is defined in the threedimensional phase space (u s , v s , t) and can be very complex. Since the unit has excitation threshold for sufficiently small values of d no response appear. In this case it displays lowamplitude subthreshold oscillations of periodic or chaotic shape. Let us consider d > d cr when spike response appears. To characterize the dynamics we introduce spiking phase map.

As the response spikes appear due to periodic stimulation we can treat the master signal as a reference (base) signal. Let a master spike appears at the time moment t m n , and the following response spike at the time moment t s n . Then, we define relative response spiking phase as

φ n = (t s n -t m n )/T, n ∈ Z ( 2
)
where T is the period of the reference signal (Fig. 3 small neighborgood of some piece-wise continuous curve Π(φ n ) that is called phase response curve (PRC).

B. Phase response curve

The dynamics of the slave unit (1) is described by two-dimensional nonautonomous system with periodic forcing in the three-dimensional phase space (u s , v s , t). The projection of its trajectories to the plane (u s , v s ) for the time series shown in FIG. 3 from neighborhood. The shape of the layer is similar to the shape of the limit cycle in the master unit (FIG. 4, dashed curve). Note also, that if response spike appear the points tend to curve interval, L, near the peak of the u s -variable. To describe the spike dynamics let us fix a half-line, M , in the phase plane (u, v) originating in some point O located inside the master limit cycle and hitting the limit cycle near the peak point of u-variable (FIG. 4 (b)). Let us set the initial conditions, t = 0, for the slave unit at the point K s corresponding to the intersection of L and M . Then, let us set the master unit at the point of the limit cycle having the angle ϕ with the half-line M accounted counterclockwise and defined in the interval [0, 2π]. This angle is proportional to the time interval, τ 0 , between the master and slave spikes. As time proceed, the trajectories of system (1) map the (u s , v s ) variables back to the half-line M in the neighborhood of the curve interval L. Let us denote this time interval by τ 1 (ϕ). Then, at t = τ 1 the master unit reaches the angle ψ(ϕ). Therefore, the system trajectories define a point map F : ϕ → ψ defined in the interval [0, 2π]. As the angles are proportional to the spiking times the dynamics can be also expressed in terms of master-slave time differences Π : τ 0 → (τ 1 + τ 0 -T ). By normalizing times on master oscillation period, φ n = τ 0 /T, φ n+1 = (τ 1 + τ 0 -T )/T we find the spiking phase map in the form Note that the approximate 1D map description of the response dynamics becomes possible due to the existence of an curve interval L approximating invariant curve on the (u s , v s ) section of non-autonomous 2D system. Note also that even if the return intersection point, K s , appears with a shift relative to L it further evolves with approximately the same phase.

φ n+1 = Π(φ n ), (3)
Such phase definition is similar to conception of isochrons recently introduced in [16].

IV. PHASE MAP ANALYSIS

.

A. Phase map attractors

Let us analyze basic spiking modes of the response and bifurcation scenarios of their appearance. It has been found that for large values of d the spike sequences from the master and slave units become synchronized with 1 : 1 spike frequency ratio. In this case slave spikes are phase-shifted relative to the master. The shift value vanishes with d → ∞.

This corresponds to the stable fixed point 0 < φ * < 1 of the map (3). (chaotic) structures alternating with periodic windows of complex spike frequency locking ratios. Let us consider how the PRC evolves with changing d. We illustrate the dynamics using normalized phase variable, φ n , mod1, defining the phase shift between the response spike and preceding master spike.

• For 1 : 1 locking mode near the boundary of chaotic behavior the map (3) has stable and unstable fixed points in the neighborhood of saddle-node bifurcation (FIG. 6 (a)). • Then, with decreasing d the saddle-node bifurcation takes place and the two fixed points disappear (FIG. 6 (b)). However, the map curve stays very close to the bisector line where trajectories can stay for a long time. This corresponds to a long lasting 1 : 1 spiking when the inter-spike phase shift is monotonically increasing. Then, trajectories jump onto another part of the map curve corresponding to 2 : 1 mode.

In one iteration this part is mapped back to saddle-node bifurcation region. It means that after missing a one spike the response dynamics returns back to 1 : 1 region with increasing phase shift. Such map behavior corresponds to intermittency mechanism of chaotic dynamics. The regular ("laminar") phase is intermitted by jumps ("turbulent phase") that breaks the "laminar" dynamics. Chaotic response signal corresponding to map attractor (FIG. 6 (b)) is illustrated in FIG. 7 (a).

• Along with chaotic trajectories the map can have periodic windows with complex spike frequency locking. Note that the present study does not cover all possible behaviors. In particular, we do not consider here various multi-stable spiking modes when two or more attractors co-exist in the map (3). Spiking modes may also co-exist with stable subthreshold responses. In this case the map (3) will have undefined intervals where no spiking response appears.

B. Spiking phase map in electronic experiments

The electronic circuit coupling the two electronic neurons is presented in Fig. 10.

V. SPIKE NUMBER ANALYSIS AND SPIKE ENCODING

We have studied the phase map attractors using normalized phase shift variable, of "1 and so on. 1 : 1 → {..., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...}

f n+1 f n mod1 f n+1 f n mod1 f n+1 f n (f) (e) (d) mod1 f n+1 f n (c) (b) mod1 f n+1 f n (a)
2 : 1 → {..., 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...}
3 : 2 → {..., 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, ...} 6 : 5 → {..., 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, ...} . . . . . . . . . In other words the sequences can be treated as a kind of "spike-communication-code" that 13 close trajectories diverge as time proceed, and the system "forgets" its initial state. However, for shorter time lengths (much less than the instability time scale) the digital route is welldefined and predictable. Thus in the chaotic case one can reset the initial condition, z 0 , and deal with finite interval routes taken from the chaotic map.

VI. DISCUSSION

We have investigated the dynamics of simple neuronal architecture comprising two interacting neurons. One of the neurons is oscillating sending to the other one a periodic spike train. The signal is transmitted through unidirectional coupling. Such an interaction can be a simplest view on a synaptic transmission when the incoming pulse causes the potential changes (postsynaptic potential) at the receiver cell. If this change exceeds the excitation threshold the cell exhibits a response pulse or spike. We have shown that the response spike sequences can be extremely diverse. Note that in such a master-slave architecture the timing lations. Moreover, when stimulated periodically the stable and unstable manifolds of the saddle point (FIG. 1 (b)), that define the excitation threshold, may form a complex homoclinic structures in the three-dimensional non-autonomous phase space. The presence of such structures can be noted in FIG. 6 where the response points form two inter-related scrolls with one corresponding to the subthreshold signal and the other to the spiking mode.

We have shown that the spiking response can be approximately described by 1D phase map. The PRC function has been obtained. PRC analysis has shown major bifurcation scenarios of chaotic response. The map displays different kinds of intermittency, complex periodic orbits and chaotic attractors. We have also described how the spiking response can be associated with digital sequences. Such sequences can be associated with a kind of "spike-communication-code" that could be used for internal representation of sensorial information in the neuronal systems. The information can be transformed in a "spiking route" that is controlled by the interaction strength and can be found by appropriate reset of initial spiking phase.

Finally, we have presented analog electronic circuit reproducing the master-slave neuronal configuration. The results obtained have confirmed the theoretical predictions. Spiking phase maps of regular and chaotic structure have been obtained experimentally. Indeed, the phenomena of coupling-controlled spiking response has robust properties and can be useful in perspective design of neurodynamic-based information processing systems.

FIG. 1 :

 1 FIG. 1: (a) Phase portrait of the master unit dynamics. The unit has three fixed points and stable limit cycle corresponding to periodic spike train. (b) Excitable dynamics of the slave unit. The saddle separatrix defines the excitation threshold of the stable rest state O 1 .

  (a)) has three fixed points and stable limit cycle for I m > I * ( ) corresponding to the periodic spike sequence. The slave unit is in its excitable mode, I s < I * ( ), where I * ( ) is a bifurcation parameter[? ]. It has stable fixed point, O 1 , corresponding to the neuron rest potential, unstable fixed point, O 3 , and saddle point O 2 . The excitation threshold is defined by the saddle incoming separatrix. If the perturbation of the rest state is strong enough, i.e. lies below the separatrix curve, then the neuron responds with excitation pulse. The MFHN unit has been recently implemented with analog electronic circuit (FIG.

6 V

 6 , leading to α = 0.5, β = 1.96. For the oscillating (master) unit we set: I = 0.22 (E 1 = 0.3876 V ) , = 0.2 (C = 1 nF ). Parameters, for the excitable (slave) unit: I = 0.19 (E 1 = 0.332 V ), = 0.2 (C = 1 nF ).

FIG. 3 :

 3 FIG. 3: Dimensionless phase plane (V,W) showing stability with excitation threshold. Parameters : η = 0.19 and ε = 0.2 (E 1 = 0.332 V and C = 1 nF ). Insets: Excitation pulse (upper inset) corresponding to (a); perturbation reaching resting state (lower inset) corresponding to (b). Parameters of insets: abscissa 20 µs/div; ordinate 225 mV /div.

FIG. 4 :

 4 FIG. 4: Dimensionless phase plane (V,W) showing bistability. Parameters : η = 0.226 and ε = 0.2 (E 1 = 0.3981 V and C = 1 nF ). Inset: Spiking train of pulses (upper inset)corresponding to (A); perturbation reaching resting state (lower inset)corresponding to (B). Parameters of inset: abscissa 20 µs/div; ordinate 500 mV /div.

FIG. 5 :

 5 FIG. 5: (a) Fragment of superimposed master and slave time series shown by dashed and solid curves, respectively. The spiking phase, φ n , is defined as a time shift between master and slave oscillations peaks. Parameter values: = 0.441, I m = 0.218, I s = 0.21. (b) Spiking phase diagramshowing the phases located in the neighborhood of some piece-wise continuous curve.

  defined in the half-line [0, +∞) by the PRC function Π(φ n ). The map (3) describes two main characteristics of the spiking sequences. (i) The phase shift between the response spike and the closest preceding master spike. In this case the spiking phase, φ n , is taken on modulus 1. (ii) The number z n = [φ n ] of integrated master spike to have the response.

  For smaller d complex spikes modes appear. The bifurcation diagram illustrating the map attractors is shown in FIG. 5 (a). The two enlarged regions in FIGs. 5(b,c) show a fine structure of the diagram for particular intervals of d. It appears that the map attractors can have complex

FIG. 6 :

 6 FIG.6: Phase definition in the (u, v)-plane. The master limit cycle is shown by dashed curve. The dots show the (u, v)-projection of nonautonomous slave dynamics. They occupy a thin layer that we approximate by invariant curve interval L (see text for details).

FIG. 7 :

 7 FIG. 7: (a) Spiking phase bifurcation diagram for normalized phase, φ n , mod1, depending on interaction strength d. For fixed d the dots show map attractors. (b), (c) Two enlarged regions of the diagram. Parameter values: = 0.441, I m = 0.218, I s = 0.21.

FIGURE 6 (••

 6 c) illustrates the example of 6 : 5 mode when the slave unit responds with 5 spikes for each 6 master spikes. With further decrease of d the number of map discontinuities increases. This corresponds to growing number of spikes to be missed in the response sequence (FIG.7(b)). In FIG.6 (d) the map has one unstable fixed point, P , with the multiplier, µ(P ) < -1. This point belongs to the map compartment corresponding to 2 : 1 spike locking. When the multiplier approaches to the bifurcation value µ(P ) = -1 the map displays another type of intermittency (FIG.6 (d)). Trajectories stay for along time rotating around the fixed point with increasing amplitude, then leave the region jumping over the discontinuity to another map compartment with different locking ratio. With changing d the points P acquires stability, |µ(P )| < 1 through the inverse period doubling bifurcation (FIG. 6 (e)). The left part of the bifurcation diagram on FIG. 5 (a) corresponds to a rare spiking response (FIG. 7 (c)). In this case has a large number of discontinuities corresponding to transitions between M : 1 spiking modes with M = 1, 2, 3, . . .. Map trajectories jumps from one mode to another forming a complex (chaotic) routes. Note that the continuous map compartment in the center corresponds to 1 : 1 ratio. The trajectory shown by solid circles hits this region at almost each cycle. Then, one can expect spike doublets in the response sequence what is clearly seen in FIG. 7(c).

φ n , mod1 .

 mod1 Let us consider another characteristics of the response sequence describing spike number dynamics, z n (d) = [φ n ] ∈ Z. This variable takes discrete values showing the number of slave spikes missed in the response sequence. In other words, it describes the number of "integrated" master spikes to get a response at the particular moment of time. So that the spiking sequences can be associated with sequences of digits. For example, simple frequency locking mode, 1 : 1, is associated with trivial sequence of zeros, 2 : 1 with periodic sequence

FIG. 8 :

 8 FIG. 8: PRC and corresponding map attractors shown by dots for different values of d. Parameter values: = 0.441, I m = 0.218, I s = 0.21. (a) Stable fixed point near saddle-node bifurcation, d = 0.07183. (b) Intermittency attractor emerged from the saddle-node bifurcation, d = 0.0716. (c) Stable periodic window with 6 : 5 spike locking ratio, d = 0.07093. (d) Intermittency attractor near period doubling bifurcation, d = 0.06878. (e) Stable fixed point corresponding to 2 : 1 spike locking mode, d = 0.068. (f) Map attractor for large number of PRC discontinuities, d = 0.06543.

( 4 )

 4 If the map (3) has one discontinuity (FIGs. 6 (a)-(c)) its dynamics can be described by binary sequences (4). Chaotic attractors (FIG. 6 (b)) of the map are described by chaotic sequences. The attractors, hence the sequences, are controlled by the value of d (FIG. 5).

FIG. 9 :FIG. 10 :

 910 FIG. 9: Response time series for different values of d (a) d = 0.0716. (b) d = 0.06878. (c) d = 0.06543.

FIG. 11 :

 11 FIG. 11: Bifurcation diagram (FIG. 5) for digital variable z.

FIG. 12 :

 12 FIG. 12: (a) A chaotic route for digital variable z in (z n , z n+1 )-plane. (b) Digital sequence {z n }. Parameter values: = 0.441, I m = 0.218, I s = 0.21, d = 0.06543.

. Thus, the dynamics of the circuit is given by Eq.[START_REF]Principles of Neural Science[END_REF] with = R 0 R 6 C/L 1 , α = R 0 /R 6 , β = R 0 (L 1 + L 2 )/(R 6 L 2 ), I = γR 0 E 1 /R 6 and time variable t = t/R 0 C.The experimental phase portrait corresponding to the two theoretical cases (A) and (B) presented in Fig.1are shown in Fig.3 and

respectively.

Acknowledgments

This research was supported in part by Russian Foundation for Basic Research (grants 03-02-17135), by grant of President of Russian Federation (MK 4586.2004.2). V.B.K. acknowledges Russian Science Support Foundation for financial support.