Anita Quiles
email: anita.quiles@gmail.com

E Aubourg

B Berthier

E Delqué-Kolic

G Pierrat-Bonnefois

M W Dee

G Audreu-Lanoë

" C Bronk

C Moreau

A Quiles Ab

E Delque-Kolic

G Andreu- Lanoë

C Bronk

Bayesian modelling of an absolute chronology

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

From the beginning of modern Egyptology, the establishment of an absolute chronology for Ancient Egypt has been an ambition which has concentrated the minds of many scholars. Pioneers in the field constructed a relative chronology by re-studying Graeco-Roman sources and deciphering ancient Egyptian writing. Firstly, Herodotus's writings were re-translated providing the first pieces of information that led up to Manethon's epitomes, which proposed lists of kings and the lengths of their reigns. These documents led to the first organization of Egyptian time into 31 groups of kings, called dynasties. The second major advance came with the studies of Champollion and his followers, who succeeded in understanding how to read hieroglyphs, the writing of the ancient Egyptians. Thus, Egyptian sources became comprehensible and archaeology brought to light documents of great interest which contained lists of the kings who reigned in Egypt. The Palermo Stone, the Abydos reliefs and the Turin Canon proposed lists of the kings and their order of succession. All these documents, to which are added a host of archaeological discoveries, have permitted Egyptologists to build a relative chronology for ancient Egypt. Henceforth, the 3000 years of dynastic history were divided into six periods. The three "stable" periods, called the Old, Middle and New Kingdoms are separated by three periods when the Egyptian state was divided between the North and the South: the First, Second and Third Intermediate Periods. Following the Manethonian tradition, these periods are organized into dynasties: the Old Kingdom (OK) consists of the Dynasties (Dyns) 3-6, the Middle Kingdom (MK) of the Dyns 11-13 and the New Kingdom (NK) of the Dyns.18-20. Dyns.1-2 precede the Old Kingdom and thus form the so-called Early Dynastic Period (EDP); Dyns.7-10 make up the First Intermediate Period (FIP); Dyns. 13-17, the Second Intermediate Period (SIP); Dyns.21-25, the Third Intermediate Period (TIP). Lastly, the Third Intermediate Period is followed by the Late Period (Dyn.26) and the Graeco-Roman Period. This study will focus on the 18 th Dynasty of the NK, for which museum held sufficient material to permit the development of relevant Bayesian models.

In Ancient Egypt, the notion of time was totally different to ours and relied on deep religious principles. Civil time was organized by a calendar based on agricultural seasons, divided into different units of time: a year contained three seasons (Akhet, Peret and Shemu) of four months and each month lasted thirty days. Five feast-days, called epagomenol days, were added at the end of the year. There was no larger unit of time than the year. Historical time was just expressed in reference to the ruling king: as examples, the construction of the Akhmenu of Karnak was described as being built in year 24 of Thutmoses III, and the battle of Qadech took place in year 5 of the reign of Ramses II.

Such an organization of time means that we would need to know the exact succession of all the kings and their precise length of reign, to re-establish an absolute chronology for Ancient Egypt. Unfortunately, our state of historical and archaeological knowledge does not permit the reconstruction of such a long period of history and Egyptologists used to date an event by citing the relevant Kingdom or Intermediate Period and the Dynasty and possibly the king under which the event occurred.

In the present study, we propose a new approach based on a statistical method, called Bayesian modelling, to establish an absolute chronology for Egypt's 18 th Dynasty. In this model, which includes different absolute astrophysical and physical methods, radiocarbon dates give the base-information, which is then constrained by results deduced from Sothic and Lunar methods calculated with a Bayesian approach, as well as Egyptian textual sources, and incorporated as a priori to the radiocarbon information.

Materials and Methods

2.1. Bayesian modelling A probabilistic approach, called Bayesian statistics, allows us to define observations taking into account our global knowledge on a studied object/event. Bayesian modelling is based on subjective hypotheses that represent the state of the established knowledge when the dating is performed. Such an approach highlights the combination of two systems of time: the first is relative time (represented by the archaeological evidence) and incorporated as a priori to a second time, which is absolute time, represented by measurement and called the likelihood. The age deduced by combining archaeological and dating information, is the a posteriori law, and this distribution has to be accepted or rejected depending on particular criteria that we have to define in advance. Such an approach seems to be particularly appropriate for dating archaeological objects/events since contextual as well as stratigraphic information is usually established before measurement. In particular, radiocarbon dates are rarely the only information we have on analyzed objects. The archaeological context in which they were found/conserved as well as the correlations between objects (inbuilt age, contemporaneity, position of one object to another, etc…) give extra information. This may improve the precision of the age distributions deduced by radiocarbon dating, by constraining the radiocarbon densities by using archaeological evidence. So, this Bayesian approach was firstly performed on Lunar dating, by writing our own C++ programme. Then, we modelled radiocarbon measurements using OxCal 4.1, and we used priors (combine, phase, sequence, gap, interval, boundary, Delta_R), agreement and convergence factors following the terminology proposed in (Bronk Ramsey, 2009).

Sothic Dating

Sothic Dating was the first method used by Egyptologists to establish the backbone of an absolute chronology for Ancient Egypt (Borchardt, 1899;Neugebauer, 1929;Parker, 1977;Krauss, 1985). This method is based on the observation of the star Sirius/Sothis, called Sopdet by the ancient Egyptians. Sopdet was an Egyptian goddess who had been known since the first Egyptian dynasty, and who was identified with the star Sothis, the brightest star in the night sky. Each year, Sothis disappears from the sky of an observer in Egypt and reappears about seventy days later. This phenomenon is called a Sothic rising. By pure coincidence, the time of reappearance of this star in the Egyptian sky corresponded with the time of the annual Nile flood. As a result, Sopdet was considered the goddess of fertility. This coincidence led Egyptians to base their calendar on the Sothic rising. When their calendar was created, the first day of the first month of the year was a day close to the Nile flood, and was also the day of the reappearance of Sothis. However, the Egyptian year lasted 365 days whereas the Sothic year was of 365 days plus one quarter. Thus, a gap of one day per four years had been created between this Egyptian civil year and the Sothic year. It meant that four years after the creation of the calendar, the Sothic rising occurred the second day of the year, eight years after, the third day etc.… So, a Sothic cycle is the time which separates two Sothic risings which have occurred the same day in the civil Egyptian calendar. It lasted 1460 Julian years (365*4=1460), which means 1461 Egyptian years. In reality, this value is not totally exact: the Sothic cycle is not fixed because of the own movement of the star and the length of Sothic year evolves in time.

Egyptians were aware of this offset between their civil calendar and the Sothic cycle but never modified it. Thus, they sometimes engraved the correspondence between the heliacal rising of Sothis's star and the date on which it occurred (king, year, month and day). These "Sothic equations" are of the utmost importance because an absolute chronological date for the year of the ruling king can be deduced from them. Sothic documents are attested, in particular a Sothic equation was found on a stone excavated at Elephantine Island, near Aswan. This block refers to an offerings' calendar, and in this document, presently exhibited at the Louvre Museum (E 3910), it is written:"3 th month of the season Shemu, day 28, day of the feast of the coming of Sopdet, …" which means that a Sothic rising occurred the 28 th day of the third month of the season Shemu (Shemu.III.28), under an unknown king. A Roman writer called Censorinus wrote in De Die Natali (Censorinus) that the 19 July 138 C.E was the first day of this Egyptian civil year. This association gave us the first piece of information to evaluate absolute dates on which Sothic risings occurred. Then, we used a method, based on the arcus visionis estimation. It is the smallest angle difference between the sun and the rising of a star that is necessary for the celestial object to be just visible at its rising by a person on the earth. Our method estimates the absolute date on which the Sothic rising occurred, according to Julian years (Aubourg, 2000).

2.3.Lunar Dating

In parallel with the civil calendar, a lunar calendar was used in ancient Egypt, especially in temples (Parker, 1950 ;Depuydt, 1997). It was based on 12 months of 29 or 30 days, and the first day of the month was the first day of the moon's invisibility. Thus, the length of the month was determined by observation of the moon on the 29 th day: if the moon was again visible, the month lasted one more day and the new one begun the day after the following day, but if it was invisible, this following day was considered the first day of the new month. Some lunar equations are attested in the texts, in particular during the reign of Thutmoses III: they consist in associating one phase of the moon with a date in the Egyptian civil calendar. Now, days and hours of new moons around the reign of Thutmoses III have been tabulated thank to the JPL-HORIZONS software (JPL Horizons) developed by the NASA. With a Bayesian approach, we have determined the more probable dates of observation of new moons attested in Egyptian texts, by synchronizing these lunar equations with the previous Sothic method.

2.4.Radiocarbon Dating

At LMC14 (CEA Saclay, France), two sets of radiocarbon analyzes were carried out on short-life objects conserved at the Louvre museum and archaeologically attributed to a specific reign or short period of the 18 th Egyptian Dynasty. Care was taken to sample objects without preservatives. Samples were treated using Saclay's routine pretreatment process for organic material (acid -base-acid procedure consisting of HCl (0.5N, 1 hour, 80°C), NaOH (0.1N, 1 hour, 80°C) and HCl (0.5N, 1 hour, 80°C)). They were then dried and combusted at 900°C with CuO and silver wire. CO2 produced was collected cryogenically before being reduced to graphite with H2 and Fe powder at 600°C. Measurements were performed with the AMS method at the ARTEMIS Facility. Calibration and Bayesian modelling were realized using the OxCal4.1 program (Bronk Ramsey, 1995), using IntCal09 curve (Reimer et al, 2009). Radiocarbon results are used to be given with one standard deviation and calibrated ages, with two sigmas deviation. In the modelling, an offset of 19 ± 5 years was added to each radiocarbon date, following the study of the Oxford Laboratory on the reservoir offset in Egypt (Dee et al, 2010) and using the Delta_R function (Bronk Ramsey, 2009).

2.4.1.Sennefer's tomb at Deir el-Medineh

