Inter-brain neural dynamics in biological and artificial intelligence systems

Authors: Xingjian Zhang, Nguyen Phi, Qin Li, Ryan Gorzek, Niklas Zwingenberger, Shan Huang, John L. Zhou, Lyle Kingsbury, Tara Raam, Ye Emily Wu, Don Wei, Jonathan C. Kao, Weizhe Hong

Abstract

Social interaction can be regarded as a dynamic feedback loop between interacting individuals as they act and react to each other^1,2^. Here, to understand the neural basis of these interactions, we investigated inter-brain neural dynamics across individuals in both mice and artificial intelligence systems. By measuring activities of molecularly defined neurons in the dorsomedial prefrontal cortex of socially interacting mice, we find that the multi-dimensional neural space within each individual can be partitioned into two distinct subspaces—a shared neural subspace that represents shared neural dynamics across animals and a unique neural subspace that represents activity unique to each animal. Notably, compared with glutamatergic neurons, GABAergic (γ-aminobutyric acid-producing) neurons in the dorsomedial prefrontal cortex contain a considerably larger shared neural subspace, which arises from behaviours of both self and others. We extended this framework to artificial intelligence agents and observed that, as social interactions emerged, so too did shared neural dynamics between interacting agents. Importantly, selectively disrupting the neural components that contribute to shared neural dynamics substantially reduces the agents’ social actions. Our findings suggest that shared neural dynamics represent a fundamental and generalizable feature of interacting neural systems present in both biological and artificial agents and highlight the functional significance of shared neural dynamics in driving social interactions.

Social interaction is an adaptive process that fundamentally supports the survival and growth of nearly all animal species. During social interaction, the behavioural decisions of participating individuals are not isolated but rather intrinsically linked—the actions and internal states of individuals are continuously influencing and adapting in response to one another^1,2^. Traditionally, studies of social interaction largely focused on a single brain and how it processes and reacts to social inputs. However, this approach may not fully capture the dynamic nature of social behaviour. To better understand how the social brain functions, an alternative approach is to treat all participants of an interaction as a single integrated system and to measure neural activities simultaneously across brains. As neural activities from multiple brains essentially share the same timeline, this approach may reveal emergent neural properties that reflect the reciprocal nature of the interaction.

Studies using non-invasive techniques have demonstrated that shared neural dynamics, such as inter-brain synchrony, emerge across participants during social engagement in humans and other primates^2–8^. Experiments using in vivo calcium imaging and electro-physiology have recently demonstrated inter-brain neural correlations in socially interacting mice^9^ and bats^10^ with single-cell resolution. Yet, the underlying neuronal components remain poorly understood. It remains unclear (1) whether and how specific neuronal cell types may differentially contribute to inter-brain neural dynamics; and (2) whether and how shared neural dynamics manifest in a higher-dimensional neural space. Although cell-type specificity and multi-dimensional state spaces are two of the most fundamental aspects in understanding any neural processes, they have not been studied with respect to inter-brain neural dynamics in any brain regions in any species.

The discovery of inter-brain neural dynamics across multiple species further raises the important question of whether this represents a fundamental property that is inherent to any interacting agents, including artificial ones. Recent advances in deep reinforcement learning have led to self-evolving artificial intelligence (AI) systems that improve task performance through interaction with an environment or other agents^11–14^. This presents an exciting opportunity for investigating whether the complex dynamics of social interactions may also emerge when multiple artificial agents interact with each other and whether this can be used to model social interactions in biological systems. Given the possibility that artificial agents may learn to interact with each other in ways that are analogous to biological organisms, we explored whether this might be driven by dynamics in their neural networks that resemble those observed in biological systems.

Inter-brain dynamics in select cell types

The dorsomedial prefrontal cortex (dmPFC) is important for encoding social information and regulating social behaviour^15–22^. Our previous studies have demonstrated that mice exhibit a correlation of aggregated neural activities in the dmPFC across interacting individuals^9^. As glutamatergic neurons account for the majority of cortical neurons^23^, this raised the hypothesis that inter-brain correlations arise primarily from the activity of glutamatergic neurons. To test this, we used in vivo microendoscopic imaging to record calcium activity simultaneously in two mice engaging in free interactions (Fig. 1a,b). We separately examined glutamatergic and GABAergic subpopulations of dmPFC neurons using the calcium sensor GCaMP6f in a cell-type-specific manner (CaMKII-GCaMP6f and mDLX-GCaMP6f, respectively) (Fig. 1c,d and Extended Data Fig. 1a–d). We recorded 9,006 CaMKII-positive glutamatergic neurons from 13 pairs of mice and 3,331 mDLX-positive GABAergic neurons from 14 pairs of mice. During these sessions, mice spent around 49% of the time engaging in active behaviours (Fig. 1e), and among these, around 48% involved social behaviours (Fig. 1f,g).

We first measured the Pearson correlation of aggregate (mean) activities of dmPFC glutamatergic neurons between two mice and found that it was greater than chance (Fig. 1h,j). This inter-brain correlation was not due to the autocorrelation, as the cross-correlation of the neural activity, which peaked at 0 s, was disrupted in phase-randomized signals (Fig. 1k and Extended Data Fig. 1g–j). Additionally, after removing autocorrelations by pre-whitening neural activities, the processed signals remained correlated across brains (Extended Data Fig. 1n–s).

Similar to glutamatergic neurons, dmPFC GABAergic neurons were also correlated across brains (Fig. 1i,m,n and Extended Data Fig. 1e,f). Surprisingly, although GABAergic neurons constitute a minor fraction of dmPFC neurons^23^, they exhibited substantially higher correlations across brains than glutamatergic neurons (Mann–Whitney two-sided test, P = 0.0023). Whereas glutamatergic neurons were mainly correlated at slower timescales, GABAergic neurons consistently displayed a higher inter-brain correlation across a wide range of timescales (Fig. 1l,o and Extended Data Fig. 1k–m). These inter-brain correlations were not due to autocorrelations (Fig. 1n,o and Extended Data Fig. 1n–s). In addition, the observed difference in inter-brain correlation between GABAergic and glutamatergic neurons was not due to differences in average firing rates (Extended Data Fig. 1t–x and Supplementary Note 1). By contrast, inter-brain correlation was significantly reduced when mice were separated by a physical divider (separation session) (Extended Data Fig. 1e–h). Moreover, when we examined the inter-brain correlation between mice from different dyads that did not directly interact, the correlation was not significantly above chance for both cell types (Extended Data Fig. 1e,g). Thus, the observed inter-brain correlations depend on direct, ongoing social interaction.

One possible explanation for the higher inter-brain correlation in GABAergic neurons is that (1) GABAergic neurons more robustly represent social behaviour; and (2) GABAergic neurons that represent social behaviour contribute to inter-brain neural correlation. Using receiver operating characteristic (ROC) analysis, we found that 45% of GABAergic neurons were responsive during social interaction in general (any social behaviours), compared to only 10% of glutamatergic neurons (Fig. 1p). This difference was not as pronounced for non-social behaviours. Moreover, we found, using support vector machine (SVM) classifiers, that although population activities in both cell types could decode individual social behaviours (for example, attack and investigation) from baseline or between different social behaviours, GABAergic populations achieved significantly higher performance (Fig. 1q–t and Extended Data Fig. 2). Additionally, GABAergic neurons could better distinguish social behaviours versus a non-social behaviour, self grooming. Thus, social information is more robustly represented in the activities of GABAergic neurons than glutamatergic neurons.

As GABAergic neurons encode both social and non-social information, we next examined whether GABAergic neurons encoding social behaviour contributed to inter-brain neural correlations. We found that inter-brain correlation was significantly reduced after computationally removing GABAergic neurons encoding social behaviour, whereas removing GABAergic neurons encoding non-social behaviour led to an opposite change in inter-brain correlation (Fig. 1u, Extended Data Fig. 1y and Supplementary Note 2). By contrast, removing glutamatergic neurons encoding social behaviour did not significantly decrease inter-brain correlation (Fig. 1v and Extended Data Fig. 1z). Together, these results support our hypothesis that GABAergic neurons are more strongly modulated by social information compared with glutamatergic neurons and that GABAergic neurons encoding social behaviours contribute significantly to inter-brain neural correlation.

Inter-brain shared neural subspace

Prior studies of inter-brain neural dynamics primarily analyse the aggregate activity of all neurons within a brain region^2^. To determine whether there is higher-dimensional structure to inter-brain neural dynamics, we used partial least squares correlation^24^ (PLSC) to identify sets of orthogonal dimensions in neural state space that maximize the cross-covariance of neural activity between the two mice (Fig. 2a,b, Methods and Supplementary Note 3). Specifically, we used the singular value decomposition to compute the left and right singular vectors of the cross-covariance matrix, which were orthogonal bases that captured the shared neural covariance of the interacting mice (Methods). We then projected the neural activity of each mouse onto their respective subspaces to compute neural activity patterns within these shared dimensions (Fig. 2b,c). We identified the neural dimensions whose activity was significantly correlated above chance (Fig. 2d) and refer to them as ‘shared neural dimensions’ and to their associated subspace as the ‘shared neural subspace’ (Fig. 2c,e). We then computed the complementary ‘unique neural subspace’ by performing principal component analysis (PCA) on neural activity projected into the null space of the shared neural subspace (Fig. 2b,c,f). Together, this partitioned the neural activity of each mouse into two neural subspaces—a shared neural subspace that captures shared neural dynamics across mice, and a unique neural subspace that captures activity unique within each mouse (Fig. 2a–c).

Using this approach, we found that both GABAergic and glutamatergic neurons contained several shared neural dimensions during interaction sessions but not separation sessions (Fig. 2e–h, Extended Data Fig. 3a and Supplementary Note 4). We confirmed that, in both male and female pairs, activity of the top shared neural dimension had higher-than-chance correlations, whereas the activity of the top unique neural dimension did not (Fig. 2i,j, Extended Data Fig. 3b–i and Supplementary Notes 5 and 6). Notably, we found that the neuronal ensembles that participated in the top shared and unique neural dimensions were largely distinct (Fig. 2k, Extended Data Fig. 3j and Supplementary Note 7). Further, largely different neuronal ensembles gave rise to different shared neural dimensions (Fig. 2l and Extended Data Fig. 3k). Thus, the shared and unique neural subspaces emerge from different neural ensembles.

Shared subspace in GABAergic neurons

We found that, compared with glutamatergic neurons, GABAergic neurons generally had a larger number of shared neural dimensions across mice (Fig. 2g,h and Extended Data Fig. 3a; for analyses using deconvolved spikes, see Extended Data Fig. 4a–f and Supplementary Note 8). This difference was not simply due to a greater number of neural dimensions in GABAergic neurons within each mouse—in fact, GABAergic population activity was lower-dimensional than that of glutamatergic neurons (Fig. 2m). This result was also not due to a different number of neurons between the two populations, as we observed consistent results after down-sampling glutamatergic neurons to match the average number of GABAergic neurons (Extended Data Fig. 4g–r). Furthermore, the number of shared neural dimensions during social interaction moments were substantially higher than during non-social moments (Extended Data Fig. 5a–l). By contrast, during separation sessions, we found no significant shared neural dimensions for either cell type (Fig. 2g,h). Together, dmPFC GABAergic neurons, compared to glutamatergic neurons, contain more shared information across mice during social interaction.

