Learning interaction rules from multi-animal trajectories via augmented behavioral models

Abstract:
Extracting the interaction rules of biological agents from movement sequences pose challenges in various domains. Granger causality is a practical framework for analyzing the interactions from observed time-series data; however, this framework ignores the structures and assumptions of the generative process in animal behaviors, which may lead to interpretational problems and sometimes erroneous assessments of causality. In this paper, we propose a new framework for learning Granger causality from multi-animal trajectories via augmented theory-based behavioral models with interpretable data-driven models. We adopt an approach for augmenting incomplete multi-agent behavioral models described by time-varying dynamical systems with neural networks. For efficient and interpretable learning, our model leverages theory-based architectures separating navigation and motion processes, and the theory-guided regularization for reliable behavioral modeling. This can provide interpretable signs of Granger-causal effects over time, i.e., when specific others cause the approach or separation. In experiments using synthetic datasets, our method achieved better performance than various baselines. We then analyzed multi-animal datasets of mice, flies, birds, and bats, which verified our method and obtained novel biological insights.

heading: 1 Introduction
text: Extracting the interaction rules of real-world agents from data is a fundamental problem in a variety of scientific and engineering fields. For example, animals, vehicles, and pedestrians observe other’s states and execute their actions in complex situations. Discovering the directed interaction rules of such agents from observed data will contribute to the understanding of the principles of biological
∗fujii@i.nagoya-u.ac.jp
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
agents’ behaviors. Among methods analyzing directed interactions within multivariate time series, Granger causality (GC) [23] is a practical framework for exploratory analysis [49] in various fields, such as neuroscience [65] and economics [3] (see Section 2). Recent methodological developments have focused on inferring GC under nonlinear dynamics (e.g., [73, 31, 81, 55, 43]).
However, the structure of the generative process in biological multi-agent trajectories, which include navigational and motion processes [54] regarded as time-varying dynamical systems (see Section 3.1), is not fully utilized in existing base models of GC including vector autoregressive [27] and recent neural models [73, 31, 81]. Ignoring the structures of such processes in animal behaviors will lead to interpretational problems and sometimes erroneous assessments of causality. That is, incorporating the structures into the base model for inferring GC, e.g., augmenting (inherently) incomplete behavioral models with interpretable data-driven models (see Section 3.2), can solve these problems. Furthermore, since data-driven models sometimes detect false causality that is counterintuitive to the user of the analysis, e.g., introducing architectures and regularization to utilize scientific knowledge (see Sections 3.2 and 4.2) will be effective for a reliable base model of a GC method.
In this paper, we propose a framework for learning GC from biological multi-agent trajectories via augmented behavioral models (ABM) using interpretable data-driven neural models. We adopt an approach for augmenting incomplete multi-agent behavioral models described by time-varying dynamical systems with neural networks (see Section 3.2). The ABM leverages theory-based architectures separating navigation and motion processes based on a well-known conceptual behavioral model [54], and the theory-guided regularization (see Section 4.2) for interpretable and reliable behavioral modeling. This framework can provide interpretable signs of Granger-causal effects over time, e.g., when specific others cause the approach or separation.
The main contributions of this paper are as follows. (1) We propose a framework for learning Granger causality via ABM, which can extract interaction rules from real-world multi-agent and multi-dimensional trajectory data. (2) Methodologically, we realized the theory-guided regularization for reliable biological behavioral modeling for the first time. The theory-guided regularization can leverage scientific knowledge such that “when this situation occurs, it would be like this” (i.e., domain experts often know an input and output pair of the prediction model). Existing methods in Granger causality did not consider the utilization of such knowledge. (3) Biologically, our methodological contributions lies in the reformulation of a well-known conceptual behavioral model [54], which did not have a numerically computable form, such that we can compute and quantitatively evaluate it. (4) In the experiments, our method achieved better performance than various baselines using synthetic datasets, and obtained new biological insights and verified our method using multiple datasets of mice, birds, bats, and flies. In the remainder of this paper, we describe the background of GC in Section 2. Next, we formulate our ABM in Section 3, and the learning and inference methods in Section 4.