Bouquets of flowers found in Sennefer's tomb at Deir el-Medineh by the French Egyptologist Bernard Bruyère (Bruyère, 1929) (Figure 1), were radiocarbon dated. 47 samples were analyzed; we sampled different short-life plants (leaves, twigs, etc..) to ensure the consistency of the dates. They came from seven different bouquets. The archaeological material found inside the tomb shows three burial phases. They occurred between the beginning of the reign of Tutankhamun and the beginning of the reign of Horemheb, which means a period of about 15 years (Bruyère, 1929). In his excavation reports, Bernard Bruyère states that he found all the bouquets at the entrance of the tomb, which means that such material precisely dates the same archaeological event: one of the three phases of burial.

2.4.2.Basketries from Deir el-Medineh

Basketries from the eastern cemetery of Deir el-Medineh (Figure 2) are held at the Louvre Museum. They were found by Bernard Bruyère in the excavations he led between 1928-39, then offered to the Louvre Museum as in recognition of the museum's support for the excavation. The cemetery was used by a relatively modest, almost exclusively female, population. It can be inferred by the study of the archaeological material found inside these tombs (ceramic, scarabs…) that burials in this cemetery occurred at the beginning of the 18 th Dynasty. In particular, scarabs with the names of queen Ahmes Nefertari, queen Hatshepsut, and king Thutmoses III are attested (Bruyère, 1937;Pierrat-Bonnefois, 2003). 19 basketries were analyzed by 78 measurements performed on short-life samples of palm, alfa, grapes and seeds, 53 by LMC14 laboratory and 25 by ORAU laboratory. Some of these results have already been incorporated in the 'Radiocarbon-based chronology for dynastic Egypt' proposed by the Oxford laboratory (Bronk Ramsey et al, 2010).

Results

Sothic Dating

The Elephantine Calendar refers a Sothic rising on Shemu III.28 at Elephantine Island, which corresponded to the 328 th Egyptian civil day. To elaborate, Figure .3 shows that around 1450 B.C.E, the Sothic rising occurred at Elephantine at about the 10.5 ± 0.75/0.5 Julian July. This 10.5 ± 0.75/0.5 Julian July has to be the 328 th or the 327 th day of the civil Egyptian calendar, since we have to consider an uncertainty of one day between the observation day and the effective day of the Sothic rising. From Censorinus text, we know that 19 Julian July 138 EC corresponded to Akhet.I.1, from what we calculated that the 11/10 Julian June 138 EC corresponded to Shemu.III.28/27. By subtraction, we deduced that Shemu.III.28/27 corresponded to the 10.5 ± 0.75/0.5 Julian July between 1439 and 1448 BCE, which enable to conclude that this heliacal rising of Sothis took place in 1443.5 ± 4.5 B.C.E (Figure 3). Unfortunately, the year of the ruling king on which this observation was made isn't mentioned, but it may firstly be shown that this block was engraved during the reign of Thutmoses III (Bommas, 2000), which lasted 53 years (Hornung, 2006). Besides, the study of the archaeological context on which the Elephantine calendar was found leads us to the conclusion that the wall on which it was located could not have been engraved before the year 33 of this king. Thus, we deduced that 1443.5 ± 4.5 B.C.E must be between the year 33 and 54 of Thutmoses III (Quiles, 2010).

In addition, for Egypt's 18 th Dynasty, successions of kings and lengths of their reigns are relatively well-known, as summarized in TAB. 3. Combining the derived equation with the known lengths of reigns, termini were calculated for each reign of this dynasty. Thus, we deduced that 1568.5 ± 4.5 B.C.E. and 1285.5 ± 4.5 B.C.E. were respectively termini post and ante quos for the beginning and the end of the 18 th Dynasty. These termini were obtained without preferring one Egyptologist's estimations to another.

Lunar Dating

