Marine Jeanjean

Krista Mcgrath

Silvia Valenzuela-Lamas

Ariadna Nieto-Espinet

Renate Schafberg

Pere Miquel Parés-Casanova

Sergio Jiménez-Manchón

Claude Guintard

Faiza Tekkouk

Rania Ridouh

Cyprien Mureau

Allowen Evin
email: allowen.evin@umontpellier.fr

ZooMS confirms geometric morphometrics species identification of ancient sheep and goat

Keywords: Caprines Identifications, GMM, Ovis aries, Capra hircus, Paleoproteins, Third lower molar

Geometric morphometrics can effectively distinguish isolated third lower molars of present-day sheep and goat, but its applicability to archaeological specimens has yet to be established. Using a modern reference collection of 743 sheep and goats and a 2D landmark based geometric morphometric (GMM) protocol, this study aimed to morphometrically identify 109 archaeological specimens, used as case studies, dating from the Late Neolithic to the modern period/era. These morphometric identifications were then compared to molecular identifications via collagen peptide mass fingerprinting, known as Zooarcheology by Mass Spectrometry (ZooMS). ZooMS confirmed the morphometric identifications for 104 specimens, with the five misidentified specimens all morphometrically identified as goat.

Modern sheep and goats have larger teeth and distinct shapes compared to their archaeological counterparts, suggesting strong differences between archaeological and modern specimens potentially linked with recent breed improvement or geographic origin of the specimens. In addition, for both species, some of the archaeological dental morphologies do not match with any of our modern references.

This study validates the applicability of geometric morphometrics for identifying isolated archaeological sheep and goat teeth. It represents a stepping stone for future, non-destructive, bioarchaeological studies of the two species.

INTRODUCTION

Zooarchaeology by mass spectrometry (ZooMS) is increasingly used in archaeology for reliable taxonomic identification of organic materials [START_REF] Richter | A primer for ZooMS applications in archaeology[END_REF] such as bones [START_REF] Buckley | Species identification by analysis of bone collagen using matrix-assisted laser desorption/ionisation time-of-flight mass spectrometry[END_REF], leather [START_REF] Ebsen | Identifying archaeological leatherdiscussing the potential of grain pattern analysis and zooarchaeology by mass spectrometry (ZooMS) through a case study involving medieval shoe parts Denmark[END_REF] or parchment [START_REF] Fiddyment | Animal origin of 13th-century uterine vellum revealed using noninvasive peptide fingerprinting[END_REF]. ZooMS, by using peptide mass fingerprinting of collagen, provides an inexpensive tool for archaeological studies and has been used to distinguish domestic caprine species and to detect their arrival in Africa [START_REF] Coutu | Palaeoproteomics confirm earliest domesticated sheep in southern Africa ca[END_REF][START_REF] Culley | Collagen fingerprinting traces the introduction of caprines to island Eastern Africa[END_REF][START_REF] Janzen | Distinguishing African bovids using Zooarchaeology by Mass Spectrometry (ZooMS): New peptide markers and insights into Iron Age economies in Zambia[END_REF][START_REF] Culley | Iron Age hunting and herding in coastal eastern Africa: ZooMS identification of domesticates and wild bovids at Panga ya Saidi, Kenya[END_REF], Asia [START_REF] Taylor | Early pastoral economies along the Ancient Silk Road: biomolecular evidence from the Alay Valley, Kyrgyzstan[END_REF][START_REF] Taylor | Evidence for early dispersal of domestic sheep into Central Asia[END_REF], and Europe [START_REF] Martins | Radiocarbon dating the beginning of the Neolithic in Iberia: New results, new problems[END_REF][START_REF] Sjögren | Evaluating caprine remains of the Swedish Funnel Beaker culture through ZooMS[END_REF], and to study distinct farming practices [START_REF] Buckley | Collagen fingerprinting of archaeological bone and teeth remains from Domuztepe, South Eastern Turkey[END_REF][START_REF] Price | Animal Management Strategies during the Chalcolithic in the Lower Galilee: New Data from Marj Rabba (Israel)[END_REF][START_REF] Vaiglova | An integrated stable isotope study of plants and animals from Kouphovouno, southern Greece: A new look at Neolithic farming[END_REF][START_REF] Birch | Combined osteomorphological, isotopic, aDNA, and ZooMS analyses of sheep and goat remains from Neolithic Ulucak, Turkey[END_REF][START_REF] Prendergast | Sorting the sheep from the goats in the Pastoral Neolithic: morphological and biomolecular approaches at Luxmanda, Tanzania[END_REF][START_REF] Gron | Separating caprine (Capra/Ovis) distal tibiae: A case study from the Polish Neolithic[END_REF][START_REF] Ananyevskaya | Specialized wool production economy of prehistoric farmstead of Chap I in the highlands of Central Tian Shan (Kyrgyzstan)[END_REF][START_REF] Bläuer | Multidisciplinary History of Goats in Finland: A Comparative Approach[END_REF].

Sheep (Ovis aries) and goat (Capra hircus) are two important species of agro-pastoral systems in the Mediterranean basin since the Neolithic and were often herded together [START_REF] Boyazoglu | Livestock research and environmental sustainability with special reference to the Mediterranean basin[END_REF][START_REF] Galaty | Pastoral Systems[END_REF][START_REF] Enne | The role of grazing in agropastoral systems in the Mediterranean region and their environmental sustainability[END_REF][START_REF] Durmus | Sustainability of Small Ruminant Production in Mediterranean Region[END_REF][START_REF] Marandet | L'homme et l'animal dans les sociétés méditerranéennes[END_REF][START_REF] Nieto-Espinet | Resilience and livestock adaptations to demographic growth and technological change: A diachronic perspective from the Late Bronze Age to Late Antiquity in NE Iberia[END_REF][START_REF] Valenzuela-Lamas | Animal husbandry across the Western Roman Empire: Changes and continuities[END_REF][START_REF] Halstead | Pastoralism or household herding? problems of scale and specialization in early Greek animal husbandry[END_REF][START_REF] Trentacoste | Systems change: Investigating climatic and environmental impacts on livestock production in lowland Italy between the Bronze Age and Late Antiquity (c. 1700 BC -AD 700)[END_REF]. However, they possess distinct ecological [START_REF] Shackleton | A Review of the Social Behavior of Feral and Wild Sheep and Goats[END_REF][START_REF] Balasse | Distinguishing sheep and goats using dental morphology and stable carbon isotopes in C4 grassland environments[END_REF] and economical properties [START_REF] Serre | Le théatre d'agriculture et ménage des champs[END_REF][START_REF] Bourrier | Les industries des abattoirs[END_REF][34][START_REF] Helmer | L'élevage des caprinés néolithiques dans le sud-est de la France : saisonnalité des abattages, relations entre grottes-bergeries et sites de plein air[END_REF][START_REF] Blaise | Economie animale et gestion des troupeaux au Néolithique final en Provence[END_REF], making it of prime interest to establish secure identifications of ancient specimens before studying the two species separately. Criteria for distinguishing sheep and goat among archaeological assemblages exist, notably from post-cranial bones [START_REF] Boessneck | Osteological Differences between Sheep (Ovis aries Linné) and Goat (Capra hircus Linné)[END_REF][START_REF] Zeder | Assessing the reliability of criteria used to identify postcranial bones in sheep, Ovis, and goats, Capra[END_REF][START_REF] Prummel | A guide for the distinction of species, sex and body side in bones of sheep and goat[END_REF] and have greatly improved the ability to differentiate the two species. On teeth, discrete criteria have been proposed [START_REF] Balasse | Distinguishing sheep and goats using dental morphology and stable carbon isotopes in C4 grassland environments[END_REF][START_REF] Halstead | Sorting the sheep from the goats: Morphological distinctions between the mandibles and mandibular teeth of adult Ovis and Capra[END_REF][START_REF] Helmer | Discrimination des genres Ovis et Capra a l'aide des premolaires inferieures 3 et 4 et interpretation des ages d'abattage: l'example de dikili tash (Grece)[END_REF][START_REF] Payne | Morphological distinctions between the mandibular teeth of young sheep, Ovis, and goats, Capra[END_REF][START_REF] Zeder | Assessing the reliability of criteria used to identify mandibles and mandibular teeth in sheep, Ovis, and goats, Capra[END_REF][START_REF] Gillis | An assessment of morphological criteria for discriminating sheep and goat mandibles on a large prehistoric archaeological assemblage (Kerma, Sudan)[END_REF][START_REF] Grine | Analysis of enamel ultrastructure in archaeology: The identification of Ovis aries and Capra hircus dental remains[END_REF], but none has been universally used for identifying isolated lower third molars that are often find among archaeological assemblages, since, unlike bones, they are better preserved [START_REF] Lyman | Bone density and differential survivorship of fossil classes[END_REF], but also often left unidentified because of the lack of secure identification criteria.