heading: 2 Granger Causality
text: GC [23] is one of the most popular and practical approaches to infer directed causal relations from observational multivariate time series data. Although the classical GC is defined by linear models, here we introduce a more recent definition of [73] for non-linear GC. Consider p stationary time-series x = {x1, ...xp} across timesteps t = {1, ..., T} and a non-linear autoregressive function gj , such that
xjt+1 = gj(x 1 ≤t, ...,x p ≤t) + ε j t , (1)
where xj≤t = (...,x j t−1,x j t ) denotes the present and past of series j and ε j t represents independent noise. We then consider that variable xi does not Granger-cause variable xj , denoted as xi 9 xj , if and only if gj(·) is constant in xi≤t. Granger causal relations are equivalent to causal relations in the underlying directed acyclic graph if all relevant variables are observed and no instantaneous (i.e., connections between two variables at the same timestep) connections exist [59]. Many methods for Granger causal discovery, including vector autoregressive [27] and recent deep learning-based approaches [73, 31, 81], can be encapsulated by the following framework. First, we define a function fθ (e.g., an multilayer perceptrons (MLP) in [73], a linear model in [27]), which learns to predict the next time-step of the test sequence x. Then, we fit fθ to x by minimizing some loss (e.g., mean squared error) L: θ? = argminθ L(x, fθ). Finally, we apply some fixed function h (e.g., thresholding) (e.g., [45]) to the learned parameters to produce a Granger causal graph estimate for x: Ĝx = h(θ?).
Furthermore, we need to differentiate between positive and negative Granger-causal effects (e.g, approaching and separating). Based on the definition of [45], we define the effect sign as follows: if gj(·) is increasing in all xi≤t, then we say that variable xi has a positive effect on xj , if gj(·) is decreasing in xi≤t, then x
i has a negative effect on xj . Note that xi can contribute both positively and negatively to the future of xj at different delays.
Overall, the causality measures, however elaborate in construction, are simply statistics estimated from a model [71]. If the model inadequately represents the system properties of interest, subsequent analyses based on the model will fail to address the question of interest. The inability of the model to represent key features of interest can cause interpretational problems and sometimes erroneous assessments of causality. Therefore, in our case, incorporating the structures of the generative process for animal behaviors (i.e., Eq.(2)) in a numerically computable form will be required. We thus propose the ABM based on a well-known conceptual model [54] in biological sciences in the next section.

heading: 3 Augmented behavioral model
text: Our motivation for developing interpretable behavior models is to obtain new insights from the results of Granger causality. In this section, we firstly formulate a well-known conceptual behavioral model [54] so that it can be computable. Second, we propose (multi-animal) ABMs with theorybased architectures based on scientific knowledge. Further, we discuss the relation to the existing explainable neural models [2]. The diagram of our method is described in Appendix C.

heading: 3.1 Formulation of a conceptual behavioral model
text: In movement ecology, which is a branch of biology concerning the spatial and temporal patterns of behaviors of organisms, a coherent framework [54] has been conceptualized to explore the causes, mechanisms, and patterns of movement. For example, two alternative structural representations [54] were proposed to model a new position pt+1 from its current location pt (for details, see Appendix A): the motion-driven case pt+1 = fU (fM (Ω, fN (Φ, rt,wt,pt), rt,wt,pt)) + εt, and the navigation-driven case pt+1 = fU (fN (Φ, fM (Ω, rt,wt,pt), rt,wt,pt)) + εt, where wt is the internal state, Ω is the motion capacity, Φ is the navigation capacity, and rt is the environmental factors (these are conceptual parameters). fM , fN , and fU are conceptual functions to represent actions of the motion (or planning), navigation, and movement progression processes, respectively.
For efficient learning of the weights in the model (i.e., coefficient of Granger causality) in this paper, we consider a simple case with homogeneous navigation and motion capacities, and internal states. Moreover, to make the contribution of fM , fN , and fU interpretable after training from the data for extracting unknown interaction rules (and for obtaining scientific new insights), one of the simplified processes for agent i is represented by
xit+1 = f i U (f i N (r i t,x i t), f i M (r i t,x i t), r i t,x i t) + ε i t, (2)
where xi ∈ Rd includes location pi and velocity for the agent i. We here consider ri ∈ R(p−1)dr including p− 1 other agents’ dr-dimensional information. This formulation does not assume either motion-driven or navigation-driven case. Such behaviors have been conventionally modeled by mathematical equations such as force- and rule-based models (e.g., reviewed by [77, 47]). Recently, these models have become more sophisticated by incorporating the models into hand-crafted functions representing anticipation (e.g., [29, 51]) and navigation (e.g., [8, 76]).
However, these conventional and recent models are sometimes too simplistic and customized for the specific animals, respectively; thus it is sometimes difficult to define the dynamics of general biological multi-agent systems (i.e., multiple species of animals). Therefore, methods for learning parameters and interaction rules of behavioral models are needed. There have been some researches to estimate specific parameters (and their distributions) of the interpretable behavior models (e.g., [84, 85, 15]), and others to model the parameters and rules in purely data-driven manners (i.e., sometimes uninterpretable) only for accurate prediction (e.g., [16, 28]). In the proposed framework, we consider flexible data-driven interpretable models to focus on inferring GC for exploratory analysis from the observed data without specific knowledge of the species and obtaining additional data.
Recently, some attempts have been made to explore flexible and interpretable models bridging theorybased and data-driven approaches. For example, a paradigm called theory-guided data science has
been proposed [30], which leverages the wealth of scientific knowledge for improving the effectiveness of data-driven models in enabling scientific discovery. For example, scientific knowledge can be used as architectures or regularization terms in learning algorithms in physical and biological sciences (e.g., [61, 21]). In biological multi-agent systems, an approach extract interpretable dynamical information based on physics-based knowledge [18] from multi-agent interaction data, and another approach made a particular module such as observation (e.g., [19]) interpretable in mostly black-box neural models. However, these data-driven models did not sufficiently utilize the above scientific knowledge of multi-animal interactions. In the next subsection, to make the model (e.g., of GC) flexible and interpretable, we propose an ABM with theory-based architectures.

