INDIVIDUALITY IN THE HIVE - LEARNING TO EMBED LIFETIME SOCIAL BEHAVIOR OF HONEY BEES
Year: 2020

Abstract:
Honey bees are a popular model for complex social systems, in which global behavior emerges from the actions and interactions of thousands of individuals. While the average life of a bee is organized as a sequence of tasks roughly determined by age, there is substantial variation at the individual level. For example, young bees can become foragers early in life, depending on the colony’s needs. Using a unique dataset containing lifetime trajectories of all individuals over multiple generations in two honey bee colonies, we propose a new temporal matrix factorization model that jointly learns the average developmental path and structured variations of individuals in the social network over their entire lives. Our method yields inherently interpretable embeddings that are biologically plausible and consistent over time, which allows comparing individuals regardless of when or in which colony they lived. Our method provides a quantitative framework for understanding behavioral heterogeneity in complex social systems applicable in fields such as behavioral biology, social sciences, neuroscience, and information science.

heading: 1 INTRODUCTION
text: Animals living in large groups often coordinate their behaviors, resulting in emergent properties at the group level, from flocking birds to democratic elections. In most animal groups, the role an individual plays in this process is thought to be reflected in the way it interacts with group members. Technological advances have made it possible to track all individuals and their interactions in animal societies, ranging from social insects to primate groups (Mersch et al., 2013; Gernat et al., 2018; Mathis et al., 2018; Graving et al., 2019; Pereira et al., 2019). These datasets have unprecedented scale and complexity, but understanding these data has emerged as a new and challenging problem in itself (Pinter-Wollman et al., 2014; Krause et al., 2015; Brask et al., 2020).
A popular approach to understand high-dimensional data is to learn semantic embeddings (Frome et al., 2013; Asgari & Mofrad, 2015; Camacho-Collados & Pilehvar, 2018; Nelson et al., 2019). Such embeddings can be learned without supervision, are interpretable, and are useful for accomplishing downstream tasks. Individuals in animal societies can be described with semantic embeddings extracted from social interaction networks using matrix factorization methods. For example, in symmetric non-negative matrix factorization (SymNMF), the dot products of any two animals’ factor vectors reconstruct the interaction matrix (Wang et al., 2011; Shi et al., 2015), see Figure 1 a and b). If the embeddings allow us to predict relevant behavioral properties, they serve our understanding as semantic representations. However, in temporal settings where the interaction matrices change over time, there is no straightforward extension of this algorithm. The interaction matrices at different time points can be factorized individually, but there is no guarantee that the embeddings stay semantically consistent over time, i.e. the prediction of relevant behavioral properties will deteriorate.
In living systems, interaction dynamics are highly variable; individuals differ in when they appear in the data and how long they live. Different non-overlapping groups of individuals, e.g. from different years, may not interact with each other at all. How can we find a common semantic embedding even in these extreme cases? How do we learn embeddings that generalize to different groups and still provide insights into each individual’s functional role? If animals take on roles partially determined by a common factor, such as age, how can we learn this dependency?
Several approaches to extend NMF to temporal settings have been proposed in a variety of problem settings. Yu et al. (2016) and Mackevicius et al. (2019) propose a factorization method for time
series analysis. Gauvin et al. (2014) focus on the analysis of communities that are determined by their temporal activity patterns. Jiao et al. (2017) consider the case of communities from graphs over time and enforce temporal consistency with an additional loss term. Yu et al. (2017) and Wang et al. (2017) represent networks as a function of time. Temporal matrix factorization can be seen as a tensor decomposition problem, for which methods with many applications have been proposed in the literature, see Kolda & Bader (2009) for a review. In particular, time-shifted tensor decomposition methods have been used in multi-neuronal spike train analysis when recordings of multiple trials from a population of neurons are available (Mørup et al., 2008; Williams, 2020). However, to our knowledge, no method yet considers the case of entities following a common trajectory depending on an observable property (e.g., the age of an individual), which we show to be a powerful inductive bias.
We approach this question using honey bees, a popular model system for studying individual and collective behavior (Elekonich & Roberts, 2005). Bees allocate tasks across thousands of individuals without central control, using an age-based system: young bees care for brood, middle-aged bees perform within-nest labor, and old bees forage outside (Seeley, 1982; Johnson, 2010). Colonies are also organized spatially: brood is reared in the center, honey and pollen are stored at the periphery, and foragers offload nectar near the exit. Therefore, an individual’s role is partially reflected in its location, which allows us to evaluate whether the embeddings our method learns are semantically meaningful.
We propose jointly learning two meaningful representations of honey bee social behavior: 1) an individuality embedding that characterizes the lifetime behavior of each individual and 2) a daily representation of the individual’s functional position in the social network that can be derived from the individuality embeddings. We show that these representations can be learned in an unsupervised fashion, using only interaction matrices of the individuals over time. We analyze a dataset obtained by tracking thousands of individually marked honey bees in two colonies, at high temporal and spatial resolution over a total of 155 days, covering entire lifespans and multiple generations. We evaluate how well the embeddings capture the semantic differences of individual honey bee development by evaluating their predictiveness for different tasks and behaviorally relevant metrics.
Our method can be used to study dynamic groups in which individual units change their latent properties over time, and group composition is not fixed due to birth, death, or migration of individuals. We provide a new perspective on how to extract fundamental factors that underlie behavioral and developmental plasticity in animal groups.