Moreover, approximately 30% of the total neural variance in GABAergic neurons was shared, whereas only around 5% was shared in glutamatergic neurons (Fig. 2n). This difference was not simply due to a larger number of shared neural dimensions—when we examined the top (the first) shared neural dimension, PLSC1, it captured 3.5 times more neural variance in GABAergic neurons than glutamatergic neurons (Fig. 2o) and these top dimensions (PLSC1) exhibited higher correlations (Fig. 2i). This difference between GABAergic and glutamatergic neurons was similarly observed in both male and female mice (Extended Data Fig. 3f,g and Supplementary Note 5). Thus, the larger amount of neural variance explained by GABAergic neurons is not only due to a greater number of shared neural dimensions but also because of a larger amount of variance explained by individual dimensions. The shared and unique neural subspace also displayed distinct levels of stability across timescales and across different interaction partners (Extended Data Figs. 6 and 7 and Supplementary Notes 9 and 10).

Using partial least squares regression (PLSR), we modelled the shared and unique neural subspaces using annotated behaviours from the corresponding mice. Notably, social behaviours accounted for a larger fraction of variance than non-social behaviours in the shared neural subspace (Fig. 2p). Conversely, non-social behaviours explained a larger fraction of variance than social behaviours in the unique neural subspace (Fig. 2q). Thus, although the shared neural dimensions were identified purely using neural activity without any behavioural information, they captured important social behaviour-related information about both mice. We found that mice engaging in more mutual social interaction had significantly higher variance in the shared neural subspace (Fig. 2r and Extended Data Fig. 5m). Mice also exhibited higher inter-brain correlations during mutual (bidirectional) social interactions than during unidirectional social moments (only one mouse displaying social behaviour) and non-social moments (Extended Data Fig. 5a–l). This was consistent with our observation that the top shared neural dimension was associated with multiple types of social behaviours, especially aggression-related ones such as attack (Fig. 2s and Extended Data Fig. 8a).

Aggressive interaction is one of the most common forms of mutual social interaction, during which both mice closely attend to each other (Extended Data Fig. 8b). In our experiments, a fraction of mice exhibited high levels of antagonistic behaviour, and these antagonistic mouse pairs displayed strong inter-brain correlations (Extended Data Fig. 8c–l and Supplementary Note 11). Inter-brain correlations in GABAergic neurons were higher during aggressive moments than during non-aggressive social behaviour (Extended Data Fig. 8m–p). Moreover, mouse pairs exhibiting higher levels of aggression had a larger shared neural subspace within GABAergic neurons (Fig. 2t). This may reflect the increased attention required during aggression, potentially leading to stronger neural coupling. Thus, the shared neural subspace is strongly influenced by mutual social interaction between mice, particularly during behaviours that involve close mutual attention (such as aggressive interaction).

Intra- versus inter-brain correlations

Shared neural activity across animals may depend on neural representations of social events as well as correlations among neurons within the same brain. To evaluate the correlation among GABAergic neurons within individual mice, we performed PCA on the neural activity of each mouse. We found that PC1 captured more neural variance than chance (Fig. 3a). Similarly, the mean correlation between each neuron and PC1 was higher than chance (Fig. 3b), suggesting substantial intra-brain correlations among GABAergic neurons. Notably, this intra-brain correlation was positively correlated with the inter-brain correlation between the top shared neural dimensions across mice (Fig. 3c). By contrast, glutamatergic neurons comparatively captured less neural variance in PC1 during social interaction and exhibited substantially lower mean correlation between each neuron and PC1 (Extended Data Fig. 9a,b), suggesting a lower level of intra-brain correlation. There was no correlation between the intra- and inter-brain correlations in glutamatergic neurons (Extended Data Fig. 9c–f).

Although the intra-brain correlation during separation sessions also exceeded chance levels, it was not positively correlated with inter-brain correlations (Fig. 3d). Moreover, within the interaction session, although intra-brain correlation was significantly higher than chance during both social and non-social moments (Extended Data Fig. 9j,k), a positive correlation between inter-brain coupling and intra-brain correlation was observed only during social, but not non-social, moments (Fig. 3e,f). Thus, intra-brain correlation alone is not sufficient to explain shared neural activity in the absence of social behaviour.

The observation of intra-brain correlation during both interaction and separation sessions raises the question of whether the same or different neuronal ensembles contribute to intra-brain correlations during social versus non-social periods. Using a PCA–ICA (independent component analysis) approach^19,25^, we identified different clusters of neuronal ensembles that exhibited intra-brain correlations in interaction versus separation sessions (Fig. 3g–i). Clusters that were correlated during interaction sessions showed reduced correlations during separation sessions, and vice versa (Fig. 3j,k). Similarly, ensembles identified during social moments within interaction sessions displayed stronger intra-brain correlations during social moments than during non-social moments, and vice versa (Fig. 3l,m and Extended Data Fig. 9g,h). Thus, intra-brain synchrony during social interaction arises from specific neuronal ensembles that are unique to the social context. Indeed, social ensembles were primarily tuned to individual social behaviours, whereas non-social ensembles were predominantly associated with non-social behaviours (Fig. 3n and Extended Data Fig. 9i).

Compared with glutamatergic neurons, a larger fraction of GABAergic neurons contributed to intra-brain ensembles (Extended Data Fig. 9m,n), even after subsampling the glutamatergic population to match the number of GABAergic neurons (Extended Data Fig. 9o,p). Within GABAergic populations, neural ensembles exhibited significantly higher intra-ensemble correlation than within glutamatergic populations (Extended Data Fig. 9q–t). Thus, GABAergic neurons contribute more prominently to intra-brain ensembles compared to glutamatergic neurons.

To further examine the relative contribution of intra-brain correlation and cross-animal temporal coupling to inter-brain correlation, we used different temporal shuffling schemes^19^ that selectively disrupted each component (Fig. 3o and Methods) and quantified the reduction of covariance in the top shared neural dimension (PLSC1). The covariance of PLSC1 was significantly reduced when either component was disrupted (Fig. 3p). By contrast, during separation sessions, PLSC1 covariance was reduced only when disrupting intra-brain correlation but not temporal coupling (Extended Data Fig. 9l). Thus, inter-brain synchrony during social interaction emerges from both within- and between-animal neural coupling.

Behavioural subspaces between animals

During social interactions, animals sometimes display simultaneous and coordinated motor actions. This raises the question of whether shared neural dimensions merely reflect coordinated behaviours. Using a deep learning-based key-point tracking algorithm^26^, we tracked the frame-by-frame poses of both mice in a pair (Fig. 4a and Extended Data Fig. 10a). We also derived additional pose and action features that captured their temporal dynamics, and measured head orientations of mice using a digital gyroscope attached to miniaturized microendoscope (Fig. 4a,b and Methods). These yielded a comprehensive set of 75 behavioural features, which required 38 principal components to capture 90% of the total variance, higher than that of the manual behaviour annotations (12 principal components) (Fig. 4c). These features could decode all major behaviour types above chance, indicating that these features captured relevant information about manually annotated behaviours (Extended Data Fig. 10b). Using generalized linear models (GLMs), we found that activities of GABAergic and glutamatergic neurons explained the majority of behavioural features above chance, and most strongly represented the derivatives of body position (speed and acceleration) and posture features (distance between body and tail) (Extended Data Fig. 10c,d and Supplementary Note 12).

We next examined whether behaviour dimensions were coordinated or uncoordinated between individuals. Unlike neural activity space, the behavioural space is manually constructed and may contain partially redundant features (Methods). Thus, we used canonical correlation analysis (CCA)^27^ instead of PLSC, because CCA removes redundant information prior to identifying correlated behavioural dimensions across mice (Methods). We refer to these significantly correlated dimensions as the ‘coordinated behavioural subspace’, and to the null space of these dimensions as the ‘uncoordinated behavioural subspace’ (Fig. 4d).

Across 27 pairs of mice, while most pairs contained at least one coordinated behavioural dimension (Fig. 4e), the coordinated dimensions captured only 9% of the total variance in the behavioural space (Fig. 4f), suggesting that mice predominantly engage in uncoordinated movements even during social interaction. Although aggressive behaviours (for example, attack or chase) by one mouse frequently co-occurred with defensive responses (such as defend or escape) by the other, during these moments only 20% of the total behavioural variance was attributed to coordinated dimensions, whereas 80% was attributed to uncoordinated dimensions (Extended Data Fig. 11a).

We found that different behavioural features contributed differentially to individual coordinated dimensions within each mouse (Fig. 4g). To further determine whether the coordinated dimensions emerge from the same behavioural features and, thus, identical behaviours across mouse pairs, we computed the cosine similarity of the significant feature weights in the top coordinated behavioural dimension (CC1) between mouse pairs (Fig. 4h). The cosine similarity was not significantly above chance for most mouse pairs, suggesting that in the coordinated behavioural subspace, mouse pairs engaged in correlated but non-identical behaviours (Fig. 4h).

Relating behavioural and neural subspaces

We next explored how the coordinated and uncoordinated behavioural subspaces contribute to activity in the shared and unique neural subspace. We found that the full behavioural space of both mice accounted for approximately 28% of the variance in the shared neural subspace in GABAergic neurons, significantly more than the approximately 13% in the full neural space (Fig. 4i). In comparison, the full behavioural space accounted for less variance in the unique neural subspace than in the full neural space. These suggest that behavioural information is enriched in the shared neural subspace (Fig. 4i and Extended Data Fig. 12a).

To examine the contribution of the coordinated or uncoordinated behavioural subspace to the shared neural subspace (Fig. 4j), we used PLSR to model the shared neural subspace using the coordinated or uncoordinated behaviours. Within the part of shared neural subspace that can be explained by the full behaviour information, only a fraction was attributed to coordinated behaviour (Fig. 4k and Extended Data Fig. 12b). This indicates that a sizeable remaining fraction of the shared neural subspace was not explained by the coordinated behaviours. Indeed, the subject’s own (self) uncoordinated behavioural subspace and the partner’s uncoordinated behavioural subspace both contributed significantly and non-redundantly to the shared neural subspace (Fig. 4l and Extended Data Fig. 12c). In particular, among the part of shared neural subspace that could be explained by behaviour, the partner’s uncoordinated behaviours accounted for 20% of the neural variance. Compared with the total neural space, this partner representation was stronger in the shared neural subspace (Extended Data Fig. 12f,g). Thus, the shared neural subspace does not merely reflect coordinated behaviour; self and partner-uncoordinated behaviour also contributed to a significant fraction of the shared neural subspace.