Thutmoses III relates two lunar equations in his Annals. The first one is attested as the day of the beginning of the Megiddo Battle. He explains: "Year 23, I st month of Shemu, day 21, day of the battle of Megiddo, day of the feast of the new moon" (Urk.IV.657.1). The following year, the king inaugurated the Akhmenu, a festival hall he built in the Karnak complex. He said: "My majesty orders to prepare the ceremony of foundation to come the day of the feast of the new moon, for the act of foundation of this monument, in the year 24, 2 nd month of the season Peret, last day" (Urk.IV, 836.2). Knowing the accession date of Thutmoses III in the Egyptian calendar (Shemu.I.4, (Urk.IV.180, 15-16)) and being aware that during the New Kingdom, the first year of a king always lasted 365 days, we deduced that these two dates are separated by 649 days. Besides, thanks to Sothic dating we established years 23 and 24 of Thutmoses III were between 1480 and 1448 B.C.E. And by calculating the correspondent day of these two celebrations in the Julian calendar, we deduced that these new moons occurred between 1 st and 20 th May for year 23, and between 8 and 27 February for year 24 (Quiles, 2010). With JPl-Horizons software, the new moons of February and May between 1480 and 1448 B.C.E. were tabulated. Texts do not however specify whether the Egyptians observed the moon in the morning or in the evening. This information is of prime interest because the day of observation of a new moon could have been different from the actual day (for instance, if Egyptians used to observe the new moon in the morning, a new moon which was sighted at about 9 pm was observed on day D+1 rather than day D.). For this reason, we have written two models of calculation in C++ language, using a Bayesian approach, to determine the more probable day of observation. As likelihood, we have taken a Gaussian law of distribution, in view of the probability of being able to observe the new moon once you know it has occurred. Two a priori laws were tested, according to an observation in the morning, or in the evening. These a priori laws were defined as Heaviside functions except around the boundaries where the probability of observation was fixed with a Gaussian law (sunrise and sunset). The difficulty was to set boundaries of the hours up to which observation was possible. They were set in reference to the hours of sunrise and sunset (consequently, the boundaries are not the same for the new moon's observation in May or in February). As a result, the a posteriori law answered the question: "knowing the effective hour of the new moon, and postulating a type of observation (morning or evening), what is the probability the observation was done on day D or day D+1?"

After having tabulated the more probable days for observation of the two new moons (February and May) between 1480 and 1448 B.C.E., we checked the ones separated by 649 days. The sum of results of the two models enables us to propose a set of 12 years that are possible as year 1 for the reign of Thutmoses III. These are summarized in Figure 4. (Quiles, 2010) 3.3. Radiocarbon dating 3.3.1.Sennefer's tomb In the 47 analyses carried out on bouquets from Sennefer's tomb (TAB.1), five outliers were identified. At first we considered samples coming from the same bouquet as dating the same archaeological object, so having the same age. These ages were combined using R_Combine function in OxCal 4.1 (Bronk Ramsey, 1995) so as to get one age per bouquet, and they were then constrained by termini post and ante quos for, respectively, the accession date of Tutankamun and Horemheb, deduced from Sothic dating (1356.5 ± 4.5 B.C.E. and 1312.5 ± 4.5 B.C.E.). The seven ages obtained are contemporaries and spread out over a period from 1420 to 1260 calBC (Figure 5).

Because the seven bouquets were found one on top of the other at the entrance of the tomb, we can proceed to consider that all the bouquets precisely date to the same archaeological event. The combination of the seven probabilities, constrained by the same termini, leads to the proposition that the bouquets were put out in Sennefer's tomb between 1358 and 1312 calBC (Figure 5).

3.3.2.Basketries from Deir el-Medineh

The 19 basketries from the eastern cemetery of Deir el-Medineh were analyzed by one to nine measures (TAB.2) and 10 outliers were identified. The validated dates have been combined to get one probability density for the age of each object. 14 are in the reign of Thutmoses III whereas four basketries (E 14487, E 14479, E 14617 and E 16397) could be older (around the beginning of the New Kingdom). For specifying the period of activity of this cemetery, we have tried to identify the original tomb of some of the basketries. The study done by Y. Gourlay and the description of the tomb detailed by B. Bruyère permit the finding of the original tomb of nine objects, but it was not possible to proceed for the last ten (Gourlay, 1981 ;Bruyère, 1937). As an example, we observe that four basketries come from tomb 1382, in which three sarcophagi were found (Bruyère, 1937). These three burials did not occur at the same time, so we cannot say that the four basketries have the same age. Moreover, we cannot consider that the basketries precisely date the tombs or a burial phase, because a person could have been buried with a basket that was 10/20 years older. So an archaeological uncertainty of one to two generations has to be applied between the age of the basketries and the ages of burial phases in these tombs.

Egypt's 18 th Dynasty model

The comparison of dates deduced from Sennefer's tomb and the ones from basketries, synchronized with Sothic and Lunar methods, allow an absolute chronology for Egypt's 18 th Dynasty to be modelled with a Bayesian approach.

Accession date for Tutankhamun

The association of the ages of the Sennefer's bouquets with a historical period is more precise than for the basketries, as the former are associated with the interval from the beginning of the reign of Tutankhamun until the beginning of the reign of Horemheb, which is a period of 15 years, whereas an uncertainty of one to two generations surround basketries. That is why we chose to favor Sennefer's probability density for the construction of our chronology. However, to establish an absolute chronology we need a starting point associated with a precise historical event, like the accession date for Tutankhamun. So, we multiplied the combined density we got for "one phase of burial in Sennefer's tomb" by a rectangular function from 1 to 15 years, to simulate a distribution of age for the beginning of the reign of Tutankhamun. This age was then constrained by the termini post and ante quos for the beginning of this reign, deduced from the Sothic method. We obtained the result that Tutankhamun became king of Egypt between 1359 and 1319 calBC (2σ). This age takes into account archaeological as well as physical uncertainties (Figure .6).

Accession date for Thutmoses III

The information from Sennefer's tomb is more precise than the attribution of the basketries to Thutmoses III. The length between the beginning of the reign of the latter and Tutankhamun is relatively well known (between 140 to 149 years, TAB.3). That is why, in the same way, we simulated an age for the beginning of the reign of Thutmoses III, from the age of Tutankhamun (1499-1463 calBC (2σ), Figure .6). Then, the ages obtained with the basketries were incorporated in a "phase".

Additionally, the years calculated by Lunar dating as possible estimates for year 1 of the reign of Thutmoses III were added in another "phase". We defined these lunar dates by Normal distributed errors of one year. These two phases and the simulated date were integrated in another bigger phase, so as to constraint the simulated age by the two first phases. Finally, the density we get for the accession date of Thutmoses III was constrained by termini post and ante quos for the beginning of this reign, deduced by the Sothic method. The youngest lunar date was rejected by calculation (agreement factor was too low) and with a set of 11 lunar dates, the model is accepted with an agreement factor of 81. It suggests we consider that the reign of Thutmoses III began between 1499 and 1471 calBC (2σ).

The reign of Thuthmose III has been the central focus of several chronological analyses on the New Kingdom, using astronomical as well as historical approaches (Krauss 1985;Leitz, 1989;Grimal, 1988;Beckerath, 1994;Shaw, 2000;Hornung et al., 2006). The aim of this paper is to provide a new account, using only the most reliable archaeological information in combination with the latest chronometric techniques.

Accession Dates for each king of Egypt's 18 th Dynasty

Consequently, we have obtained two dates for two different reigns of the 18 th Dynasty and we have termini post and ante quos for the beginning and the end of this period, thanks to Sothic dating. That is, we have four corner points. For this period, the succession of the kings and the length of their reigns are well attested in history and archaeology. The majority of them are known with an uncertainty of about one to two years, but three are not as well-established and their uncertainty is of about 10 years. For the first group, we chose to define a probability distribution for the length of the reign by a rectangular function. This affords equal probability to each year of the function. Even if this kind of distribution law is restrictive for the calculation of the model, it represents our state of knowledge. On the other hand, the length of Thutmoses I st , Thutmoses II and Horemheb's reigns are not as well established. For instance, some Egyptologists affirm that Horemheb's reign lasted 13 years whereas others claim it reached 27 years. That is why we chose to define this probability density using the Before and After OxCal functions (Bronk Ramsey, 2009). Because the reign lasted at last 13 years, we used the restrictive After(13) function. Then, because it could have lasted more that 13, and until 27 years, we combined the former with a Before(27+T( 5)) function. Such a probability density makes the distribution more flexible. Lengths of Thutmoses I st and Thutmoses II's reigns were defined using the same functions.

As a result, we have simulated an interval of age per king by the multiplication of these relationships on the length of the reigns, to one of the two simulated radiocarbon's densities we got for Thutmoses III and Tutankhamun (Figure 7 and TAB. 3).

Implications on our model

Our model calculated an age distribution for the accession date of each king of Egypt's 18 th Dynasty, which synchronizes radiocarbon, Sothic, lunar, historical and archaeological information. Intervals are about 32 (Thutmoses III) to 45 (Horemheb) years with 95% probability (TAB. 3). They are shorter for the beginning of the dynasty which means the terminus post quem we calculated is relatively precise. For the end of the dynasty, the model is more flexible because intervals are larger. This can be explained by the fact the terminus ante quem we use is lower. It was calculated thanks to Sothic method, combined with historical information on the lengths of reigns. It means we have summed the maximum lengths for each reign after Thutmoses III, and taken as a starting point the first possible year possible for estimating the year of observation of the Sothic rising engraved on the Elephantine Calendar. Anyway, given our state of knowledge, we could not calculate a higher terminus ante quem without preferring one Egyptologist's interpretation over another. But, that could show that a couple of reigns were shorter that what we used, in particular for reigns for which the length is not well established. Another explanation would be the Sothic rising's observation could have occurred closer to the end of Thutmoses III's reign, rather than around year 33. But, this can not be definitively established.

This model may be further improved if we get information on another constraint (like a more precise knowledge on the length of a reign, another astrophysical observation, etc…). But, at present, modelled intervals are strongly correlated, which means new information on just one reign will significantly modify each interval. Getting more dates for different reigns that Thutmoses III and Tutankhamun seems to be the best solution to get an "asymptotic model". The more constraints we establish, the less intervals will vary.

5.

Conclusions Radiocarbon dating is a physical method which can be used to determine probability densities for the age of archaeological events. Thanks to historical and archaeological evidence, the succession of kings and lengths of their reigns for the 18 th Egyptian Dynasty are relatively well-established. Besides, we have astrophysical equations attested in Egyptians texts which allow us to calculate anchor points in Egyptian chronology. The first one is Sothic dating and in this paper we have shown, firstly, that these equations allow us to minimize uncertainties by using arcus visionis method. As a result, we achieved a relatively precise anchor point for the reign of Thutmoses III. Combined with our knowledge on the succession of the kings and the lengths of their reigns, we were able to fix termini for each of the reigns of the 18 th Dynasty. Lunar dates were then a second source of information and we developed a Bayesian model which calculates the more probable days of observations of new moons, given the type of observation (morning or evening). By incorporating the two lunar equations attested in the Annals of Thutmoses III, we deduced a set of 12 years as possible estimates for the accession date of this king. Then, two radiocarbon studies were carried out at LMC14 and ORAU on samples archaeologically attributed to a period Tutankhamun-Horemheb and Thutmoses III. The first leads to precisely dating one phase of burial in Sennefer's tomb, and from that result we simulate a date for the beginning of Tutankhamun, then for Thutmoses III. The second was constrained by densities we obtained on basketries from Deir el-Medineh, and by lunar dates. According to our knowledge on the length of each reign of the 18 th Dynasty, we finally determined a probability distribution for the length of each reign of this period, which allows the simulation of a period for the accession date of each king of this dynasty. These temporal probability densities incorporate Sothic, Lunar, radiocarbon, archaeological and historical information. Our model was applied on a reliable Egyptian period and results we got are in perfect agreement with dates previously deduced by some Egyptologists from methods unused here. This allows us to conclude that such an approach should certainly be applied to older Egyptian periods, for which the textual sources are less extensive.

Acknowledgements

This project was realized in the framework of a Physics PhD led at the LMC14 lab (AQ), supervised by B. Berthier (BB), and funded by the University Paris Diderot, and a Master of Egyptology at the University Paris IV-Sorbonne (AQ), supervised by E. Aubourg (EA) and D. Valbelle to whom we are grateful. Radiocarbon measurements were mainly performed at the Laboratoire de Mesure du Carbone 14 (CEA Saclay, France), funded by the CEA, CNRS, Ministère de la Culture et de la Communication, IRSN and IRD (AQ, EDQ, CM). We are indebted to all the LMC14 team for its rigorous work and its involvement in this project. Oxford dates were realized at the Oxford Radiocarbon Accelerator Unit (University of Oxford) in the framework of the Egypt Chronology Project, funded by the Leverhulme Trust (MD, CBR). Archaeological samples analysed in this study are held at the Département des Antiquités égyptiennes from the Louvre Museum which generously authorized us to collect them (GPB, GAL). We are also grateful to Dominique Bavay and Nicolas Durand for their thorough comments during the writing of this paper.

  E., 2000, Sirius et le cycle sothiaque, BIFAO 100, 37-46.  Beckerath, J.Van, 1994, Chronologie des ägyptischen Neuen Reiches. Gerstenberg Verlag, Hidlesheim.  Bommas, M., 2000, Der Tempel des Chnum der 18. Dyn. Auf Elephantine, Inauguraldissertation zur Erlangung der Doktorwürde vorgelegt der Fakultät für Orientalistik und Altertumswissenschaft der Universität Heidelberg.  Borchardt, L, 1899, Der zweite papyrusfund von Kahun und die zeitlische festlegung des mittleren reiches der ägyptischen geschichte, Zeitschrift für Ägyptische Sprache und Altertumskunde (ZÄS) 37, 99.  Bronk Ramsey, C., 1995, Radiocarbon Calibration and analysis of stratigraphy: the OxCal program Radiocarbon 37 (2), 425-430.  Bronk Ramsey, C., 2009, Bayesian analysis of radiocarbon dates, Radiocarbon 51 (1), 337-360.  Bronk Ramsey, C, Dee, M. W., Rowland, J. M., Higham, T. F. G., Harris, S. A., Brock, F. , Quiles, A., Wild, E. M., Marcus, E. S., Shortland, A. J., 2010, Radiocarbon-Based Chronology for Dynastic Egypt, Science 328, 1154-1557.  Bruyère, B., 1929, Rapport sur les fouilles de Deir el-Medineh, La tombe de Sennefer.  Bruyère, B., 1937, Les tombes du cimetière Est de Deir-el-Medineh, FIFAO XV.  Censorinus, De Die Natali, (French Translation): J. Mangeart, Livre de Censorinus sur le Jour Natal, Paris, 1843.  Dee, M.W., Brock, F., Harris, S.A, Bronk Ramsey, C., Shortland, A. J., Higham, T. F.G., Rowland, J. M., 2010, Investigating the likelihood of a reservoir offset in the radiocarbon record for ancient Egypt, JAS 37, 687-693.  Depuydt, L. 1997, Civil calendar and Lunar calendar in Ancient Egypt, Orientalia Lovaniensia Analecta 77.  Grimal, N., 1988, Histoire de l'Egypte ancienne, Fayard.  Gourlay, Y., 1981, Les sparteries de Deir el-Medineh, XVIII-XX dynasties. Documents de fouilles de l'IFAO, Le Caire.  Hornung, E., 2006, The New Kingdom, in Ancient Egyptian Chronology, Handbook of Oriental Studies, Brill, pp. 197-217.  Hornung, E., Krauss, R., Warburton, D. A., 2006, Ancient Egyptian Chronology, Handbook of Oriental Studies, Brill.  JPL Horizons, http://ssd.jpl.nasa.gov/?horizons.  Krauss, R., 1985, Sothis-und Monddaten, Studies zur astronomischen und technischen Chronologie altägyptens, Hildesheimer Ägyptologische Studien 20, Hildesheim: Gerstenberg.  Leitz, Ch., 1989, Studien zur ägyptischen Astronomie. Ägyptologische Abhandlungen 49, Wiesbaden: Harrasowitz.  Neugebauer O., 1929, Astronomische Chronologie, Berlin.  Parker, R. 1950, The calendars of Ancient Egypt. Studies in ancient oriental civilization n° 26, University of Chicago Press.  Parker, R. 1977, The Sothic date of the twelfth and eighteenth dynasties. In Studies in honor of G. R. Hughes.  Pierrat-Bonnefois, G., 2003, Cimetière est du village ou cimetière à l'est de Deir el-Medineh., In Deir el Medineh et la Vallée des rois, (ed) G. Andreu-Lanoë, Editions Khéops, pp. 49-65.  Quiles, A., 2010, Le règne de Thoutmosis III, mise en perspective des données participant à l'établissement d'une chronologie absolue et restitution du Calendrier d'Eléphantine dans son contexte d'origine. Master in Egyptology, University Paris IV-Sorbonne.  Reimer et al, 2009, Intcal09 and Marine09 radiocarbon age calibration curves, 0-50,000 years cal BP., Radiocarbon 51(4), 1111-1150.  Shaw I., 2000, Ed., The Oxford History of Ancient Egypt, Oxford University press, Oxford.  Urkunden des ägyptischen Altertums.

Figure 1 :

 1 Figure Captions:

Figure 2 :

 2 Figure 2: Basketry E 16391, held at the Louvre Museum. Two 14 C dates were performed on palm and halfa samples and the deduced combined density shows that this object was fabricated during the reign of Thutmoses III.

Figure 3 :

 3 Figure 3: Determination of the Sothic rising's day in Julian July, from (Aubourg, 2000, p. 46). We can observe on this graphic that around 1450 B.C.E, Sothic rising occurred at 24°N, at 10.5 ± 0.75/0.5 Julian July. Because the 19 Julian July 138 EC corresponded to Akhet.I.1, we deduced that 10.5 ± 0.75/0.5 corresponded to Shemu.II. 28/27 from 1448 to 1339 B.C.E.

Figure 4 :

 4 Figure 4: Results of Bayesian modelling of lunar dates attested in Thutmoses III Annals. A set of 12 estimates are given for the accession date of this king (in dark).

Figure 5 :

 5 Figure 5: Modelled age for each bouquet from Sennefer's tomb (upper). The seven bouquets combined age stands for the event "one phase of burial in the Sennefer's tomb", which occurred between 1358 and 1312 calBC (2σ). Termini Post and Ante Quos are deduced from Sothic method. A posteriori laws are in dark, likelihood in grey.

Figure 6 :

 6 Figure 6: Simulations of Tutankhamun (upper) and Thutmoses III (down) Accession Dates, from the Sennefer density obtained Figure 5. We deduced that Tutankhamun begun king in the period 1359-1319 calBC (2σ), and Thutmoses III between 1499 and 1463 calBC (2σ).

Figure 7 :

 7 Figure 7: Simulations of kings Accession dates for Egypt's 18 th Dynasty. Deduced intervals are about 32 (Thutmoses III) to 45 (Horemheb) with 95% probability.

  Results of 78 basketries samples from the eastern cemetery of Deir el-Medineh, undergoing AMS radiocarbon dating at LMC14 (SacA-53) and ORAU (OxA and OxA X (small carbon mass sample)-25) laboratories. ( * Département des Antiquités égyptiennes, Musée du Louvre). Outliers are identified in grey.

	Sample No. DAE* No. SacA 11161	Type Palm frond	mg C 0.50	δ 13 C -29	Radiocarbon age (BP) Radiocarbon age (BP) Radiocarbon age (BP) 3207 ± 31
	Sample No. DAE* No. SacA 15966 E 14000 Bouquet 1 SacA 18758 E 14000 Bouquet 1 SacA 15967 E 14000 Bouquet 1 SacA 18759 E 14000 Bouquet 1 SacA 15968 E 14000 Bouquet 1 SacA 18760 E 14000 Bouquet 1 SacA 15969 E 14000 Bouquet 2 SacA 18761 E 14000 Bouquet 2 SacA 11129 E 14477 E 14483 OxA 19485 E 14483 E 14483 E 14483 SacA 11130 E 14477 OxA 19152 SacA 11131 E 14477 OxA 19486 Cyperus Papyrus Type halfa grass halfa grass leaf leaf reed reed leaf leaf Palm frond Halfa grass Palm frond Palm frond Textile SacA 11132 E 14477 Date SacA 11162 Cyperus Papyrus OxA 19448 E 14477 Palm frond OxA 19449 E 14477 Palm frond E 14477 E 14477 E 16401 SacA 11163 E 16401 Palm frond SacA 11164 E 16401 Halfa grass E 16401 OxA 19450 Textile SacA 16407 Halfa grass OxA 19451 Textile SacA 11134 Palm frond SacA 11165 E 16393 Halfa grass SacA 11166 E 16393 Palm frond SacA 11167 Grape SacA 15970 E 14000 Bouquet 2 twig SacA 18762 E 14000 Bouquet 2 twig E 14489 SacA 11135 E 14489 Textile SacA 16397 E 14489 Halfa grass SacA 11168 Plant remains E 14489 SacA 11169 Halfa grass SacA 15971 E 14000 Bouquet 2 halfa grass or reed SacA 18763 E 14000 Bouquet 2 halfa grass or reed SacA 15972 E 14000 Bouquet 2 reed SacA 18764 E 14000 Bouquet 2 reed SacA 15973 twig SacA 18765 SacA 16398 Palm frond SacA 11137 E 16390 Palm frond E 14617 E 14617 E 14617 SacA 11170 E 14617 Textile E 14617 E 14617 E 14617 SacA 11138 E 16390 Halfa grass OxA X-2287 Plant remains SacA 11139 E 16394 Halfa grass OxA 19153 Halfa grass SacA 16399 E 16394 Halfa grass OxA 19452 OxA 19154 Textile E 16394 Halfa grass E Grape E twig SacA 11140 Halfa grass Palm	mg C 0.70 0.22 1.30 0.22 0.48 0.25 1.20 0.38 0.50 1.90 0.60 0.80 0.70 2.10 0.90 0.60 2.10 1.80 0.70 0.60 2.14 0.40 1.83 0.77 0.18 0.55 0.72 0.79 0.22 1.01 0.39 0.20 1.00 0.77 0.22 1.20 0.71 1.10 0.44 0.84 1.00 0.85 0.40 0.80 1.10 0.63 1.60 1.30 0.42 0.22 0.95 0.72	δ 13 C -28.4 -26.3 -29.5 -29.1 -25.8 -25.2 -30.3 -34.1 -28.6 -10.6 -26.8 -24.4 -22.7 -24.3 -16.8 -26.8 -27.1 -26.1 -25.1 -10.5 -23.7 -13.3 -23.9 -28 -15.4 -25.2 -27.8 -30.4 -32.1 -27.5 -10.1 9.5 -21.8 -26.6 -25.5 -30.0 -30.6 -28.6 -29,6 -29.3 -24.2 -11.3 -22.4 -10.8 -22.1 -10.7 -24.5 -10.7 -27.8 -25.8 -27.1 -13.2 -25.8	Radiocarbon age (BP) Radiocarbon age (BP) Radiocarbon age (BP) 3101 ± 34 3128 ± 26 3123 ± 39 3089 ± 26 3047 ± 36 3042 ± 29 3031 ± 36 2990 ± 3164 ± 3262 ± 3202 ± 3249 ± 3231 ± 182 ± 3171 ± 178 ± 3245 ± 3275 ± 3283 ± 3161 ± 3291 ± 3234 ± 3237 ± 3226 ± 3197 ± 3155 ± 3161 ± 44 3024 ± 37 3085 ± 3097 ± 3202 ± 3756 ± 112 ± 26 3079 ± 35 3170 ± 26 166 ± 30 -94 ± 43 3029 ± 35 3232 ± 3170 ± 3285 ± 3260 ± 3333 ± 3224 ± 110 ± 3188 ± 3209 ± 3227 ± 3241 ± 3192 ± 3060 ± 26 3250 ± 3234 ±	25 28 28 23 22 22 21 36 30 22 22 31 31 31 30 28 25 22 25 21 61 25 18 26 22 25 33 30 23 24 23 28 30 26 23 25 25
	SacA 15974 SacA 18766 SacA 11141 SacA 16403 SacA 16400 SacA 16404	E 14000 Bouquet 3 E 14000 Bouquet 3 E 14000 Bouquet 3 E 14000 Bouquet 3 E 14000 E 14000 E 16396 E 16396 E 16396 Bouquet 3 Bouquet 3 E 14665 E 14665 E 14665 E 14665 E 14665	leaf Palm frond Palm frond leaf Halfa grass Palm frond	1.40 1.00 0.53 0.23 0.52 0.50	-30.8 -25.6 -28.7 -32.4 -12.5 -28.3	3040 170 3200 2990 3191 3232	± ±	36 ± ± 26 ± ±	25 23 24 25
	SacA 15975 SacA 18767 SacA 11143 OxA 19147 SacA 16401 OxA 19148	E 16397 E 16397 E 16397	reed Halfa grass Halfa grass reed Halfa grass Palm frond	1.00 0.80 0.50 0.21 0.38 0.90	-25.4 -26.0 -10.6 -22.4 -13.4 -23.4	3088 165 3261 3113 3271 3186	± ±	33 ± ± 25 ± ±	21 32 29 28
	SacA 15976 SacA 18768 OxA 19146 SacA 16408 SacA 15977 SacA 18769 SacA 15978 SacA 11149 SacA 11144 SacA 11145 SacA 11148 TAB 2:	E 14000 Bouquet 4 E 14000 Bouquet 4 E 14487 E 14000 Bouquet 4 E 14000 Bouquet 4 E 14000 Bouquet 4 E 14000 Bouquet 4 E 14000 Bouquet 4 E 16391 E 16391	MO Halfa grass Halfa grass MO reed reed reed Palm frond Halfa grass Palm frond Halfa grass	0.91 1.60 0.45 0.86 0.83 0.79 0.31 0.80 0.08 0.46 0.25	-28.3 -10.7 -12.9 -25.7 -25.2 -24.6 -25.5 -19.4 -10.6 -26.2 -17.6	3015 127 3243 3137 3044 3131 2999 3280 3377 3196 3089	± ± ± ± ±	34 ± ± 26 37 28 33 ± ± ± ±	24 25 26 53 22 31
	SacA 15979 SacA 18770 SacA 16405 SacA 16406	E 14488 E 14488 E 14488 E 14488 E 14488 E 14488	leaf leaf Halfa grass Plant	0.90 0.60 0.86 0.42	-30.2 -30.5 -10.4 -15.8	3046 3120 3219 3163	± ±	36 27 ± ±	25 28
	SacA 15980 SacA 18771 OxA 19453 OxA 19480		OM Halfa grass OM Palm frond	1.10 1.70 0.86 1.90	-26.9 -12 -25.9 -24.4	3100 3264 3170 3251	± ±	41 ± 25 ±	29 26
	SacA 15981 SacA 18772 SacA 11150 SacA 11151		reed Halfa grass reed Halfa grass	0.75 0.30 1.20 0.90	-27.6 -8.2 -27,0 -17.3	3076 3185 2924 260	± ±	35 ± 27 ±	28 24
	SacA 15982 SacA 18773 OxA X-2287 SacA 16411	E 14000 Bouquet 5 E 14000 Bouquet 5 E 14000 Bouquet 5 E 14000 Bouquet 5 E 14000 Bouquet 5 E 14000 Bouquet 5 Bouquet 5 Bouquet 5 Bouquet 5 Bouquet 5 E 14000 E 14000 E 14000 E 14000 E 14491 E 14491 E 14491 E 14491 E 14491	leaf Halfa grass leaf Halfa grass	1.50 0.50 0.40 0.44	-30.3 -33.2 -8.7 -13.4	3135 3049 3283 3193	± ±	39 26 ± ±	33 27
	SacA 15983 SacA 18774 OxA 19149 SacA 11152		twig Halfa grass twig Halfa grass	0.60 0.90 0.23 0.60	-26.9 -22.1 -27.6 -10.2	1587 165 1178 3178	± ±	34 ± 29 ±	23 25
	SacA 15984 SacA 18775 SacA 11153 OxA 19481 SacA 15985 OxA 19150 SacA 18776 SacA 11154 SacA 15986 SacA 11155 SacA 15987 OxA 19482 SacA 18778 SacA 11156 SacA 15988 SacA 11158 SacA 18779 SacA 11159 SacA 15989 SacA 16409 SacA 18780 SacA 15990 OxA 19581 OxA 19151 SacA 16410	E 14480 E 14480 E 14480 E 14480 E 14000 Bouquet 6 E 14000 Bouquet 6 E 14000 Bouquet 6 E 14487 E 14487 E 14487 E 14000 Bouquet 6 E 14000 E 14000 E 14000 Bouquet 6 Bouquet 6 Bouquet 6 E 14000 Bouquet 7 E 14000 Bouquet 7 Bouquet 7 E 14000 E 14479 E 14479 E 14479 E 14479 E 14479 E 14479 E 14479 E 14479 E 14479	twig twig Palm frond Halfa grass OM Palm frond OM Palm frond OM Halfa grass wood ? Palm frond wood ? Halfa grass leaf Textile leaf Palm frond twig Halfa grass twig leaf Halfa grass Halfa grass Halfa grass	1.22 0.86 0.50 1.70 1.09 1.00 0.24 1.00 1.19 0.60 1.05 1.50 0.24 0.70 1.35 0.80 0.33 0.80 1.28 0.69 0.23 1.47 2.10 1.20 0.59	-27.9 -29.8 -28.1 -10.1 -30.4 -26.3 -25.9 -23.9 -28.3 -12.4 -28.4 -25.3 -27.3 -11.8 -29.2 -25.2 -36.9 -24.8 -27.9 -10.3 -25.3 -27.7 -9.9 -10.3 -10.5	3044 3116 3261 3233 3004 3153 2984 3304 3104 3041 3055 3277 3089 3198 3041 3284 2864 3305 3086 3259 3190 3036 3258 3107 3214	± ± ± ± ± ± ± ± ± ± ± ±	36 ± ± 29 ± 34 ± 28 ± 34 ± 33 ± 28 ± 34 ± 38 ± 35 ± 37 ± 28 ±	28 25 27 25 25 26 25 22 23 24 26 27 24
	OxA 19483		Textile	1.70	-24.6	3226		±	26
	OxA 19484		Palm frond	2,00	-25.2	3257		±	26
	SacA 11160		Halfa grass	0.70	-10.7	3200		±	22
		E 14483

TAB 1: Results of 47 samples from Sennefer Tomb undergoing AMS radiocarbon dating. Outliers are identified in grey ( * Département des Antiquités égyptiennes, Musée du Louvre).

Bayesian modelling of an absolute chronology for Egypt's 18 th Dynasty by astrophysical and radiocarbon methods


Zoran Tiganj
email: zoran.tiganj@gmail.com

Mamadou Mboup
email: mamadou.mboup@univ-reims.fr

Neural spike sorting using iterative ICA and deflation based approach

We propose a spike sorting method for multi-channel recordings. When applied in neural recordings, the performance of the Independent Component Analysis (ICA) algorithm is known to be limited, since the number of recording sites is much lower than the number of neurons. The proposed method uses an iterative application of ICA and a deflation technique in two nested loops. In each iteration of the external loop, the spiking activity of one neuron is singled out and then deflated from the recordings. The internal loop implements a sequence of ICA and sorting for removing the noise and all the spikes that are not fired by the targeted neuron. Then a final step is appended to the two nested loops in order to separate simultaneously fired spikes. We solve this problem by taking all possible pairs of the sorted neurons and apply ICA only on the segments of the signal during which at least one of the neurons in a given pair was active. We validate the performance of the proposed method on simulated recordings, but also on a specific type of real recordings: simultaneous extracellularintracellular. We quantify the sorting results on the extracellular recordings for the spikes that come from the neurons recorded intracellularly. The results suggest that the proposed solution significantly improves the performance of ICA in spike sorting.

Introduction

Information between neurons is transmitted using electrochemical signaling, which includes communication through action potentials (AP).

Multi-site electrodes with an increasing number of channel recordings are being more and more used for accessing neural activity (see, e.g., [START_REF] Simeral | Neural control of cursor trajectory and click by a human with tetraplegia 1000 days after implant of an intracortical microelectrode array[END_REF] where 100-channels recording system, described in [START_REF]Cerebus Data Acquisition System -users manual[END_REF], has been used). Each recording site will capture a mixture of activities from a number of neurons around.

Finding the firing instants of individual neurons form extracellular recordings is a very challenging problem in neuroscience, known as spike sorting (when recorded with extracellular electrodes action potentials are usually called spikes). For reviews and comparison of the existing spike sorting methods, we refer to [START_REF] Lewicki | A review of methods for spike sorting: the detection and classification of neural action potentials[END_REF], [START_REF] Gibson | Comparison of spike-sorting algorithms for future hardware implementation[END_REF] and [START_REF] Brown | Multiple neural spike train data analysis: state-of-theart and future challenges[END_REF] to name a few.

A preliminary step before actually sorting the spikes is to find and extract them from the original noisy measurement. This is the spike detection problem. For recordings with good Signal to Noise Ratio (SNR) spike detection can be solved by a simple thresholding [START_REF] Pouzat | Using noise signature to optimize spike-sorting and to assess neuronal classification quality[END_REF]. When the SNR is not good enough, different spike detection algorithms can be applied, e.g.: [START_REF] Kim | Neural Spike Sorting Under Nearly 0-dB Signal-to-Noise Ratio Using Nonlinear Energy Operator and Artificial Neural-Network Classifier[END_REF], [START_REF] Nenadic | Spike detection using the continious wavelet transform[END_REF] and [START_REF] Tiganj | Spike detection and sorting: Combining algebraic differentiations with ica[END_REF]. In general, only spikes with high enough amplitude can be detected. These correspond to the activity of the neurons closest to the recording sites. We mention however that, since the medium around the electrodes is not isotropic in general, the term closest does not always translate in the sense of minimum distance. So here and in all the remaining, closest neurons means neurons that are most visible in terms of current flow at the recording points during APs. The cumulative activity of all remote (in the sense above) neurons is considered as an undesired perturbation in the recorded signal, forming the background noise. And even for the neurons close to the electrodes, not all of their spikes are detected: some are inevitably lost due to the background noise corruption. However, analyzing the performance in spike sorting, we can say that missing some spikes (not detected spikes or spikes left unsorted) is less problematic than assigning spikes into a wrong cluster. We refer to [START_REF] Ilan | Information theoretic bounds on neural prosthesis effectiveness: The importance of spike sorting[END_REF] for detailed argumentation of this statement.

Commonly, spike sorting techniques use spike shape and amplitude as discriminative factors: even though all the neurons fire almost the same action potential waveform, due to the propagation and the velocity effects, spikes recorded with an extracellular electrode have different shapes and amplitudes if they are coming from different neurons [START_REF] Henze | Intracellular features predicted by extracellular recordings in the hippocampus in vivo[END_REF]. Thus, after spike detection, feature extraction techniques such as Principal Component Analysis (PCA) [START_REF] Adamos | Performance evaluation of PCA-based spike sorting algorithms[END_REF], wavelet based feature extraction [START_REF] Quiroga | Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering[END_REF], [START_REF] Hulata | A method for spike sorting and detection based on wavelet packets and Shannon's mutual information[END_REF] or some others [START_REF] Ghanbari | A graph-Laplacian-based feature extraction algorithm for neural spike sorting[END_REF], [START_REF] Chah | Automated spike sorting algorithm based on Laplacian eigenmaps and k-means clustering[END_REF] are commonly applied on the detected spikes.

Spikes are then sorted into different clusters, depending on the extracted features. A large number of clustering algorithms have been proposed for spike sorting. K-means [START_REF] Takahashi | A new approach to spike sorting for multi-neuronal activities recorded with a tetrode-how ICA can be practical[END_REF], EM [START_REF] Wood | Automatic spike sorting for neural decoding[END_REF], template matching [START_REF] Vargas-Irwin | Automated spike sorting using density grid contour clustering and subtractive waveform decomposition[END_REF], Bayesian clustering [START_REF] Prentice | Fast, scalable, Bayesian spike identification for multi-electrode arrays[END_REF] and super-paramagnetic clustering [START_REF] Quiroga | Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering[END_REF] are just some popular/recent examples. The number of clusters is often determined manually. However, sometimes some unsupervised procedures are used. As examples we mention Basysian Information Criterion (BIC) (e.g. [START_REF] Novak | Identifying number of neurons in extracellular recording[END_REF] and [START_REF] Fraley | How many clusters? Which clustering method? Answers via modelbased cluster analysis[END_REF]) as well as some simple ad-hoc procedures, such as finding the maximal density of the spikes in the feature vector space [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF]. The essential problem in spike sorting arises from the fact that the activities of neurons close to the electrode are often significantly destroyed by the background noise. Thus, very often different spikes fired by the same neuron appear in a recording as very different waveforms and it is very hard to recognize that they should be sorted into the same cluster.

In multi-channel (multi-site) recordings, the activity of each neuron in the neighborhood is captured by several recording sites. The principle of many spike sorting algorithms is to exploit such information redundancy. This is actually a Blind Source Separation (BSS) problem which is commonly approached via Independent Component Analysis (ICA). It allows one to find a new basis of data representation in which the mutual information between the induced virtual recording sites is minimized [START_REF] Shiraishi | Multi-neuron action potentials recorded with tetrode are not instantaneous mixtures of single neuronal action potentials[END_REF], [START_REF] Brown | Independent component analysis at the neural cocktail party[END_REF], [START_REF] Takahashi | A new approach to spike sorting for multi-neuronal activities recorded with a tetrode-how ICA can be practical[END_REF] and [START_REF] Takahashi | Automatic sorting for multi-neuronal activity recorded with tetrodes in the presence of overlapping spikes[END_REF]. For further reading, we refer to [START_REF] Comon | Handbook of Blind Source Separation, Independent Component Analysis and Applications[END_REF] for a nice tutorial on ICA (see also [START_REF] Hyvärinen | Independent component analysis: algorithms and applications[END_REF]). Unfortunately, the contribution of ICA is limited by several characteristics of neural recordings. The foremost is that the number of electrodes is much smaller than the number of neurons around. Also, the activities of different neurons are not independent. Consequently, ICA alone can not fully separate neural activities, but it can often transform the recorded signal so that the sorting is easier. Thus in general ICA is followed by some classical spike sorting algorithm (such as those we have mentioned in the previous paragraph).

In this paper we present a new algorithm for spike sorting from multi-channel extracellular recordings. The algorithm is iterative and consists of two nested loops. We use a deflation technique to improve the performance of ICA: after we separate the firing instants of a single neuron (using ICA and some spike sorting method), we remove them from the original recording and repeat the procedure until the algorithm becomes unable to separate any more neurons or until a prescribed number of neurons is obtained. Moreover, within each iteration we implement, in an internal loop, another iterative algorithm for removing the noise and all the spikes that are not coming from the neuron which is the closest to (the most visible by) the electrode. This neuron is identified as the one whose activity has the largest projection on the recording sites.

More precisely, the iteration of the internal loop consists of a sequence of ICA, and spike sorting. For the classical sorting algorithm we devised a simple greedy method, described in [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF]. The algorithm is based on the observation that the distribution of a neural signal deviates significantly from the uniform distribution and is rather unimodal. The detected spikes to be sorted are first processed with some feature extraction technique and then represented in a space with reduced dimension by keeping only a few most important features. The resulting space is next filtered in order to emphasis the differences between the centers and the borders of the clusters. Using some prior knowledge on the lowest level activity of a neuron, as, e.g., the minimal firing rate, we find the number of clusters and the center of each cluster. The spikes are then sorted using a simple greedy algorithm which grabs the nearest neighbors. This sorting method is selected for two reasons: 1) It does not require manual supervision, which would not be suited for an iterative algorithm. 2) It gives a good measure of distance between the clusters which, as we will see later, helps to determine which spikes can be deflated.

