A deep learning approach for the analysis of birdsong

Many birds included in this study were previously recorded and analyzed in Garcia-Oscos et al., 2021. This includes seven birds which were injected with a pscAAV-GFP-shFoxP1 virus before exposure to a song tutor, leading to disrupted songs (referred to as FP1 KD birds in this manuscript and FP1-KD SE in Garcia-Oscos et al., 2021). Eight birds were included in this group in the previous paper. One was omitted from this manuscript because it exhibited completely typical song, likely due to weak viral expression. A further 10 birds which were injected with a control virus before exposure to a song tutor (Ctrl SE in Garcia-Oscos et al., 2021), and 10 which were injected with the pscAAV-GFP-shFoxP1 after tutor song exposure (FP1-KD BI in Garcia-Oscos et al., 2021) are included in the current study. Both of these groups exhibit species-typical song production and are included in the ‘typical’ group in this study. Finally, 8 additional birds which were raised in isolation from an adult song model until they were at least 90 days post-hatch (‘Full isolates’ in the FP1 paper and ‘Isolates’ in the current study) were included in this study, for a total of 35 birds. See Garcia-Oscos et al., 2021 for more information on viral injections and rearing conditions.

An additional 8 adult isolate birds, 37 typically reared adult birds, and 10 juvenile birds were recorded for this study. For isolated birds, fathers were removed from breeding cages before the young reached 12 dph. These young remained housed with their mother and siblings in a room containing only other isolate breeding cages, until they were weaned between 40 and 60 dph, at which point they were housed individually, with auditory but no visual access to other isolate reared males. Typically, reared and juvenile birds were raised in our main colony room which contained about 55 breeding pairs. They had unlimited access to their father’s song in their home cages, until they were weaned between 40 and 60 dph, at which point they were housed in group cages with other males.

Recordings were obtained by individually housing birds in sound-attenuating chambers. They were continuously recorded using Sound Analysis Pro 2011 (Tchernichovski et al., 2000). All birds were placed on a 14 hr:10 hr day:night cycle and provided ad libitum access to food, water, and grit. All procedures were performed in accordance with protocols approved by the Animal Care and Use Committee at UT Southwestern Medical Center.

In addition to the birds recorded at UTSW, this study includes recordings of 25 pupils and 6 tutors from the Rockefeller University Field Research Center Song Library (Tchernichovski et al., 2021), 5 juvenile birds from Duke University (Brudner et al., 2023), and 5 early deafened and 4 sham deafened birds from Hokkaido University (Mori and Wada, 2015). See citations for more information on rearing and recording conditions.

A random subset of 30 song files from a single day of recording was annotated for each of 35 adult birds from UTSW, and 15 song files were annotated for each of 25 adult birds from the Rockefeller Song Library. Manual annotation was performed using the evsonganaly application in MATLAB (Tumer and Brainard, 2007), and involved (1) amplitude threshold syllable segmentation with a threshold selected for each song file based on visual inspection of the amplitude trace and spectrogram, (2) manual correction of erroneous syllable onsets or offsets, and (3) assignment of syllable labels to each syllable based on visual inspection of the spectrogram. All annotations were prepared by one of two expert annotators. These annotators consulted with each other in cases of ambiguous syllable segmentation or labeling and used the same existing set of labeled songs as a reference. Both annotators agreed to adopt the more conservative approach of splitting syllables into smaller segments in cases of ambiguous segmentation and assigning different syllable labels in cases of ambiguous labeling.

For all applications except training TweetyNet (Cohen et al., 2022), segments that reflect cage noise were dropped from the annotations based on visual inspection of the spectrograms. A second set of annotations was made which retained noise segments, labeling them as such, with all syllables and calls labeled as simply ‘vocalizations’ for the purpose of training TweetyNet.

Amplitude segmentation was performed using the ‘RMSEDerivative’ class in the AVN Python package’s segmentation module. Each song file is bandpass filtered between 200 Hz and 9000 Hz, then the root mean square energy (RMSE) of each audio frame is computed with a hop length of 512 samples and a frame length of 2048 samples. The RMSE values of each song file are normalized, then the RMSE’s first derivative is compared against user-specified thresholds. A syllable onset is identified as a positive crossing of the ‘onset’ threshold. Syllable offsets tend to be marked by more gradual changes in RMSE compared to syllable onsets, making it difficult to identify them consistently. To mitigate this, we perform onset to onset segmentation with this method, meaning each segment included a song syllable and the silent gap that immediately followed it. If a syllable onset is not followed by another onset within 300 ms (as in the end of a song bout), the offset is set as the first negative crossing of an ‘offset’ threshold after the syllable onset.