heading: 3.2 Augmented behavioral model with theory-based architectures
text: In this subsection, we propose a ABM using interpretable neural models with theory-based architectures for learning GC from multi-animal trajectories. In general, it is scientifically beneficial if a model mimics the data-generating process well, e.g., because existing scientific insights can be leveraged or revalidated. In our case of GC, additionally, it is expected to eliminate obvious erroneous causality by utilizing existing knowledge, we thus propose a theory-based ABM for learning GC.
Generally, scientific knowledge can be used to influence the architecture of data-driven scientific models. Most design considerations are mainly motivated to simplify the learning procedure, minimize the training loss, and ensure robust generalization performance [30]. In some cases, domain knowledge can be used designing neural models by decomposing the overall problem into modular sub-problems. For example, in our problem, to describe the overall process of multi-animal behaviors, modular neural models can be learned for different sub-processes, such as the navigation, planning, and movement processes (f iN , f i M , and f i U , respectively) described in Section 3.1. This will help in using the power of learning frameworks while following a high-level organization in the architecture that is motivated by domain knowledge [30]. Specifically, to accurately model the relationships between agents (finally interpreted as causal relationships) with limited information in usual GC settings, we explicitly formulate the f iN , f i M , and f i U and estimate f i N and f i M from data.
In summary, our base ABM can be expressed as
xit = K∑ k=1 ( F i,t,kN (h i t−k) F i,t,k M (h i t−k) ) hit−k + ε i t, (3)
where hit−k ∈ Rdh is a vector concatenating the self state xit−k ∈ Rd and all others’ state rit−k ∈ R(p−1)dr , and denotes a element-wise multiplication. K is the order of the autoregressive model. F i,t,kN , F i,t,k M : Rdh → Rd×dh are matrix-valued functions that represent navigation and motion functions, which are implemented by MLPs. For brevity, we omit the intercept term here and in the following equations. The value of the element of F i,kN is [−1, 1] is like a switching function value, i.e., a positive or negative sign to represent the approach and separation from others. The value of the element of F i,kM is a positive value or zero, which changes continuously and represents coefficients of time-varying dynamics. Relationships between agents x1, ...,xp and their variability throughout time can be examined by inspecting coefficient matrices Ψiθt,k = ( F i,t,kN (h i t−k) F i,t,k M (h i t−k) ) . We separate Ψiθt,k into F i,t,k N (h i t−k) and F i,t,k M (h i t−k) for two reasons: interpretability and efficient use of scientific knowledge. The interpretability of two coefficients F i,kN and F i,k M contributes to the understanding of navigation and motion planning processes of animals (i.e., signs and amplitudes in the GC effects), respectively. The efficient use of scientific knowledge in the learning of a model enables us to incorporate the knowledge into the model. The effectiveness was shown in the ablation studies in the experiments. Specific forms of Eq. (3) are described in Appendices E.2 and G.2. The formulation of the model via linear combinations of the interpretable feature hit−k for an explainable neural model is related to the self-explanatory neural network (SENN) [2].

