Topological Data Analysis and Multiple Kernel Learning for Species Identification of Modern and Archaeological Small Ruminants

Research aim

The main contribution of this paper is twofold. Firstly, a machine learning routine is proposed to automatically identify morphologically related animal species, like small ruminant herbivores such as sheep, goat, Alpine ibex, roe deer and gazelle. This routine relies on topological data analysis. Secondly, an original statistical framework is presented, allowing to weight the extracted topological features in such a way as to exploit each one of them (i.e. multiple kernel learning).

Three additional remarks are needed. First, although the classification is automatic, the expertise of the zooarchaeologist is required for the analysis of the results to construct an analytical framework and gain an understanding of the way statistics work. Second, the above-mentioned species have been chosen for three main reasons: 1) their morphological proximity; 2) their simultaneous presence in certain archaeological contexts; 3) the large number of taxonomic criteria available in the literature to differentiate some of them, such as sheep and goats, compared to the lack of data for others, such as roe deer and gazelle. Third, the model described here is based on whole bones in a good state of preservation.

2.1. Astragalus

The anatomical part selected for this study is a short bone, the astragalus, from the tarsal joint of the foot (Supplementary data 1). It preserves very well in archaeological faunas because it is a small bone, particularly compact and robust and rarely broken intentionally due to its low nutritional value (Barone 1976; Popkin et al. 2012). These bones presented several anatomical criteria for identifying wild and domestic ungulates discussed by multiple scholars for almost sixty years. The distinction between sheep and goats is well documented (Boessneck et al. 1964; Clutton-Brock et al. 1990; Fernandez 2001; Prummel & Frish 1986; Salvagno & Albarella 2017; Zeder & Lapham 2010; Zeder & Pilaar 2010), but still poses problems (Sipilä et al. 2023) and represent a challenge as demonstrated by recent GMM (Gaastra, 2023; Jeanjean et al. 2022; Lloveras et al. 2022; Pöllath et al. 2019; Pöllath et al. 2018; Vuillien 2020) and molecular studies (Jeanjean et al. 2023; Le Meillour et al. 2020). In addition, there are few criteria for distinguishing between roe deer, gazelle and ibex (Buitenhuis 1988; Crégut-Bonnoure, 2020; Fernandez 2001; Gudea & Stan 2012; Lavocat 1966; Peters 1989).

2.2. 3D models dataset

The dataset included 150 3D complete astragali belonging to five taxa: Alpine ibex (Capra ibex), sheep (Ovis aries), goat (Capra hircus), roe deer (Capreolus capreolus) and gazelle (Gazella cuvieri, Gazella dorcas, Gazella spekei and Gazella sp.) (Supplementary data 2). This dataset does not consider the specimens’ geographical origin or provenance which do not concern our research topic. Gazelle species are also grouped at the genus level for statistical reasons. The specimens belong to National Museum of Natural History Mammalian and Birds collection of Paris, CEPAM (UMR 7264) and Archéorient (UMR 5133) labs zooarchaeological reference collection, modern sheep and goat collected for the EvoSheep collection (ANR-17-CE27-0004), modern goats form BALUT Laboratory Iran, archaeological site of “Grotte de l’Observatoire” Museum of Prehistoric Anthropology of Monaco and archaeological site of “Tell Sheikh Hassan” Archeorient lab (UMR 5133) (Table 1). In order to provide a homogeneous group for ML analysis, each taxon is represented by 30 specimens (3D astragali).

The astragali were scanned using the Artec Spider blue LED surface scanner and Artec Studio reconstruction software (version 16) and EinScan Pro 2X (Figure 1). The 3D models are reconstructed at a resolution between 0.3 and 0.1 mm using a textured polygonal mesh. Meshes are exported in “ASCII.ply” and “obj” archiving format (Vergnieux et al. 2017).

Figure 1

2.3. Topological data analysis (TDA)

Although an in-depth presentation of TDA clearly is outside of the scope of this paper (the interested reader is referred to Chazal & Michel 2021), in this section are sketched the main ideas TDA relies on.

