Inter-brain neural dynamics in biological and artificial intelligence systems

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][3][4][5][6][7][8] . Experiments using in vivo calcium imaging and electrophysiology 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][12][13][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][16][17][18][19][20][21][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 showing behaviours of two mice in an interaction session. h,i, Example traces of three individual glutamatergic (h) or GABAergic (i) neurons (top) and aggregate (mean) activity of all recorded neurons in a mouse (bottom) in an interaction session. j,m, Inter-brain correlation of mean activity in dmPFC glutamatergic ( j) or GABAergic (m) neurons during interaction sessions compared with the corresponding shuffled data. k,n, Cross-correlation of mean activity of glutamatergic (k) or GABAergic (n) neurons across mice during interaction sessions, compared with separation sessions and phaserandomized signals. l,o, Inter-brain correlation of glutamatergic (l) or GABAergic (o) neurons across different timescales, compared with phaserandomized signals. The timescales are the period (1/frequency) of the frequency bands used for the fast Fourier transform (FFT) decomposition of mean activity. p, Fraction of social, non-social, mixed or neutral responsive neurons in the GABAergic or glutamatergic population. q-t, Performance of a SVM decoder trained to distinguish between pairs of behaviours using population activity from each cell type comparing attack versus sniff (q), attack versus approach (r), attack versus self groom (s) and sniff versus self groom. u,v, Inter-brain correlation after computational removal of socialresponsive, non-social-responsive, and random neurons from the GABAergic (u) or glutamatergic (v) population. k,l,n,o, Data are mean ± s.e.m. In box plots, the centre line is the median, box edges delineate interquartile ranges (IQRs) and whiskers extend from minimum to maximum values (j,m,q-t) or 1.5× IQR (u,v). See Supplementary Table 1 for details of statistical analyses. ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05; NS, not significant (P > 0.05).

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

Article

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

-0.3 -0.1 0.1 0.3 -0.4 -0.2 0 0.2 0.4 0 10 20 30 0.4 0.5 0.6 0.7 0.8 0.9 0 0.02 0.04 0.06 0.08 0.10 Probability PLSC1 correlation -0.2 0 0.2 0.4 0.6 0.8 1.0 **** NS -0.2 0 0.2 0.4 0.6 0.8 1.0 0 10 20 0 10 20 *** *** 30 0 10 20 **** 30 Neural variance explained (%) 0 20 40 60 80 100 Neural variance explained (%) 0 10 20 30 40 *** 50 60 Nonredundant neural variance explained (%) 0 10 20 30 **** 40 50 Nonredundant neural variance explained (%) 50 100 150 200 250 300 350 5 4 3 2 1 0.06 0 -0.09 0.07 0.11 50 100 150 200 250 300 350 5 4 3 2 1 0.94 0.84 0.83 0.73 0.78 R 2 = 0.681, P = 0.0003 Unique subspace Shared subspace Covariance matrix Mouse 2 Animal 2 neurons Animal 1 neurons Mouse 1 SVD a c Mouse 1 Mouse 2 Shared neural subspace Unique neural subspace Unique neural subspace p o m n q r s t PLSC1 PLSC2 Both b Attack Escape Defend General sniffs Tussling Chasing Face sniffs Climb Dig Stand Approach Self groom Follow Biteobj Flinch Exploreobj Genital sniffs Threaten 0 0.5 1.0 1.5 Absolute coefficients S o c ia l N o n -s o c ia l S o c ia l N o n -s o c ia l Glutamatergic GABAergic ** Inter Inter Sep Sep 0 20 40 60 80 Number of PCs NS **** **** High Low Mutual social time Unique dimensions Shared dimensions 100 s k l 0 2 4 6 0 20 40 60 80 GABAergic Unique subspace Shared subspace Glutamatergic 30% 70% 5% 95% d Time (s) Time (s) e f |WPLSC1| -|WU1| |WPLSC1| -|WPLSC2| Cell count 0 10 20 30 Cell count Approach Sniffing Attack Chase Escape Defend Self groom Dig Shared dimensions (PLSC) Unique dimensions (U) g Glutamatergic neurons GABAergic neurons j i h Percentage of attack behaviour PCC PCC Mouse 1 Mouse 2 Mouse 1 Mouse 1 behaviour Mouse 2 behaviour Mouse 1 behaviour Mouse 2 behaviour Mouse 2 Number of shared dimensions S e p a r a t io n s e s s io n In t e r a c t io n s e s s io n S e p a r a t io n s e s s io n In t e r a c t io n s e s s io n G lu t a m a t e r g ic G A B A e r g ic ΔPCC of PLSC1 ΔPCC of U1 G lu t a m a t e r g ic G A B A e r g ic G lu t a m a t e r g ic G A B A e r g ic Neural variance explained by the shared space (%) Mouse 1 Mouse 2 Shared subspace Unique subspace Mouse 1 Mouse 2 = WXΔW T Y Projection Projection Time Time Shared : WY Shared : WX Unique : PCs of ((W ⊥ Y ) T Y) Unique : PCs of ((W ⊥ X ) T X) WX Neuron 1 Neuron 2 Neuron 3 Neuron 1 Neuron 2 Neuron 3 WY WY WX PLSC1 U1 Both 97.5% Observation

Fig.

2 | Shared and unique neural subspace in socially interacting mice. a, Schematic showing partitioning neural activity of interacting mice into shared and unique neural subspaces. Created in BioRender. Phi, N. (2025) https://BioRender.com/fxtjeqz. b, PLSC is used to identify shared neural dimensions that maximize cross-covariance between two mice. The unique neural subspace is derived from non-significant dimensions. SVD, singular value decomposition. Variables are defined in the Methods. c, Example shared and unique neural dimensions. d, Example null distribution of correlation for the top shared neural dimension (PLSC1) using temporally permuted data. e,f, Example traces of shared (PLSC) (e) and unique (U) (f) neural dimensions, aligned with the behaviours of both mice. Right, PCC for each PLSC or U. g,h, Number of significant PLSCs in GABAergic (g) or glutamatergic (h) neurons during interaction and separation sessions. i,j, ΔPCC of the top shared (i) or unique ( j) neural dimensions in GABAergic or glutamatergic neurons (that is, the PCC of observed data minus the PCC of shuffled controls). U1, PC1 of unique neural subspace. k,l, Distributions of significant GABAergic neurons in the top shared versus unique neural dimensions (k) or the first versus second shared neural dimension (l), based on their weights. W U1 , W PLSC1 and W PLSC2 correspond to loadings of neurons in U1, PLSC1 and PLSC2, respectively. m, Number of principal components that capture 50% neural variance during interaction (Inter) or separation (Sep) sessions. n, Percentage of neural variance explained by the shared and unique neural subspace in GABAergic or glutamatergic neurons during interaction sessions, averaged across mice. o, Percentage of neural variance explained by the top shared neural dimension (PLSC1) in each mouse. p,q, Percentage of non-redundant neural variance explained in either the shared (p) or unique (q) neural subspace by social and non-social behaviours in GABAergic neurons. r, The shared neural subspace explains more variance in mice showing more mutual social interaction in GABAergic neurons. s, Contribution of annotated behaviour to the top shared neural dimension PLSC1 (Methods) in GABAergic neurons. Data are mean ± s.e.m. Biteobj, bite object; Exploreobj, explore object. t, Percentage of time showing attack positively correlates with the size of shared neural subspace in GABAergic neurons across mouse pairs. Each data point represents an average of two mice within a pair. In box plots, the centre line is the median, box edges delineate IQRs and whiskers extend from minimum to maximum values (g-j) or 1.5× IQR (m,o-r). See Supplementary Table 1 for details of statistical analyses.

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).

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