heading: 3.3 Relation to self-explanatory neural network
text: SENN [2] was introduced as a class of intrinsically interpretable models motivated by explicitness, faithfulness, and stability properties. A SENN with a link function g(·) and interpretable basis concepts h(x) : Rp → Ru follows the form
f(x) = g(θ(x)1h(x)1, ..., θ(x)uh(x)u), (4)
where x ∈ Rp are predictors; and θ(·) is a neural network with u outputs (here, we consider the simple case of d = 1 and dr = 1). We refer to θ(x) as generalized coefficients for data point x and use them to explain contributions of individual basis concepts to predictions. In the case of g(·) being sum and concepts being raw inputs, Eq. (4) simplifies to f(x) = ∑p i=1 θ(x)ixi. In this paper, we regard the movement function f iU as g(·) and the function of f iN and f iM as θ for the following interpretable modeling of f iU , f i N , and f i M . Appendix B presents additional properties SENNs need to satisfy and the learning algorithm, as defined by [2]. Note that our model does not always satisfy the requirements of SENN [2, 45] due to the modeling of time-varying dynamics (see Appendix B). SENN was first applied to GC [45] via generalized vector autoregression model (GVAR): xt = ∑K k=1 Ψθk(xt−k)xt−k + εt, where Ψθk : Rp → Rp×p is a neural network parameterized by θk. Ψθk(xt−k) is a matrix whose components correspond to the generalized coefficients for lag k at timestep t. The component (i, j) of Ψθk(xt−k) corresponds to the influence of x j t−k on x i t. However, the SENN model did not use scientific knowledge of multi-element interactions and may cause interpretational problems and sometimes erroneous assessments of causality.

heading: 4 Learning with theory-guided regularization and inference
text: Here, we describe the learning method of the ABM including theory-guided regularization. We first overview the learning method and define the objective function. We then explain the theory-guided regularization for incorporating scientific knowledge into the learning of the model. Finally, we describe the inference of GC by our method. Again, the overview of our method is described in Appendix C.

heading: 4.1 Overview
text: To mitigate the inference in multivariate time series, Eq. (3) for each agent is summarized as the following expression:
xt = K∑ k=1 [( F 1,t,kN (h 1 t−k) F 1,t,kM (h 1 t−k) ) h1t−k, . . . , ( F p,t,kN (h p t−k) F p,t,k M (h p t−k) ) hpt−k ] + εt, (5)
where xt and εt concatenate the original variables for all p agents (various F s and Psis are learned parameters). We train our model by minimizing the following penalized loss function with the mini-batch gradient descent
T∑ t=K+1 ( Lpred(x̂t,xt) + λLsparsity(Ψt) + γLTG(Ψt,ΨTGt ) ) + T−1∑ t=K+1 βLsmooth(Ψt+1,Ψt), (6)
where {xt}Tt=1 is a single observed time series of length T with d-dimensions and p-agents; x̂t is the one-step forecast for the t-th time point based on Eq. (5); Ψt ∈ Rpd×Kdh is defined as a concatenated matrix of [Ψ1θt,K , . . . ,Ψ 1 θt,1 ], . . . , [Ψpθt,K , . . . ,Ψ p θt,1
] in a row; ΨTGt is a coefficient determined by the following theory-guided regularization; and λ, β, γ ≥ 0 are regularization parameters. The loss function in Eq. (6) consists of four terms: (i) the mean squared error (MSE) prediction loss, (ii) a sparsity-inducing penalty term, (iii) theory-guided regularization, and (iv) the smoothing penalty term. The sparsity-inducing term Lsparsity is an appropriate penalty on the norm of Ψt. Among possible various regularization terms, in our implementation, we employ the elastic-net-style penalty term [86, 56] Lsparsity(Ψt) = 1T−K ( α ‖Ψt‖1 + (1− α) ‖Ψt‖ 2 F ) , with α = 0.5, based on [45]. Note that other penalties can be also easily adapted to our model. The smoothing penalty term, given by Lsmooth(Ψt+1,Ψt) = 1T−K−1 ‖Ψt+1 −Ψt‖ 2 F , is the average norm of the difference between generalized coefficient matrices for two consecutive time points. This penalty term encourages smoothness in the evolution of coefficients with respect to time [45]. To avoid overfitting and model selection problems, we eliminate unused factors based on prior knowledge (for details, see Appendices E.2 and G.2).