By contrast, we hypothesized that the unique neural subspace primarily emerges from self, rather than partner, behaviour (Fig. 4m). Indeed, self-uncoordinated behaviour contributed substantially more to the unique neural subspace than coordinated behaviour or partner behaviour, which was near chance level (Fig. 4n,o and Extended Data Fig. 12d,e). Thus, whereas the shared neural subspace emerges from both self and partner behaviours, the unique neural subspace emerges largely from one’s own behaviours. In particular, the representation of partner-uncoordinated behaviour in the shared, but not unique, neural subspace provides conclusive evidence that dmPFC neurons can specifically encode the behaviours of others that are not correlated with one’s own behaviour and that this representation is enriched in the shared neural subspace (with a 2.3-fold increase; P < 0.0001; Wilcoxon matched-pairs signed rank test).

Social interaction in multi-agent systems

To study interactions between AI systems, we trained artificial agents using multi-agent reinforcement learning^13,14^ (MARL) in a ‘social’ environment. Using RLlib^28^, we created two agents, referred to as the explorer and chaser. The explorer’s goal was to explore new tiles, whereas the chaser aimed to block the explorer by colliding with the explorer (Fig. 5a,b). Each agent had partial vision of the environment (Fig. 5b). Both agents were trained to maximize their own rewards with both non-social (exploring new tiles) and social (collision) goals (Methods). Each agent was implemented using an actor–critic architecture with a 256-unit recurrent neural network (RNN). The explorer and chaser did not share any network parameters. Agents were trained using proximal policy optimization^29^ (PPO) for 20,000 training iterations, each with 4,000 environment steps.

Agents started by making random decisions and successfully learned to perform the task through training. Chasers’ rewards increased, commensurate with an increase in collisions and the percentage of time it kept the explorer within its vision (Fig. 5c,e,f), indicating the emergence of following behaviour. Explorers exhibited a steep increase in reward in the initial phase of training (Fig. 5d) with increased new tile explorations (Fig. 5h), but as chasers improved, explorers received less reward at later training epochs (Fig. 5d). As a control, we created a ‘non-social’ environment by removing rewards for social interaction (collision) and partner vision, while keeping the rewards for exploring new tiles (Methods). In this environment, agents learned to acquire an increasing number of tiles (Fig. 5g), but social collisions remained similar to the untrained network (Fig. 5e). The chaser did not actively approach the explorer, as the time it kept its partner in vision remained as low as the untrained network (Fig. 5f). Thus, agents did not engage in social interaction in a non-social environment.

The performance of one agent is directly coupled to the performance of their co-trained partner agent. To have a consistent benchmark across all performing agents and across different training stages, we also evaluated agent performance in social and non-social environments against random agents (Fig. 5i–n). When playing against a random agent, social chasers (chasers trained in the social environment) consistently outperformed non-social chasers (chasers trained in the non-social environment) (Fig. 5i–k). Meanwhile, social explorers outperformed non-social explorers in increasing their distance from the random chaser, and non-social explorers were better at exploring new tiles (Fig. 5m,n). By contrast, non-social chasers outperformed social chasers in new tile exploration (Fig. 5l). Together, these results show that agents exhibited social behaviour when trained in a social environment, but not when trained in a non-social environment.

To further characterize agents’ behaviour, we measured the chaser’s moving directions relative to the position of the partner (Methods). In the initial training phase, chaser agents exhibited largely random movement directions relative to the other agent (Fig. 5o). After training, chasers displayed an increased tendency to move towards their partner (that is, within 60°; Fig. 5p–r). This suggests that through training, MARL-trained agents developed strategies based on their roles and optimized rewards.

Inter-brain dynamics in artificial agents

We next analysed the activity dynamics within agents’ neural networks that underlie their behavioural strategies (Fig. 6a). We examined an agent’s internal representation of their partner by decoding social events (collision) and partner behaviours (approach and escape) using an agent’s RNN activity. We found, using SVM classifiers, that a social agent’s RNN activity could decode both social events (collisions) and partner behaviours (approach or escape) above chance (Fig. 6b–e). By contrast, a non-social agent’s activity could not decode the other agent’s behaviours.

We next tested whether social agents exhibited shared neural dynamics resembling socially interacting mice. Indeed, we found that trained social agents had significantly more shared neural dimensions and a higher PLSC1 correlation than non-social agents (Fig. 6f–h). The emergence of a shared neural subspace among social agents is not due to a similar vision of the input space—when we provided the non-social agents with the same vision input as in the social task or the full vision of the arena, they contained nearly zero shared dimensions and significantly less correlated PLSC1 compared to social agents (Fig. 6g,h and Extended Data Fig. 13a,b). Shared neural dimensions also did not simply arise from identical actions of the two agents (Extended Data Fig. 13c,d). Thus, shared neural dynamics predominantly arise from the social interaction of the task rather than shared inputs of the same environment or identical actions. Inter-agent shared neural dimensions were also observed in three additional social AI environments, which featured varying collaborative and competitive goals, richer social environments with complex visual inputs, and a more diverse behavioural space (Supplementary Note 13).

We showed earlier that activity in different neuronal ensembles gives rise to the shared and unique neural subspaces in mice. When we examined RNN units that exhibited a significant contribution to PLSC1 and U1, we found that largely different units contributed to shared and unique neural dimensions (Fig. 6i). Similarly, largely different units gave rise to different shared neural subspaces in social agents (Fig. 6j).

Similar to interacting mice, disrupting intra-agent correlations or temporal coupling resulted in a substantial reduction in the covariance of PLSC1 (Fig. 6k). In comparison, disrupting only intra-agent correlations, but not temporal coupling, reduced the covariance of PLSC1 in non-social agents (Extended Data Fig. 13e), similar to that observed in separation sessions of mice (Extended Data Fig. 9l). This suggests that in both rodents and artificial agents, shared neural dynamics emerge from both intra-brain or intra-network correlation and temporal coupling across agents during social interaction.

We next explored what behavioural information is encoded in the shared neural subspace of trained agents. Partner behaviour was represented in social agents, but not in non-social agents (Fig. 6l), and this representation was stronger in the top shared neural dimension (Fig. 6m), resembling the findings in rodents (Extended Data Fig. 12f,g). Notably, the chaser exhibited a stronger partner representation compared with the explorer (Fig. 6n and Extended Data Fig. 13f), suggesting that partner representations are learned to serve specific behavioural strategies. Remarkably, when we examined mice that were more aggressive and exhibited a higher tendency to chase the other mouse in a dyadic interaction, these ‘chaser’ mice also exhibited a stronger partner representation (Fig. 6o), echoing the asymmetric neural properties observed in artificial agents.

We hypothesized that chasers with a stronger representation of the explorer’s behaviour may also exhibit a higher frequency of collisions. To investigate this, we quantified chasers’ neural variance explained by the explorer’s behaviour in their neural action subspace (Methods). Our results revealed a substantial increase in partner representation during training for social chasers, while non-social chasers displayed a decrease (Fig. 6p). This increase in partner representation was strongly correlated with improved chaser performance, with chasers who had a higher partner representation achieving more collisions (Fig. 6q). This correlation was not observed in chasers trained in the non-social environment (Fig. 6r). These findings suggest that a better representation of one’s partner correlates with enhanced task performance.

Finally, we examined whether the neural components that underlie shared neural dynamics causally contribute to social behaviours. Although it is technically challenging to specifically remove high-dimensional shared neural dynamics in animals, we can computationally remove shared neural dimensions in artificial agents and assess their effect on social interactions. We removed the top 10 shared neural dimensions from the chaser by projecting the agent’s RNN activation into the null space of these PLSCs prior to computing the agent’s action. This removed the contribution of these shared neural dimensions to the policy decision-making process while leaving the remaining neural dynamics unique to one individual intact. This perturbation resulted in a significant decrease in both the number of collisions and the per cent time chasers kept partners within their field of vision compared with the original agent pairs and a significant increase in the average distance to the other agent (Fig. 6s–u). Therefore, removing shared neural dimensions significantly impaired social behaviours. As a control, removing a comparable but random amount of variance from chasers did not result in a significant decrease in social behaviour compared to the original agent pairs (Fig. 6s–u and Extended Data Fig. 13g). Thus, the neural components underlying the shared neural subspace have a crucial role in the emergence of social interaction in agents.

Discussion

Although neural correlations between individuals have been observed in human subjects and animal studies^2,5,9,10^, the underlying neural components are poorly understood. Here we examined inter-brain dynamics in specific, molecularly marked neuronal subpopulations. We found that GABAergic neurons in the dmPFC were considerably more correlated across individual mice compared with glutamatergic neurons. Such distinctions may stem from different response properties, including coding properties, of GABAergic and glutamatergic neurons. Indeed, in many cortical areas, glutamatergic and GABAergic neurons are differentially modulated by sensory inputs and brain states and exhibit different response profiles^30–32^. In particular, mPFC GABAergic neurons have been shown to regulate social behaviours^15,17,18^. Our finding expands our cellular-level understanding of inter-brain dynamics, and provides a foundation for future investigations of other neuronal subpopulations defined by molecular markers or connectivity.

Although previous studies showed that inter-brain dynamics may occur in certain frequency bands^33,34^, these measurements were typically derived from aggregated, unidimensional activity. Using PLSC analysis, we identified a high-dimensional neural subspace that is shared between two brains. Compared to unidimensional measurements, the multi-dimensional characterization of the unique and shared neural subspace can be used as an effective framework to represent the relationship of neural dynamics across social animals. Using this approach, we show that GABAergic neurons contain a substantially larger shared neural subspace compared to glutamatergic neurons. This reveals a previously unappreciated distinction in the structure of inter-brain synchrony between the two neuronal populations. Although the shared and unique neural subspaces were identified without behavioural information, unique neural dimensions captured non-social information within each mouse, whereas the shared neural dimensions captured social behaviour-related information about both mice. In particular, shared neural dynamics are not simply attributable to temporally coordinated behaviour between two individuals, but also emerge from the representation of others’ unique (uncoordinated) behaviour. The shared and unique neural subspaces may serve as a population-level mechanism for selectively routing activity from the mPFC to downstream areas^35^. An intriguing question is whether and how downstream neurons form selective connections to access information from these subspaces to modulate social and non-social behaviours. In humans, inter-brain correlation is observed across many distinct brain regions, depending on social context and task^2–4^. The shared neural subspace that we identified probably reflects a fundamental principle of neural systems that exists in other brain regions and other species including humans.

By modelling social interactions in AI agents, we found that artificial agents developed behavioural strategies that allowed them to interact with each other, and that, analogous to interacting animals, cross-individual shared neural dynamics emerge from interactions between artificial agents. Remarkably, many characteristics of shared neural dynamics also appear to be similar between mice and agents. As observed in interacting mice, shared neural dynamics between agents were not simply due to shared inputs or coordinated action outputs, but arose when the two agents engaged in distinct actions. In both mice and agents, shared neural dynamics arise from neuronal correlations within an individual and interaction-induced temporal coupling. We recognize that AI agents and environments do not fully capture all aspects of real-world animal interactions, such as the detailed behavioural features in mice. Nevertheless, artificial agents are conceptually similar to biological agents in that they (1) receive social inputs; (2) make social decisions that affect other agents; and (3) possess a neural network encoding these inputs and decisions. In addition, although the biological and artificial systems we examined have differences in network (1) in biological systems, we can only examine a small fraction of the entire brain, whereas in artificial agents, we have access to the complete neural network involved in decision-making; and (2) specific cell types in mice do not correspond directly to neural networks in artificial agents—they share common organizational principles, such as representations of environmental stimuli and the actions of the interaction partner. These shared principles are likely to contribute to the notable similarities in the neural dynamics between biological and artificial systems. Such parallels suggest that shared neural dynamics represent a fundamental, generalizable property of interacting neural systems.