-10 0 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 0 10 -10 -5 0 5 10 0 5 10 15 20 **** ** c g PC1 **** **** ** ** R 2 = 0.

32, P = 0.037 Ensembles identified in interaction sessions **** **** Ensembles identified in separation sessions Interaction session Separation session b a e j t-SNE 1 t-SNE 2 t-SNE 1 t-SNE 1 t-SNE 1 Ensembles (identified in interaction sessions) Ensembles (identified in separation sessions) During interaction session During separation session i h p o Perturbing temporal coupling but maintaining intra-brain correlation Perturbing intra-brain correlation but maintaining temporal coupling Temporal coupling Intra-brain correlation Inter-brain coupling + Neurons Neurons Neurons Neurons 0 10 20 30 Interaction Separation O b s e r v e d S h u f fl e d O b s e r v e d S h u f fl e d Interaction Separation O b s e r v e d S h u f fl e d O b s e r v e d S h u f fl e d 0 0.1 0.2 0.3 0.4 PCC (cells, PC1) 0.1 0.2 0.3 0.2 0.4 0.6 PCC (cells, PC1) 0.1 0.2 -0.05 0 0.05 0.10 PCC (cells, PC1) R 2 = 0.19, P = 0.11 Social moments in interaction sessions R 2 = 0.56, P = 0.002 R 2 = 0.005, P = 0.83 0.1 0.2 0.3 0.2 0.4 0.6 PCC (cells, PC1) Non-social moments in interaction sessions 0.1 0.2 0.3 0.2 0.4 0.6 PCC (cells, PC1) l 0 0.1 0.2 0.3 PCC (pairwise) 0 0.05 0.10 0.15 0.20 0.25 PCC (pairwise) 0 0.1 0.2 0.3 0.4 Ensembles identified in social moments Ensembles identified in non-social moments 0 0.1 0.2 0.3 t-SNE 2 D u r in g in t e r a c t io n D u r in g s e p a r a t io n Social Non-social Mixed Ensembles (social) Ensembles (non-social) 63% 15% 22% 39% 20% 41% ΔPCC of PLSC1 ΔPCC of PLSC1 ΔPCC of PLSC1 ΔPCC of PLSC1 O r ig b, Intra-brain neural correlations, measured by average PCC between individual GABAergic neurons and PC1 within a mouse, during interaction or separation sessions (subsampled to 50 neurons). c-f, Scatter plots of intra-brain correlation versus inter-brain coupling (difference in PCC of observed PLSC1 and that of shuffled controls) in GABAergic neurons during the interaction session (c), separation session (d), social moments in interaction sessions (e) and non-social moments in interaction sessions (f). g,h, Example t-distributed stochastic neighbour embedding (t-SNE) plot showing the loadings of the GABAergic neuron ensembles that are identified in the interaction (g) or separation (h) session, during interaction (left) or separation sessions (right). Coloured dots represent single neurons in the intra-brain neural ensembles identified in interaction (g) or separation (h) sessions. i, Example raster of the ensemble neural activity of individual ensembles during interaction sessions. j,k, Average pairwise PCC between neurons within the same GABAergic ensembles that are identified in interaction ( j) or separation (k) sessions, during interaction or separation sessions. l,m, Average pairwise PCC between neurons within the same GABAergic ensembles that are identified in social (l) or non-social (m) moments in interaction sessions, during social or non-social moments. n, Behaviour tuning of GABAergic neuron ensembles identified in social or non-social moments within interaction sessions. o, Methods used to disrupt intra-brain correlation or temporal coupling. To disrupt temporal coupling, the activity of one mouse was shuffled with respect to that of the other. To disrupt intra-brain correlation, spikes were randomly shuffled across neurons of the same mouse within each time bin prior to being re-convolved into calcium traces. p, Covariance of PLSC1 in GABAergic neurons between interacting mice after disrupting intra-brain coupling or temporal coupling, relative to the shuffle baseline. j-m, Data are mean ± s.e.m. In box plots, the centre line is the median, box edges delineate IQRs and whiskers extend to 1.5× IQR (a,b,p). See Supplementary Table 1 for details of statistical analyses.

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).