heading: 4.2 Theory-guided regularization
text: The third term in Eq. (6) is the theory-guided regularization for reliable Granger causal discovery by leveraging regularization with scientific knowledge. Here we utilize theory-based and data-driven prediction results and impose penalties in the appropriate situations as described below. Again, let x̂t be the prediction from the data. In addition to the data, we prepare some input-output
pairs (x̃t−k≤t, x̃t) based on scientific knowledge. We call them pairs of theory-guided feature and prediction, respectively. In this case, we assume that the theory-guided cause or weight of the ABM ΨTGt is uniquely determined. When the difference between x̂t and x̃t is below a certain threshold, we assume that the cause (weight) of x̂t is equivalent to the cause of x̃t.
In animal behaviors, the theory-guided prediction utilizes the intuitive prior knowledge such that the agents go straight from the current state if there is no interaction. In this case, x̃t includes the same velocity as the previous step and the corresponding positions after going straight. The penalty is expressed as LTG(Ψt,ΨTGt ) = 1T−K exp(‖xt−x̃t‖ 2 2/σ)‖Ψ′t‖2F , where Ψ′t ∈ Rpd×K(p−1)dr is the weight matrix regarding others’ information (i.e., eliminating the information of the agents themselves from Ψt) and σ is a parameter regarding the threshold. Note that here the matrix Ψ ′TG t corresponding to ΨTGt is a zero matrix representing no interaction with others (i.e., ‖Ψ′t −Ψ ′TG t ‖2F = ‖Ψ′t‖2F ).
Next, we can consider the general cases. All possible combinations of the pairs are denoted as the direct product H0 := L × M × · · · × M = {(l,m1, . . . ,mp) | l ∈ L ∧ m1 ∈ M ∧ · · · ∧ mp ∈ M}, where L = {1, . . . , p} and M = {−1, 0, 1} if we consider the sign of Granger causal effects (otherwise, M = {0, 1}). However, if we consider the pairs (x̃t−k≤t, x̃t) uniquely determined, it will be a considerably fewer number of combinations by avoiding underdetermined problems. We denote the set of the uniquely-determined combinations as H1 ⊂ H0. We can then impose penalties on the weights: LTG(Ψt,ΨTGt ) =
1 |H1|(T−K) ∑ l,m1,...,mp∈H1 ( exp(‖xt − x̃t‖2F /σ)‖Ψ′t −Ψ ′TG t,l,m1,...,mp ‖2F ) , where Ψ ′TG t,l,m1,...,mp
∈ Rpd×K(p−1)dr is the weight matrix regarding others’ information in Ψt. In animal behaviors, due to unknown terms, such as inertia and other biological factors, the theory-guided prediction utilizes the only intuitive prior knowledge such that the agents go straight from the current state if there are no interactions (i.e., |H1| = 1).

heading: 4.3 Inference of Granger causality
text: Once Ψt is trained, we quantify strengths of Granger-causal relationships between variables by aggregating matrices Ψt across all K, d, dr, t into summary statistics. Although most neural GC methods [73, 55, 31, 81, 45] did not provide an obvious way for handling multi-dimensional time series (i.e., d > 1), our main problems include two- or three-dimensional positional and velocity data for each animal. Therefore, we compute the norm with respect to spatial dimensions d, dr, and the sign of the GC separately. That is, we aggregate the obtained generalized coefficients into matrix S ∈ Rp×p as follows:
Si,j = signmax K+1≤t≤T signmax1≤k≤K median
q=1,...,dr u=1,...,d
(Ψi,j)  maxK+1≤t≤T ( max 1≤k≤K ( ‖ (Ψi,j)t,k ‖F )) , (7)
where Ψi,j ∈ R(T−K)×K×d×dr is computed by reshaping and concatenating Ψt over K + 1 ≤ t ≤ T . ‖ (Ψi,j)t,k ‖F is the Frobenius norm of the matrix (Ψi,j)t,k ∈ R
d×dr for each t, k. The signmax is an original function to output the sign of the larger value of the absolute value of the maximum and minimum values (e.g., signmax({1, 2,−3}) = −1). If we do not consider the sign of Granger causal effects, we ignore the coefficient of the signed function. If we investigate the GC effects over time, we eliminate max function among t. Note that we only consider off-diagonal elements of adjacency matrices and ignore self-causal relationships. Intuitively, Si,j are statistics that quantify the strength of the Granger-causal effect of xi on xj using magnitudes of generalized coefficients. We expect Si,j to be close to 0 for non-causal relationships and Si,j 0 if xi → xj . Note that in practice Si,j is not binary-valued, as opposed to the ground truth, which we want to infer, because the outputs of Ψi,j are not shrunk to exact zeros. Therefore, we need a procedure deciding for which variable pairs Si,j are significantly different from 0. To infer a binary matrix of GC relationships, we use a heuristic threshold. For the detail, see Appendix D.1.