heading: 2 METHODS

heading: 2.1 TEMPORAL NMF ALGORITHM
text: SymNMF factorizes a matrix A ∈ RN×N+ such that it can be approximated by the product FF T , where F ∈ RN×M+ and M N :
F̂ = argmin F≥0 ∥∥∥A− FF T∥∥∥2 Ai,j ≈ f(i) · f(j)T f(i) = Fi,: f(i) ∈ RM+ (1) When applied to social networks, f(i) can represent the role of an entity within the social network A (Wang et al., 2011; Shi et al., 2015) - however, in temporal settings, factorizing the matrices for different times separately will result in semantically inconsistent factors F .
Here we present a novel temporal NMF algorithm (TNMF) which extends SymNMF to temporal settings in which A ∈ RT×N×N+ changes over time t. We assume that the entities i ∈ {0, 1, . . . , N} follow to some extent a common trajectory depending on an observable property (for example the age of an individual). We represent an entity at a specific point in time t using a factor vector f+(t, i) such that
Ât,i,j = f+(t, i) · f+(t, j)T Â ∈ RT×N×N+ f+(t, i) ∈ RM+ (2)
In contrast to SymNMF, we do not directly factorize At to find the optimal factors that reconstruct the matrices. Instead, we decompose the problem into learning an average trajectory of factors m(c(t, i)) and structured variations from this trajectory o(t, i) that depend on the observable property c(t, i):
f(t, i) = m(c(t, i)) + o(t, i) f+(t, i) = max(0,f(t, i)) (3)
c : NT×N → N m : N→ RM+ o : NT×N → RM
This decomposition is an inductive bias that allows the model to learn semantically consistent factors for entities, even if they do not share any data points (e.g., there is no overlap in their interaction partners), as long as the relationship between functional role and c(t, i) is stable. Note that in the simplest case c(t, i) = t, TNMF can be seen as a tensor decomposition model, i.e. the trajectory of all entities is aligned with the temporal dimension t of A. In our case, c(t, i) maps to the age of individual i at time t.
While many parameterizations for the function o(t, i) are possible, we only consider one particular case in this work: We learn a set of individuality basis functions b(c(t, i)) (shared among all entities) that define a coordinate system of possible individual variations and the individuality embeddings φ, which capture to what extent each basis function applies to an entity:
o(t, i) = K∑ k=0 φi,k · bk(c(t, i)) φ : RN×K bk : NT → R (4)
where K is the number of learned basis functions. This parameterization allows us to disentangle the forms of individual variability (individuality basis functions) and the distribution of this variability (individuality embeddings) in the data.
We implement the functions m(c(t, i)) and b(c(t, i)) with small fully connected neural networks with non-linearities and several hidden layers. The parameters θ of these functions and the entities’ embeddings φ are learned jointly using minibatch stochastic gradient descent:
θ̂, φ̂ = argmin θ,φ ∥∥∥A− Â∥∥∥2 (5) Note that non-negativity is not strictly necessary, but we only consider the non-negative case in this work for consistency with prior work (Wang et al., 2011; Shi et al., 2015). Furthermore, instead of one common property with discrete time steps, the factors could depend on multiple continuous properties, i.e. c : RT×N → RP , e.g. the day and time in a intraday analysis of social networks.