Precise perturbation of high-dimensional neuronal components that are specifically involved in shared neural dynamics is technically challenging in animals, making it difficult to determine how inter-brain synchrony functionally contributes to social interaction. An important advantage of artificial agents is that we have full access to the neural network of each agent and can easily perturb the system^36,37^. We found that selectively disrupting the neural components that contributed to shared neural dynamics substantially reduced agents’ social actions, demonstrating the significance of shared neural dynamics in social interaction. Thus, multi-agent systems may serve as a platform not only for understanding emergent behavioural strategies during social interaction but also for testing the causal role of specific neural components, which could be otherwise challenging to examine in animal models. A deeper understanding of AI social interaction may also enhance understanding and designing future social AI models, an important endeavour as AI increasingly becomes a part of everyday life (Supplementary Note 14).

Methods

Animals

All experiments were carried out in accordance with the NIH guidelines and approved by the UCLA Institutional Animal Care and Use Committee (IACUC). Male or female C57BL/6J mice (Jackson Laboratories 000664) at 8–10 weeks of age and 25–30 g of weight were used for microendoscopic calcium imaging and behavioural experiments. Male Vgat^cre/+^ mice (Vgat is also known as Slc32a1) (Jackson Laboratories 028862) were used for the anatomical co-localization experiment in Extended Data Fig. 1a. Vgat-ires-cre mice^39^ were first purchased from Jackson Laboratories (028862) and backcrossed to C57BL/6J to generate a breeding colony. Mice were maintained in a 12 12 h light/dark cycle (light 00–09:00) with food and water ad libitum, and were housed under controlled environmental conditions at a temperature of 21–23 °C and relative humidity of 30–70%. Mice were individually housed for three weeks prior to imaging and behaviour experiments. All experiments were performed during the dark cycle of the mice.

Stereotaxic surgeries

Mice were anaesthetized with 1.5 to 2.0% isoflurane delivered through the SomnoSuite low-flow anaesthesia system during virus injection, GRIN lens implantation, and baseplate mounting. For microendoscopic calcium imaging experiments, we unilaterally injected 300 nl of AAV1-CaMKII-GCaMP6f virus (Addgene) or 500 nl AAV1-mDlx-GCaMP6f^40^ (Addgene) at 40 nl min^−1^ into the dmPFC of the right hemisphere using the stereotactic coordinates (anterior–posterior (AP): +2.0 mm, medial–lateral (ML): +0.3 mm, dorsal–ventral (DV): −1.8 mm from bregma). The mice were given 3 to 5 days for recovery so that the virus was well diffused before the GRIN lens implantation. To plant the 1.0 mm diameter GRIN lens (Inscopix), we created a 1.1–1.2 mm diameter circular craniotomy centred above the virus injection site (AP: +2.0 mm, ML: +0.3 mm from bregma). The GRIN lens was implanted through the craniotomy (DV: −1.65 mm from bregma) and secured to the skull using super glue and dental cement. Mice were given one subcutaneous injection of ketoprofen (4 mg kg^−1^) on the same day of the surgery and ibuprofen in drinking water (30 mg kg^−1^) starting on the surgery day for 4 days. Mice with implants were individually housed. After three weeks of recovery, the mice were ready for imaging. We adjusted the miniscope position while imaging through the GRIN lens to find the ideal imaging plane and mounted the metal baseplate for the miniscope on the top of the skull accordingly. The baseplates were secured using dental cement. After mounting the baseplate, the mice were given two to three days for full recovery before being used for experiments. All mice were handled and habituated to the experimental setup for at least 4 days before experiments. For the mDLX–Vgat co-localization experiment (Extended Data Fig. 1a), 500 nl of AAV1-mDlx-GCaMP6f and 300 nl of AAV5-hSyn-DIO-mCherry (Addgene) were co-injected into the dmPFC of Vgat^cre/+^ mice (AP: +2.0 mm, ML: +0.3 mm, DV: −1.8 mm from bregma). Fourteen days following the injections, brain sections were collected and processed.

Histology

After imaging experiments, the mice were transcardially perfused with 4% paraformaldehyde (PFA). The brain was extracted and further fixed in the same solution for 24 h. Then it was transferred into the 20% sucrose PBS solution for 24 h or until the brain sank for cryoprotection. Finally, 30-μm coronal sections were obtained using the cryostat, mounted on slides, stained with DAPI (1:5,000), and imaged using a fluorescence microscope (Leica DM6 B) to confirm the expression of the virus and the location of the lens implant. To measure the co-localization between mDlx-GCaMP6f-expressing and Vgat::mCherry-expressing neurons, 30-μm coronal sections were obtained with a standard histology procedure described above. Images were acquired using a confocal microscope (Zeiss LSM880). The number of fluorescence-positive neurons and their co-localization were manually quantified.

Behaviour assays

For each of the free social interaction experiments, we paired two stranger male or female mice. The mice were pre-habituated to the miniscopes and the open arena for at least four days consecutively before the experiments. On the experiment day, they were transferred to the test room at least 1 h before for habituation. The room was illuminated by dark red light to minimize anxiety. A new cage (26.5 cm × 11.4 cm × 14.8 cm) with fresh woodchip bedding was used as the arena. The recording setup included 2 miniscopes recording at 15 frames per second (fps) and 1 webcam at 30 fps for top-view behaviour recording, all connected to one computer. The miniscopes were connected to the data acquisition device (DAQ) through a long flexible silicone rubber mini coax cable (Cooner Wire CW2040-3650SR), which allowed the mice to move freely. We used the Miniscope DAQ software to drive and control all three devices, so the timestamps were automatically aligned.

Before each experiment, we mounted miniscopes on the mice, placed them in the neutral arena divided by a solid opaque board in the middle, and let them further habituate for an extra 20 min. We started experiments with a 10-min recording with the divider as a control session, followed by a 20-min free interaction session without the divider. During the free interaction session, we managed the cables carefully so that the mice were not affected by the tangled cables. We used a total of 27 unique males and 12 unique females—14 males were used to form 14 pairs for GABAergic neurons, 13 males were used to form 13 pairs for glutamatergic neurons, 7 females were used to form 8 pairs for GABAergic neurons, and 5 females were used to form 9 pairs for glutamatergic pairs.

Behaviour analysis

Annotation of mouse behaviour.