heading: 5 Related work
text: Methods for nonlinear GC. Initial work for nonlinear GC methods focused on time-varying dynamic Bayesian networks [70], regularized logistic regression with time-varying coefficients [34], and kernelbased regression models [46, 69, 39]. Recent approaches to inferring Granger-causal relationships leverage the expressive power of neural networks [50, 79, 73, 55, 31, 43, 81] and are often based on
regularized autoregressive models. Methods using sparse-input MLPs and long short-term memory to model nonlinear autoregressive relationships have been proposed [73], followed by a more sample efficient economy statistical recurrent unit (eSRU) architecture [31]. Other researchers proposed a temporal causal discovery framework that leverages attention-based convolutional neural networks [55] and a framework to interpret signs of GC effects and their variability through time building on SENN [2]. However, the structure of time-varying dynamical systems in multi-animal trajectories was not fully utilized in the above models.
Information-theoretic analysis for multi-animal motions. In this topic, most researchers have adopted transfer entropy (TE) and its variants and have analyzed them in terms of e.g., information cascades rather than causal discovery among animals. In the pioneering work, [78] analyzed information cascades among artificial collective motions using (conditional) TE [40, 41]. [64] applied variants of conditional mutual information to identify dynamical coupling between the trajectories of foraging meerkats. TE has been used to study the response of schools of zebrafish to a robotic replica of the animal [9, 36], to infer leadership in pairs of bats [58] and simulated zebrafish [10], and to identify interactions in a swarm of insects (Chironomus riparius) [42]. Local TE (or pointwise TE) [67, 40] has been used to detect local dependencies at specific time points in a swarm of soldier crabs [75], teams of simulated RoboCup agents [11], and a school of fish [13]. Since biological collective motions are intrinsically time-varying dynamical systems, we compared our methods and local TE in our experiments.
Other Biological multi-agent motion analysis. Previous studies have investigated leader-follower relationships. For example, the existences of the leadership have been investigated via the correlation in movement with time delay (e.g., [53, 66]) and via global physical (e.g., [4]) and statistical properties [57]. Meanwhile, methods for data-driven biological multi-agent motion modeling have been intensively investigated for pedestrian (e.g., [1, 25]), vehicles (e.g., [5, 63, 72]), animals [16, 28], and athletes (e.g., [83, 37]). In most of these methods, the agents are assumed to have the full observation of other agents to achieve accurate prediction. In contrast, some researches have modeled partial observation in real-world multi-agent systems [26, 38, 24, 17, 18, 19, 22]. However, the above approaches required a large amount of training data and would not be suitable for application to the multi-animal trajectory datasets that are measured in small quantities.

heading: 6 Experiments
text: The purpose of our experiments is to validate the proposed methods for application to real-world multianimal movement trajectories, which have usually a smaller amount of sequences and no ground truth of the interaction rules. Thus, for verification of our methods, we first compared their performances to infer the Granger causality to those in various baselines using two synthetic datasets with ground truth: nonlinear oscillator (Kuramoto model) and boid model simulation datasets. We used the same ABM as applied to real-world multi-animal trajectory datasets: mice, birds, bats, and flies. To demonstrate the applicability to the multi-element dynamics other than multi-animal trajectories, we validated our method using the Kuramoto dataset (the results are shown in Appendix F). Each method was trained only on one sequence according to most neural GC frameworks [73, 31, 45]. The hyperparameters of the models were determined by validation datasets in the synthetic data experiments (for the details, see Appendices E and G). The common training details, (binary) inference methods, computational resources, and the amount of computation are described in Appendix D. Our code is available at https://github.com/keisuke198619/ABM.