heading: 2.2 REGULARIZATION
text: Due to the inductive bias, the model performs well without additional explicit regularization in the datasets considered in this work (see Table 2). However, we find that the model’s interpretability can be improved using additional regularization terms without significantly affecting its performance. We encourage sparsity in both the number of used factors and individuality basis functions by adding L1 penalties of the mean absolute magnitude of the factors f(t, i) and basis functions b(c(t, i)) to the objective. We encourage individuals’ lifetimes to be represented with a sparse embedding using an L1 penalty of the learned individuality embeddings φ.
We also introduce an optional adversarial loss term to encourage the model to learn embeddings that are semantically consistent over time, i.e. to only represent two entities that were present in the dataset at different times with different embeddings if this is strictly necessary to factorize the matrices A. We jointly train a discriminative network d(φi) that tries to classify the time of the first occurrence of all entities based on their individuality embeddings φ. The negative cross-entropy loss of this model is added as a regularization term to equation 5 in a training regime similar to generative adversarial networks (Goodfellow et al., 2014). See appendix A.2.1 for more details and A.4 for an ablation study of the effect the individual regularization terms have on the results of the model.

heading: 2.3 DATA
text: 2.3.1 HONEY BEE DATA
Two colonies of honey bees were continuously recorded over a total of 155 days. Each individual was manually tagged at emergence, so the date of birth is known for each bee. Locations and identities of all honey bees (N=9286) were extracted from the raw images and used to con-
struct daily aggregated temporal interaction networks based on the counts of spatial proximity events. The dataset also contains labels that can be used in proxy tasks (see section 2.5) to quantify if the learned embeddings and factors are semantically meaningful and temporally consistent. See table 1, and appendix A.1 for details. We publish the data together with this paper1.

heading: 2.3.2 SYNTHETIC DATA
text: Although some ground truth labels exist for the honey bee datasets (section 2.5), factorizing these data is still fundamentally an unsupervised learning problem. We generated synthetic datasets to evaluate whether the model can identify groups of individuals with common latent factors that determine their interaction frequencies. We model a common lifetime trajectory and groups of individual variation of factors using smoothed Gaussian random walks. We randomly assign individuals to a group with random times of emergence and disappearance from the data and compute interaction matrices by calculating the dot product of the factors of 1024 individuals for 100 simulated days. We then measure how well the individuality embeddings φ of a fitted model correspond to the truth groups using the adjusted mutual information score (Vinh et al., 2009), and the mean squared error between the ground truth factors and the best permutation of the factors f+. We evaluate the model on 128 different random synthetic datasets with increasing Gaussian noise levels in the interaction tensor. See appendix A.1.2 for more details on the generation process.
In both datasets, we define c(t, i) as the age in days of an individual i at time t.

heading: 2.4 BASELINE MODELS
text: Symmetric NMF: We compute the factors that optimally reconstruct the original interaction matrices using the standard symmetric NMF algorithm (Shi et al., 2015; Kuang et al., 2015), for each day separately, using the same number of factors as in the TNMF model.
Aligned symmetric NMF: We consider a simple extension of the standard SymNMF algorithm that aligns the factors to be more consistent over time. For each pair of subsequent days, we consider all combinatorial reorderings of the factors computed for the second day. For each reordering, we compute the mean L2 distance of all individuals that were alive on both days. We then select the reordering that minimizes those pairwise L2 distances and greedily continue with the next pair of days until all factors are aligned. Furthermore, we align the factors across colonies (where individuals cannot overlap) as follows: we run this algorithm for both datasets separately and align the resulting factors by first computing the mean embedding for all individuals grouped by their ages. As before, we now select from all combinatorial possibilities the reordering that minimizes the L2 distance between the embeddings obtained from both datasets. See section A.5.1 for pseudo code.
Tensor decomposition: We also compare against a constrained non-negative tensor decomposition model with symmetric factors F ∈ RN×M+ and temporal dynamics constrained to the diagonals, i.e. D ∈ RT×M×M+ and Dt = diag(dt), dt ∈ RM+ .
Ât = FDtF T (6)
F̂ , D̂ = argmin F ,D T−1 T∑ t=0 ∥∥∥At − Ât∥∥∥2 (7) Temporal NMF models: We evaluate variants of the temporal symmetric matrix factorization algorithms proposed by Jiao et al. (2017) and Yu et al. (2017).
1See Anonymous (2020) for dataset
For the tensor decomposition and temporal NMF baselines, we follow the procedure given above for the Aligned symmetric NMF to find the optimal reordering to align the factors obtained by applying models to the two datasets separately.