At the end of the external loop, we sort the overlapped spikes. This is done by taking all the possible pairs of the separated neurons and applying ICA only on the segments of the signal during which at least one of the neurons in a given pair was active.

The paper is organized as follows. In section 2 we describe two types of signal that we use in the simulations to come: 1) An artificially created neural mixture, where we simulate a four-channel extracellular recording in population of 1000 neurons distributed around the electrode. 2) A real simultaneous intracellular-extracellular recording [START_REF] Henze | Intracellular features predicted by extracellular recordings in the hippocampus in vivo[END_REF] we test the sorting algorithms on the extracellular recording and use the intracellular to verify the results for the intracellularly recorded neuron. In section 3 we give a step-bystep description of the proposed algorithm. In section 4 we give the complete results of the simulations for each step of the proposed algorithm. Concluding discussion is given in section 5.

Simulation settings

In this section we will describe the simulated and the real signal that we will use for the demonstration of the proposed algorithm in section 3 and for quantifying the results in section 4. With the simulated recordings we can validate the sorting accuracy for all the separated neurons, while with the real recordings we can validate the accuracy only for a single neuron -the one which was simultaneously recorded with the intracellular electrode.

Simulated neural recordings

To describe the proposed algorithm in an illustrative way we will use the simulated extracellular neural recording, in order to have the precise knowledge of which neuron fired which spikes. In this way we are able to clearly demonstrate the performance of the proposed algorithm. Also, using such signal, we are able to precisely measure the sorting accuracy and the sorting detection and thus compare the proposed approach with the other approaches.

