Deep Neural Network Identification of Limnonectes Species and New Class Detection Using Image Data

We are facing a biodiversity crisis. At present we have discovered, identified, and described only a very small fraction of the organisms with which we share this planet (e.g., \citeNPSweetlove2011). Unfortunately, because of habitat loss, climate change, and other anthropogenic influences, it is likely that we will lose this diversity faster than we can recognize it (e.g., \shortciteNPCowieetal2022). Our ability to accurately identify and delineate biological diversity is critically important because it influences our understanding of true biological richness and informs conservation efforts to safeguard these natural resources. Two critical tasks involved in understanding biological diversity are (1) the ability to correctly identify an individual organism to the taxonomic rank of species, and (2) the ability to detect individuals that may represents new species that are undescribed and therefore unrecognized by the scientific community. Traditionally, we have relied on the expertise of taxonomic specialists who drew from an understanding of morphology and natural history to identify new species. In recent decades, tools such as DNA sequencing and phylogenetic analyses have been employed by systematists in this pursuit. These new methods have advanced our rate of species discovery (e.g., \shortciteNPYAMASAKI2017177) and revealed species complexes in which morphological similarity has obfuscated true diversity (e.g., \shortciteNPFreudensteinetal2016).

We consider a species complex to be a group of two or more phylogenetically related species-level taxa that share such morphological similarity that group members cannot easily be distinguished from one another (\shortciteNPScherzetal2019 and \shortciteNPdeSousa-Paulaetal2021). Examples of species complexes can be found across all domains of life and present unique challenges to taxonomists and conservation biologists alike who would benefit from a new suite of tools that facilitate the recognition and detection of both known and unknown species.

In this paper, we focus on one such species complex comprising a group of Southeast Asian stream frogs. The frogs allied to the Limnonectes kuhlii (abbreviated as L. kuhlii) complex were considered to be a single species for about 200 years. Broadly distributed throughout Southeast Asia these frogs share a tremendous amount of morphological and ecological similarities. On the basis of DNA evidence alone, molecular phylogenetic studies (e.g., \shortciteNPYeetal2007, \citeNPMcLeod2008, \citeNPMCLEOD2010, \shortciteNPMatsuietal2010, and \shortciteNPStuartetal2020) recovered more than 22 unique evolutionary lineages representing potentially new species within the L. kuhlii complex. Using these data in combination with morphological evidence, subsequent taxonomic studies have described several new species from within the L. kuhlii complex (e.g., L. megastomias in \citeNPMcLeod2008, L. nguyenorum in \shortciteNPMcLeodetal2012, L. isanensis in \shortciteNPSuwannapoometal2017, and L. fastigatus in \shortciteNPStuartetal2020). Figure 1 shows four examples of frogs from four different clades (i.e., lineages).

Despite these advances, delineating species boundaries and describing additional new species within this complex has been confounded by the high degree of morphological similarity among the genetically distinct lineages. One of the diagnostic morphological characters that has proven useful in distinguishing between species in the L. kuhlii complex is the texture of the skin on the frog’s hind legs. This texture (created by the patterning of raised tubercles) is traditionally described qualitatively on the spectra of “smooth to rough” or “dense to sparse”, but without quantification of any kind. The pattern of tuberculation is consistent within a species but varies between species (\shortciteNPMcLeodetal2011, \shortciteNPMcLeodetal2012, and \shortciteNPStuartetal2020). The consistent yet subjectively qualitative nature of this character opens the door for applying machine-learning-based methods to standardize and automate species recognition and the detection of new species. Thus our specific research questions in this paper are: 1) whether machine-learning methods can use a qualitative character, such as the pattern of tuberculation on frog legs, to classify individuals into known classes (viz., species), and 2) whether a machine-learning-based tool can be developed to detect new classes (viz., new species). Figure 2 shows examples of the texture of the frog legs from the four clades for the four samples as shown in Figure 1, which highlights visible differences in tubercle structure on the hind limbs.

The first research question can be formulated as an image classification problem. Thus, deep-learning techniques can be applied. Deep learning has proven to be a powerful tool in resolving image classification problems across diverse applications (e.g., \shortciteNPKrizhevskyetal2012, and \citeNPRawatetal2017). \shortciteNNorouzzadehetal2018 applied a deep-learning method to monitor wildlife using camera-trap images, and \shortciteNWangetal2021 applied deep-learning algorithm on CT images to screen for COVID-19 disease. Whereas the application of machine learning is commonly used in species identification (e.g., iNaturalist, www.inaturalist.org, and PictureThis, www.picturethisai.com) and there has been some application of machine-learning and deep-learning methods to detect species from audio files of frog calls (e.g., \shortciteNPHuangetal2009, \shortciteNPHassanetal2017, and \shortciteNPXieetal2017), no study has applied these methods to problems associated with species complexes and new species detection. In this paper, we develop a convolutional neural network (CNN) model that can serve as an automatic identifier or classifier based on image/pattern recognition. This tool works in such a way that a photograph of an individual frog can be classified into one of the existing species using an existing set of images.