The input data here is a collection of 150 3D point clouds, and the aim is to extract some useful information (or features) from each point cloud in order to use it to assign the cloud/bone to a species. In order to extract the features, a “continuous” shape is built from the point cloud by progressively connecting data points that are closer to each other with an edge. Here, the reader can safely consider that two points are close if their Euclidean distance in the 3D space is smaller than a given threshold ϵ > 0. As far as ϵ grows, more and more points are connected and the shape that appears is a collection of simplicial complexes (Figure 2). A simplicial complex can be seen as a higher-dimensional generalization of a neighbouring graph: whereas the latter only includes vertices and edges, the former also contains faces, namely triangles and tetrahedrons. The nested family of simplicial complexes that add to each other is called filtration and as long as the family grows, relevant topological features such as connected components, loops and voids are collected via specific methods, such as persistent homology (PH, Otter et al. 2017; Zomorodian & Carlsson 2004) and stored into the so-called persistence diagrams (PDs), that is described in some detail in the next section. It has been shown (Chazal & Michel 2021) that the topological features provide insights into the underlying structure of the data, making TDA particularly useful for analysing highly dimensional and noisy datasets.

Figure 2

In order to provide the reader with an intuition of what TDA is, the exposition of this pipeline is simplified. However, some remarks are needed:

1. Since the way simplicial complexes are built mainly relies on the relative distance between the points in the cloud, the orientation of the 3D shapes is irrelevant and there is no longer a need to register the collection.

2. Several notions of distance between points can be chosen and several ways of aggregating simplicial complexes (i.e. filtrations) exist.

3. The topological features are collected in PDs, but alternatives exist (e.g. persistence images, barcodes, etc.).

These few remarks should help the reader to figure out the vastness and richness of topological data analysis.

2.4. Discrete optimal transport

To quantify similarities between PDs, several metrics have been developed (Biasotti et al. 2011; Efrat et al. 2001). It should be stressed that the notions of similarity and distance between PDs are two sides of the same coin. Indeed, in general given two data points x and y if the distance d(x, y) between them is known, a measure of similarity is obtained via

for any real positive λ. That said, the study of distances between discrete probability distributions, using the Wasserstein distance (Lacombe et al., 2018) from optimal transport has been resorted (OT, Cuturi 2013). Optimal transport provides an intuitive way to quantify the similarity between two probability distributions by considering the minimum cost required to transform one distribution into the other. Here, discrete probability distributions are given consideration. In more details, two persistent diagrams $P_X\hspace{0.17em}=\hspace{0.17em}\{x_1\cdots ,x_N\}$ and $P_Y\hspace{0.17em}=\hspace{0.17em}\{y_1,\cdots ,y_M\},$ with N and M denoting the number of points in each diagram are considered. One can easily define a probability distribution μ_X (respectively μ_Y) by putting mass 1/N (1/M) over each point in P_X (P_Y).

Definition 1:

The Wasserstein distance of order p between μ_X and μ_Y is given by

2

$W_p\left({\mu }_X,{\mu }_Y\right)\hspace{0.17em}=\hspace{0.17em}{\left(\underset{\pi \in \Gamma \left({\mu }_X,{\mu }_Y\right)}{min}{\displaystyle \sum _{i,j}}d\left(x_i,y_j\right){\pi }_{ij}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$p$}\right.}$

where Γ(μ_X, μ_Y) denotes the set of joint probability mass functions (also called couplings or transport plans) on P_X × P_Y with marginals μ_X and μ_Y, respectively, and d(x,y) is the distance between points x and y in the underlying metric space.

Here, d(x, y) are considered to be the Euclidean distance and set p = 2. Intuitively, one has to imagine a total mass of 1 is split into N equal portions and distributed over the points of P_X. Now the aim is to carry all this mass from P_X to P_Y under the constraint that, at the end of the day, the total mass is uniformly distributed over the points of P_Y. The Wasserstein distance quantifies the optimal cost of such an operation, and the optimal transport plan (the one minimising Eq. (2)) tells how much weight one should move from where to where in order to minimise the effort.

However, the PDs gives an additional issue: a green point (for instance) on P_X could be partially or totally transported to a red or a blue point on P_Y, whereas the aim is to keep different topological features (connected components, loop and voids) well separated. Thus, a PD is “splitted” into three PDs, one for each colour (see Figure 4) and compute three Wasserstein distances for each pair of bones, one for each topological feature represented in the PDs.