To create the recording we first generate spatial positions for 1000 artificial neurons. We distribute the neurons randomly in a 3D space, respecting some constraints about their minimal and maximal mutual distances. We place an electrode with four virtual recordings sites at four different locations around the middle of the space. Distance between two neighboring sites is set to be equal to an average distance between two neighboring neurons. To emulate the effect of extracellular medium on propagating action potentials we simulate attenuation of the amplitude, but also a moderate waveform smoothing which arises from a lowpass filtering properties of the medium [START_REF] Bedard | Modeling extracellular field potentials and the Frequency-Filtering properties of extracellular space[END_REF]. We consider the smoothing between the recording sites to be negligible and simulate only the smoothing which happens during the action potential propagation between neurons and the electrode. We assume that the action potentials fired by the same neuron are always recorded as spikes of the same shape. Moreover, since the smoothing is generally moderate, we assume that the 1000 simulated neurons exhibit only 32 different spike shapes.. The 32 shapes are in fact templates generated by averaging spikes coming from 32 different neurons, recorded by real extracellular recordings (the recordings will be described in the next subsection). The length of each template is 60 samples, corresponding to 4ms real time. For each neuron we first generate a random firing rate between 5Hz and 10Hz and then generate Poisson distributed firing times.

To make the simulated signal more realistic, we have to simulate the sampling effects. Indeed, in real recordings spikes coming from a particular neuron generally look different, because they are sampled at different time instants with respect to the beginning of each spike. We simulate this effect by upsampling the spike templates by a factor of 100. Then we associate the spike templates with the corresponding time locations in the signal. Finally we subsample back to the original recording rate.