We manually annotated behavioural videos on a frame-by-frame basis with the assistance of a MATLAB toolkit (https://github.com/pdollar/toolbox) and a custom behaviour annotation software that we developed (https://github.com/hongw-lab/behavior_annotator). We annotated 12 social behaviours (approach, follow, chase, escape, attack, defend, tussle, flinch, threaten, general sniff, sniff face, sniff genital) and 6 non-social behaviours (self groom, dig, climb, explore object, bite object, stand (rearing)). ‘Object’ refers to parts of the cage or experimental setup that mice may interact with, such as the cage ventilation outlet; no additional object was introduced in the experiment. We excluded frames in which mouse behaviours were interrupted by cables or the human experimenter. Annotations were organized in a binary matrix (session duration × 18 behaviours), where each cell contained 1 if the behaviour occurred during that time frame, or 0 if not. Behaviours were mutually exclusive, with only one category active per frame. The fraction of time behaving was calculated as the fraction of the total time that the mice were engaged in any of the 18 behaviours. The fraction of social and non-social time was calculated as the fraction of the amount of time the mice were engaged in social or non-social behaviour out of the total behavioural time. Attack, threaten and chase were categorized as offensive behaviours, while defend, flinch and escape were catagorized as defensive behaviours. The aggressor and submissor within each mouse pair were determined by comparing offensive versus defensive behaviour durations. The mouse with more time spent in offensive than defensive behaviours was classified as the aggressor, while its partner was designated as the submissor.

Posture tracking of mouse behaviour.

We used SLEAP^26^, a deep learning algorithm, to track the posture and position of both mice in our experiments. Our tracking model defined seven nodes (nose, miniscope base, miniscope top (LED), left ear, right ear, body and tail base) and six edges (miniscope base to nose, miniscope base to miniscope top (LED), miniscope base to left ear, miniscope base to right ear, miniscope base to body and body to tail base). Our model was trained in the default top-down multi-animal mode. The model first identified the mice and then estimated the pose of each identified mouse. The model was trained on more than 9,500 manually labelled frames (approximately 10–20% of frames per video, pooled across all mice). The trained model reliably estimated the position of all the nodes of both mice on more than 90% of frames. We manually corrected the incorrectly estimated frames to ensure the posture estimation outputs were accurate. A total of 75 features were derived from the tracking data, and these features are summarized in Supplementary Table 2. Manual annotation and posture tracking provide complementary information about mouse behaviours. SLEAP tracking provides feature-rich context-free objective measurements of animal kinematics. Manual annotation, on the other hand, provides interpretable behavioural events and carries contextual information.

Decoding annotated behaviours from tracking features.

We evaluated the tracking features against the behavioural annotation using a binary SVM classifier. The SVM input was all 75 tracking features, and the output was whether the target behaviour occurred or not. Each annotated behaviour type was classified against all other behaviours. The positive and negative classes were balanced in training and validation, and only behaviours with at least 50 frame occurrences were tested. The performance of SVM classifiers was the average fivefold cross-validation accuracy. Chance-level decoding accuracy was obtained by averaging the fivefold cross-validation accuracies of 100 random permutations.

Microendoscopic calcium imaging analysis

Calcium signal preprocessing.

Calcium fluorescence videos of both mice were simultaneously recorded at 15 fps through Miniscopes (https://github.com/aharoni-lab/miniscope-v4). Raw videos were first fed into the motion-correction algorithm NoRMCorre^41^ to eliminate motion artefacts. We then used the bandpass filter function in ImageJ (filterLarge = 40, filterSmall = 3, percentage of the image size) to remove the fluorescent background from the corrected videos. From the filtered videos, we applied CNMF-E^42^ (constrained nonnegative matrix factorization) to automatically detect and extract regions of interest (ROIs). All ROIs were manually inspected to remove duplicated ROIs and ROIs that did not represent cell bodies. Unless otherwise mentioned, we used the CNMF-E denoised traces for all analyses. Frames of calcium activity were rarely dropped. When ten or fewer consecutive frames were dropped, we linearly interpolated the calcium traces.

While aggressive behaviours (for example, attack and tussling) were associated with increased neural activity in GABAergic neurons (Extended Data Fig. 11c), there was no correlation between neural activity levels and the magnitude of mouse movement (measured by movement speed) (Extended Data Fig. 11b,e). Thus, the observed neural activity differences reflect specific behavioural contexts rather than motion artefacts. Furthermore, we found that aggressive behaviours were not associated with increased activity in glutamatergic neurons (Extended Data Fig. 11d,f). If motion artefacts, rather than the specific behavioural contexts, were affecting neural activities, we would expect such artefacts to also appear in glutamatergic neurons as well. The absence of increased activity in these neurons supports the conclusion that the observed neural activity increases during aggressive behaviour are unlikely to be due to motion artefacts.

Behaviour modulation of single-cell response.

Single-neuron correlates to individual types of behavioural events (such as attack or escape) were analysed through ROC analysis, as previously described^43,44^. These behavioural events were determined from the manually labelled behavioural data. At every time point, a behavioural event either happened (positive label) or did not happen (negative label). This is a binary classification problem. ROC analysis quantifies how well a neuron can detect a behavioural event across several discrimination thresholds. For every threshold, true positive rates and true negative rates were calculated and plotted against each other to produce the ROC curve. The area under the ROC curve (auROC) was used to quantify how strongly a single neuron could predict a specific behaviour. To determine significance, the observed auROC for each neuron was compared to its own null distribution based on circularly permuting neural activity for 2,000 random time shifts. We used a circular permutation (circshift function in MATLAB) to preserve the temporal dynamics of the calcium signals while disrupting the temporal relationship between the calcium signal and behavioural event occurrences. A neuron’s calcium response was considered significant if its auROC value exceeded the 95th percentile of the null distribution (inhibited response if auROC <2.5th percentile, excited response if auROC >97.5th percentile).

Decoding behaviours from neural population activity.

We quantified how behaviours are represented in the neural population activity through binary classification with an SVM classifier. At each time point, the SVM input was a vector of the population calcium activity (a vector whose dimensionality is the number of neurons), and its output was whether the behaviour occurred or not (0 or 1). We balanced negative and positive class labels in the training and validation data. The performance of the SVM was the average cross-validation accuracy across 10-fold cross-validation. Chance levels were computed by calculating the average tenfold cross-validated accuracy values after circularly permuted calcium signals. We constructed a null distribution using 2,000 random time shifts. The neural population representation of behaviour was considered significant when its decoding accuracy exceeded the 95th percentile of the null distribution. Only behaviours with at least 50 frame occurrences were used in this analysis.

For the GABAergic cell type, we used all the neurons recorded per mouse (that is, 119 ± 8 neurons per mouse) to decode behaviours (Fig. 1q–t and Extended Data Fig. 2a–e). For glutamatergic neurons, we down-sampled all recorded neurons (346 ± 23 neurons per mouse) to match the average number of neurons recorded in the GABAergic population in the behaviour decoding analysis (Fig. 1q–t and Extended Data Fig. 2a–e). As an additional control, we also performed an N-matched decoding analysis by down-sampling both cell types to 50 neurons per mouse (Extended Data Fig. 2f–n).

Analysis of the anatomical distribution of neurons in the shared and unique neural subspace.

We investigated the anatomical distribution of neurons that contributed significantly to the top shared (PLSC1) and unique (U1) neural dimensions. Significant neurons were defined as those with weights exceeding 1.5× s.d. from the mean in either the PLSC1 or U1 dimensions. For each group, we calculated the centre of mass of these significant neurons and compared across groups.

Analysis of inter-brain neural dynamics

Pearson correlation of neural activity across brains.

To compute the Pearson correlation of neural activity across brains, we first averaged the neural activity of all recorded neurons for each mouse. We then calculated the Pearson correlation coefficient of these two average traces across the experimental session. We also constructed a null distribution by circularly shifting the average population activity of one mouse relative to the other across 100 random time lags.

Cross-correlation of neural activity across brains.

To quantify the temporal dynamics of neural correlations across brains, we computed the cross-correlation of the average neural activity in paired mice. This involved measuring the correlation between two average population activity traces at time lags ranging from −60 to +60 s. We then plotted these correlation values at their corresponding time lags. We performed this analysis for both interaction and separation sessions. To control for cross-correlations that may arise due to autocorrelation in the signal, we generated a surrogate signal using phase randomization. In brief, we first transformed the original signal into the frequency domain using the FFT, and then randomized the phases across the signal while maintaining its amplitudes. We then transformed it back to the time domain to create the surrogate data and performed cross-correlation on the surrogate data to obtain the phase-randomized control. Phase randomization transforms the input signal into a randomized signal while preserving the power spectrum and autocorrelation. If autocorrelation contributes to inter-brain correlation, phase-randomized signals would still be correlated. We showed that the randomized neural activity, despite having the same autocorrelation (Extended Data Fig. 1i,j), had chance-level inter-brain correlation (Fig. 1j–o). This demonstrates that autocorrelation does not cause inter-brain correlation.

As an alternative control to rule out autocorrelations, we pre-whitened the calcium activity traces to remove autocorrelation before computing inter-brain correlations. For each neuron, we fit an autoregressive integrated moving average model ARIMA(1,0,1)^45^, removed the model-inferred activity, and retained the residual. We confirmed that this effectively removed autocorrelated components of the neural activity. We found that pre-whitened neural signals from two brains were still correlated in both GABAergic and glutamatergic neurons, with the levels of correlation higher in GABAergic neurons. Moreover, we performed our analyses using deconvolved calcium spikes that do not have a slow transient decay, and also observed significant inter-brain correlation. Together, these multiple lines of evidence exclude the possibility that autocorrelation is the source of inter-brain correlation.

Cross-correlation of neural activity across brains at different frequency bands.

To explore neural correlations across brains at different frequency bands, we used the FFT to decompose neural activity into distinct frequency bands. The period (1/frequency) of the frequency bands used for the FFT decomposition of the population average activity was indicated on the x axis in Fig. 1l,o and Extended Data Fig. 1k–m,r,s. This neural activity may refer to the average population activity or other projections of the neural activity depending on the analysis. Cross-correlations were computed across interacting mice for each frequency band and compared with the phase-randomized control as well as cross-pair controls (computing the cross-correlation between neural activity of mice in different pairs).

Analysis of shared neural dynamics using PLSC.

We sought to identify the neural population structure of inter-brain neural coupling beyond average activity. We therefore used a dimensionality reduction method to identify shared dimensions of neural population activity between mice as well as dimensions unique to individual mice. In particular, we used PLSC^24^, a factorization method, to decompose the neural activity of the brain of each mouse into two orthogonal one capturing activity correlated across mice (referred to as the shared neural subspace) and the other capturing activity unique to each individual mouse (referred to as the unique neural subspace).

PLSC performs an optimization to identify two projections of the neural population activity, represented by weight matrices WX and WY, one for each mouse in a pair. PLSC optimizes these weight vectors to maximize the covariance between the neural activities of both mice, as expressed argmaxWX,WYcovWX⊤X,WY⊤Y where X∈Rn1×T is a z-scored matrix containing the neural activity of mouse 1 with n1 neurons across the length of the entire experimental session, T. Each neuron has zero mean and a s.d. of 1. Similarly, Y∈Rn2×T is the neural activity of mouse 2. The solution to this optimization can be found via singular value decomposition. The singular value decomposition of the cross-covariance matrix C of neural activity from both mice C=XY⊤=WXΔWY⊤ where C∈Rn1×n2 is a cross-covariance matrix, Δ represents the diagonal matrix of singular values (in descending order), and WX⊤WX=WY⊤WY=I so that the columns are orthonormal.

We quantified the cross-correlation of the ith PLSC dimension by computing corrWX,i⊤X,WY,i⊤Y, where WX,i is the ith column of the matrix WX. To test whether this dimension was significantly shared between mice, we performed a permutation test. We computed two different null distributions, one of the singular value Δ and one of the cross-correlation of the ith dimension, by temporally permuting the data for 10,000 randomly generated time lags. Dimensions exceeding the 97.5th percentile in both distributions were considered significant. Only the top sequentially significant dimensions were counted—for example, if the dimensions 1, 2, 3 and 5 passed the significant test while the 4th dimension did not, the mouse pairs were considered to only have 3 shared dimensions. For the rest of the dimensions that did not pass the significance threshold, we projected the neural activity of the mouse onto this subspace and performed PCA to establish the unique neural subspace. The unique neural subspace is the null space of the shared neural subspace (that is, the space orthogonal to the shared neural subspace). Both shared and unique neural subspaces are computed independently of behavioural information.

PLSC1 correlation refers to corrWX,1⊤X,WY,1⊤Y, the correlation of the top shared neural dimension. ΔPLSC1 is the PLSC1 correlation minus the chance correlation level computed from the null distribution. U1 correlation refers to the correlation of the top unique neural dimension (PC1 of the unique neural subspace). ΔU1 is the U1 correlation minus the chance correlation level computed from the null distribution.

We used all recorded neurons from both mice to construct the cross-covariance matrix for PLSC (that is, 119 ± 8 neurons in GABAergic neurons and 346 ± 23 neurons in glutamatergic neurons per mouse). To compare GABAergic and glutamatergic neurons, we also down-sampled glutamatergic neurons to match the average number of neurons recorded in the GABAergic population (Extended Data Fig. 4m–r). As an additional control, we performed an N-matched PLSC analysis by down-sampling both cell types to 50 neurons per mouse (Extended Data Fig. 4g–l).

The timescale required for neural data depends on the behavioural context. For our study of free social interaction, which lacks a fixed trial structure, we empirically find that a minimum time window of 30–60 s is needed to generate shared dimensions, while longer time windows produce a more stable neural space. However, for other behavioural contexts, such as those with a repeated trial structure, shorter time windows may be sufficient.

Temporal lags of shared neural activity across interacting mice.

To examine if there are temporal lags in neural signals between interacting mice, we computed both the covariance and correlation of the top shared neural dimension (PLSC1) across lags from −30 to +30 s, using 100 ms intervals. Here, the lag represents the timing neural activity of the more aggressive mouse relative to that of the more submissive mouse within each pair. A positive lag indicates that the neural activity of the more aggressive mouse precedes that of the submissive mouse.

Evaluating the similarity of neural dimensions

Cosine similarity score.

The similarity between neural dimensions (for example, PLSC1, U1) was measured by the absolute value of cosine c1⋅c2c1c2, where c1 and c2 are coefficients of neurons in the neural dimensions of interest. The higher the score, the more similar the two neural dimensions.

Compare neural dimensions from different behavioural moments.

To evaluate the similarity of neural dimensions across different behaviours, we performed PLSC analyses on neural activity during specific behaviours. Specifically, we generated behavioural segments by randomly sampling 8 s of neural activity during a given behaviour and 8 s of neural activity when mice were not engaged in any defined behaviour. Shared and unique dimensions were established from these behavioural segments, and we extracted the coefficients of individual neurons contributing to PLSC1 and U1. For a shuffle control, we temporally permuted these segments 200 times, obtained the PLSC1 and U1 coefficients, and computed the cosine similarity scores between different behaviours. To assess the similarity of neural dimensions within a specific behaviour across time, we used a similar approach, while ensuring that the two sampled segments were at least 10 s apart.

Given that neural activities are dynamically modulated by distinct behaviours, we expect that the shared neural subspace of both neuron types change dynamically but may retain some similarity across select behaviour pairs compared to random distributions. Indeed, we found that in both GABAergic and glutamatergic populations, the weights of these dimensions were distinct across behaviours, while cosine similarities of these weights exceeded those of random distributions for a subset of behaviour pairs (Extended Data Fig. 7c–f). It is important to note that these higher-than-chance similarities do not imply that these shared neural dimensions are identical across behaviours—rather, they change dynamically across behaviours, consistent with the rich dynamics of moment-by-moment behaviours in animals (Supplementary Note 9).

In addition, we used a similar approach to evaluate the inter-brain correlation and the shared neural subspace during social and non-social moments within interaction sessions (Extended Data Fig. 5a–h). Here, we generated behavioural segments by sampling 20 s of neural activity during the social or non-social moments and 20 s of neural activity when mice did not engage in any defined behaviour. These segments were then used to establish the shared neural dimensions.

To explore how the shared neural subspace is modulated by different behavioural feedback from interaction partners, we extracted neural activity during instances when the subject mouse engaged in the same behaviour but the partner responded with different behaviour. We then applied the same procedure as above to examine the resulting shared neural subspace.

Compare neural dimensions when the same mouse interacts with different partners.

We evaluated the stability of the neural dimensions when the same mouse interacted with different partners. Using CellReg^46^, we registered neurons from the same mice across two interaction sessions, each with a different partner. We then calculated the coefficients of these registered neurons that contributed to the shared neural dimensions and used these coefficients to compute the cosine similarity scores, which quantified the similarity of shared neural dimensions across different interaction partners. We found that GABAergic neurons had a similarity score consistently higher than the temporally shuffled controls (Extended Data Fig. 7a,b) and higher than glutamatergic neurons (P = 0.0381, two-tailed Mann–Whitney test). This suggests that, across different pairs, the shared neural subspace of GABAergic neurons is more aligned than for glutamatergic neurons.

Compare neural dimensions at different timescales.

To evaluate the similarity of shared neural dimensions on faster timescales, we applied a sliding window of one minute epoch length. We included only epochs containing at least one significant shared neural dimension and calculated the cosine similarity of the PLSC1 weights between epochs separated by at least one minute. To create a null distribution, we circularly shifted the neural activity of one mouse relative to the other within each pair, using 1,000 randomly generated time lags, and quantified the cosine similarity of the corresponding epoch pairs. For slower timescales, we examined two 10-min segments from each interaction session and computed the cosine similarity of the PLSC1 weights between the segments. Chance levels were determined by temporally shifting neural activity, following the same approach as in the faster timescale analysis.

Considerations of different methods.

While CCA^27^, PLSC and reduced rank regression (RRR) are all statistical methods for analysing relationships between two sets of variables, there are important differences between them. Compared to RRR, PLSC and CCA are somewhat similar to each other; the only difference is that CCA identifies the linear combination of variables that maximizes their correlation, whereas PLSC identifies the linear combination of variables that maximizes their covariance. As CCA may identify correlated neural signals that account for a small fraction of the variance, it is more susceptible to noise. Compared to CCA, PLSC takes into account the variance of the shared signals, making PLSC more robust for noisy, trial-free social interaction data like ours. Nevertheless, when we used CCA to analyse our neural data, it yielded similar conclusions as our PLSC analysis (Supplementary Note 3).

To identify significant shared neural dimensions using CCA, we performed CCA on the z-scored neural activity of both mice within each pair. We then determined the significance of the top shared neural dimension (CC1) by comparing their correlations to the 97.5th percentile of the null distribution generated from 20,000 random temporal shuffles of the data. We found that the top shared dimensions identified using CCA were highly correlated with those identified using PLSC in both GABAergic and glutamatergic populations (Extended Data Fig. 4s). Additionally, in the CCA analysis, the top neural dimension captured approximately 3.5 times more neural variance in the GABAergic subpopulation than in the glutamatergic population (Extended Data Fig. 4t). Thus, our overall conclusions that shared neural dimensions arise during social interaction across cell types, and that the GABAergic subpopulation consistently contains a larger shared neural subspace than the glutamatergic subpopulation, are the same for both methods. Nevertheless, we should note that we do not expect the exact shared subspace to be identical between the two methods, as they are mathematically different—CCA maximizes the projection correlation, whereas PLSC maximizes the projection covariance. CCA may therefore identify dimensions that have correlated activity but account for only a small fraction of the total neural variance (Extended Data Fig. 4v, w), making it more susceptible to noise.

RRR, on the other hand, is typically used for asymmetric analyses, as it models the linear relationship between inputs and outputs by reducing the dimensionality of the regressors, focusing on predictive power. As RRR does not identify correlated signal in a symmetric manner, it is inherently not suited for our analysis.

Neural space analysis using deconvolved spikes.

We performed additional analyses using spike events extracted through CNMF-E and found that our conclusions remain unchanged. Consistent with our original analysis, we found that significantly more shared neural dimensions emerged during interaction sessions than during separation sessions. Additionally, GABAergic neurons overall showed a greater number of shared neural dimensions than glutamatergic neurons across mice. Approximately 28% of the total neural variance in GABAergic neurons was shared, compared to only 2% in glutamatergic neurons, with the top neural dimension capturing approximately 5 times more neural variance in the GABAergic subpopulation than in the glutamatergic population. These results are consistent with those observed using fluorescence signals.

In all of our analyses, the activity of each neuron is z-scored, such that all neurons (both GABAergic and glutamatergic) are normalized to have the same activity mean and s.d. To further examine whether differences in firing rates between the neuron subtypes might contribute to the observed inter-brain correlation difference between GABAergic and glutamatergic populations, we subsampled the top 20% of glutamatergic neurons with the highest firing rates, aligning their average firing rates more closely with those of GABAergic neurons. We found that the inter-brain correlation of GABAergic neurons remained higher than this select group of glutamatergic neurons. As an alternative approach, we also subsampled the top 50% of glutamatergic neurons and the bottom 50% of GABAergic neurons based on average firing rates to achieve more similar firing rates. Similarly, we found that the inter-brain correlation of this selected GABAergic population was higher than in the corresponding glutamatergic population. Together, these results suggest that differences in activity levels between neuron subtypes do not account for the observed inter-brain correlation differences across cell types.

Pearson correlation coefficient analyses following cell deletion.

We computed the contribution of neurons to inter-brain correlation by removing neurons that encode either social or non-social behaviours and quantifying the reduction in the Pearson correlation coefficient. Behaviour encoding was determined using the previously described single-cell ROC analysis. For each behavioural group (social or non-social), the top 20% of neurons, ranked by their auROC values, were selectively excluded. We subsequently recalculated the Pearson correlation coefficient of the average population activity between mice. To control for chance effects, we randomly permuted the auROC values among neurons and removed the top 20% of these permuted neurons before calculating the corresponding Pearson correlation coefficient.

Comparison of neural ensembles that gave rise to different neural dimensions.

To examine whether different shared neural dimensions emerge from the same neural ensembles, we identified significant neurons for each of the top 2 shared dimensions (PLSC1 and PLSC2); these were neurons whose absolute weights were 1.5× s.d. above the weight distributions. We then identified neuronal populations that contributed to only one dimension or to both and plotted their difference in absolute weight across the dimensions (absolute weight in PLSC1 minus absolute weight in PLSC2) as histograms. The same analysis was applied to examine whether the top shared and unique neural dimensions (PLSC1 and U1) emerged from different neural ensembles.

Identification of intra-brain neural ensembles.

We identified intra-brain neural ensembles across sessions using a PCA–ICA method^19,25^. In brief, the number of neural ensembles present was identified as the number of significant principal components, defined as those whose eigenvalues were greater than the 97.5th percentile of the null distribution. This null distribution was constructed by computing the maximum eigenvalue from PCA applied to each of the 1,000 permutations, generated by temporally shifting population activity by a random lag. ICA was then performed on these significant principal components to identify independent correlated neural ensembles. Ensemble activity was computed by projecting the neural activity of all neurons onto ensemble weight vectors derived using the PCA–ICA method (by multiplying the time-by-neuron activity matrix with the neuron-by-ensemble weight matrix). Code can be obtained from https://github.com/tortlab/Cell-assembly-detection/blob/master/.

Contribution of temporal coupling and intra-brain correlation to inter-brain coupling.

To estimate the contribution of intra-brain correlation and temporal coupling via social interaction to inter-brain coupling, we computed the shared neural dimension covariance after disrupting each factor. To disrupt intra-brain correlation, we used the inferred spikes from the calcium activity (produced in CNMF-E) and binned them into 300-frame bins (that is, 20-s). Spikes within each bin were then permuted across neurons^19^. This disrupted intra-brain correlation while maintaining average population activity. We then re-convolved the spikes with the MATLAB function conv with the default kernel used in the CNMF-E process. We then binned the data into 1 s bins. After this, we re-computed the shared neural dimensions and their covariance. To disrupt temporal coupling, we shifted the neural activity of one mouse relative to the other. We then re-computed the shared neural dimensions and their covariance. Covariance values from the original data, from disrupted intra-brain correlation, and from disrupted temporal coupling were each baseline-corrected by subtracting the covariance computed from a shuffle control. In this control, neural activity was first randomized across neurons at each time point, followed by neuron-wise circular shifts with randomly assigned time lags. Each permutation was repeated 100 times to compute an average covariance.

Contribution of state and state-transition to inter-brain coupling.

We defined five behavioural states, each of which covers a group of previously defined aggression (attack, chase, tussle, threaten), defence (escape, defend, flinch), social (general sniff, sniff face, sniff genital, approach, follow), non-social (dig, climb, explore object, bite object, stand), and idle (self groom, other). To estimate the fraction of PLSC1 or U1 variance explained by behavioural states alone, we fitted the exponentially filtered behavioural states matrix to PLSC1 or U1, respectively, using linear models. We defined state transitions as the moments when mice changed their behavioural states from one to another. To estimate the unique fraction of variance explained by behavioural state transitions, we fitted both the exponentially filtered state matrix and the exponentially filtered transition matrix to PLSC1 or U1 using linear models (full models). We then calculated the difference between the variance explained by full models and the variance explained by corresponding state-only models, which reflected the unique contribution of behavioural transitions.

Analysis of inter-animal behaviour coupling

High-dimensional behaviour coupling analysis using CCA.

To explore higher-dimensional behaviour coupling during social interaction, we used CCA using the canoncorr MATLAB function. This method introduces an additional step of matrix whitening on behaviour data before performing singular value decomposition. Unlike neural data, where each neuron carries a distinct physical meaning, behaviour features are arbitrarily defined, and their covariance predominantly reflects redundancy. Consequently, CCA was chosen to identify pairs of eigenvectors maximizing correlation among behaviour features across mice.

We applied the technique to identify shared dimensions relating the behaviour of the mice, producing two orthogonal subspaces. The first captures behaviour correlated across mice, termed the coordinated behaviour subspace, while the second captures behaviour not correlated to each individual mouse, termed the uncoordinated behaviour subspace.

The significance of these dimensions was assessed by comparing their correlation with the corresponding null distribution generated using temporally permuted data (2,000 randomly generated time lags). Dimensions exceeding the 97.5th percentile of the null distributions were considered significant. Similar to the shared neural subspace, only the top consecutive significant dimensions were counted as described in ‘Analysis of shared neural dynamics using PLSC’. The top coordinated behaviour dimension (CC1) refers to the most correlated pair of behaviour dimensions identified by CCA.

We note that correlation in behaviour features reflects temporal coordination in the behaviours of the mice, as opposed to identical movements. These correlations may arise from different features across mice. A more in-depth analysis and description of the methods used for this is discussed in ‘Behaviour symmetry analysis’.

Behaviour symmetry analysis

Within animals.

Because we defined the behavioural features, they may be highly correlated. In our CCA analysis, we sought to identify how uncorrelated behavioural features are shared. We therefore identified orthogonal behaviour dimensions by applying PCA to the behavioural data. We then computed CCA on the PCA-transformed data. Although CCA on the PCA-transformed data will yield the same correlation as on the raw data, we performed PCA to identify whether the weights of orthogonal behavioural features (principal components) were significant or not. We defined significant features in CC1 and CC2 as those with absolute weights at least 1.5× s.d. above the weight distribution for either of the dimensions. We term the vector of weights of these significant features the ‘loading vectors’. We then computed the cosine similarity of the loading vectors for CC1 and CC2. To control for chance, we temporally shifted the transformed behaviour spaces of two mice relative to each other and constructed a null distribution of cosine similarity values. Observed values were compared to percentiles (95th) to determine whether CC1 and CC2 dimensions emerged from significantly similar features.

Across animals.

To explore whether behaviour synchrony mostly comprises identical movement or temporally coordinated actions, we concatenated the behaviour spaces of both mice within a pair. We then orthogonalized using PCA to remove feature redundancy and provide shared coordinates for projection. Coordinated dimensions were re-computed using these principal components, and significant features from CC1 were identified. The cosine similarity of loading vectors from both mice was measured to determine if CC1 emerges from similar features across mice by comparing the observed value to the 95th percentile of the null distribution.

Temporal lags of coordinated behaviour across interacting animals.

We computed the correlation of the top coordinated behaviour dimension (CC1) across lags from −30 to +30 s, using 100 ms intervals, to assess whether temporal lags exist between behaviours of the two interacting mice. The lag represents the timing of the aggressor’s behaviour relative to the submissor’s within each pair. A positive lag indicates that the aggressor’s behaviour precedes that of the submissor mouse.

We measured the cross-correlation of CC1 between aggressor and submissor mice. When examining the timing of behaviours in these mouse pairs, we observed a distinct peak at 0-s lag, with no significant difference in correlation within the 1-s window before and after (Extended Data Fig. 11k,l). This indicates that coordinated behaviours between interacting mice also occur without an observable lag in our setup, though a small lag may exist at a timescale much faster than what our calcium imaging and video recording rate (at 15–30 fps) can detect.

Analysis of the relationship between neural activity and behaviour

Modelling a linear relationship between annotated behaviours and the top shared neural dimension.

To explore the behavioural context of the top shared neural dimension (PLSC1), we smoothed the binary annotated behaviours with an exponential kernel and z-scored them. These were used as predictors in a GLM model with a normal distribution to model PLSC1 activity. We then obtained the absolute coefficients of the annotated behaviours to examine their contribution to PLSC1 (Fig. 2s, Extended Data Fig. 8a).

Variance explained between multi-dimensional behavioural and neural spaces

Modelling a linear relationship between the high-dimensional behavioural and neural space.

To address the multicollinearity and high-dimensional nature of our behaviour and neural space, we performed partial least squares regression (PLSR)^24^, a latent variable modelling technique, using the plsregress function in MATLAB. In brief, PLSR simultaneously decomposes the input and output spaces into different sets of latent variables that maximize their covariance. To identify significance, we computed the neural variance explained by the observed model and subtracted the chance neural variance obtained by temporally shifting the data using 1,000 randomly generated time lags.

Calculating neural variance explained by high-dimensional behavioural space in various neural spaces.

For each of the neural spaces of interest (shared neural subspace, unique neural subspace, or full neural space), we combined the high-dimensional behaviour space of the subject and partner as predictors and modelled the neural space activity using PLSR. We then subtracted the chance level, calculated by temporally shuffling the predictors relative to the neural activity using 1,000 randomly generated time lags.

Calculating non-redundant neural variance explained.

To calculate the non-redundant variance explained by the model, we initially computed the full model variance and subtracted the chance level^47^, described above. Next, we adjusted the model by introducing a new set of predictors where the variables of interest were temporally shifted relative to the rest. This disrupts the temporal relationship of the variables of interest. The resulting new model variance, which accounts for chance, was then subtracted from the full model variance to determine the non-redundant variance. This process yields a lower bound for the variance captured by the variables of interest.

To compute the non-redundant neural variance of social and non-social behaviours in the shared neural subspace, we applied an exponential decay to binary vectors indicating social and non-social behaviours. We fit all social and non-social behaviours to the shared neural subspace, quantifying the variance explained by the model and subtracting chance-level variance. To quantify the non-redundant neural variance of social behaviours, we temporally shifted the inputs of the PLSR (all social behaviours) relative to the outputs (the neural space), refit the model, and quantified the neural variance explained by this permuted model. The chancel-level variance explained is estimated by the average of the variance explained by 100 permuted models. We then subtracted this chance-level variance from the original model variance. The residual variance in the shared neural subspace therefore reflects variance explained non-redundantly by social behaviours. The same process was repeated for non-social behaviours. We also repeated the process to assess the contribution of social and non-social behaviours in the unique neural subspace. Given that the unique neural subspace is generally larger than the shared neural subspace, we limit the unique subspace dimensionality by only using the top principal components in the unique neural subspace that account for approximately the same amount of neural variance as the shared neural subspace.

To compute the non-redundant neural variance explained by coordinated behaviours, the predictors were the coordinated behaviours, self-uncoordinated behaviours, and partner-uncoordinated behaviours. To quantify the non-redundant contribution of coordinated behaviours to the shared neural subspace, we computed the decrease in variance explained after permuting the coordinated behaviours. To compute the non-redundant neural variance explained by the self and partner’s uncoordinated behaviours, we used the self and partner’s original behaviour space (unpartitioned) as predictors for the full model. This allowed all possible redundancy to be accounted for without underestimating the uncoordinated behaviours. To quantify the non-redundant contribution of self/partner behaviours to the shared neural subspace, we computed the decrease in variance after permuting the unpartitioned behaviour space of self/partner (while keeping other predictors unchanged). We presented the variance explained by coordinated behaviours as well as uncoordinated behaviours of self and partner as a ratio to the total variance possibly explained by both self and partner behaviours in the relevant neural space.

To compute the non-redundant neural variance explained by the partner’s behaviours in AI agents, the predictors comprised collision events, x and y positions of self and partner, actions (up, down, left and right) of self and partner, and behaviours of self and partner. These behaviours were new fields for both agents, approach for chasers (indicating whether the chaser made a move to reduce the distance from the last time step), and escape for explorers (further details explained in ‘Task design’). For the artificial neural activity, we projected the RNN activation (ht) of each agent into their ‘neural action subspace’ by using the trained output weight (Waction), computing WactionTWactionht. We then fit the predictors to the neural action subspace activity in the full model, and assessed the non-redundant contribution of the partner’s behaviours by permuting them and measuring the resulting decrease in variance explained. Note that because we fit a linear model, we did not require WactionTWaction to be a projection matrix.

Modelling tracking behavioural features with neural activity.

We assessed the representation of each of the 75 tracking behavioural features in the neural space. In brief, we fit the full neural activity to individual tracking features using a GLM and computed the variance explained of the model. The chance-level variance was then subtracted to obtain the variance explained of behavioural features by the neural activity. In order to estimate the chance-level variance, we temporally shifted the neural activity relative to the behavioural feature and refit the GLM. The chance-level variance explained was averaged across 100 random permutations.

Multi-agent reinforcement learning

Task design.

We modelled competitive social interactions in rodents within a 10 × 10 grid world using RLlib^28^. In this task, there was a chaser and explorer agent, each equipped with a 7 × 7 field of vision. We modified the rewards of this task to define both a social and a non-social task. Agents in the non-social task did not have their partner’s location as input. This non-social task acted as a control for our analyses. The rewards are described in Supplementary Table 3. Note that in the non-social task, the reward for collision is the same as the reward received every time step if no new fields or collision occurs.

In addition to ‘collision’ and ‘new fields explored’, we defined event types for subsequent neural-behavioural analysis. For the chaser, we defined an approach event as when the agent takes an action that reduces its distance to the explorer’s last position. For the explorer, we defined an escape event as when the agent takes an action that increases its distance from the chaser’s last position. Within escape events, we further defined ‘escape far’ when the updated distance between agents is more than five units, ‘escape near’ when the updated distance is three to five units, and ‘escape close’ when the updated distance is less than three units. These escape events are mutually exclusive. The initial locations of agents were randomly initialized at the beginning of each episode.

In the analysis examining whether asymmetric partner representations emerged independently from reward definitions (Extended Data Fig. 13f), we modified the following rewards while keeping the rest of the task the same (Supplementary Table 4). New agents were trained using these rewards. This result was consistent across environments with different reward values, suggesting that this asymmetric relationship was not due to specific values of the rewards, but was inherent to the distinct role of chasers and explorers.

Agent architecture.

The input to each agent was the concatenation of two vectors. Each vector had a length of 100, corresponding to the 100 tiles in the 10 × 10 grid world. One vector was the one-hot location of the agent (self). The other vector was either the one-hot location of the opponent if it was within the agent’s vision, or else the zero vector.

We used a recurrent actor–critic architecture for our agents. The inputs were fed into a vanilla RNN with 256 hidden units using the ReLU activation. The RNN units were mapped to two parallel linear layers to compute the action logits (actor) and the state value (critic). The architecture is,

RNN: ht=ReLUWinputxt+binput+Wrecht-1+brec, where Winput∈R256×200,xt∈R200×1,Wrec∈R256×256,ht∈R256×1,binput∈R256×1,brec∈R256×1, and ht=ReLU(xt). The value readout is defined as vt=Wvalueht+bvalue, where Wvalue∈R1×256,ht∈R256×1,bvalue∈R,vt∈R.

The action readout is defined as at=Wactionht+baction, where Waction∈R4×256,ht∈R256×1,baction∈R4×1,at∈R4×1.

Agent training and evaluation.

We trained agents with PPO^29^ for 20,000 epochs, each epoch consisting of 4,000 environment steps sampled from 40 full episodes (each of length 100 timesteps). This training length allowed the agents’ performance against a random agent to plateau. We used PPO with a Kullback–Leibler divergence penalty on policy updates using default parameters from RLlib. We also applied L2 regularization to the model weights with loss function weight λ = 0.3. Other parameters were set to RLlib defaults. For each task (social versus non-social), we trained a total of ten pairs of agents from different initial seeds.

To analyse the relationship between neural activation and agent behaviours, we evaluated each agent pair in 25 episodes, each lasting 500 timesteps. Because MARL is a challenging learning setting where, from the perspective of one agent, the environment is constantly changing, we observed that it was at times possible for agents to exhibit degenerate behaviour. We therefore excluded episodes where agents ceased moving in the same state or repeated movements within the same tile for more than 1% of the episode length.

To assess agents’ performance at various training iterations, we paired each trained agent with a standardized random agent. The random agent’s actions at any timestep were randomly sampled from a uniform distribution. These agent pairs engaged in 100 episodes, each spanning 100 timesteps. This ensured that agents were evaluated against a common opponent they were not trained with.

For chaser agents, we evaluated criteria including the number of collisions, the percentage of time agents kept their partner in vision, the average distance between agents, and the number of new fields. The first three criteria reflect how well the chaser exhibited chasing behaviour. For explorer agents, we evaluated the number of new fields and the average distance between agents. These reflect how well the explorer simultaneously explored new tiles while maintaining distance from the chaser.

Quantification of action flow field and angle.

To characterize the actions of trained agents at the conclusion of their training, we examined their overall movement with respect to their opponent’s relative position. We calculated the angle between the agent’s action (up, down, left and right) and the vector pointing from the agent to the opponent in an egocentric perspective. We only performed this analysis when the opponent was in the agent’s field of vision.

We generated a matrix representing this angle across their relative x and y distances and averaged the results across all 10 agent pairs. This averaged matrix was then used to construct the flow field, polar plots and bar plots reflecting the agent’s overall behavioural strategy.

Decoding social events and partner behaviour in MARL.

In order to evaluate the partner (opponent) representation in the artificial neural network, we used the RNN activity to decode collision events or partner behaviours including approaching and escaping. For decoding chasing or approaching, we fit RNN activity during time steps when the partner was in vision. The performance of the SVM was assessed using balanced accuracy, Balanced accuracy=12TPTP+FN+TNTN+FP, where TP, TN, FP, and FN stand for true positive, true negative, false positive and false negative, respectively. For comparison, we also fit SVM classifiers using the RNN activity of non-social agents.

Shared neural space in MARL.

To examine whether inter-agent neural coupling emerged in social agents, we performed the analyses described in ‘Analysis of shared neural dynamics using PLSC’ for MARL agents. To explore whether the shared neural subspace emerged solely from shared visual input, we provided the non-social agents either with partial vision (7 × 7) or full vision of the partner and re-computed the shared neural subspace.

MARL partner representation in neural action space.

To assess how partner representation impacts agents’ actions across training, we investigated the extent to which neural variance in the neural action space (the 4D neural subspace spanned by WactionTWaction) could be attributed to partner behaviours after subtracting the variance already accounted for by self behaviours. For each agent, we projected the neural activation of the RNN onto the neural action space using the trained weights by computing WactionTWactionht. Subsequently, we used the technique described in ‘Calculating non-redundant neural variance explained’, to compute the unique contribution of partner behaviours to the agents’ neural activity.

MARL perturbation.

To perform perturbations, we concatenated all rollout episodes and identified the top 10 PLSCs for each agent pair. We then projected out the top 10 PLSCs by projecting the RNN activation for each chaser agent into the null space (I − PP^T^) of the PLSCs. We used the null space activity to compute the agent’s action logits in new online episodes. This method, therefore, evaluates how agents interact when PLSC activity is removed. We evaluated agent performance by quantifying the number of collisions, the percentage time they kept the partner in vision, and the average distance between agents across 100 episodes, each 100 timesteps. As a control, we also regressed out the top 25 random principal components that did not overlap with the top 10 PLSCs. To do this, we computed the null space of the top 10 PLSCs and then projected the RNN activation into this null space space. We then temporally permuted the activation of all units with different time lags to disrupt structure, and performed PCA on this permuted space. These 25 random principal components account for approximately the same amount of variance as the top 10 PLSCs. Finally, we projected the RNN activation onto the null space of the top 25 principal components prior to computing the agent’s action logits. We also collected 100 episodes for this set of agents.

Additional MARL tasks

To examine how agents’ goals may influence neural correlations across interacting agents, we examined three different tasks with varying collaborative and competitive goals and complexities (Extended Data Fig. 13).

Mutual interaction task.

We used the same 10 × 10 grid world environment and 7 × 7 field of vision (Extended Data Fig. 13h) as in the chaser-explorer task. Rewards were defined such that cooperative mutual interaction between agents resulted in a positive reward (+1) for each agent, while non-social exploration received a small positive reward (+0.1) and all other events resulted in a small baseline penalty (−0.1). Training and evaluation were carried out using the same procedures outlined for the chaser-explorer task. Agents were trained with PPO for 20,000 epochs, each comprising 40 episodes with 100 timesteps per episode. 10 pairs of agents were trained. Agent pairs were evaluated in 25 episodes, each with 500 timesteps.

Game theory-based coins task.

The Coins task is a multi-agent RL environment that can be viewed as a higher-dimensional spatially extended version of the iterative Prisoner’s Dilemma involving both cooperation and competition^48^. A blue and a red agent are placed in a 3 × 3 environment where they can pick up blue and red coins (Extended Data Fig. 13m). Agents can pick up coins of any colour. Anytime an agent picks up a coin, it receives a +1 reward. However, picking up a coin of the opposing colour penalizes the other agent with a −2 reward. By analogy to the Prisoner’s Dilemma, if both agents pick up each other’s coins, they both receive negative rewards. Agents can also exploit the other agent by picking up their coins, even as the other agent picks up coins of their own colour. Finally, if both agents pick up their own coins, they will both receive positive rewards. Agents may therefore display varying levels of cooperation (picking up their own coins and not their partner’s coins) or selfishness (picking up both their and the partner’s coins). The rewards are described in Supplementary Table 5.

Training cooperating agents in the Coins task is empirically challenging. Agents trained naively with PPO do not cooperate above chance level, but pick up any available coin^49,50^. A cooperating agent should learn that picking up the partner’s coin has less state value than its own coin—even though the agent receives the same +1 reward for picking it up—because its actions will negatively impact its partner. This negative impact on its partner will modify the partner’s policy, and may lead to a negative influence on itself. This presents a technical agents must learn cooperating policies that may or may not lead to future negative rewards, depending on whether their partner also learns a cooperating policy. Because an agent cannot directly control its partner’s policy, we emphasize that this learning has to occur in both agents to achieve cooperation and avoid exploitation.

To train agents, we used a gated recurrent unit network instead of a vanilla RNN (original and MI task), which helps agents retain task observations and memory over longer timescales, facilitating training^51^. Using a gated recurrent unit network also enables the evaluation of whether shared dynamics emerge in a different RNN architecture.

Complex social task.

We incorporated a more complex social environment using the DeepMind Melting Pot environment^52^. which incorporates more complex visual inputs, dynamic internal states, and a more diverse behavioural space. In brief, in this task, two agents have distinct the gatherer earns rewards by retrieving and collecting objects (apples and acorns), while the capturer earns rewards by intercepting the gatherer (Extended Data Fig. 13p). Gatherers receive an immediate +1 reward for collecting an apple. Gatherers may receive up to +18 delayed rewards for collecting an acorn by ‘eating’ the acorn, a process that involves additional complexity. To eat an acorn, the gatherer must perform an ‘interact’ action to eat 1/3 of the acorn, for which it receives +6 reward. During eating, the gatherer cannot move for 5 timesteps, increasing its susceptibility to capture. Gatherers are provided a 3 × 3 safe zone, where capturers cannot enter; gatherers can take advantage of this safe zone by retrieving acorns to this area before consuming it. Each agent has a dynamic amount of stamina and can choose from one of the eight possible actions at each time no-operation, move forward, move backward, move left, move right, turn left, turn right, and interact. Each agent’s stamina capacity is 18 units. Every action, except no-operation, reduces stamina, while no-operation increases stamina. Eating an acorn always reduces 6 stamina units and freezes the explorer for 5 timesteps. For all other actions, the amount of stamina reduced or increased depends on the current stamina when stamina is between (a) 7−18 units, (b) 1–6 units, or (c) 0 units, an action will reduce or increase stamina by (a) 2 units, (b) 3 units, and (c) 5 units in the gatherer, and (a) 1 unit, (b) 2 units, and (c) 7 units in the capturer. Further, each action induces freezing for (a) 1 timestep, (b) 2 timesteps, or (c) 4 timesteps in the gatherer and (a) 0 timesteps, (b) 1 timestep, or (c) 6 timesteps in the capturer.

We adopted previously trained actor–critic baseline agents^52^. Each agent has a field of view of 5 to the left, 5 to the right, 9 forward, and 1 backward. This vision input, represented in RGB pixel values, was processed through a convolutional neural network (CNN). The CNN contained 2 convolutional layers followed by two fully connected layers, providing features to a long-short term memory network (LSTM) with 128 units. The input to the LSTM was, therefore, CNN-derived visual features. The LSTM then outputs the policy and value function for actor–critic training. Agents were trained using IMPALA^53^ with a contrastive predictive coding auxiliary objective. Further details of the environment, training, and network architecture were previously described^52^. We chose to analyse these agents because they exhibited generalized and complex behaviours as a result of large-scale training across the Melting Pot suite. Because our goal was to evaluate whether shared neural dimensions emerge in more complex MARL agents, we selected 5 well-trained gatherers and capturers and performed new rollouts with pairs of gatherers and capturers. Gatherers and capturers displayed complex apple and acorn gathering, as well as capturing strategies, in our rollout environment. Finally, analysing the LSTM enables the evaluation of whether shared dynamics also emerge in another RNN architecture. The rewards are described in Supplementary Table 6.

Statistics and reproducibility

All statistical analyses were conducted using Prism (v10, GraphPad), MATLAB (R2022a, MathWorks) or Python (3.10). Details about the types of statistical tests used, test statistics, P values, and sample sizes are provided in Supplementary Table 1. When parametric tests were used, data normality was confirmed using the D’Agostino–Pearson test. P values were corrected for multiple comparisons when necessary. Bar plots show mean ± s.e.m. In box plots, the centre line indicates the median, the box limits indicate the upper and lower quartiles (IQR), and the whiskers indicate the minimum and maximum values or 1.5× IQR. The significance threshold was held at α = 0.05 (NS, not significant (P > 0.05); *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). All experiments were replicated in several subject mice with similar results. Sample sizes were not predetermined using statistical methods. Experiments were randomized whenever possible. Experimenters were not blind to group allocation.

Extended Data

Supplementary Material

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41586-025-09196-4.