Examples of pairwise kernel similarity matrices can be seen in Figures 3 and 5. Kernels were obtained from the Wasserstein distances via Eq. (1), where λ was set equal to the standard deviation of the corresponding distance matrix. The rows and columns of the kernel matrices correspond to the 150 bones and the i-th row and j-th column elements denote the Wasserstein similarity between the PDs associated with the i-th and j-th bones, for the corresponding topological features. The main diagonal is the brightest region of each kernel matrix since each bone is trivially at similarity one (and zero distance) from itself. Darker regions in a matrix correspond to higher distances. To narrow down the focus, the 1 and 2-dimensional features (respectively blue and green points in Figure 3) are only considered, which proved to be more informative.

Figure 3

2.5. Supervised multiple kernel learning

The blueprint of supervised machine learning can be briefly described as follows. Assume a training dataset of $\left\{x_1\cdots ,x_N\right\}$ observations is given, together with labels $\left\{y_1\cdots ,y_N\right\}$ Here, the i-th observation x_i is a 3D scan of a bone, in the form of a point cloud, and its label y_i can be seen as an integer ranging from 1 to Q and labelling the species of the bone. As seen in Section 2.2, N = 150 and Q = 5 for us. Next ingredient is a function of the data (the classifier), say fθ, depending on some parameters θ and associating to each observation xⁱ a predicted label $f_{\theta }\left(x_i\right)\hspace{0.17em}=\hspace{0.17em}{\widehat{y}}_i$, i.e. another integer between 1 and Q. then fθ is “trained” by solving the following minimisation problem

where L(·, ·) is a loss function. So, roughly speaking, the value of θ should be such that the mismatch between the predicted and the actual labels is minimal, on average. Once θ is optimised, the final aim is to be able to correctly predict the label y^* of a new test data point x^*, via fθ(x^*). In order to check that it is actually the case, in these experiments, all dataset (N = 150) has been split into train (Ntrain = 120) and test (Ntest = 30) and the accuracy (i.e. the proportion of correctly identified specimens) reported on the test dataset. To assess the robustness of the test accuracy, 50 random train/test splits are performed allowing the computation of a mean accuracy and a standard deviation.

Now, one way to define fθ is to pass through kernel matrices. Several ML classifiers are specifically designed to leverage the representation of the data through kernel matrices. Popular methods include support vector machines (SVM), kernel discriminant analysis (KDA) and kernel logistic regression (KLR). An in-depth description of such methods is outside the scope of this paper and the interested reader is referred to (Hastie et al. 2009, Chapters 4, 5 and 12).

In this work, the focus is on multi-class KLR with multiple kernels. Indeed, as mentioned above, different types of filtrations and/or topological features lead to different similarity matrices, each of which could contain different discriminant information and relying on a single kernel matrix chosen by the user could be not the optimal strategy. Multiple kernel learning (MKL) approaches [47] are specifically designed to manage this kind of situation: they allow one to properly weigh the input kernels and possibly discard the useless ones. Formally, D different kernel matrices $\left\{K^{\left(1\right)}, \dots , K^{\left(D\right)}\right\}$ with $K^{\left(d\right)}\in {ℝ}^{N\times N}$ were considered for each d. The aim is for an optimal convex combination of such matrices, defined by

By adopting a hybrid Bayesian formulation, the kernel weights β_d are treated as random variables, following an exponential prior distribution, whereas the other weights intervening in the KLR are treated as parameters to optimise. This approach allows to revisit MKL in an original way, being a com- promise between pure optimization strategies (Rakotomamonjy et al. 2008) and fully Bayesian ones (Damoulas & Girolami 2008; Gonen 2012). Moreover, recent developments in importance-weighted stochastic variational inference (Sobolev & Vetrov 2019) have enabled the optimization problem in a fully differentiable manner, via Stochastic Gradient Descent (SGD, Bottou 2010), and perform posterior inference on β₁, …, β_D.

3.1 Distance matrices (TDA) and optimal transport distances