heading: 6.1 Synthetic datasets
text: For verification of our method, we compared the performances to infer the GC to those in various baselines using two synthetic datasets with ground truth. To compare with various baselines of GC methods, we tackled problems where the true causality is not changed over time. Here, we compared our methods (ABM) to 5 baseline methods: economic statistical recurrent unit (eSRU) [31]; amortized causal discovery (ACD) [43]; GVAR [45] (this is the most appropriate baseline); and simple baselines such as linear GC and local TE modified from [81, 43]. Except for ACD [43], most baselines did not provide an obvious way for handling multi-dimensional time series, whereas our main problems include two- or three-dimensional trajectories for each animal. Therefore, we modified the baselines except for ACD to compute norms with respect to spatial dimensions (2 or 3) for comparability with the proposed method. Note that the interpretations of the relationships estimated by ACD and Local TE are different from other methods, thus the sign of the relationship could not be investigated (we denote such by N/A in Table 1).
We investigated the continuously-valued inference results: the norms of relevant weights, scores, and strengths of GC relationships. We compared these scores against the true structures using areas under receiver operating characteristic (AUROC) and precision-recall (AUPRC) curves. We also evaluated thresholded inference results: accuracy (Acc) and balanced accuracy (BA) scores. For the inference methods of a binary matrix of GC relationships, see Appendix D.1. For all evaluation metrics, we only considered off-diagonal elements of adjacency matrices, ignoring self-causal relationships.
Boid model. Here, we evaluated the interpretability and validity of our method on the simulation data using the boid model, which contains five agents movement trajectories (for details, see Appendix G). In this experiment, we set the boids (agents) directed preferences: we randomly set the ground truth relationships 1, 0, and −1 as the rules of attraction, no interaction, and repulsion, respectively. Figure 1 illustrates that e.g., boid #5 was
attracted to boid #1 (i.e., true relationship: 1) and boid #3 avoided boid #1 (true relationship: −1). In this figure, our method detected the changes in the signed relationships whereas the GVAR [45] did not.
The performances were evaluated using Si,j in Eq. (7) throughout time because our method and ground truth were sensitive the sign as shown in Figure 1. Table 1 (upper) shows that our method achieved better performance than various baselines. The ablation studies shown in Table 1 (lower) reveal that the main two contributions of this work, the theory-guided regularization LTG and learning navigation function F kN and motion function F kM separately, improved the performance greatly. These suggest that the utilization of scientific knowledge via the regularization and architectures efficiently worked in the limited data situations. Similarly, the results of the Kuramoto dataset are shown in Appendix F, indicating that our method achieved much better performance than these baselines. Therefore, our method can effectively infer the GC in multi-agent (or multi-element) systems with partially known structures.

heading: 6.2 Multi-animal trajectory datasets
text: We here analyzed biological multi-agent trajectory datasets of bats, birds, mice, and flies and obtained new biological insights using our framework (for the results of flies, see Appendix H). We used the same ABM as used in the boid dataset. In real-world data, since there was no ground truth, we used the hyperparameters of the boid simulation dataset. As a possible verification method, our method can be verified by investigating whether the GC result follows the hypothesis using mice and flies (Appendix H) datasets, which were controlled based on scientific knowledge. Next, we show that
our methods as analytical tools can obtain new insights from birds and bats datasets based on the quantitative results. Compared with the methodologies mentioned in Section 5 (i.e., uninterpretable information-theoretic approaches or using non-causal features), our method has advantages for providing local interactions: interpretable signs of Granger-causal effects over time (i.e., our findings are all new).
Mice for verification. As an application to a hypothesis-driven study, for example, we show the effectiveness of our method using three mice raised in different environments. The hypothesis is that, as is well known (e.g., [52]), when grown in different cages, they are more socially novel to others, thus more frequently attractive and repulsive movements will be observed. In this experiment, we regarded the same/different cage as a group pseudo-label, and confirmed that our method could extract its features in an unsupervised manner. We analyzed the trajectories of three mice in the same and different cages measured at 30 Hz for 5 min each (see also Appendix H). As shown in Figure 2A, our method extracted a significantly larger duration in the different cage than that in the same cages for both movements (p < 0.05; p is a statistical p-value), whereas GVAR [45] did not in repulsive movements (p > 0.05) but did in attractive movements (p < 0.05) and extracted too much interaction. The main reason for the too much interaction in GVAR was the overdetection of the attraction and repulsion as shown in Figure 1 right (black break line). The statistical analysis and videos are presented in Appendix H and the supplementary materials. Our methods characterized the movement behaviors before the contacts with others, which have been previously evaluated (e.g., [74]).
Growing birds. Animals grow while interacting with other individuals, but the directed interaction between young individuals has not been fully investigated as longitudinal (i.e., long period) studies. Here, as an example, we analyzed the flight GPS (two-dimensional) trajectories of three juvenile brown boobies Sula leucogaster over six times for 34 days (11.91± 0.09 [h] for each day), which were recorded at 1 Hz. We segmented two or three bird trajectories within 1 km and during moving (over 1 km/s) each other, in which interactions were considered to exist, and obtained 25 sequences of length 367 ± 278 [s] (for details, see Appendix H). Results of inferring GC in Figure 2B show that on the first measurement day, the most frequent directed interactions were observed between ID 1 and 2 (particularly two individuals had more repulsions and ID 1 had fewer attractions). On the
other hand, in the second and subsequent measurements, it was observed that the most interacting individuals changed every measurement day. One possible factor of the decrease in the duration of interactions (especially repulsion) may be the habituation with the same individuals. Measurements and analyses over longer periods will reveal the acquisition of social behavior in young individuals.
Wild bats. As an example of an exploratory analysis, we applied our method to three-dimensional trajectories of eastern bent-wing bats Miniopterus fuliginosus that left a cave (some bats returned to the cave). Although some multi-animal studies have investigated leader-follower relationships (see also Section 5), those in wild bats are unknown. We used two sequences with 7 and 27 bats of length 237 and 296 frames, respectively, which were obtained via digitizing the videos at 30 Hz. Details of the dataset are given in Appendix H. As a result, among 138 interactions of all 34 individuals within the leaving and returning groups, there were 46 interactions where the locationally-leading (i.e., flying forward) bats repelled the following bats in the same direction, 27 interactions where the leading bats were attracted from the following bats, 65 ones with no interactions (the results of the following bats were not discussed because it was obvious; see also example results of 27 bats in Figure 2C). Since bats can echo-locate other bats in all directions up to a range of approximately 20 m [44, 6], the locationally-leading bats can be influenced by the locationally-following bats in the same direction (if no perception, they cannot be influenced). The results suggest that the groups of flying bats would not show simple leader-follower relationships.