Assuming an isotropic setting, the extracellular potential typically decays as 1/r k , where r is the Euclidian distance between the firing neuron and the recording site. Depending on the position and orientation, the power k can take the value 1 (monopole model, near the soma), 2 or higher (dipole model or multipole model, far-field). For sake of simplicity, we consider k = 1 in our simulation. The first column on figure 2 shows an example of the simulated recording obtained in this way. We refer to [START_REF] Harris | Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements[END_REF] and [START_REF] Pettersen | Extracellular spikes and current-source density[END_REF] for further discussions on the propagation of action potentials through the extracellular medium.

Real neural recording

A drawback of the validation using the simulated signal lies in the complexity of a real neural mixture, which is very hard to reproduce artificially. Real recordings, depending on the tissue in which the recording is taken, consist of superposed activity of up to millions of neurons. The neurons interact through the large number of dendrites and axons, what makes the spatial configuration extremely complex. Thus, the simulated signal is often not a truthful counterpart of the real signal.

Apart from the simulated signal we examine the proposed algorithm using the real neural recordings. We use the simultaneous intracellular-extracellular recordings from a rat hippocampal area CA1 done by Henze et al., and described in [START_REF] Henze | Intracellular features predicted by extracellular recordings in the hippocampus in vivo[END_REF]. When the extracellular electrode is close to the intracellular, the activity of the neuron recorded intracellularly will be also visible in the extracellular recording. Of course, with the extracellular electrode activities of a vast number of neurons are recorded, so to separate the activity of individual neurons we have to do spike sorting. Figure 1 shows a part of such simultaneous recording. The four top plots show the extracellular recording with tetrode and the bottom plot shows the corresponding intracellular recording. The action potentials fired by the neuron recorded intracellularly are clearly visible on all four extracellularly recorded channels. This indicates that the recording sites are located very close to each other. On the plots of the extracellular channels we can also see the activity of other neurons that are around the extracellular electrode. Notice that many spikes from the other neurons have amplitudes similar to those of the spikes that come from the neuron recorded intracellularly. Thus, the task of spike sorting is obviously not trivial in this case.

The algorithm description

In this section we describe the proposed algorithm for multi-channel spike sorting. The algorithm iteractively removes the activities of distant neurons from the original recording as well as the activities of the neurons close to the electrode, once they are separated from the mixture. By doing so we reduce the number of sources, which brings the application of ICA in a more comfortable setting.

Apart from the recorded signal, to use the algorithm we only need to provide a lower bound for the firing rate, call it G, and the spike detection threshold level. This lower bound is usually easy to set since, generally, we are not interested in separating the activities of neurons that fired only a few spikes. A more detailed discussion on how to chose G is given in our previous paper on spike sorting [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF]. The choice of the threshold level, on the other side, is not very critical, since we will simply neglect all the spikes that belong to the cluster which has the lowest average peak-to-peak amplitude of the spikes. This will make the algorithm more reliable, since the spikes with the lowest amplitudes usually come from many different neurons, whose activity is only partially detected.

We assume that a four-channel recording is given, but the algorithm is the same for any number of channels. The four recorded signals are labeled as E1, E2, E3 and E4. An example of a part of the simulated four-channel recording is shown on figure 2 and of the four-channel real recording on figure 3. Spikes from the neuron recorded intracellularly are marked with the red stars. . 5000 samples from the simulated four-channel extracellular recording (first column). 2D feature vector space with the extracted spikes -positive peak amplitude vs negative peak amplitude (second column). Results of sorting algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] are given in the third column (only one cluster was detected for each channel).

To demonstrate that spike sorting is very difficult for such recording we apply directly the greedy based spike sorting method described in [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF]. First, we detect the spikes using the spike detection method described in [START_REF] Tiganj | Spike detection and sorting: Combining algebraic differentiations with ica[END_REF] and [START_REF] Mboup | A Volterra filter for neuronal spike detection[END_REF]. We then project Figure 3. 5000 samples from the real four-channel extracellular recording (first column). 2D feature vector space with the extracted spikes -positive peak amplitude vs negative peak amplitude (second column). Spikes that come from the neuron recorded intracellularly are labeled with the red stars, while the rest of the spikes are labeled with the blue dots. Results of sorting algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] are given in the third column, where the black dots marks are the detected centers of the clusters.