The topological dimension 1 (loops) without normalisation indicates that Alpine ibex, roe deer, and gazelle respectively form dense clusters, in contrast to sheep and goat (Figure 3 Up). Sheep and goat exhibit an intra-specific significant topological heterogeneity, revealed by a marked colour gradient from dark blue to yellow. This is mainly due to both the breed factor and the geographical heterogeneity of the specimens studied. For sheep, specimens from Ethiopian breeds (in dark green and yellow) are distinct from specimens from French breeds (in dark blue) (Figure 4). The same is true for goat. Specimens attributed to breeds of Iranian origin (in dark blue) differ from specimens of French breeds (in green and yellow) except for one specimen (number 57) (Supplementary data 4). Furthermore, although ibex form a homogeneous cluster, one individual stands out (dark blue): this is the only modern ibex in the dataset. When comparing species, the same topological dimension reveals a high degree of separation between wild and domestic species. In more details: ibex form a single community; roe dee look quite similar to gazelle (two wild taxa) and sheep often look indistinguishable from goat (two domestic species) despite no apparent link to their breed or geographical origin. Nevertheless, there is a certain topological proximity between roe deer and gazelle, on one side, and specimens of French sheep and goat, on the other (whereas the morphology of gazelle and roe deer is clearly distinguishable from that of African sheep and Asian goat).

Figure 4

Also, the topological dimension 2 (holes) without normalisation highlights differences between Alpine ibexes and other species (Figure 3 down). These differences could be explained by the size of the astragalus. Ibexes are taller than the other species in the dataset. However, if the size effect were crucial, gazelle, the smallest taxon in the dataset, should be clearly distinguishable from the other taxa. This is not the case. Consequently, although this parameter is undoubtedly important, it does not appear to be the sole factor responsible for the structuring of the dataset. Interestingly, whereas in the previous matrix roe deer and gazelle looked quite similar, here the difference between them is more accentuated (although an increased similarity with goat appears). This point clearly illustrates why the adoption of several kernel matrices is beneficial to the taxonomic identification at the species level. The topological dimension 1 (loops) with bone normalisation allows for the clear distinction of Alpine ibex, roe deer and gazelle from groups of domestic caprine (sheep and goat) (Figure 5 Up). This outcome indicates that normalisation does not directly impact classification at the inter-specific level.

Figure 5

Finally, the topological dimension 2 (holes) with bone normalisation might seem of no particular interest (Figure 5 Down) since the similarities between specimens of each species render it impossible to distinguish between the species in question (except for roe deer, partially). However, as demonstrated in the next section, supervised MKL remarkably benefits from this matrix since it is the only one highlighting similarities between specimens of the same species far from each other in previous representations, especially sheep and goat.

3.2 Classification: Supervised MKL and automatic taxonomic identification of bones

For the supervised MKL part, as mentioned, the whole dataset (N = 150 bones) was randomly divided into 120 (80%) as train dataset and 30 (20%) as test dataset. It is recalled that such a random split is repeated 50 times. Table 2 reports the average test accuracy together with its standard deviation (in small characters). The second column reports the global test accuracy, whereas the remaining columns show the test accuracies for each species.

In order to better inspect how the test accuracy behaves as a function of the train/test data split, a box-and-whisker plot is provided in Figure 6a. From the bottom to the top, each box-and-whisker reports for each species the minimum test accuracy, the lower quartile, the median accuracy, the upper quartile, and the maximum test accuracy. Outliers are represented as single points. The only modern specimen present in the Alpine ibex dataset is the outlier. The gazelle outlier corresponds to two specimens: a modern specimen (specimen number 1992 1844) and an archaeological specimen from Tell Sheikh Hassan (specimen number TSH 7Z23).

Figure 6a

The average test accuracy of 81,1% (sensibly higher than 15% that one would obtain from a random assignment) is boosted by the accuracy on the wild species, which our approach can successfully identify. Figure 6a shows that the median test accuracy on wild species is 100%. The mean/median accuracies for sheep and goat are sensibly lower and come with a higher variability. However, as pointed out when describing the kernel matrices in Section 3.1, the non-homogeneous intra-specific blocks for these two species correspond to different breeds.