Article

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

40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 6.5 7.0 7.5 8.0 Training iteration Average distance to opponent (units) 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 4 6 8 2 Training iteration Average distance to opponent (units) 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000

40 60 80 100 20 Training iteration 0 2 4 6 8 10 0 2 4 6 8 0 20 40 60 80 0 20 40 60 80 100 Exploration Collision +1 +0.1 In vision Out of vision +1 -1 Behaviour Rewards a b e f g h i j k l m n Chaser Explorer Chaser Chaser Chaser Explorer c d Social agents Non-social agents **** Relative x distance Relative x distance o p 0 0.1 0.2 Early training **** **** **** **** **** * Late training r q

Early stage of training Late stage of training 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 0 10 20 Training iteration Rewards 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 0 20 40 Training iteration Reward 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000

40 60 80 100 20 Training iteration 0 5 10 15 20 25 S o c ia l N o n -s o c ia l S o c ia l N o n -s o c ia l S o c ia l N o n -s o c ia l 0 20 40 60 80 100 Social environment Non-social environment Chaser Explorer Chaser Chaser Explorer Chaser 1 -3 0 3 1 -6 0 6 1 -9 0 9 1 -1 2 0 1 2 1 -1 5 0 1 5 1 -1 8 0 0 0.2 0.4 0.6 Angle (degrees) Probability Social environment Non-social environment