Geometric morphometrics (GMM) is more and more used in bioarchaeology [START_REF] Evin | Archaeophenomics of ancient domestic plants and animals using geometric morphometrics : a review[END_REF]. This series of analytical tools are based on the multivariate statistical analysis of set of cartesian coordinate data allowing to explore in-depth morphometric variation [START_REF] Rohlf | A Revolution in Morphometrics[END_REF]. Recently, such approaches have been used to quantify sheep and goat morphometric variation based on the size and shape of their lower third molar [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. This study demonstrated the ability to identify adult modern specimens of the two species with a probability of 93.3%. In addition, the identification of Middle Ages specimens from the site of Missignac (Aimargues, southern France), used as a first case study, appeared congruent with the initial zooarchaeological faunal analyses and interpretations [START_REF] Maufras | Genèse, évolution et désertion de Missignac (Aimargues, Gard), villa des Ve-XIIIe siècles. Archéologie du Midi Médiéval[END_REF]51]. However, in order to assess whether the same geometric morphometrics (GMM) protocol can be used to identify ancient sheep and goat teeth originating from a broader spatio-temporal scale, it is necessary to validate the approach using securely identified archaeological specimens. Archaeological teeth can be found either enclosed in mandibular bone that possess distinct features between the two species [START_REF] Halstead | Sorting the sheep from the goats: Morphological distinctions between the mandibles and mandibular teeth of adult Ovis and Capra[END_REF][START_REF] Zeder | Assessing the reliability of criteria used to identify mandibles and mandibular teeth in sheep, Ovis, and goats, Capra[END_REF], or isolated (i.e. without the bone), which make their identification particularly challenging even if some discrete criteria have been proposed for both mandible and teeth [START_REF] Balasse | Distinguishing sheep and goats using dental morphology and stable carbon isotopes in C4 grassland environments[END_REF][START_REF] Halstead | Sorting the sheep from the goats: Morphological distinctions between the mandibles and mandibular teeth of adult Ovis and Capra[END_REF][START_REF] Helmer | Discrimination des genres Ovis et Capra a l'aide des premolaires inferieures 3 et 4 et interpretation des ages d'abattage: l'example de dikili tash (Grece)[END_REF][START_REF] Payne | Morphological distinctions between the mandibular teeth of young sheep, Ovis, and goats, Capra[END_REF][START_REF] Zeder | Assessing the reliability of criteria used to identify mandibles and mandibular teeth in sheep, Ovis, and goats, Capra[END_REF][START_REF] Gillis | An assessment of morphological criteria for discriminating sheep and goat mandibles on a large prehistoric archaeological assemblage (Kerma, Sudan)[END_REF][START_REF] Grine | Analysis of enamel ultrastructure in archaeology: The identification of Ovis aries and Capra hircus dental remains[END_REF]. Being able to provide reliable identification is therefore of prime interest to bioarchaeological studies.

Here we aim to (1) identify archaeological specimens based on their third molar shape using geometric morphometrics (GMM), (2) compare those morphometric identifications with newly generated molecular Mass Spectrometry (ZooMS) identifications, and (3) compare the size and shape morphometric variation of modern sheep and goat with their securely identified archaeological counterparts.

MATERIAL AND METHODS