The second research question can be formulated as an out-of-distribution detection (OOD) problem (e.g., \shortciteNPhendrycks2017a, and \shortciteNPliang2018enhancing). The developed new class detection (NCD) tool allows for an image to be “rejected”, suggesting the detection of a potential new species. CNN and OOD techniques are well developed in the machine-learning discipline. For example, \shortciteNLeeetal2018 presents a simple unified framework for detecting OOD samples under the setting of adversarial attacks to neural network. The applications of CNN and OOD tools to biological research in the context of species complexes are new, and our approach provides quantitative tools to biological scientists working in the area. Our NCD methods extend the framework in \shortciteNLeeetal2018 and customize it for the frog NCD problem.

The rest of the paper is organized as follows. Section 2 describes the collection and preprocessing of the frog image data used in this study. Section 3 describes the CNN model for classification and describes a new method for NCD. Section 4 presents the classification and NCD results for the frog datasets and the NCD results for the MNIST dataset (\citeNPDeng2012) for further validation. Section 5 contains some concluding remarks.

In this section, we develop a classifier that can predict the species based on frog leg images. We use on deep neural networks (e.g., \citeNPGoodfellowetal2016). In particular, we use the convolutional neural networks (CNN) for the prediction that can provide classification based directly on the images of frog legs. In the literature, the CNN has been found to be a powerful tool in image classification. A CNN consists of an input layer, an output layer, and multiple hidden layers. The hidden layers consist of convolutional layers, pooling layers, and fully connected layers. Sometimes normalization layers and dropout layers are also used. More details can be found in Chapter 9 of \citeNGoodfellowetal2016.

Here we introduce the network architecture used for the frog images. Figure 4 illustrates the structure of the CNN used in this paper, which is a relatively simple structure due to the limited number of images, compared to the magnitude of tens of thousands of images in some other applications. We also tried CNN with more complicated structures. However, the prediction results are not getting better. For the CNN structure in Figure 4, the input image is ${\boldsymbol{X}}_{i}$, which is a $(3,512,512)$ tensor. The kernel size for the convolution is $(3,3)$. After the convolution, we do maximum pooling. The kernel size for the maximum pooling is $(2,2)$ with a stride of $(2,2)$ to reduce the image size. Between each convolutional layer, batch normalization is applied. For each fully connected layer, the random dropout rate is set to 0.8.

As shown in Figure 4, we have three convolutional layers after the input layer, which are labeled as Conv1, Conv2, and Conv3, respectively. The number of filters we use in Conv1, Conv2, and Conv3 are 32, 16, and 4, respectively. After application of Conv3, we obtain a $(4,62,62)$ tensor and it is flattened to a $15{,}376\times 1$ vector. The flattened vector then goes through four fully connected layers, labeled by FC1, FC2, FC3, and FC4, respectively. The lengths of those four output vectors are 256, 128, 64, and 4, respectively. For convenience, we denote the outputs of those fully connected layers by

$\displaystyle{\boldsymbol{z}}_{i}^{{\text{(1)}}},{\boldsymbol{z}}_{i}^{{\text{(2)}}},{\boldsymbol{z}}_{i}^{{\text{(3)}}},\text{ and }{\boldsymbol{z}}_{i}^{{\text{(4)}}},$

(1)

respectively. Here, $i$ is the index for the input image. In the output layer, as labeled “OP” in Figure 4, we use the softmax function to map ${\boldsymbol{z}}_{i}^{{\text{(4)}}}$ to a probability. That is

$\displaystyle p_{ij}=\frac{\exp(z_{ij}^{{\text{(4)}}})}{\sum_{k=1}^{m}\exp(z_{ik}^{{\text{(4)}}})},\quad j=1,\dots,m,$

(2)

where ${\boldsymbol{z}}_{i}^{{\text{(4)}}}=(z_{i1}^{{\text{(4)}}},\dots,z_{im}^{{\text{(4)}}})^{\top}$. Here, $j$ is the index for class label.

Let ${\mathsf{S}}$ be the set of indexes for samples used in model training and let ${\mathsf{T}}$ be the set of indexes for samples used in testing. Note that we do not have a validation set because of the small sample size and the fact that we do not have hyper-parameters in the model training. Let ${\boldsymbol{\theta}}$ be the set of all parameters in the CNN model. The categorical entropy is used as the loss function, that is,