Table 3 shows the weights β₁, …, β₄, introduced in Eq. (2) and estimated via importance-weighted stochastic variational inference (IW-SVI). In more detail, our Bayesian inference strategy allows us to sample from the approximate posterior distribution of β₁, …, β₄, given the data and the other model parameters. Table 3 reports the posterior mean of each weight with its standard deviation (in small character). Kernel density estimates of weights’ posterior distributions can be seen in Figure 6b. The mode of each estimated density (i.e. the maximum a posteriori estimate of the corresponding weight) is away from zero, meaning that multiple KLR exploits all the topological features considered, either for normalized or unnormalized point clouds. Interestingly, the most important weight is put on β₄, corresponding to the apparently less expressive kernel matrix of Figure 6b. This happens precisely because this kernel matrix is the only one allowing the algorithm to densify the clusters of sheep and goat, respectively. Furthermore, the weight given to each matrix indicates the influence of size on classification. If size was the more discriminating factor, then β₁ (non-standardised matrix) should have a greater weight than β₃ (standardised matrix). The results show the opposite (Table 3), indicating that bone size is not the more discriminating factor for the final classification.

Figure 6b

This can clearly be seen in Figure 7 which shows the learned kernel matrix K introduced in Eq. (3) and obtained as a linear combination of the Wasserstein kernels, weighted via the optimal betas shown in Table 3. K exhibits brighter diagonal blocks (especially for wild species) and the blocks corresponding to sheep and goats look more similar to dense blocks with respect to what Figures 3 and 5 show.

Figure 7

4.1. Inter-specific identification: wild vs. domestic

About the main research question of the inter-specific identification of wild and domestic small ruminants, TDA proved to be very proficient at correctly identifying wild specimens from a 3D scan of their astragalus, whereas it suffers with the identification of domestic specimens.

In other words, TDA makes it possible to clearly distinguish ibex, roe deer, and gazelle astragali from sheep and goat astragali, which underlines its value. Indeed, it is often difficult to distinguish these wild animals from sheep and goats in archaeological contexts where all these species are represented and where anatomical criteria are not discriminating. In contrast, the use of TDA does not address the challenges faced by the zooarchaeological community in distinguishing sheep from goats. Here, the result is dependent on the discriminative capacity of TDA and the dataset; the obtained classification exhibits the intraspecific morphological variability of the selected sheep and goat breeds.

Nevertheless, the fact that the median accuracy is higher than the mean for the wild species and the presence of outliers in Alpine ibex and gazelle is explained by the heterogeneity of our dataset: a high intra-specific variability of a few specimens lowers the mean accuracy. For instance, the single modern Alpine ibex in the collection turns out to be dissimilar from the other ibexes. As we saw in the figures in Section 3, it is dissimilar from the other ibexes; with regard to the train dataset, multiple KLR can safely identify all the archaeological ibex in the test. Conversely, when the modern specimen is included in the test dataset, the classifier fails to identify it (since trained on archaeological ibex dated to the Upper Pleistocene), and turns it into an outlier in the Alpine ibex group. This outcome is consistent with the difficulties encountered by zooarchaeologists in using modern datasets to identify ancient animal populations. This is due to the contrast between the number of current species and subspecies compared to fossil species and their morphological diversity (Crégut-Bonnoure 2020; Crégut-Bonnoure & Fernandez 2018; Urena et al. 2018). However, our findings indicate that this approach should be employed with great caution when attempting to address the question of the evolution of morphology. Indeed, upon noting the presence of two outliers within the gazelle group, a more detailed examination reveals that interpreting these differences in relation to the available temporal origin information is challenging. The dissimilarities between modern and archaeological gazelles (Holocene, between 9000 and 6000 BC) are not immediately apparent: TDA makes it difficult to identify variability among the gazelle species. This is likely due to the dataset itself, which does not accurately reflect the morphological variability among the modern species (Gazella cuvieri, Gazella dorcas, Gazella spekei). However, we intend to compare available proteomic data (Culley et al. 2021; Janzen et al. 2021; Le Meillour et al. 2020; Le Meillour et al. 2023) and future morphological and morphometric studies (Vuillien 2024).

4.2. Beyond the classification problem: is astragalus a good taxonomic marker?