Third lower molars of 743 modern specimens, corresponding to 521 sheep and 222 goats of various breeds and geographical origins from Western Eurasia, including western Europe and north Africa, were used as a reference for the geometric morphometrics analysis (SI-table 1). Part of these specimens were already included in a previous study [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. Only one tooth per specimen was considered. All studied specimens were adults, with fully erupted third lower molars. Tooth wear stages were used to estimate the age at death of the specimens following the classification established by [START_REF] Payne | Kill-off Patterns in Sheep and Goats : The Mandibles from Aşvan Kale[END_REF][START_REF] Payne | Morphological distinctions between the mandibular teeth of young sheep, Ovis, and goats, Capra[END_REF] and they were ranging from stages 0 to 17G (i.e. from 1-2 years to more than 8 years). Sheep and goat display distinct sexual dimorphism and changes through age based on their third molar morphometry, but both factors were found to have little impact on the species differenciation [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF].

In addition, 109 individuals (SI-table 2) from 13 archaeological sites, selected for their chronological and spatial diversity, dated from the Neolithic to Modern period and located in the north-western part of the Mediterranean basin (present-day Catalonia region in Spain and Occitanie region in Southern France), were used to compare GMM and ZooMS identifications (figure 1, table 1). The selected archaeological specimens consisted of third lower molars with at least a small portion of adherent mandibular bone. For each specimen, the size and shape of the third lower molar was quantified through geometric morphometrics and used for morphometric identification, and a fragment of mandibular bone (normally between 10 to 50 mg) was sampled for molecular ZooMS identification.

Geometric morphometrics

The occlusal view of the third lower molars was photographed using a Nikon d90 LSR camera paired with a 60 mm macro lens (AF-S Micro NIKKOR) attached to a photographic arm (manfrotto 244RC) following a previously published protocol [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. Right molars were preferentially photographed, and when necessary, photographs of left molars were mirrored prior to GMM data acquisition. Third molar size and shape were quantified using two dimensional landmarks and sliding landmarks based geometric morphometrics approaches following the protocol developed in our previous study [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF] (figure 2). The coordinates of seven landmarks and of 48 sliding semi-landmarks distributed along six curves (eight equidistant points in each) were acquired using TPSdig2 [START_REF] Rohlf | [END_REF] (for a formal description of the landmark positions see SI table 2 of [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]) (figure 2). This protocol does not include points along the mesial part of the teeth that is often affected by lateral wear [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. All data are available in SI table 3.

Figure. 2

Geometric morphometrics protocol: Position of the seven landmarks (in blue) and 48 sliding semi-landmarks (in red) whose coordinates were acquired on a right third lower molar of a sheep. Landmark positions are described in SI table 2 of [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF].

Morphometrics identification

All modern and archaeological specimens were superimposed together using a Generalized Procrustes Analysis (GPA) during which the sliding semi-landmarks were allowed to slide by minimizing the sum of the Procrustes distances between each individual and the mean conformation [START_REF] Perez | Differences between sliding semi-landmark methods in geometric morphometrics, with an application to human craniofacial and dental variation[END_REF][START_REF] Sampson | Eigenshape analysis of left ventricular outlines from contrast ventriculograms[END_REF][START_REF] Sheets | A combined landmark and outline-based approach to ontogenetic shape change in the Ordovician trilobite Triarthrus becki[END_REF]. The Procrustes residuals were first analysed using a Principal Component Analysis (PCA) [START_REF] Gower | Generalized procrustes analysis[END_REF][START_REF] Rohlf | Extensions of the procrustes method for the optimal superimposition of landmarks[END_REF][START_REF] Goodall | Procrustes methods in the statistical analysis of shape revisited[END_REF] (SI figure 1). A predictive linear discriminant analysis (pLDA), computed on the 20 first PCA scores that maximise the differences between groups [START_REF] Baylac | Fourier descriptors, Procrustes superimposition, and data dimensionality: an example of cranial shape analysis in modern human populations[END_REF], was then used to predict to which species the archaeological specimens belonged. Because unbalanced sample size can have a profound effect on discriminant analyses [START_REF] Kovarovic | Discriminant function analyses in archaeology: Are classification rates too good to be true?[END_REF], and the fact that sheep greatly outnumber goat in the reference dataset, we used a resampling procedure to down sample the number of sheep to the smaller number of goats, as recommended by previous analysis [START_REF] Evin | The long and winding road : identifying pig domestication through molar size and shape[END_REF]. The morphometrics identification was therefore based on 100 predictive discriminant analyses (pLDA) performed on two equal size samples (here 222 in each group) whose results were summarized by the percentage of times a specimen was identified to its correct group [START_REF] Evin | Unravelling the complexity of domestication: A case study using morphometrics and ancient DNA analyses of archaeological pigs from Romania[END_REF]. All discriminant analyses were retained and no threshold was fixed for posterior probabilities (i.e. all identifications were based on the 100 pLDA). Only LDA with a cross validation percentage (CVP) above the 3 rd quartile of all CVP were retained in order to select the most discriminant analyses. Morphometric identifications were only based on molar shape data, and not size, since it has been previously identified as the most discriminative criteria between sheep and goat [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF].

ZooMS

Between 10 and 50 mg of bone was weighed out for each specimen. 250 µL of 0.6 M hydrochloric acid (HCl) was added to each sample then stored at 4°C until demineralised. Once demineralised, the acid was discarded and three washes of 200 µL each of ammonium bicarbonate (AmBic, NH₄HCO₃, pH8) were performed to remove any remaining acid. 100 µL of AmBic was added and the samples were gelatinised at 65°C for 1 hour. Samples were centrifuged for 1 minute and 50 µL of the gelatinised supernatant was transferred to a new tube and 0.4 µg/µL of trypsin was added. The trypsinated samples were digested at 37°C overnight, after which the samples were acidified with 1 µL of 5% trifluoroacetic acid to stop the trypsin. C18 zip tips (Pierce™, Thermo Scientific™) were used to purify the samples which were eluted in a 50 µL volume.

Comparison of ancient and modern sheep and goat

The overall percentage of correct identification was computed comparing the geometric morphometrics and ZooMS identifications. Securely ZooMS-identified archaeological sheep and goat were then compared with their modern counterparts using Wilcox-tests for size and Procrustes Anova (procD.lm function from geomorph package [START_REF] Adams | Geomorph: Software for geometric morphometric analyses[END_REF]) for shape. Size and shape variation were visualised using boxplot and PCA, respectively. Homogeneity of size and shape difference between species and between archaeological vs. modern specimens were tested using 2-way Anova and 2-way Procrustes Anova, respectively.

Morphometrics proximities between modern sheep and goat and the archaeological specimens identified via ZooMS were visualized by a neighbour-joining network computed on Mahalanobis distances. Graphics were created in R using the 'ggplot2' package [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF] Change in molar shape and size variance between archaeological and modern specimens of sheep and goat was assessed, respectively, using the disparity test of the geomorph R package [START_REF] Collyer | RRPP: An R package for fitting linear models to high-dimensional data using residual randomization[END_REF][START_REF] Collyer | RRPP: Linear Model Evaluation with Randomized Residuals in a Permutation Procedure[END_REF] and Fligner test (stats package, [START_REF] Core | R: a language and environment for statistical computing[END_REF]).

All the analysis were performed in R version 6.1.524 [START_REF] Core | R: a language and environment for statistical computing[END_REF] through Rstudio version 4.2.2 [START_REF]RStudio: Integrated Development Environment for R[END_REF].

RESULTS

Identification for archaeological specimens

Based on their third molar shape, 59 archaeological specimens were identified as sheep (54.1%) and 50 as goat (45.9%) (figure 3).

Morphometrics comparison with modern specimens

The ancient samples studied were selected to represent a diversity of geographic origin (Catalonia and the South of France) and chrono-cultural periods (Late Neolithic to modern period) but certainly do not cover the entire past sheep and goat diversity in the North Western Mediterranean basin. Sample sizes per region and chrono-cultural period were too small to envision detailed morphometrics comparison and all archaeological specimens were analysed jointly (See SI-figure-2 for a visualization of size differences between sites and chrono-cultural periods). Overall, the archaeological sheep and goat possessed smaller third lower molars than their modern relatives (figure 3, table 2). Modern sheep also appear larger than goat, contrary to the pattern observed for archaeological specimens (figure 3, table 2). A PCA of all archaeological and modern sheep and goat shows that the modern specimens studied show new morphologies that didn't exist in the past (figure 4). Shape differences between sheep and goat are present for both modern and archaeological specimens (respectively p-value=0.013 and p-value=0.001). Mean shape visualisation and dissimilarity network on shape confirm that there is a higher proximity within species than within periods, that is to say the archaeological goat (identified through ZooMS) reassemble more to modern goat than to archaeological sheep (and conversely) (figure 5), even if there are shape differences between archaeological and modern sheep (p-value= 0.001) and goat (p-value=0.004).

In addition, differences between sheep and goat are homogeneous between modern and archaeological specimens for both shape (interaction term of a 2-way Procrustes ANOVA, p-value= 0.463) and size (p-value= 0.138)

The only difference in terms of size or shape variance is seen between modern sheep and goat, with modern sheep showing more variation than modern goat (table 3).

DISCUSSION

We demonstrate, for the first time, that geometric morphometrics can be an efficient alternative to molecular identification with 95.4% of the morphometric identifications being congruent with those provide by the molecular ZooMS identification.

Identification of paleontological or archaeological remains is often done by actualism, defined as 'the methodology of inferring the nature of past events by analogy with processes observable and in action in the present' [START_REF] Rudwick | Uniformity and progression: reflections on the structure of geological theory in the age of Lyell[END_REF]. In order to fully explain the morphological diversity in the past, a better understanding of current diversity and its incorporation in the analyses is needed [START_REF] Kitts | The Structure of Geology[END_REF][START_REF] Pobiner | Applying Actualism: Considerations for Future[END_REF]. The hypothesis based on these identifications assumed that modern specimens are relevant for studying remains from the past, with negligible influence of extrinsic (e.g. temporal or geographical) variation. Modern specimens are often the only available reference of securely identified species, however, it has been demonstrated that modern references can lead to misidentifications. For example, modern dog breeds are not a good proxy to identify archaeological specimens [START_REF] Brassard | Unexpected morphological diversity in ancient dogs compared to modern relatives[END_REF]. This is also the case for other domestic species for which it is well known that direct selection over the last two centuries was intensified, followed by standardization of morphology and performance with the emergence of the breed concept [START_REF] Taberlet | Are cattle, sheep, and goats endangered species?[END_REF][START_REF] Taberlet | Conservation genetics of cattle, sheep, and goats[END_REF][START_REF] Denis | Races Bovins: Histoire, Aptitudes, Situation actuelle[END_REF][START_REF] Denis | Les Races d'animaux domestiques en France[END_REF]. This also applies to sheep and goat, whose breeds became genetically more and more uniformized [START_REF] Taberlet | Conservation genetics of cattle, sheep, and goats[END_REF], even if our morphometric results do not detect increase in diversity between our ancient and modern datasets. However, we detected morphologies that exist today that were not present in the studied archaeological dataset. This result will have to be explored further using more archaeological specimens.

We found that a large majority (95.4%) of the specimens were correctly identified based on their third molar shape, suggesting that the morphometric integrity, and thus the distinction of the species, are maintained through time, at least for the third lower molar. This means that, based on their third molar shape, the differences between sheep and goat remain the same through time. On the other hand, modern specimens of both species possess larger teeth than their archaeological counterparts with a strong overlap between sheep and goat, confirming the impossibility of using size alone for species identification. This size difference between archaeological and modern specimens, although little studied, has also been suggested for cattle [START_REF] Albarella | Shape variation of cattle metapodials: age, sex or breed? Some examples from Medieval and Postmedieval sites[END_REF] and pigs [START_REF] White | 2014 Shipping sheep or creating cattle: Domesticate size changes with Greek colonisation in Magna Graecia[END_REF].

The five misaligned specimens, those identified as sheep by ZooMS and goat by geometric morphometrics suggest that some ancient sheep dental morphologies were more similar to the herecaptured present-day goat dental morphologies. It seems important to note that the percentage of correct identification of 95.4% observed here is similar to the 93.3 % percentage of correct cross-validation between modern sheep and goat obtained in a previous study [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. At least four hypotheses could be envisaged to explain the misidentified specimens.

First, even if our modern reference collection includes 521 sheep and 222 goats belonging to a range of breeds originating from Europe and North Africa, the entire existing modern diversity is certainly not represented. Furthermore, this reference collection was purposely built to study the north western Mediterranean region, and so may not be directly suitable for other geographic areas. Similarly, two waves of exotic foods are known during Iron Age [100], the Roman and medieval periods [101,102] that may have influenced the diversity we observe in our archaeological dataset. Consequently, we cannot exclude that these five individuals may reflect animal translocation across regions, from areas where sheep dental morphology was similar to that of present-day goats.

Second, these animals could correspond to taxa not included in our analysese.g. the mouflon (Ovis orientalis musimon, a feral form of sheep ), the alpine (Capra ibex) or Iberian ibex (Capra pyrenaica) or the Barbary sheep (Ammotragus lervia)-whose dental morphology and geographical distributions in the past are largely unknown, and whose distinction based on collagen peptide sequences is limited. Due to the highly conserved nature of the amino acid sequence of collagen, ZooMS cannot be used to distinguish between domestic and wild specimens of the same species, and sometimes cannot distinguish between closely related taxa, as their collagen sequences are identical or near identical [START_REF] Buckley | Distinguishing between archaeological sheep and goat bones using a single collagen peptide[END_REF]. For example, domestic sheep, mouflon and Barbary sheep cannot be distinguished from another using ZooMS as their collagen peptide markers are the same. It is possible that these five samples belong to taxa that are closely related to domestic sheep however have dentition more similar to that of domesticated goats. This hypothesis could be tested in the future through e.g. palaeogenomics.

Third, age and sex are known to have a significant impact on sheep and goat third molar morphology, but that both factors have limited impact on the ability to discriminate the two species [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. Indeed, even if age (estimated through wear stages) affects the third lower molar shape of both species, a closer proximity between species than between age categories has been found [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. In addition, it was demonstrated that the impact of age on the between species differenciation was mostly due to the oldest specimens (more than 8 years) [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF]. Because the mis-indentitied specimens had an estimated age of 4-8 years, it seems unlikely that age caused the mismatch between the morphometric and molecular identifications. Sexual dimorphism in teeth size and shape has been identified in both species; but as for age, a greater proximity has been found between species than sex, and the sexual dimorphism has no impact on modern caprine identification [START_REF] Jeanjean | Sorting the flock: Quantitative identification of sheep and goat from isolated third lower molars and mandibles through geometric morphometrics[END_REF].

Finally, even if very rare, sheep-goat hybrids exist (e.g. [103,104]). Such hybrids would likely be difficult, if not impossible, to identify using either ZooMS or GMM, as presumably the results for both analyses would simply indicate either sheep or goat, rather than a mix.

Nevertheless, with a ~95% success rate, this newly generated data allows to build a secure reference collection of ancient morphologies for the two species from the Neolithic to Modern periods in the North Western Mediterranean basin. This will enable non-molecular, and thus non-destructive, methods of analysis of dentition for larger zooarchaeological assemblages to improve our understanding of the husbandry history of these two species in space and time. Because this study aims only at validating the methodological approach, the archaeological interpretations require further consideration. While the number of specimens studied here did not allow detailed analyses, adding more specimens will allow further exploration of the environmental and socio-cultural factors that shaped the past agrobiodiversity and its evolution. While it is known that modern [START_REF] Balasse | Distinguishing sheep and goats using dental morphology and stable carbon isotopes in C4 grassland environments[END_REF] and ancient goat and sheep have different diet[105], the impact of diet in tooth morphology in those species remain to be explored, but in any case did not have a major impact on the taxonomic identification.

The observed size and shape differences between archaeological and modern specimens of both sheep and goats are in agreement with previous zooarchaeological results. Here we observed that archaeological specimens have on average smaller molars than their modern counterparts. Sheep and goat are known to have increased in size during Modern period [106-108]. Caprines were particularly improved upon over the last centuries [108-110], which could explain the differences noticed between archaeological and modern groups. However, despite the presence of both primitive and standardised breeds in our referential, we did not find that modern specimens were more or less diverse than the archaeological groups. The only observed variance difference is between modern sheep and goat, where modern sheep are, for both size and shape, more variable than modern goat. This result could be due to the greater number of sheep compared to goat in the modern referential. Moreover, today officially recognised sheep breeds are also more numerous than goat, with for example around 50 different French sheep breeds compared to only 15 breeds of goat [START_REF] Denis | Les Races d'animaux domestiques en France[END_REF]111]. Moreover goats were mainly selected for the purpose of milk, while sheep have been selected for both meat, milk and wool [START_REF] Denis | Les Races d'animaux domestiques en France[END_REF]111], which could explain this greater morphological diversity, and the outnumber of sheep breeds compare to goat.

Archaeological sheep and goat dental remains can prove difficult to identify to the species level [START_REF] Zeder | Assessing the reliability of criteria used to identify postcranial bones in sheep, Ovis, and goats, Capra[END_REF][START_REF] Halstead | Sorting the sheep from the goats: Morphological distinctions between the mandibles and mandibular teeth of adult Ovis and Capra[END_REF][START_REF] Payne | Morphological distinctions between the mandibular teeth of young sheep, Ovis, and goats, Capra[END_REF][START_REF] Zeder | Assessing the reliability of criteria used to identify mandibles and mandibular teeth in sheep, Ovis, and goats, Capra[END_REF]112,113]. This study is a new contribution toward analyses dedicated to understanding the evolution of both sheep and goat in archaeological records using non-destructive morphometric identification of dental specimens. Being able to identify isolated third molars, that often remain unidentified in archaeozoology, allows new perspectives for further research. It is now possible to use the powerful and non-invasive quantitative tools of geometric morphometrics, via archaeophenomics [START_REF] Evin | Archaeophenomics of ancient domestic plants and animals using geometric morphometrics : a review[END_REF], to explore in more depth the separate evolution of sheep and goat through time, and follow their socio-economic and cultural roles in past animal husbandry. Such studies require secure identification and quantification which, particularly for large scale assemblages, is best accomplished through integrating molecular and non-molecular methodologies. Integrative approaches that combine GMM and biomolecular methods, such as ZooMS, are in high demand, with such integration leading to clearer perceptions of domestic species evolution through time and space.

Figure 3

 3 Figure 3 Third lower molar size variation: Boxplot showing variation in log centroid size between modern (dark) and archaeological (light) sheep (green) and goat (blue) identified by ZooMS. For differences between sites and chronologies, see SI figure 2.

Figure 4 :

 4 Figure 4: Two first axes of the PCA showing variation between modern (dark) and archaeological (light) sheep (green) and goat (blue) identified by ZooMS. PCA for each species are available in SI figure 3.

Figure 5 :

 5 Figure 5: Visualisation of the mean shape (A) and dissimilarity network (B) between archaeological (light) and modern (dark) sheep (green) and goat (blue). Archaeological specimens are identified by ZooMS.

Table 1

 1 Description of the 13 archaeological sites sampled, indicating geographic origin, chronocultural period of occupation, and number of specimens studied by geometric morphometrics and ZooMS. Location of the 13 archaeological sites of the sampled specimens. The numbers refer to the column ID of table 1.

	ID Country City	Name of the Site	Chrono-cultural Occupation	Number specimens	of	References
	1 Spain	Sant Quirze del Vallès	Bòbila Mas Duran Madurell-	Neolithic		1	[53]
	2 Spain	Maldà	Cantorella	Neolithic		1	[54]
	2 Spain	Maldà	Cantorella	Bronze Age		1	[54]
	3 Spain	Sabadell	Bòbila Madurell -Can Gambús 1	Bronze Age		2	[55]
	1 Spain	Sant Quirze del Vallès	Bòbila Mas Duran Madurell-	Iron Age		2	[53]
	4 France Sigean	Pech Maho	Iron Age		15	[56]
	5 Spain	Arbeca	Vilars	Iron Age		12	[57]
	6 Spain	Calafell	Alorda Park	Iron Age		1	[58]
	7 Spain	L'Escala	Empúries	Iron Age		1	[59]
	8 Spain	Sabadell	Can Feu	Iron Age		1	[60]
				Iron	Age/
	4 France Sigean	Pech Maho	Antiquity		1	[56]
				transition

  The 109 identified were based on 31 LDA with cross-validation percentages above 94.22% [90% Confidence Interval, CI: 93.01-95.50].For the same 109 specimens, ZooMS identified 64 (58.7%) specimens as sheep and 45 (41.3%) as goat (all spectra supporting the results are available in the Dryad Digital Repository https://doi.org/10.5061/dryad.fj6q5741m[START_REF] Jeanjean | ZooMS Spectra from: ZooMS confirms geometric morphometrics species identification of ancient sheep and goat[END_REF]).The two approaches provided the same species identification for 104 (95.4%) of the archaeological specimens corresponding to 59 sheep and 45 goats. The five (4.6%) remaining specimens are from the sites of Can Gambús (Catalonia, Late Antiquity) (n=3) and Missignac (France, Late Antiquity and Middle Age) (n=2). These misidentified specimens are all morphometrically identified as goat but molecularly identified as sheep. All of these misidentified specimens fell into the most common age category of 4-8 years old (i.e. 11G) (SI Table2). Looking more closely at these misidentifications, it appears that three out of the five specimens were identified 100% of the time to goat via GMM (i.e. all LDA give 'goat' as identification, no matter the probability of each of the 100 LDAs), while the two last ones were attributed respectively to goat only in 74.2% and 81.8% of the cases (See SI figure1for probability of each LDA).

Table 2 :

 2 Differences in size between modern sheep and goat, and the archaeological specimens identified by ZooMS. Pairwise Wilcoxon test were done between each group. The p-values in bold are significant (p-value < 0.05) after Bonferroni correction for multiple comparisons.

		Archaeological Goat Archaeological Sheep Modern Goat
	Archaeological Sheep 1
	Modern Goat	3.9x10 -6	1.4x10 -7
	Modern Sheep	9.6x10 -11	2.9x10 -14	9.1x10 -8

Table 3

 3 Comparison of size and shape variance between modern and archaeological sheep and goat. Size variance was tested with Fligner test, with a Khi 2 test statistic, and shape variance with a disparity approach (pairwise test on Procrustes variance, distance between each Procrustes variance).

	Size		Shape
	Khi²	p-value	distance	p-value

ACKNOWLEDGEMENTS

This project received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 852573). We would like to thank the Laboratori de Proteòmica CSIC-UAB at the Universitat Autònoma de Barcelona for access to the Bruker Ultraflex III MALDI-TOF-MS, and the Centres Cientifics i Tecnologics of the Universitat de Barcelona for access to the AB SCIEX 4700 MALDI-TOF-MS. We acknowledge the Service Regional de l'Archéologie of the Occitanie Region for permission of studying and destructive sampling of the archaeological material, as well as the directors of the different excavations included in this study.

We also acknowledge Museum National d'Histoire Naturelle (Paris), Museum für Haustierkunde, Oniris Nantes and the Musée d'Histoire Naturelle de Nîmes for allowing to consult their collections.

CREDIT

SI Figure 1: Identification probability of each LDAs for each archaeological specimen, sorted by archaeological sites numbers. The upper part of the graph, in green, corresponds to the probability of each LDAs for specimens identified as sheep by the geometric morphometrics, while the lower part, in blue, corresponds to the goats identified by the geometric morphometrics. The colour of each boxplot correspond to the ZooMS identification. Specimens with a red star are misidentified by geometric morphometrics. Specimens with a black star are those which are both identify as goat or sheep depending on the LDAs. The most numerous identification was kept. SI Figure 2: Boxplot showing differences of log centroid size for goat and sheep between archaeological sites and chrono-cultural occupation. Refers table 1


Esteban Restrepo
email: esteban.restrepo@inria.fr

Senior Member, IEEE Paolo Robuffo Giordano

Distributed Biconnecitvity Achievement and Preservation in Multi-agent Systems

Keywords: Multi-robot systems, distributed control, network connectivity

We propose a distributed control law to increase the robustness of a multi-agent network to node failure or removal. More precisely, our approach is able to maintain biconnectivity of an initially biconnected graph. Remarkably, if the graph is not initially biconnected or if the property is lost after a node removal or failure, our approach is also able to render the graph biconnected after a time instant which can be tuned by the user. The proposed control algorithm is completely distributed using only locally available information and requires the estimation of a single global parameter akin to existing connectivitymaintenance algorithms in the literature. Numerical simulations illustrate the effectiveness of our approach.

I. INTRODUCTION

I T is well known that connectivity of the graph representing the interconnections in a multi-agent system, that is, the existence of a path from any node to any other node in the network, is a necessary and sufficient condition for convergence to the agreement subspace in consensus-based coordination algorithms [START_REF] Mesbahi | Graph theoretic methods in multiagent networks[END_REF]. In light of this, a great number of algorithms have been proposed in order to maintain the connectivity of the graph in networks of autonomous multi-agent systems where the interaction among the nodes is limited to a finite region of the space around each agent-see e.g. [START_REF] Sabattini | Decentralized Estimation and Control for Preserving the Strong Connectivity of Directed Graphs[END_REF]- [START_REF] Ajorlou | A class of bounded distributed control strategies for connectivity preservation in multi-agent systems[END_REF].

Mobile agents in realistic scenarios, such as, e.g., robotic vehicles, are prone to multiple kinds of failures or may be required to leave the network in order to fulfill some specific task. These situations, however, may lead to disconnection of the graph and subsequently to the failure of the systemwise coordination mission. Therefore, it is necessary to design algorithms that not only maintain connectivity but that render the graph robust to node removals. In other words, the graph must remain connected even after one (or several) of the nodes and all its incident edges are removed. Such a property is the so-called biconnectivity [START_REF] Golumbic | Algorithmic graph theory and perfect graphs[END_REF]. In the present paper we propose a control law that is able to render the graph biconnected, possibly starting from a non biconnected graph, and to maintain said property for all time.

Biconnectivity maintenance has been much less studied than the simple connectivity counterpart. In [START_REF] Zareh | Enforcing biconnectivity in multi-robot systems[END_REF] a control law is proposed that enforces biconnectivity of the graph by using the gradient of the third smallest eigenvalue of the perturbed Laplacian. However, the third smallest eigenvalue is difficult to estimate in a distributed way and each agent would need to estimate a different eigenvalue for every critical node in the graph, increasing the amount of computations each agent needs to perform. Moreover, this approach is able to render a graph biconnected but no guarantees are given on the preservation of the biconnectivity in the presence of, e.g. disturbances or other control objectives. In [START_REF] Murayama | Distributed control for bi-connectivity of multi-robot network[END_REF] a biconnectivity-maintenance-andrestoration approach is proposed. However, it needs that each node i continuously estimates the values of the second smallest eigenvalue of the perturbed Laplacian at every node, that is, each node needs to estimate N different second smallest eigenvalues, making this approach not scalable. An extension of the latter is proposed in [START_REF]Dynamical improvement of network connectivity for multi-robot robustness using perturbed algebraic connectivity[END_REF] where the authors define a set of frangible nodes, i.e., critical nodes that may soon fail or be removed from the graph. A perturbed Laplacian is computed multiplying the edges of all the frangible nodes by a small constant. Then, the so-called perturbed algebraic connectivity is calculated. Based on a gradient control of the perturbed algebraic connectivity the robustness of the graph with respect to the removal of the frangible nodes is improved. However, it only increases the number of edges of the critical nodes and holds only if the initial condition on the perturbed algebraic connectivity is satisfied. That is, no guarantees are given on the maintenance or achievement of biconnectivity.

In this letter we propose a gradient-based biconnectivity controller that solves two of the open challenges not addressed by the previous works in the literature mentioned above, namely, scalability and guarantees for biconnectivity maintenance and recovery. Indeed, on one hand, in contrast to the works mentioned above, we guarantee both that the graph becomes and remains biconnected, even if it is not initially so or if biconnectivity is lost after, e.g., a node failure. Moreover, biconnectivity is achieved after a specified time instant, which can be tuned by the user. On the other hand, our approach relies only on locally-available information and only one estimation needs to be performed by the agents, that of the perturbed algebraic connectivity, akin to other global-connectivity-maintenance approaches in the literature, which greatly reduces the number of computations needed with respect to the approaches mentioned above, making our approach more scalable. This letter is organized as follows. In Section II the model and the problem statement are presented. The biconnectivity control is presented in Section III. Finally, the main results are illustrated via numerical simulations in Section IV and some concluding remarks are given in Section V.

II. MODEL AND GRAPH THEORY

We consider N agents modeled by the single-integrator dynamics ṗi = u i ,

where p i and u i denote, respectively, the position and control input of agent i. The agents communicate with each other and this communication is described by an undirected graph G(t) = (V, E(t)), where V = {1, 2, . . . , N } is the index set and E(t) ⊆ V 2 is the edge set characterizing the information exchange between agents, that is, an edge e k := (i, j) ∈ E(t), k = {1, . . . , M (t)}, is an ordered pair indicating that node j has access to information from node i.

We assume that the communication between the agents is proximity-based, that is, and edge e ij connecting robots i and j exists if and only if ||p i -p j || ≤ R, with R denoting the upper-bound on the distance under which the communication is reliable. Hence, the set of neighbors of agent i is given by

N i := {j ∈ V, | j ̸ = i and ||p i -p j || < R}. (2)
To model this proximity-based communication we use the adjacency matrix A ∈ R N ×N , where the element a ij , represents the weight of edge e ij . Namely, a ij > 0 if and only if j ∈ N i , and a ij = 0 otherwise. More precisely, given r < R and denoting d ij := ||p i -p j ||, we let

a ij =        1, if d ij ≤ r 1 2 1 + cos π(dij -r) R-r , if r < d ij ≤ R 0, R < d ij . (3)
Since G(t) is undirected, it follows that a ij = a ji . The Laplacian matrix L ∈ R N ×N is a symmetric positive semi-definite matrix given by L = diag{A1} -A, where 1 ∈ R N is the vector of all ones. It is well known that some fundamental properties of the graph are associated with the Laplacian matrix. Specifically, denoting λ 2 as the second smallest eigenvalue of L, we have that λ 2 > 0 if and only if G(t) is connected and λ 2 = 0 otherwise-see [START_REF] Godsil | Algebraic graph theory[END_REF] for more details on graph theory. Furthermore, we let ν 2 denote the unit-norm eigenvector of L associated with λ 2 . Under the weights definition (3), λ 2 becomes a connectivity measure of the proximity-based graph and, in particular, a function of the system state (e.g. of the agents relative positions).

Although maintaining the connectivity of the graph, i.e. λ 2 (t) > 0 for all t, is a sufficient condition for most coordination tasks in multi-agent systems, it may prove lacking in certain scenarios where individual agents disconnect from the group due to, e.g., failures, attacks, recharging or user commands. In such scenarios, the removal of one of the robots may lead the remaining graph to become disconnected.

Biconnectivity is the property of a graph to be robust to a node removal [START_REF] Golumbic | Algorithmic graph theory and perfect graphs[END_REF]. More precisely, let G -i be the graph remaining after the removal of node i. Then, G is said to be biconnected if, for any i ∈ V, G -i is connected. If G -i is disconnected, node i is called an articulation point. Therefore, an equivalent definition for a biconnected graph is for it to be a connected graph with no articulation points.

A. Perturbed algebraic connectivity

In order to consider the robustness of a graph with respect to a node removal we define the perturbed adjacency matrix A i (ε), whose elements are given by

a i kj (ε) = εa kj , if k = i or j = i a kj , otherwise, (4)
where ε is a small positive constant. The perturbed adjacency matrix A i (ε) captures the case in which the node i is almost removed from the graph. Indeed, note that ε → 0 would correspond to the case where node i is effectively removed. Now, the perturbed Laplacian matrix (at node i) is defined as

L i (ε) = diag{A i (ε)1} -A i (ε).
The second smallest eigenvalue of the perturbed Laplacian at node i is denoted λ i 2 (ε) and its associated eigenvector

ν i 2 (ε). If node i is an articulation point, then the graph G -i is composed of M i remaining connected components C i m , m ∈ {1, . . . , M i }.
In [START_REF] Murayama | Articulation node importance estimation and its correctness for robustness of multi-robot network[END_REF] it is shown that the connectivity of G -i and the size of C i m is related to the value

λi 2 = lim ε→0 + λ i 2 (ε) ε , (5)
which is called the perturbed algebraic connectivity. Indeed, letting ν i 2,i (ε) be the ith component of ν i 2 (ε) and denoting by C i m the number of nodes in C i m and by D i m the sum of their weights, i.e., D i m = j∈Ni C i m a ij , it is shown in [START_REF] Murayama | Articulation node importance estimation and its correctness for robustness of multi-robot network[END_REF] that

λi 2 = N j=1 a ij (ν i 2,i (0) -ν i 2,j (0)) 2 ||ν i 2 (0)|| , ( 6
)
λi 2 ν i 2,j = D i m C i m (ν i 2,j -ν i 2,i ), ∀j ∈ C i m . (7)
From ( 7) and the definition of ν i 2 , namely j ν i 2 = 0, if ν i 2,i ̸ = 0, we get the algebraic equation

M i m=1 D i m C i m D i m -C i m λi 2 + 1 = 0. (8)
Based on the algebraic equation ( 8), it is shown in [START_REF]Dynamical improvement of network connectivity for multi-robot robustness using perturbed algebraic connectivity[END_REF] that there exists a lower bound λ such that if the perturbed algebraic connectivity satisfies λi 2 > λ then the remaining graph G -i is connected.

III. BICONNECTIVITY CONTROL

A. Perturbed graph

In this paper, the objective is to design a control law that renders the graph biconnected and maintains said property for all time. For this purpose we make use of the perturbed Laplacian and the perturbed algebraic connectivity. However, instead of considering the perturbed graph only at a single critical node, we set a general perturbed graph that accounts simultaneously for all the articulation points in the graph.

In light of the above, we associate each agent with a dynamic parameter, which will be defined later, denoted ρ i ∈ [0, 1 -δ], with δ > 0 a small constant, and we let ε i = 1 -ρ i . Then, we define the perturbed adjacency matrix Ã, whose elements are given by

ãij (d ij ) = min{ε i , ε j }a ij (d ij ). (9)
Note that if either node i or j are critical points, then the edge-weights a ij and a ji are multiplied by a small constant to express that this node is almost disconnected. Now, the perturbed Laplacian matrix is given by L = diag{ Ã1}-Ã. Let λ2 be the second smallest eigenvalue of L and ν2 its associated eigenvector. Then, with a slight abuse of notation, we denote the perturbed algebraic connectivity of L by λ2 , which is approximated by

λ2 ≈ λ2 δ . ( 10
)
The main idea is to set ρ i = 1 -δ (equivalently ε i = δ) if node i is an articulation point, and to apply a control law based on the gradient of the perturbed algebraic connectivity to maintain λ2 above a threshold λ, defined in Section III-B below, in order to guarantee the biconnectivity of the graph. Then, when the node stops being an articulation point, ρ i smoothly decreases to 0, so that it is no longer a "critical node" and the normal edge weights are recovered with ε i = 1.

The authors in [START_REF] Zareh | Decentralized biconnectivity conditions in multi-robot systems[END_REF] proposed a way in which each node can verify if it is an articulation node based on the value of the third smallest eigenvalue λ i 3 of the perturbed Laplacian L i (ε) (at a single node i). This approach, however, requires the knowledge, by every robot, of the third smallest eigenvalue (and its associated eigenvector) of L i (ε) at each articulation point. However, the third smallest eigenvalue is difficult to estimate in a distributed way and, moreover, each agent would need to perform a different estimation of λ i 3 for every articulation node in the graph, increasing the amount of computations performed by each agent.

In contrast, in this paper we propose to check if a node is locally biconnected. Let G l i ⊂ G denote the local subgraph centered at node i and formed by the neighbors of node i, without itself. That is,

G l i = (V l i , E l i )
, where V l i = N i and an edge e kj ∈ E l i exits if and only if e kj ∈ E, with k, j ∈ N i . Then, a node is called locally biconnected if the second smallest eigenvalue λ l 2,i of the local graph G l i is positive. Remark 1: In an undirected graph, to characterize the local subgraph, each node only needs to receive the positions of its neighbors. Then, based on these, the local Laplacian matrices can be determined.

• A sufficient condition for a graph to be biconnected is that all its nodes are locally biconnected [START_REF] Zareh | Decentralized biconnectivity conditions in multi-robot systems[END_REF]. Therefore, we define a dynamic law for the parameter ρ i based on the connectivity of the local graph. More precisely we set

ρi = -κ 1 ρ i + κ 2 2 1 + sign(σ λ -λ l 2,i ) . (11)
The dynamic system ( 11) is an exponentially stable system with an additive disturbance which is different from zero only when λ l 2,i is smaller than a small threshold σ λ . This means that if a node i is not locally biconnected, i.e. λ l 2,i = 0, the value of ρ i increases and, equivalently, ϵ i decreases, thereby "perturbing" the graph. We set ρ i (t 0 ) = 0 if node i is locally biconnected at t 0 , and ρ i (t 0 ) = 1 -δ otherwise.

Remark 2: From (11) we have that ρ i (t) ∈ [0, 1 -δ] for all t ≥ t 0 . To see this, note that the second term on the right-hand side of ( 11) is always positive. Therefore, ρi (t) ≥ -κ 1 ρ i (t), together with ρ i (t 0 ) ∈ [0, 1 -δ], means that ρ i (t) ≥ 0 for all t. Then, the derivative of V ρ = 0.5ρ 2 i yields V ≤ 0 for κ 1 ρ i ≥ κ 2 , hence, the set R i := {ρ i ∈ R : 0 ≤ ρ i ≤ κ2 κ1 } is forward invariant. Choosing κ 1 and κ 2 such that κ2 κ1 = 1 -δ we get the result.

• Remark 3: Designing a control law based on the perturbation parameter ρ i given in [START_REF] Murayama | Articulation node importance estimation and its correctness for robustness of multi-robot network[END_REF] would render every node locally biconnected, which is more conservative than just achieving biconnectivity, since the graph may become biconnected even if every node is not locally so. However, using λ l 2,i instead of, e.g., λ i 3 , offers a good trade-off between conservativeness and complexity as it requires less computations and can be calculated locally.

•

B. Lower bound of the perturbed algebraic connectivity

In order to derive a lower bound for λ2 that guarantees that the graph remains biconnected we use the algebraic equation [START_REF] Murayama | Distributed control for bi-connectivity of multi-robot network[END_REF]. The following result, originally proposed in [8, Th. 3.1], is adapted hereafter to better fit the notations of this paper.

First let us assume that there is a single agent i such that ρ i > 0, that is, there is a single "critical" node. Assume further that the graph G -i is disconnected with two connected components. Let us define α, β ∈ (0, 1), then, without loss of generality, the number of nodes in each connected component and the sum of their weights are given, respectively, by

C i 1 = β(N -1), C i 2 = (1 -β)(N -1), and D i 1 = αD i , D i 2 = (1 -α)D i .
Then, from (8), we have

D i αβ(N -1) β(N -1) λ2 -D i α + D i (1 -α)(1 -β)(N -1) (1 -β)(N -1) λ2 -D i (1 -α) = 1, (12)
where

β ∈ 1 N -1 , 2 N -1 , . . . , N -2 N -1 (13)
since C i 1 = β(N -1) should be an integer. Hence, from (12) we have that the perturbed algebraic connectivity satisfies

αD i β(N -1) ≤ λ2 ≤ (1 -α)D i (1 -β)(N -1) . ( 14
)
The maximal algebraic connectivity is then found by computing the gradient of λ2 with respect to the parameters α and β (see [START_REF] Murayama | Distributed control for bi-connectivity of multi-robot network[END_REF]Th. 3.1] for more details), and is given by

λ2 = num(D i , N ) den(D i , N ) , ( 15
) num(D i , N ) = 2D 2 i N (N -1) -D i N (N -6) -9D i +D i (N -3) 4D 2 i N 2 -4D i N (N + 1) + N 2 -2N + 9, den(D i , N ) = 8D 2 i N (N -2) + 4D i N (N -6) + 36D i -2N 2 + 12N -18).
The value in (15) corresponds to the maximum algebraic connectivity if the node i is an articulation point. Moreover, since in the presence of multiple critical nodes the perturbed algebraic connectivity is smaller or equal than for a single agent, if we make λ2 > λ, where λ is set to the right-hand side of (15), then the graph remains biconnected.

Note that in this paper, unlike in [START_REF] Murayama | Distributed control for bi-connectivity of multi-robot network[END_REF], the expression in ( 15) is a function of the degree of node i, D i , since we do not consider the normalized Laplacian. However, from the definition of the weights (3), the value of the degree is time-varying and depends on the number of neighbors and robot's i proximity to them. In light of the latter, we simplify (15) by taking the maximum value for the degree, given by D i = N -1 when the node i is completely connected and sufficiently close to its neighbors, and we obtain

λ = num(N ) 8N 4 -28N 3 + 11N 2 + 50N -45 , ( 16
)
num(N ) := (N -1) 2N 3 -5N 2 + 8N -9 +(N -3) 4N 4 -12N 3 + 5N 2 + 2N + 9 .
The bound ( 16) is more conservative than (15) and, in principle, would drive the critical nodes to become completely connected. However, the latter does not occur since, as soon as a node i becomes locally biconnected, the last term on the right-hand side of ( 11) vanishes and ρ i converges to 0 exponentially, and the node stops being a critical node. Moreover, as can be seen in Fig. 1, λ as a function of the number of agents N ranges from 0.47 to 0.5 with an asymptote at 0.5. Therefore, without loss of generality, one could take λ = 0.5. This makes our approach, compared to the works in the literature, more scalable since it does not depend on the number of agents in the system and is useful in situations where the number of agents changes or is unknown.

C. Control law for biconnectivity

The design of the control law is inspired by the connectivity preserving controller proposed in [START_REF] Sabattini | Distributed Control of Multirobot Systems With Global Connectivity Maintenance[END_REF]. Thus, we introduce the scalar function V λ : D → R ≥0 , λ2 → V λ ( λ2 ), which is C 1 over its domain and has the property that V λ ( λ2 ) → ∞ as λ2 → ∂D. Then, the gradient controller is set to be

u c i = -κ ∂V λ ( λ2 ) ∂p i = -κ ∂V λ ∂ λ2 λ2 ∂p i , κ > 0. (17)
From [START_REF] Sabattini | Distributed Control of Multirobot Systems With Global Connectivity Maintenance[END_REF] and recalling [START_REF] Godsil | Algebraic graph theory[END_REF], we have that

∂ λ2 ∂p i = 1 δ ν⊤ 2 ∂ L ∂p i ν2 = 1 δ N j=1 (ν i -νj ) 2 ∂ã ij ∂p i . (18)
Replacing ( 18) into (17) we obtain

u c i = -κ ∂V λ ∂ λ2 N j=1 (ν i -νj ) 2 δ ∂ã ij ∂p i . (19)
Moreover choosing

V λ ( λ2 ) = coth( λ2 -λ * ) -1, (20)
where λ * is a properly chosen bound, the gradient control law becomes

u c i = -κ csch 2 ( λ2 -λ * ) N j=1 (ν i -νj ) 2 δ ∂ã ij ∂p i . (21)
The objective of biconnectivity maintenance is equivalent to rendering forward invariant the set

D := { λ2 ∈ R ≥0 | λ2 > λ * }. (22)
It can be shown that under the gradient control (17), with λ * = λ in ( 16), the multi-agent system remains biconnected for all time. However, biconnectivity may be lost after a node is removed due to, e.g., a failure. Hence, the initial graph may not be biconnected, i.e., there is at least one articulation point, and it may be possible that λ2 ≤ λ = λ * . But, since V λ is only defined in the set (22), the controller (17) would not be able to render the graph biconnected, indeed it may not be well defined. In order to overcome this problem, we set instead λ * as the solution to the dynamic system

λ * (t) = -κ λ (λ * (t) -( λ + ϵ)), (23)
where ϵ > 0 is a small constant, λ is given by ( 16), and λ * (t 0 ) = 0. Note that, under (23) and λ * (t 0 ) = 0, λ * (t) is positive for all t and exponentially converges to λ + ϵ with a rate determined by κ λ > 0. Since we assume that the initial graph is at least connected and ε i ≥ δ, for all i, we have that λ2 (t 0 ) > 0 and therefore, given λ * (t 0 ) = 0, we have that λ2 (t 0 ) ∈ D so that V λ and the input (21) are well defined. Under the control law (21), setting λ * (t) as the solution to (23), we guarantee that there exists a (tunable) time instant T (κ λ ) > t 0 such that for all t ≥ T (κ λ ) the graph is biconnected. This is stated in the following result.

Proposition 1 (Main result): Assume the initial network G(t 0 ) is connected. Then, under the gradient control law (17) with λ(t) * as the solution to (23) with λ(t 0 ) * = 0 and λ as in ( 16), the graph G(t) becomes and remains biconnected. More precisely, the ( 22) is rendered forward invariant and there exists a (tunable) time instant T (κ λ ) > t 0 such that G(t) is biconnected for all t ≥ T (κ λ ).

□ Proof: Assume that λ2 (t) ∈ D (this assumption will be relaxed later). Define the potential function V λ : R ≥0 × D → R >0 , given by V λ (t, λ2 ) = coth( λ2 -λ * (t)) -1, which is positive, i.e., V λ (t, λ2 ) > 0, in its domain of definition. Its derivative yields

Vλ (t, λ2 ) = ∂V λ ∂ λ2 λ2 ∂p ⊤ ṗ + ∂V λ ∂t =csch 2 ( λ2 -λ * (t))   -κ λ2 ∂p 2 -λ * (t)   . ( 24
)
Given ( 23) and λ * (t 0 ), we have that λ * (t) ≥ 0 for all t ≥ t 0 . Therefore, we have

Vλ (t, λ2 ) ≤ -κ csch 2 ( λ2 -λ * (t)) λ2 ∂p 2 ≤ 0. (25)
From (25) we conclude that V λ is bounded along the trajectories of p(t). Now, in order to prove the forward invariance of D we proceed by contradiction. Let us assume that there exists a τ such that λ2 (t) ∈ D for all t ∈ [t 0 , τ ) and λ2 (τ ) / ∈ D. Therefore, from continuity of the solutions we have that λ2 (t) → ∂D as t → τ . From the definition of V λ this implies that V λ (t, λ2 (t)) → ∞. However, the latter is in contradiction with (25) which implies that V λ (t, λ2 (t)) is bounded. Therefore, we conclude that D is forward invariant.

Finally, since from (23) λ * exponentially converges to λ + ϵ with a rate determined by κ λ , there exists a (tunable) time instant T (κ λ ) > t 0 such that λ * (T ) = λ. From the forward invariance of D, λ2 (t) > λ * (t) holds for all t ≥ t 0 . Therefore, λ2 (t) > λ for all t ≥ T (κ λ ) implying that the graph becomes and remains biconnected for all t ≥ T (κ λ ).

IV. NUMERICAL EXAMPLES

In this section we present some numerical examples to validate our theoretical results. The multi-agent system consists of N = 21 agents that interact over a proximity-based undirected graph where the edge weights are given by (3) with R = 2.2 m and r = 1.2 m. The initial configuration with the initial graph are represented in Fig. 2. From the latter, it is clear that the initial graph is not biconnected, therefore λ2 (t 0 ) < λ. However, since λ * (t) is the solution of (23), we have that λ2 (t 0 ) > λ * (t 0 )-see Fig. 4. 4 that the value of the perturbed algebraic connectivity λ2 remains above the imposed time-varying lower bound λ(t) * . The latter guarantees that the graph becomes biconnected as can be seen in Figs. 3 and5. Indeed note that the parameters ρ i become 0 after 3s.

-5 -3 -1 1 3 5 7 -6 -4 -2 0 2 4 x [m] y [m]
For comparison Fig. 6 shows the results when the agents' inputs are given by the connectivity maintenance control in [2]. Although the connectivity is maintained, it is clear from Fig. 6 that the graph is not biconnected.

Finally, to demonstrate the applicability of the proposed approach in a more meaningful example we studied the scenario where agents may join and leave the network. In this scenario, the agents start from the same initial configuration as before (Fig. 2). As the agents evolve under the biconnectivity law described above, three new agents (initial positions marked with a triangle in Fig. 7) join the network at around 6s, 10s, and 17s. Moreover, at 20s and 24s, respectively, two agents (initial positions marked with a square in Fig. 7), leave the network. Figure 8 depicts the evolution of the number of agents over time. Whenever an agent leaves or joins, the value of λ * is reset to zero in order to account for the possible loss of biconnectivity. This could be distributedly achieved in practice using a flooding algorithm [START_REF] Oliva | Distributed finite-time average-consensus with limited computational and storage capability[END_REF] to update the value of λ * . Then, biconnectivity is recovered under our proposed controller. This can be seen in Fig. 9 where the perturbed algebraic connectivity remains at all times above the prescribed bound. Figure 10 reports the values of the parameters ρ i , where it can be seen that the local biconnectivity is recovered after each addition or removal.

V. CONCLUSIONS

In realistic multi-agent scenarios with limited-range communication, it is necessary to account for potential disconnections of the network due to the removal of one or multiple nodes following, e.g., a failure or an attack. In such cases, the more robust property of biconnectivity is preferred to guarantee robustness with respect to node removal. We propose a gradient-based control law that is able, using only locally available information, to not only maintain biconnectivity but to also recover it, after it is lost, in tunable pre-specified time. Current and future work focuses, on one hand, on a quantitative and experimental comparison between the localbiconnectivity-based approach proposed in this letter and the approaches based on computing λ 3 , and on the other hand, considering field-of-view and collision-avoidance constraints in "generalized biconnectivity" framework, and applying the results to the context of open multi-robot systems.

Fig. 1 :

 1 Fig. 1: Variation of λ with respect to N .

Fig. 2 :

 2 Fig. 2: Initial configuration and initial graph.The parameters of the dynamic law[START_REF] Murayama | Articulation node importance estimation and its correctness for robustness of multi-robot network[END_REF] are set to δ = 0.02, κ 1 = 10, κ 2 = 9.8, and σ λ = 0.005, and the gain of the controller (21) is set to κ = 0.25. Moreover we set κ λ = 0.3 in (23) and λ ≈ 0.48, obtained from (16) by replacing N = 21. Figs.3-5show the results of applying the controller (23). It can be seen from Fig.4that the value of the perturbed algebraic connectivity λ2 remains above the imposed time-varying lower bound λ(t) * . The latter guarantees that the graph becomes biconnected as can be seen in Figs.3 and 5. Indeed note that the parameters ρ i become 0 after 3s.For comparison Fig.6shows the results when the agents' inputs are given by the connectivity maintenance control in

Fig. 3 :

 3 Fig. 3: Trajectories of the agents and final graph.

Fig. 4 :

 4 Fig. 4: Perturbed algebraic connectivity. The dashed red line represents the bound λ * (t).

Fig. 5 :

 5 Fig. 5: Dynamic indicator of local biconnectivity. ρ i = 0 means that node i is locally biconnected.

Fig. 6 :

 6 Fig. 6: Trajectories of the and final graph.

Fig. 7 :

 7 Fig. 7: Trajectories of the agents final graph.

Fig. 8 :

 8 Fig. 8: Evolution of the number of agents in the network.

Fig. 9 :

 9 Fig. 9: Perturbed algebraic connectivity. The dashed red line represents the bound λ * (t).

Fig. 10 :

 10 Fig. 10: Dynamic indicator of local biconnectivity ρ i .

This work was supported by the ANR-20-CHIA-0017 project "MULTI-SHARED"