Social agents Non-social agents Social agents Non-social agents 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 0 10 20 Training iteration 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 40 60 80 20 Training iteration Number of new fields 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 40 60 80 20 Training iteration Number of new fields 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 5 10 15 20 0 Training iteration Number of collisions 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 40 60 80 20 Training iteration Number of new fields 0 20 40 60 80 100 250 500 750 1,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 40 60 80 20 Training iteration Number of new fields Number of collisions Time partner in vision (%) Time partner in vision (%) Sampling RNN Action distribution Chaser Winput Waction Sampling RNN Action distribution Explorer Winput Waction Chaser Explorer Environment Relative y distance Multi-agent reinforcement learning P(up) P(down) P(left) P(right) Action (e.g. up) Action (e.g. left)

Chaser and explorer position (in vision) new fields explored by chaser (g) or explorer (h). i-n, Evaluation of artificial agents, trained in social or non-social tasks, when acting with an agent taking random actions. The chaser's collision (i), the percentage of time the chaser keeps the explorer in vision ( j), chaser's (k) or explorer's (n) average distance to the partner, and new fields explored by chaser (l) or explorer (m). All metrics were evaluated using the first 100 timesteps in rollouts. o,p, Flow field showing chasers' moving direction towards the explorer in vision during early (o) or late (p) stages of training. The blue dot represents the explorer's location, and x and y axes indicate the relative distance of the chaser to the explorer in x and y directions, respectively. q, Polar plot showing the chasers' movement towards the explorer in vision in early (grey) or late (orange) stages of training. r, Angles between chasers' moving direction and their orientation towards the partner during early (grey) or late (orange) stages of training. c-h,i-m,n (left), r, Data are mean ± s.e.m. In box plots, the centre line is the median, box edges delineate IQRs and whiskers extend to 1.5× IQR (i-n). See Supplementary Table 1 for details of statistical analyses. 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 approach (e) from population activity in trained explorers. f, Example shared neural dimensions in interacting agents in a social environment, aligned with their behaviours. g, Number of significantly shared neural dimensions in trained social agents, compared with trained non-social agents with and without vision of the partner. h, ΔPCC of the top shared neural dimensions in trained social agents, compared with trained non-social agents with and without vision of the opponent (that is, the PCC of observed data minus the PCC of shuffled controls). i,j, Distributions of significant neurons in the top shared (PLSC1) versus unique (U1) neural dimensions in trained social agents (i) or the first versus second shared neural dimensions ( j), based on the absolute weights within these dimensions. k, Cross-covariance of PLSC1 in social agents after computationally disrupting intra-brain correlation or temporal coupling. l,m, Neural variance in trained social and non-social agents explained by partner's behaviours, after discounting self behaviours, in the full neural space (l) or in the top shared neural dimension (PLSC1) (m). n,o, Non-redundant neural variance explained by partner's behaviours in opposing agents (n) or mice (o). p, Partner representation (neural variance explained by partner) in the neural action space of social or non-social chasers, which is the four-dimensional neural subspace spanned by the RNN output matrix (Methods). q,r, Scatter plots showing the relationship between partner representation and the percentage of time with collisions in social (q) or non-social (r) agents. s-u, Removing activities in the top ten shared neural dimensions in chasers reduces collision events (s) and the percentage of time chasers keep the explorer in vision (t) while increasing the average distance between chasers and explorers (u), compared to removing a similar amount of random non-overlapping variance. n,o, Data are mean ± s.e.m. In box plots, the centre line is the median, box edges delineate IQRs and whiskers extend to 1.5× IQR (b-e,g,h,k-m,s-u). See Supplementary Table 1 for details of statistical analyses.