heading: 2.5 EVALUATION METRICS
text: Reconstruction: We measure how well the original interaction matrices A can be reconstructed from the factors. We do not require the model to reconstruct the interaction matrices as well as possible because we only use the reconstruction as a proxy objective to learn a meaningful representation. Still, a high reconstruction loss could indicate problems with the model, such as excessive regularization.
Consistency: We measure to what extent the individuality embeddings φ change over time. For each model, we train a multinomial logistic regression model to predict the source cohort (date of birth) and calculate the area under the ROC curve (AUCcohort) using a stratified 100-fold cross-validation with scikit-learn (Pedregosa et al., 2011). The baseline models do not learn an individuality embedding; therefore we compute how well the model can predict the cohort using the mean factor representation of the individuals over their lives. We define consistency as 1− AUCcohort of this linear model. Note that a very low temporal consistency would indicate that the development of individual bees changes strongly between cohorts and colonies, which we know not to be true.
Mortality and Rhythmicity: We evaluate how well a linear regression model can predict the mortality (number of days until death) and circadian rhythmicity of the movement (R2 score of a sine with a period of 24 h fitted to the velocity over a three-day window). These metrics are strongly correlated with an individual’s behavior (e.g. foragers exhibit strong circadian rhythms because they can only forage during the daytime; foragers also have a high mortality). We follow the procedure given in Wild et al. (2020) and report the 100-fold cross-validated R2 scores for these regression tasks.
For this data, we expect the factors f+ and individuality embeddings φ to be semantically meaningful and temporally consistent if they reflect an individual’s behavioral metrics (mortality and rhythmicity) and if they do not change strongly over time (measured in the consistency metric).

heading: 3 RESULTS
text: We implemented the model using PyTorch (Paszke et al., 2019) and trained it in minibatches of 256 individuals for 200 000 iterations with the Adam optimizer (Kingma & Ba, 2015). See appendix A.2.3 for the architecture of the learned functions, a precise description of the regularization losses, and further hyperparameters. The code is publicly available 2.

heading: 3.1 SYNTHETIC DATA
text: We factorize the interaction matrices of the 128 synthetic datasets with varying levels of Gaussian noise. We confirmed that our model converges in all datasets and evaluate whether we can distinguish
2https://anonymous.4open.science/r/b2b7e2fc-aa04-4cf8-85c7-646d8dc46400
the individuals’ ground truth group assignments. To that end, we extract the individuality embeddings φ from the models and measure how well they correspond to ground truth data using the adjusted mutual information (AMI) score. Furthermore, we measure the mean squared error between the best permutation of learned factors f+ and the ground truth factors. We find that for low levels of noise, our model can identify the truth group assignments with high accuracy, and are still significantly better than random assignments even at very high levels of noise (see figure 3). Note that for this experiment, we evaluated a model with the same hyperparameters as used in all plots in the results section (see Table 2) and a variant without explicit regularization except the L1 penalty of the learned individuality embeddings φ (λembeddings, because this regularization is required to meaningfully extract clusters), which was set to 0.1. See appendix A.1.2 for more details on the synthetic datasets.