the spikes in a 2D vector space by keeping only two features per spike: the positive and the negative peak amplitudes. Notice that one could also use, e.g., the two first principal components. In the particular case, better (generally more reliable) results were obtained when the positive and the negative peak amplitude were used as the features, instead of the commonly applied PCA. The plots of the features are shown next to the corresponding signals, in the second column on figures 2 and 3. In the third column on figures 2 and 3 we give the outputs of the algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF]. The spikes that are left unsorted by the algorithm are automatically eliminated and not plotted. Note that the colors/symbols chosen for the representation of the clusters for the plots in the second column are not related to the one in the third column. It is evident that from such feature vector space it is not possible to do an accurate spike sorting.

We will now demonstrate the performance of the proposed algorithm using stepby-step description. The flowchart that describes the proposed algorithm is shown on figure 4.

• Step 1: ICA. We process the input signal with the FastICA algorithm described

ICA E1 E2 E3 E4
Remove everything but spikes from E1, E2, E3 and E4. Set:

T rackOneN euron=false. IC1 IC2 IC3 IC4 ICA E1 * E2 * E3 * E4 *
Find which IC is the most suitable for the further processing.

IC1 * IC2 * IC3 * IC4 * IC *
The IC * contains at least two clusters?

Remove spikes from the cluster which is the furthest from the cluster that contains the spikes with the largest amplitudes; Track-OneNeuron=true. in [START_REF] Hyvärinen | Independent component analysis: algorithms and applications[END_REF] (we used the default parameters given within the FastICA MATLAB package) and obtain IC1, IC2, IC3 and IC4 as the outputs (figures 5 and 6). The second and the third column on the both figures are again shown only to give an illustration of what the results would be if the feature extraction based sorting method was directly applied on the outputs of the ICA algorithm. The difference between the clusters is more obvious than on figures 2 and 3.

Yes

• Step 2: Remove everything but spikes from E1, E2, E3 and E4. Set: k=false. We first want to find which one of the four ICs contains the largest projection of the spikes from neurons close to the electrode. The criterion we set to determine this IC is spike dynamics. As a measure of the spike dynamics we use the average difference between the positive and the negative peak. Since spikes naturally contain a large peak, such difference should be a good measure of the dynamics. We first find all the spikes from all four IC using the method described in [START_REF] Tiganj | Spike detection and sorting: Combining algebraic differentiations with ica[END_REF] and [START_REF] Mboup | A Volterra filter for neuronal spike detection[END_REF]. Then, we calculate the corresponding spike dynamics for each IC and consider as actual spikes only those detected from the IC with the largest spike dynamics. Next, from all the 4 input recordings we remove (set to zero) all but the actual spikes. We assume that the removed parts of the signal correspond only to the noise and the cumulative activity of the neurons that are distant from the electrode. The outputs are called E1 * , E2 * , E3 * and E4 * .

• Step 3: ICA. We apply ICA on E1 * , E2 * , E3 * and E4 * . The resulting ICs, called IC1 * , IC2 * , IC3 * and IC4 * , are shown on figure 7 for the simulated signal and on figure 8 for the real signal. To demonstrate the results as before, we plot the spikes in the feature vector space (second column) and output of the greedy sorting algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] (third column). As the number of the sources is now lower, the ICA algorithm performed a little better and a little bigger difference between the clusters is visible, in comparison with the results from figures 5 and 6. • Step 4: Find which IC is the most suitable for the further processing.

In the further steps we will need only one IC. Thus in this step we want to chose the one which contains well preserved spike waveforms. It is likely that a good choice would be the one for which the detected spikes have the largest dynamics. We calculate the spike dynamics in the same way as in the step 2. For the further processing we chose the IC which has the largest value of the dynamics and call it IC * .

• Condition 1: The IC * contains at least two clusters?

• If Condition 1 is true: Remove spikes from the cluster which is the furthest from the cluster that contains the spikes with the largest amplitudes; TrackOneNeuron=true. In this step we have to apply the spike sorting algorithm proposed in [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF], this time it is not only for demonstration, but to find which spikes we can remove from the signal. For each cluster we compute the spike dynamics as the average of peak-to-peak amplitudes from the first G spikes. The spikes that have been classified in the cluster furthest from the one with the largest spike dynamics are now removed from E1 * , E2 * , E3 * and E4 * . Next, the algorithm goes back to step Result of applying ICA on the simulated recording after all but spikes was removed from the signal (set to zero) -5000 samples are displayed. The four ICs are shown in the first column. 2D feature vector spaces with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] on the vector space from the second column. IC that has been chosen for the further processing is marked with a red rectangle. The black triangle, with the tip pointing up, marks the center of the cluster that contains the spikes with the largest dynamics. Opposite, the black triangle with the tip pointing down marks the center of the cluster that contains the spikes with the smallest dynamics.

3. We show on figures 9 and 10 the outputs after applying ICA on the new E1 * , E2 * , E3 * and E4 * for the simulated and the real recording respectively. From the simulated recording (figure 9) we have detected only one cluster on each of four IC. We use the same criterion as before to chose which IC we should keep. Since the condition 1 is now false, the algorithm will continue to the condition 2. Notice that almost all the remaining spikes are fired by the same neuron. From the real recording (figure 10) we have detected more than one cluster on IC1 * and IC2 * ; for the further processing we will use IC1 * . The cluster that will be removed is the one labeled with a black triangle whose tip is pointing down. After that cluster is removed, spikes from the neuron recorded intracellularly will be almost the only one remaining in the signal. Only a Figure 8. Result of applying ICA on the real recording after all but spikes was removed from the signal (set to zero) -5000 samples are displayed. The four ICs are shown in the first column. 2D feature vector spaces with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from [START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] on the vector space from the second column. All the labels are the same as for figure 7.

few spikes from the intracellularly recorded neuron were removed. We set in this step the flag T rackOneN euron, in order to track if the algorithm converged just after a new neuron is separated, what would indicate generally that no more neurons can be separated from the given recording.

• If Condition 1 is false: Condition 2: TrackOneNeuron = true? • If Condition 2 is true: Neuron(i)=IC * ; Remove the spikes fired by the i th neuron; Set: i = i + 1. If T rackOneN euron is true that means that we have detected more than one cluster in the previous iteration so now, since only one cluster is left (the condition 1 was false), we assume that the remaining spikes are coming from the same neuron. This is the i th separated neuron. To continue our iterative algorithm we set to zero all the samples from the recording (from E1, E2, E3 and E4) when any of the already detected neurons (Neuron(1), Neuron(2),...,Neuron(i)) was active and go back to step 1. By doing this we removed an important high amplitude sources.

• If Condition 2 is false: Find the overlapped spikes and update Neuron. This is the last step of the proposed algorithm. We address now a very important problem in spike sorting: simultaneous firing. The proposed deflation based algorithm is obviously not able to detect simultaneous firing, since once a spike has been sorted the corresponding part of the signal is set to zero for the following iterations. Thus, if two or more neurons fired at approximately the same time, the spike will be assigned only to one of these neurons. To find the simultaneously fired spikes we analyze the activity of the sorted neurons two by two, for all possible pairs of the detected neurons. When analyzing the activity of each pair, we keep in the recording only the segments when one of the neurons in the pair fired a spike. Each such recording, call it si,j , where i = 1, ..., S -1 and j = i + 1, ..., S can be expressed ‡ as:

si,j = s i + s j + (1 -δ(|s i | + |s j |)) S k=1,k =i,j s k ,
where S is the number of sorted neurons and where δ(u) takes the value 1 for u = 0 and 0 elsewhere. Notice that we use absolute values of s i and s j to exclude the pairs i, j such that s i = -s j = 0. It is obvious that most of the energy in any si,j comes from s i and s j . Thus, applying ICA on si,j leads to the separation of s i and s j : one IC will represent mostly s i , another one s j and the remaining ICs will represent (1 -δ(|s i | + |s j |)) S k=1,k =i,j s k . Since the activity of s i and s j will be separated, any ‡ For notational simplicity, we drop the time argument and write s k instead of s k (t) for the recorded activity of neuron #k. Note that s k (t) is a sparse signal: We have s k (t) = 0 except for intervals [τ, τ + T ] when an action potential is emitted from time τ . The length T of the interval corresponds to the duration of an action potential.

spikes fired simultaneously by these two neurons should be visible on both of the ICs that represent the activity of these neurons. We demonstrate this by an example shown in figure 11. In the first two rows we plot part of the activity of two simulated neurons that fired three spikes simultaneously. The next four plots show four channels of the simulated recording. The activity of the two neurons is visible, but hard to recognize, especially to separate. We remove from the four channels everything but the segments when any of the two neurons was active and apply ICA on such sparse signal. The result of ICA is shown on four bottom plots on the same figure. Since this is the last step of the proposed algorithm, we already have the estimation of the activity of all the separated neurons. Thus, it is easy to find, from the result of ICA, that in this example IC2 represents the activity of Neuron 1 and that IC3 represents the activity of Neuron 2. All the overlapped spikes are visible on both IC2 and IC3. They are now very easy to localize and to update the detected activity of the separated neurons.

Results

We apply the algorithm described in the previous section on 5 real recordings and on 10 simulated recordings. The results are given in tables 1 and 2. There are four key-operations in the proposed algorithm: 1) removing of the noise (Noise rem.); 2) removing of the spikes from the clusters far from the one that contains the spikes which show the largest average peak-to-peak amplitude (Clus. rem.); 3) removing of the spikes from already separated neurons (Sorted rem); 4) detection of the overlapped spikes (Over. det.). We analyze contribution of each of these steps by simply comparing the proposed algorithm to the algorithms which perform none, only one or more of these steps. Namely, we compare the proposed algorithm with the following: (i) Only the spike sorting algorithm, without performing ICA (Spike sor.).

(ii) ICA and spike sorting (iii) ICA, spike sorting and removal of the spikes from already separated neurons.

(iv) As 3, but with removal of everything but spikes (step 2 of the algorithm).

(v) As 4, but with removal of clusters that generally contain spikes small in the amplitude (clusters which are far from the one that contains the spikes which show the largest average peak-to-peak amplitude).