Article

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][31][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][3][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 structures: (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).

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.

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.

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 prewhitened 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. 1km,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 subspaces: 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 W X and W Y , 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 by:

is a z-scored matrix containing the neural activity of mouse 1 with n 1 neurons across the length of the entire experimental session, T. Each neuron has zero mean and a s.d. of 1. Similarly,

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 satisfies:

is a cross-covariance matrix, Δ represents the diagonal matrix of singular values (in descending order), and

so that the columns are orthonormal. We quantified the cross-correlation of the ith PLSC dimension by computing corr( , )

, where W X i , is the ith column of the matrix W X . 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

, 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.

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.

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 behaviours: 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.

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 ( t h ) of each agent into their 'neural action subspace' by using the trained output weight (W action ), computing t action T action W W h . 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 W W action T action 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:

The value readout is defined as

W h b a

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,

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.

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 challenge: 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 objectives: 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 step: 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 level: 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. 1c. e, Interbrain correlation of aggregate (mean) activity in dmPFC GABAergic neurons during interaction sessions compared to separation sessions and cross-pair controls. f, Cross-correlation of aggregate (mean) activity of GABAergic neurons during interaction and separation sessions, and the phase-randomized neural activity of interaction sessions. This panel shows a replot of the data from Fig. 1n using a smaller lag window. g, Inter-brain correlation of aggregate (mean) activity in dmPFC glutamatergic neurons during interaction sessions compared to separation sessions and cross-pair controls. h, Cross-correlation of aggregate (mean) activity of glutamatergic neurons in interaction and separation sessions, and the phase-randomized neural activity of interaction sessions. This panel shows a replot of the data from Fig. 1k using a smaller lag window. i, j, Example trace (top) and autocorrelation (bottom) of aggregate (mean) activity (i) and phase randomized control (j). k-m, Inter-brain correlation (PCC) of individual components from the Fourier transform of the aggregate (mean) neural traces (Mean ± SEM). Data from interaction sessions shown in k are replotted in l and m alongside data from the corresponding separation sessions. n, Autocorrelation of an example GABAergic neuron activity before and after pre-whitening. o, Inter-brain correlation of pre-whitened GABAergic (left) and glutamatergic (right) neurons in interaction sessions, separation sessions, and shuffled controls. p, q, Cross-correlation of pre-whitened neural activity in GABAergic (p) or glutamatergic (q) neurons across animals during interaction sessions, compared to separation sessions and shuffled controls (Mean ± SEM). r-s, Inter-brain correlation (PCC) of individual components from the Fourier transform of the pre-whitened aggregate (mean) activity of GABAergic (r) or glutamatergic (s) neurons (Mean ± SEM). Of note, prewhitening partially removes real signals, potentially reducing the observed difference in inter-brain correlation between interaction and separation sessions. t, Firing rates in GABAergic and glutamatergic neurons. u, v, Firing rates (u) and inter-brain correlation (v) in different cell types after subsampling the top 20% glutamatergic neurons based on firing rates. w-x, Firing rates (w) and inter-brain correlation (x) in different cell types after subsampling the top 50% glutamatergic neurons and the bottom 50% GABAergic neurons based on firing rates. t-x, mean ± SEM. y, z, Inter-brain correlation (PCC) from only socialresponsive or non-social-responsive neurons in GABAergic (y) or glutamatergic (z) population. p**** < 0.0001, p*** < 0.001; p** < 0.01; p > 0.05, n.s. Boxplots: medians, interquartile ranges (IQR), and whiskers indicating min-max (y, z) or 1.5× IQR (e, g, o). See Supplementary Table 1 for details of statistical analyses.

Extended Data Fig. 2 | Encoding of social behaviors in GABAergic and glutamatergic neurons. a-e, Performance of SVM decoders trained to distinguish different behaviors from baseline (non-behavioral moments) using population activity from each cell type, with glutamatergic neurons downsampled to match the average number of GABAergic neurons. f-i, Performance of SVM decoders trained to distinguish pairs of behaviors using population activity with both neuronal populations down-sampled to 50 neurons per animal. j-n, Performance of SVM decoder trained to decode different behaviors from baseline (non-behavioral moments) using population activity with both neuronal populations down-sampled to 50 neurons per animal. Boxplots: medians, interquartile ranges (IQR), and whiskers indicating min-max (a-n). See Supplementary Table 1 for details of statistical analyses. 1 for details of statistical analyses. Extended Data Fig. 5 | Shared neural subspace during social and non-social moments. a-h Inter-brain correlation (PCC), number of significant shared neural dimensions, chance-level-subtracted PCC of PLSC1 (the PCC of observed data minus the PCC computed from temporally shuffled data), and fraction of variance explained by shared neural subspace during mutual social, unidirectional social, and non-social moments within interaction sessions in GABAergic (a-d) or glutamatergic neurons (e-h). Comparing mutual social conditions between GABAergic and glutamatergic neurons in a and e, p = 0.0003; b and f, p = 0.0644; c and g, p < 0.0001; d and h, p = 0.0008 (two-way ANOVA followed by Sidák multiple comparisons). i-l: Inter-brain correlation (PCC), chance-level subtracted PCC of PLSC1, fraction of variance explained by shared dimensions, and number of shared dimensions in GABAergic or glutamatergic neurons during mutual social versus concurrent rearing (mutual non-social) behaviors within interaction sessions. In i-l, all mutual social data are replotted from a-h. m, Higher neural variance explained by the glutamatergic shared neural subspace in animals that engaged in a higher level of mutual social interaction. Boxplots: medians, interquartile ranges (IQR), and whiskers indicating 1.5× IQR (a-m). See Supplementary Table 1 for details of statistical analyses. Behavioral categories are: attack includes attack, tussle, chase, and threaten; defense includes defend, escape, and flinch; sniff includes general-sniffs, face-sniffs, and genital-sniffs; approach includes approach and follow; non-social includes self-groom, dig, explore object, bite object, climb, and stand. l, m, Fraction of variance in GABAergic (l) or glutamatergic (m) neural dimensions PLSC1 or U1 explained by behavioral states and state transitions, behavioral states alone, and their differences, which reflect the unique contribution of transitions. n-q, Differences in x-and y-coordinates between the center of mass of neurons that significantly contribute to PLSC1 versus U1 in the GABAergic or glutamatergic population. 1 pixel corresponds to approximately 2.5 µm. Boxplots: medians, interquartile ranges (IQR), and whiskers indicating 1.5× IQR (a-q). See Supplementary Table 1 for details of statistical analyses.

Extended

Extended Data Fig. 7 | Similarity of neural subspaces across different partners and behaviors. a, b. Cosine similarity between PLSC1 coefficients when one animal interacts with different partners (Mean ± SEM). c-f, Chancelevel subtracted cosine similarity between GABAergic PLSC1s (c), GABAergic U1s (d), glutamatergic PLSC1s (e), or glutamatergic U1s (f) established during different categories of behaviors. Asterisks mark the observation groups that are significantly higher than their shuffle controls. Behavioral categories include: attack includes attack, tussle, chase, and threaten; defense includes defend, escape, and flinch; sniff includes general-sniffs, face-sniffs, and genital-sniffs; approach includes approach and follow; non-social includes self-groom, dig, explore object, bite object, climb, and stand. g, Schematics of computing cosine similarity as a measurement of PLSC1 stability during different behavioral feedback from the interaction partner. Created in BioRender. Phi, N. (2025) https://BioRender.com/r9xshzl. h, i, Chance-level subtracted cosine similarity between PLSC1s when a subject animal displays a particular behavior (indicated by the color of the bars) while the partner responds with two different types of behavioral feedback. Asterisks mark the pairs of responses that are significantly more similar than shuffle controls. The results showed that the weights of the shared GABAergic neural subspaces of the subject animals during their non-social moments were more similar between the moment the partner animal was sniffing or approaching and the moment the partner was conducting non-social behaviors (grey bars in h). As animals in both situations were not involved in mutual social interaction, the weights of the shared neural subspace were more similar between the two conditions. However, this does not mean that the shared neural subspace is bigger-the shared neural subspace could be more similar even when it is small. Barplots: mean ± SEM (a-f,h,i). See Supplementary Table 1 for details of statistical analyses. 1 for details of statistical analyses.

Extended Data

Extended Data Fig. 9 | Additional analyses on intra-brain correlations and neuronal ensembles. a, Neural variance explained by the first principal component (PC1) during interaction sessions, social or non-social moments within interaction sessions, and during separation sessions in glutamatergic neurons, compared to temporally shuffled data (subsampled to 50 neurons). b, Intra-brain neural correlation, measured by average PCC between individual neurons and PC1 during interaction sessions, social or non-social moments within interaction sessions, and during separation sessions in glutamatergic neurons, compared to temporally shuffled data (subsampled to 50 neurons). c-f, Scatter plots of intra-brain neural correlation versus inter-brain coupling above chance (PCC of observed PLSC1 minus the PCC of temporally shuffled data) during interaction sessions (c), social (d) and non-social (e) moments within interaction sessions, and during separation sessions (f) in glutamatergic neurons. g, h, Average pairwise PCC between neurons within the same neural ensembles identified during social (g) or non-social (h) moments within interaction sessions in glutamatergic neurons. i, Behavior tunings of glutamatergic ensembles identified during social and non-social moments within interaction sessions. j, Neural variance explained by the first principal component (PC1) during social or non-social moments within interaction sessions, and during separation sessions in GABAergic neurons, compared to temporally shuffled data (subsampled to 50 neurons). k, Average PCC values between individual neurons and their corresponding PC1 during social or non-social moments within interaction sessions, and during separation sessions in the GABAergic population, compared to temporally shuffled data (subsampled to 50 neurons). l, Cross-covariance of PLSC1 in GABAergic neurons in separation sessions after disrupting intra-brain coupling or temporal coupling, relative to the shuffle baseline. m, n, Percentage of significant neurons contributing to GABAergic and glutamatergic ensembles identified during interaction (m) or separation (n) sessions. o, p, Percentage of significant neurons contributing to GABAergic or glutamatergic ensembles identified during interaction (o) or separation (p) sessions in animals that contain at least one ensemble (GABAergic or glutamatergic neurons are subsampled to 50 neurons). One animal (glutamatergic) does not contain any ensembles during separation sessions. q, r, Average pairwise PCC between neurons within the same glutamatergic ensembles during interaction (q) or separation (r) sessions after being subsampled to match the average GABAergic neuron numbers (Mean ± SEM). s, t, Average pairwise PCC of GABAergic and glutamatergic ensembles identified during interaction (s) or separation (t) sessions after subsampling glutamatergic neurons to match GABAergic neuron numbers (Mean ± SEM). Boxplots: medians, interquartile ranges (IQR), and whiskers indicating 1.5× IQR (a, b, j-p). Barplots: mean ± SEM (g, h, q-t). See Supplementary Table 1 for details of statistical analyses.