heading: 3.2 HONEY BEES
text: Mean lifetime model: The model learns a sparse representation of the developmental trajectory of a honey bee in the space of social interactions. Only two factors are effectively used (they exceed the threshold value of 0.01). These factors show a clear trend over the life of a bee, indicating that the model captures the temporal aspects of the honey bee division of labor (See Figure 4).
Interpretability of factors: To understand the relationship between the factors and division of labor, we calculate how the factors map to the fraction of time an individual spent on the brood area, honey storage, or dance floor (where foragers aggregate). Time spent on these different substrates is a strong indicator of an individual’s task. The factor f1, which peaks at young age (Figure 4), correlates with the proportion of time spent in the brood area, while a high f0 indicates increased time spent on the dance floor. Therefore, the model learned to map biologically relevant processes.
Individuality basis functions and individuality embeddings: Due to the regularization of the embeddings, the model learns a sparse set of individuality basis functions. As encouraged by the model, most individuals can predominantly be described by a single basis function. That means that while each honey bee can collect a unique set of experiences, most can be described with a few common individuality embeddings which are consistent across cohorts and colonies. In the context of honey bee division of labor, the basis functions are interpretable because the factors correspond to different task groups. For example, b12(c(t, i)) (accounting for≈ 10.7% of the individuals) describes workers that occupy nursing tasks much longer than most bees. As the individuality embeddings φ only scale the magnitude of the basis functions, they can be interpreted in the same way. Individual lifetime trajectories in the factor space can be computed based on the mean lifetime trajectories (m), individuality basis functions (b(c(t, i))) and individuality embeddings (φ). See figure 5 for examples
of individual lifetime trajectories from workers that most strongly corresponded to the common individuality basis functions.
Evaluation: We verify that the learned representations of the individuals are meaningful (i.e., they relate to other properties of the individuals, not just their interaction matrices) and semantically consistent over time and across datasets using the metrics described in the section Evaluation metrics. We compare variants of our model with different adversarial loss scaling factors and factor L1 regularizations, the baseline models, and the individuals’ ages. We expect a good model to be temporally consistent and semantically meaningful. All variants of our model outperform the baselines in terms of the semantic metrics Mortality and Rhythmicity, except for the Yu et al. (2017) model, which performs comparably well in the Mortality metric. The adversarial loss term further increases the Consistency metric without negatively affecting the other metrics. A very strong adversarial regularization (see row with λadv = 1 in Table 2) prevents the model from learning a good representation of the data. See Table 2 for an overview of the results. We also evaluate the tradeoff between the different metrics using a grid search over the hyperparameters (see appendix A.3).
Scalability: The functions m(c(t, i)) and b(c(t, i)) are learned neural networks with non-linearities. The objective is non-convex and we learn the model parameters using stochastic gradient descent. Optimization is therefore slower than the standard NMF algorithms that can be fitted using algorithms such as Alternating Least Squares (Kim et al., 2014). We found that the model converges faster if the reconstruction loss of the age based model m(c(t, i)) is additionally minimized with the main objective in equation 5. Due to the minibatch training regime, our method should scale well in larger datasets. Small neural networks were sufficient to learn the functions m(c(t, i)) and b(c(t, i)) in our experiments. Most of the runtime during training is spent on the matrix multiplication f+(t, i) · f+(t, j)T and the corresponding backwards pass.

heading: 4 CONCLUSION
text: Temporal NMF factorizes temporal matrices with overlapping and even disjointed communities by learning a common embedding of the lifetime development of the individual entities. In the context of honey bees, this embedding is biologically meaningful, consistent over time, and shows that interaction patterns follow a common lifetime trajectory. Differences from the mean are described in a coordinate system of individual variability. Honey bee colonies are known to exhibit suites of correlated traits (Wray et al., 2011), but this method provides a perspective at the individual level and opens the possibility for understanding how colony-level traits respond to biotic and abiotic pressures. The basis functions are interpretable with respect to the division of labor within colonies and offer a valuable tool to understand and quantify the influence of experimental manipulations (e.g. pesticides or increased predatory pressure) on an individual, even between experimental trials on completely different sets of individuals. While we applied our method to honey bees as an exemplary system
with many individuals that exhibit an entangled, non-overlapping social structure, our method can be applied to any setting in which some interaction structure follows a general pattern over an observable (such as time) to detect structured deviations at the individual level.

Review summaries:
The authors introduce a novel method for non-negative matrix factorization for timeseries and apply it to longitudinal honey bee interaction data.  The model leverages consistency of individuals over time by forcing the factors (or rather, the residuals of the factors with respect to a global trajectory) to be linear combinations of a small set of temporal basis functions.  These temporal basis functions are functions of the bee’s age.  In other words, the factor embedding of each bee is a vector of linear combinations of 16 learned basis functions of time.  Since all bees use the same 16 temporal basis functions, given these basis functions the lifetime embedding of each bee is encapsulated by a small matrix of numbers, namely the coefficients for the temporal basis functions for each factor (and in practice only two factors were significant, so each bee’s life is embedded in 32-dimensional space).  There are a number of regularizations on the temporal basis functions and the embedding coefficients.

The authors present a matrix factorization model to jointly characterize the lifetime interactions of thousands of bees over generations. The problem is fascinating as both a technical and scientific question and the modeling framework appears novel. Although not directly addressed, the authors appear to be trying to solve a *tensor* factorization problem, not just the special case of a matrix (which is of course a 2-d tensor). It would have been interesting to see results comparing their method with a non-negative variant of, say, PARAFAC/CANDECOMP or generalizations thereof. It would have at least have been appropriate to explain why or why not existing tensor methods are not appropriate.

This paper proposes a NMF formulation ||A-FF^T||^2 where A and F are different types of information extracted from social datasets. In the honeybee example the authors highlight, A represents the spatial relationship between bees, and F encodes the age of the bees. The authors setup F to be decomposable into two types of embeddings, one which characterizes the group activity and the other which characterizes the individual activity.