heading: 7 Conclusions
text: We proposed a framework for learning GC from multi-animal trajectories via a theory-based ABM with interpretable neural models. In our framework, as shown in Figure 1, the duration of interaction and non-interaction, attraction and repulsion, their amplitudes (or strength), and their timings can be interpretable. In the experiments, our method can analyze the biological movement sequence of mice, birds, and bats, and obtained novel biological insights. One possible future research direction is to incorporate other scientific knowledge into the models such as body inertia (or visuo-motor delay). Real-world animals have certain visuo-motor delays, but they also predict the others’ movements (i.e., the visuo-motor delays may be smaller). This is an inherently ill-posed and challenging problem, which will be our future work.
For societal impact, our method can be utilized as real-world multi-agent analyses to estimate interaction rules such as in animals, pedestrians, vehicles, and athletes in sports. On the other hand, there are some concerns in our method from the perspectives of negative impact when applied to human data. One is a privacy problem by the tracking of groups of individuals to detect their activities and potential interactions over time. This topic has been discussed such as in [60]. Although we did not apply our method to human data, solutions for such a problem will improve the applicability of the proposed method in our society.

heading: Acknowledgments
text: This work was supported by JSPS KAKENHI (Grant Numbers 19H04941, 20H04075, 16H06541, 25281056, 21H05296, 18H03786, 21H05295, 19H04939, JP18H03287, 19H04940, and 21H05300), JST PRESTO (JPMJPR20CA), and JST CREST (JPMJCR1913). For obtaining flies data, we would like to thank Ryota Nishimura at Nagoya University.

Review summaries:
This paper proposes a framework for learning Granger causality with the goal of finding interaction rules from trajectory data of multiple animals. The framework has an "augmented behavior model" based on a conceptual behavior model studied in [53]. The experiments are on a range of datasets, from synthetic to real data of mice, birds, and bats. The main baseline of comparison is GVAR [44], which this method compares against for synthetic data and on real mice data. Based on the experiments, the interaction rules found are in the form of attraction, no interaction, or repulsion between different agents.

The proposed paper introduces a novel approach for learning the Granger causality based on augmented behavioral models. The objective of the discussed method is to extract interaction rules from real-world multi-agent trajectory data with neural networks. The main contributions of this work are a framework for learning Granger causality, methods for efficient and interpretable learning, and the decomposition of motion and navigation. The method is evaluated based on a number of experiments and based on datasets of mice, flies, birds, and bats.    

The manuscript proposes a hybrid method for learning Granger causality that combines data-driven learning methods and theory-based behavioral models for the goal of modeling multi-agent behavior.   The main motivation is that current models either do not consider the data generation process behind the observed time series or if they do, they are tailored to specific animal species, thus having generalization limitations.  More specifically, a Self-Explainable Neural Network (SENN) -like model is used where scientific knowledge on multi-agent interactions is integrated. In this regard, key contributions of the proposed method include: a theory-guided loss function to regularize training and a motion and navigation function Experiments on synthetic and realistic data show the effectiveness of the proposed method.

In this paper, the authors propose a method for learning rules of interaction between individuals from　trajectory data of animals. Their method is based on the Granger causality methodology for time series data but differs from previous studies in that it introduces a new network architecture and a regularization term based on the theory of animal behavior, which enables more accurate learning. The method has been validated using synthetic datasets and then applied to real animal behavior data to gain new knowledge about animal behavior.