The results obtained for species classification prompt the question of whether the astragalus can be employed as an interspecific identification marker. The geographical provenance of domestic and wild specimens appears to exert a more pronounced impact than the distinction between species. As we demonstrated, TDA cannot see sheep and goats as uniform and separated clusters (within the limits of the explored filtrations) and this fact has some significant consequences for the weighting of the kernels. The misidentification of sheep and goat species is attributed to the morphological variability of domestic breeds and their geographical origin. In addition, the difference in the ibex group observed between the single modern specimen and the other archaeological specimens may be correlated with the morphological evolution of the Alpine ibex over time and would point to the potential of the astragalus as an ecomorphological marker (Barr 2014; DeGusta and Vrba, 2003; Plummer et al. 2008). The biometrical and biomolecular studies carried out for this species show its morphotypic diversity over time, linked to environmental changes, the rocky areas frequented and human pressure (Crégut-Bonnoure 2020). However, it is still difficult to identify these factors precisely from faunal remains. TDA could be a valuable tool for exploring morphological variability at the intra-specific level for domestic and wild species and their adaptation across time and the environment, such as recent GMM studies applied to the same modern sheep breeds (Bader et al. 2022) and archaeological Alpine ibex populations. In archaeozoology, identifying and classifying of domestic sheep and goat morphotypes is of great importance. Such an analysis can provide insight into the evolution of zootechnical practices, economic development, and human society (Vila et al. 2021).

4.3. More robust assessment of the result obtained via TDA

To correctly estimate the potential of the method proposed in this paper, it will be crucial in the future to analyse the same dataset with other approaches either directly based on human expertise to have a human accuracy relying on anatomical criteria or automatic, such as GMMs. This paper aimed to test TDA and supervised MKL on a zooarchaeological dataset, nevertheless, TDA features would certainly express the best of their potential in conjunction with other features, such as anatomical criteria and/or GMM patterns.

Another crucial aspect of machine learning routines in general is explainability: to know whether a classifier can be trusted, we need to understand how it works (on this topic, see Rudin 2019). Unfortunately, whereas TDA keeps track of the lifelong of topological features such as connected components, loops, and holes, it is currently not possible to know where these features are located on the surface of a 3D bone. This makes it very difficult, for instance, to assess whether the similarities that our approach detected between French sheep and roe deer/gazelle are based on meaningful and previously unexplored morphological patterns or if they are due to some other geometric components. Although some model-agnostic explaining techniques exist (Lundberg & Lee 2017; Ribeiro et al. 2016), their use in the context of 3D point clouds, in conjunction with TDA, is not immediate at all. This issue will be addressed in future research to propose new anatomical features. Finally, it will be of interest to test deep learning models. Although requiring a huge amount of training data, such models can obtain impressive results, and their behaviours could be more easily captured than the one of TDA, either via explaining techniques or ad hoc architectures.

Conclusion and perspectives

This paper mainly focuses on the taxonomic identification of wild and domestic small ruminants from 3D scans of complete astragalus of modern and archaeological specimens. The problem was framed as a supervised learning problem and addressed using TDA, optimal transport and an original inference routine. From one side, the topological features extracted with TDA proved to be very discriminant in classifying wild species (Alpine ibex, roe deer, gazelle). The main strengths of the proposed approach are a median test accuracy of 100% for these species and the fact that our routine is entirely automated (the expert’s intervention is only required to analyse the results). On the other hand, TDA/MKL partly failed to identify the modern domestic species goat and sheep. However, an in-depth analysis of the reasons for such difficulty revealed that TDA might be better suited for the intra-specific classification of such species, for which our method seems likely to perceive a lot of detail. Another drawback of TDA methods is their lack of explainability: it is not possible to know which part of the bone contributed the most to the identification. In light of the above remarks, a few avenues for future research can be outlined: (i) create ad hoc datasets of 3D scans to assess the capabilities of TDA/MKL intra-specific classifiers; (ii) combine TDA features to anatomical criteria and GMM morphological patterns; (iii) test deep learning methods on an increasing dataset to test others classification approaches. In conclusion, this research demonstrates the effectiveness of the TDA/MKL methods in extracting authentic biological information that can be interpreted by archaeozoologists, as well as novel information that cannot be detected by traditional anatomical criteria. In this sense, these approaches represent a milestone in the dialogue between two scientific disciplines, mathematics and bioarchaeology, by providing information comparable to that obtained by palaeogenetics.