In keeping with the TweetyNet and WhisperSeg segmentation methods which don’t require per-bird parameter adjustments, the same ‘onset’ and ‘offset’ thresholds were used for all birds in the dataset. These thresholds were selected using AVN’s segmentation.Utils.threshold_optimization_many_birds() function, which compares the F1 scores relative to manual segmentation obtained with multiple different threshold values to identify the threshold value that results in the lowest mean F1 score across all 35 UTSW birds used for amplitude segmentation validation. The same thresholds selected based on the UTSW birds were used to segment the 25 Rockefeller Song Library birds, as a test of the generalization of this method without the need for manual segmentations.

The vak Python package was used to prepare datasets for, train, and generate segmentation predictions with the TweetyNet model (Cohen et al., 2022). TweetyNet is a deep neural network consisting of a block of convolutional layers followed by a bidirectional long short-term memory (LSTM) layer. The model takes a 1 s spectrogram of song as input, and labels each frame within that spectrogram. TweetyNet was designed for simultaneous syllable labeling and segmentation, in which case it would label each frame of the spectrogram with a syllable label or as silence. However, to make this model generalize to new birds without any additional training data, we instead trained TweetyNet to label each frame as a vocalization, silence, or noise (common sources of noise include the bird hopping around its cage and flapping its wings), rather than a more specific syllable type. When trained with such data from many birds, it can learn to distinguish vocalizations from noise and silence in a sufficiently general manner that the model can be applied to previously unseen individuals.

Manually annotated song files with label classes ‘noise’ and ‘vocalization’ were used to train TweetyNet in a leave-one-out cross-validation scheme, meaning the model was trained with data from all but one bird and tested on the withheld bird for each of the 35 birds in the UTSW dataset. A model trained with data from all 35 UTSW birds was used to segment the 25 validation birds from the Rockefeller Song Library dataset, to test the model’s ability to generalize to new colonies. Full model training and prediction procedures can be found in this paper’s accompanying GitHub repository. For more information on the TweetyNet model itself, see Cohen et al., 2022.

WhisperSeg is an instance of the Whisper Transformer model which was pre-trained for automatic human speech recognition and fine-tuned for animal voice activity detection with a multi-species animal vocalization dataset (Gu et al., 2023). It takes a spectrogram representation of up to 2.5 s of song as input and outputs the indices of vocalization onsets and offsets in the spectrogram. These indices are then converted to timestamps, and a consistent labeling scheme for an entire song file is achieved through a ‘majority-vote’ post-processing step across overlapping 2.5 s song segments. Syllable segmentation was performed using the whisperseg-large-ms-ct2 model, with hyperparameters optimized for zebra finch song segmentation, based on Gu et al., 2023. This was done using scripts shared in this paper’s accompanying GitHub repository, which must be installed separately from AVN according to the instructions provided in the repository. For more information on the WhisperSeg model and its training, see Gu et al., 2023.

Segmentation methods were compared on the basis of their precision, recall, and F1 scores relative to manual annotations.