$\displaystyle L({\boldsymbol{\theta}})=-\sum_{i\in{\mathsf{S}}}\delta_{ij}\log(p_{ij}),$

(3)

where $\delta_{ij}={\mathds{1}}(y_{i}=j)$ and ${\mathds{1}}(\cdot)$ is the indicator function. PyTorch (\shortciteNPPyTorch) is used for training the model to find the ${\boldsymbol{\theta}}$ that minimizes the loss function.

To measure the model prediction accuracy, we use $K$-fold cross validation. We have an unbalanced classification problem. If we split the data randomly across all four species, the training set may not have all four species. To ensure that all four species are in our training set and the testing set consists of equal proportions from each of the four species for each fold, we use a balanced cross validation. For each species with class label $j$, we denote the index set by ${\mathsf{C}}_{j},j=1,\dots,m$. The process for balanced cross validation is described as follows.

1. We randomly split ${\mathsf{C}}_{j}$ in to $K$ equal folds, which are denoted by ${\mathsf{C}}_{jk},k=1,\dots,K$. 2. For a given $k$, construct the training set as ${\mathsf{S}}=\cup_{j}\cup_{l\neq k}{\mathsf{C}}_{jl}$ and the testing set is as ${\mathsf{T}}=\cup_{j}{\mathsf{C}}_{jk}$. Then train the model based on set ${\mathsf{S}}$ and make predictions for samples in set ${\mathsf{T}}$. 3. Repeat Step 2 for $k=1,\dots,K$. 4. Repeat Steps 1 to 3 to obtain different random splits.

By using a balanced cross validation, we can ensure that each class always appears in both the training and testing sets.

In this section, we describe our algorithms for new class detection (NCD) for frog species. To evaluate the performance of our NCD algorithm, we take one of the classes in frog clades, say class $k$, as the samples from the out-of-distribution (OOD) population. Then the rest of the three classes are used to build the model.

For a trained deep-learning model, the outputs of the intermediate layers, such as those in (1), represent the extracted features from the sample. If a new sample does not belong to one of existing classes in the training set, then one would expect those extracted features would be different from those features from samples from existing classes. Thus, detecting new classes based on the outputs of intermediate layers can be treated as detecting outliers based on these multivariate outputs. Similar ideas have been used in the literature (e.g., \shortciteNPLeeetal2018). In this paper, we customize this idea for the our frog species detection problem and use parametric discriminant analyses for NCD.

Suppose we already have a trained deep neural network classification model, such as the one in Figure 4 (but with 3 outputs for FC4 and OP layers to match the dimension). Let ${\boldsymbol{z}}_{i}$ be one of those outputs in (1). For those three classes that used for modeling for the frog data, we split them into the training set (${\mathsf{S}}$) and the testing set (${\mathsf{T}}$). The samples that we set aside for OOD testing is denoted by ${\mathsf{O}}$. We prepare three datasets for the NCD problem.

• For the training set, we collect $\{{\boldsymbol{z}}_{i}\}$ for those $i\in{\mathsf{S}}$, which is the output of a fully connected layer. • For the in-distribution (ID) testing set, we feed the neural network with image ${\boldsymbol{X}}_{i}$ for those $i\in{\mathsf{T}}$ and then we collect $\{{\boldsymbol{z}}_{i},u_{i}\}$ for those $i\in{\mathsf{T}}$. Here $u_{i}$ is the OOD sample indicator. For those $i\in{\mathsf{T}}$, we have $u_{i}=0$. • For the OOD testing set, we feed the neural network with image ${\boldsymbol{X}}_{i}$ for those $i\in{\mathsf{O}}$ and then we collect $\{{\boldsymbol{z}}_{i},u_{i}\}$ for those $i\in{\mathsf{O}}$. For those $i\in{\mathsf{O}}$, we have $u_{i}=1$.

Then we assume ${\boldsymbol{z}}_{i}$ given class $j$ follows a multivariate normal distribution. That is, the distribution of $({\boldsymbol{z}}|y=j)$ is modeled as ${\rm N}({\boldsymbol{z}}|{\boldsymbol{\mu}}_{j},{\boldsymbol{\Sigma}}_{j})$, where the mean ${\boldsymbol{\mu}}_{j}$ is estimated by

$\displaystyle\widehat{{\boldsymbol{\mu}}}_{j}=$

$\displaystyle\frac{1}{n_{j}}\sum_{i\in{\mathsf{C}}_{j}\cap{\mathsf{S}}}{\boldsymbol{z}}_{i}.$