The proposed algorithm is actually an extension of the algorithm (v), with a sequence for the detection of the overlapped spikes. We express the results through two parameters: 1) sorting accuracy (SA(%)) which we define as ratio of accurately sorted spikes and total number of detected spikes and 2) sorting detection (SD(%)) which we define as ratio of accurately sorted spikes and total number of spikes in the signal.

Results of the spike sorting on the real signal

For the real recordings we can quantify the results only for one neuron per recording, the one which was simultaneously recorded intracellularly. Thus SA and SD are very simple to calculate: SA(%) = 100 * C/(F + C) and SD(%) = 100 * C/T . Here C is the largest number of the spikes from the neuron recorded intracellularly, which are sorted into the same cluster (correct detections). F is the number of spikes from other neurons which are wrongly placed into the same cluster for which C is calculated (false spikes). T is the total number of spikes fired by the intracellularly recorded neuron.

In table 1 and on figure 12 we give SA and SD for each of the real recordings. Names of the recordings are kept from [START_REF] Henze | Intracellular features predicted by extracellular recordings in the hippocampus in vivo[END_REF]. Durations of the recordings were from 60 to 120 seconds with the sampling frequency of 10kHz for some recordings and 20kHz for others. The total number of the detected spikes per recording was from 800 to 2000, while the number of the spikes fired by the intracellularly recorded neuron was from 200 to 500.

The results suggest that the proposed approach improves significantly the sorting results in comparison with the basic spike sorting algorithm and the ICA algorithm. It generally decreases the number of falsely sorted spikes and increases the number of spikes that are left unsorted (missed spikes). As we mentioned in the introduction, in the trade-off between the number of missed spikes and the number of false spikes it is generally highly preferable to have more missed spikes.

When the intracellularly recorded neuron is sorted first (recordings d561102, d533101 and d533102) there is no benefit from the removal of the spikes from already sorted neurons (results in the third and the fourth column of table 1 are the same in such cases). Also, for those recordings contribution from detection of the overlapped spikes is very small, since most of the overlapped spikes were initially correctly assigned to the intracellularly recorded neuron (results in the last two columns of table 1 are the same/similar in such cases).

Results of the spike sorting on the simulated signal

Using simulated recordings we can quantify the results for more than one neuron and also we can easily perform several simulations and average their results. For each simulated recording we set the maximum number of outer loops (sorted neurons) to 6. However, in some situations (especially for simpler methods like the basic spike sorting and combination of the basic spike sorting with ICA) less than 6 neurons have been found.

Here we do a computation of SA and SD in a slightly different way than for the real signal. We compute the average SA and the average SD over 10 simulations for each j th sorted neuron, where j goes from 1 to 6, which is the maximum number of the sorted neurons:

SA j = N j i=1 SA i,j N j , SD j = N j i=1 SD i,j N j
where N j is the number of iterations in which at least j neurons were detected. We give also a standard deviation (sd) of SA and SD in order to quantify the robustness of each algorithm.

The results are given in table 2 and on figure 13. Even though the results vary significantly from neuron to neuron (spikes from the neurons close to the electrode are, as expected, sorted much more accurately that from the more distant ones), we can say again that overall each step of the proposed algorithm leads to significant improvement of the final results.

An important property of the proposed algorithm is that increasing SA is paid by decreasing SD as well, as is even more evident than on the figure 12. However, the last step (detection of overlapped spikes) increases SD, especially as j increases, since for the larger j (more distant neurons) many spikes were not detected because they overlapped spikes from the closer neurons (which correspond to the lower j).

The standard deviation is larger for the neurons further from the electrode and tends to decrease as a more complex algorithm is used. The standard deviation illustrates the

Discussion

In this paper we proposed a new algorithm for spike sorting from multi-channel recordings. We use iterative application of ICA and deflation to sort the spikes. Activity at the recording sites is mutually highly dependant, while the activity of single neurons is mutually much less dependant. This naturally suggests application of ICA to separate the neural mixture.

A condition that the number of sensors should be equal or greater than the number of sources appears as the major problem and limitation for application of ICA on neural recordings. The proposed algorithm uses deflation to eliminate or at least to reduce the contribution of different sources. In particular, we remove: 1) everything but spikes, 2) spikes which have already been sorted and 3) spikes which belong to the clusters far from the one which contains the spikes of the largest average peak-to-peak amplitude.

Such deflation would result in inability of the algorithm to detect spikes fired by different neurons at approximately the same time -simultaneous firing. Only one spike per simultaneous firing could be detected. The solution to this problem is implemented as the last step of the proposed algorithm. We take two of the separated neurons and keep in the recording only the segments in which one of these two neurons was active. Then we apply ICA on such sparse signal. Since we have only two high-energy sources in the signal, ICA will, in most situations, be able to separate the activity of the two neurons, including simultaneously fired spikes. We repeat this procedure for all the possible pairs of the separated neurons.

The results given in tables 1 and 2 suggest that each step of the algorithm improves the final sorting results. Overall, the algorithm significantly increases the number of correctly sorted spikes while increasing the number of missed spikes by a relatively small amount. ICA alone was not sufficient to obtain reliable results. It is important to mention that if a recording is very noisy, the proposed algorithm will not be able to detect any cluster.

In some particular cases, such as multi-electrode recordings from peripheral nerve interfaces, the spike waveforms can arrive to the different recording sites with time delays up to 0.1ms [START_REF] Citi | Characterization of tflife neural response for the control of a cybernetic hand[END_REF], [START_REF] Citi | On the use of wavelet denoising and spike sorting techniques to process electroneurographic signals recorded using intraneural electrodes[END_REF]. FastICA algorithm is not adapted for such mixtures. In its present form, the proposed algorithm is therefore not suitable for those cases. Nonetheless, the proposed deflation methodology would still be useful if FastICA is replaced by some convolutive version of ICA. In the convolutive case [START_REF] Thomas | Time-domain fast fixed-point algorithms for convolutive ICA[END_REF], [START_REF] Dyrholm | Model selection for convolutive ICA with an application to spatio-temporal analysis of EEG[END_REF], the recorded mixture is described by an FIR filter model of the mixing process. Such model could be beneficial in this particular case of spike sorting for two reasons: 1) Ability to account for different delays of the source components; 2) Ability to account for echoes coming from different parts of a nerve around the electrode.

As spike features we used simply positive and negative peak amplitude. However, we could use practically any feature extraction method. Methods that can make the difference between the spikes fired by different neurons more evident will course lead to better performance of the proposed algorithm, in terms of the number of sorted neurons and the sorting accuracy. Nevertheless, no matter how good the chosen feature extraction technique is, the proposed algorithm should always improve the final sorting results when they are compared to the results obtained with only ICA algorithm. Finally we can say that the obtained results suggest that the proposed algorithm gives better sorting results as compared to ICA.

Figure 1 .

 1 Figure 1. Simultaneous extracellular-intracellular recording. The four recording sites are very close, so the action potentials recorded with the intracellular electrode are clearly visible on all four extracellularly recorded channels.

Figure 2

 2 Figure2. 5000 samples from the simulated four-channel extracellular recording (first column). 2D feature vector space with the extracted spikes -positive peak amplitude vs negative peak amplitude (second column). Results of sorting algorithm from[START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] are given in the third column (only one cluster was detected for each channel).

TFigure 4 .

 4 Figure 4. Flowchart of the proposed multi-channel spike sorting algorithm.

Figure 5 .

 5 Figure 5. Result of applying ICA on the recording shown on figure 2. The four ICs are shown in the first column (5000 samples). 2D feature vector space with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from [23] on the vector space from the second column. The red rectangle marks the IC that has the largest spike dynamics.

Figure 6 .

 6 Figure 6. Result of applying ICA on the recording shown on figure 3. The four ICs are shown in the first column (5000 samples). 2D feature vector space with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from [23] on the vector space from the second column. The red rectangle marks the IC that has the largest spike dynamics.

Figure 7 .

 7 Figure 7. Result of applying ICA on the simulated recording after all but spikes was removed from the signal (set to zero) -5000 samples are displayed. The four ICs are shown in the first column. 2D feature vector spaces with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from[START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] on the vector space from the second column. IC that has been chosen for the further processing is marked with a red rectangle. The black triangle, with the tip pointing up, marks the center of the cluster that contains the spikes with the largest dynamics. Opposite, the black triangle with the tip pointing down marks the center of the cluster that contains the spikes with the smallest dynamics.

Figure 9 .

 9 Figure 9. Result of applying ICA on the simulated recording, after removing (setting to zero) of the low amplitude spikes (see figure7) from E1 * , E2 * , E3 * and E4 * (5000 samples). The four ICs are shown in the first column. 2D feature vector spaces with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from[START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] on the vector spaces from the second column. All the labels are the same as for figure7.

Figure 10 .

 10 Figure10. Result of applying ICA on the real recording, after removing (setting to zero) of the low amplitude spikes (see figure8) from E1 * , E2 * , E3 * and E4 * (5000 samples). The four ICs are shown in the first column. 2D feature vector spaces with the spikes extracted from the ICs are in the second column. The plots in the third column display the results of applying the algorithm from[START_REF] Tiganj | A non-parametric method for automatic neural spikes clustering based on the non-uniform distribution of the data[END_REF] on the vector spaces from the second column. All the labels are the same as for figure7.

Figure 11 .

 11 Figure 11. The first two plots show part of the activity of the two simulated neurons. The following four plots show the part of the four-channel simulated recording. The four plots at the bottom show the results of applying ICA on the four channels. Before applying ICA the recording was modified by setting to zero all but the segments where one of the two neurons was active. It appears that the activity of Neuron 1 is represented by IC3 and the activity of Neuron 2 is represented by IC1. The three arrows on the top of the figure point out the overlapped spikes from the two neurons. All of the three overlapped spikes are successfully sorted.

Figure 12 .

 12 Figure 12. Comparison of different spike sorting algorithms for 5 real recordings, in terms of sorting accuracy (SA) and sorting detection (SD).

Table 1 .

 1 Comparison of different spike sorting algorithms for 5 real recordings, in terms of sorting accuracy (SA) and sorting detection (SD).

Acknowledgments

This work was supported by the PHC COGITO project.