where a true positive is a syllable onset in the automatic segmentation that is within 10 ms of a syllable onset in the manual annotation, a false positive is a syllable onset in the automatic segmentation that doesn’t match with a syllable onset in the manual annotation within 10 ms, and a false negative is a syllable onset that is present in the manual annotation which doesn’t match with an automatic segmentation onset within 10 ms. When determining onset alignments, we ensure that each syllable onset in the manual annotation can only ‘match’ with a single syllable onset in the automatic segmentation, and vis-versa. This was done using AVN’s ‘segmentation.Metrics.calc_F1()` function. Across all three metrics, scores closer to 1 indicate better agreement between automatic and manual segmentations. These same features were calculated for syllable offsets as well, but with an allowance of 20 ms rather than 10 ms, to account for the greater variability in exact offset segmentation across all methods tested.

To further examine the temporal precision of each method relative to manual annotation, we also calculated the time difference in milliseconds between matched syllable onsets and offsets between automatic segmentations and manual annotations. This was done using AVN’s ‘segmentation.Metrics.get_time_delta_df()’ function.

UMAP dimensionality reduction (McInnes et al., 2018) is performed on spectrograms of syllables, prior to HDBSCAN clustering for label assignment (McInnes and Healy, 2017). First, spectrograms of each segmented syllable are produced. The audio is first bandpass filtered between 500 Hz and 15 KHz, then amplitude normalized independently for each syllable rendition. The short-term Fourier transform of the normalized audio for each syllable is computed with a window length of 512 samples and a hop length of 128 samples. The resulting amplitude spectrogram is converted to decibels using the librosa ‘amplitude_to_db()’ function (McFee et al., 2023). The db-scaled spectrogram is then padded to match the dimensions of the longest syllable in a given bird’s dataset, or clipped to 870 ms if it exceeds that duration. This limit is very generous, and only ever applies to segmentation errors, but it is necessary to avoid memory issues during UMAP computation. Normalized spectrograms are then flattened from an array to a single long vector, and the vectors corresponding to each spectrogram are concatenated into an array separately for each bird. This spectrogram array is used to calculate the UMAP embeddings of each syllable using the UMAP Python package’s ‘UMAP()’ function. This approach is based on Sainburg et al., 2020.

At a high level, UMAP dimensionality reduction involves constructing a graphical representation of the syllable set, where each syllable spectrogram can be thought of as a point in high dimensional space which is connected to other syllables near it by edges. These edges are weighted based on the distance between points and the local density of those data points. The high-dimensional graph is then projected into lower dimensions in a way that best preserves its overall structure.

UMAP dimensionality reduction can be a useful initial step when attempting to cluster high-dimensional data points because many clustering algorithms, especially density-based clustering algorithms such as HDBSCAN, can suffer from the ‘curse of dimensionality’. When clustering spectrograms directly, each pixel in the spectrogram is a dimension, meaning each spectrogram exists as a point in a space with thousands of dimensions. In such a high-dimensional space, points will be very sparsely distributed, even if the spectrograms appear largely very similar. As a result, it is very difficult to detect regions of higher point density to serve as the basis of clusters. Reducing the dimensionality of the dataset forces points closer together, such that regions of high density separated by lower density can be more easily detected. UMAP is particularly adept at emphasizing local clusters in high-dimensional data because of how its initial embedding graph is constructed.

UMAP parameters were selected based on suggestions in the UMAP-learn documentation for clustering using UMAP embeddings and based on visual inspection of plots and labeling outcomes compared to manual annotations for birds from the UTSW dataset.

The ‘Hierarchical Density-based Spatial Clustering of Applications with Noise’ (HDBSCAN) clustering algorithm (McInnes and Healy, 2017) was applied to the UMAP embeddings of syllable spectrograms for each bird independently, in order to assign syllable labels. This method was selected based on the results in Sainburg et al., 2020 for clustering Bengalese finch syllables, a species closely related to zebra finches.

Essentially, HDBSCAN works by calculating the ‘mutual reachability’ distance between points in the UMAP space, based on the distance between them and their local densities. These mutual reachability distances serve as edges connecting nodes (points representing individual song syllables) in a graph, which are then pruned to obtain a minimum spanning tree (a graph using the minimum number of total edges to connect all points). The minimum spanning tree is then converted to a hierarchy by sorting the edges on the basis of their mutual reachability scores. Clusters of points are identified by defining a minimum cluster size and selecting the clusters that persist over the longest span of the hierarchy. As with the UMAP parameters, the same HDBSCAN hyperparameter set was used for all birds. The hyperparameter values were selected based on v-measure scores and visual inspection of confusion matrices for WhisperSeg segments compared to manual annotations for birds from the UTSW colony.

Syllable labeling was assessed by comparing automatically assigned syllable labels to manual annotations. Automatically labeled segments first had to be aligned to the manual annotations to identify pairs of labels in the automatic clustering and manual annotation that referred to the same vocalization. This was done using the same method described in the Segmentation Validation section, in which syllable onsets are uniquely matched to their closest counterpart across segmentation methods, this time up to a maximum distance of 100 ms. False positive syllable detections (i.e. syllable present in the automatic segmentation without a manual annotation counterpart) are assigned to their own manual annotation category (‘x’), and false negative syllable detections (i.e. syllables present in the manual annotation without an automatic segmentation counterpart) are assigned to their own cluster (‘1000’) for the purposes of visualization and quantification.

Once syllables have been aligned between automatic segmentation and manual annotations, the HDBSCAN cluster labels are compared to manual labels for each bird to construct a confusion matrix, which gives the number of syllables in each HDBSCAN cluster that carry each of the possible manual labels. The confusion matrix values can then be used to compute homogeneity, completeness, and v-measure scores, to evaluate the correspondence between HDBSCAN labels and manual annotations for each bird. Homogeneity measures the extent to which syllables with the same AVN label also carry the same manual annotation label, completeness measures the extent to which syllables with the same manual annotation label also carry the same AVN label, and the V-measure is the harmonic mean of these two scores.

where $H\left(\text{manual labels}|\text{clusters}\right)$ is the conditional entropy of the manual labels given the cluster labels, $H\left(\text{manual labels}\right)$ is the entropy of the manual labels, and vis-versa. In all cases, a higher score indicates better correspondence between clusters and manual labels, with a maximum possible score of 1 and minimum score of 0.

As UMAP embeddings are inherently stochastic, the subsequent HDBSCAN labels can vary over repeated runs with the exact same syllable dataset, depending on the random initialization of the UMAP embedding. To examine how consistent the labels are across different random UMAP initializations, we calculated the v-measure scores for 35 birds from the UTSW colony across 30 different random initializations and calculated the standard deviation of the v-measure scores for each bird, as a measure of the outcome variability due to UMAP’s stochasticity.

Beginning with a table of all AVN labels, syllables that are preceded and followed by a period of silence longer than 200 ms are removed, as they likely reflect calls produced outside of song. Song bouts are then identified as sequences of at least two syllables that are separated by silent gaps no longer than 200 ms. These bouts are aligned based on a user-specified alignment syllable, such that the first instance of the alignment syllable is in the same position across all bouts. This alignment is important as bouts typically begin with a variable number of introductory notes, which will obscure patterns in syllable sequence across bouts when they are not aligned to the first non-introductory note syllable. After alignment, bouts are ordered such that bouts with similar sequences after the alignment syllable are together in the final plot, which also helps emphasize patterns across bouts. This is done using AVN’s syntax. Syntax_Data.make_syntax_raster() function. See avn’s documentation for additional information and examples.

As with syntax raster plots, syllables that are preceded and followed by more than 200 ms of silence are dropped from the AVN labels as they likely reflect calls produced outside of song. Silent gaps longer than 200 ms and file bounds are then added as states to the AVN label sequence. All syllable transitions between AVN labels are then counted, including transitions to and from periods of over 200 ms of silence; meaningful transitions as they reflect the beginnings or ends of song bouts. Transitions to and from file bounds are ignored, as these are artifacts of the recording and don’t reflect meaningful behavioral states. The transition counts are then divided by the total number of renditions of the first syllable type in the transition to get the conditional probability of the second syllable, given the first syllable. This is done using AVN’s syntax.Syntax_Data.make_transition_matrix() function.

Syntax stereotypy is quantified using the entropy rate of the syllable transition matrix.  

where $p_{i,k}$ is the probability of transitioning from initial syllable type i to following syllable type k, and ${\pi }_i$ the probability of syllable type i occurring, regardless of what syllable precedes or follows it. An entropy rate approaching 0 indicates that all transitions are highly predictable. The maximum possible entropy rate score is ${log}_2\left(N\right)$ where N is the number of syllable types in the bird’s repertoire plus one to account for silence as a possible state. To directly compare scores between birds without being biased by the number of syllable types in their song (a feature which depends strongly on the number of syllable types present in their tutor’s song), we divide the entropy rate score by ${log}_2\left(N\right)$ such that it is now bounded between 0 and 1. This is done using AVN’s syntax.Syntax_Data.get_entropy_rate() function.

A repetition bout refers to every instance in which a syllable is produced, either a single time or multiple times in a row. For example, in the syllable sequence abcaaabc, syllable a has two repetition bouts, one of length one, meaning the syllable was produced without being repeated, and one of length 3, meaning the syllable was produced three times in a row. The number and length of repetition bouts are calculated for each syllable type in a bird’s repertoire. The mean repetition bout length and coefficient of variation (CV) of repetition bout length are then calculated for each syllable type.

To facilitate comparisons across birds, which have different numbers of syllable types, the mean repetition bout length and CV of repetition bout length for the syllable type with the highest mean repetition bout length are selected to represent the bird’s overall tendency to repeat syllables, excluding syllable types that reflect putative calls or introductory notes. Typical zebra finches often repeat calls or introductory notes, but rarely repeat song syllables, so looking at the repetition of song syllables is more informative when detecting or comparing birds with abnormal syntax. That said, in certain experiments, repetition bout features of calls or introductory notes may be of greater interest, in which case they can also be specifically identified using AVN.

An AVN syllable type is considered a putative introductory note if it (1) is no less than 5% less likely to be transitioned to from silence than the syllable type most commonly transitioned to from silence, meaning it tends to occur at the start of a vocalization bout and (2) it has a single dominant transition to a syllable type other than itself which is not silence, meaning that after a number of repetitions, it is eventually followed by a predictable next syllable type, which should reflect the start of a motif. These criteria were determined based on inspection of the syntax properties of introductory notes in manual song annotations. An AVN syllable type is considered a putative call if it is (1) not a putative introductory note and (2) produced in a bout of one or two syllables preceded and followed by at least 200 ms of silence in more than ⅓ of all utterances. This criterion was again determined based on visual inspection of manual song annotations.

A syllable duration distribution is constructed based on the segment durations output by WhisperSeg for each bird. A histogram of the ${\log }_{10}$ of syllable durations is calculated, with 50 evenly spaced bins ranging from –2.5 to 0. As in Goldberg and Fee, 2011, the log of syllable durations is used because the syllable duration distributions of juvenile birds are roughly exponential and therefore linear in log space. Histograms are normalized to produce a probability density function across syllable durations. The entropy of this distribution is then calculated as

where $p_i$ is the density the ith bin in the histogram, and N is the total number of bins (50, in this case). The resulting entropy can range from 0 to 1, with higher scores indicating less predictable syllable durations, consistent with the songs of immature birds.

Entropy is calculated similarly for silent gap durations, where gap durations are defined as the time difference between a syllable offset and the immediately following syllable onset, up to a maximum duration of 200 ms. A log transform was not applied to gap durations before constructing a histogram with 20 10 ms bins.

Rhythm spectrograms are a visualization of the strength and stereotypy of rhythmic patterns in a bird’s song, generated by concatenating the rhythm spectra of multiple song bouts, as first proposed in Saar and Mitra, 2008. A song bout’s rhythm spectrum is the spectrum of the first derivative of its amplitude. If a song’s amplitude has consistent repeating fluctuation patterns (as we expect for a bout composed of multiple repetitions of the same stereotyped motif), then its spectrum will exhibit harmonic banding patterns. If, by contrast, there are no repeating rhythms in the song’s amplitude, the rhythm spectrum will have a more even spread of energy across frequency bands. To detect these harmonic patterns more easily in the rhythm spectrogram, all rhythm spectrograms are plotted as the rolling average of 10 song bouts, smoothing out some bout-to-bout variation in the spectra to make harmonic bands more obvious.

To ensure consistent dimensions and resolution across song bouts, the rhythm spectrum is calculated for segments of song of a fixed duration, rather than complete song bouts. This also eliminates the need for any segmentation and labeling of song files to identify bouts, making this timing analysis method completely independent of possible segmentation and labeling errors. Each .wav file is broken into multiple 3-s-long frames, with a hop length of 0.2 s. The three frames with the highest total amplitude (i.e. the three windows containing the most vocalizations) from each file have their rhythm spectra calculated, and the mean of their spectra is taken as the rhythm spectrum for that file. Because of this windowing system, only files at least 3+3 ×0.2 s in duration can be windowed this way, so shorter .wav files are ignored.

The derivative of the amplitude of each frame is centered at 0 by subtracting the mean value, then multiplied by a Hanning window to reduce spectral leakage when calculating the spectrum. The transformed amplitude derivative is then padded to a total length of 100000 frames, resulting in a smoother spectrum with more interpolated values. A bandpass filter is then applied, keeping only frequency components above 1 Hz and below 500 Hz, as these are the frequencies consistent with typical zebra finch motif and syllable periods. Finally, the real component of the Fourier transform is calculated, constituting the frame’s ‘rhythm spectrum’. Only portions of the rhythm spectrum corresponding to frequencies between 0 and 30 Hz are included in rhythm spectrograms and downstream feature calculations, as this is the range with the strongest harmonic banding for typical zebra finches. This is all done using AVN’s avn.timing.RhythmAnlysis.make_rhythm_spectrogram() function.

We quantify the strength of the harmonic content of a bird’s rhythm spectrum (i.e. the strength of its rhythm) by calculating the Wiener Entropy of the mean rhythm spectrum across bouts. Wiener entropy is a common acoustic feature used to assess the harmonic nature of zebra finch syllables, with scores near 0 reflecting signals with little harmonic structure, and scores ranging to negative infinity for signals with more harmonic structure.

This is calculated using AVN’s avn.timing.RhythmAnalysis.calc_rhythm_spectrogram_entropy() function.

Whereas the rhythm spectrum entropy measures the overall strength of the rhythms in a set of songs, the peak frequency variability reflects the consistency of the rhythm across multiple song renditions. The exact spacing of the harmonics in a rhythm spectrogram depends on the shape of the amplitude trace of the bird’s motif. It isn’t obvious how different motifs and motif lengths affect the banding pattern of the rhythm spectrogram, so it doesn’t make sense to compare the appearance of different birds’ rhythm spectra beyond the prominence of harmonic bands. Likewise, the frequency of the harmonic band with the highest magnitude doesn’t carry any special meaning. However, in a very stereotyped bird, that harmonic band will be consistent across songs. If the bird sings its song slightly faster or slower, the band can shift slightly in the frequency domain. So, we measure a bird’s rhythm stereotypy by looking at the variability of the frequency with the highest magnitude across song files (the peak frequency).

In practice, the frequency band with the highest energy can jump between harmonic bands across files, even while the overall timing is largely unchanged, so to truly capture fluctuations in the underlying rhythms in a bird’s song, we restrict the range of ‘peak frequency’ values to a 3 Hz band about the median peak frequency across bouts. The CV of the peak frequency within this range is calculated as the peak frequency variability. This is done using AVN’s avn.timing.RhythmAnalysis.calc_peak_freq_cv() function.

Spectrograms of manually segmented and labeled syllables from 21 adult zebra finches from the UTSW colony were used for model training. All validation was performed with spectrograms of WhisperSeg segmented syllables from a test set of 30 tutor-pupil pairs from UTSW and 25 tutor-pupil pairs from the Rockefeller Song Library (Tchernichovski et al., 2021), none of which were included in training. For the triplet loss model, these spectrograms were normalized for amplitude, then clipped or padded to a uniform duration of 180 ms. To reduce computational costs, all frequency bands below 2 kHz and above 6 kHz are discarded. For the variational autoencoder model, spectrograms were computed for each syllable then interpolated to consistent time and frequency representations; 128 frequency bands ranging from 400 Hz to 10 kHz, and 128 time bins. Syllables longer than 300 ms were discarded. Short syllables were stretched in time by a factor of $\surd \frac{300ms}{t}$, then symmetrically zero-padded to a size of 128 time bins, such that all syllables’ spectrograms had the same dimensions.

The proposed neural network model is composed of five convolutional layers alternating with four ‘Multiscale Analysis Modules’ (MAMs), followed by a global pooling layer and three fully connected linear layers (Figure 6—figure supplement 3). This architecture is based on the model proposed in Thakur et al., 2019, which was used for species classification with field recordings of bird, frog, and toad vocalizations. The first convolutional layer consists of 32 3x3 kernels, with the 4 subsequent convolutional layers consisting of 64 3x3 kernels with a stride length of 2 along the frequency axis, resulting in downsampling by a factor of 2 along that dimension. Each MAM is composed of four parallel strands, each processing the data at different scales. The first strand consists of a single convolutional layer with 32 1x1 kernels, the second, third, and fourth strands start with a 32 filter 1x1 kernel convolutional layer, followed by a convolutional layer with 32 3x3, 5x5, and 7x7 kernels, respectively. The output of each of these strands is concatenated channel-wise, resulting in a 128-channel representation of the data which is passed to the next layer. This parallel strand organization allows the model to perform feature extraction at multiple different scales without increasing the depth of the model, saving computational cost and limiting potential overfitting. This approach was first proposed in Szegedy et al., 2015. The ReLU activation function is used after every layer (Nair and Hinton, 2010). The output of the final linear layer is an 8-dimensional vector, which represents an input syllable’s embedding. These vectors are normalized to have a length of 1, such that all embeddings lie on a unit eight-dimensional hypersphere.

The model is trained using dynamic triplet loss with triplet mining. Triplet loss involves presenting the model with batches of triplets, where each triplet consists of an anchor, a positive, and negative syllable. The anchor and positive syllables carry the same manual annotation label from the same bird, and the anchor and negative syllables carry different labels, either from the same bird or from an unrelated bird. The loss function to be minimized is:

where N is the total number of possible triplets in training, $f\left(A_i\right)$ is the embedding of the anchor in triplet i, $f(P_i)$ is the embedding of the positive, $f(N_i)$ is the embedding of the negative, and α is a margin parameter. If the positive is closer to the anchor than the negative by at least α, the loss for that triplet is 0.

During training, many randomly sampled triplets will already yield a loss of 0, and therefore will not lead to any change in the model. As a result, such triplets are not presented to the model during training. The remaining triplets that result in a positive loss value can be divided into two groups: hard triplets and semi-hard triplets. Hard triplets are cases where $||f\left(A\right)-f\left(P\right)||>||f\left(A\right)-f\left(N\right)||$, and semi-hard triplets are cases where $||f\left(A\right)-f\left(P\right)||<||f\left(A\right)-f\left(N\right)||$ and $||f\left(A\right)-f\left(P\right)||-||f\left(A\right)-f\left(N\right)||>\alpha$. Hard triplets result in higher loss values and, therefore, larger model weight updates, which can result in unstable training. Previous studies have shown that, as a result, training with semi-hard triplets alone can lead to faster and more stable model convergence (Schroff et al., 2015). We found that we achieved the best model performance when training with a ratio of 75 semi-hard triplets: 25 hard triplets, so this ratio was used to train the final model. A forward pass of the model is performed for each mini-batch in training to determine the ‘hardness’ of all possible triplets from that batch. Triplets are then sampled according to the specified ratio of semi-hard to hard triplets and are presented to the model for training.

Our approach is said to be dynamic triplet loss because the value of the margin parameter $\alpha$ is updated dynamically over the course of training. The value is initially set to 0.1 and increased stepwise by 0.2 every time a training epoch had fewer than 2500 hard or semi-hard triplets per batch on average, up to a maximum value of 0.7. This allows the model to begin learning an easier task, of separating syllables of different classes by a smaller margin. As the margin increases, the task gradually becomes more difficult, leading to more stable model convergence compared to starting with a higher margin value. Each time the margin parameter is increased, triplets that previously resulted in a loss of 0 can become semi-hard triplets, meaning the set of triplets presented to the model also expands over the course of training. The initial margin value, maximum margin value, margin step size, and non-zero triplet threshold were all determined empirically. The weight optimization was performed with the Adam optimizer, with weight decay of 0.0001 (Kingma and Ba, 2014).

Ultimately, this dynamic triplet loss training constitutes a form of deep metric learning, where high-dimensional inputs (e.g. spectrograms of syllables) are mapped onto a lower-dimensional space where the similarity between samples is proportional to the distance between them. Training with triplets is advantageous in a context where the total amount of labeled training data is limited, as the number of possible triplets is proportional to the cube of the number of training samples.

In addition to using this triplet loss model, we also calculated similarity score using syllable embeddings generated by a variational autoencoder (VAE) model, as described in Goffinet et al., 2021. Briefly, using the autoencoded-vocal-analysis Python package (Goffinet et al., 2021), we trained the base VAE model with an eight-dimensional latent embedding using syllable spectrograms from the 21 UTSW zebra finches used to train the triplet loss model.

Syllable embeddings were obtained by running a forward pass of the trained triplet loss or VAE model with all segmented syllables that an individual bird produced on a given day. This results in a set of thousands of syllables in the 8-dimensional embedding space. Two bird songs are compared by calculating the Earth Mover’s Distance (EMD) between their syllable embeddings using the PyEMD package (Doran, 2014). If one were to imagine the syllable embedding empirical distributions as piles of dirt, the earth mover’s distance is the minimum cost of moving the earth from one distribution to match the other, where cost is defined as the amount of dirt moved multiplied by the distance over which it is moved. This value can range from 0 for identical distributions to positive infinity for distributions that are infinitely far apart. As our triplet loss model’s embedding space is limited to points lying on a unit radius 8-dimensional hypersphere, the maximum possible value for EMD is 1.41 for distributions that are each concentrated on a single point, maximally separated within the constraints of the embedding space.

The EMD score considers many song renditions from each bird being compared, allowing a better overall comparison of the similarity between two birds’ song production as compared to multiple pairwise comparisons between renditions. One limitation of EMD, however, is that it is completely agnostic to syllable sequencing, so a pupil that imitated all syllables from his tutor but sings them in a completely different order will have a similar EMD score to a pupil that imitated all syllables from a tutor and produces them in the same order. That said, there is new evidence that zebra finches recognize songs independently of syllable order, raising questions about the importance of syllable order in song perception (Ning et al., 2023). The EMD score is also symmetrical, meaning that the presence of syllables in the tutor’s song that weren’t imitated by the pupil will have the same impact on the EMD score as new, improvised syllables present in the pupil’s song but not in the tutor’s song.

In addition to EMD, we calculate the dissimilarity between distributions of embedded syllables using maximum mean discrepancy (MMD), as in Goffinet et al., 2021. MMD compares two distributions of syllables by computing the squared distance between their mean features mapped in a reproducing kernel Hilbert space (RKHS) using a Gaussian radial basis function (RBF) kernel, with a kernel bandwidth equal to half the median pairwise distance between 1000 randomly sampled syllable embeddings from our validation dataset. The result is a measure of the dissimilarity between distributions, where a larger score indicates greater dissimilarity between the two syllable distributions. As with EMD, this is a symmetrical score that doesn’t consider the syllable sequencing similarity between the two birds being compared.

To validate the performance of our models and EMD or MMD scores for similarity scoring, we compute EMD and MMD scores between a pupil and itself, a pupil and its tutor, pupil and a ‘sibling’, and a pupil and an unrelated bird. For comparisons between a pupil and itself, embeddings of all recorded syllables from a day of song are computed. If the bird has fewer than 4000 syllables in this dataset, the dataset is randomly split in half and the two halves are compared. If a bird has more than 4000 syllables, two sets of 2000 syllables are randomly sampled and compared. For comparisons between a pupil and its tutor, up to 4000 syllables are sampled from each of the tutor and the pupil and compared. In the case of pupil and ‘sibling’ comparisons, a sibling is defined as another bird sharing the same song tutor. We expect that typical zebra finches that learned from the same tutor will have similar songs, but that these will generally be less similar than a pupil compared to its tutor directly. Each pupil is compared to up to three ‘siblings’, depending on availability in our dataset. Finally, each pupil is compared to three randomly selected pupils who don’t share their song tutor. For each of these comparisons, up to 4000 syllables are randomly sampled from each bird as well, to help reduce the compute time for EMD and MMD.

As in Mandelblat-Cerf and Fee, 2014, contrast index is calculated as:

where self-similarity is the EMD or MMD score between pupil and itself, calculated as described in the previous section, and cross similarity is the mean similarity between a pupil and 3 unrelated birds, again calculated as described in the previous section. As EMD and MMD are dissimilarity scores, rather than similarity scores, more negative values reflect better contrast between comparisons, so we report the absolute value for ease of comparison to existing similarity scoring methods.

Similarly to the contrast index score, the tutor contrast index is calculated as:

where tutor similarity is the EMD or MMD score between a pupil and its tutor, calculated as described previously, and cross similarity is the mean similarity between a pupil and three unrelated birds. As with the contrast index, absolute values are reported, with higher values indicating larger contrast in score between a pupil vs. tutor comparison and pupil vs. unrelated bird comparison for a given pupil bird.

A panel of 11 expert human annotators was each presented with 126 pairs of spectrograms and was instructed to rate their similarity on a scale from 1 (not similar) to 10 (very similar). The spectrograms were generated using Sound Analysis Pro 2011 (Tchernichovski et al., 2000), began at the beginning of a song bout, and included at least one full motif when motif structure was present. Raters were presented with 4 pupil-tutor spectrogram pairs per pupil, 8 comparisons between a tutor and itself to ensure that raters were using the full rating scale, and 10 duplicated tutor-pupil comparisons to ensure that the scorers were internally consistent. No individual scorer differed from the mean score by more than an average of 2 standard deviations, none differed by more than 2 points on the duplicated comparisons, and all but one made use of the full scale (this scorer never assigned a perfect 10/10 score, so their scores were rescaled such that they spanned the full range). Similarity scores for each tutor-pupil pair were obtained by taking the mean similarity score across their 4 spectrogram pairs, across all scorers. The full scoring set is available at https://forms.gle/9TDu1fwGGYXWKhgB6. These scores were previously generated for and published in Garcia-Oscos et al., 2021. Of the birds evaluated, 15 were also in the similarity scoring validation set, so the correlation between these 15 birds’ mean human similarity scores and tutor-pupil EMD and MMD scores was used to evaluate the agreement between methods.

For comparison to the expert human similarity scores, Sound Analysis Pro 2011 % similarity scores were calculated for the same set of 15 pupils. A representative motif from the tutor song was selected and compared to between 30 and 60 motif renditions from the pupil bird when the pupil was over 90 dph, using the asymmetric time-courses similarity tool. The final reported scores are the mean % similarity across all comparisons for a given pupil.

For the six birds from UTSW and five birds from Duke University (Brudner et al., 2023) from which we had recordings at 90-100d ph and earlier time points, we computed the MMD between their juvenile and adult songs. For each bird, 4000 WhisperSeg segmented syllables were sampled from a full day of song recordings when the birds were between 90 and 100 dph, to serve as the mature song distribution. Up to 4000 WhisperSeg segmented syllables were sampled from each day of available recordings prior to or shortly following the mature song date for comparison. MMD scores were calculated using embeddings from the triplet loss model as described previously. As the scores appeared to follow an exponential pattern, where the rate of song dissimilarity change slowed over development, we fit an exponential function to the data using the scipy.optimize curve_fit() function (Virtanen et al., 2020), and plotted this function alongside the data.

• Todd F Roberts

• Therese MI Koch

• Ethan S Marks