(4)

Here $n_{j}$ is the number of samples in ${\mathsf{C}}_{j}\cap{\mathsf{S}}$ and $n$ is the number of samples in ${\mathsf{S}}$.

For the OOD detection, we consider both linear discriminate analysis (LDA) and quadratic discriminant analysis (QDA). For LDA, the variance-covariance matrix $\Sigma_{j}$ are estimated by

$\displaystyle\widehat{{\boldsymbol{\Sigma}}}_{j}=\widehat{{\boldsymbol{\Sigma}}}=$

$\displaystyle\frac{1}{n-1}\sum_{j}\sum_{i\in{\mathsf{C}}_{j}\cap{\mathsf{S}}}({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})^{\top},$

(5)

which is common across all class labels. For QDA, each class has different $\widehat{{\boldsymbol{\Sigma}}}_{j}$, which is

$\displaystyle\widehat{{\boldsymbol{\Sigma}}}_{j}=$

$\displaystyle\frac{1}{n_{j}-1}\sum_{i\in{\mathsf{C}}_{j}\cap{\mathsf{S}}}({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})^{\top}.$

Note that the estimate for the variance-covariance matrix, $\widehat{{\boldsymbol{\Sigma}}}_{j}$, may not be positive definite. In such cases, we use the estimator from \citeNLEDOIT2004365, which works for non-positive definite cases by modifying the sample covariance matrix as,

$\displaystyle\widehat{{\boldsymbol{\Sigma}}}_{j}^{\text{LW}}=\rho_{1}{\boldsymbol{I}}_{n_{j}}+\rho_{2}\widehat{{\boldsymbol{\Sigma}}}_{j}.$

(6)

Here ${\boldsymbol{I}}_{n_{j}}$ in (6) is the $n_{j}\times n_{j}$ identical matrix, and $\rho_{1}$ are $\rho_{2}$ are scalars. Note that $\rho_{2}\widehat{{\boldsymbol{\Sigma}}}_{j}$ is a nonnegative definite matrix and adding $\rho_{1}{\boldsymbol{I}}_{n_{j}}$ to its diagonal terms will make $\widehat{{\boldsymbol{\Sigma}}}_{j}^{\text{LW}}$ positive definite.

Specifically, $\rho_{1}$ are $\rho_{2}$ are determined by the following optimization problem,

$\displaystyle(\rho_{1},\rho_{2})^{\top}=\operatorname*{argmin}_{\rho_{1},\rho_{2}}\mathds{E}\|\widehat{{\boldsymbol{\Sigma}}}_{j}^{\text{LW}}-{\boldsymbol{\Sigma}}_{j}\|^{2}.$

(7)

Here $\|\cdot\|^{2}$ is the squared Frobenius norm for matrices, and the expectation $\mathds{E}(\cdot)$ is taken with respect to the distribution ${\rm N}({\boldsymbol{z}}|{\boldsymbol{\mu}}_{j},{\boldsymbol{\Sigma}}_{j})$. The idea behind (7) is to find a matrix that is guaranteed to be positive definite and it is as close to ${\boldsymbol{\Sigma}}_{j}$ as possible.

Intuitively, those ${\boldsymbol{z}}_{i}$’s from the OOD samples should be far away from the center of distribution of ${\boldsymbol{z}}$, while those ${\boldsymbol{z}}_{i}$’s from the ID samples should be close to the center of distribution of ${\boldsymbol{z}}$, Thus, we use the Mahalanobis distance-based confidence score of ${\boldsymbol{z}}_{i}$ to measure the distance of ${\boldsymbol{z}}_{i}$ to its nearest class (e.g., \shortciteNPDenouden2018ImprovingRA). That is, for LDA, we compute

$\displaystyle w_{i}=\min_{j}\;({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})^{\top}\widehat{{\boldsymbol{\Sigma}}}^{-1}({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j}),i\in{\mathsf{T}}\cup{\mathsf{O}},$

and for QDA, we compute

$\displaystyle w_{i}=\min_{j}\;({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j})^{\top}\widehat{{\boldsymbol{\Sigma}}}_{j}^{-1}({\boldsymbol{z}}_{i}-\widehat{{\boldsymbol{\mu}}}_{j}),i\in{\mathsf{T}}\cup{\mathsf{O}}.$

For a sample with $w_{i}$ that is beyond a threshold, it will classified as a new class. So the data for the NCD problem is denoted as $\{w_{i},u_{i}\},i\in{\mathsf{T}}\cup{\mathsf{O}}$. We train a logistic regression model as

$\displaystyle u_{i}\sim\text{Bernoulli}(\xi_{i}),\quad\xi_{i}={\rm Pr}(u_{i}=1)=\frac{\exp(\alpha_{0}+\alpha_{1}w_{i})}{1+\exp(\alpha_{0}+\alpha_{1}w_{i})},$

(8)

where $\xi_{i}$ is the probability. With the logistic regression and the area under the receiver operating characteristic curve (AUC), a threshold can be determined. Note that to learn the logistic model, training samples from the OOD set are needed.

We first fit the CNN model as in Figure 4 to all the frog image data as summarized in Table 1. We use balanced cross validation as described in Section 3.1 with $K=5$, which leads to 80% for training and 20% for testing. Table 2 shows the accuracy of each fold and the overall accuracy for the testing set is around 73.1%, which is reasonably good, given that the size of training set is only around 150. To gain insights on mis-classifications, Table 3 shows the overall confusion matrix. We can see that the model mainly mis-classifies clades 4 and 5 to other clades, which may be caused by limited sample size.

We also plot some output of intermediate layers of the CNN model. Figure 5 plots an input image from one frog and the output of the last two convolution layers (i.e., Conv2 and Conv3) in the CNN model from the input frog image. We can see that different filters pick up different features from the image. Especially from Figure 5(c), we can see that the tuberculation of frog leg are activated in the CNN model, which provides insight that the tuberculation can contribute to classification of frog clades.

Using the CNN model built in Section 4.1, in this section, we conduct some OOD detections for the frog data with one OOD clade. The data splitting strategy is as described in Section 3.2. That is, we set one clade as the OOD sample and the rest to build the CNN model. Table 4 shows the fitting metrics based on the LDA on the frog data, with each species left out as OOD in each column. The accuracy of detecting OOD samples is high but also depends on the OOD sample size. Generally, more OOD samples will have satisfactory OOD accuracy. Figure 6 shows the visualization of the NCD for clade 5. We can see that the distribution of known species and new species is different. But for the layers with 3 and 64 units, there are no major gaps between the two distributions, which shows it is challenging to do NCD. One reason could be due to the small sample sizes of the OOD sets. Because of the limited sample size, we do not fit the QDA model because the covariance matrix can not be robustly estimated for each class.

In this section, we present additional results on OOD detection because some new samples were gathered, especially for clades 10, 11, 16, 18, 20, and 21. Because of the availability of additional images on multiple new clades other than clades 4, 5, 8, and 12, we can form an OOD class from multiple OOD clades from those new clades. In particular, the leg image counts for those new clades are 24, 9, 40, 61, 28, and 27, respectively. In total, we have 189 leg images for the OOD class. This section uses images from clades 4, 5, 8, and 12 to train the CNN model, as shown in Figure 4. Figure 7 shows the cumulative distribution of Mahalanobis distance of different fully-connected layers for the LDA and QDA methods. We see that the QDA procedure has slightly better classification power. Table 5 shows the metrics for LDA and QDA for OOD detection based on the frog data. We can see that the QDA method has better overall performance when considering AUC, true positive rate (TPR), and true negative rate (TNR). Overall, the OOD detection ability is good based on the multiple OOD clade data, though there is still room for improvement.

To better study the performance of our NCD methods, we further test our method by using the much larger MNIST image dataset in \citeNDeng2012, which has tens of thousands of images. The MNIST data are images of handwritten numbers. Although the frog images and MNIST images look different, the general idea of using the features extracted from CNN to do the OOD detection can be applied to different types of image data. Because the input image of the MNIST data is different (much simpler to some extent) from the frog images, we need to slightly modify our CNN model. In particular, we use the CNN architecture as shown in Figure 9, with fewer layers than the CNN for the frog data.

For NCD of the MNIST, we now can consider both classification methods, LDA and QDA. Figures 10 and 11 shows the distribution of the Mahalanobis distance using the MNIST data based on LDA and QDA where the number “9” is left out as OOD. The distribution of the layer’s output of OOD is significantly different from the training and testing data. We see that the QDA procedure has better classification power on layers with 320 and 50 outputs while LDA is better on the final fully connected layer with 8 outputs. One simple idea is when building a logistic model, use the QDA Mahalanobis distance from the layers with 320 and 50 outputs, and the LDA distance from the layer with 8 outputs. Table 6 shows the accuracy results for the LDA and QDA in terms of AUC, TPR, and TNR, based on the MNIST data. We can see that both LDA and QDA have good accuracy for detecting new classes, and QDA has better prediction accuracy than the LDA especially as measured by TPR. In addition, QDA’s performance is more stable when the OOD classes are changing.